- N /-'1.1- -- I Wave Propagation and Scattering in Dense Random Media by Paul Robert Siqueira A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering: Electromagnetics) in The University of Michigan 1996 Doctoral Committee: Associate Professor Kamal Sarabandi, Chair Professor Anthony England Professor Fawwaz T. Ulaby Professor John Volakis Associate Professor Andrew Yagle RL-941 = RL-941

,/?/- - 9y Wave Propagation and Scattering in Dense Random Media by Paul Robert Siqueira A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering: Electromagnetics) in The University of Michigan 1996

I

ABSTRACT Wave Propagation and Scattering in Dense Random Media by Paul Robert Siqueira Chair: Kamal Sarabandi This dissertation addresses the important problem of electromagnetic field propagation through and scattering from random media. The mathematical parameter that characterizes the coherent interaction of an electromagnetic field with a random medium is embodied in the fundamental constant of effective permittivity. Effective permittivity relates the small scale interaction of fields in an inhomogeneous medium to the macroscopic behavior of the aggregate. In granular or aerosol media, these inhomogeneities are represented as discontinuities or inclusions, separate from a homogeneous background such as free space. The ability to relate the theoretically understood microscopic behavior of fields in random media to the observed macroscopic behavior is a fundamental problem of applied physics and remote sensing. This dissertation is the culmination of a thorough investigation into this phenomenon and presents a unique, rigorous and consistent method of determining the fundamental parameter of effective permittivity for both two- and three-dimensional random media. In the development of this work, a number of new tools were implemented to model electromagnetic fields in random media. These were (i) a packing algorithm to simulate particle arrangements in gravity-deposited random media such as sand, snow or soils, and (ii) a numerical method for determining effective permittivity independent of particle density, shape, dielectric contrast and arrangement method. The packing algorithm referenced to here, provides essential unknowns such as the correlation function and/or the pair distribution function for use with commonly applied theoretical methods like the Born approximation and the quasi-crystalline approximation. Additionally, the packing algorithm may be used as the first step in the full numerical method for determining effective permittivity. This numerical method surpasses existing theoretical techniques by directly solving the integral form of Maxwell's equations, thereby eliminating approximations employed by theoretical methods to make them tractable. Thus, the numerical method is capable of modeling complex problems often found in nature. As a result, this work represents a solid step forward in our basic understanding of how electromagnetic fields propagate through, and scatter from, random media.

~ Paul Robert Siqueira 1996 All Rights Reserved

To see a world in a grain of sand and heaven in a wild flower, hold infinity in the palm of you hand, and eternity in an hour. William Blake, 1757-1827 English Poet/Painter/Engraver Science and art belong to the whole world, and the barriers of nationality vanish before them. Johann Wolfgang van Goethe, 1749-1832 German Poet/Dramatist/Novelist

This thesis is dedicated to Edir and Mary Christine, Kevin and Thais and in Memory of Eleanor and Joseph. The sheer quantity of words that fill this page, and those that follow, cannot describe all they have given me.

ACKNOWLEDGEMENTS The work that is contained within these pages is the result of the encouragement and friendship of many individuals. While it is impossible to express my gratitude to all of those whose company I have enjoyed and benefited from throughout my graduate career, it would be remiss of me not to recognize the contributions of some of these individuals within the pages of the document that represents the culmination of my academic studies. The first and foremost of my thanks goes to Professor Kamal Sarabandi who has served many functions ranging from an advisor, collaborator and friend. The creative components of my research that are sprinkled throughout the following pages are in a large part due to Kamal's encouragement. His inspiration, persistence and trust gave me the opportunity me to delve into complicated research problems as deeply and completely as I saw fit. His intuition, insight and understanding were always there to greet me at the other end. I can think of no better way for pursuing excellence in graduate education than to follow the example that Kamal has set forth. I would also like to thank and recognize Professor F.T. Ulaby for his early guidance in my graduate career and for inviting me to assist him in teaching an introductory electromagnetics course. His tenure in the field of remote sensing has laid the path for me and many others to pursue interesting and fulfilling careers. During the five years I spent working at the Radiation Lab I have had the chance to encounter and interact with many remarkable individuals. My special thanks go towards Ruben de la Sierra for his undying will, to Tayfun Ozdemir for his unique perspective, to Stephane Legault for his insightful discussions, to Jong Gwan Yook for his graceful manner, to John Kendra for his intellect, to Eric Li for his outgoing nature, to Paul Dahl for his wisdom, to Kathleen Bergen for her encouragement, and to Andrew Oliver for being Andy. In the hallways, offices and SIR-C missions, I have enjoyed many meaningful encounters with Craig Dobson, Leland Pierce, Jasmeet Judge, Josef Kellndorfer, Ron Hartikka, Roger DeRoo, Jim Stiles, Jim Ahne, Edward Kim, Joe Landry, Sunil Bindiganavale, Jian Gong, Youssry Botros and Mike Nurenberger. Much of what I have learned and enjoyed in the field of remote sensing was obtained through thoughtful discussions and experiences with Adib Nashashibi, a respected friend, experimentalist and educator. My sincere wish for him, and all of my acquaintances, is for all of their dreams to come true. The early years of my graduate education were spurned forward by the example of several outstanding individuals in Ames, Iowa. These are Professor John Basart, Professor Phil Appleton, Professor Bill Lord, and my good friend, Ali Safaeinili. As no education can exist in a vacuum from the aesthetic aspects of life, I wish to recog iii

nize Anna Stephanopoulou who has been a paradigm of success and a cherished companion. For the qualities of my character that have inspired me to achieve a doctorate degree I wish to thank the members of my family, beginning with my mother Mary Christine who encouraged my curiosity and my father, Edir Barros, whose many talents have been an example of all that an applied mind can achieve. My brother, Kevin Jay, for his individuality and intelligence, my sister, Thais Maria, for her humanity, and her son, Lake for representing the future. If tears could be preserved within these pages, they would be for the memory of my maternal grandparents, Josef and Eleanor, who will continue to live in all my accomplishments and in the fabric of my life. Finally, I wish to thank the Radiation Laboratory for its diversity and excellence in the field of electromagnetics and remote sensing. All of my expectations have been fulfilled and exceeded within the walls that define this dynamic and respectable organization. iv

TABLE OF CONTENTS DEDICATION............................. ACKNOWLEDGEMENTS..................... LIST OF TABLES.......................... LIST OF FIGURES...................... CHAPTERS 1 Introduction....................... 1.1 Background................... 1.2 Outline...................... 11 iv iv V 1 1 5 2 Measurement of Sand at 35 GHz 2.1 Introduction....... 2.2 Theory.......... 2.3 Experimental Setup... 2.3.1 35 GHz System. 2.3.2 Ground Truth.. 2.4 Calibration........ 2.5 Measurements...... 2.6 Conclusion........ A Motivating *..................... *............... *....... Example -........ 11................ 11................ 13................ 15............... 15................ 19............... 23............... 24............... 3 1 3 Theoretical Methods.............................. 3.1 Introduction............................... 3.2 Coherent Fields............................. 3.2.1 Mixing Formulas........................ 3.2.2 Foldy's Approximation..................... 3.2.3 Quasi-Crystalline Approximation............... 3.2.4 Variations............................ 3.2.5 Pair Distribution Function................... 3.3 Incoherent Fields............................ 3.3.1 Born Approximation...................... 3.3.2 Radiative Transfer....................... 4 Packing Algorithm for Determining Particle Arrangements......... 4.1 2. Packing Algorithm.......................... 4.2 Permittivity Correlation and Pair Distribution Functions....... 33 33 34 34 37 40 48 48 53 54 60 63 67 74 ii

4.3 Conclusions............................... 5 T-Matrix..................................... 5.1 Form ulation............................... 5.2 Limitations of the Recursive T-Matrix Algorithm.......... 5.3 Exam ples................................ 5.4 Variations................................ 5.5 Sum m ary................................ 6 Determination of 6eff in Two and Three Dimensions............. 6.1 Introduction............................... 6.2 Coherent Method for Determining Effective Permittivity...... 6.3 Two-dimensional techique....................... 6.3.1 Initial Results.......................... 6.3.2 Effective Permittivity Dependence on Particle Size..... 6.3.3 Convergence Considerations.................. 6.3.4 Two-Dimensional Particle Arrangement Methods...... 6.3.5 Evaluation of Theoretical Methods.............. 6.3.6 Conclusions........................... 6.4 Three-dimensional Technique..................... 6.4.1 Convergence Considerations.................. 6.4.2 Sensitivity Analysis....................... 6.4.3 Calculation of Effective Permittivity in Three Dimensions 6.4.4 Dependence on Volume Fraction and Particle Size...... 6.4.5 Limitations of the Three-Dimensional Technique...... 6.4.6 Conclusions and Recommendations for the Three-Dimensional Technique................. 6.5 Incoherent Method of Determining Extinction............ 6.6 Summary and Conclusions....................... 86 87 88 95 98 105 107 109 109 110 111 111 116 123 127 130 139 141 142 147 147 150 157 160 161 169 171 171 172 174 177 178 180 185 190 191 195 196 197 199 7 Scattering from Wheat Grain Heads...... 7.1 Introduction............... 7.2 Theoretical Construction........ 7.2.1 Truncated Cylinder Model... 7.2.2 Physical Grain Head Model.. 7.2.3 Electromagnetic Analysis of the 7.3 A Semi-Empirical Grain Head Model. 7.4 Extinction Coefficient........... 7.5 Conclusion................ 7.6 Semi-Empirical Model Coefficients... 8 Conclusions and Directions for Future Work 8.1 Summary of Accomplishments..... 8.2 Directions for Future Work........................... Physical............................................ Grain Head............ *................. Model e............................... BIBLIOGRAPHY.................................... iii

LIST OF TABLES Table 2.1 Ratio of incoherent to total power in dB............................. 31 7.1 Measured physical parameters of wheat grain heads on different dates [Stiles, 1995]........................................ 176 7.2 Dielectric values for the wheat grain head and wheat grains for varying moisture content. The dielectric model is based upon [El-Rayes and Ulaby, 1987]. [77 7.3 Forward-scatter coefficients for vv-polarization....................... 191 7.4 Forward-scatter coefficients for hh-polarization.......................... 192 7.5 Back-scatter coefficients for vv-polarization............................. 192 7.6 Back-scatter coefficients for hh-polarization............................. 93 7.7 Specular-scatter coefficients for vv-polarization...................... 93 7.8 Specular-scatter coefficients for hh-polarization.......................... 194 iv

LIST OF FIGURES Figure 1.1 Random medium volume scattering model. Only the energy scattered from the random collection of particles is observed by the radar in this model. However, the energy incident upon the collection and the energy received by the radar must pass through the homogeneous medium with permittivity ceff. Embodied within this parameter is both the absorptive and scattering losses from the random components that actually make up the homogeneous dielectric............................................................ 7 2.1 Transmission line equivalent model.......................................... 14 2.2 Sensitivity of AO to c"/c' for varying slab thickness....................... 15 2.3 Estimated near-field antenna patterns of the transmit antenna for vertical polarization at three different ranges.............................. 17 2.4 Sand box and network analyzer setup.............................. 18 2.5 Bistatic measurement autocorrelation function...................... 19 2.6 Determination of particle grain size distributions, modeled as prolate ellipsoids. Each grain was modeled as an ellipsoid....................... 21 2.7 Particle grain size distributions for major and minor axes............. 21 2.8 A photograph of the sand surface........................ 22 2.9 Correlation length of 50-70 sand surface............................ 22 2.10 Measured and theoretical AO as a function of incidence angle.......... 25 2.11 Measured and theoretical Ad as a function of incidence angle for a 4 mm sand depth.............................................................. 26 2.12 Measured and theoretical A!~ as a function of incidence angle for a 7mm sand depth.............................................................. 27 2.13 Measured and theoretical Ai as a function of incidence angle for a 13 mm sand depth................................................... 28 2.14 Measured and theoretical AO as a function of incidence angle for a 57 mm sand depth............................. 29 3.1 Electromagnetic model for Polder-Van Santen mixing formula......... 35 3.2 Single particle interaction model for the effective field approximation 37 3.3 Two particle interaction model for the quasi-crystalline approximation... 40 3.4 Domain of integration for the quasi-crystalline approximation.......... 43 3.5 Percus-Yevick radial distribution function for hard spheres. Horizontal axis represents radial axis normalized to the sphere diameter, b........... 50 v

3.6 Particle arrangement method for a classical fluid. (a) 30% volume fraction using sequential addition (dotted circles) and 10% volume fraction resulting from particle extraction (solid circles) (b) 10% volume fraction using sequential addition....................................... 52 3.7 Comparison of the different pair distribution functions obtained from simulations demonstrated in Fig. 3.6............................... 52 3.8 Volume scattering model. Observed power in the backscatter direction arises from the incoherent field scattered from dielectric discontinuities within the volum e....................................... 53 3.9 Elliptical segment traces (0(0i)) of the power spectra] density sampled by observing the backscatter RCS of a random medium. Shown are a number of traces for different frequencies............................. 58 3.10 Sample data and correlation functions of snow [from Vallese and Kong, 1981]. (a) sample cross section of snow (dark regions represent ice, light regions represent air), (b) two-dimensional correlation function of snow sample, (c) horizontal cut of the two-dimensional correlation function (shown is the horizontal correlation length used to fit an exponential function, (d) vertical cut of the two-dimensional correlation function (shown is the vertical correlation length used to fit an exponential function)......................... 59 3.11 An elemental volume in radiative transfer. Scattering of incident power is accounted for by an extinction coefficient, K, and a source function, F.... 60 4.1 Illustration of the two-dimensional packing algorithm. Particles are fit to the surface Sp to find the minimum height, z. Once a particle has been fit to the surface, it becomes part of the surface for the next iteration of the process. 68 4.2 Illustration of the lower surface of a particle, Mlower. The particle is disp cretized into a horizontal grid and the maximum and minimum heights of the particle at points on the grid define the upper and lower surfaces.... 69 4.3 Illustration of equation (4.4) for two different points (il,j2) on the surface S(j). The dilation of the surface is shown for each point on the surface where the maximum is highlighted by a circle. This is the height of where the center of the ellipse would lie for a given horizontal location, j. The minimum of these heights defines the point of minimum potential energy................ 70 4.4 Demonstration of the two-dimensional packing algorithm. Particle radius is 2mm ~ 0.4mm. In this case the volume fraction is 0.8. This particular simulation took 23 seconds on a Sun 10 workstation.....................71 4.5 Demonstration of the three-dimensional packing algorithm. Shown is one slice taken out of a 50mm3 cube filled with ellipsoidal particles with radii of 2mm + 0.4 mm. This simulation using 3300 particles took approximately 17 minutes using 4 processors of an IBM SP Parallel computer. Particles of arbitrary shape may be modeled at minimal additional computation cost.. 72 4.6 Demonstration of the packing algorithm for two-dimensional rocks. Rock surfaces are generated by a random walk about a uniform radius 2mm ~ 0.4mm. This particular simulation took approximately 30 seconds on a Sun 10 W orkstation.................................... 73 vi

4.7 Geometry for a two-dimensional random medium. Medium 0 is considered free space and medium 1 is a random media consisting of particles as described in Figure 4.4................................................. 75 4.8 Two-dimensional correlation function for the particles in Figure 4.4. The correlation function is not axially symmetric................................... 77 4.9 The power spectral density derived from data shown in Figure 4.4. Strong peaks along the diagonal highlight the periodicity in the correlation function seen along this axis. Arcs drawn along a constant radius are proportional to the backscatter intensity as a function of incidence angle....................... 78 4.10 Gaussian (solid line, -) and truncated Rayleigh (dot dashed line,.-) probability distribution functions of particle radii used in the packing algorithm. Both distributions have the same mean.................................. 79 4.11 Slices of Figs. 4.8 and 4.9 taken along the principle axes of (a) 0 0~, (b) 9 450, and (c) 0 = 90~ (next page). The solid lines (-) and the dotdashed lines (.-) show the numerically derived values for the Gaussian and Rayleigh particle size distributions respectively. The dotted lines (...) show the equivalent correlation function and spectrum for the best fit exponential. 81 4.12 Illustration of the visible region of the Born approximation at two different frequencies shown overlaid on a contour map of the power spectral density (Figure 4.8) for the Gaussian particle size distribution. Note that the greatest angular variation due to the packing structure occurs near resonance at f = 21GHz............................................................. 82 4.13 Results of the two-dimensional Born approximation for both 10 and 35 GHz for the Gaussian particle size distribution. Solid lines (-) are for vertical polarization, dashed (- -) are for horizontal polarization and the symbols (o and x) indicate vertically polarized results of the Born approximation using the best fit exponential (Figure 4.11)...................... 83 4.14 Results of the two-dimensional Born approximation at particle resonance (f = 21GHz) for both numerically derived (Gaussian particle size distribution) and exponential correlation functions. Solid lines (-) and circles (o) indicate vertical and dashed lines (- -) and crosses (x) indicate horizontal polarizations. The symbols represent the theoretical exponential curves in both cases.................................................... 83 4.15 The pair distribution function, p(rjj-i) for the packed array shown in Figure 4.4 (Gaussian size distribution). Dark areas indicate low probability and bright areas indicate a high probability of finding another particle at rj given a particle at -r. The inset graph gives the quantitative pair distribution function along the vertical, horizontal and 45 degree axes.................. 85 5.1 Illustration of spherical wave vector translation in the multiple scattering T-matrix equation...................................................90 5.2 (a) scattering from n spheres in the presence of n + n' spheres, (b) scattering from sphere j in the presence of n + n' spheres........................... 92 5.3 Illustration of (5.14).................................................. 93 vii

5.4 Limitations of RATMA. (a) Imaginary component of the diagonal elements of a21 computed using the matrix product Q20-,301 as in (5.24) for three different values of 1r21. For each 1r2l1 (A/10, 2A/10, and 3A/10), the exact value of a21 ( —) is given alongside the approximate values obtained by the matrix product with different values for the maximum number of dipole moments used (see legend). Pmax is given by (5.21). (b) Physical geometry used in the analysis. Shown are the three different particle positions (open circles) in relation to the fixed particle (shaded circle), and the distance, A for calculating Pmax. Particle diameter is A/10............................... 97 5.5 Illustration of the incident and scattered electric field directions in reference to the principle scattering plane defined by the incident vertical electric field, E.........................n................ 9 9 5.6 Magnitude and phase of a single sphere offset from the origin by koT =- 5y. In this example, k0a = 0.5 and c = 6.93 + i0.1, Pmax = 9, Mmax = 1....... 99 5.7 Geometry for two spheres offset from the origin........................100 5.8 Co-polarized vertical and horizontal electric fields scattered from two spheres offset from the origin by k0o = 5. In this example, koa = 0.5 and c = 6.93 + i0.1, Pmax - 9, Mmax = 1.............................. 101 5.9 Two spheres aligned along the z-axis. The position of the inner sphere may be adjacent to the outer sphere, or opposite of it (as shown)...... 101 5.10 Scattered field (vv-pol) from two offset spheres centered along the z-axis. Both plots show results for both the traditional T-matrix approach (labeled "exact") and the recursive algorithm, where Pmax is given................103 5.11 Average error (in dB) over observation angles between the scattered field of two interacting spheres calculated using the recursive T-matrix algorithm (RATMA) and the direct T-matrix approach for vertically and horizontally polarized fields. The variable, Pmax refers to the number of dipole moments kept in the spherical wave expansion as specified by [Bohren and Huffman, 1983]. The inset illustration (center left) describes the geometry of the simulation which specifies d, the diameter of the spheres, and D, the distance of the inner sphere from the origin................................. 104 5.12 Geometry of the 57 sphere composite. Spheres are allowed to overlap so that the total volume enclosed by the individual spheres is the same as the volume of the composite sphere........................................ 105 5.13 Co-polarized vertical (solid lines and circles) and horizontal fields (dashed lines and stars). The Mie series solution is shown as continuous lines and the results from the composite sphere calculated using RATMA are shown as symbols. For koA = 1.26, k0a = 0.33, and e = 2.0 + iO.0, values of Pmax = 9 and Mmax = 1 were used...................................106 5.14 Rotation of coordinate axes from local coordinates (x0, yo, zo) to global coordinates (x, y/, z) by the Eulerian angles of rotation (,/3, 7)..............107 6.1 Model for numerically determining ceff for a random medium...........11 6.2 Sample of a random medium with a volume fraction of 20%. Dimensions for this particular example are given in millimeters for a 48 x 16 mm rectangular slab with 1.6 mm diameter inclusions. The dashed line indicates the boundary that confines the sample of the random medium........... 112 viii

6.3 TM and TE polarized scattered fields for 100 realizations of the random medium illustrated in Figure 6.2. In this example, the observing frequency is 10 GHz and the permittivity of the inclusions is 3.6 + iO.l.............. 113 6.4 TM and TE polarized scattered fields averaged over 100 realizations of the random medium. In this example, the observing frequency is 10 GHz and the permittivity of the inclusions is 3.6 + iO.l.......................114 6.5 Forward and backscatter field magnitude distributions. Numerical simulations (bar graph) follow closely the Rice-Kakagami distribution (solid line). This example is for a square boundary of dimension 1.5A x 1.5A, 40% vollume fraction, TE polarization............................. 115 6.6 Algorithm results for 10% volume fraction. Shown here is a.) a sample collection of particles confined within an imaginary boundary, b.) average scattered TM field (solid line) and the scattered field from a homogeneous dielectric with permittivity, ceff (dashed line), c.) average scattered TE field (solid line) and the scattered field from a homogeneous dielectric with permittivity, Ceff (dashed line).......................... 117 6.7 Algorithm results for 40% volume fraction. See Figure 6.6 for description of individual sub-figures............................... 118 6.8 Algorithm results for 80% volume fraction. See Figure 6.6 for description of individual sub-figures............................. 119 6.9 Error function over the space of realistic ceff. Differences between the mean scattered field of the random medium and a homogeneous dielectric slab are minimized to determine the effective permittivity. To aid in illustration, the difference curve was inverted so that the peak of the displayed surface represents the best fit dielectric.......................... 120 6.10 Simulation and mixing formula results for TM polarization. Real ( ) and imaginary (- - -) components of the effective permittivity derived from the Polder-Van Santen mixing formula compared with effective permittivity obtained by numerical simulation. Symbols indicate different particle diameters of A-/10 (x), 2A-/10 (o), and 3AU/10 (*). Inclusion permittivity is ci = 3.6+il.0121 6.11 Simulation and mixing formula results for TE polarization. See Figure 6.10 for more details................................... 122 6.12 Solution dependence on boundary dimension for TM and TE polarizations. Different lines on the graph represent the average permittivity for a given volume fraction of f = 0.1, 0.2, 0.3, 0.4 or 0.5........................124 6.13 Standard deviation of backscattered field magnitude for three different volume fractions (0.1, 0.3 and 0.5), and both TM and TE polarizations. Similar behavior was observed for the forward scattered field as well.............. 125 6.14 Effective permittivity comparison for differing boundary shapes. Shown are the real and imaginary permitivities for three different boundary shapes: o - circular disk (diameter = 2A-), + - rectangular slab (3A- x Ai) and and x - square box (1.5A, x 1.5A-). The imaginary part of the effective permittivity is multiplied by ten so that both parts of the effective permittivity may be displayed on the same graph. Lines indicate results from the two-dimensional Polder-Van Santen mixing formula.................. 126 ix

6.15 Real part of effective permittivity versus boundary dimension for three different particle sizes: Ai/10 - x, 2A/10 - o, 3,A/10 - * and a volume fraction of 40%, TM polarization.............................. 127 6.16 Simulation of a classical fluid for a volume fraction of 30%. Shown are results from Monte-Carlo numerical simulations and the PY equation in two dimensions solved by Lado........................... 129 6.17 Particle arrangement simulation using particle extraction from a near perfect lattice for a volume fraction of 30%. Shown is the average pair distribution function over all angles and the pair distribution function from the vertical and horizontal directions which accentuate the azimuthal asymmetry. Because the basic structure of the lattice remains unchanged for different volume fractions, the pair distribution function does not change with particle number density........................................ 130 6.18 Simulation of a two-dimensional snow-type medium for a volume fraction of 30%. Shown is the average pair distribution function over all angles and the pair distribution function from the vertical and horizontal directions which accentuate the azimuthal asymmetry....................... 131 6.19 Comparison between five methods of determining effective permittivity for a TM polarized field: EFA (-.), PVS (...), QCA-PY (o), QCA-HC (+) and NUM-PY (*). Plots (a) through (c) illustrate Re(eeff) and (d) through (e) illustrate lm(6eff) as a function of volume fraction using a model of particles (ei = 3.6 + i0.1) suspended in a classical fluid. Particle diameter ranges from /lO10(kd = 0.33) (a and d), 2Ai/10(kd = 0.67) (b and e), to 3A,/10(kd =: 1.0) (c and f)........................................ 132 6.20 Comparison between five methods (QCA-PY, QCA-HC, EFA, PVS and NUMPY) of determining effective permittivity for a TE polarized field. See the caption of Figure 6.19 for details......................... 133 6.21 Comparison between five methods (QCA-PY, QCA-HC, EFA, PVS and NUMPY) of determining the effective permittivity for TM and TE polarized fields. The permittivity of the included particles is ci = 3.6 + iO.01. The real component (not shown) is essentially unchanged from the previous example when i = 3.6 + i0.1. See the caption in Figure 6.19 for details............ 134 6.22 Performance of QCA-HC at very low frequencies. Shown is Im(ceff) for particles with diameter A-/800, A-/80, A-/20, and A-/10 for both TM and TE polarizations (E = 3.6 + i0.1)............................ 137 6.23 Validity regions based on 20% differences between imaginary Ceff for the numerical method of determining ~eff (NUM-PY) and the theoretical methods of QCA-PY (light gray), EFA (hatched) and PVS (dark gray) for (a) TM, ci = 3.6 + i0.1, (b) TM, ei = 3.6 + i0.01, (c) TE,e, = 3.6 + iO.1, (d) TE, ~i = 3.6 + iO.01. The plain white regions indicate where none of the theoretical models yield satisfactory results....................... 138 6.24 Effective permittivity comparison between particle arrangement methods. Volume fraction = 30%, ei = 3.6 + i0.01, a.) TM Polarization and b.) TE Polarization.................................... 140 x

6.25 Vertically polarized scattered fields (magnitude and phase) calculated for different orientations of the incident field and the plane containing the observed scattered field. Shown are 50 realizations for a two wavelength diameter sphere with 1/10 th free space wavelength dielectric spheres (c = 6.0 + iO.1) and a volume fraction of 20%.................................144 6.26 Vertically polarized averaged scattered field using (i) 330 different orientations of the incident field and plane of observation for a single realization of the random medium, (ii) 30 independent realizations of the random medium and (iii) 50 independent realizations of the random medium. See Figure 6.25 for further details........................................ 145 6.27 Scattered field magnitude from 50 realizations of a 20% volume fraction random medium calculated using different values of n' (number of spheres added at each iteration) using the recursive aggregate T-matrix algorithm. When n = N, RATMA is the equivalent of the direct T-matrix algorithm (in this case N = 115).................................... 146 6.28 Scattered fields (vv-polarization) from homogeneous dielectric spheres of varying diameter. Each separate plot displays the scattered field for a dielectric sphere with a real permittivity of 2.1 and imaginary components (going from top to bottom) of 0.0, 0.02, 0.04, 0.06, 0.08, and 0.10 (these are typical values for a 40% volume fraction of dielectric spheres with a dielectric of 6.93 + iO.1)................................... 148 6.29 Scattered fields (hh-polarization) from homogeneous dielectric spheres of varying diameters and imaginary permittivity. See figure 6.28 for further details............................................ 149 6.30 Comparison between the average scattered fields for a random medium with 10% volume fraction (solid line) and the best fit Mie solution for a homogeneous dielectric sphere (dashed line) with ceff = 1.24 + ZiO.02...............151 6.31 Comparison between the average scattered fields for a random medium with 20% volume fraction (solid line) and the best fit Mie solution for a homogeneous dielectric sphere (dashed line) with Ceff = 1.49 + iO.032............. 152 6.32 Comparison between the average scattered fields for a random medium with 30% volume fraction (solid line) and the best fit Mie solution for a homogeneous dielectric sphere (dashed line) with 6eff = 1.73 + iO.044............153 6.33 Comparison between the average scattered fields for a random medium with 40% volume fraction (solid line) and the best fit Mie solution for a homogeneous dielectric sphere (dashed line) with ceff = 1.90 + iO.048............. 154 6.34 Contour and gray-scale plot of the differences between the Mie solution for a dielectric sphere and the averaged scattered field obtained from the random medium as a function of real and imaginary effective permittivity. In this example the volume fraction is 20% and the permittivity of the inclusions is e = 6.93 + iO.1. Dark areas represent areas where the difference is the least between the two solutions............................. 155 xi

6.35 Real refractive index and normalized extinction coefficient comparison between the numerical method (symbols) and the theoretical techniques of quasi-crystalline approximation with coherent potential (solid lines) and the effective field approximation (dashed lines) for three particle radii ranging from 3 to 3.5 mm at 10 GHz. Normalized extinction (shown in the lower plot) consistently increases with particle size for all three methods shown.. 158 6.36 Illustration of the volume fraction and particle radius limits for the determination of effective permittivity using the coherent method in three dimensions. X's illustrate those points that have been calculated in the previous sub-sections (i.e. radius = 3mm to 3.5mm and volume fraction ranging from 10% to 40%). The observing frequency is 10 GHz (A =3 cm)........ 159 6.37 Illustration of a volume current due to dielectric fluctuation in a random medium. Shown are some spherical components of the random medium... 162 6.38 Normalized incoherent power as a function of observation angle for three different square boundary dimensions (0.75A) (), Ai- - ), and 1.5A.(... )) and an inclusion permittivity of i = 3.6 + iO.01. The diameter of the inclusions is A,/10.......................................... 165 6.39 Normalized coherent field as a function of observation angle for three different square boundary dimensions (0.75AX(-),A X( ---), and 1.5Ai(...)) and an inclusion permittivity of qi = 3.6 + iO.01. The diameter of the inclusions is A/10......................................... 165 6.40 Normalized extinction coefficient comparison between the coherent method ( —), the incoherent method ( ---), and the Polder-VanSanten mixing formula (.. ) for two different inclusion permittivities: low-loss, ei = 3.6 + iO.01(*) and lossy, 6i = 3.6 + iO.l(o)................................ 1.66 6.41 Comparison between the coherent (solid lines) and incoherent (dashed lines) numerical techniques for determining normalized extinction coefficient tK,/k for different volume fractions and inclusion volumes at 10 GHz. The three volumes used in the comparison relate to particle radii of 3, 3.25 and 3.5 mm. The dotted line indicates comparable results obtained from the quasicrystalline approximation with coherent potential.................. 168 7.1 Two models for the grain head of wheat. Individual grains may contribute more strongly to the backscatter field than the "smooth" cylinder model most commonly employed................................173 7.2 Schematic diagram of the scattering components of interest for the wheat grain head. Note that the backscatter from the grain layer is unattenuated by other layers in the vegetation.......................... 174 7.3 Electromagnetic comparison of the grain head model with the cylinder model at L-, C- and X-band. Example shown is for Mv = 14%, L = 7cm, andD = l.lcm........................................ 179 7.4 Comparison between RATMA (solid) and the direct T-matrix algorithm (dashed) method of determining the scattered from a wheat grain head at C-Band. Shown is a typical example of a 7cm long head with volumetric moisture, Mv = 14%, consisting of 35 grains whose dimensions are 6.4 x 3.0 x 3.0mm. 181 7.5 Scattering terms of interest for the wheat grain head. Shown are the validity regions most critical for the semi-empirical model................ 182 xii

7.6 Overall comparison of the grain head model (solid lines) with the cylinder model (dotted lines) and semi-empirical model (x) at L-, C- and X-band. Example shown is for Mv = 14%, L= 7cm, andD = 1.1cm............... 186 7.7 Agreement between the semi-emperical model and numerical model for grain head scattering at L-band. Shown are 600 points for back-, forward-, and specular-scattering over the full range of incidence angles (20 < 0 < 70) and grain permittivities (5 < c' < 25, 0 < c" < 10)................ 187 7.8 Agreement between the semi-emperical model and numerical model for grain head scattering at C-band. See Figure 7.7 for details.................. 187 7.9 Agreement between the semi-emperical model and numerical model for grain head scattering at X-band. See Figure 7.7 for details.................. 188 7.10 Imaginary component of Spp at L-Band to be used in the calculation of extinction coefficient. This figure shows results for the grain head model (solid), cylinder model (dashed) and semi-empirical model (x); vv- and hhpolarizations. In this example, L = 7cm, A = 1.3cm2, c = 20 + i6.25(M, = 14%)............................................................. 189 7.11 Imaginary component of Spp at C-Band to be used in the calculation of extinction coefficient. This figure shows results for the grain head model (solid), cylinder model (dashed) and semi-empirical model (x); vv- and hhpolarizations. In this example, L - 7cm, A = 1.3cm2, - = 15+i5.0(Mv = 14%).189 7.12 Imaginary component of Spp at X-Band to be used in the calculation of extinction coefficient. This figure shows results for the grain head model (solid), cylinder model (dashed) and semi-empirical model (x); vv- and hhpolarizations. In this example, L = 7cm, A = 1.3cm2, c = 15+i5.0(Mv = 14%).190 xiii

CHAPTER 1 Introduction 1.1 Background The topic of this thesis is volume scattering in very dense media with discontinuous structure on a similar scale to the observing wavelength. To further explain and understand the necessity of studying such a subject, it is useful to look at three basic spectral bands of electromagnetic waves suitable for remote sensing applications for which the study has direct relevance. These three spectral bands along with a rough wavelength scale are: Optical/Infrared (IOnm - 100,um), Millimeter wave (100/m - 1cm), and Microwave/Radio (1cm - 1O0m). Out of the three regions, optical remote sensing has the longest history, beginning with the use of hot-air balloons replete with an artist and sketchpad flying behind enemy lines during the American civil war. Eventually the eye of the artist gave way to a more complete information gathering system with the advent of the camera. Photographs are analyzed today to determine anything from forest tree density to treaty compliance. The ability of the camera to record differences in color creates a record of the sensitivity of optical wavelengths to molecular resonances and surface roughness. Ancient footpaths can be distinguished with aerial infrared photos because of the recordable differences in temperature caused by organic residues along the path relative to the surrounding area. At optical wavelengths the resolution is typically much greater than the variations we are trying to detect and therefore reflection from regions such as the footpath can be considered to be independent from the reflection from the surrounding area. Optical illumination is typically from an incoherent source such as the sun and interference effects that are familiar feature of lasers is not a consideration in analyzing optical photographs (coherent illumination which is common in microwave and millimeterwave remote sensing will play an important role in the upcoming discussion). The relative shortness of optical wavelengths in comparison to the region under 1

2 study and the fact that optical remote sensing is almost exclusively done incoherently allows approximations to be made to make theoretical analysis relatively simple. Furthermore, because optical processing and imaging is a familiar format to human nature, it is the most developed and mature of the regions to work in and optical sensing is typically accomplished incoherently and does not suffer from the coherent problems of volume scattering that will be discussed in this dissertation. A more recent area, stemming back to the days of World War II, that utilizes coherent fields, is the microwave region. Since the advent of radio at the turn of the century, equipment has been developed to generate and measure energy at radio and microwave frequencies. In these regions the signals can be used to penetrate through clouds and ground cover to detect objects that may not be visible at optical wavelengths. Microwaves can be used to determine bulk dielectric properties and surface roughness on a similar scale to the observing wavelength. Early uses of microwaves were in radar detection and identification of an enemy posing a threat, hence the common use of the term "target" even in today's more peaceful applications of microwave technology. Modern applications of microwaves are used in everything ranging from locating avalanche victims buried in snow to determining chemical composition and surface roughness properties of distant planets. The longer wavelengths used by microwaves allows for measurements of both magnitude and phase of a received signal (i.e. coherent fields), thus doubling the amount of information available. Coherent phase effects also pose the problem of signal fading in radar applications where the received field may vary wildly in a seemingly unpredictable manner. The solution of fading is to average together large numbers of independent measurements of the target under study, thus reducing resolution and information content. Microwave equipment is relatively easy to build with physically simple dimensions and tolerances that can be realized in a manufacturing process, thus reducing the cost and complexity of microwave equipment in comparison to optical instruments. Alternatively, this can prove to be a weakness when it becomes necessary to reduce instrument size, as is the case for satellite applications. Millimeterwave remote sensing bridges the gap between optical and microwave frequencies. Higher energies, compactness and increased resolution make millimeterwave technology an attractive alternative to microwaves; ruggedness, ability to measure phase and reliability are three basic advantages of millimeterwaves over optical wavelengths. Despite these advantages, millimeterwave technology is the youngest of three remote sensing regions mentioned so far. Millimeterwaves suffer from the problems of theoretical complexity of data analysis, lack of available technology and narrow atmospheric windows through which mi

3 crowaves can pass through the atmosphere with little attenuation. As technology advances, the theoretical problems of analysis become more important. The problem of theoretical complexity is not unique to the millimeterwave wave region but it is in this area where the problem becomes its greatest. The depths of this problem will be explained more fully in a short time, but first it should be pointed out that all of the advantages and disadvantages of using one of the three remote sensing regions pales in comparison to the fact that each one of the three regions yields different information about the target under study. The sulfuric acid clouds of Venus provide a good example of the uses of all three of the remote sensing regions. The optical reflectance of sunlight off of the planet's atrnosphere allows us to infer the chemical composition of the clouds that obscure our view of the planets surface. Microwaves, which penetrate these clouds, provide information such as surface temperature which is deduced from observed black body radiation. Reflection measurements made by the Magellan spacecraft at 13cm has provided a complete map of the planets surface and volcanos at resolutions on the order of tens of meters (largely -unexplored from the Magellan data are echoes occurring from the volume laying below the surface). Millimeterwaves can be used to explore atmospheric convection and droplet sizes, information useful to the study of erosion and climatology. All three remote sensing regions provide important pieces of information in understanding the complete picture of Venus as a planet. Similarly, these three regions are necessary in obtaining a complete picture of any "target" and thus each one is important in its own right. As mentioned previously, millimeterwave remote sensing has some theoretical difficulties that limit our ability to explore all of the benefits that this region has to offer. The difficulty arises because for many different classes of targets, particularly natural ones, physical characteristics vary on a similar scale as the wavelength. The concept of inhomogeneity is important in this context. If we consider the response of a deep layer of snow to the three remote sensing regions, we would see that the snow appears as a uniform or layered medium to the lower microwave frequencies. On the opposite side of the spectrum, optical wavelengths do not penetrate the snow layer very deeply; much of what we see is the incoherent reflection of light off ice crystals at the snow surface. Individually, these crystals are much larger than the wavelength of the illuminating field. At millimeter wavelengths, the crystals of ice are no longer at one extreme of the relative wavelength scale, and the incident energy passes into the volume of the snow layer. Thus, penetration depth and grain size become an important factor in the scattering equation. This fact makes the millimeterwave region well suited to explore the crystal structure in a snow pack or the droplet size in a rain storm.

4 Similar to millimeterwaves, microwaves are sensitive to structure on a larger wavelength scale. Microwaves are sensitive to targets such as leaves on a tree or stalks of corn in a farmer's field. The volume scattering phenomena of Bragg diffraction is an example of the sensitivity of coherent light to crystal structure. For an example such as microwaves in vegetation it is common to find that discontinuities or individual scatterers do not make up an appreciable part of the target's volume. Leaves for instance take up only about 10% of the total area occupied by a tree, the tree trunk and branches take up even less; much of the space occupied by a tree is empty air. Similarly in the millimeterwave region, water droplets in clouds and in rain storms are all sparsely scattered in a larger domain. Analysis of this class of problems is made simple because individual scatterers can be considered to react independently of one another to the incident field. The total response of the target can be considered a net sum of the individual scatterers responses to the incident field. As the density of the medium increases, so does the complexity of the volume scattering analysis. If we were to use microwaves to explore the volume of a quarry of rocks (a densely packed volume) the analysis problem is much more complicated than for a sparse distribution of scatterers. Each one of the rocks interacts strongly with the incident field and reflects energy that in turn reacts with neighboring rocks. The use of incident and reflected fields in this situation is no longer straightforward. The problem is exacerbated at millimeterwave frequencies because much of the structure of typical targets varies on the same scale as the observing wavelength. Soils and snowpacks can no longer be considered to be homogeneous nor sparsely distributed. The response of such media to a coherent field in this class of targets is not simple to derive theoretically because basic material properties such as permittivity are not easily determined. The level of complexity incurred in the theoretical analysis increases by a full order of magnitude. The difficulties involved with analyzing the interaction of coherent fields with a dense media returns us to the central topic of this dissertation: analysis of scattering from a very dense media. The purpose will be to provide practical answers to physical problems in dense media volume scattering that are not readily determined by analytic theory. These applications are worthy of illustration through several examples. The first of these may be in the detection of targets buried within a random medium background, such as might be the case in land mine detection. In this instance, it is important to know the scattering and propagation characteristics of the medium surrounding the target so that detection schemes can be optimized to separate the target from the clutter. Modeling

5 efforts have the ability of deepening our insight into a problem and developing a fuller understanding of the physical principles involved. The second of these examples may be in the characterization of a snow field with a millimeterwave radar system mounted on an airplane or part of a satellite system. Because the backscatter response of the radar at these frequencies is dependent on the physical characteristics of the snow field such as density, particle size and wetness, it may be possible to use the radar to estimate these parameters over a large geographic area. In turn, knowledge of these parameters could then be used to predict water resources or as a component in earth biosphere models such as global heat budget analysis. 1.2 Outline This thesis begins with the description of a controlled experiment whose purpose was to determine the effective permittivity of a layer of fine sand at a frequency of 35 GHz (A = 8.54mm). At this frequency, the individual granules of sand exist on a similar dimensional scale as the observing wavelength and therefore it is expected that individual grains interact strongly with the incident field and with each other within the sand layer. Aside from the immediate problem of determining the fundamental parameter of permittivity for the layer of sand, the study of the interaction of electromagnetic fields with this medium poses an interesting question: "How do electromagnetic fields propagate through a highly inhomogeneous medium that is discontinuous on the same dimensional scale as the observing wavelength?" It is this question that has motivated the work that makes up the theme of the dissertation that follows. In response tohat has been stated above, a number of the question that has been stated above, a number of theoretical methods have been developed by others to determine the parameter of effective permittivity (the parameter that characterizes how a coherent electromagnetic field propagates through random media). Embodied within this parameter is the phase delay that the field encounters (the real component of effective permittivity) and the amount of losses encountered by the coherent field, both absorptive and scattering (the imaginary component of effective permittivity). Chapter 3 of this dissertation addresses a variety of these theoretical methods with the objective of developing the theories that will be referred to in future chapters. Chapter 3 is split into two domains, one that addresses the coherent interaction of electromagnetic fields with random media and another that addresses the incoherent fields. Coherent interaction refers to the phase delay and loss of an electromagnetic field as it propagates

6 from one end of a random medium to another. The characterization of this phenomenon is fundamental to relating the phenomenon of volume scattering to surface scattering and is also the vital first step in developing the subject of the second half of Chapter 3, that of characterizing incoherent fields in a random medium. It is the incoherent component of a field that a radar actually observes from the volume of a scatterer. This energy is termed incoherent because the random nature of the scatterer volume has the effect of randomizing the phase of the reflected signal. Thus, while the average field scattered from this medium may be zero, the power, or the variance of the signal, is non-zero and is proportional to the physical and electromagnetic characteristics of the random medium. This concept is worthy of further illustration. To begin, we take an idealized situation of a random collection of scatterers buried within a large, homogeneous dielectric and ask ourselves what would be the scattered field observed by the radar (Figure 1.1). Because the surface of the homogeneous medium is flat, at any angle other than the nadir direction, the surface would have no contribution to the backscattered field. Some of the energy incident on the medium however would be reflected from the surface in the specular direction and some would be propagate into the homogeneous dielectric (a process theoretically described by Fresnel reflection/transmission coefficients). Because the homogeneous dielectric has no discontinuities, there would be no energy reflected back to the radar until the field reached the buried random collection of scatterers. At this point, we might consider the interaction of the incident field with the random medium as a separate problem and try to determine what part of the incident field would interact with the medium and how much of that part would be reflected back towards the radar. Once this has been determined, if the permittivity of the homogeneous dielectric were known, the energy observed by the radar could be determined directly. At this point, we can step away from the idealized model of the random collection of scatterers buried within a homogeneous dielectric, and ask ourselves what does this homogeneous dielectric actually consist of? The answer to this is that it consists of random components just like the one that was embedded in the simplified model. Much of the energy (but not all) that was incident upon the random medium sample actually passes through into the forward direction and continues to propagate. The energy that is lost by scattering the incident fields into other directions appears directly in the imaginary component of the effective permittivity that characterizes the behavior of coherent fields as they propagate through the medium. Thus, knowledge of the effective permittivity Ceff encompasses both the coherent and to some extent the incoherent behavior of fields within

7 Figure 1.1: Random medium volume scattering model. Only the energy scattered from the random collection of particles is observed by the radar in this model. However, the energy incident upon the collection and the energy received by the radar must pass through the homogeneous medium with permittivity ceff. Embodied within this parameter is both the absorptive and scattering losses from the random components that actually make up the homogeneous dielectric.

8 the medium. Furthermore, accurate knowledge of what the value of effective permittivity is for a given medium (i.e. accurate characterization of coherent fields) is fundamentally important to understanding how incoherent fields are excited by the incident field and how they behave within the medium as they propagate back to the radar. The theoretical methods that are described in Chapter 3 address this issue, but because of the complexity of the electromagnetic interactions involved, it is difficult to ascertain what the limitations of the different theoretical methods are. In this light, Chapters 4 through 6 endeavor to address the problem of developing an independent, numerically exact solution, for determining how fields propagate within random media. Chapter 3 addresses the physical aspect of this problem by first developing a method for simulating the manner in which discrete particles may be arranged together in a dense medium. This method is referred to as the "Packing algorithm" that was developed initially to mimic the sand experiment described in Chapter 2. The algorithm however has a variety of other applications, such as providing correlation and pair distribution functions (in both two and three dimensions) for use with some of the theoretical methods developed in Chapter 3. Knowledge of discontinuity locations in a dense media also has the application of being the first step in a full Monte-Carlo simulation of electromagnetic fields in dense random media. To this end, the method of moments may be used to efficiently analyze the behavior of electromagnetic fields in two dimensions, or for three dimensions, the T-matrix (transition matrix) algorithm can be used in three dimensions. This approach, less common than the method of moments, is described in depth by Chapter 5. Given a random arrangement of particles (obtained either by the packing algorithm or a different method), and employing a numerical solution technique to determine the scattered fields, it is shown in Chapter 6 that the fundamental parameter of effective permittivity can be obtained. This is accomplished by comparing the average scattered field obtained from the samples of a random medium with that of a uniform dielectric body whose shape and size is the same as the random medium samples. It is shown that this technique provides reliable and consistent results in both two and three dimensions. Given the technique as it has been developed, it is thus possible to return back to the theoretical methods of characterizing coherent fields within a random medium and to evaluate these methods based on results obtained from the numerical technique. The techniques and methods that have been developed in the previous chapters can be used for addressing a variety of different problems in electromagnetics and environmental remote sensing. The penultimate chapter of this dissertation (Chapter 7) illustrates one such problem: that of how to characterize the scattering from wheat grain heads. This

9 problem is particularly important at the microwave frequencies of C- and X-band, where the observed backscatter from mature wheat plants can be shown to be strongly influenced by attenuation and contributions from the grain component of the vegetation layer. Chapter 8 reviews the work that has been accomplished by this dissertation and gives suggestions and directions for future work. In summary, the contributions described by this thesis are: 1. Characterization and bistatic measurement of fine sand at 35 GHz. Measurements are used to determine the real component of effective permittivity for the medium.. 2. Development of the two-dimensional version of the quasi-crystalline approximation. This treatment is instructive in its approach and provides a basis/alternative view for the three dimensional problem developed by other investigators. 3. Development of the two-dimensional version of the Born approximation and graphic relation of the predicted backscatter to the power spectral density of the permittivity fluctuation correlation function. This treatment illustrates how the backscatter and particle size distribution are related. Further discussion addresses the problem of obtaining the correlation function and the variety of correlation/power spectral density functions that are valid. 4. The development of the packing algorithm for efficiently determining particle arrangements in two and three dimensions. Results are used to illustrate applications to the Born approximation and the numerical evaluation of the quasi-crystalline approxirnation. 5. Implementation and evaluation of the recursive aggregate T-matrix algorithm [Wang and Chew 1993]. This algorithm was needed for the three-dimensional treatment of effective permittivity. Independent implementation of this algorithm led to several important discoveries such as errors in the original publications of this material and limitations of the algorithm not discussed in these publications. This algorithm served as the driving numerical method for determining effective permittivity in three dimensions described by the following item. 6. Numerical determination of effective permittivity in both two and three dimensions. As with any method, the numerical method has its limitations, but insofar as computational resources are made available, a reliable and consistent method is set forth for determining this fundamental quantity.

10 7. Evaluation of the quasi-crystalline approximation, the effective field approximation and the Polder-VanSanten mixing formula in both two and three dimensions. The numerical method for determining effective permittivity is used to analyze similar results predicted by these three fundamental theoretical techniques. 8. Development of a semi-empirical model for the scattering and extinction due to wheat grain heads. This model is illustrated as an application for some of the tools that have been developed for other purposes in this dissertation and highlights a technique for addressing other similar problems that may be encountered in electromagnetic analysis in environmental remote sensing.

CHAPTER 2 Measurement of Sand at 35 GHz - A Motivating Example - 2.1 Introduction The following chapter gives a quantitative description of a set of bistatic measurements that were used to determine the dielectric constant of a layer of sand at 35 GHz. At these frequencies, scattering from and interaction between individual grains of sand becomes a significant factor. The motivation behind this work is to gain a better understanding of mechanisms of surface and volume scattering for distributed targets, with the ultimate goal of relating electromagnetic analysis to radar observations in the backscatter direction. The focus of this chapter, and much of the dissertation that follows, will be on the effective permittivity. For surface scattering, the effective permittivities of two adjacent media are used to calculate the Fresnel reflection coefficient and thus determines the magnitude and phase of the reflected wave. In volume scattering the scattering of a wave is caused mainly by discontinuities within the volume. For a mixture of discrete particles, the theoretical determination of effective permittivity is a statistical problem where assumptions must be made about interactions between fields and particles within the volume. Low frequency equations relating the effective permittivity of a mixture to the dielectric constants of its components are called mixing formulas, for which there are a variety (Tinga et al., 1973; van Beek, 1967). The experimental determination of effective permittivity provides a method of checking which of these formulas is most appropriate and in this way provides information about the nature of electromagnetic interactions and scattering mechanisms present within the volume. The permittivity of a material describes a fundamental relationship between the applied electric field and the behavior of charges under the influence of this field. Because this relationship is so fundamental, there are in fact a large number of methods for measuring this quantity and it is the equipment availability and simplicity of measurement that ul 11

12 timately determines the method used to find the permittivity at a specific frequency. As an experimental target, sand is very useful because electromagnetically important characteristics about the medium can be easily controlled (i.e. grain size, wetness, and purity). The method for determining the dielectric constant described in this chapter is designed to be a simple extension to these measurements. Aside from the determination of dielectric constant for a statistical target, this measurement also has benefits such as providing checks for the operation and calibraations of the sand on a submillimeter scale. Polarimetric bistatic measurements themselves are not often found because of the complexity and detail involved in the setup of the problem. One recent bistatic measurement of sand was made by Ulaby et al. (1987) where the bistatic scattering coefficient was measured as a function of azimuth angle for a fixed angle of incidence (0i = 66~) at 35 GHz. The experiment used a homogeneous mixture of common sand with a reported dielectric constant of c = 2.5 + i0.03, measured at 10 GHz. Measurements were made of the scattering coefficient for three types of surfaces: smooth, rough, and very rough. Results showed that the measured reflection coefficient from a smooth sand surface agreed well with the predicted value for a semi-infinite medium and that surface surface scattering increased as a function of surface roughness. In contrast, this chapter describes a controlled experiment of forward scattering from a layer of sand whose surface has been made electromagnetically smooth. Observations are always made in the specular direction as a function of sand thickness and angle of incidence. Efforts were made to assure minimal contributions of incoherent surface and volume scattering in the forward direction. A refined sand (US Silica ASTM 50-70) of uniform size distribution and 99.8% pure silica (Si02) was used. The 50-70 designation indicates that all particles passed through a wire grid of 50 squares to the inch, but did not pass a similar grid of 70 squares to the inch (i.e. sand particles are constrained to have diameters between 0.36mm and 0.51mm). The smoothness of the sand surface, the grain sizes and correlation length are physically quantified to submillimeter accuracy (such accuracy is necessary at 35GHz where the free space wavelength is 8.57mm). The sand layer is modeled as a mixture of sand particles surrounded by an air medium and using a measured volume fraction of the sand particles, the Polder-Van Santen mixing formula is used to determine the effective permittivity to fields within the sand layer. This mixing formula consists of two unknown permitivities, the effective permittivity of the distributed mixture (which is the objective of this experiment) and the apparent permittivity of the

13 mixture to fields within the sand layer. Typically for sparse mixtures, this latter permittivity is set to that of free space and for dense mixtures it is set to the effective permittivity. The results of this study however indicate that the permittivity apparent to fields within the layer is actually closer to that of the sand particles themselves and in this sense implies strong electromagnetic interactions between particles. This chapter begins by deriving a function relating the phase of the scattered electric field to the permittivity and thickness of a dielectric slab. An equivalence is made between the plane wave model of the the physical situation and an equivalent model using transmission line theory. This is followed by a description of the 35 GHz bistatic system that was used in the measurement, along with a method for calculating the co-polarized system constants. It is then shown that the sand being used as a distributed target possesses the physical characteristics necessary that allow it to be equivalently modeled as a dielectric slab. A description of the least squares model fitting routine is given along with calculated results for the effective dielectric constant of the sand. 2.2 Theory Transmission line theory can be used to draw an equivalent model for a plane wave incident at an oblique angle onto a dielectric surface (Ulaby et al., 1981; Ramo et al., 1984). Through the model, a reflection coefficient can be calculated using standard transmission line techniques. Specifically, a model will be derived here for a dielectric slab placed on top of a perfectly conducting plane. From the experimental setup shown in Figure 2.1, it is clear that the magnitude of the reflected wave will be near unity for low loss dielectrics (i.e. energy is not dissipated in the system). For a solid plate of silica glass, the complex permittivity at 10 GHz is reported to be Ec = 3.78 + iO.O01 (Holloway, 1973), giving an effective depth of penetration for the glass plate of 6p = Ax/~/27re" = 2.65m (where, c = c' + ic"). Assuming the dielectric losses of the sand and glass plate are similar, it is then expected that most of the energy incident on the sand layer will not be dissipated and therefore the phase of the reflected wave will be the most sensitive to the dielectric constant of the sand layer. Depending on the polarization of the incident wave, the impedances of the equivalent model can take on one of two values: Z rfn cos 0n vertical polarization (2.1) rn sec On horizontal polarization

14 grouna plane Z3 Figure 2.1: Transmission line equivalent model for a coherent wave scattering from a dielectric surface over a ground plane. where On is the intrinsic impedance of the nth medium and 0n is the angle of incidence in medium n. When n 5 1, t, for low loss dielectrics can be obtained from Snell's law: On sin-1 (- n a-. (2.2) The reflection coefficient at each interface is calculated by n kZn+ - Zn here f 1 vertical polarization (2.3) Zn+i + Zn 0 horizontal polarization. The reflection coefficient for the entire system is then Er R1 + R2-_2t(a2 sec 02+j32 cos 02) Re Ei 1 R1 R2e-2t(a2sec02+j32cos02) (2.4) where a2 = 2fjm{ /c}f and I32 = je{/} are the dielectric loss and propagation factors, respectively, t is the thickness of the dielectric slab and 02 is the angle of incidence in the sand layer. For a reference incident electric field E{ = Eo(l + iO), the reflected wave can be written as o -|RV V or =|Rh |Zh. (2.5) Eo Eo The subscripts v and h indicate vertical and horizontal polarizations and the angle X refers to the phase of the complex reflection coefficient, Re. Now using Xv as the reference phase, the difference in phase between the two polarizations becomes Aq = ( &h - v (2.6) where Aq$(c, t, 0) is a function of complex dielectric constant (cc), slab thickness (t) and incidence angle (i). Computer analysis has shown that AO is not sensitive to values of the

15 imaginary component of ~c when c' > 100l". This can be seen in Figure 2 where a plot is made of the sensitivity of the phase difference function (2.6) to the ratio of "/c'. lFor a given c' = 2.6, and varying slab thickness (t = 5, 10, and 20mm), e" was increased in logarithmically spaced steps from 0.001 to 10. The average difference between successive Ad functions calculated for each c"/c' is plotted to give a measure of the ability of the method to detect changes in the imaginary component of the dielectric constant. As can be seen from the figure, the ratio c"/c' should be greater than 0.01 to detect differences in '" for sand depths less than 20mm. Average phase change as a function of increasing ~".: ct CA 0;o Uc E 0 ct c) Q) c3 60 50 40 30 20 10 0 L — 0.001 0.01 0.1 1 10 ratio: E"/E' Figure 2.2: Sensitivity of AO to c"/c' for varying slab thickness (t = 5, 10 and 20mm). The vertical scale is a measure of the average difference of AO for increasing values of c". For this plot, ' = 2.6 remained constant. 2.3 Experimental Setup 2.3.1 35 GHz System The sand measurements were made using a 35 GHz bistatic network analyzer based scatterometer. The output port of a HP8753C Vector Network Analyzer was used to mod

16 ulate an IF signal from 2 to 3 GHz which was in turn used to modulate a 32 GHz RF signal from a x3 harmonic mixer (10.67 GHz => 32 GHz), the resulting signal varying between 34 and 35 GHz. A 1 GHz bandwidth gives a range-resolution capability of the system of Ar - = 15 cm. (2.7) 2B The actual resolution used for the measurement was actually larger than this because software gating on the network analyzer was used to isolate the response of the target from system and room noise. A typical gate of 15nS was used which gives the equivalent resolution of 2.3 meters. All sand depths used in the measurements were much less than this, making the measurement an integrated effect of the reflected wave from the top surface of the dielectric and the multiply reflected wave from the sand layer. This was implicitly assumed in the preceding theoretical derivation. The transmitter antenna (circular aperture, d = 15.25 cm) half power beamwidth at 2.5 meters was measured to be 5.4 degrees and the receiver (rectangular aperture, 5.1 x 5.1cm), has a half power beamwidth greater than 10 degrees (/301/2 > A/lx,y). With the receiver having a larger beamwidth than the transmitter, pointing of the two antennas becomes less critical and we can often make the assumption that the receiver gain is constant over the HPBW of the transmitter. At a distance of 2.5 meters, the transmitter is being operated close to the near field region (D2/A = 2.7m). Tests have shown however that within the range of 2.3 to 2.7 meters, the phase and amplitude of the transmitter pattern remain fairly consistent (Figure 2.3). These patterns were made using a receive antenna of the same dimensions as the transmitter to detect reflected fields from a small trihedral placed at three ranges. The horizontal axis angles were doubled to adjust for the multiplication of the transmit and receive antenna patterns and the received power was corrected to account for free space propagation loss. The similarity of the magnitude and phase patterns for the three ranges indicates that targets in the range of 2.3 to 2.7 meters may be treated as being in the far field. A large 1.8 x 1.8m rotatable box was used to hold the ground plane and sand layer. Once the transmitter and receiver were set at their appropriate locations, the received signal from the ground plane was maximized by making fine adjustments to the angle of the receiving antenna. A brief illustration of the entire setup is shown in Figure 2.4. Illustrated also is the HP8753 network analyzer where the RF output port provides a signal to the transmitter, from which a sample of that signal is returned to the reference port of the network analyzer. The received vertical and horizontal polarization signals are separated by an ortho-mode transducer placed directly behind the receiving antenna, and the down-mixed IF signals

17 Near field magnitude, vertical polarization -15 —.... =- ---- 2.3m -2 0 2.7m ----- -25 - ', -40 - -45 - -50 -55 I. -2 0 2 4 6 8 10 Azimutth (degrees) (a) Figure 2.3: Estimated near-field antenna patterns of the transmit antenna for vertical polarization at three different ranges (2.3, 2.5, and 2.7m). Plot (a) is the magnitude and (b) on the following page, represents the phase. are sent to the A and B ports of the HP8753 so that both received polarizations may be detected simultaneously. The measured signal from the sand is modeled as the summation of a coherent reflected field, Ecoh- from the smooth specular surface of the sand and an incoherent field, Einci resulting primarily from scattering within the sand volume: Em = Ecoh + E-in (2.8) The nature of the incoherent field is such that it is a random process with zero mean ((Einc) = 0) and variance proportional to the power of the incoherent field (Pinc) (Einc)2 = ((Em - Ecoh)2) (2.9) The coherent field is deterministic as given by (2.4) and (2.5) and can be estimated by finding the sample mean of the measured fields Ecoh = (Em) = (Ecoh + Einc) Ecoh (2.10) The coherent electric field is estimated by taking a large number of samples of the target and averaging the measurements. Each sample must be statistically independent of the other

18 Near field phase, vertical polarization U) G) a1) a1) Q) U) rd 180 170 160 150 140 130 / 2.3m phase 2.5m phase 2.7m phase ----, I i// ', - /f:," IA: -~.I A.. 120 - 110 - 100 -2 0 2 4 6 Azimuth (degrees) 8 10 (b) Figure 2.3 cont'd: Estimated near-field antenna patterns of the transmit antenna for vertical polarization at three different ranges (2.3, 2.5, and 2.7m). Plot (a) from the previous page, is the magnitude and (b) represents the phase. for tterero. Since the incoherent field to averagdue to zero. Since the incoherent field is due to the random locations of scatterers in the antennas field of view, the target must be moved a sufficient amount such that successive measurements of the target are uncorrelated. This movement is achieved by rotating the sand table a predetermined angle between each measurement. The angle of rotation is determined by autocorrelating a sequence of measurements taken of the received signal for small angle increments over a full rotation of the table. The change in angle at which the autocorrelation function drops below a value of 0.2 is taken as a sufficient Figure 2.4: Sand box and network analyzer setup.

19 Bistatic correlation 0.3 4-i C1) 0) 0 u 0 H 4.) 0) 0 U 0.25 0.2 0.15 0.1 0.05 0 2 4 6 8 10 12 14 16 18 20 table rotation angle Figure 2.5: Bistatic measurement autocorrelation function. A set of 270 measurements were taken at one degree increments of the table and the magnitude of the electric fields were autocorrelated to determine the minimum angle of rotation required to achieve statistical independence. increment angle to provide the required degree of decorrelation between measurements. This was done for a sand depth of 7 mm and incidence angle of 21 degrees using one degree increments of table rotation. The resulting autocorrelation function is shown in Figure 2.5 where it can be seen that the statistical correlation between adjacent measurements dropped to less than 0.2 for 4 degree angles of rotation. Thus each complete measurement of the sand for a given incidence angle and sand depth consisted of 90 measurements. The number of samples actually used in determining the coherent field for any given incidence angle and sand depth was slightly less than this because it was determined that a seam in the ground plane was causing non-stationarities in the statistics over a small range of table angles. Sample sizes more typically consisted from 70 to 80 samples for this reason. 2.3.2 Ground Truth In order to model the electromagnetic scattering from sand, the following physical characteristics must be considered: 1. Correlation length and rms height of the rough surface.

20 2. Particle size and shape distributions. 3. Volume fractions of the constituent particles making up the distributed target (this includes solids, liquids and free space) and their respective dielectric constants. Because the sand layer is being modeled as a dielectric slab, volume scattering contributions in the specular direction must be much smaller than the coherent component. Volume scattering is dependent upon the particle size (relative to wavelength) and the dielectric contrast (permittivity of the particle relative to the background material). For particles whose radius, a, is much smaller than the wavelength, A, the Rayleigh approximation to Mie scattering for dielectric spheres can be used (Ulaby et al., 1981). The scattering cross section, Qs is the ratio of scattered power, Ps, to incident power density, Si. Using the Rayleigh approximation gives: PS 2A2 (27ra6 C-12 Q Si r( A) c+2 (2.11) The sand was characterized by photographing a random sample of 143 sand particles and modeling each of the particles as an ellipsoid (Figure 2.6). Measurements of the major and minor axes (Figure 2.7) showed that the mean major axis length was 0.37 mm (std. dev. = 0.05 mm) and the mean minor axis length was 0.28mm (std. dev. = 0.04 mm), which gives a mean eccentricity for the ellipsoids of e = (1- (mi2)2)1/2 0.65. Using a worst case particle size of 0.5mm in (2.11) gives an upper limit to the scattering cross section of Qmax = O.lxO-9m2, or an efficiency factor of smax = Qs/Tra2 = 0.5x10-3. Thus it is likely that scattering by the individual particles (i.e., volume scattering) although present, will not be a significant factor in the forward scattering direction. The surface of the sand was made smooth by using a straight metallic edge so that the sand layer may be modeled as a dielectric slab with smooth boundaries. Because of the small wavelength being used, it was necessary to quantize the smoothness of the surface on the order of 1/32 of a wavelength, 0.3 mm (the Fraunhofer criterion). To do this, a sample of sand was placed between two glass plates and the surface was flattened using a paper blade. Photographs taken of the sample's surface under a microscope showed that the mean variation in the surface smoothness was better than 0.3 mm when sampled every 0.17 mm on the surface (Figure 2.8). The conclusion from this measurement therefore is that the surface can be considered electromagnetically smooth at 35 GHz. The correlation from on from one point on the surface to another, the correlation length 1, is also an important parameter when considering the scattering from a surface. For this

21 Figure 2.6: Determination of particle grain size distributions, modeled as prolate ellipsoids. Each grain was modeled as an ellipsoid Maior axis size distribution for 50-70 sand (143 grains) ct 0 e0 o (O 40 40 - 30 20 10 -10 0.15 0.2 0.25 1 4 —l 1 1 1 1 1 1 1 0.45 0.5 0.55 0.6 0.3 0.35 0.4 length in mm Minor axis size distribution for 50-70 sand (143 grains) 40 30 20 10 0 0.15 v I i ----— I --- —-....... 0.2 0.25 0.3 0.35 0.4 length in mm 0.45 0.5 0.55 0.6 Figure 2.7: Particle grain size distributions for major and minor axes.

22 figure 2.8: A photograph of a sand surface sandwiched between two glass plates and flattened by a paper blade. experiment, it is important to minimize the scattering contribution from the surface. Using the sampled surface discussed above, an autocorrelation plot of the surface was made (Figure 2.9). The result shows that the correlation length is smaller than the rate at which the surface was sampled (0.17 mm) which further supports the assumption that the surface is very flat. 1 C 4U u C U - r. 4 —j ct t: 0 I) 0.5 0 Correlation length for 50-70 sand surface........:...:.....: i ~~~ I -......... 0 0.5 1 1.5 2 2.5 3 3.5 Distance (mm) 4 4.5 5 Figure 2.9: Correlation length of 50-70 sand surface.

23 2.4 Calibration Looking in the specular direction at a reflected signal from a flat surface with low volume scattering, it is expected that most of the energy received is the same polarization as the transmitted signal. Any cross-polarized response from the distributed target will give a measure of the combined effect of 1. Misallignment of the transmit/receive polarization vectors 2. Cross talk between the vertical and horizontal channels in the antenna structure, and 3. Scattering by the volume of the target. The system contributions to the cross-polarized signal (items 1 & 2) were measured to be approximately 15 dB down from the transmitted polarized wave while the cross-polarized signal due to the target alone (item 3 from above) is better than 20 dB less than the transmitted signal (i.e., measured volume scattering is negligibly small). The difference between the v & h channels for the transmit and receive combination need to be calibrated so that only the target characteristics are measured. Because for this experiment we are only interested in the co-polarized response (namely the phase differences between h and v), a simple calibration technique can be employed using a ground plane [Sarabandi et al, 1990]. Here, the measured signal for both co-polarized responses is modeled as K~ - e-k(rt+rr)R T and h -K ik(rt+rr)RTS~ rtrr rtrr (2.1 2) where S~ and Soh are the scattering coefficients for a ground plane, rt and Tr are the transmitter/receiver distances to the target, k is the wave number, K is a system constant and the RpTp pairs represent transmit and receive constants for each set of co-polarized responses, p = {v, h}. By measuring the reflected fields from a ground plane (EI E'j), calibration constants can be calculated by K - K -jk(rt+rr)R T- VV and Ahh K je-k(rt+rr) RhT -_ hh rt rr S~v r Sth (2.13) where So = 1 for vertical polarization and Sh = -1 for horizontal polarization. Measurements of an unknown target (EVV or Ehh) are calibrated to give the unknown target's scattering coefficients S'v and Sh by Sv ~ and Sh - E~ h (2.14) vv Khh

24 which assumes the distances are not changed. A very strong advantage of using the phase difference function (2.6) to determine the effective permittivity is that it removes the necessity for the system to be phase stable over the period of observation. What must remain constant is the phase difference between the vertical and horizontal polarization channels. The transmitted polarization is set by phase plates at the antenna aperture where the mechanical position of the plates determines the transmitted polarization. If the positioning of the plates remains consistent for the duration of the experiment, the calibration will correct for phase differences caused by the positioning of the plates and the measurement will not be dependent on phase and amplitude drifts of the transmitter. 2.5 Measurements Measurements were made of the sand using the setup previously described (section 2.3). Sand depth was varied between 4mm and 6cm. At each depth, a set of six, 90 sample measurements were taken for incidence angles of 21, 28, 44, 56, 61 and 70 degrees from vertical. These angles were determined by constraints on the physical positioning and pointing of the radar. From these data, the sample mean of the phase difference A(c, t, Oi) was found for each depth and angle of incidence. A least mean square fitting routine was then used to fit the model to the measured data. It was found that the sample variance increased as a function of incidence angle (i.e. the uncertainty of the measurement was different for each angle of incidence). It would be wrong then to equally weigh differences between measurement and theory for all angles of incidence. Thus, the least squares algorithm was made sensitive to these differences by incorporating a weighting factor that made the model/measurement errors linearly related to the standard deviation of the measurement: 2 /measurement - model a2 error2 = (2.15) std. dev. The sum of these errors constitutes the least squares measure, which was minimized over a range of dielectric constants min E error2(c, t, (2.16) where the summations are over angle of incidence and different sand thicknesses for which measurements were taken. The dielectric constant giving the least error is then the best estimate of the effective dielectric constant of the sand. Using a set of five sand depths

25 17mm measured and theoretical phase difference -H,I, 4 4,ill 4 -H' 4 0 r., (I) ci) H a) (ci 04 360 340 320 300 280 260 240 220 200 180 160 140 120 100 0 10 20 30 40 50 60 Angle of incidence (degrees) 70 80 90 Figure 2.10: Measured and theoretical AO as a function of incidence angle for two values of dielectric constant and a sand depth of 17 mm. The solid line represents a value of c = 2.64 + i0.0 (the best fit model) and the dashed line represents the mixing formula (2.17) value of e = 2.4 + i0.0; individual points are measured values of phase difference for which the error bars indicate ~d 1 standard deviation from the measured mean. (4 mm, 7 mm, 13 mm, 17 mm and 57 mm) a value of c, = 2.64 + iO.0 gave the best fit between the model and the measured data. The results comparing measured phase difference with theoretical phase difference (2.6) for a sand depth of 17 mm is shown for an example in Figure 2.10; Figures (2.11) through (2.14) contain similar figures for the remaining depths measured. Note that the dielectric loss factor c", while likely present to some degree, was not detectable by this method (see section 2). The measured value of dielectric constant discussed above compares well with dielectric constants theoretically calculated from a mixing formula. The mixing formula calculates the effective dielectric constant in terms of the measured physical parameters of the sand (ground truth). For randomly oriented ellipsoidal inclusions, the most common mixing

26 I.H 04 1 H a a) - U) 04 (D -11 (D Ul a) 44 44 360 340 320 300 280 260 240 220 200 180 160 140 120 1 nn 4mm measured and theoretical phase differences I I I I I I I I ----------------- I I I I I I I I II I... I. Jl u U 0 10 20 30 40 50 60 Angle of incidence (degrees) 70 80 90 angle for a 4 mm sand Figure 2.11: Measured and theoretical AO as a function of incidence depth

27 7mm measured and theoretical phase differences 1 -a) a) rd r.1 360 340 320 300 280 260 240 220 200 180 160 140 120 1 nn {/ - - /. i I " --- —----- -— iE___ " < """"" J ~ i I i! I I I I I I -l v v 0 10 20 30 40 50 60 Angle of incidence (degrees) 70 80 90 Figure 2.12: Measured and theoretical A> as a function of incidence angle for depth a 7 mm sand

28 13mm measured and theoretical phase differences I.H 12 a, a) 260 240 220 200 180 160 140 120 100 80 60 40 20 0 0 10 20 30 40 50 60 Angle of incidence (degrees) 70 80 90 Figure 2.13: sand depth Measured and theoretical AOq as a function of incidence angle for a 13 mm

29 57mm measured and theoretical phase differences I I I I I I I -H -ri a) U 0) -rl 0) U) Q) 280 260 240 220 200 180 160 140 120 100 80 60 40 20 ----------— +!-, / / I ' I II! I \ \, / ',x II I I I I 0 10 20 30 40 50 60 70 80 Angle of incidence (degrees) 90 Figure 2.14: Measured and theoretical A/ as a function of incidence sand depth angle for a 57 rmm

30 formula comes from Polder-Van Santen (1946): m = 6h + - C h)E[( - ) (2.17) where Au is the depolarization factor of an ellipse along the u axis (a, b and c are the three axes of the ellipse, u C {a, b, c}), cm is the effective (measured) dielectric constant, 6h and ci are the dielectric constants of the host and sand particles (inclusions) respectively, and the volume fraction, V, is the fraction of the sand volume to the total volume of the mixture (i.e. vi = Vsand/vtotal). For prolate spheroids, the depolarization factors are given by: AC = [n - ) 2e] and Aa = Ab = -(1 - AC). (2.18) 2e [l n(+ e ) 2 The eccentricity, e, as discussed in section 2.3.2, has a value of e 0.65. The mixing formula was derived by isolating an individual particle of the dielectric, calculating an effective field incident on the particle and then calculating a "reaction field" from the particle itself. The variable, *, is the effective dielectric constant of the medium immediately surrounding the particle which is representative of the apparent "reaction fields" from neighboring particles. We know that ~h < C* < (C, and usually e* = m is chosen for dense materials. While it is common practice to use ce = cm for dense mixtures (Ulaby et al., 1981; De Loor, 1973), it was found that c* = t = 3.6 (a value close to e = 3.78, the dielectric constant of a solid slab of silica) gave a much better agreement between measurement and theory. Using a measured volume fraction of 0.6 (ratio of sand volume to total volume) and * = cm gives cm = 2.4 while using c* = 3.6 gives the measured value of cm = 2.64. This suggests that for the very dense mixture of the sand, individual particles of silica have a dominant, if not complete, effect on the fields internal to the sand layer. Finally, as a last part of measurement analysis, it is of interest to quantize the amount of incoherent received power (due to volume scattering) to the total power received. Using (2.9) to find (Pinc) and (Pm) + (Pin,) to find the total power, (Ptot), the ratio of incoherent to total transmitted power can be found. This is tabulated and expressed in decibels in Table 2.1 for each sand thickness and incidence angle measured. The average power ratio of incoherent to total power is -15 dB with some ratios becoming as large as -8 dB. This is likely due to positioning differences of the radar between the calibration and measurement procedures and possible strong reflections from the table edges for certain incidence angles (this may be occurring at the incidence angle of 63 degrees). Using the phase difference function (2.6) should alleviate the radar positioning problem and using a least squares measure that is proportional to the variance of the received signal (2.15) allows for the

31 sand depth - 4mm 7mm 13mm 17mm 57mm polarization - vv hh vv hh vv hh vv hh vv hh incidence 20.5 -17.46 -17.98 -20.35 -20.48 -17.30 -17.43 -14.39 -13.84 -14.65 -14.98 angle 27.5 -17.69 -18.13 -21.65 -21.60 -16.90 -16.47 -15.91 -15.59 -12.95 -12.45 44.0 -15.12 -14.62 -20.90 -21.60 -15.32 -14.45 -17.50 -18.46 -13.94 -14.10 55.8 -13.77 -10.60 -19.26 -14.97 -16.56 -20.15 -16.20 -15.19 -13.36 -12.64 61.5 -14.47 -9.36 -19.00 -12.29 -15.62 -23.56 -16.43 -9.77 -12.64 -12.94 70.0 -14.43 -9.09 -18.45 -11.25 -12.09 -24.23 -16.16 -15.90 -8.33 -16.33 Table 2.1: Ratio of incoherent to total power in dB measurement uncertainty to be weighed into the fitting routine. Table 2.1 does however reflect that the lower angles of incidence introduce a greater amount of received incoherent power, as is expected because the transmitted signal must pass through more sand than it would for high angles of incidence and the projection of the antenna beam onto the sand surface is larger at low angles of incidence. 2.6 Conclusion In this chapter, a detailed set of measurements of very fine sand at 35 GIIHz is described. The purpose of this experiment was to characterize the effective permittivity (or equivalently the propagation constant) of the medium. Through this analysis it was learned that * A strong electromagnetic interaction exists between neighboring particles; a stronger interaction than that predicted by the standard form of the Polder-VanSanten mixing formula * The total power reflected from the sand surface is dominated by the coherent component. While an incoherent component is present, it is not sufficient to be detected by the presented method of determining effective permittivity. The execution of this experiment raised a number of interesting research questions, which have been the motivating theme for much of the work presented in this dissertation. That theme is: how do we relate the propagation of electromagnetic fields through a macroscopic medium to the physical properties of its microscopic constituents? This chapter has provided a physical feel for the problem as it exists when the observing wavelength and discontinuity dimension exist on similar scales. We now step into the realm of theory and numerical

32 simulation to further our understanding of this complicated, yet fundamental process in electromagnetics.

CHAPTER 3 Theoretical Methods 3.1 Introduction The preceding chapter described a careful set of measurements used to determine the effective permittivity of a dense, inhomogeneous medium. For natural media, the effective permittivity is the parameter that mathematically characterizes the manner in which coherent fields propagate through the medium. Often however, a radar observes a combination of coherent fields (scattered from a surface for instance) and incoherent fields (scattering from a volume). Of these two phenomena, it is the coherent fields that are the most fundamental, because it is the coherent field in a random medium that excites the incoherent field. If the coherent field is not understood accurately, then it stands that we cannot fully understand where the power in the incoherent field arises from. In this chapter, commonly used theoretical methods for modeling coherent and incoherent fields are presented. There are a number of components that should be considered first before delving into the formulation. To render a theory tractable, simplifications are often made to reduce the problem into a simpler form. Physical knowledge encompassing the following factors are an important first concern in developing and choosing the correct theoretical method: 1. density of the discontinuities (volume fraction) 2. dielectric contrast between the discontinuities and the background 3. dimension of the discontinuities relative to the wavelength 4. statistical positioning of the discontinuities with respect to one another (internal microstructure) 5. discontinuity shape and orientation (preferred or uniformly distributed) 33

34 Of these, the first three, volume fraction, contrast and dimension, have the greatest effect on electromagnetic fields. Low density, low contrast and small dimensioned discontinuities do not perturb the mean field appreciably, and thus, single interactions dominate much of the observed behavior of the medium. The development oftheory to describe the theory to describe the phenomena of electromagnetic wave propagation through these media is fairly straightforward. As these parameters change so that the interaction for a given discontinuity is strongly interacting with the incident field, the analysis becomes increasingly complex. Under these circumstances, it is necessary to account for more than just single particle interactions with the incident field. Often, neighboring discontinuities will form a single, statistically coherent scatterer, that becomes the dominant scattering element for the medium. The theoretical challenge in these situations is to understand and characterize this multiple discontinuity interaction with the propagating fields. In the following sections, the coherent field problem (i.e. determination of effective permittivity) is addressed first and then followed by theoretical methods of determining incoherent scattered power. For both the coherent and incoherent subjects, the simplest situation of low frequency, low density and low dielectric contrast is developed first before approaching the more rigorous, multiple interaction models that are commonly used for dense random media. 3.2 Coherent Fields 3.2.1 Mixing Formulas Mixing formulas relate the macroscopic propagation constant to the microscopic properties of the constituents. Most mixing formulas are inherently low frequency in that they do not account for losses due to scattering from individual discontinuities. Of the variety of mixing formulas that exist, the one given by Polder and VanSanten in 1946 is the most prevalent [Polder & VanSanten, 1946; de Loor 1956]. The framework of the Polder-Van Santen mixing formula is the same for both the twoand three-dimensional problems. In this context, a volume may refer to a two-dimensional area or the more traditional three-dimensional quantity. The derivation of the Polder-Van Santen mixing formula begins with an arbitrary volume where the average fields along the boundary are everywhere known (Figure 3.1). The average electric flux density within the

35 Figure 3.1: Electromagnetic model for Polder-Van Santen mixing formula volume may be written as (D)=eff (E)= + f()+( - -h) jf (E'n dv (3.1) where f is the volume fraction of inclusions whose permittivity is given by Ci, and 6h is the permittivity of the host material (nominally free space). The electric field within the inclusion at low frequencies is related to the incident field through the polarizability tensor, - int Gk, such that En = t E, and axx axy axz e aOyz Oyy a yz (3.2) zX O azy Ozz If the orientation of the inclusions is randomly distributed about all possible angles, then the average induced field within the inclusions may be written as (Eint)= () (E) (3.3) In three dimensions (c, takes on the form of (axX + ayy + azz)/3; for two dimensions it is either arzz for TM polarization or (aXX + ayy)/2 for TE polarization. In general, the polarizability tensor is a function of the inclusion shape and dielectric contrast. Two Dimensions For the special two-dimensional case of TM polarization, the electric field always couples into the i-direction, and the integral in (3.1) is only over the included volume, making the

36 polarizability tensor equal to unity. Thus, the mixing formula becomes a simple linear interpolation between the two extremes of zero and unity volume fraction: |eff = h + f(Ci - Eh). (3.4) For TE polarization, the electric field lies in the plane of the discontinuity and can have components in both the x- and a-directions. The polarizability tensor then does become a function of the shape and contrast of the discontinuity. This quantity may be determined analytically for circular inclusions and numerically for inclusions with other shapes [Sarabandi and Senior, 1990]. The diagonal elements of the tensor are written as shape Co C1 c2 1 1 (i +c 1 4 0 0 Cc0^) = where (3.5) 2 i + ~h Ci + c2 A 5.28 4.17 5.95 D 4.32 3.38 3.76 The mixing formula for TE polarization is 1 E6 - Eh Ei + C1 Ceff = Ch + + -- Co (3.6) 2 ec + 2 h Ci + C2 Three Dimensions In three dimensions, the diagonal elements of the polarizability tensor are given by 1uu (3.7) i+A e -1) where Au is the depolarization factor along the u-axis, which is typically chosen to lie along one of the principal axes of the discontinuity. For spherical particles, Au = 1/3; values of A, for ellipsoids may be found in [Ulaby 1986, pg. 2040] or in general by the elliptical integral given in [Polder and Van Santen 1946]. The resulting mixing formula, which is the classic result given by Polder and Van Santen, is 11 =eff = h + f (e- h) + (3.8) Variations In the formulas that have been derived so far, the dielectric contrast between the host and included materials plays the role of modifying the magnitude of the induced electric field (via the polarizability tensor). In reality, the mean field surrounding an inclusion may

37 E sca E inc B- M Figure 3.2: Single particle interaction model for the effective field approximation not be the mean field as it would appear in the host medium. Rather, it may be more accurate to model this field as if it existed in a medium with a dielectric constant of C* that may take on a range of values [De Loor 1973]. For sparse media, this dielectric constant would be that of the host material (i.e. E' = Ch). As the density increases, C* may be either eff or even c, for a very dense medium. To include this variation into the mixing formulas so far derived, the constant Ch in the polarizability expressions given by (3.5) and (3.7) is replaced by e*. In three dimensions, the result would be identical to (2.17). 3.2.2 Foldy's Approximation In the previous section, a mixing formula was derived based upon the average electric fields within a random medium. The electric fields outside of the discontinuities were that of the incident field; the fields within the discontinuities was determined through a lowfrequency polarizability tensor. In contrast, Foldy's approximation, or the effective field approximation (EFA), determines the effective permittivity via a dispersion relation dlerived from a single interaction of a discontinuity with the incident field (Figure 3.2). Thus, the effect of scattering is taken into account on a, first order basis. We begin the general treatment of this theory by assuming a uniform plane wave propagating in a medium with an effective propagation constant of, such that K -= W/btoeff and (E= Eoei' (3.9) In the above equation, () represents the fields outside of the discontinuities; because this theory only accounts for single interactions, this uniform plane wave is everywhere valid. As a consequence, the average internal fields for each discontinuity are also plane waves [Brown 1980], and the average field within the scatterers may be written as 2-f - _)( -int ( 2+k ( ) (3.10)

38 where (Eint) is the average field internal to the discontinuities, f is the volume fraction, and ci is the relative permittivity of the discontinuities (inclusions). The scattered fields are related to the internal fields of the scatterer by the free space dyadic Green's function, G and ( sca k(c - 1) (Eint) Gdv'. (3.11) Substituting in the expression from (3.10) above, we have - E sca) (-K2 + ) J (E) Gdv'. (3.12) Two Dimensions In two dimensions, the scalar for of the Green's function is [Tai, 1994] Go = - Ho (kop) a e-~ ' (3.13) 4 4 V/p k.(18 In the above equation, the far field approximation is included for the zeroth order Hankel function of the first kind. For a TM polarized field, the electric field is solely in the z-direction and (3.13) becomes 4ifTT = (-,2 + ko) f ciko kks.p'dvI, (3.14) where Ak = K - ko, k = k, Vp is the volume enclosed by a single discontinuity, and TTM is the average of the forward scattering amplitude of a single scatterer. Following [Brown 1980], the integrand is expanded into a power series and integrated term by term, the resulting integral is equal to the volume of the particle, and so we arrive at the sought after dispersion relation 4inoTTM = -_2 + kO. (3.15) where no = f/Vp is the number of discontinuities per unit volume. In terms of effective permittivity, (3.15) is i4noTTM eff 1 - (3.16) A similar derivation may be derived for the TE polarized, two-dimensional field. In this approach, it is easiest to solve for an Hz directed field, which would leave much of the

39 derivation unaltered save the expression for the forward scattered field TTE. In this case, the effective permittivity becomes i4nO TE eff- 1 TTE (3.17) For circular cylinders, the forward scattering amplitude may be determined analytically [Ruck 1970, pg. 205]. For the TM and TE polarizations, these expressions are 00 00 TTM = TM and TTE = TTE (3.18) m=0 m=0 TTM _ -kiJm(koa)Jm(kja) + koJ'(koa)Jm(kia) T^M m (3.19) kiH()(koa)J (kia) - koH(l1)(koa)Jm(kia) TE -ki/iJm(koa)J (kia) + ko/eoJm (koa)Jm(kia) In the above equations, Jm and H() are Bessel and Hankel functions of the mth order, the prime indicates a total derivative, and the factor 6m is equal to unity for m = 0 and two for 7n > 1. Three Dimensions In three dimensions we use the scalar free-space Green's function given by -ikok.(r-') Go = 47rR (3.21) in (3.12) where R = r- r'|. After performing similar mathematics as was done for the two dimensional case, the expression for the dispersion relation becomes [Brown, 1980] 47rnoF = -2 + k2 (3.22) which gives for the effective permittivity 47rn0 Ceff -1 + ko2F|. (3.23) Here, F is the forward scattering amplitude of the discontinuities. For spherical inclusions, the expression for F is derived from the Mie series given by [Ruck et al. 1970, pg. 141] F (i)n-l (n 1)(A - iB)] (3.24) n=1 l 2

40 E B sca E inc Figure 3.3: Two particle interaction model for the quasi-crystalline approximation and ( -) 2n + 1 j(koa) [kiajk(kia)]' - jn(kia) [koajn(koa)]' A )n(n + 1) jn(]kia) [koah(1o)(k' (ka)] [kiaj n(a)]' (3.25) = ( n+l 2n + 1 [j(koa) [kiajn(kia)]' - i jn(kia) [koaj'(koa)]' n(n ) h(1)(koa) [kiajn(kia)]' - -ijn(kia) [koah()(oa)] 3.2.3 Quasi-Crystalline Approximation In the previous section, an expression for the effective permittivity of a random medium was derived based on single scattering interactions with the incident field. As a result, reff becomes a linear expression related to the forward scattered field from a representative discontinuity, and the number density of the discontinuities. In this section, the next level of complexity is approached by taking into account all two-particle interactions (Figure 3.3). Note that this does not imply a "double-scattering" interaction as might be the case if we were to model an incident field upon one discontinuity, alone, and then to use this scattered field and the incident field to determine the scattered field from a second particle. Rather, a transition-matrix (T-Matrix) method is employed that accounts for all interactions between two particles simultaneously, thus determining the scattered field from the pair of scatterers. Because of the pair-wise interaction that is modeled by this theory, it becomes necessary to statistically describe the position of two particles with respect to one another in the random medium. This additional complexity will be addressed as a separate issue in a future section. The general and three-dimensional derivation of the quasi-crystalline approximation (QCA) follows closely that of [Tsang et al. 1986], whereas, the two-dimensional presentation of this theorem is unique to this dissertation and has been published in [Siqueira and

41 Sarabandi, 1996]. As has been the case in previous sections, the term "volume" may refer to a two-dimensional area, or a true, three-dimensional volume. We begin with the multiple scattering equation for plane wave incidence using the Tmatrix approach. Consider N particles whose centers are denoted by 4T (l - 1,..., N), randomly distributed over an volume V. Expanding the fields in terms of cylindrical or spherical basis functions, it can be shown that [Tsang, 1985; pg. 454] N N ( /) i. i = A 3(kT)T ) + zeikrljinc (3.27) j=l where W() is the vector of exciting field coefficients for the 1th particle, ainc is the vector of coefficients for the incident field, cr(krjrj) is the translation matrix from a coordinate system centered on the jth particle to one centered on the 1th particle, and T is the transition matrix which translates the exciting field of the jth particle into the scattered field from the jth particle. The incident field is written in terms of the basis functions as E int(r) a incRgtn (kr) = tainc. (3.28) n=O where, Rg4 represents a column vector of basis functions that is regular (finite) at the origin. The superscript t indicates a vector transpose, making the column vector a row vector in this instance. In general, the basis function will be designated as either '1 -for outgoing fields or RgI for standing waves and incident fields. By determining the expected value of (3.27) with respect to the 1th particle, we arrive at El [w-(] ] El [(kri-)F)E [w(J)]] + - iki-. (3.29) j3l Under the conditions that the number of particles are large and particle density is not very high, we make the approximation fundamental to QCA that Eji [w )] - Ej [W ] = w(() (3.30) which states that the expected value of the exciting field on the jth particle, which will excite the 1th particle, given the location of the jth and th particles, is equivalent to the expected value of the exciting field given only the jth particle. In other words, the set of exciting field for the particles which illuminate the 1th particle are not themselves dependent on the

42 positions of other particles in the medium. Thus, three particle interactions are excluded, and using the approximation in (3.30), (3.29) becomes w(rl) = (N - l) ] k(krlrj)T w(rj)p(rjrT) drj + eikain (3.31) If we normalize the conditional probability p(Tjlr1) to the volume, V, we have the pair distribution function, g(r) such that p(r = (3.32) which has the asymptotic property of approaching unity when Tj - ri is large. For a large system of particles, the fraction NV- is approximately the particle density, no. Letting V be the positive half-space xj > 0 [Tsang, 1985; pg. 491] w(Tl) = no I g(rjlrl)a(kTrlr)T()w(r(j) drj + e ikzanc (3.33) 37.>0 To evaluate the integral in (3.33) we use a trial solution for w(rT) such that it expresses a field traveling in a medium with effective constant of propagation, Kt w(Tl) - aEeiX. (3.34) After substituting (3.34) in (3.33) and making a simple coordinate transform (rjrl = r l = r), the integral in (3.33) becomes s eci't no ~ g(r)a(kr)eiKxd-rTaE = noeiir (j (kr)e dr + no 0 [g(r) - 1] W(kr)ei dri TaE = ei/i (1 + I2)TaE. (3.35) where V is the positive half-space excluding the volume occupied by the 1th particle. To evaluate II, we note that the vector addition matrix a(rjri), solves the wave equation, because it is composed of basis functions of order p for either the two- or three-dimensional problems. Here, rjrl is the translation vector between coordinate systems, 0j1 is the angle that the vector makes with the x-axis, and m and n (p = n - m) denote the harmonics in the jth and the th coordinate systems respectively. Additionally, a plane wave propagating in an effective medium with a wave-number of t satisfies the wave equation. Thus, we have V2'p +2 k2p - 0 (3.36) V2eix + e2 ei X 0 (3.37)

43 I'....-.:........S S Figure 3.4: Domain of integration for the quasi-crystalline approximation. Shown are the three surfaces (Sd, Se, and Sx,) and their inward pointing surface normals for the integral in (3.38). By multiplying (3.36) by e&' and (3.37) by ap, subtracting (3.37) from (3.36) and applying Green's second identity, I, becomes a surface/contour integral [L]nm -- k2 "2 Vo' - 'pVeI] dS (3.t38) as shown in Fig. 3.4. I1 may be written in terms of the three components of the surface where inward pointing normals were used for convenience. By employing the far-field condition, Ioo = 0, we are left with an integration over the exclusion volume, IC and the half space surface, Id. While both IC and Id are evaluated differently based on the dimensionality of the problem, it will be shown that Id has a ei(k-)x' dependence for both the twoand three-dimensional problems. For the second integral in (3.35), we have [i]nm - n > b - 1] H( )(kp)e-'PeKd- (3.40) x > -xi where we note that 12 is dependent on the propagation constant I, the pair distribution function g(r), the particle diameter, and the particle location x1. If we notice however

44 that [g(r) - 1] is nearly zero for Irl greater than a few particle diameters, we can treat the integration as if it were over an unbounded space, thus ignoring the boundary at x -= -xi (and particles near that boundary) [I2] no j j [g() - 1] eiPcOs dp dQ (3.41) nm Q b where Q represents integration over the two- or three-dimensional solid angles (note: Q represents the solid angle for a disc/sphere of radius b). Furthermore, if we assume that the pair distribution function is axially symmetric, (3.41) may be written as a single integral which can be calculated numerically. The exact form of this integral is unique to both the two- and three-dimensional problems, and is specified in (3.61) and (3.69). Depending on the particular problem, the appropriate form of I2 given by (3.40), (3.41), or (3.61,3.69) may be used. Referring back to (3.35), we can multiply I7 and I2 by eiNx to get a new expression s = [Sleixl + s2eik1] TaE (3.42) where the explicit dependence of (3.35) on the constant of propagation (k or K) is made clear (note: s2 is a function related only to Id from above). In (3.42) we have 1 = ^2 -k2I + I2 (3.43) and - ~-i(k-)x -no (3.44) SK2 - kI Using (3.42) in (3.33), the multiple scattering equation may be written aeei = [ 1eixl + 2eikxl] TaE + eikainc (3.45) If we balance the exponential terms of eikx and e&'X we then have two independent equations 0 = s2TaE + ainc (3.46) aE = TaE (3.47) which comprise the Ewald-Oseen extinction theorem (3.46) and the Lorentz-Lorenz law (3.47). To find the constant of propagation, we solve (3.47) by noticing that the determinant of the matrix Q = [iT - I] (3.48)

45 must be equal to zero for the non-trivial solution. Thus the solution for QCA rests in minimizing the determinant of a matrix whose elements are related to si and the single particle transition matrix, T. To implement QCA, the matrix Q is constructed from sI (3.43), which employs Ie and 12 which will be determined below. The effective constant of propagation of the random medium is min eff [detQ] (3.49) The user defined variables in (3.49) are: 1.) the pair distribution function, g(r) 2.) the particle density, no, 3.) the particle diameter, b, and 4.) the maximum order of the cylindrical/spherical wave expansion, pmax such that p C [O,pmax]. From a practical point of view, the pair distribution function is the most difficult of these variables to determine. Two Dimensions The basis functions for the field expansions in two dimensions are the cylindrical functions [Chew, Wang, and Gurel, 1992] n= Hn)(k>)e (3.50) Rgn = J(kr<)ein (3.51) where r> = larger(T', ) (3.52) r< = lesser(r, r') (3.53) when r is the point of observation and r' is the source point. The translation matrix, a, is given by the Hankel function addition theorem [Chew et al., 1992] [((krrl)l]nm = H1) l(krlrjr)ei(m-n)kJ, (3.54) where rjrJ is the translation vector between coordinate systems, jl is the angle that the vector makes with the x-axis, and m and n denote the harmonics in the jth and the!th coordinate systems respectively.

46 Explicitly, the integral IC in (3.39) and (3.43) may be written as [Ie] = 'I ei PCOSe-iP0H'(l)(kp)kpd0 nm J p - H(1)(kp)e-iP ( eiKPCOS p d kpH'(l)(kp)Yp - pH()(kp)-Y (3.55) where Yp = - iP eiKPcos' d: -T= 27riJp(Kp). (3.56) Substituting (3.56) into (3.55) and evaluating at p = b we have [Ie] = 27riP [kbH(1)(kb)Jp(b) - KbH(l)(kb)J;(lb)]. (3.57) The surface integral Id in (3.39) and (3.44) over the planer surface x =-xl is [Id] = / e aYp - Yp e dy \1~0 ~ X dy _ ei XP- X i^e X (3.58) oxwhere k1-2ik Xp = H(l)(kp)e-idy = -2e (3.59) -00 Substituting (3.59) into (3.58) and evaluating at x = -xi gives |[d =mm 2ie(!-t) [I + -| (3.60) and [72] =m no iP' [g(p) - 1] H(l)(kp)Jp(p)p dp. (3.61) In the general and two-dimensional version of the derivation, the field quantity being solved may be for either E, (TM polarization) or H, (TE polarization). Of the equations used to determine K, only the transition matrix in (3.48) is dependent on the field polarization. For a circularly symmetric cylinder, the T-matrix is diagonal and is given by (3.19), for TM polarization and (3.20) for the TE polarized fields.

47 Three Dimensions In three dimensions, the basis functions for the field expansions are the spherical functions [Wang and Chew 1993] = [= N] (3.6(2) RgT = [RgM RgN] (3.63) where 1 M = V x Tr and N =-V x M (3.64) ko, h) (kor)YIm(, ) and Rgilm = ji(kor)Ylm(O, ). (3.65) In the above equations, ji and h() are the spherical Bessel and Hankel functions of order, 1, and Yim is the associated Legendre polynomial, defined in this dissertation as Yi(, ~) ( () (211) (1- m)! m( )ei (3.66) 47r ( + m)! Pi(cs)e where P (cos S) are the ordinary Legendre polynomials given in [Abramowitz and Stegun 1964]. It is necessary to exactly specify which form of the associated Legendre polynomial is being used because the convention varies from application to application (see [C(hew 1990, pg. 395] vs. [Wang and Chew 1993]). The translation matrix formulas for CT are found in [Tsang et al. 1986, pg. 448]. Using these equations, the exclusion integral I, in (3.39) and (3.43) is e = 4iPb [kbh(1)(kb)jp(b) - Kbhlp)(kb)j(I(b)] (3.67) Additionally, the surface integral Id in (3.39) and (3.44) over the planer surface x =-xl is Id = -27rip+ei(k) [i + ]. (3.68) The three dimensional form of the exclusion integral given by (3.41) is [I2] no / (- i [g(T) - 1] hp(kr)jp(Kr)rdr. (3.69) For a medium composed of dielectric spheres, the T-matrix which is used in (3.48) is For a medium composed of dielectric spheres, the T-matrix which is used in (3.48) is

48 diagonal and is given by An 0 (3.70) 0 Bn where An and Bn are components of the Mie series expansion for a dielectric sphere. These constants are the same ones used in the previously derived three-dimensional form of Foldy's approximation (and, incidentally, these are the same coefficients that will be used in the T-matrix electromagnetic solution technique), and are given in (3.25) and (3.26). 3.2.4 Variations The theory that has been presented in two and three dimensions is termed the quasicrystalline approximation. It is termed as such because it assumes that a pair distribution function exists that statistically describes the positions of discontinuities within a random medium, but that at distances greater than a few particle diameters, the medium is essentially orderless. This physical representation of the medium implies that the electromagnetic field interactaction between two particles takes place through free space. This interaction is formulated through the use of a free space Green's function. If however, we say that the field is traveling through a medium whose permittivity is not unity (i.e. free space) but has a permittivity that is equal to that of the effective permittivity of the medium, the theory is termed the quasi-crystalline approximation with coherent potential (abbreviated QCA-CP). Details of this theory are given in [Tsang et al. 1985]. The effect however is essentially to replace the argument of the cylindrical and spherical wave functions discussed above with the propagation constant, for the effective medium. This approach has the added advantage over QCA of being self-consistent (i.e. it is expected that the behavior towards unity volume fraction would give an effective permittivity equivalent to that of the inclusions). 3.2.5 Pair Distribution Function The pair distribuution fnction arises from the attempt of QCA to model two-particle interactions within a random medium. The role of the pair distribution function is to

49 provide a statistical description of the position of one scatter with respect to another. This is a problem well studied in the statistical theory of liquids, where molecules display a quasi-crystalline order on a short range scale, and is disordered over long ranges [Fisher 1961]. Formulation The probability density function which describes the likelihood of possible particle locations (or translational states) of a set of s particles can be written as p(li2i...,I-s ) (3.71) which is the probability of finding s particles within small regions, dr, centered around the locations ri, r2,.o. The particle distribution function (or the somewhat misnamed correlation function), g(...), is defined such that (nr2...,r) = 1 (n 2...,) (3.72) If the particle locations may be located anywhere within the bounding volume, V, independent of one another, the probability of finding a particle at any one given location is p(r) = I/V and thus ne-7 p(rlr2,...,)= p(Yr)= vs (3.73) 1<i<s which means g(T1, r 2,..., 7) = 1. For a pair of particles, the probability of finding a particle at T2, given the location of a particle at r1 is p(r2K1) - p(2, 1) 1 =() - r(3.74) which may be considered as the definition of the pair distribution function used in QCA. In the absence of external forces, it is expected that the pair distribution function is radially symmetric. In this case, the probability density of finding a particle at a given distance (p or r) from another particle is determined by integrating (3.74) over 0 and A, giving P(p) g(p) -~ for two dimensions (3.75) p(r) g(r) V in three dimensions. (3.76) Ideal Gas Pair Distribution Function In the limiting case of an ideal gas composed of hard spheres, the pair distribution function takes on a simple form. For a sparse distribution of spherical particles, given the

50 f f= 0.50 ~ - - f = 0.40 | ~ -' f= 0.30 23 -.. f= 0.20 Radial distance, r/b Figure 3.5: Percus-Yevick radial distribution function for hard spheres. Horizontal axis represents radial axis normalized to the sphere diameter, b. position of one particle, another particle may be found at any distance greater than one particle diameter. This pair distribution function is referred to the hole-correction formula and is written as 0 r<2a g(r)- = (3.77) 1 r0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 The use of this pair distribution function in the quasi-crystalline approximation has the effect of eliminating the I2 term in (3.43) and thus the propagation constant is dependent only upon the exclusion integral given in (3.57) or (3.67). Percus-Yevick Pair Distribution Function For a medium composed of impenetrable particles, the presence of one particle precludes the presence of another. In the case of an ideal gas, this effect was taken into account only on a first order basis, through the hole correction formula. As density increases however, this process of exclusion becomes iterative, making the medium appear crystal-like on a scale of a few particle diameters. The most accepted model of this phenomena is given by [Percus and Yevick 1958], which was solved numerically by [Wertheim 1963] in three dimensions written afor two dimensions. Often the pair distribution function is thought of as being a robust quantity and the pair distribution function of one medium is applied to another with the hopes that the function distribution function of one medium is applied to another with the hopes that the function

51 would still be applicable. For instance, it is very common to apply the Percus-Yevick pair distribution function for classical fluids to the problem of snow [Wen et al., 1990; Tsang and Kong, 1992; Tsang, 1992]. In fact, the pair distribution is very dependent on the method of particle arrangement as will be illustrated by example here and in the next chapter. One interesting aspect of the pair distribution is that it remains unchanged under the process of particle extraction (assuming that particles are not allowed to rearrange themselves). The example given here highlights how similar arrangements of particles may result in strikingly different pair distribution functions. In general, the pair distribution function may be written as [McQuarrie, 1976, pg. 258] g(r2, T ) = P(72 r1 )/p(r1), (3.78) where p(r2\r) is the probability of a particle being located at T2 given a particle located at r1, and p(rI) is equal to the particle number density, no. Note that when 1r2 - rl I> b (where b is the average grain diameter), p(r2lrl) = p(r1), and then g(72, 1) =1, as expected. Now, given a particular particle arrangement method (say, a classical fluid), we have the probability function po(r2IT1). (3.79) If we reduce the particle density in our system by the factor f, we have P fo(1) - fpo(i1) (3.80) and similarly (because particle removal is independent of particle position), Po(r21) = fpo(72\1). (3.81) Then by (3.78), the pair distribution function remains unchanged. This result can be easily demonstrated numerically by generating random arrangements of particles. A common method of simulating a classical fluid is to randomly introduce particles into an empty container (Figure 3.6). In Figure 3.6a, the particles represented by dotted lines illustrate a classical fluid simulation for a volume fraction of 30%. The particles represented by solid lines in Figure 3.6a show the remaining particles after reducing the number of particles from a volume fraction of 30% to a volume fraction of 10% by random extraction. Figure 3.6b gives a classical fluid distribution of particles for a volume fraction of 10% where, as for the 30% case, particles have been randomly placed in the box but not removed. The pair distribution function for these different methods is given in Figure 3.7.

52 100. I,,. ' 'N....o. O.:: 0.....: 100 50 50.. U O 0 oCP O0 0 To o.0 0 0( '""' 0 r) -50;.;... " C -100 — -100 -50 0 50 100 (a) U1 u0 D 0 ( -5 0 0 0 -100 --100 -50 0 50 100 (b) Figure 3.6: Particle arrangement method for a classical fluid. (a) 30% volume fraction using sequential addition (dotted circles) and 10% volume fraction resulting from particle extraction (solid circles) (b) 10% volume fraction using sequential addition. 1.5 03 -Q CC. 0L 0. 0.5 1 1.5 2 2.5 3 Radial Distance (particle diameters) Figure 3.7: Comparison of the different pair distribution functions obtained from simulations demonstrated in Fig. 3.6.

53 radar backscatter specular medium 0 Figure 3.8: Volume scattering model. Observed power in the backscatter direction arises from the incoherent field scattered from dielectric discontinuities within the volume. Clearly, the pair distribution of the 10% volume fraction obtained from extraction yields pair distribution functions are distinctly different from the function obtained for the 10% classical fluid simulation. The reason that the two simulations for 10% volume fraction differ is that in the particle extraction case, particles are randomly removed, regardless of their position. When creating the 30% volume fraction simulation (from which the extraction example was derived), particles are occasionally excluded if they overlap, and it is this exclusion process which makes the classical fluid pair distribution functions different for the 30% and the 10% volume fractions. Reducing the number of particles does not alter the effect of this exclusion process. Thus, the fact that the pair distribution function does not change when particles are extracted has been shown both theoretically and by example. 3.3 Incoherent Fields Incoherent methods refer to electromagnetic analysis techniques that are used to characterize the incoherent nature of volume scattering. Incoherent energy is the component of volume scattering that a radar actually observes. Specifically, this is apparent for a flat interface where only the specular direction is dominated by the coherent field. In all other directions of observation, the radar signal observed is the incoherent power. This incoherent power arises from dielectric discontinuities that lie below the surface, the mean scattered electric field being zero.

54 3.3.1 Born Approximation Consider Figure 3.8 where a coherent incident field is incident upon a dielectric half space with permittivity fluctuations specified by E1(r). We separate the fluctuating permittivity and wave number into a coherent component and a fluctuating component as (lf(T) - 1() - Elm (3.82) kf(r) - k(r- (3.83) Often, im =-< E1(r) >, but the apparent permittivity may also be rlm = Eeff where ceff is derived from a mixing formula or comparable method from the previous section (use of ceff derived by this manner is referred to as the distorted Born approximation). The wave equation within medium 1 may be written as V x V x E1-2E1 = kf2(r)El (3.84) where k 2(Y)El behaves as a source term for the electric field. The solution to (3.84) may be written as E = E) + Go (rl). Elklf(r)d (3.85) where E is the electric field observed in medium zero, E is the incident field and the integral term represents the scattered field. The scattered field is composed of a source term, k2f(r)El, and the dyadic Green's function for a single interface with a source in medium one and observation point in medium zero, G01 (rr7). At this stage in the formulation, the equations that have been discussed so far are exact. The difficulty in solving these equations arises from two sources: (i) the effective propagation constant, ceff in the random medium is not known, and (ii) the electric field in medium one is unknown. The Born series (from which the Born approximation arises) circumvents the latter of these two problems by approximating the electric field in the random medium by the highest order known electric field. For the first term of the series, this is the zeroth order electric field induced by the incident field crossing interface between the two media, E)(Tl). Thus, the observed scattered field in medium zero, due to the zeroth order electric field in medium one, is E(1)Ol (,) (f ) lf()d-0. (3.86) In the above, the far-field dyadic, single interface, Green's function may be written as Go) = KICo - ik1.r kzs e( )l(klOzs)(] ) + Trl1 h(kozs)hl(klz )] (3.87)

55 where Srkr 2-dimensions C: - v rkr (3.88) 3-dimensions Here, rl and rl1~ are the transmission coefficients across the boundary from medium one to medium zero. These are specified by k01 = 2kl kz rls perpendicular polarization (3.89) kOz + klz kOz 0 _ 2Elkoz k2 koz 10 01 I parallel polarization. (3.90) 711s -ok lz+ koz - k klIz Is The incident electric field in medium one can be written in general as, (r) = EO eikll(-kli) + E011 1eikrhl (-k1z). (3.91) Equations (3.87) through (3.91) utilize a coordinate system that is referenced to the direction of propagation, k. Following [Tsang et al., 1986], e is related to perpendicular polarization (electric field is perpendicular to the plane of incidence) and h is related to parallel polarization (magnetic field is perpendicular to the plane of incidence). These unit vectors can be specified as follows: e(kz) = k- ykx) (3.92) kp h(kz) (Ak, +!ky) + (3.93) where kx = k cos sin 0 kz = k cos ( (3.94) ky = k s sin sin 0 k = k cos 0. Furthermore, for the calculation of (3.86), it will be necessary to calculate the dot product of various combinations of e and h. The four possible combinations are el (lzs) el (-zi) = cos s (3.95) e(I(klzs). hi(-k1i) = sin (3.96) k1 klzs hi(klzs). l(-kz) = sin (3.97) kl 1 /l(k~klzs) klk-l cosq) = (3.98) /~ ~iCs- kzkz 1 s where, i = 0 was used for the reference direction of the incident field.

56 To find the average scattered power from the dielectric discontinuities in the random medium (medium one), we multiply E(l)() by its conjugate and determine the ensemble average, Jr 2 _(2 2 ( )2) = (r ) O) (l) (G (i, T2) E (2)) (klf()kl ) r2)) drldr2. (3.99) For a perpendicularly polarized incident field, the perpendicular component of the scattered field is (M(_2) 12 E0_ [2 ] 10 10 12 2 J2 21 i- - (ll)l 2) E= loil2 ITTi cos2 / i(k 8)-i(k ) (f(rl)klf(r2)) drldr2 (3.100) The integral term in the above equation is common to all combinations of incident and received polarizations and can be written as B J ei(k s)*i(kk)2 kf ( )2) ) drl d2, (3.101) where A is the area illuminated by the radar. By first noticing that for a stationary medium, Kk f(Yl)k(r)) is related to the correlation function of the permittivity fluctuations as C( r - 2) = (t f(r )), (3.102) we can utilize the power spectral density, X to simplify (3.101). The power spectral density is written as k-r (k) - (2) fC(rl - r2)eik'dnr (3.103) where n refers to the dimensionality of the problem (i.e. n=2 or 3). Substituting for (kif(T1 )k2f(72)) in (3.101), it can be shown that 8'2k 4A(kx- - kx -k - kz ) B r z-zs ZI ) (3.104) kzi -'1 k\zs for two dimensions, and for three dimensions, 47r3ko4AO(kxi- kcsI ky - kys -k zizs) (- k. 1) B ' 1 -.k 1k- (3.105) kIz + lzs By substituting (3.101) in (3.100), an expression arises for the co-polarized perpendicular scattered field for the random medium in both two and three dimensions. The bistatic

57 coefficient, 7v, is determined by dividing (3.100) by the incident power density, Pi, which is given by f 27rr 2-D P = IEol2AcosOi/S S = (3.106) 47rr2 3-D Combining the general forms of (3.100), (3.101) and (3.95), the bistatic coefficient can be written for transmit polarization, t and receive polarization, r being one of v or h for vertical and horizontal polarization respectively. + C ( 27r2 k (2-D) (3.07) 7rt(e ) I ) r (3-D) I701 T01 12i1oz2 k cos2 COS (hh) I S I i klzs kzt s Iso r1 O12I koz k1i 2 sin2, (hv),01 kl k sin2 (vh) Is 1Zi k 101 01121 (k pips - kIziklzs cos,)12 (vv) Backscatter Radar Cross Section The backscattering radar cross section is determined by setting 0s = 0 in (3.107). In this instance, there are no cross-polarized terms, and we have 23 \ 0/ 1 0 114 k| z4 14 2r = 0 1co h (Trt 2 -osk or or (3.108) 7r24 I |0114 \ 0 151fIl / Correlation Function The backscatter RCS determined via the Born approximation is dependent upon the angle of incidence, the effective permittivity, and the correlation function of the permittivity fluctuations (via the power spectral density). By the power spectral density, the observed backscattered power as a function of observation angle and frequency can be regarded as a method of probing the internal microstructure of the random medium. For a given frequency, (3.108) samples a elliptical segment (semi-axes of 2k1 and 2To) in r- space, where X refers to the power spectral density function given in (3.103). For given angles of incidence and varying frequency, we sample along radial lines from the origin; at normal incidence (0 = 0~), the power spectral density is sampled along the kx = 0 axis, for grazing incidence (0i = 90~), the power spectral density is sampled along the diagonal defined by a line tan- ) - from the negative kx-axis (Figure 3.9). This relationship will be illustrated

58 k~i k x l 0 =90 0=0 Figure 3.9: Elliptical segment traces (0(0j)) of the power spectral density sampled by observing the backscatter RCS of a random medium. Shown are a number of traces for different frequencies. by example in the following chapter, where it will be shown that this relationship may be used to detect the mean grain size. There exists a specific class of functions that may qualify as valid correlation functions. A necessary condition for these functions is that the power spectral density, b must be positive for all values of k, and kF. Among the functions that meet this criterion are the radially symmetric exponential and Gaussian functions. While the Gaussian correlation function can be considered heuristic in nature, Debye [1957] derived a canonical exponential form of the correlation function based on the assumption that discontinuities are arranged randomly in space. Given a set of experimental data, it is common to fit an exponential function to the sampled correlation function [Vallese and Kong, 1981]. This approach however would occlude the asymmetry that might be expected between the k, and k, directions due to external effects such as gravity. Additionally, for a semi-crystalline structure, it is expected that the presence of a particle at one location makes it likely that there will not be a particle within a short distance of that particle. This implies a negative correlation, a factor that neither the Gaussian nor exponential functions are capable of modeling. Both of these effects are visible in the snow data from [Vallese and Kong, 1981], the essential results of which are summarized if Figure 3.10.

59 ' - -, 4-' (a) (b).... '............., 1:,..... I'/I |-................................ 2 42 5 8 14 S 8 0 4.^, 8 4 i mm 'r~m (c) (d) Figure 3.10: Sample data and correlation functions of snow [from Vallese and Kong, 1981]. (a) sample cross section of snow (dark regions represent ice, light regions represent air), (b) two-dimensional correlation function of snow sample, (c) horizontal cut of the twodimensional correlation function (shown is the horizontal correlation length used to fit an exponential function, (d) vertical cut of the two-dimensional correlation function (shown is the vertical correlation length used to fit an exponential function).

60 1+dl ds I 97 \' Figure 3.11: An elemental volume in radiative transfer. Scattering of incident power is accounted for by an extinction coefficient, is, and a source function,.F. 3.3. 2 Radiative Transfer Radiative transfer refers to a well known theory of energy transport proposed by [Chandreshakar 1960]. Instead of starting with Maxwell's equations as do the analytic models, radiative transfer relies on the conservation of energy and the incoherent addition of scattered power performed through an integral operation. Thus radiative transfer avoids the mathematical complexity of the analytic models which makes them difficult to use and impractical to extend to higher orders. Solution to the integral equation provides the radiated incoherent intensity as a function of position of the observer and source. The incoherent addition for the scattered energy does limit the density for which radiative transfer is applicable because of the lack of a method for determining the extinction coefficient and difficulty in defining an elemental scatterer (which may consist of one or more particles). As mentioned previously radiative transfer provides the basis for dense media radiative wave theory which has been made to account for phase coherence between scatterers using a the analytic quasi- crystalline approximation. The quasi-crystalline approximation accounts for phase coherence between scatterers on the microscopic scale and radiative transfer is used to model the behavior of energy transfer on the macroscopic scale. Radiative transfer then is presented here in its classic form which can be later modified to account for phase coherence in a dense media. The theory begins by analyzing a change in field intensity, represented by a four corn

61 ponent Stokes vector, I, as it passes through an elemental volume, ds (Figure 3.11): di - |ds | (3.109) ds The extinction matrix, K, consists of absorptive losses, Ka, and scattering losses, Ks. The scattering losses are due to the scattering of incident intensity from the direction s into other directions. These scattered fields can act as sources in other elemental volumes and is accounted for in the source function, f which is given by f - P(, S ')(r, s )dQ'. (3.110) The phase function P(s, S') accounts for the scattering of energy from the A' direction into the s direction by an elemental volume. Solution to the radiative transfer equation (3.109) for a layered media is accomplished by splitting the intensity vector into an upward and downward directed intensity vectors. In this case the radiative transfer integro-differential equation becomes the set of coupled integro-differential equations of the following form di n ds+ ds HI+(~,,z)+~f+ (3.1]1) dI __,-(- z) + -)+n- (3.112) ds L where [ = cos 0 (,u is positive for upward traveling energy and negative for downward traveling). The new source functions, F+ and F- account for scattering into the upward and downward directions respectively. By imposing boundary conditions for the energy at the upper and lower interfaces of the layered medium, a formal solution for the intensity components can be obtained. The solution depends on the unknowns of the extinction matrix, aK, and the phase function, P. The extinction matrix can be related to the effective propagation constant in the scattering media and the phase function can be related to the basic scattering element (either a single scatterer or a collection of scatterers) in the medium. For sparse media, the Foldy approximation and single scatterer phase functions can be used for JK and P respectively. For dense media, the quasi-crystalline approximation is used in conjunction with the Percus-Yevick pair distribution function to determine the extinction matrix, and either a T-matrix approach or the distorted Born approximation is typically used to determine the phase function. In all instances it is necessary to reduce theoretical complexity by assuming that the media consists of spherical particles. Alternatively, both

62 the extinction matrix and the phase fnnction can be determined either experimentally [Thesis proposal by John Kendra, 1993] or by numerical methods. It is these numerical methods that are a distinct part of this proposal and will be developed in the following sections. The pnrpose of the numerical methods however is not to fulfill the unknown gaps of radiative wave theory, rather they are presented here as an alternative formulation to the volume scattering problem. Intermediate results from the numerical methods may be used as tools for providing the unknowns in the theoretical methods discussed so far or alternatively as a method in its own right.

CHAPTER 4 Packing Algorithm for Determining Particle Arrangements The study of electromagnetic wave interaction with a collection of random particles (volume scattering) is of importance because of its application to a variety of radar remote sensing problems. Volume scattering theories are developed to determine basic electromagnetic properties of the medium such as the effective propagation constant, the attenuation constant, and the incoherent scattered power. Modeling efforts for volume scattering can be categorized into two groups: incoherent approaches and coherent approaches. In incoherent volume scattering theories, such as radiative transfer [Chandrasekhar, 1960], the effect of the phases of the fields scattered between neighboring particles is ignored. These methods are usually applied to sparse media where single scattering properties of constituent particles are used to formulate the volume scattering problem. Intensity approaches implicitly make the sparse medium assumption that individual scatterers are randomly positioned with respect to one another. As the scatterer density increases, multiple scattering between particles becomes, significant and the scattering solution based on the incoherent approach can become prohibitively complex [Tsang and Ishimaru, 1987]. Coherent approaches, on the other hand, such as the Born approximation and the quasi-crystalline approximation (QCA) [Tsang et al., 1985], account for the interaction between particles through the inclusion of a permittivity fluctuation correlation function or a scatterer center pair distribution function, both of which provide statistical descriptions of the location of scatterers with respect to each other. The need for the correlation or pair distribution functions adds another cornplexity to volume scattering theories, namely that associated with how to determine these functions for the medium under consideration. These functions play an important role in determining the scattering behavior of the random medium and thus must be characterized accurately. The importance of this characterization has been somewhat overlooked in the literature. Simple Gaussian and exponential functions are usually used as an approximation for the correlation function since they are amenable for algebraic manipulation, but 63

64 there is little evidence in support of the hypothesis that these are accurate representations of natural media. It is the purpose of this chapter to (i) present a numerical technique to determine the correlation and pair distribution functions accurately and efficiently based on physical modeling of particle arrangements in a random medium, (ii) demonstrate the application of the correlation function to the two-dimensional Born approximation and make a comparison with theoretically derived results, (iii) to use the derived correlation function to demonstrate limitations in the use of the exponential correlation function, and (iv) to highlight the use of packing algorithms as a method of scattering analysis that can be used to analyze or enhance existing scattering theories. The methods of determining the correlation function or the pair distribution function can be categorized into experimental theoretical, and numerical approaches. Experimental determination involves capturing an undisturbed sample of the volume under study and analyzing it for the desired information. This approach is quite difficult, very time consuming, and its accuracy depends on the particle size and measurement method. A classical example of such a process is given by [Vallese and Kong, 1981] where a layer of snow was infused with liquid plastic. The final result was the two-dimensional correlation function of the permittivity fluctuations due to the two dissimilar dielectrics of air and ice (see Figure 3.10). The observed correlation function for one sample was then fit to a theoretical model (exponential) from which further calculations could be carried out. Further assumptions, such as azimuthal symmetry and validity of the theoretical model, were required to estimate the three-dimensional correlation function. Theoretical methods of determining the correlation and pair distribution functions, while benefiting from a greater generality than experimental methods, must use simplifying assumptions to make the theories tractable and easy to handle. For the correlation function, it is common to use an exponential function [Debye et al., 1957], which is derived by assuming that individual particles are positioned randomly with respect to one another. Such an assumption is valid in the limit when particle sizes are distributed over a wide range and/or when the particle density is low. This function will be the used as a basis of comparison further in the paper. The pair distribution function of [Percus and Yevick, 1958] is another such theoretically derived function commonly used in conjunction with QCA and dense media radiative transfer [Tsang and Ishimaru, 1987]. This particular distribution is tailored for a medium consisting of discrete-sized spherical particles and in which there are no interparticle forces except exclusion, as in a classical fluid; external forces, such as gravity, are not taken into

65 account. For macroscopic media such as snow, soils, and sand, the Percus-Yevick distribution may not be applicable because particle positions are dominated by forces such as friction, gravity, and inertia. It is not difficult to come up with examples where the Percus-Yevick pair distribution function would give considerably different results than would be expected for macroscopic particles under the influence of these conditions. For instance, if we were to create a dense assembly of uniform spheres in their lowest energy state, a crystal, and then randomly remove 10-20% of them, the pair distribution function would not change because the basic structure of the crystal still remains (see Section 3.2.4). In a classical fluid however, particles are allowed to rearrange themselves to reach the lowest energy state, which is dependent on their relation to one another. Thus the pair distribution function in a classical fluid would change as the Percus-Yevick distribution predicts, but there is no reason to expect that the pair distribution functions of macroscopic particles under a different set of forces would be similar. Thus it is important to explore other avenues of determining these functions to either verify the validity of employing these theoretical techniques to specific applications or generate the functions themselves. Numerical methods refer to algorithms for determining the correlation function and pair distribution functions through the use of a computer to simulate the scatterer positions in a given space. Some of these methods are used to validate already existing theoretical distributions [Ding et al., 1992; Broyles et al., 1962], while others are used to determine the desired correlation or pair distribution functions directly [Buchalter and Bradley, 1992, 1994]. The class of packing algorithms can be split into two groups, those that physically model the pouring deposition of particles and those that generate possible particle arrangements. Of the physical models of pouring, a general trend is to model the interaction of hard grains as they "pour" downward into a given volume. Models such as this are used to investigate local particle ordering and the effects of shaking packed sets of particles. The algorithm of Buchalter and Bradley is the most elaborate of these methods, it models collisions and rotations of uniform two-dimensional [Buchalter and Bradley, 1992] or threedimensional [Buchalter and Bradley, 1994] monosized ellipsoids as they pour into a twoor three-dimensional volume. Nevertheless, modeling of the physical pouring of particles is complex and time consuming; as this set of models more closely attempts to match the physical processes involved in the pouring process, the computation time will increase proportionally. For this reason, the numerical pouring model of physical deposition may not be appropriate for Monte Carlo simulations and may not be amenable to the real world

66 situation of varying particle shape and size. For the scattering problem, we may be more interested in computational efficiency in determining possible arrangements of particles rather than physically modeling their deposition if the resulting arrangement appears reasonable. In the arrangement set of algorithms, some methods, such as in the sequential addition method [Ding et al., 1992], simply find open spaces for particles to fit within a confined area or incorporate an external force such as gravity, from which local potential energy can be minimized and particle stability can be maximized [Visscher and Bolsterli, 1972]. While the existing algorithms of this type can be considered to be much more computationally efficient than the physical deposition models, up until now they have been lacking in flexibility to model arrangements of particles of different shapes, particularly when it comes to the three-dimensional problem. The method described here should be considered a method of generating particle arrangements rather than a numerical pouring model. It takes into account the three important factors of determining particle arrangements: (i) nonpenetration of neighboring particles, (ii) gravity, and (iii) particle shape. This presented method is capable of determining particle distributions in both two and three dimensions and is computationally efficient enough to provide a large sample size for Monte Carlo simulations. The new method is not restricted to spherical/elliptically shaped particles and, as will be described, may be used in conjunction with coherent volume scattering theories to explore the scattering from dense random media. This chapter is arranged as follows: first the packing algorithm is described in Section 4.1, and examples for both the two- and three-dimensional problems are given. Methods describing accurate calculation of the correlation function and pair distribution function from the computer-generated particle positions are given in Section 4.2. In order to demonstrate the importance of the packing algorithm to the problem of volume scattering from a very dense medium, scattering solutions based on the Born approximation using the correlation function computed from the packing algorithm and an exponential function having the same correlation length are compared as a function of incidence angle and frequency. Because the angular backscatter response is proportional to the two-dimensional Fourier transform of the correlation function, it is possible to introduce the concept of a "visible region" where we can gain insight of how the choice of observation frequency samples the Fourier space of the correlation function. This is done in the second half of section 3, where the angular dependence of the backscatter response is demonstrated for a particular example at two different frequencies, one below and one above the resonance of the mean particle size. In

67 the Appendix, we offer a brief discussion of the dependence of the pair distribution function on the particle arrangement method. 4.1 Packing Algorithm We begin by introducing a user-defined distribution of particle states, Pp, where the state of a particle is a vector representing the particle size, shape, orientation, and dielectric constant. For the purposes of demonstration, elliptically shaped particles will be used, however, this packing algorithm is not restricted to the simplified elliptical shapes. In two dimensions the state of such a particle is represented by Pp = (ap, bp, p, ~p) (4.1) where ap, bp are the lengths of the two principle axes, 0p describes the orientation, and Ep is the dielectric constant of the particle. To simulate the particle positions in two dimensions, a rectangular region with j being the discrete horizontal coordinate and z being the vertical coordinate is considered (Figure 4.1). Coordinates of the intersection of the principle axes of the particle p denote the position of the particle, which, together with the particle state, completely specifies the particle in the medium. The vector state of a particle is chosen by a random number generator with a prescribed distribution for the different states ap,bp, (p, and cp. In this paper, it is assumed that the principle axes of the particle have a normal distribution (N(a, Ta)) and the orientation angle of particles has a uniform distribution (Uj(-7r,7r)). The exact distributions of the particle states may be chosen according to physical measurements or an empirical model. The particle is then laid down individually via the packing algorithm into the region to be filled. By making the limits of the region large compared to the sizes of the individual particles and by using the periodic boundary conditions, a semi-infinite layer of particles can be simulated. Periodic boundary conditions refer to making the spillover of a particle on one boundary appear on the opposite boundary of the region. The packing algorithm starts by setting down an initial layer on the bottom of the rectangular region so that the lower surface is covered by a single layer of particles. This set of particles creates a bumpy solid surface, Sp. Particles are then sequentially added to the region in the following manner: 1. A particle state Pp is selected from the distribution, thus setting the particle shape and orientation.

68 new particle, p oko < ~j Figure 4.1: Illustration of the two-dimensional packing algorithm. Particles are fit to the surface Sp to find the minimum height, z. Once a particle has been fit to the surface, it becomes part of the surface for the next iteration of the process. 2. The surface of the particle is discretized and then categorized into a set of lower and upper surface points, as shown in Figure 4.2. The points of the lower surface will be used in matching the surface of the particle to the current surface, Sp. 3. A fitting method is used to find the lowest height, zp, as a function of the horizontal coordinates, that the particle will fit to the surface. This is equivalent to minimizing the potential energy of the particle. 4. The particle's position, 7p, and state, Pp, are stored in a file, and the surface is updated to Sp+, using the upper surface points determined in step two. 5. The first step is returned to until the desired layer thickness is reached. The algorithm in its current state assumes that, averaged over a large sample, orientation is uniformly distributed between 0 and 2wr. An additional complexity may be added to the algorithm to minimize zp as a function of particle orientation also (i.e. 0 is no longer one of the preset states given in (1)). While this may provide a more realistic variation to the algorithm, it is unlikely that the additional computational costs would outweigh the benefits. Biases in the orientation distribution that the user wishes to impose can be more efficiently introduced via the probability distribution chosen for the orientation states. The fitting method described in step three is one of the key components that makes this algorithm computationally efficient. While any fitting method will be sufficient in this

69 I I I I I I I I +m -m i +m Figure 4.2: Illustration of the lower surface of a particle, Mlower. The particle is discretized into a horizontal grid and the maximum and minimum heights of the particle at points on the grid define the upper and lower surfaces. treatment, the one used here borrows a concept from the image processing technique of gray scale morphology [Serra, 1982; Giardina and Dougherty, 1988; Appleton et al., 1993]. In this packing algorithm, a gray scale dilation is used to find the minimum height that a given shape fits to the surface. To describe this procedure, consider a particle such as the one shown in Figure 4.2. The surface of a particle is represented by an upper (MUPPCr) and a lower surface (Mlower). The lower part of the surface determines how the particle will fit the the current surface, Sp, and the upper part determines how the particle will contribute to the new surface, 5p+,, once the particle has been deposited. Given a parametric equation for the surface of a particle, zp(i), these surfaces can be described by MuPPer(i) {max, z(i): -m < i < m} (4.2) Mw (ri) = {mmin zp(i): -a < i < m} (4.3) where for a normally oriented ellipse, zp(i) ~bp I- -), with ap and bp given in (4.1) and m being the horizontal limit of the extent of the ellipse (m > max(ap,bp) Vp). Note that the form of zp(i) is used here only for demonstration purposes and does not constrain the process once (4.2) and (4.3) have been determined. The height of where a particle will rest for a given horizontal position, j, along the

70 I S(il + i)- Mlw (i) S(j2 + i)- Mweri) max ( ----.~ rmn I i I I I I I I I I I I I I I l l I m I ii i i2 J=J1 J=J2 Figure 4.3: Illustration of equation (4.4) for two different points (jl,j2) on the surface S(j). The dilation of the surface is shown for each point on the surface where the maximum is highlighted by a circle. This is the height of where the center of the ellipse would lie for a given horizontal location, j. The minimum of these heights defines the point of minimum potential energy. surface Sp is given by the morphological dilation (Figure 4.3) height(j) = {maxi (S (j + i) -Mower()) -m < i < +m} (4.4) The minimum height attainable used by step three in the procedure is zp = min (height(j)). (4.5) Successful algorithms have been implemented for both two- and three-dimensional elliptical particle distributions, and examples are shown in Figures 4.4 and 4.5 (to make the three-dimensional problem more efficient, a local minimization over an area of several particle diameters was performed rather than global minimization, indicated by (4.5)). The two-dimensional example clearly shows the algorithm's ability to generate very dense arrangements of particles in a short amount of computer time, two characteristics that make the algorithm amenable to Monte Carlo dense media volume scattering analysis. Furthermore, an analogous dense medium three- dimensional example is given for nonspherical particles with a Gaussian size distribution. Results from the packing algorithm can be used to determine physical characteristics such as volume fraction, correlation length and function, and pair distribution function, of a IW

71 80 0 20 40 60 N._ 40 20 0 20 40 60 80 100 Horizontal axis, j (mm) Figure 4.4: Demonstration of the two-dimensional packing algorithm. Particle radius is 2mm ~ 0.4mm. In this case the volume fraction is 0.8. This particular simulation took 23 seconds on a Sun 10 workstation.

72 2O::iij 15i:::~"~~~~'i~'j::~:~i ~ s E i:::::::::::::: E~:.:,j:~~:~:, V ~:::::::::::::::a ~ ~~~2~~ N~:~::~::~::~::j::._f~ i::~~~::~:~::: x ~ ~~~~~~~~~~~~~~~~~ l 0 ~ ~ ~ ~:~ ~~ i::::i:::~:iii:::::x._o~:::~~~~:':~"::''~:: IZr q)'""~ >:: 5 r ~ ~ ~ ~:~:;~~:~~~~ -~~ I'fo: 0r ~:~~~~~~~~~~:::::::~: ~ x~~ 0~~ 5~ 1 0152 Horionta axi s j (mm takn utof 5mm cbe iledwih elisoda p a tce wth ai f2m4.4 m an IBM SP Parallel computer.55 Parices of arbitrary- ~shpmybeodldtmima additional~:~~ ~~: computat~:~:ion cost.~::~::

73 50 0 10 20 30 40 50 Horizontal axis, j (mm) Figure 4.6: Demonstration of the packing algorithm for two-dimensional rocks. Rock surfaces are generated by a random walk about a uniform radius 2mm ~t 0.4mm. rThis particular simulation took approximately 30 seconds on a Sun 10 Workstation. volume under study. Future work in this area may be to include impurities such as water and determine the coating of individual grains through the implementation of simple physical processes such as surface tension and entropy maximization. Other natural substances, such as snow, may also be simulated by reducing the minimization of potential energy constraint, which can be accomplished by making individual particles stick once they reach the surface of the packed layer. The simplicity with which the algorithm can be changed for particles with a general cross section is demonstrated in Figure 4.6 where packing of two-dimensional rocks is simulated. In this figure the standard packing algorithm is used with a random walk modification of the individual particle surfaces. E N ~20 10 0 0 10 20 30 40 50 Horizontal axis, j (mm) Figure 4.6: Demonstration of the packing algorithm for two-dimensional rocks. Rock si~rfaces are generated by a random walk about a uniform radius 2mm ~ 0.4mm. This particular simulation took approximately 30 seconds on a Sun 10 Workstation. volume under study. Future work in this area may be to include impurities such as water arid determine the coating of individual grains through the implementation of simple physical processes such as surface tension and entropy maximization. Other natural substances, such as snow, may also be simulated by reducing the minimization of potential energy constraint, which can be accomplished by making individual particles stick once they reach the surface of the packed layer. The simplicity with which the algorithm can be changed for particles with a general cross section is demonstrated in Figure 4.6 where packing of two-dimensional rocks is simulated. In this figure the standard packing algorithm is used with a random walk modification of the individual particle surfaces.

74 4.2 Permittivity Correlation and Pair Distribution Functions One of the principal uses of the packing algorithm is in the calculation of the permittivity correlation and pair distribution functions which can be used in conjunction with coherent theoretical approaches to compute scattering in inhomogeneous random media. One such classical approach is the Born approximation (Figure 4.7) described in Chapter 3, where the permittivity fluctuations act as distributed sources in an effective homogeneous medium. The scattering solution is then determined using a perturbation series. It was shown that for a random medium where the variation of the fluctuating permittivity is relatively small, the scattering coefficient from the medium can be directly calculated from the Fourier transform of the correlation function given by C(T- r') - (f (r)f()) (4.6) where jf(j) is the fluctuating part of the permittivity function. The most commonly accepted form for the correlation function is the exponential function which was derived originally by [Debye et al., 1957] for a sparse random collection of spherical particles and is given by C()- e-I/ro. (4.7) The corresponding two-dimensional power spectral density is 2 0k2- / (4-8) Here, the parameter To is related to the mean diameter of particles in the random medium. For dense random media, experimental and numerical methods have been attempted to determine C(r - T') from the recorded samples of the medium under study. For most practical cases, the behavior of the correlation function when r' is in the near vicinity of r does indeed resemble that of an exponential or Gaussian function. Because of the ease with which these functions can be manipulated algebraically, their use has become widely prevalent in the literature. A difficulty that arises however, is that the power spectral density is proportional to the integral of the correlation function over all space, and thus estimation of the correlation function for only small values of r - 4 may not be sufficient for accurate estimation of the power spectral density for different ranges of observing frequencies. The packing algorithm described in the previous section is an ideal tool for determining the validity and range of applicability of assumptions made in the determination of this correlation function. Furthermore, since the packing algorithm is capable of making a full

75 gki \ kS r = 23 (iandom med'i-al zi) kz 2 klzi Figure 4.7: Geometry for a two-dimensional random medium. Medium 0 is considered free space and medium 1 is a random media consisting of particles as described in Figure 4.4. Monte Carlo realization of a random medium, it may be used to enhance the understanding of the physics behind the scattering mechanisms in a very dense medium. Computationa]ly, it is much faster to perform simulations in two dimensions to first develop an understanding and then in the future to extend the results to three dimensions. In two dimensions, it was shown in Chapter 3 (see (3.108)) that the radar cross section using the Born approximation is given by,V(O) = - Xo;14C hh(1) = kIYo 141C 1 K 2 27k2g(2kX, -2k1zi) I 2 2k'i cos 0i kzi k = ko sin Oi klzi = klcosji (4..9) where 90 is the angle of incidence, Xolj and kYoli are the transmission coefficients from medium 1 to medium 0, kim is the mean field propagation constant due to the average dielectric, (c), k'" is the imaginary part of the propagation constant in the scattering medium and the power spectral density b is ( k, kz) = 1J C( r-r )eik Fd2r. (4.10) The unknowns in this formulation are the propagation constant and the correlation function for the random medium. The propagation constant can be either the spatial average of the dielectric constant within the random medium (Born approximation) or derived from a

76 mixing formula (distorted Born approximation). In what follows, the Born approximation will be developed for a prototypical example in which the described packing algorithm will be used to provide the correlation function. After illustrating the two-dimensional correlation function with the associated power spectral density, a comparison will be made with similar results obtained by using an exponential correlation function (4.7) having the same correlation length. In the calculation of the permittivity fluctuation correlation function from the collected samples of the medium, it is common to generate the autocorrelation function of each sample and then average the resultant autocorrelation functions. However, a difficulty may arise from this procedure: the Fourier transform of the correlation function may become negative over some partial frequency range. This is a result of both a poor estimation of the tail of the correlation function and possible in the sample cross section, the former being a consequence of the finite size of the sampled medium. To circumvent this difficulty while keeping the size of the sampled medium within realistic dimensions, the power spectral density should be computed directly from the Fourier transform of a profiled sample of the random medium. Since the correlation function of the permittivity fluctuation process is continuous, it implies that the process is continuous in the root mean squared sense, which in turn implies that the power spectral density can be directly computed from )(kA, k z) = lim { X-oo XZ Z - oo X Z 2 E f I (x, z)ei(kxx+kz)dxdz. (4.11) The example considered here simulates the observation of a two-dimensional layer of sand particles with permittivity cs = 3.15 + j0.005, (e) = 2.54 + j0.004, (CE) = 0.56, and volume fraction f = 0.8. Individual realizations of Gaussian-distributed packed particles were created using a mean radius r =2 mm with a standard deviation of 0.4 mm for the two principle axes and a uniform distribution for the orientation angle of the ellipses. Observation frequencies used for this discussion will be at a below-resonance frequency of 10 GHz (AS/2 = 8.5mm) and an above-resonance frequency of 35 GHz (A,/2 = 2.4mm), where resonance refers to the average half wavelength (AS/2) dimension of the scattering particles. Over three hundred packing realizations of the dimensions 0.2 x 0.2 m were used to determine the power spectral density directly from the Fourier transform of the packed profiles given by the packing algorithm. The correlation function and the power spectral density are shown in Figures 4.8 and 4.9, respectively, where the nonaxial symmetry

77 1 I 0.8-, 0.6 -0.4 0.2, 0 ' 04 0e 15 50 -5 (,~~ -10 ~ s~ Figure 4.8: Two-dimensional correlation function for the particles in Figure 4.4. The correlation function is not axially symmetric. patterns are clearly displayed. The strongest peak in the power spectral density lies along the diagonal, as would be expected in a tightly packed array of two-dimensional spheres with a narrow size distribution. A wider distribution would give a more axially symmetrical two-dimensional spectrum. This can be demonstrated by using a Rayleigh distribution to determine particle radius rather than using a Gaussian distribution. Because the packing algorithm has a limited resolution due to the discretization of the particle ((4.2) and (4.3)), the Rayleigh distribution is truncated at the very small scale. Furthermore, for the sake of comparison with previous theoretical and experimental work, particle shapes used with the Rayleigh distribution were made circular rather than elliptical, this has the effect of reducing the randomness of particle locations, but this is an effect that is ameliorated by the wider particle size distribution. The probability density functions of both distributions are shown in Figure 4.10, where it can be seen that the Rayleigh distribution is wider by about a factor of three. Thus for the wider particle size distribution we would expect that the correlation function and

78 -55 - -80 4000 40000 4000 t -200000 2.00 0 t -4000 -4000 Figure 4.9 The power spectral density derived from data shown in Figure 4.4. Strong Fipea l he i go l higlight t odt in the correlation function seen along this axis. Arcs drawn along a constant radius are p as a function of incidence angle.

79 Figure 4.10: Gaussian (solid line, -) and truncated Rayleigh (dot dashed line,.-) probability distribution functions of particle radii used in the packing algorithm. Both distributions have the same mean. spectral density to more closely approximate that of the exponential function (which would be correct for completely random particle placement) and this is indeed shown to be the case. Plots in Figures 4.1-1a through 4.11c illustrate the correlation function and its power spectral density along three principle axes of the two-dimensional volume for both particle distributions and the exponential correlation function. It can be seen that at both short and long distances, the correlation function of both particle distributions and the theoretical exponential function agree very well, but it is at the midrange distances where the functions can be seen to differ. These differences become more apparent in the low-frequency region of the power spectral density for the derived correlation functions. Note again a better agreement of the Rayleigh-distributed particle result with the theoretical exponential function, but still there is an appreciable error at the low-frequency region. Referring to the Gaussian distribution function, the effect of the strong peaks in the power spectral density function along the XZ axis in Figure 4.11b can be clearly seen in the extended periodicity of the correlation function along this axis. Again, this is indicative of the natural ordering of the system for tightly packed arrays of particles with a narrow particle size distribution, and we note that it is less of a feature when the particle sizes follow a Rayleigh distribution. An important consequence of this lack of symmetry, however, is that it would be inappropriate to calculate the correlation function along one axis of the distribution (such as the function in Figure 4.11a) and to enforce axial symmetry. The effect of such an action would be that the resulting power calculated through (4.10) would

80 not be positive for all frequencies because the chosen correlation function along the one axis is not, and cannot be, the correlation function for an axially symmetrical medium. In effect, then, if an exponential or Gaussian function is not to be used, the only appropriate way to calculate the correlation function is to completely sample the fluctuations in twoor three-dimensional space. Furthermore, and more importantly, as is shown by Figures 4.11a through 4.11c, the deviation of the numerically derived power spectral density from the exponential power spectral density is as high as an order of magnitude at the lower frequencies, which is precisely the region of validity for the Born approximation. As indicated by (4.9) and demonstrated at the end of Chapter 3, only a portion of /(kX, k,) is "visible" at a particular frequency. The loci of the spatial frequencies as the incidence angle varies from 0~ to 90~ is an ellipse with semi-axes of 2ko and 2k1m as illustrated in Figure 4.12 for this particular example. In this figure it can be seen that while the spectral density is not axially symmetric, the visible regions at each of the observing frequencies are nearly so. The effects of the asymmetry are most strongly observed at frequencies near the resonance of the particles within the medium. This figure can be used to understand the relationship between observing frequency and physical structure of a random medium. Decreasing particle sizes would have the effect of moving features radially outward in the figure, the same effect as decreasing frequency. Narrowing the particle distribution would have the effect of increasing the magnitude of the peaks shown in the figure and further disturb the high-frequency agreement of the observed backscatter with the backscatter response expected for a medium whose particles follow an exponential function. Another interpretation of the figure shows that the particle size distribution of a random medium may be explored by frequency sweeping near the particle resonance at angles slightly greater than 45~. Given the derived correlation function for a specific particle size distribution we can carry out the remaining mathematical operations in (4.9) to get the observed backscattering cross section, which is shown in Figures 4.13 and 4.14. In Figure 4.13 the backscattering coefficients of the medium for both the vertical and horizontal polarizations at 10 GHz and 35 GHz are shown and compared with those derived based on an exponential correlation function for a Gaussian particle size distribution. Figure 4.14 is the same as Figure 4.13 for an observation frequency at particle resonance (21 GHz), where differences between observed and predicted backscatter response are expected to be the greatest. Referring to Figures 4.8 through 4.14, the following observations are in order: (1) differences between the backscattering coefficients estimated from the correlation function

81 c 0 0 0.8 0.6 0.4 0.2 0 I I I I I I A I I I I I I I -r 7 - -%... 0 -15 -10 -5 0 5 10 Distance along the horizontal X-plane (mm) 15 20 10-6 10-7 Q_ 10-in -8 -4000 -3000 -2000 -1000 0 1000 2000 Wavenumber along the horizontal X-plane 3000 4000 (a) I 0 a1) 0 0 0.8 0.6 0.4 0.2 0 A.,., I T I I I I I I- -- I I -0I.2 I I I I I I 0 -15 -10 -5 0 5 10 Distance along the diagonal XZ-plane (mm) 15 2 0 (b) Figure 4.11: Slices of Figs. 4.8 and 4.9 taken along the principle axes of (a) 0 = 0~, (b) 0 = 45~, and (c) 0 - 90~ (next page). The solid lines (-) and the dot-dashed lines (.-) show the numerically derived values for the Gaussian and Rayleigh particle size distributions respectively. The dotted lines (...) show the equivalent correlation function and spectrum for the best fit exponential.

82 c 0 0 0 0.8 0.6 0.4 0.2 0 A 1) I I I I I I I -0.,~ I 0 -15 -10 -5 0 5 10 Distance along the vertical Z-plane (mm) 15 2 10-6F........ - '........... -... *'-\ /." -10 - 10-8 I8-II I I -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 Wavenumber along the vertical Z-plane Figure 4.11 cont'd: Slices of Figs. 4.8 and 4.9 taken along the (previous page), (b) 0 = 450 (previous page), and (c) 0 = 90~ previous page for details (c) principle axes of (a) 0 = 0 (this page). See caption on Spectral Intensity, Phi 2 y 0 Kx Figure 4.12: Illustration of the visible region of the Born approximation at two different frequencies shown overlaid on a contour map of the power spectral density (Figure 4.8) for the Gaussian particle size distribution. Note that the greatest angular variation due to the packing structure occurs near resonance at f = 21GHz.

83 O 10GHz \ 35GHz -6............. -66 Angle of Incidence (degrees) Figure 4.13: Results of the two-dimensional Born approximation for both 10 and 35 GHz for the Gaussian particle size distribution. Solid lines (-) are for vertical polarization, dashed (- -) are for horizontal polarization and the symbols (o and x8) indicate vertically polarized results of the Born approximation using the best fit exponential (Figure 4.1 1). 12-'-'-'-' ---' 21 GHz - 6- ' ''' '' 0 \ -4 10 2 3 4 5 6 7 8.o..; \ \0 -1 i J i ' \ i 0 10 20 30 40 50 60 70 80 90 Angle of Incidence (degrees) (C — Figure 4.14: Results of the two-dimensional Born approximation at particle resonance (f = 21GHz) for both numerically derived (Gaussian particle size distribution) and exponential correlation functions. Solid lines (-) and circles (o) indicate vertical and dashed lines (-) and crosses (x) indicate horizontal polarizations. The symbols represent the theoretical exponential curves in both cases.

84 derived from the Monte Carlo simulation and the exponential function can be as high as 10 dB; (2) angular variation may depend strongly on the observation frequency and particle size distribution (see Figure 4.14); and (3) the correlation function and the respective power spectral density are not likely to be radially symmetrical, but this may not have a strong effect on the observation with the exception of frequencies near resonance. Another statistical function relating to particle positions that is commonly used in the study of random media is that of the pair distribution function, p(Fj1r,), which is the probability of finding a particle located at a position Tj, given a particle located at position Y, [McQuarrie, 1976] (see Section 3.2.4). The knowledge of the pair distribution function becomes necessary where the effective propagation constant of a dense random media is to be determined [Tsang et al., 1985; Lax, 1952]. As the packing algorithm can be used to determine the correlation function, it can also be used to determine the pair distribution function. Figure 4.15 illustrates the aggregate pair distribution function for the particle arrangement consisting of two-dimensional ellipses shown in Figure 4.4. In this figure, the pair distribution function is shown as a gray scale image in two dimensions where the brightness is directly proportional to the probability. As expected, a dark area covers the central region indicating the zero probability of having two particles intersecting one another. More remarkable is the ringing effect along the diagonals of the image, indicating periods of high and low probability, which is reflective of the results obtained by the correlation function. This effect is highlighted in the graph inset in Figure 4.15, where the pair distribution function is shown quantitatively along the horizontal, vertical, and diagonal axes. The strong peaks along the diagonal axis are a result of the dense packing with a narrow size distribution. In comparison, the Percus-Yevick pair distribution function predicts an axially symmetrical function for single-sized spherical particles. The example given here clearly is not so and cannot be so for monosized particles under the influence of gravity. This does not disprove the Percus-Yevick pair distribution function, but it does call into question its application to macroscopic granular media under the influence of external forces. Additional treatment of the sensitivity of the pair distribution function to particle extraction was given in Section 3.2.4. In this treatment it is shown how the pair distribution function does not change when particles are randomly removed from a simulated arrangement of particles. The central assumption to this observation is that when particles are removed, the remaining particles are not allowed to physically rearrange themselves. This is an arguably unphysical phenomenon, but it emphasizes how similar arrangements of par

85 -25 -20 -15 — 10 E E -5 N '~ 0.( _ 5 10 15 25 -20 -10 0 10 20 Horizontal Axis, j (mm) Figure 4.15: The pair distribution function, p(arj xi) for the packed array shown in Figure 4.4 (Gaussian size distribution). Dark areas indicate low probability and bright areas indicate a high probability of finding another particle at Tj given a particle at Tr. The inset graph gives the quantitative pair distribution function along the vertical, horizontal and 45 degree axes.

86 tides may result in considerably different pair distribution functions. This point is especially important for those who use the Percus-Yevick pair distribution function in the modeling of snow. 4.3 Conclusions In this chapter,a new approach for determining particle arrangements in a dense medium has been introduced. The approach is computationally efficient, flexible, and was demonstrated to work in two dimensions as well as three for nonspherical particles. A twodimensional "rocks" example was also given that demonstrated the ability of the algorithm to work with particles of arbitrary shape. The algorithm was then used to compute the correlation function for two different size distributions of particles, one Gaussian and the other Rayleigh, and it was highlighted that the proper way to calculate the correlation function was through an averaged periodogram in the spectral domain rather than averaging individual correlation functions or assuming radial symmetry from the outset. The correlation functions of the two particle size distributions were compared with a theoretical exponential correlation function, where it was shown that the derived functions deviated from the theoretical at low frequencies but agreed well at high frequencies. It was also shown that the wider Rayleigh distribution gave a better fit overall to an exponential function, as would be expected. Derived correlation functions from the packing algorithm were used to compute radar backscatter via the Born approximation, where it was shown that an increase in the backscatter response as a function of observation angle can be related to angular asymmetry in the correlation function. This asymmetry may be due to the natural ordering of particles when arranged under the influence of gravity. Finally, another application of the packing algorithm was given as it applies to the determination of the particle pair distribution function, which is a key unknown in the quasi-crystalline approximation, and an argument was presented as to why the more common Percus-Yevick pair distribution function is not appropriate for use with granular media.

CHAPTER 5 T-Matrix Much of the work that follows this chapter relies on solution techniques to determine scattered fields from collection of discrete scatterers confined to a volume situated in free space. The purpose of this is to utilize numerical solutions to Maxwell's equations to provide an in depth knowledge of how fields are behaving within different varieties of random media. Given this detailed knowledge, a technique will be presented (in the following chapter) of how to determine a fundamental electromagnetic characteristic about random media, that is, how to determine the propagation constant, or effective permittivity. In two dimensions, the method of moments [Richmond, 1965; Richmond 1966] can be used to determine volume currents induced by an incident field on a collection of scatterers with known positions and permittivities. By applying the free space Green's function in two-dimensions, the scattered field may determined. In the two-dimensional approach, the incident field is split into TM and TE components, and solved separately. This is an exact solution to Maxwell's equations insofar as integrations can be carried out accurately (either numerically or analytically) over each of the components of the medium. In its own right, the method of moments is well developed, explored and documented [Harrington, 1966; Balanis, 1989; Chew, 1995] and need not be reviewed here. In three dimensions however, application of the method of moments is much more cumnbersome due to the large number of unknowns (vector volume currents) inherent to the problem. As a consequence, it is unrealistic to solve any large scale problems using the method of moments. To circumvent this difficulty, the T-matrix (transition matrix) method, first proposed by [Waterman 1971; Peterson and Strom 1973], has been further developed by [Wang and Chew 1993], to utilize recursive techniques to solve for the scattering die to a large number of scatterers. Because the development and application of the recursive T-matrix solution technique is fairly recent, it is reviewed here for clarity, completeness, and also to highlight its uses and limitations. Outside of the original work published by Chew, 87

88 this is the first full treatment and implementation of the generalized recursive aggregate T-matrix algorithm [Wang and Chew 1993]. This in itself is a significant accomplishment because (i) it offers an alternative perspective of the recursive T-matrix technique, (ii) the implementation revealed several errors in the original publication of the three-dimensional work, and (iii) limitations and abilities of the algorithm can be explored in a more objective format. What follows is first a formulation of the generalized Recursive Aggregate T-Matrix Agorithm (RATMA) and then a number of examples that are either fundamental examples scattering from a collection of spheres or examples taken directly out of the work published by Professor Chew. While many of the examples illustrate the very strong capabilities of RATMA, one example in particular was constructed to demonstrate a fundamental limitation of the algorithm. This limitation itself has not been discussed or presented in the literature and in this light is an important contribution. 5.1 Formulation The formulation of this technique is based upon the spherical wave expansion of electromagnetic fields. These functions, composed of the spatial vectors M and N given by (3.64) are ordered into a vector i given by (3.62). In this context, the term vector may refer to field directions in space or to the ordered vector of basis functions T. To minimize the confusion that may result, vector field quantities normally specified by bold face type and an overbar will revert to non-emphasized characters. Thus, a given incident field, Einc may be specified as nc -t-inc E = Rg4 a (5.1) indicating that the transpose (indicated by a superscript t), or horizontal, basis function vector is multiplied by a vertical vector of incident field coefficients specified by ai'. For plane wave incidence, the incident field coefficients are given by: inc = [M aN (5.2) a) = (-1)ie- 1) [Oit-(cos(0)) - b-(cos())] ainc = (-l)q-iP e-i 1) [5 q(coS(0)) -itm(cos())]. The functions, s and t are related to the associated Legendre polynomial and its derivative and can be found in [Tsang et al. 1986, pg. 173]. In the above expressions, p refers to the

89 order of the dipole moment, and q is the dipole number such that q = {-p, -p + 1,..., p} and p = 1,.., Pmax, where Pmax is the truncation number for the number dipole moments used in the calculation. For a given Pmax, the number of elements in a jnc is 2Pmax (Pmax + 2). Similarly, the scattered electric field can be expressed in terms of the scattered field coefficients, aca, such that Escata sca The T-matrix (or transition matrix) is then defined by the relation sca _ Tainc sca T (5.4) For a single sphere, centered at the origin, the T-Matrix reverts to the Mie series and is given by (3.70). For the more general problem, it will be necessary to sum together the effects of a number of spheres, displaced from the origin by the vector Tri. For a single sphere, this problem can be easily addressed using the vector translation theorem, as in E a = t301(rr1). T /30(rlr)ai. (5.5) In the above equation, the vector translation matrix Pij indicates a translation of spherical wave functions centered around the jth coordinate system to the ith coordinate system. Equivalently, /3j is the same as Rga found in (3.27). The notation is changed slightly here to be consistent with notations that differ between the authors [Tsang et al., 1986] and [Chew 1995]. In (5.5) it can be seen that the incident field vector in global coordinates is transformed to the local coordinates of the scatterer, multiplied by the T-matrix, transformed back to the global coordinate system and finally multiplied by the spherical wave basis functions in global coordinates. Use of the vector translation theorem is dependent on the point of observation with respect to the vector that describes the translation of coordinates (thlis will be discussed shortly). The fundamental equation for the direct T-matrix algorithm relies on the continued application of this transformation via the application of the different forms of the vector addition theorem. Figure 5.1 illustrates the multiple scattering equation given by [Peterson and Strom 1973; Tsang et al., 1986, pg. 452] N(5.6) w(J - * ( ' i ( a1' (5.6) i —1

90 s - Rg P(k,ro).ainc o 0 7 ' Figure 5.1: Illustration of spherical wave vector translation in the multiple scattering Tmatrix equation from the remaining N - 1 particles plus the incident field. The exciting field coefficients, w^i are unknown and must be solved for by the matrix equation implied by (5.6). Once the exciting fields are known for each particle, the scattered field can be determined by summing up all of the contributions from all of t he particles together. The mathematical rules for the application of the vector addition theorem for spherical waves are often not well illustrated in the literature. In reference to Figure 5.1, for an observation point on the surface s and the translation of spherical wave functions centered on the ithcoordinate system to the jthcoordinate system, the vector translation theorem is given by -(*) = V(kj) i Ijl > Tj (5.8) Rg~(k7i) = Rg~(kFj) '3^i V|T^jl (5.9) The incident field in (5.6) is expanded in terms of regular spherical wave functions, Rg centered at the global origin. Thus, (5.9) is used in (5.6) to translate these functions to a coordinate system centered at the jthparticle. Similarly, the product of Ti(l) w(i) represents outgoing spherical waves, c, whose origin is the center of the ithparticle, and in turn, (5.7) is used to translate the spherical wave functions to a coordinate system whose origin is the

91 center of the jthparticle. Note that the observation point, Tsj will always be less than the distance between the origins of the two coordinate systems, rji. To develop the concept of a recursive algorithm, the T-matrix of a collection of spheres can be written in terms of the T-matrices of the individual components. This new T-matrix is often referred to as the aggregate T-Matrix, designated as T. Extensive work on the recursive technique for determining this quantity has been done by Prof. Chew of the University of Illinois [Chew 1990; Chew 1992; Chew, Wang and Gurel 1992; Chew et al., 1992; Wang and Chew 1993]. The theory behind the recursive algorithm is as follows (Figure 5.2). We begin by assuming that the aggregate scattering matrix from n spheres, 7(n)i has been determined (at the beginning of the recursion, when n = 0, the aggregate T-matrix, 7(n):= 0). The aggregate T-matrix of these n spheres in the presence of an additional n' spheres can be written as n+n Tn(n+n) = t(n) + T(n) '* S,oj ~ Tj(3+n) /j0 (5.10) j=n+l the jthparticle itself), a and /3 are vector translation matrices as described above (a is a vector translation matrix similar to /3 with the exception that spherical Hankel functions are used instead of spherical Bessel functions). Figure 5.2(a) illustrates the different component terms of (5.10). Similarly, the T-matrix of the jthsphere, referenced to the origin, may be written as n+n' Tj(n+-,) /jo = Tj(l) ' jo + T() 'o ~ T(+ j() ~ a i Ti(n+n') ~ Po i=l (5.11) Zj3 The components of (5.11) are illustrated in Figure 5.2(b). In (5.11), Tj(1) refers to the single scattering matrix of the jthparticle in the absence of other particles. In (5.10) and (5.11) the principle unknown is the scattering matrix of the jthparticle in the presence of n + n' scatterers (i.e. Tj(n+n,)/3o). By substituting (5.10) into (5.11), an algebraic expression can be written for these unknowns 5.12)

92 - Tj (n+n') I /i inc - ' (n) (a) Einc, IB, G '' t11 I (b) Figure 5.2: (a) scattering from n spheres in the presence of n + n' spheres, (b) scattering from sphere j in the presence of n + n' spheres. where -j I -. T,(n) ' CfO' i = I;= T ).to + T3.(l) - Ot '' 1 4 i (5.13) Given a solution for Tj(n+,n).-jo, the aggregate T-matrix of n + n' scatterers is written as (see Figure 5.3) n+ —n' (nm+n') = n(n+n') + t30j ' T j(n+n) ' 3jo j=n+l (5.14) After substituting (5.10) into (5.14), an expression can be written for the aggregate T-matrix in terms of known quantities, ~T() and Tj(n+n') '/3j0 n+nf - -(5 = = = -(n+nT) = T(n) + 0 (3oj + 7(n). Oj ) T (n+n, ).j (5.15) j=n+l The critical components of the recursive algorithm are (5.12), (5.13), and (5.15). The

93 A Tj (n+n') IR I I E ~l I I I I I I - I I.1 F I I T I -W T I n(n+n') 'I I AM, \\ // Mr I I first of these, A l,n+l (5.12) may Arn+ri',n+2 Figure 5.3: Illustration of (5.14). be written as a matrix of matrices, as in... A lnn Tn+l,(nj+n') f3n+1,... Ani+2,n+njI Trn+2,(n2+n) )n20(5.1.. A~n,nl+n'i LT n,(n+n') )nn, -Tn+1(1) (n+1,0 + i&n+1,O -rn - -Tn+2(i) (n+2,0 ~ arn+2,O 7r(n)) L +n' (1) ( n+n',O + nri,. T(~) L6) which can be more concisely written as AC - B. (5.17) Furthermore, (5.15) can be written in matrix form as 7'n +n') T (n) + [(3O,n+l ~ =r(rn) a Co,n+i)... (30,n+nI + 7n) 'Cf,n+n/) (.8 -Tn+1(1) (n+1,0 ~ Ctn+1,O 1 ~n [- Trn+r/(l) (3nr+n',0 + Cin+n',O Thin)) I or Trnn)=-(n) ~ V C. (5.19)

94 From the above, (5.17) can be used to solve for C which can then be substituted into (5.19) to determine the aggregate T-matrix. The recursive form of the T-matrix just derived is a superset of direct T-matrix method given by [Waterman 1971; Peterson and Strom 1973; Tsang et al., 1985]. If, at the beginning of the algorithm, n' = N, the total number of spheres to be added, (5.12), (5.14) and (5.15) reduce to a single matrix inversion which is equivalent to the direct T-matrix method. The equivalence between the multiple scattering equation for the direct T-matrix described by (5.6) and the T-matrix algorithm developed in this chapter, given by (5.11), can be seen by right multiplying (5.11) by ainc and making the substitution T (1)' w() =- Tj(N) jo * ac. (5.20) which would make the initial starting point of the two methods identical. The dimension of the arrays utilized in the recursive T-matrix method is an important parameter. Array dimensions are governed by the order of the dipole moment used in calculating the T-matrices of individual scatterers, denoted as Mmax, the order of the dipole moment used for the aggregate T-matrix, denoted as Pmax, and the number of new scatterers added at each iteration, n'. For a given value of Mmax (or Pmax), there are a total number of M = 2Mmax(Mmax + 2) (or P = 2Mmax(Mmax + 2)) components of the electric field (the factor of two coming from a combination of M and N vector fields). The following table lists the dimensions of the utilized by the recursive T-matrix method in terms of these components. variable: Ai- T() 3ji, -ji o 3 joo |oj oj (n) A B C | ) rows: M M M M P P n'M n'M n'M P cols: M M M P M P n'M P P n'M Note that the one matrix inversion that must be performed for each iteration, is on A, which has dimensions of n'M x n'M. For Rayleigh sized spheres, Mmax =1, and therefore a 6n' x 6n' matrix must be inverted. Values of n' on the order of 30 spheres per iteration will maintain the matrix inversion within reasonable limits, but this will be at the sacrifice of some accuracy (which will be demonstrated in the next section). The number of dipole moments required for the aggregate T-matrix is best determined by the formula given by [Bohren and Huffman, 1983] as Pmax = koA + 4(koA)1/3 + 2 (5.21) where ko is the free space wave number and A is the radius of the sphere enclosing the total of N particles comprising the medium under analysis.

95 5.2 Limitations of the Recursive T-Matrix Algorithm The T-matrix algorithm and RATMA are numerically exact in that if a sufficient number of terms of the spherical wave expansion are retained, and if the machine precision is sufficient, an exact solution to Maxwell's equations will result. Computational limits however require a more practical solution to these equations, and limits must be put on PIax and Mfmax as was alluded to above. Experience has shown that the direct T-matrix method converges to a unique solution by implementing (5.21) from above for the value of Pmax and Mmax = 1 to 3 for the component spheres (assuming they are less than one wavelength). Reducing Pmax below the value specified in (5.21) will have the effect of giving erroneous results, while using an insufficient Mmax will constitute a loss of power in the system (i.e. some energy is lost in the higher order terms of the spherical wave expansion). The recursive T-matrix approach (RATMA) however has another set of limitations with respect to the choice of Pmax. This limitation is not governed by (5.21) and thus requires special consideration. This limitation is best illustrated by comparing (5.12) and (5.13) derived using the traditional and recursive T-matrix algorithms for the two-sphere problem. Using these equations, it can be shown that total scattered field from the second (outer) sphere, calculated using the direct T-matrix approach, can be expressed by T2(2) * f320 = T2(1) * ^21 * Tl(l) * a12 * T2(2320 + a21 * Tl(l) * C10o + T20 1 (5.22) Similarly, for the same quantity, the recursive algorithm gives T2(2) * 20 = T2(1) * 20 * 301 * T (1) *,1 2 * T2(2)320 + 21 *(I) * ( o + 32o (5.23) The difference between these two equations lies in the vector translation relationships given by C21 = a20 An (5.24) ^12 = 10 * 02 (5.25) in matrix notation, or AUv,nm (kT12) = Z1 E=-i A (,ijT2o0)RgAij,nm (k-ol) +B,,j (k2)RgBij,nm (krl) (5.26)

96 Bv,nm(krl2) = C-1 E}=-i B,3 (kr20)RgAjnm(krol) +A, (k2o)RgBijnm (kro ) (5.27) in summation notation. In the above [r1o < jr201 and (5.27) is the expanded form of (5.24) in terms of the vector addition theorem coefficients given by [Chew 1990, pg. 591]. This form of the addition theorem is shown explicitly in summation form here to highlight the fact that the summation is truncated after Pmax terms instead of the infinite bound given in (5.27). Returning to matrix notation of (5.24) and (5.25), we see that the right-hand-side performs the operation of transferring vector spherical wave functions centered upon one of the spheres to the origin and then translating to the remaining sphere. Given an infinite number of spherical harmonics (i.e. Pmax = xo), this would be equivalent to directly transforming the vector spherical wave functions from one sphere center to the other (i.e. the left hand side of (5.24)). When Pmax is finite however, the relations given by (5.24) through (5.27) are approximate, and the convergence of this relation with respect to Pmax is a function of both the distance of the spheres from the origin and the proximity of the spheres to one another. This is due to the singular nature of Hankel functions at or near the origin. The effect of this limitation can be observed by calculating the diagonal elements of a21 directly and comparing them with the matrix product 2Q -* p0oi as in (5.24). This is done if Figure 5.4 for three different values of 21 (A/10, 2A/10, and 3A/10) and three different limits for the number of dipole moments used in the spherical wave expansion (Pmax, Pmax + 2, Pmax +4), where the reference Pmax is determined by (5.21) and is related to the magnitude of r2O which is fixed at A/2 in this example. As [r21| gets smaller, the imaginary component of aV21 increases due to the singular nature of Hankel functions. This is illustrated for the first three dipole moments illustrated in 5.4 for the three different distances of r21. As the distance r21 decreases, more terms in the global coordinate dipole expansion are required to reconstruct the large imaginary component of ~a21 from the matrix multiplication between ak20 and 3o01. When the value of r21 is relatively large (3A/10), there are a sufficient number of terms in the global dipole expansion to reconstruct 621. This can be seen in the lower set of three plots of Figure 5.4. As the scattering centers get closer however, and r21 gets smaller, the imaginary component of (X21 increases and a larger number of dipole moments in the global expansion are required to maintain accuracy of the dipole expansion. Thus, for Rayleigh sized scatterers and smaller, spaced closely together, the recursive T-matrix algorithm requires an unrealistically large number of dipole moments to accurately account for the interaction between scatterer centers.

97 106 -exact 10 - - Pmax+4 - Pmax+2 Pmax 104 _ / ~ I: 103 101 10-1 __ 1 2 3 Dipole moment, n origin (a) (b) Figure 5.4: Limitations of RATMA. (a) Imaginary component of the diagonal elements of ca21 computed using the matrix product C20 /0o1 as in (5.24) for three different values of r2I i. For each r2i) (A/10,2A/10, and 3A/10), the exact value of 21 (-) is given alongside the approximate values obtained by the matrix product with different values for the maximum number of dipole moments used (see legend). Pmax is given by (5.21). (b) Physical geometry used in the analysis. Shown are the three different particle positions (open circles) in relation to the fixed particle (shaded circle), and the distance, A for calculating Pmax. Particle diameter is A/10.

98 The consequence of this limitation is that the accuracy of the recursive form of the Tmatrix algorithm in accounting for strong interactions between neighboring spheres depends not only upon the value of Pmax specified in (5.21), but also on how closely together the spheres are packed together and how much they are expected to interact. For the two-sphere problem, this error can be on the order of a one dB uniform loss (i.e. offset of power over all observation angles, as will be shown). For multiple spheres, or random media, the error is less predictable, but overall it should be noted that the error reduces with increasing Pmax and decreasing electromagnetic interaction between particles. 5.3 Examples The following section is a series of examples utilizing the recursive aggregate T-matrix (RATMA) described in the previous sections. The electric field quantities that will be shown are elements of the scattering matrix, S, in the far field such that ~ikor Esca = - SE inc(5.28) kor The far field components of the electric field, M and N, comprising the spherical basis functions, 4 are written as F F = e1- f (cos 0) + is (cos 9)] MN - - kor (5.29) oFF [e s(kcos 0) + it (cos 0)] Npq c kor V2/7r The geometry of the setup is shown in Figure 5.5 where the principle observing plane is shown in reference to the incident electric field vector (i.e. the principle is defined as having its normal in the direction of the vertical electric field, Enc). Single offset sphere Given the scattered electric field for a single sphere centered at the origin [Ruck et al., 1970] as Eo(0, 0), the observed field due to a sphere offset from the origin by a vector r~ can be written as E(0, ) Eo(0, b)6e-ko.(ki-s). (5.31) If kor7 = 5y, the phase of the observed field is shifted by lkoTrs sin qa, and the electric field magnitude and phase appear as they do in Figure 5.6, where both the theoretical and T-matrix results can be seen to match very well.

99 inc v AL kinc inc h ksca Figure 5.5: Illustration of the incident and scattered electric field directions in reference to inc the principle scattering plane defined by the incident vertical electric field, EV. -15 o -2.. -25 - ) -30 / CZ E ~ ' -35 \.a~\ I o-40 0I a -45; -50 0 30 60 90 120 150 180 Observation angle (deg.) 4 r 2 -C a) 0o 0 ~ -2 w -4 0 30 60 90 120 150 180 Observation angle (deg.) Figure 5.6: Magnitude and phase of a single sphere offset from the origin by kor6 = 5y. In this example, k0a = 0.5 and c = 6.93 + iO.l, Pmax = 9, Mmax = 1.

100 z 5' by; radius = a X v inc - sca v cinc Vr h inc sca Figure 5.7: Geometry for two spheres offset from the origin Two non-interacting spheres For the case of two non-interacting spheres, the electric fields determined by (5.31) for a two single spheres, offset from the origin, may be added directly to determine the total field. For two spheres, offset by 6y along the y-axis as in Figure 5.7, the total electric field is given by E(O, ~) = 2Eo(0, a) cos [ko6y sin 0sca] (5.32) As an example, the single sphere solution shown in Figure 5.6 can be adapted by adding another sphere on the other side of the y = 0 plane. Using the same parameters as in the single sphere example, it can be seen in Figure 5.8 that there is excellent agreement between the non-interacting analytical solution and the T-matrix solution. Two interacting spheres In this section, a simple example is set up to explore three effects of interaction between two spheres offset from the origin (Figure 5.9). In this figure, three spheres are shown lying along the z-axis. The first and second spheres are positioned at equal distances from the origin in opposite directions; only one of these spheres will be used at a time in the ensuing two simulations, one where the spheres are separated by a large distance and are not interacting, and the other where the spheres are closely spaced and strongly interact. This setup guarantees that the magnitude of the translation formulas used in both of these examples is the same, and therefore the accuracy of the translation formulas will be the

101 m -10 a) (D = -20 r0) E -30 "0 a -40 o -50 I -rI -60 C 7 I I I I K KI \ /E I I I I I I I I I if I I 11 I X I I' Il II *1 I I I II I I II 11 I I I I - 1 I 1 30 60 90 120 150 180 Observation angle (deg.) ) 30 60 90 120 150 180 Observation angle (deg.) I * L * 'i II*. Figure 5.8: Co-polarized vertical and horizontal electric fields scattered from two spheres offset from the origin by koS = 5. In this example, k0a = 0.5 and e = 6.93 + iO.1, Pmax - 9,Mmax= 1. z-axis A ) interacting / 0o 0 non-interacting Figure 5.9: Two spheres aligned along the z-axis. The position of the inner sphere may be adjacent to the outer sphere, or opposite of it (as shown).

102 same between the two examples (the only difference between the examples is the degree of interaction between the two spheres). The solution for each configuration is performed using the recursive algorithm and the direct T-matrix inversion techniques. The radii of the spheres are ka = 0.63 and the distance of the outer sphere to the origin is kd30 = 6.3 (thus kd31 = 1.3, kd32 = 11.3, and Pmax = 16). As shown in the lower plot of Figure 5.10, the two non-interacting sphere simulation gives equivalent results for both the recursive and traditional T-matrix approaches. These results were previously shown to agree very well with theoretical results for non-interacting spheres (the value of Pmax was calculated by (5.21)). When the interaction between spheres becomes significant however, we see a shift between the result obtained from the direct and recursive T-matrix algorithms. Furthermore, the recursive method more closely approximates the result given by the traditional method as the value of Pmax is increased, even though the minimal value of Pmax specified by (5.21) has been surpassed. The geometry used in this study precludes the possibility of this error being due to a violation of the vector addition theorem (as is stated to be the cause by [Chew and Lu, 1995]) and illustrates the impact of the approximation described by (5.24) through (5.27). The results shown by the two-sphere example illustrated by Figure 5.10 indicates a loss of energy in the system when using the recursive T-matrix algorithm. This power loss can be directly be attributed to truncation of the slowly convergent infinite series expansion implied by (5.24) and is a limitation of the recursive T-matrix algorithm. This limitation is only mentioned in a cursory manner by the literature, and when it is mentioned, the causes are inappropriately labeled as being due to a violation of the conditions for the vector addition theorem (specified in (5.7) through (5.9)). A violation of the vector addition theorem occurs when a source is located within the boundary of the aggregated T-matrix. Clearly, this example rules this possibility out and highlights the fact that errors are due to a truncation of the infinite series specified by (5.27). As a final component of this interacting sphere example, the average error over observation angle between the direct and recursive algorithm over all angles of observation is shown in Figure 5.11 as a function of distance from the origin (normalized with respect to the distance between two adjacent spheres), polarization, and the number of terms used in the spherical wave expansion calculated using RATMA. In this example, Pmax refers to the limit determined by using (5.21), which increases as the distance from the origin increases. As can be seen in the figure, the error decreases steadily (but slowly) as the number of terms in the spherical wave expansion is increased. Furthermore, the errors in the hori

103 Interacting Spheres (vv pol) Co 0 L-6 0 a) cQ cn - Non-interacting Spheres (vv pol) -7 0 -9 0) "-10 0 cO -11 0 30 60 90 120 150 Observation Angle (deg) 180 Figure 5.10: Scattered field (vv-pol) from two offset spheres centered along the z-axis. Both plots show results for both the traditional T-matrix approach (labeled "exact") and the recursive algorithm, where Pmax is given.

104 1 l l I l 0.8 - VV Pmax+2 co Pmax +4\. ^ 0) |0.4 0-i 0) D 0.2 origin 01 0 0.5 1 1.5 2 2.5 3 3.5 Distance from the origin (D/d) Figure 5.11: Average error (in dB) over observation angles between the scattered field of two interacting spheres calculated using the recursive T-matrix algorithm (RATMA) and the direct T-matrix approach for vertically and horizontally polarized fields. The variable, Pmax refers to the number of dipole moments kept in the spherical wave expansion as specified by [Bohren and Huffman, 1983]. The inset illustration (center left) describes the geometry of the simulation which specifies d, the diameter of the spheres, and D, the distance of the inner sphere from the origin. zontally polarized field are much less than for the vertically polarized field because in this arrangement, the coupling between the two spheres is much less. Composite sphere As a final example, a composite sphere may be created out of an arranged set of individual spheres (this example is taken directly out of [Wang and Chew 1993]). This example requires that (i) the volume of the composite sphere to be the same as the large, single sphere and (ii) the individual composite spheres to be electrically small and to have a low dielectric contrast (i.e. to reduce volume scattering) and (iii) Pmax to be larger than the value specified by (5.21). Using a small sphere radius of koa = 0.33, a large sphere radius of koA = 1.26, and sphere permittivity of c = 2.0 + iO.0, the composite sphere can be built out of 57 individual spheres (Figure 5.12). The Mie series solution and numerical solutions for the co-polarized fields are shown to be in very good agreement in Figure 5.13 with a

Figure 5.12: Geometry of the 57 sphere composite. Spheres are allowed to overlap so that the total volume enclosed by the individual spheres is the same as the volume of the composite sphere. best fit uniform dielectric Mie sphere which demonstrates that the recursive T-matrix can be used to some degree for high density arrangements when the interactions between neighboring spheres is not strong. It is expected however, due to the limitations in the recursive T-matrix algorithm, that the imaginary component in this composite sphere is overstated. 5.4 Variations So far the T-matrix method has been demonstrated in it's most basic form, as a method of summing together the scattering effects of a group of spheres. Given the spherical wave expansion of the scattered fields from any scatterer or group of scatterers, their effect may be included into a larger scattering solution via RATMA. Additionally, the T-rnatrix can be derived using the extended boundary condition technique for objects whose surfaces are expressible as mathematical functions (i.e. spheroids, etc.) [see Tsang et al., 1986; Chew 1990]. This technique is based upon extending the boundary condition at the physical surface of the scatterer to a virtual surface that is described by the largest sphere that may be inscribed within the scatterer. Thus, the extended boundary condition technique becomes less accurate as the volume of the scatterer outside of the inscribing sphere increases. For such instances, it may be more appropriate to split the scatterer into a group of subscatterers as was done in the composite sphere example.

106 O,.,,, I l l l l l E-20, 1 '- 0 -V0 V 0 ') -L Observation angle (deg.) Observation angle (deg.) Figure 5.13: Co-polarized vertical (solid lines and circles) and horizontal fields (dashed lines and stars). The Mie series solution is shown as continuous lines and the results from the composite sphere calculated using RATMA are shown as symbols. For koA = 1.26, khoa = 0.33, and e = 2.0 + iO.0, values of Pmax = 9 and Mmax = 1 were used. Elliptical scatterers The T-matrix for spheroids whose surface is defined by the parametric equation -(5.33) 30 + -2 I moments are retained and the T-matrix may be written as it and T it (5.35) 2 -1) 2 kkc^c- 1) -40.31 0 30 60 and t = 0 10 ( 1)A150 Obs 1 [nl (-eg)- 2e1 prolate spheroids(c > a) {2 Cl - (1e2 sin1 e] oblate spheroids(c < a) where e = /1- (min(a, c)/max(a, c))2 is the eccentricity and Aa = (1 - AC)/2. (5.37) is the depolarization factor used in the low frequency mixing formula (3.8). =T 0 0 Tell - - = —N and Tell -- 0 To 0 (5-34 To - ito - t 2 an d T 1 - itl - tl2(.3 2 k3a2c(e~-l) 2 l) 3 In 1+e - 2e prolite spheroids(c > a ) Ac - 2e TI e (5.37) le2 sin oblate spheroids(c < a) the depolarization factor used in the low frequency mixing formula (3.8).

107 Zo Z Z,Z 'y 'Yo Y' 0 x 0X ' X't X"' O ', a Figure 5.14: Rotation of coordinate axes from local coordinates (xo, Yo, zo) to global coordinates (x, y, z) by the Eulerian angles of rotation (a,,, y) Rotation of the T-matrix For non-spherically shaped scatterers, it may convenient to determine the the T-matrix of the scatterer in a local coordinate system and then to rotate that T-matrix relative to the global coordinate system. Such a procedure is accomplished by multiplying by the rotation matrix D given by [Tsang et al., 1986, pg. 193] - - -1 T = DTo D (5.38) The rotation takes place in three steps i) rotation of an angle 7 about the z0-axis in local coordinates to give an intermediate, primed coordinate system, ii) rotation by P about the y'-axis to give an intermediate double-primed coordinate system and iii) rotation of a about the z"-axis to situate the scatterer in the global coordinate system. This series of rotations is illustrated in Figure 5.14. Note that the first rotation would have no effect on a scatterer that is axially symmetric about the local zo-axis (such as the two types of ellipsoids discussed in the prior section). 5.5 Summary The T-matrix method and recursive aggregate T-matrix (RATMA) of computing the scattered fields from a collection of subscatterers has been presented. It was shown that the method provides excellent results for a variety of problems ranging from a single offset sphere, two non-interacting spheres and a composite sphere composed of 57 subscatterers. Limitations of the recursive T-matrix algorithm for the problem of two interacting spheres were illustrated and explained. Additional treatment was given to the determination of the T-matrix for elliptically shaped particles and rotation of the T-matrix from a local

108 to a global coordinate system. The T-matrix algorithm that has been presented will be used extensively in the upcoming sections for addressing three-dimensional random media problems in electromagnetics.

CHAPTER 6 Determination of eff in Two and Three Dimensions 6.1 Introduction In this chapter a new numerical method for determining effective permittivity of dense random media in two and three dimensions is presented. The core idea of the method is to compare the average scattered field of a random collection of scatterers confined within an imaginary boundary with the scattered field from a homogeneous dielectric of the same shape and size as the imaginary boundary. The presentation of the material follows the chronological sequence in which the research was performed. The two-dimensional random media problem is first developed to provide insight into the dependence of the method's convergence on particle size, boundary shape and boundary dimension. The numerical analysis is performed using the method of moments numerical solution to Maxwell's equations where it is shown that the mean scattered field from a random collection of scatterers does converge to the scattered field that would be observed from a similar scatterer composed of a homogeneous dielectric. Results from this analysis are then used to evaluate the two-dimensional version of the effective permittivity theories discussed in Chapter 3 of this thesis. This portion of the work has been published by [Siqueira and Sarabandi, 1996]. The three-dimensional version of this same technique has been more recently developed by implementing the T-matrix algorithm described in Chapter 5 of this dissertation for a spherically shaped confining boundary. Results indicate that similar success can be demonstrated, however, the three-dimensional approach is much more limited due to cormputational limits imposed by the T-matrix algorithm. 109

110 6.2 Coherent Method for Determining Effective Permittivity In this section the numerical procedure for characterizing the effective permittivity of a random medium containing discrete scatterers is outlined. The procedure for determining the effective permittivity is a four step process: (1) a collection of particles with a specified volume fraction and particle arrangement contained within an imaginary boundary is generated, (2) given an incident field, a numerical electromagnetic solution technique (such as the method of moments or the T-matrix algorithm) is used to solve for the scattered field, (3) scattered fields for each realization of the collection of particles are averaged over a set of observation angles to determine the coherent scattered field which is related to the shape and size of the imaginary boundary, as well as the effective permittivity and (4) the effective permittivity of the medium is found by finding the best fit between the coherent field from the Monte-Carlo simulations and the field scattered from a homogeneous medium with permittivity, Ceff having the same boundary as the imaginary boundary. Occasionally, the packing algorithm described in Chapter 4 is used to simulate the arrangement of particles with arbitrary shape, size and orientation for "natural", dense media. Alternatively, a more simple algorithm which randomly places particles into an empty space can be used to mimick the particle arrangement for an ideal fluid. In either case, the arrangement of particles becomes the input step for the numerical solution for the scattered fields. Once this is obtained, the process is repeated many times in order to generate the statistics of the scattered field. This Monte-Carlo simulation determines the average (coherent) scattered field as a function of observation angle, 0, for a random dense medium, (Erandom(O)). By using the coherent scattered field obtained in this manner, it is possible to determine the effective propagation constant which can in turn be used to compare with QCA or even to create an empirical mixing formula for the simulated media. To characterize the propagation constant of the mean field, we search for the dielectric constant of a uniform, homogeneous dielectric whose shape and size are the same as that used to bound the collection of random scatterers and exhibits the same scattering pattern. An example of the general concept of the technique is shown if Figure 6.1 for a rectangular dielectric slab. A numerical or analytical method can be used to determine the scattered field from this uniform dielectric as a function of observation angle and dielectric constant, Euniform((, 0) which is used in conjunction with a minimization technique to find the best fit permittivity for the observed scattered field of the uniform dielectric body to the average

111 o 0 Ei,. -0 aL -:- 0 - * 0... Figure 6.1: Model for numerically determining 6eff for a random medium. field from the random scatterers: 6eff mm [IS niform(E 0) - (Erandom(s))l] E (6.1) To find the minimum of (6.1) by iteration requires calculation of the scattered field for each new trial. For this reason, a simple boundary shape whose analytical solution is known (such as the Mie series solution for a dielectric sphere) makes this calculation simple to perform. Alternatively, in two dimensions, a method of moments eigenalysis procedure [Sarabandi 1994] can be used to convert the inversion of the impedence matrix operation in the method of moments to a simple matrix multiplication. This technique is specific for a dielectric body with fixed dimensions but unkwnown permittivity and is an ideal tool for this analysis. 6.3 Two-dimensional techique 6.3.1 Initial Results The two-dimensional treatment of this method begins with Figure 6.2 which illustrates a 20% volume fraction rectangular sample of a random medium with circular inclusions whose diameter is 1/10th of the wavelength within the inclusions. For an incident field in the broadside direction (from above) a method of moments technique is used to calculate the scattered field for each realization. A sample of the TM and TE polarized scattered fields of 100 realizations of the random medium is illustrated in Figure 6.3 where it can be seen that the variation in the scattered field magnitude between each sample is considerable and nonsymmetric about the forward scattering direction (180 degrees). The random behavior of the scattered field is a clear consequence of the random arrangement of inclusions between each trial. The average of this same coherent scattered field however has characteristics

112 10 - -- - --- - i * * -— P -1 0 rI I I I,, l -25 -20 -15 -10 -5 0 5 10 15 20 25 Figure 6.2: Sample of a random medium with a volume fraction of 20%. Dimensions for this particular example are given in millimeters for a 48 x 16 mm rectangular slab with 1.6 mm diameter inclusions. The dashed line indicates the boundary that confines the sample of the random medium. that are very representative of a uniform medium with a symmetrical shape (Figure 6.4). The statistics of the forward and backscatter field magnitude accurately follow a RiceNakagami distribution as shown in Figure 6.5 [Stark and Woods 1986, pg. 94]. This is expected for the magnitudes of complex vectors whose real and imaginary components are normally distributed, independent random variables. It is expected that scattered fields at the other observation angles follow a similar distribution and that by the central limit theorem, the mean value of the scattered field may be obtained by averaging a sufficient number of samples together. Demonstration of the ability of the method to find a true effective permittivity (comparable to that of a homogeneous dielectric) is given in Figures 6.6 through 6.8 for volume fractions of 10%, 40% and 80% respectively (these volume fractions are illustrated here to demonstrate the the technique works over the entire range of volume fractions; similar results have been obtained for all volume fractions falling within this broad range). It can be seen in the demonstrated cases that the bistatic scattered field from the homogeneous dielectric gives an excellent fit to the entire average bistatic scattered field from the random medium even though only two points in the forward and backscatter direction are used in determining the effective permittivity. This is particularly remarkable in Figure 6.6 (10% volume fraction) where the locations of the scatterers for any particular realization is not o— 5 -fractio~~~~~~~~~~~~~~~~~~n ot 0,4%ad8%rsetvl tee ouefatosaeilsrtdhr

113 TM-Polarization.... I4 r L lu I. (1) -o 0) C E -10. - o -20 a) W -30 F 0 60 120 180 240 300 360 TE-Polarization. I I IIII In I V ' 0 v m O 1 — = -10 'c 0a) E -20 h -30 0 0 Ia -40 L K -50 L 0 60 120 180 240 300 360 Angle of observation (deg) Figure 6.3: TM and TE polarized scattered fields for 100 realizations of the random medium illustrated in Figure 6.2. In this example, the observing frequency is 10 GHz and the permittivity of the inclusions is 3.6 + i0.1

114 TM-Polarization I A m 10 -co 0 a) "0. 0 -C: 0c) E -u -10 0 0 -20 - a) LU -30 -0 60 120 180 240 300 360 TE-Polarization 0 a) -0.-o c -10 E -20 0 o -30 LU W -40 L 0 60 120 180 240 300 Angle of observation (deg) 360 Figure 6.4: TM and TE polarized scattered fields averaged over 100 realizations of the random medium. In this example, the observing frequency is 10 GHz and the permittivity of the inclusions is 3.6 + iO.l

115 Magnitude distributions for forward and backscatter directions 4 E: * - CA rz <U cT3 -~ *4-J ""C IC 0 CL.4 >^ 4-) * V-1 cn >1 *4-> IE, *;C)> 3.1 3.2 3.3 Field Magnitude 3.7 Figure 6.5: Forward and backscatter field magnitude distributions. Numerical simulations (bar graph) follow closely the Rice-Kakagami distribution (solid line). This example is for a square boundary of dimension 1.5A x 1.5A, 40% vollume fraction, TE polarization.

116 reflective of the shape and size of the imaginary boundary. In Figure 6.9 it is demonstrates that the solution to (6.1) is unique. This observation was made by plotting the error function in (6.1) over the space of realistic Ceff (i.e. 1.0 < ff < nclusion and 0.0 < eff < 10 clusion) whereby it was noticed that only one local and global minimum exists for the two-dimensional problem. 6.3.2 Effective Permittivity Dependence on Particle Size Multiple trials of the presented method have been shown to provide values of effective permittivity independent of the boundary shape and independent of the boundary size once the size has passed some critical limit. This following subsection describes the analysis of the dependence of effective permittivity for varying particle size. The shape considered in this example is a rectangular 3A, x lA, dielectric slab where Ai is the wavelength in a dielectric slab with permittivity sc = 3.6 + i1.0 which reflects the real part of the inclusion permitivities. The size is chosen in this example so that at a volume fraction of 100%, the discretization of the homogeneous slab of 10 samples per wavelength will still be valid. The mean particle diameter is chosen to be one of Ai/10,2Ai/10, or 3A,/10 and follows a normal distribution with unity mean and 25% standard deviation (N(1, 0.25)). A minimum of 300 realizations were performed for each volume fraction and the mean scattered field was determined. From this mean scattered field, the forward and backscatter directions were used as sample points in the minimization algorithm given by (6.1). Other points might have been used, but empirical testing has shown that this was not necessary, convergence to a unique solution occurred for all volume fractions and all particles sizes used for both TM and TE polarizations of the incident and scattered fields. Results from the simulations with varying volume fraction and particle size can be compiled onto a single plot that details the dependence of effective permittivity (real and imaginary components) on these parameters. Comparison can then be made with the theoretically derived mixing formula of Polder and Van Santen given in the previous section. This comparison is shown in Figures 6.10 and 6.11 for the TM and TE polarizations. For the TM polarization it is interesting to note the exact agreement with the mixing formula model for small particle size for real permittivity. This behavior indicates that the mean field approach taken by the mixing formula is indeed effective for small particles illuminated by an E-polarized field. As particle size increases, we note a deviation from the mixing formula model. The imaginary component of the effective permittivity is especially illuminating as we note a measured loss term greater than the one predicted by the mixing

117 5 0 -5 0 0 0 L 0 10% Volume Fraction 0 00..0 0....... C O..........0 Oo 0 0 I I '-,1I I -20 -10 0 10 20 TM Polarization TE Polarization 360 Figure 6.6: Algorithm results for 10% volume fraction. Shown here is a.) a sample collection of particles confined within an imaginary boundary, b.) average scattered TM field (solid line) and the scattered field from a homogeneous dielectric with permittivity, 6eff (dashed line), c.) average scattered TE field (solid line) and the scattered field from a homogeneous dielectric with permittivity, ceff (dashed line).

118 40% Volume Fraction TM Polarization /... 10 0 -10 -20......................... 0 60 120 180 240 300 36C TE Polarization 20 i l l 10 -10.... -20 -.q 0 60 120 180 240 300 360 Figure 6.7: Algorithm individual sub-figures. results for 40% volume fraction. See Figure 6.6 for description of

119 80% Volume Fraction TM Polarization,,, 10 0................................................... -10 0 60 120 180 240 300 36 TE Polarization one..... iO 10 0...................................................................................................... -10 -Of) L -;/t....... ) 60 120 180 240 300 360 Figure 6.8: Algorithm results for individual sub-figures. 80% volume fraction. See Figure 6.6 for description of

120 3.5 / 1 2.5^' 0.2 Real permitivity 2 0 Imaginary permitivity Figure 6.9: Error function over the space of realistic 6eff. Differences between the mean scattered field of the random medium and a homogeneous dielectric slab are minimized to determine the effective permittivity. To aid in illustration, the difference curve was inverted so that the peak of the displayed surface represents the best fit dielectric.

121 Effective Permittivity for TM - Polarization with Varying Particle Size 4 x i/10 3.5 - 0 2/10 32X/10 3 3 /10 3 - - model ~' - - model ~" / 4!.>.t'Z E 2.5 a. UJ -a) 2.-, <x 1.5 E c 0.5 cc 0.5 A I/ I // X limiting values <, " -.K X.. - -- - x -.1 - - - A X, X- - 0 0.1 0.2 0.3 0.4 0.5 Volume Fraction 0.6 0.7 0.8 0.9 1 Figure 6.10: Simulation and mixing formula results for TM polarization. Real (-) and imaginary (- - -) components of the effective permittivity derived from the Polder-Van Santen mixing formula compared with effective permittivity obtained by numerical simulation. Symbols indicate different particle diameters of A-/10 (x), 2A-/10 (o), and 3Ai/10 (*). Inclusion permittivity is ~c = 3.6 + il.0

122 Effective Permittivity for TE -Polarization with Varying Particle Size 4 r 3.5 C. 3 E ]. a) 2.5 a) a) I r CZ E 1.5 "0 c m 1 a, w: x ki/10 0 2ki/10 3 3 ki/10 model ~' --- mnrol p" w it,,/ I A I / in u 6 limiting values ' &.l t WI 0.5 -?x I ----- - = — I - -- --- - - I..... I I i! - 0 0.1 0.2 0.3 0.4 0.5 Volume Fraction 0.6 0.7 0.8 0.9 1 Figure 6.11: Simulation and mixing formula results for TE polarization. See Figure 6.10 for more details.

123 formula, a result that indicates that the larger particles are contributing to the scattering process, a factor not accounted for in the mixing formula model. For the TE polarization we see a significant deviation from the mixing formula results. The deviation of the numerical model however is an improvement over the mixing formula in that it follows the expected trend towards the limiting case of unity volume fraction. Similar to the case of TM polarization, the loss term for the H-polarized field demonstrates increased scattering losses due to the larger particle sizes as is expected. This section has demonstrated how the presented method can be used to determine the effective permittivity for a random medium. The behavior of the effective permittivity was shown to behave in a manner consistent with what would be expected for changing volume fraction and particle size. An important question is whether or not the permittivity calculated is dependent on the shape and size of the bounding area used. The next section addresses this problem and offers answers as to how small the imaginary boundary may be to still reflect the large scale behavior of the material parameters. 6.3.3 Convergence Considerations This section seeks to offer evidence for the convergence of the solution for effective permittivity and to test independence of the method from the shape of the imaginary boundary used. To begin this controlled study, small particles of uniform diameter A-/ 10 are used, where, as before, Ai is the field wavelength within the included material whose permittivity is chosen to be ~i = 3.6 + iO.1, the background material being free space. There are some important differences between the results shown in this section and those given in the previous section. In this section, the imaginary part of the scatterer permittivity is chosen to be ten times smaller than that used in the previous section to assure that multiple scattering is allowed to take place more significantly. As a consequence however, since the imaginary component of the effective permittivity is a factor of approximately thirty-six times less than the real component, numerical errors in the estimate of the imaginary component may be much more evident than those for the real component. Another difference of concern between this and the previous section is that scatterers of uniform diameter are used here. Thus for high volume fractions in the range of 50% to the limit of 91%, the resulting collection of scatterer positions will be nearly crystalline and therefore anisotropic. To avoid this additional complexity, only volume fractions ranging from ten to fifty percent are shown here. Already we have seen the ability of the method to achieve an excellent fit between the

124 Dependence of Real Effective Permitivity on Boundary Dimension 2.5 - TM polarization ~ -e - ~~ f=0.5 + - f=0.4 2 -1X 5 x-x f=0.3 co 1.5 - A: -- — 3( - 3( -- f=0.2 -~ - e - o f=0.1 1~ 0 0.5 1 1.5 2 2.5 2 - TE polarization 1 o - - ---- f=0.5 1.6 - + + - f=0.4 1.4 - X f=0.3 - K I f=0.2 1.2 - --- --- ~ ------------------- f=o. 1 1 ------- i ------- i --- —------— i ------- 0 0.5 1 1.5 2 2.5 Boundary Dimension (wavelengths) Figure 6.12: Solution dependence on boundary dimension for TM and TE polarizations. Different lines on the graph represent the average permittivity for a given volume fraction off= 0.1, 0.2, 0.3, 0.4 or 0.5. average scattered field from the random medium and the scattered field of a homogeneous dielectric. There are two questions that will be addressed in this section: 1.) How large should the scattering boundary be? and 2.) How many realizations are required to determine the forward and backscatter field accurately? As it happens, the answer to these two questions are related; the larger the confining boundary gets, the larger is the variation of the electric field in the forward and backscatter directions. Thus as the bounded area increases, so does the uncertainty of the calculated field. To determine the convergence of the method for changing boundary size, we refer to Figure 6.12 which displays the real part of the effective permittivity for the TM and TE polarizations as a function of boundary size for differing volume fractions of particles. We note in the figure that the solution is stable for a wide range of boundary sizes (i.e. horizontal lines). As boundary size decreases farther from the displayed range, it is expected that errors will increase due to an insufficient number of scatterers located within the boundary. At the other extreme of a large boundary size, it is expected that the numbers will become error prone due to the increased uncertainty in measuring the forward and backscattered field (see Figure 6.13 and the following paragraph). The statistics of the forward and backscatter field magnitude was shown in Figure 6.5 to follow a Rice-Nakagami distribution where the mean value and standard deviation are

125 Standard Deviation of Back Scattered Magnitude 0.8.0 = 0.6 **=.. - - a 0.4 -._o — X 0.2 0.5 1 1.5 2 0.5 0.4 N 0.2. c.................O,.. 0.5 1 1.5 2 Boundary Dimension (wavelengths) Figure 6.13: Standard deviation of backscattered field magnitude for three different volume fractions (0.1, 0.3 and 0.5), and both TM and TE polarizations. Similar behavior was observed for the forward scattered field as well. dependent on the boundary shape and size. To explore the size dependence of the standard deviation, Figure 6.13 shows the standard deviation of backscatter magnitude for changing boundary size and volume fraction for both the TM and TE polarization. It can be seen from the graphs that the variation in observed field is not strongly dependent on volume fraction, but does increase steadily as the boundary size increases. If we make the observation that the standard deviation increases by a factor of two between the boundary sizes of 1.5A, and 2.0Ai, to achieve a given accuracy of field magnitude, the larger boundary would require four times as many simulations. The ideal boundary size would then appear to fall somewhere between 1.OAi and 1.5A-. There are three basic geometries that have been tested: a circular disk, a rectangular slab and a square. Real and imaginary effective permitivities for these shapes are shown in Figure 6.14 (note that the imaginary component is multiplied by a factor of ten) as a function of volume fraction for both TM and TE polarizations. Given also for comparison are the 2-D mixing formula results derived in Chapter 3. From this figure it can be seen

126 Comparison Between Different Boundary Shapes 3 0.L I1 - jeff 0 -----, -— " --- —------ ----------------- 0 0.1 0.2 0.3 0.4 0.5 0.6 2 i; 1 5 — c1.5 N 0 ------ - - -------- T --- —------------ -------------— I 0 0.1 0.2 0.3 0.4 0.5 0.6 Volume Fraction Figure 6.14: Effective permittivity comparison for differing boundary shapes. Shown are the real and imaginary permitivities for three different boundary shapes: o - circular disk (diameter = 2Ai), + - rectangular slab (3Ai x A) and and x - square box (1.5Ai x 1.5A2). The imaginary part of the effective permittivity is multiplied by ten so that both parts of the effective permittivity may be displayed on the same graph. Lines indicate results from the two-dimensional Polder-Van Santen mixing formula. that there is essentially an exact match between the real permitivities for all three shapes with some small deviation being noticed for the imaginary part of the effective permittivity. These differences are likely due to difficulties in estimating an imaginary component that is much smaller than the real component of the permittivity. So far, in this section it has been demonstrated that the proposed technique yields stable results independent of boundary shape and size. It was also shown that the forward and backscattered field magnitudes accurately follow a Rice-Nakagami distribution whose variance increases as a function of boundary dimension. One final test is to show the dependence of minimum boundary size on particle diameter. As particle size increases, it is expected that increased scattering and reduction in the number scatterers for a given volume fraction will increase the minimum boundary dimension for which the algorithm

127 inclusion diameter d=3Xi/10 m I A I I II 0.5 1 1.5 2 o 0 d = ki/5 ------— ~ ------ 0 ---- ~9 2.5 0.5 1 1.5 2 2.5 d=ki/10 V v A X\ A /E 0.5 1 1 Boundary Dimension (Xi) 1.5 2 2.5 Figure 6.15: Real part of effective permittivity versus boundary dimension for three different particle sizes: A/10 - x,2A/10 - o,3AA/10 - * and a volume fraction of 40%, TM polarization. will converge. Figure 6.15 demonstrates this behavior by plotting the effective permittivity as a function of boundary dimension for three different particle sizes (A-/10,2A-/10, and 3Ai/10) at a volume fraction of 40%, for a TM polarized field. From this figure it is evident that the Ai/10 particles converge to a stable value even at the smallest of boundary sizes used. As particle size increases, so does the required boundary dimension. Particles of diameter 3A-/10 require a boundary equal to, or larger than, two wavelengths before reaching convergence. Thus, care must be taken to insure that 1.) boundary size is large enough to insure convergence and 2.) a sufficient number of realizations is used to sufficiently calculate the coherent field. 6.3.4 Two-Dimensional Particle Arrangement Methods In Chapter 3, the t hoery was developed for the two-dimensional quasi-crystalline approximation where it was noted that the pair distribution function was the unknown which was most difficult to characterize. In reality the pair distribution function may be es

128 timated in one of three ways: 1.) experimentally, 2.) theoretically or 3.) numerically. Experimental estimation of the pair distribution function is extremely difficult due to the fact that it entails a detailed study of the particle structure for each random medium being considered. A theoretical approach may be used to model the physical distribution of particles in a medium but we may be limited to only a certain class of media which may be described by the mathematical model. Alternatively, a numerical approach may be used to directly model the arrangements of particles via Monte-Carlo simulations based on a computer model the particle deposition. The numerical model, has the advantage of increased flexibility in addressing a variety of physical methods of particle deposition and interaction. In this section a set of four distinct pair distribution functions is presented with the goal of using these to explore the dependence of effective permittivity on the pair distribution function. Each of these functions has been given a simple label (in parentheses) that will be used for reference in a later section. Hole Correction (HC). The hole correction formula is the simplest of the theoretical pair distribution functions and can be used to model an ideal gas consisting of mono-sized particles of diameter b (see Chapter 3). g() p= < (6.2) 1 p> b HC has the effect of eliminating the 12 term in (3.43) and thus it is possible to use QCA-HC as a measure to gauge the effect of the pair distribution function on the QCA formulation. Classical fluid medium (PY). Of the theoretical pair distribution functions, the most common is given by Percus and Yevick [1958] which can be used to model particle positions in a classical fluid of hard spheres. Three-dimensional closed form solutions of PY can be found in Wertheim [1963] and have been shown to agree closely with numerical simulations of particle arrangements [Ding et al., 1992]. In two dimensions, the PY equations have been solved numerically by Lado [1968] (Figure 6.16). The assumptions of hard-sphere PY are 1.) there are no external forces acting on the particles and 2.) the potential energy between particles is infinite when they overlap and zero when their centers are separated by a distance greater than one particle diameter. The first assumption allows the system energy to be directly related to the particle arrangement and the second assumption creates a basis from which this energy may be calculated. These assumptions are the physical parallel of placing particles randomly within a confining volume until a given volume fraction or particle density has been reached.

129 I I I 0 OD 00 6 - Percus-Yevic/Ladc O 0 1 0 Monte-Carlo 00 0 00 4 —, o o 0OOO jI3 2 1 2 3 4 Radial distance, r (particle diameters) Figure 6.16: Simulation of a classical fluid for a volume fraction of 30%. Shown are results from Monte-Carlo numerical simulations and the PY equation in two dimensions solved by Lado. In performing two-dimensional numerical simulations, it was found that volume fractions from 0 to approximately 48% were attainable without resorting to arrangement, initialization techniques (such as given by Metropolis et al. [1953]). Initialization techniques allow high particle densities to be reached by initializing particle positions to a totally packed crystalline array. Individual realizations are performed by perturbing particles from their original positions by a random vector. It should be noted that by resorting to this method, we can no longer guarantee an axially symmetric pair distribution function and for this reason, only volume fractions below 45% are generally used for the analysis component of this paper. It may be however that it is possible to determine a mean pair distribution function which calculates probable particle locations based on radial distance only and ignores the angular dependence by averaging over all angles. While this violates the assumption made in determining (3.61), we can examine the effects of the approximation here. Particle extraction method (extract). Particles are arranged together in a near crystalline fashion using a packing algorithm (PA) described in Siqueira et al. [1995]. Particles are removed one by one without disturbing the lattice until the desired volume fraction is achieved (Figure 6.17). Because of the effect of gravity, the pair distribution function is not axially symmetric. The pair distribution function is shown in the vertical direction (o), the horizontal direction(x) and for the azimuthal average (solid line). Snow-type medium (snow). This two-dimensional simulation utilizes PA to simulate the process of falling snow (Figure 6.18). In this simulation, deposited particles are made to stick near their initial starting points after which the volume fraction is adjusted using

130 1 2 3 4 Radial distance, r (particle diameters) Figure 6.17: Particle arrangement simulation using particle extraction from a near perfect lattice for a volume fraction of 30%. Shown is the average pair distribution function over all angles and the pair distribution function from the vertical and horizontal directions which accentuate the azimuthal asymmetry. Because the basic structure of the lattice remains unchanged for different volume fractions, the pair distribution function does not change with particle number density. the extraction method. Particles are not allowed to find their state of minimum potential energy. This stacking process is reflected in the pair distribution function which illustrates an ordered array in the vertical direction (o) and a disordered one in the horizontal direction( x). In this section, we have discussed four different methods of arranging particles in a volume or area. Each of the different methods were presented to illustrate the variety of pair distribution functions that can be achieved by altering the method of particle deposition in the simulations. The presentation accentuated the fact that the pair distribution functions are not necessarily axially symmetric nor similar, effects that will be further explored in the paper. The resulting pair distribution functions can be used directly in the QCA theory developed in Chapter 2 to analyze the sensitivity of QCA to the pair distribution function, g(4), or we can extend our evaluation one step farther by using the particle position simulations to directly determine the effective constant of propagations for a random medium. 6.3.5 Evaluation of Theoretical Methods In this section, a comparison between the behavior of numerical technique compared to the theoretical methods developed in Chapter 3 as a function of particle size, polarization, volume fraction, scattering/dielectric losses and particle arrangement method is

131 6o 4 x horizontal Radial distance, r (particle diameters) Figure 6.18: Simulation of a two-dimensional snow-type medium for a volume fraction of 30%. Shown is the average pair distribution function over all angles and the pair distribution function from the vertical and horizontal directions which accentuate the azimuthal asymmetry. explored. To begin, a set three sets of six plots (Figures 6.19, 6.20 and 6.21) are presented which demonstrate the theoretical and numerical methods dependence on particle size, volume fraction, polarization and dielectric loss. Figure 6.19 (TM polarization) and Figure 6.20 (TEpolarization) demonstrate real and imaginary (eff performance for particles with a modest amount of dielectric loss (i = 3.6 + i0.1) while Figure 6.21 gives the imaginary component of Reff for particles with low dielectric loss ( = 3.6 + iO.01) for both TM and TE polarizations. The real component of seff is not included in Figure 6.21 because it was found that the behavior of Re(eff) did not change appreciably from the examples given in the demonstrations range from AX/10(kd= 0.33), 2Ai/10(kd= 0.67) to 3Ai/10(kd= 1.0) as the plots go from left to right (k is the wavenumber of the included particles). Referring to Figures 6.19 through 6.21, the following observations are in order: 1.) Re(eeff), T2VI polarization. The real part of the effective permittivity is practically the same for the theoretical methods of PVS, EFA, QCA-PY, QCA-HC and NUM-PY for all volume fractions, particle size and both values of particle dielectric loss (Figure 6.19(a) through (c)). The performance of EFA and PVS degrades as particle size increases, an trackat e behavior of the numerical method well in this aspect because of the inclusion of multiple scattering terms in the QCA formulation. Note that PVS may be used as a reference as particle size varies because particle size is in general not a factor in mixing

132 kd = 0.33 2.5 2.! 1.5 Q) U1) 3r kd = 0.67 5 " /ioi' /, I0 0.1 0.2 0.3 0.4 0.5 (b) 2.5 2 kd = 1.0 9/ la ) 0.1 0.2 0.3 0.4 0.5 (C) 1.5 1.5 1 4 1 0 0.1 0.2 0.3 0.4 0.5 (a) 1 1 0.2 0.2 0.15 7 / / / 0.2 0.15 0.15 0.1 X' 00 @ / a) E 0.1 0.1 0. / 05 mA' #-D + 0 0.1 0.2 0.3 0.4 0.5 Volume Fraction (d) 0. 05 + o. '', o 0 I 0 0.1 0.2 0.3 0.4 0.5 Volume Fraction (e) + 0.05 0 0 0.1 0.2 0.3 0.4 0.5 Volume Fraction (f) Figure 6.19: Comparison between five methods of determining effective permittivity for a TM polarized field: EFA (-.), PVS (...), QCA-PY (o), QCA-HC (+) and NUM-PY (*). Plots (a) through (c) illustrate Re(ceff) and (d) through (e) illustrate Im(cff) as a function of volume fraction using a model of particles (ei = 3.6 + iO.1) suspended in a classical fluid. Particle diameter ranges from Ai/10(kd = 0.33) (a and d), 2Ai/10(kd = 0.67) (b and e), to 3AX/lO(kd= 1.0) (c and f).

133 1.8i kd = 0.33 * * 1.8 1.6 1.6 1.4 K.2 0 0.1 0.2 0.3 0.4 0.5 1 o o. 1 0.2 0.3 0.4 0.5 kd = 0.67 )K X ~/ 1.8 1.6 1.4 1.2 kd = 1.0 X.. /* * Aw 1.4 1.2 1 1 0 0.1 0.2 0.3 0.4 0.5 1 0 0.1 0.2 0.3 0.4 0.5 (a) (b) (c) 0.04 0.03 0.02 E o.C 31 ). - *,3 I -. - - * / o 6 'e Io, 0 0.1 0.2 0.3 0.4 0.5 Volume Fraction (d) 0.04 / 0.03 / X / 0.02 / / 0.01 /)K 0 0.1 0.2 0.3 0.4 0.5 Volume Fraction (e) O.C O.C O.C O.C 04 / * 03 / X 02 /) /X 01 / e s...-.. 0 0 0.1 0.2 0.3 0.4 0.5 Volume Fraction (f) Figure 6.20: Comparison between five methods (QCA-PY, QCA-HC, EFA, PVS and NUMPY) of determining effective permittivity for a TE polarized field. See the caption of Figure 6.19 for details.

134 kd = 0.33 0.2 0 QI I a, E 0.15 0.1 0.05 0.2 0.15 0.1 kd = 0.67 / / 0.2 0.15 0.1 kd = 1.0 I i ) A* +0 X 0 X i + O0 / 0.( 05 /mR A 0O........ I I 0 0.1 0.2 0.3 0.4 0.5 (b) 0.05 + 0 '......... I 0 0.1 0.2 0.3 0.4 0.5 (C) 0 a. w E I a) E 0.04 0.03 0.02 0.01 0 0.1 0.2 0.3 0.4 0.5 Volume Fraction (d) 0.04 0.041 ' I / 0.( 0.( O.( )3 no 0.03 )K * 3* / / / )1 m. /' X A / K v o 0.e.,.,,,,,...... 0 0.1 0.2 0.3 0.4 0.5 Volume Fraction (e) 0. 0. 02 / ** / X 01 /* o ~ e......... 0 0.1 0.2 0.3 0.4 0.5 Volume Fraction (f) Figure 6.21: Comparison between five methods (QCA-PY, QCA-HC, EFA, PVS and NUMPY) of determining the effective permittivity for TM and TE polarized fields. The permittivity of the included particles is ~i = 3.6 + i0.01. The real component (not shown) is essentially unchanged from the previous example when ci = 3.6 + iO.l. See the caption in Figure 6.19 for details.

135 formulae. We may notice therefore, via PVS, that QCA-PY, QCA-HC and NUM-PY agree that phase velocity decreases with increasing particle size. The similarities between the results of Re(eeff) for both low loss and moderate loss inclusions indicate that dielectric loss in the particles does not appreciably effect the phase velocity. This may also highlight the effect that only near particle interactions effect the real part of ceff. 2.) Im(ceff), TMA polarization. The imaginary part of the effective permittivity displays a more complex behavior than the real part as a function of volume fraction and particle size. The theoretical and numerical methods agree well at low volume fractions with the exception of small particle sizes for moderately lossy inclusions (Figure 6.19, part (d)). In this respect, while it is known in the low frequency limit that all the methods agree, the losses indicated by NUM-PY are significantly larger than results given by the theoretical methods of QCA-HC and QCA-PY. We believe that these differences may be caused by an overestimation of absorption losses with respect to the multiple scattering losses for the methods of QCA-PY and QCA-HC. That is, numerically, we detect more multiple scattering than QCA does for moderately lossy particles with sizes that are relatively small with respect to wavelength. We may also note that QCA-PY does an excellent job of tracking the losses predicted by NUM-PY for the larger particle sizes shown in parts (e) and (f) of Figure 6.19. At a volume fraction of approximately 20% the methods of QCA-PY and QCA-HC diverge significantly with a trend reversal seen in QCA-PY due to the contribution from the pair distribution function. In both cases for the larger particles we also see a significant difference between EFA, QCA-PY, QCA-HC, NUM-PY and that of PVS, an effect that is expected because of the scattering losses incurred by the larger particles are not accounted for by the mixing formula. The change in Im(ceff) is considerably more noticeable between Figure 6.19 parts (d) through (f) and Figure 6.21, parts (a) through (c). Proportional to the PVS mixing formula the losses are significantly larger for ~i = 3.6 + iO.01 than for q = 3.6 + 0.1, indicating that multiple scattering is a much more dominant factor in the low loss case. The multiple scattering however does not translate into greater losses overall as we can see that the numerical method predicts lower losses in general than the losses predicted for particles with higher dielectric loss. What this means is that although the field in the = 3.6 + iO.Ol case interacts with more particles, the field does not suffer enough additional absorption losses due to the multiple interactions to make up for the absorption losses incurred with the lossier c = 3.6 + iO.1 particles. We also make a final note that QCA does a comparably well

136 job at predicting Im(ceff) in this case as it did in the previous TM case when c = 3.6+iO.1. 3.) Re(ceff), TE polarization. The real parts of ceff for all methods and particle sizes agree well at the low volume fractions but begin to deviate at a volume fraction of about 20% (Figure 6.20, parts (a) through (c)). This effect is expected because of the necessary non-linear contribution from electrical dipoles to the total field for the TE polarization. The trends of EFA, QCA-PY, and QCA-HC all follow that of the PVS model relatively well despite an even poorer performance for the theoretical methods in the prediction of scattering losses (discussed below). We note that NUM displays the physically expected trend of ceff approaching ci as volume fraction increases towards unity. 4.) Im(eff), TE polarization. We may first note that all methods agree well with respect to losses at very low volume fractions irrespective of particle size. In the zero to twenty percent volume fraction range it is evident that QCA-PY and QCA-HC reflect the correct trend of increased scattering loss for increased particle size. However, for this polarization, the dominant term in (3.43) is Ie (the particle exclusion integral) and not 12 (the integral dependent on the pair distribution function). Thus Im(6eff) becomes unphysically negative be seen through the effect of QCA-PY following QCA-HC into the non-physical domain of negative Im(eeff). QCA-PY and QCA-HC do not adequately describe the multiple scattering losses at fractional volumes greater than 5% and in fact the scattering losses under-predicted by QCA-PY and QCA-HC significantly fall below the lower limit given by the PVS mixing formula. In Figure 6.21, parts (d) through (f) we illustrate the behavior of Im(6eff) for low-loss particles whose permittivity is Ei = 3.6 + i0.01. Without the additional support of dielectric losses, QCA gives realistic values for Im(eeff) for only the smallest of volume fractions. As with the TM case, the losses predicted by NUM-PY are larger in proportion to the losses predicted by PVS but are smaller overall than the losses shown in Figure 6.20 (d) through (f). It may be disturbing to see the poor performance of QCA in its ability to predict losses for TE polarized fields (Figures 6.20 and 6.21). To answer questions related to this effect, it is useful to examine the low frequency behavior of QCA (or more specifically QCA-HC) for both the TM and TE polarizations (Figure 6.22). Only the hole correction formula with QCA is used in this demonstration to examine the effect of the exclusion integral, IC, term in (3.43). The results of QCA-HC are compared with the PVS mixing formula (which is essentially the low frequency limit of EFA) so that a single curve can be used as a reference for different particle sizes. The real component of ceff for both polarizations follows that

137 Comparison of QCA-hc with PVS mixing formula 0.05 __ PVS 0.04- o kd=0.004 TM Polarization x kd=0.04 -0.03- * kd=0.17 + + ~.02+ kd=0.34 + E 0.02 0.01 - + + 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 10x 10-3 8 TE Polarization 0 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Volume Fraction Figure 6.22: Performance of QCA-HC at very low frequencies. Shown is Im(eff) for particles with diameter Ai/800, Ai/80, Ai/20, and Ai/10 for both TM and TE polarizations (si = 3.6 + iO.1). of PVS very closely for all particle sizes and therefore is not shown. Only for particle diameters of Ai/20 and larger does Im(eff) deviate significantly from PVS for the two polarizations, but the deviation for the TM polarization is much milder than that shown for a TE polarized field in which (eff actually becomes negative. We attribute these differences to the higher order singularity found in the Green's function for a TE polarized field (i.e. dipole interactions instead of mono-pole interactions as is the case for TM polarization). The negative trend seen for the TE polarized field is much more difficult for the 12 term in (3.43) to overcome and it is for this reason that we see the degraded performance of the TE polarization results shown in Figures 6.20 and 6.21. Given the extensive analysis presented so far with respect to the performance of QCA, EFA, PVS and NUM, it is convenient to present the results in a more compact form. To this end, Figure 6.23 is presented to graphically illustrate the differences between the three theoretical method's (QCA, EFA and PVS) and the numerical method's (NUM) estimation of 6eff in terms of validity regions. The differences used to determine the validity regions were computed in terms of errors with respect to the numerical method given by theory - NUM err = NU x 100%. (6.3)

138 (a) ~ = 3.6 + iO.1 (c) Ei = 3.6 + i0.1 - 0.67 Q A|m 0. 67 g — x 0.33Xz-Adt03 00 -0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5 Volume Fraction Volume Fraction (b) c = 3.6 + i0.01 (d), = 3.6 + i0.01 (b) F-i = 3.6 + iO.01 (d) ej = 3.6 + iO.01......' ''' 1.0 '~ 1.0 c3 0.67 03 0.67 0.33 r.0.33 0- IA 0 Ol 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5 Volume Fraction Volume Fraction Figure 6.23: Validity regions based on 20% differences between imaginary 6eff for the numerical method of determining ceff (NUM-PY) and the theoretical methods of QCA-PY (light gray), EFA (hatched) and PVS (dark gray) for (a) TM, eq = 3.6 + i0.1, (b) TM, ei = 3.6 + i0.01, (c) TE,Ei = 3.6 + i0.1, (d) TE, c, = 3.6 + iO.01. The plain white regions indicate where none of the theoretical models yield satisfactory results. The numerical method is used as a reference because there are no current exact formulations to determine 6eff in the regions under consideration and the numerical method is the only method presented thus far that displays the expected behavior of ~eff as a function of volume fraction, particle size and polarization. The validity regions were drawn to enclose areas where differences between the numerical and theoretical methods were less than 20%. The validity regions shown are based on a grid of numerical simulations that varied the volume fraction from 5% to 45% in 5% steps and varied particle sizes from kd = 0.33 to kd = 1.0 in steps of one third. It should be noted that because NUM is used as a reference, errors in NUM will manifest themselves as a systematic offset in the presented validity regions of QCA, EFA and PVS. Additional regions of validity have also been drawn for the limiting cases of low frequency and low fractional volume where all of the different methods employed are known to converge. Because of computational limitations, it was not possible to test the high frequency limits of the theoretical methods. Figure 6.23 illustrates well the behavior of the theoretical methods as a function of

139 volume fraction, particle size and polarization. QCA-PY for the TM polarization does well with the exception of kd = 0.33 sized particles (a region which was discussed earlier) and the region of volume fractions extending between 30% and 40%. In this latter region, the reduced performance of QCA is most likely due to an enhanced component of multiple scattering with respect to the other volume fractions. This behavior has been observed both numerically and experimentally [Nashashibi and Sarabandi, 1995]. As a final analysis of QCA, we compare the numerical method's evaluation of Ceff with that of QCA for the array of different particle arrangement methods presented earlier (Figure 6.24). To simplify the treatment, only the angular averaged pair distribution functions were used. At a volume fraction of 30% and q = 3.6 + 0.01 we expect multiple scattering to play an important role in determining ceff. It can be noticed that the exact form of the pair distribution function has little effect on the performance of QCA or NUM for Re(Eeff) for both the TM and TE polarizations, with the possible exception of the snow simulation. More interesting however is the behavior of Im(eeff) to the different methods of particle arrangement. By examining Figures 6.16 through 6.18, it can be seen that the PY simulation yields a loose collection of scatterers whereas the extraction technique generally yields sets of compact groups. The snow-type medium simulation falls somewhere between these two extremes. The enhancement of scattering losses when particles are clumped together is reflected in the results given by NUM; Im(ceff) is least for PY and greatest for the extraction method. Physically, this makes sense, the grouping together of particles acts as a coherent collection of scatterers whose cross section is greater than what it would be for all of the scatterers acting independently. The implications of this result argues strongly against the use of the Percus-Yevick pair distribution function (which models a classic fluid) in the determination of extinction for granular media deposited under the influence of gravity. Snow, sand and soils are good examples of the type of media which fall into this category. 6.3.6 Conclusions A two-dimensional numerical technique based on the method of moments was discussed and developed which may to determine effective permittivity of a random medium. It was shown that the behavior of effective permittivity followed expected trends from low to high density and from low to high frequency. Results of the technique were used to evaluate threoretical methods of determining this same quantity. One of the theoretical methods in common use today is the quasi-crystalline approximation (QCA). Because QCA depends on the form of pair distribution function used, a variety of different particle arrangement

140 a.) TM Polarization b.) TE Polarization x rn ^ r" I 0) EC 1 0.0 o0.0 _Eo0 X NUM O QCA.8 o o.7 PY snow extract )6 o 4 - o0 )4 )2 -X n(I XI X NUM ^ a,> fl A V I.1 I -- U0UA * - 1.4 -I 0 0 0 1.3.. PY snow extract 0.010 6 X "0.005 E 0.0 _ nn; Q I PY snow extract PY snow extract Figure 6.24: Effective permittivity comparison between particle arrangement methods. Volume fraction = 30%, i == 3.6 + i0.01, a.) TM Polarization and b.) TE Polarization. methods were presented. The behavior of the methods and the general behavior of effective permittivity were then studied for a variety of different situations, from which it was possible to determine regions of validity for the theoretical methods based on the presented numerical method. Finally, a comparison was made between the effective permittivities found by the numerical method and the quasi-crystalline approximation as a function of particle arrangement method. From this analysis, the following conclusions have been made: 1.) The real component of effective permittivity is essentially the same for all methods up to a value of kd 1. The TM polarization in this respect is valid for all volume fractions and the TE polarization is valid up to a volume fraction of 20%. 2.) Given 20% error bounds, the losses predicted by QCA for the TM polarization agree quite well with the numerical method (Figure 6.23a,b). It is important however to use the correct form of pair distribution function (Figure 6.24a). 3.) The losses predicted by QCA for the TE polarization agree very poorly with the numerical method (Figure 6.23c,d). At volume fractions greater than 15% it is not uncommon for Im(eeff) to become unphysically negative. We believe this poor behavior is due to the strong singularity found in the Green's function for this polarization. 4.) The effect of decreasing absorption losses of the inclusions has the expected effect of increasing multiple scattering but this does not translate into higher total losses. 5.) The pair distribution function (or similarly, the particle arrangement method) does play an important role in determining extinction (Figure 6.24). It is believed that this

141 is due to the effect of two neighboring particles acting together as a single larger particle. Thus it is inappropriate to use the Percus-Yevick pair distribution when analyzing a random medium whose particle locations do not resemble a classical fluid. The work presented in the paper opens up a number of avenues for future work. In two dimensions the effect of the angular asymmetry in the pair distribution function on QCA can be explored by carrying out the two-dimensional integration in (3.41). It would also be informative to analyze the effect of varying Re(,-) on ceff and to determine an empirical method of determining the effective permittivity based on a combination of QCA and generalized results found by NUM. The most important and useful future work however, will be to bring the presented analysis and methods into three dimensions. In this sense, the study presented here offers an important set of methods and background that can be used to conduct; that investigation. 6.4 Three-dimensional Technique In the last section, the technique for determining the effective permittivity of a random medium was demonstrated in two dimensions to give reasonable, reliable results. The purpose of this section is to bring the technique that was so well explored in two dimensions into the three dimensional domain. The transformation of this process into three dimensions incurs a significant amount of complexity due to the added dimension of the problem. Issues that must be addressed, or issues that limit the ability of the algorithm to explore the full three-dimensional space can be enumerated as follows: 1. The number of particles that fill the bounding space has increased by an order of dimension (i.e. N2 -* N3). The use of shapes other than spheres makes the particle arrangement a separate task within itself. While the packing algorithm can be used to address this problem efficiently, it is not explored here because of the added degree of complexity that this would add to the problem. 2. The variety of boundary shapes that can be used are limited to the small number that can analytically be solved for (i.e. a sphere is the most practical). Solving for the scattered field of other homogeneous bodies with variable permittivity would be very time consuming. 3. Three dimensional numerical methods that can address this problem in an accurate and efficient manner are limited. Initially, it was thought that the recursive T-matrix

142 algorithm (RATMA), discussed in Chapter 5 would be appropriate for this task, but it was later discovered that the recursive aspect of this algorithm has limitations which make the algorithm unusable for this purpose. 4. The number of electromagnetic unknowns in the problem has grown considerably, thus simple trials of the method can be time consuming and computationally intensive. 6.4.1 Convergence Considerations The computational overhead incurred by the three-dimensional problem is significantly larger than that for the two-dimensional problem. Thus it is necessary to judiciously explore the limits of the method of determining effective permittivity so that the computation costs are minimized. The first consideration of interest is to determine how many realizations of the random medium are necessary to achieve convergence. In this case, one realization refers to a random medium sample that is confined within a spherical boundary (whose ideal dimension will be determined later). Within each realization of the random medium it is possible to choose an incident field direction and the orientation of the observing plane over which the scattered field is computed (see Figure 5.5). For different directions of incidence and different orientations of the observation plane, separated by a fixed distance of solid arc, it is expected that the scattered fields would be uncorrelated. In this way, the number of "independent" samples can be increased by changing the geometry of the incident field and the plane for calculating the scattered fields. This type of operation is the equivalent of rotating the sample of the random medium to generate independent samples for measurement. Furthermore, even if the scattered fields derived in this manner are correlated to some degree, the use of multiple independent realizations ensures that the correlation will not bias the answer. This is an important point because the calculation of the T-matrix for a collection of scatterers is computationally expensive compared to the cost of altering the incident and scattered field geometries. With the above considerations in mind, it is informative to display the variation in the observed scattered fields for this problem, similar to what was demonstrated for the twodimensional method. In Figure 6.25, 50 scattered fields calculated from a single realization of a random medium with a twenty percent volume fraction is shown. The variation in the signals between each calculated field is considerable, but when these fields are averaged together, a scattered field pattern that is indicative of a uniform dielectric sphere emerges (see Figure 6.26). This particular example has been for a single realization of the random

143 medium, and in of itself, is not appropriate to be used to determine effective permittivity. The scattered fields between realizations themselves need to be averaged to converge on the final coherent scattered field for the random medium, and this field is in turn used to determine the effective permittivity. Figure 6.26 displays both the vertically polarized averaged field for the scattered fields obtained from a single realization and the averaged scattered field determined using 30 and 50 realizations. From the figure it can be seen that between 30 and 50 realizations of the random medium are sufficient to converge on a single scattered field. When the approach to the three-dimensional problem was first being designed, it was thought that the recursive aggregate T-matrix algorithm (RATMA) developed in Chapter 5 would be an ideal tool for the analysis of effective permittivity. What made this algorithm very appealing was its ability to calculate the scattering from a large number of spheres by recursively adding the effects of spheres in circumferential sets of shells. Thus, a large problem could be tackled with an efficient use of the computer. As was discussed in Chapter 5 however, RATMA is limited in its ability to efficiently calculate the multiple interactions between two neighboring spheres. This limitation has the effect of increasing the losses of the homogeneous sphere because terms in the spherical wave expansion (i.e. energy) are eliminated in the vector translation formulas. This problem is exacerbated as the volume fraction increases because multiple particle interactions also increase with volume fraction. Figure 6.27 demonstrates this effect by showing the scattered field from a 20% volume fraction medium using different quantities of spheres that are added by RATMA at each iteration. The direct T-matrix solution occurs when all of the spheres are added during the first iteration (i.e. n' = N) and is assumed to be numerically exact. The differences seen in the scattered field patterns calculated by RATMA compared to the direct T-matrix approach, illustrated in Figure 6.27, are large enough to eliminate the possibility of using RATMA for the purpose of determining effective permittivity. The consequence of this limitation is that only the direct T-matrix algorithm can be used for solving for the scattered field from a collection of spheres. The direct T-matrix however is computationally limited because it requires the inversion of a large matrix whose dimension is directly proportional to the number of spheres that compose the composite scatterer. I[n general, it is not feasible to accommodate a large number of spheres (greater than 150) into the direct T-matrix algorithm; thus, an upper limit is placed on the diameter of the bounding sphere. This limit is proportional to (i) the volume fraction of the scatterers and (ii) the diameter of the subscatterers. Because of these limitations it is necessary to perform

144 VV-Polarization CD 0) E a).0 c 0. 0) - - LU 30 20 10 0 -10 -20L 0 4 2 -a) C, 3 O C o 0 a) LL 30 60 90 120 150 180 -4 0 30 60 90 120 150 Angle of observation (deg) 180 Figure 6.25: Vertically polarized scattered fields (magnitude and phase) calculated for different orientations of the incident field and the plane containing the observed scattered field. Shown are 50 realizations for a two wavelength diameter sphere with 1/10 th free space wavelength dielectric spheres (c = 6.0 + i0.1) and a volume fraction of 20%.

145 VV-Polarization 30 CrD m () A_0:3 c E 'C) (.LU) 20 10 0 -10 0 30 60 90 120 150 180 C, v0 - 2 Q0 o3 a).4-' 02 C). -2 0l 30 60 90 120 150 Angle of observation (deg) 180 Figure 6.26: Vertically polarized averaged scattered field using (i) 330 different orientations of the incident field and plane of observation for a single realization of the random medium, (ii) 30 independent realizations of the random medium and (iii) 50 independent realizations of the random medium. See Figure 6.25 for further details.

146 VV-Polarization -a 3 r 0) C) "0 E -0 0 C. a) U1 Observing angle (deg) HH-Polarization 30 CD v- 20 a) 0.-o c 10 cn O 0 L0 o -10 a) W 60 90 120 Observing angle (deg) 180 Figure 6.27: Scattered field magnitude from 50 realizations of a 20%o volume fraction random medium calculated using different values of n' (number of spheres added at each iteration) using the recursive aggregate T-matrix algorithm. When n' N, RATMA is the equivalent of the direct T-matrix algorithm (in this case N = 115).

147 an in-depth sensitivity analysis to determine if the coherent scattered field is sufficiently sensitive to accurately determine the effective permittivity. 6.4.2 Sensitivity Analysis A Mie series solution can be used to generate the scattered field patterns of a dielectric sphere with uniform permittivity. Generating a series of curves for different values of permittivity shows that the scattered field pattern is mostly sensitive to the diameter of the sphere and the real component of the permittivity. The imaginary component, for spheres whose dimension is less than several wavelengths, has a secondary effect on the overall pattern and it is expected that this component will be the most difficult to detect. With this in mind, two sets of three plots are displayed in Figures 6.28 and 6.29 to illustrate the sensitivity of the scattered fields on the imaginary component of permittivity for varying diameters of the bounding sphere (0.67A, 1.3A and 2A). From the set of plots in Figures 6.28 and 6.29, the following conclusions can be drawn: 1. The overall shape of the scattered field is sensitive to the real part of effective permittivity while the backscatter direction is the most sensitive to the imaginary component of effective permittivity. 2. To detect meaningful differences in the backscatter field accurate enough to estimate the imaginary component of effective permittivity to within ~0.01, accuracies of approximately 0.05 dB (D = 0.67A), 0.5 dB (D = 1.3A), and 0.6 dB (D = 2A) are required. This translates into 8300, 90, and 60 independent samples, respectively for the different diameters, averaged together. These numbers were obtained assuming a Rayleigh distribution for the backscatter field magnitude. By setting 30 degrees of arc between different incident field directions and orientations of the plane of observation, each individual realization of the random medium provides a total of 336 scattered field patterns, some portion of which are likely correlated; 30 independent realizations give a total of 10,000 scattered field patterns. Thus the accuracies required for the two larger sphere diameters are realizable, but the number required for the small sphere diameter may be unrealistic. 6.4.3 Calculation of Effective Permittivity in Three Dimensions For the example that is about to be developed, glass spheres (e = 6.93 + iO.1) with diameter of A/5 are used as inclusions. Glass beads are chosen for this demonstration for

148 Diameter = 0.67 wavelengths Il 15 10 5 0 -10 -15 -20 0 25 20 15 10 5 0 -5 -10 I I I I I am 1 -,I I.. 30 60 90 120 Diameter = 1.3 wavelengths 150 180 Diameter = 2 wavelengths 60 90 120 Angle of observation (deg) 180 Figure 6.28: Scattered fields (vv-polarization) from homogeneous dielectric spheres of varying diameter. Each separate plot displays the scattered field for a dielectric sphere with a real permittivity of 2.1 and imaginary components (going from top to bottom) of 0.0, 0.02, 0.04, 0.06, 0.08, and 0.10 (these are typical values for a 40% volume fraction of dielectric spheres with a dielectric of 6.93 + iO.l).

149 Diameter = 0.67 wavelengths..... 20 15 10 5 0 -5 -10 C 25 20 15 10 5 0 -5 C 30 25 20 15 10 5 0 C - I - I I I I U 30 60 90 120 Diameter = 1.3 wavelengths 150 180 30 60 90 120 150 Diameter = 2 wavelengths 180 30 60 90 120 Angle of observation (deg) 150 180 Figure 6.29: Scattered fields (hh-polarization) from homogeneous dielectric spheres of varying diameters and imaginary permittivity. See figure 6.28 for further details.

150 several reasons, (i) this material has been measured experimentally by [Nashashibi and Sarabandi 1995], and (ii) the dielectric contrast between the glass and the background (free space) is sufficient to ensure that the beads interact strongly with the incident field. In the first set of examples, the beads are placed randomly within a confining volume, without overlapping, until the desired volume fraction is obtained (volume fractions up to 40% may be obtained in this manner). This method of arranging particles within a space is equivalent to the ideal fluid numerical solution given by [Percus and Yevick 1958] (see Chapter 3) and is spherically symmetric. The scattered field patterns and the best fit dielectric sphere (diameter, A = 1.3A) solutions for 10%, 20%, 30% and 40% volume fractions are shown in Figures 6.30 through 6.33 respectively. With the exception of the 10% volume fraction example, it can be seen that there is an excellent fit between the dielectric sphere solution and the averaged scattered fields obtained from the Monte-Carlo simulations. The slight differences seen in the 10% volume fraction example are likely caused by an insufficient quantity of independent samples for estimating the coherent scattered field from the random medium. The lower volume fraction poses a problem because of the relatively small size of the bounding sphere (compared to a wavelength) and the small number of scatterers involved in each realization (at 10%, N = 30) may not be sufficient to randomize the scattered field enough to ensure proper fading statistics (i.e. Rayleigh distribution) for the field magnitude. Thus, more realizations may be necessary to obtain convergence for this particular volume fraction. A plot of the differences between the Mie scattering solution and the averaged scattered fields from the random medium as a function of both the real and imaginary components of the homogeneous sphere's permittivity is given in Figure 6.34. It can be seen from this plot that there is a local minimum within the realistic range of real and imaginary permittivity for the volume fraction of 20%. Similar behavior is observed at the 1OO, 30%, and 40% volume fractions studied in this section. 6.4.4 Dependence on Volume Fraction and Particle Size The results of fitting the Mie solution for the homogeneous dielectric sphere to the coherent observed fields from the random medium sample at different volume fractions and particle sizes can be compiled into a single plot (similar to what was done in two dimensions) and the results compared to theoretical methods discussed in Chapter 3. In three dimensions, the scope of this comparison is currently limited to a narrow range of particle sizes and volume fractions due to computational considerations. In this context, it

151 20 20 '0 5 -a) 3 - 4-' Co L -20 -0 30 60 90 120 150 180 m 20 -- 0 0 0 ao " 0 -, -20 LL ( 0 30 60 90 120 Angle of observation (deg) 150 180 Figure 6.30: Comparison between the average scattered fields for a random medium with 10% volume fraction (solid line) and the best fit Mie solution for a homogeneous dielectric sphere (dashed line) with 6eff = 1.24 + iO.02.

152 _" 20 -a) 10. - CZ 0 - -10 U.0 30 60 90 120 150 180 CD r 20 2 a) 3 a) is -20 0 30 60 90 120 Angle of observation (deg) 150 180 Figure 6.31: Comparison between the average scattered fields for a random medium with 20% volume fraction (solid line) and the best fit Mie solution for a homogeneous dielectric sphere (dashed line) with ceff = 1.49 + iO.032.

153 2- 20 (1) m 20 1 10 rcd 2 0 iI -10 -0 co =M O S- 20 a) 10 (.0-. -10 ) 0 0 - 30 60 90 120 150 180 30 60 90 120 Angle of observation (deg) 150 180 Figure 6.32: Comparison between the average scattered fields for a random medium with 30% volume fraction (solid line) and the best fit Mie solution for a homogeneous dielectric sphere (dashed line) with 6eff = 1.73 + iO.044.

154 s 30 -D ' 20.-' 0)10 ~ 0 U- I ( 0 30 60 90 120 150 180 0 30 60 90 120 150 Angle of observation (deg) 180 Figure 6.33: Comparison between the average scattered fields for a random medium with 40% volume fraction (solid line) and the best fit Mie solution for a homogeneous dielectric sphere (dashed line) with ceff = 1.90 + iO.048.

155 0.1 0.2 0.3,.0.4 0.9 1.5 1.6 1.7 1.8 1.9 2 2.1 Re( eff) Figure 6.34: Contour and gray-scale plot of the differences between the Mie solution for a dielectric sphere and the averaged scattered field obtained from the random medium as a function of real and imaginary effective permittivity. In this example the volume fraction is 20% and the permittivity of the inclusions is c = 6.93 + iO.1. Dark areas represent areas where the difference is the least between the two solutions.

156 is possible to explore volume fractions ranging from 10% to 40% and particle sizes whose radius at 10 GHz varies from 3 mm to 3.5 mm (i.e. diameter varies from A/5 to A/4.3). This range is sufficient enough to explore the effects of increased scattering losses due to larger particle sizes and the behavior of both the real and imaginary components of effective permittivity as a function of particle density. In contrast to the results shown for two dimensions, the quantities of interest will be given in terms of the effective refractive index (neff = /eff) and the normalized extinction coefficient for the medium (see Section 3.3.2 on radiative transfer) which is given by r|e/k = 21m(e) | (6.4) The real part of neff is directly proportional to the phase delay of an electromagnetic field as it passes through the random medium and sIe/k is a measure of the power lost into the incoherent field. Figure 6.35 illustrates the numerically calculated results and provides a comparison to several of the theoretical models discussed in Chapter 3. The theoretical methods compared are the quasi-crystalline approximation with coherent potential (QCA, solid lines), the effective field approximation (EFA, dashed lines), and the Polder Van Santen mixing formula (PVS, dotted line). The first plot of Figure 6.35 compares the real component of neff, which is not expected to vary strongly as a function of particle size. This is a characteristic that is common between the theoretical and numerical method results. The three separate dashed lines in this plot indicate theoretical results from the effective field approximation for the three different particle sizes, with the effect of increasing particle size causing an increase in real refractive index. The numerical results obtained fall closely between the effective field approximation and the Polder-Van Santen mixing formula using a background dielectric of free space. The numerical calculations seem to indicate that the behavior of the fields near the inclusions, at these volume fractions, experience the full dielectric contrast between the inclusions and the host material rather than the dielectric contrast between the included material and the effective permittivity of the random medium. At some point, not modeled here, it is expected that this behavior should change as volume fraction increases and the medium becomes one dominated by inclusions rather than the spaces between inclusions. In this case, the background dielectric would no longer be similar to that of free space (as the numerical technique in this example would indicate) rather it would be closer to that of the effective permittivity, Ceff. The manner in which particles are arranged within the volume may have an effect on this behavior, as in this example, where the method of particle arrangement maximizes the average distance between inclusions.

157 The second plot of Figure 6.35 displays the normalized extinction coefficient given by (6.4) which reflects the power lost to the incoherent field. The numerical method is compared with the theoretical methods of the quasi-crystalline approximation with coherent potential and the effective field approximation (the supposed upper limit). Each of the method results is given in sets of three for each of the three particle sizes used (i.e. radius, a = 3, 3.25 and 3.5 mm at 10 GHz), with the larger particle sizes having the effect of increasing Ke/k for all three methods. It is evident from this illustration that at 10% volume fraction, that the numerical method and QCA-CP agree closely for the three particle sizes, while as particle density increases, the two methods diverge with QCA-CP predicting lower scattering losses than the numerical method. This likely reflects the natural limitation of QCA-CP to address particle interactions that occur in groups larger than two (i.e. QCA models multiple interactions between pairs of particles and excludes the effect three or greater particle interactions). This change occurs most predominantly for volume fractions of 25% and greater and is more pronounced for larger particle sizes. 6.4.5 Limitations of the Three-Dimensional Technique In the previous sub-section, the coherent numerical method presented in this chapter was used to compare with similar results obtained using theoretical methods for determining the extinction coefficient. In this manner, the numerical technique was used to explore the limitations of the other methods. The three-dimensional determination of effective permittivity however, has its own set of limitations, that are imposed due to computational concerns. The limits that are imposed are difficult to asses for a new technique such as the coherent method for determining effective permittivity because it is difficult to gauge when the technique is succeeding or failing. Figure 6.36 attempts to graphically display these limits as they have been explored thus far and to illustrate those points that have been demonstrated in the previous sub-section (indicated by x's). The lower limit of the bounded region is labeled as Nmax in reference to the fact that the direct T-matrix cannot simultaneously solve for more than 150 Rayleigh sized spheres on a Convex MPP-1000 parallel computer with approximately 500 mega-bytes of random access memory. Using less spheres would require the use of a smaller bounding sphere, thus decreasing the sensitivity of the coherent method to the imaginary component of effective permittivity (see Section 6.4.2). The volume fraction limit, labeled as "Packing Method" refers to an upper limit on the volume fraction that can be obtained using the random introduction of spheres into an

158 x 1 NUM C 1.8 -.: QCA 1.6- — EFA > 1.62.... PVS ' 1.4 - r 1.2 -0 0.1 0.2 0.3 0.4 0.5 Volume fraction 0.2 (D m ~X NUM 0.15- QCA - - * - EFA 0 - -. - N- -- * C 0.05- - - 0 Z 0 "'I 0 0.1 0.2 0.3 0.4 0.5 Volume fraction Figure 6.35: Real refractive index and normalized extinction coefficient comparison between the numerical method (symbols) and the theoretical techniques of quasi-crystalline approximation with coherent potential (solid lines) and the effective field approximation (dashed lines) for three particle radii ranging from 3 to 3.5 mm at 10 GHz. Normalized extinction (shown in the lower plot) consistently increases with particle size for all three methods shown.

159 Mmax Pmax 4mm Packing 0.1 0.2 0.3 0.4 0.5 volume fraction Figure 6.36: Illustration of the volume fraction and particle radius limits for the determina(A -3 cm). empty space (i.e. simulation of the Percus-Yevick pair distribution function for an ideal fluid). While other packing algorithms may be used to surpass this limit, the isotropy of number of independent samples. This would have the effect of increasing the number of not be representative of a homogeneous dielectric. The upper limit which dictates the maximum radius allowed for the inclusions is imposed because larger inclusions require more terms in the spherical wave expansion ('Mmax) addition to the fact that a large boundary dimension may be required to enclose a sufficient number of inclusions to accurately represent the random medium. This large boundary is limited by the maximum number of spherical waves (Pmax) allowed to represent the FThue.limits that have been discussed in this sub-se ction are a best gueslimits for the litmitations of the emethod that has been presented. It is critical to poresent them in some form to putit the technique, as it hase been adeveloped, into a greater context as a method for determining o effective permittivity for random media. In comparison to the technigque demonstrated in two dimensions, these limits are very constraining, but at present they are mostly due to the

160 computational technique that has been utilized and it is likely that features of the technique such as the shape of the confining boundary may be optimized in some way to make the technique more robust. Furthermore, the T-matrix algorithm in general has proved tricky to use in this context because the truncation of a spherical wave expansion implies an overall power loss for the system. This power loss has the potential of directly appearing in the imaginary component of effective permittivity, and thus would be the source of a larger imaginary component than is appropriate for the random medium under study (for the examples generated in the previous sub-sections, care was taken to avoid this possibility). For this reason, while it is still believed that the T-matrix is an appropriate technique for this application, it would be useful to back up the results that have been presented thus far with an alternative numerical technique such as the Finite Element Method or the Method of Moments in three dimensions. 6.4.6 Conclusions and Recommendations for the Three-Dimensional Technique In this section, the three-dimensional coherent field technique for determining effective permittivity for random media was developed. In the initial treatment of the method, it was shown that between 30 and 50 independent realizations of the random medium were required to converge on a representative coherent field. This coherent field can then be used to determine effective permittivity. A complete example was illustrated using a particle diameter of A/5 (radius, a = 3 mm at 10 GHz) for volume fractions ranging from 10% to 40% where it was shown that the average scattered field obtained from the Monte-Carlo simulations agreed well with a Mie series solution for a homogeneous dielectric sphere. Additionally, the differences between the Mie series solution and the averaged scattered field were shown to have an isolated local minimum within the physical range of ~eff for a representative volume fraction of 20%. This example was then expanded to include two other particle diameters of A/4.6 (a = 3.25 mm) and A/4.3 (a = 3.5 mm) where a comparison was made between the numerical method and the theoretical methods of the quasi-crystalline approximation with coherent potential and the effective field approximation. It was shown that the numerical method predicted a lower real component of the index of refraction than these theoretical models and predicted an extinction coefficient that agreed well with QCA-CP at low volume fractions and the smaller of the particle diameters. At higher volume fractions and large particle diameters, the numerical technique indicates larger scattering losses, a possible consequence of a higher order of particle interaction than QCA-CP theoretically accounts for

161 (i.e. two-particle interactions). A discussion was then provided regarding the current limitations of the coherent technique in three dimensions. These limitations are imposed by the computational technique that has been implemented (direct T-matrix approach) and may be broadened as the technique is used to explore a wider range of problems. It is recommended that the technique is expanded/verified in the future by (i) exploring a variety of different boundary shapes and sizes to maximize the sensitivity of the coherent field technique to the imaginary component of effective permittivity (the dielectric slab used in two dimensions seems to be one such possibility), (ii) implementing and comparing alternative numerical techniques for solving scattered fields and (iii) using the results that have been obtained in the this chapter to explore theoretical/numerical methods of reducing the computational burden that a full Monte-Carlo simulation that the problem entails. 6.5 Incoherent Method of Determining Extinction The coherent method for determining Eeff in this chapter can be contrasted against the incoherent method for determining the extinction coefficient, tK given by [Tsang et al., 1992; Chew et al., 1995]. The purpose of this section is to make the distinction between the two numerical methods clear and to emphasize potential difficulties in the application of the incoherent method. These difficulties can be summarized as follows: * The incoherent field is dependent on boundary shape and size because the random medium sample is isolated in a free-space background. * The conservation of energy approach used by the incoherent method cannot be applied when the included material exhibits absorption losses. * Due to the existence of a boundary, the nature of wave propagation changes. As the mean field trapped within the boundary generates more incoherent power than would be observed if the boundary were submerged in a homogeneous dielectric with permittivity of ceff. To begin, both methods examine the scattered field from a random collection of scatterers situated in free space (Figure 6.37). For a given incident field, numerical calculations are used to determine the scattered field which is composed of a coherent and incoherent component (the coherent component does not vary between each realization, n): -sca -incoh - coh -n En + E.(6.5)

162 4Eca Figure 6.37: Illustration of a volume current due to dielectric fluctuation in a random medium. Shown are some spherical components of the random medium. By performing many realizations, N, of the random arrangement of scatterers, the coherent field may be determined by averaging the scattered field -coh 1 sca (6.6) E ^ Enr i(6.6) n=l The incoherent field for each realization may be found by solving (6.5) or by calculating the sum of the squares of the scattered field in addition to the square of the coherent field n=l 1a -r 9ie -incohe2 e t te-hcaq. 2 1 N ~nsca 2 coh(~,)1 (2 6\ n=1 The coherent field is related to the size and shape of the confining boundary and arises from the random collection of scatterers behaving in part as a coherent structure with permittivity defined as ceff. The incoherent field arises from dielectric fluctuations between the inclusions and the background and it is this field which provides the basis for the incoherent technique. Following the incoherent technique, the extinction coefficient is found by summing the power scattered into the incoherent field. This is accomplished via the integration N =RJ (K,incoh 2) (6.8) in two dimensions and N = R2 ( Enc~h 2sin0 dsds (6.9) i J \' /6.9

163 for three dimensions. The R factor normalizes the scattered power with respect to the distance from the origin. The extinction coefficient is determined by Ke = aN/V, where V is the volume (or equivalently the area in two dimensions) containing the random medium. Thus the extinction coefficient is directly related to the incoherent field scattered from the random medium. Because the source and the observer are located outside of the random medium a potential problem with the incoherent technique arises. That is, the dielectric contrast between the random medium sample and that of free space has the potential of altering the behavior of fields from the way they would appear if the source and observer had existed in a medium whose permittivity was ~eff. Because the incoherent field is excited by the coherent field, the existence of a dielectric boundary on the incoherent fields would fundamentally alter the problem. The coherent method presented in this chapter circumvents this difficulty by directly taking into account the effective boundary of the random medium sample. To highlight the difficulties of the incoherent technique, either a theoretical or numerical approach can be taken. The theoretical aspect of this problem goes as follows. Fluctuations in the scattered field between different particle arrangement realizations are attributed to spatial variations between the permittivity of the inclusions/spaces, c(x, y, z), and that of the background, ceff. These fluctuations act as a volumetric current distribution, Jf, which can be concisely described by the wave equation as V x V x E1 - -2 = w2P(C - ceff)EI = iWUJf (6.10) where E1 is the electric field within the boundary containing the random medium and K is the effective propagation constant of the medium (as in Chapter 3). Two questions may be posed: (i) What is the contribution of a volume current Jf to the coherent field? and (ii) What is the contribution of a volume current Jf to the incoherent field? The answer to the first question is straightforward. By definition, the dielectric fluctuations do not contribute to the coherent field. The second question is addressed by relating the scattered field to the current fluctuation via the dyadic Green's function, G: Esca G-JfdV. (6.11) Here, G is the dyadic Green's function in the presence of a homogeneous dielectric scatterer with the prescribed boundary. For example, in the case of a spherical cluster of particles, the Green's function is ik, G = 47r E Cmn CnM(ko0R)M(^R) + DnN(koR)N(iR)J. (6.12) rm,n

164 350,,,,, o 300 - ~ 150 - 100 different — 0- ), and.5A )) and an inclusion where M and N are vector spherical wave functions (further details can be found in [Ruck o // et al. 1970]).. 100 --------- From the above discussion, it can bgle of observen that the scattered field is a function of bo(deg)th thFigure 6.38:effective peNormaittivity and tincoherent power as a fbounctiondary, as is expected. In contrast, three diffncoherent square boundary dimensionspace (0.75A(reen's function - -), and.5does not take into an inclusionnt these pefactors (a sparse medium assumption). This assumption rusions iscounter to the incoherent/10. method's original premis vector spherical wave funandom collctions (further details can behave fou in part as auck coherent structure with permittivity, ff.1970]). From the above discllustrasion, it can be seen that the incoscatterednt field is dependenta function thof boundary shape, the effetive permittivity and by (6.7) is plotted as a funidary, as is expected. In contrast, the 6.38 fincoherent method uses the fsquaree space Grey dimens functions rangingd does not take into 1.5account thesean infactors (a sparse medium assumptof = 3.6+i.1 for the two-dimensional runs counter to the incoherently shows that the incoherent field is a function of the boundary dimension; a small boundary dimension demonstrates a slowly varying function with respect to observation angle and for a large boundary dimension, the incoherent field is peaked in the forward direction, as would be expected. This set of plots can be contrasted with the coherent field observed for the same set of boundaries in Figure 6.39.

165.C 140 -C 120 // \ 0) ~/ - 100 / l 0 V 80 - N CZ 60 -z I 4~0 /. \.-' 20_>,..Y^ 0 60 120 180 240 300 360 Anale of observation (dea) Figure 6.39: Normalized coherent field as a function of observation angle for three different square boundary dimensions (0.75AA(-), A(- - -), and 1.5A-(... )) and an inclusion permittivity of c = 3.6 + iO.01. The diameter of the inclusions is A-/10. The similarity in the shape of the incoherent power (Figure 6.38) and the coherent field plots (Figure 6.39) is a strong indication that the contribution of the current sources to the incoherent field is a function of the boundary. This dependence was theoretically described by (6.11) and (6.12) and demonstrated by Figures 6.38 and 6.39. These results indicate that the incoherent energy generated from random fluctuations in the medium must pass through an effective boundary that occurs because the random medium is confined in a background of free space (i.e. there is a statistical likelihood that the incoherent energy will interact with one or more discontinuities at or near the confining boundary). Accurate implementation of the incoherent method would occur only if the random medium sample was located within a background medium of permittivity 6eff (i.e. the wavenumber outside of the particles should not be the same as free space, but the effective wavenumber). In this instance the incoherent field would not be a function of the observation angle. In Figure 6.40 it is demonstrated that the conservation of energy approach, that the incoherent method relies upon, is not appropriate for inclusions that demonstrate absorption losses. This figure illustrates the normalized extinction coefficient, he/k, calculated by the coherent method (-), the incoherent method ( — -), and the Polder-VanSanten mixing formula ( ) as a function of volume fraction for two different particle permittivities, c =

166 TM Polarization 0.05 a,0.04 o.0 0.03.0 0 0....0'' 1 0.02 - N E ' z 0 0.01 - '....u.u /.....)K........ 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Volume fraction Figure 6.40: Normalized extinction coefficient comparison between the coherent method (), the incoherent method (- -), and the Polder-VanSanten mixing formula (...) for two different inclusion permittivities: low-loss, E = 3.6 + i0.01(*) and lossy, c = 3.6 + iO.l(o). 3.6 + i0.01(*) and qi = 3.6 + iO.l(o). It is expected that the Polder-VanSanten mixing formula provides the lower limit for the extinction coefficient because it only accounts for losses due to absorption (inclusion diameter in this example is AX/10). For the lossy case, indicated by o's, the extinction coefficient calculated by the incoherent method becomes unphysically smaller than the extinction coefficient predicted by the mixing formula and is even less than the extinction coefficient predicted by the incoherent method for the low loss inclusions. This clearly indicates that the incoherent method is incapable of predicting the absorption loss in evaluating the extinction coefficient. It is interesting to note that by adding the absorption loss calculated by the Polder-VanSanten formula to the extinction calculated by the incoherent method, the result calculated by the coherent method can be reproduced approximately. However, this is not the case for the low loss example indicated by *'s in Figure 6.40. In this case the mean-field trapped within the boundary bounces many times which gives rise to an overestimation of extinction. In three dimensions, Figure 6.41 illustrates the behavior of the coherent technique (solid lines), the incoherent technique (dashed lines) and the quasi-crystalline approximation with coherent potential (dotted lines) for determining the normalized extinction coefficient for

167 volume fractions ranging from 10% to 40% and particle radii of 3.0 mm, 3.25 mm and 3.5 mm at 10 GHz. In this set of figures, the inclusion radii are converted to volumes because it is expected that the extinction is linearly proportional to volume for small variations in inclusion radius. At the low volume fractions of 10% and 20% (bottom two graphs of Figure 6.41) very good agreement between the QCA-CP, the incoherent method and the coherent method of determining he/kk can be seen. At these volume fractions the coherent and incoherent method are expected to agree closely because the inclusion density is relatively sparse and therefore the boundary between the random medium sample and free space has a low dielectric contrast and is effectively unseen. At the low volume fractions of 10% and 20% the numerical methods agree quite well with the theoretical method of QCA-CP. This is likely due to the fact that the scattering mechanism of the random medium is dominated by single- and two-particle interactions in this relatively sparse arrangement. At the higher volume fractions of 30% and 40% (upper two plots of Figure 6.41), the behavior of the three methods with respect to one another diverges. Error bars are provided for the coherent method to illustrate the range of Im(Ceff) which the coherent method gives satisfactory results. This is a result of the coherent method relying on precise values of the scattered fields, and thus the exact value obtained for Ke/k may be a function of the fitting method used. In contrast, the incoherent method is an integrated quantity as reflected in (6.9) and therefore provides a single, unique answer that varies slightly as a function of boundary diameter and the number of independent realizations. These variations are within 10% of the values shown if Figure 6.41. While the uncertainty for the coherent method in comparison to the incoherent method of determining Ke/k is significant, the general differences between two methods is revealing. The coherent method consistently provides results that indicate scattering losses are greater than those calculated by the incoherent method. This is consistent with the previous discussion regarding possible limitations of the incoherent method when absorption losses are present in the included material (in this case, ci = 6.93 + iO.1). The difference between the two methods is directly due to the phenomenon of absorptive losses, a process not accounted for by the incoherent method. In this section, the similarities and differences between the coherent and incoherent method of determining the extinction coefficient for a random medium were outlined. It was demonstrated both theoretically and numerically that the incoherent method may have shortcomings for the three reasons that (i) the incoherent field is dependent on the boundary

168 0.06 0.05 0.04 a) I la) 0o c o o4 C - o. - x a, c N E 0 z 0.06 0.05 0.04 0.03 110 120 130 140 150 160 170 180 0.05 0.04 0.03 0.02 1 0.03 0.02 0.01 1 20% v.f. I..I:,.... 10 120 130 140 150 160 170 180 I I I I I I I 10 120 130 140 150 160 170 Inclusion volume (mmA3) 180 Figure 6.41: Comparison between the coherent (solid lines) and incoherent (dashed lines) numerical techniques for determning normalized extinction coefficient oe/k for different volume fractions and inclusion volumes at 10 GHz. The three volumes used in the comparison relate to particle radii of 3, 3.25 and 3.5 mm. The dotted line indicates comparable results obtained from the quasi-crystalline approximation with coherent potential.

169 shape and size because the random medium sample is isolated in a free-space background, (ii) the conservation of energy approach used by the incoherent method cannot be applied when the included material exhibits absorption losses and (iii) due to the existence of a boundary, the nature of wave propagation in the random medium changes. For these reasons, it is thought that the incoherent method of determining the extinction coefficient may in general provide unreliable results. 6.6 Summary and Conclusions In this chapter a Monte-Carlo numerical technique for determining effective permittivity based on the coherent scattered fields from bounded realisations of random media for both two- and three-dimensional problems was developed. In two dimensions the method of moments was used to determine effective permittivity for a wide variety of particle sizes and volume fractions. These results were compared with calculations provided by the quasi-crystalline approximation, the effective field approximation and the Polder-VanSanten mixing formula for TM and TE polarized fields. Based upon this comparison, it was possible to dilineate regions of applicability for the theoretical techniques. In three dimensions, similar calculations were performed by utilizing the T-matrix method for determining the scattered fields from a collection of spheres. It was shown that the recursive version of this algorithm (RATMA) could not be used in this treatement, and thus it was not possible to explore as wide range of situations such as was done with the two-dimensional technique. For a small range of particle sizes and volume fractions, it was possible to show that the coherent technique provided meaningful results. These results were used to compare with the theoretical methods of the quasi-crsytalline approximation and the effective field approximation. The treatment of the three dimensional technique ended with a discussion of the limitations of the technique and gave suggested directions for future research. The final section of this chapter was devoted to a comparison between the two alternative numerical techniques of the coherent and incoherent methods for determining the extinction coefficient for random media. Limitations of the incoherent method were discussed and demonstrated with numerical examples in both two and three dimensions. It was concluded that the incoherent method may not provide accurate results when the included material exhibits absorption losses or when the volume fraction becomes larger than 30%. This chapter represents the cummulative efforts of the previous chapters of this disser

170 tation. The objective has been to provide a reliable, consistent method for determining the fundamental parameter of effective permittivity (the behavior of coherent fields) for random media. In general, the two-dimensional implementation of this technique has proved to be applicable to a wider range of problems than its three-dimensional counterpart. In both sets of data however, it was possible to compare results obtained by the numerical method with other, established theoretical methods in common use. The impact has been to provide an independent method through which these theoretical methods can be gauged or evaluated and also provides a basis that the numerical method can be used to explore other random media problems of interest.

CHAPTER 7 Scattering from Wheat Grain Heads 7.1 Introduction Up to this stage, a very broad basis has been established to address theoretical problems in wave propagation encountered in remote sensing of random media. The tools and techniques employed to this end can be applied to a number of more simple problems that occur in electromagnetic analysis. This chapter addresses one such problem that has been highlighted by a coherent wheat scattering model at L-, C-, and X-band frequencies [Stiles, 1995]. In his dissertation, Dr. Stiles noted that a first order coherent scattering model of a wheat field could accurately predict experimental observations at these frequencies, that is, until the plant had reached full maturity (i.e. the grain head of the wheat plant became developed). In the first order coherent model, a circular cylinder was used to model the scattering properties of the grain head. It is believed, particularly at C- and X-bands that the cylinder model is not appropriate to correctly model the scattering behavior of the grain head. The material covered in this chapter addresses this problem by combining a mcre elaborate physical model of the grain head with the T-matrix electromagnetic solution technique to improve the fidelity of the grain head scattering model. Because it is impractical in a general sense to use this numerical technique in a complete analysis of the scattering from a wheat field, a semi-empirical model is developed based on the numerical scattering model. Thus, the grain head model fidelity is increased at no cost to the model utility. The approach to obtain the electromagnetic model follows a five step process that can be applied to a range of problems. The approach, as it applies to the specific problem of scattering from a grain head goes as follows: 1. The zero order scattering terms of a simple geometric model are derived. In the case of a grain head, a dielectric cylinder is used and the dependencies on the observing parameters of the radar (frequency, polarization and angle of incidence) and the physical 171

172 grain head parameters (moisture content, length and width) are made specific. 2. The full order scattering cross section of a dielectric cylinder is computed numerically for a variety of radar and grain head parameters. These calculations will be used for comparison with the more elaborate, high fidelity grain head model. 3. An elaborate physical model is constructed for the grain head structure using components that can be analyzed with an electromagnetic numerical code. In this case, individual grains are modeled as dielectric ellipsoids which can be ordered into a full grain head structure. 4. A high fidelity numerical electromagnetics code is used to analyze the scattering behavior of the physical structure. For the grain head problem, a T-matrix algorithm is used. 5. The functional forms of the zero order scattering terms calculated in the first step of this process are used as basis functions for a semi-empirical model that is fit to the high fidelity numerical calculations obtained in step four. By relying on the basic functional forms, the dependencies of the semi-empirical model on parameters such as moisture content and observation angle are more robust than would be the case if simple polynomials were used to perform the fit. This approach is developed in detail here so as to provide a usable model for the wheat grain head and to demonstrate a methodology for future problems in need of similar rigor. 7.2 Theoretical Construction The goal of this section is to relate the physically observable parameters of a wheat grain head to an electromagnetic model that may be used in a full scattering analysis for a wheat field. Figure 7.1 displays a photograph of a typical wheat grain head along with two electromagnetic models that will be explored here. These are the truncated dielectric cylinder model and the wheat grain model (where contributions from individual grains are taken into account). The overall goal of this chapter will be to relate the scattered field to the incident field for a wheat grain head, with known physical parameters, via the scattering matrix defined by F ESc eikr S S, i Eic 1 E kr Shv Shh L E (7.1) F~ vShEfC in

173 60 ------------- 50 Ko 1X D 450 '- -- 35C-'OR 5 0 5 Figure 7.1: Two models for the grain head of wheat. Individual grains may contribute more strongly to the backscatter field than the "smooth" cylinder model most commonly employed. the context of a first order scattering model. The incident field in this problem is defined by Ei - — eik (7.2) and ah -- for horizontal polarization. Here, the incident field is propagating in the -10xs) -5 0 5 Figure 7.1: Two models for the grain head of wheat. Individual grains may contribute more strongly to the backscatter field than the "smooth" cylinder model most commonly employed. In this treatment, the cross polarized scattering terms, Sbh and Shv are not of concern infor a the context of a first order scattering model. The incident field in this problem is defined by Et = aeikk' * (7.2) where the pilot vector', a^ may be broken up into a, = x cos 0 + z sin 0 for vertical polarization and ah = -y for horizontal polarization. Here, the incident field is propagating in the positive x and negative z directions so that ki = x sin 8 - z cos 0 where 0 is the incidence angle measured from nadir (the z-axis). In general, three components of the grain head scattering behavior are of interest for a first order coherent wheat field scattering model. These components are (i) direct backscatter, (ii) specular scatter (for the ground bounce term), and (iii) forward scatter (or the imaginary part of the forward scatter to be used for determining the extinction coefficient).

174 incident field V direct backscatter ground bounce ^ extinction Figure 7.2: Schematic diagram of the scattering components of interest for the wheat grain head. Note that the backscatter from the grain layer is unattenuated by other layers in the vegetation. These directions are illustrated in relation to the grain head by Figure 7.2 and are given by backy back -sin 0 + z cos0 (7.3) forw = x sin 0 + cos (7.4) kspec = -x sin 0 - cos 0 (7.5) respectively. Given the physical parameters of the grain head length L, cross-sectional area A (or equivalently diameter D), volumetric moisture content My, and volume fraction f, it is desired to calculate the co-polarized terms of the scattering matrix for the three principle scattering directions described above. 7.2.1 Truncated Cylinder Model Of these two models that will be addressed, the truncated cylinder model is the most simple and straight forward (Figure 7.1). In the low frequency limit, the scattering matrix can be determined by approximating the currents on a truncated cylinder with those of an infinite cylinder [Sarabandi and Senior 1990]. Integration over the finite length of the actual cylinder introduces a sin U/U or sincU function into the scattered field, where U is related

175 to the electrical length of the cylinder. Specifically, the scattering matrix is given by S = -27r A3i k, x ks x P Pa] sincU (7.6) and U = -(cos9 s- cos 9). (7.7) For an infinitely long cylinder, sincU = 5(s-0i) and the solution given by (7.6) is exact. For a truncated cylinder, the currents at either end of the cylinder are not correctly terminated, thus the accuracy of the truncated cylinder model is best in the broadside incidence and scattering directions, and becomes less accurate towards the axis of the cylinder (where the effect of the end caps will eventually dominate). In (7.6), the polarizability tensor P relates the incident field to induced currents within the cylinder at low frequencies. This is a diagonal tensor whose elements for a circular cross section are given by PxX = Pyy = 2(e - 1)/(e + 1) and Pzz = - 1 (see [Sarabandi and Senior, 1990] for details). Substituting values for the incident and scattered field directions for both v and h polarizations yields the co-polarized scattering matrix components for the low frequency model. These are given by sbk — 7r2LA 2 ] Sback = 27[2P — Pzz sin2 0 + P~ cos2 0] sinc [kL cos ] (7.8) Sck = 27T2 LAPyysinc [kL cos ] (7.9) oA3 sforw = 27r2LA [Pzz sin2 0 + P cos2 ] (7.10) srw = -22LAyy (7.11) T0 spec = 27r2LA Pzz sin20-Pcos2 (7.12) Sse = 2r A (7.13) In all cases, the magnitude of the scattering matrix is proportional to the product of the length, L, and cross-sectional area, A, of the cylinder. The elements of the polarizability tensor are solely dependent upon the permittivity of the cylinder. For vegetation type media, the dielectric constant can be related to the volumetric moisture content using formulas developed by [El Rayes and Uby 1987]. The volumetric moisture content is determined by first physically weighing grain head samples before and after drying to determine the total water content of the grain head and then dividing by the volume of the grain head. In this case, the volume is determined by the

176 date MvYl Mgrain" diameter (cm) height (cm) 6/4 0.19 0.38 0.79 6.8 6/10 0.11 0.22 0.97 7.2 6/15 0.08 0.16 1.05 7.4 6/18 0.12 0.24 1.05 7.4 6/22 0.16 0.32 1.05 7.4 6/25 0.16 0.32 1.05 7.4 6/29 0.16 0.32 1.05 7.4 7/2 0.16 0.32 1.05 7.4 7/7 0.14 0.28 1.05 7.4 7/15 0.14 0.28 1.05 7.4 Table 7.1: Measured physical parameters of wheat grain heads on different dates [Stiles, 1995]. boundaries of the cylinder model. Alternatively, a volumetric moisture content can be calculated for the individual grains making up the grain head, and the dielectric constant of the complete grain head related to the dielectric constant of the grain heads by the Polder-VanSanten mixing formula discussed in Chapter 3. In the treatment that follows, the cylinder permittivity is calculated by the former of these two approaches. Measured physical parameters from a wheat field under study have been performed by [Stiles 1995]. The volumetric moisture content of the grain head and the individual grains (assuming a 50% volume fraction) and the physical dimensions of the grain head, are given in Table 7.1 over the period of one full growing season in the state of Michigan. The range of dielectric values dependent on the volumetric moisture content calculated by the formulas developed by [El-Rayes and Ulaby 1987] are summarized in Table 7.2 at L-, C- and X-band frequencies for both the dielectric cylinder and the grains that compose the complete grain head. The low frequency version of the truncated cylinder model just developed is in general not appropriate for modeling the scattering from a grain head at the higher frequencies of C- and X-band. For this reason, it is necessary to use more terms in the Bessel and Hankel function expansions to represent the incident and scattered fields. The complete expansion for an infinite length dielectric cylinder can be found in [Ruck et al. 1970, pg. 205; Ulaby and Elachi 1990, pg. 92]. The low frequency model described above consists of the zeroth

177 Mv L-Band C-Band X-Band min wheat 0.08 5.8 + il.7 4.7 + il.2 4.3 + i1.2 grain 0.16 11.4 + i4.0 8.8 + i2.8 8.0 + i2.8 mean wheat 0.14 10.0+ i3.4 7.8 + i2.4 7.1 + i2.4 grain 0.28 18.8 + i6.3 15.2 + i4.9 13.7 + i5.3 max wheat 0.19 13.2 + i4.5 10.4 + i3.3 9.3 + i3.4 grain 0.38 25.6 + i8.3 21.4 + i6.7 19.2 + i7.7 Table 7.2: Dielectric values for the wheat grain head and wheat grains for varying moisture content. The dielectric model is based upon [El-Rayes and Ulaby, 1987]. order terms of this expansion and will be used to aid in the development of a semi-empirical model. This model is based upon calculations performed using the T-matrix algorithm in conjunction with a physical model for the arrangement of individual grains making up the complete grain head. The truncated cylinder model used to represent the grain head is expected to be deficient for four reasons: (i) the dielectric cylinder model is based upon the solution for an infinite length cylinder, thus away from the broadside direction the model becomes inaccurate, (ii) a unique equivalent dielectric constant for the cylinder cannot be obtained due to the small number of particles that constitute the cylinder, (iii) the permittivity of the grain head is anisotropic due to the vertical alignment of the individual grains (modeled as ellipsoids), and (iv) the surface of the grain head is not smooth and thus it is expected that backscatter away from broadside incidence would be underestimated by the cylinder model. 7.2.2 Physical Grain Head Model While the differences between the grain head and cylinder model may appear at first glance to be second order, especially since the total biomass volume occupied by the grain heads is relatively small, in a coherent model, the effect of the grain heads dominates the observed backscatter. This can be attributed to two reasons (i) illumination of and direct backscatter from the layer of grain heads is unattenuated by other components of the vegetation layer, and (ii) all the remaining backscatter components of the vegetation layer are attenuated by the the grain head layer. Differences in the scattering components of the grain head of a few dB may have a very strong effect on the predicted backscatter from a

178 mature wheat field. This concern is especially relevant for observations at incidence angles near the nadir direction where the cylinder model of the grain heads itself is of questionable accuracy. Physical inspection of a wheat grain head reveals that it can be considered to be composed of a collection of ellipsoidal grains centered around a vertical axis (Figure 7.1). A physical model of the wheat head can be constructed by considering the center of each ellipsoid to be positioned along one of four helical coils that wrap around this vertical axis at intervals of 90 degrees (the number of coils may be adjusted for alternate grain models that need to be studied). The pitch of the coil and number of grains per unit length may be directly related to the physical dimension of individual grains (measured to have semi-axis lengths of 3.2mm, 1.5mm and 1.5mm). Thus, given the average dimensions of individual grains, a physical model may be created as shown in Figure 7.1. 7.2.3 Electromagnetic Analysis of the Physical Grain Head Model The recursive aggregate T-matrix algorithm was used to determine the principle scattered field components for the three frequencies (L-,C-, and X-band), the range of grain permittivities (obtained from Table 7.2) extending from 5 + iO.O to 25 + i10.0, grain height from 6-8cm and nominal grain diameter of 1.2cm. A representative example of these results are shown in Figure 7.3. For each of the three frequency bands, the back-scatter, forward-scatter and specular-scatter is shown for both vv- and hh-polarizations along with the results obtained from the cylinder model. Inspection of this set of figures indicates that * The general shape of the scattering functions are very similar between the cylinder and grain head models for all three frequencies and all three scattering components. * Agreement between the cylinder and grain head model is best at broadside incidence and deteriorates towards the nadir direction. This behavior is expected because the cylinder model itself becomes inaccurate at these incidence angles. * Back-scatter nulls are absent or minimized in the grain head model at C- and X-band. This is due to the "volume scattering" nature of the grain head model. * At nadir incidence, both polarizations at all three frequencies saturate to a minimum value. At these angles, the top of the grain head is contributing to the backscatter, an effect not accounted for in the dielectric cylinder model of the grain head.

CDt t-11 II o. CD -KJ CD 0 CD 0 0 0 CD CD S 0 CD CD CD 0 CD vv-pol (T-matrix) ihh-pol (T-mvatrix) \7V...IM,)0le~lnir L-Band C-Band X-Band = -40 I-o +- -60 C-3C cn 2 -40 0 -o -50 -60 "A IA 0 -20 -40 II 20 0 -20 -40 \/ \, ( N V "'' 0 10 20 30 40 50 60 70 80 90 -60' 0 10 20 30 40 50 60 70 80 9( 1 (. 0 10 20 30 40 50 60 70 80 90 0 -......l..1. — ' I......... I -.. ---------------. - I I I ]1 I --- —--— L. -- 20 10 0 20 10 0 - -.1 - --....... 11 -. -............... 0 10 20 30 40 50 60 70 80 90 I ki )10 20 30 40 50 60 70 80 90 ) 10 20 30 40 50 60 70 80 9( -50.. - -.-...... ---06 07 0 10 20 30 40 5060 70 80 90 angle from broadside (deg.) I (. 0 ---------------- -10k II....... -,-,- -...:. I -20' -10, 0 10 20 30 40 50 60 70 80 90 angle from broadside (deg.) 0 10 20 30 40 50 60 70 80 90 angle from broadside (deg.)

180 At the initial time that this model was being developed, the limitation with RATMA discussed in Chapter 5 had not yet been fully explored. Instead, for each frequency and grain head length, the value of Pmax was determined by increasing this value until the scattered field patterns converged. Once additional insight about the limitations of the recursive aggregate T-matrix algorithm was further understood, the direct T-matrix method was used to analyze the grain head for a typical length and volumetric moisture content. Shown in Figure 7.4 is a comparison between the RATMA derived result and the direct T-matrix result. As can be seen from the figure, although minor differences between the two algorithms occur in the backscatter pattern, for the most part the agreement is very good. The good agreement between the two forms of the T-matrix algorithm can be attributed to three reasons: (i) a large number of terms in the spherical wave expansion was used, (ii) dielectric losses in this problem are high with respect to the real component, and thus multiple interactions are not a strong factor in the overall pattern, and (iii) only the principle scattering directions are investigated which implies that most of the energy in these directions is from single grain interactions with the incident field. 7.3 A Semi-Empirical Grain Head Model Up to this point, the details of developing a more accurate grain head model have been presented. The accuracy obtained has been at the expense of utility however (i.e. the T-matrix algorithm is complicated, cumbersome and relatively slow). In this section, the scattering calculations obtained by the T-matrix algorithm are used as the driving component behind a semi-empirical model whose functionality is based upon the low frequency truncated cylinder model described in Section 7.2.1. This is not to say that the dielectric cylinder model will be reconstructed, rather, the functional form of the truncated cylinder model will provide basis functions for which to construct the semi-empirical model. By relying on the same machinery that makes up the dielectric cylinder model, but tailoring the fit to match the T-matrix grain head calculations, the overall fit as a function of grain permittivity (volumetric moisture), length and cross section will be much more stable. Additionally, to make the task of creating the semi-empirical model more simple, only the critical incidence angles ranging from 20 to 70 degrees angle from broadside are matched. This relaxes restrictions on the fit of the semi-empirical model while emphasizing the conditions under which observations are most likely to be made. The validity regions and the scattering components that make up the semi-empirical model are illustrated in Figure 7.5.

181 Grain head direct T-matrix vs. RATMA (C-Band) 0 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 80 90 angle from broadside (deg) Figure 7.4: Comparison between RATMA (solid) and the direct T-matrix algorithm (dashed) method of determining the scattered from a wheat grain head at C-Band. Shown is a typical example of a 7cm long head with volumetric moisture, M6 = 14%, consisting of 35 grains whose dimensions are 6.4 x 3.0 x 3.0mm.

182 Figure 7.5: Scattering terms of interest for the wheat grain head. Shown are the validity regions most critical for the semi-empirical model. The scattered field for the principle components of back-, forward-, and specular-scatter are a minimum at nadir incidence (0 = 0) and maximum at broadside (0 = 7r/2). This general behavior should be well represented by a Fourier sine series given by 00 Esca = a2sin2nO for 0 < 0 < 7/2 (7.14) n=0 The coefficients a2n may vary as a function of grain head length, cross-sectional diameter, and grain permittivity (volumetric moisture). Based on the low frequency analysis for a truncated dielectric cylinder, it is expected that the dependence on the physical dimensions of length and cross-sectional area will be linear. The dependence on permittivity, however, may be quite complicated. To overcome this difficulty, we may use our understanding of electromagnetic principles to render this dependence in terms of a Taylor series that utilizes the polarizability tensor at normal incidence for the dielectric cylinder as the basis function. The co-polarized components of the scattering matrix for the semi-empirical model are given by the equation (derived by using (7.6) and (7.7)) LA Spp = 3 [ao + a2cos(20)] K'p (7.15) 0 an = eClnPpp(C) + c2nPpp(C) + e3nPpp() In (7.15), L and A are the grain head length and cross sectional area respectively, p is the polarization, K7 is a function that is dependent on the three different types of scatter

183 ing mechanisms being modeled (m=back for backscatter, m=forw for forward-scatter and m=spec for specular-scatter) and Ppp(E) is the functional form of the polarizability tensor at normal incidence defined by P (e-I) or Phh = (- 1)/(+ 1) (7.16) For the backscatter direction, back = a + sin 2L cos ) (7.1]7) v = Ao0/.1 (7.18) Ahh = A0 (7.19) where the electrical length of the cylinder is observed to be polarization dependent due to an anisotropic permittivity that reflects the preferential alignment of grains in the vertical direction. The backscatter field model contains an additional term, a;, that is due to the incoherent contribution of energy to the total field by individual grains. This is because the electric fields scattered from individual grains are not likely to cancel out one another. Thus, the deep nulls that are observed in the cylinder model will not be present at C- and X-band grain model. Furthermore, the incoherent term shows up only in the backscatter direction because in this direction, the signal is the most sensitive to the phenomenon of volume scattering. So while it is necessary to maintain an additive incoherent term here, it is eliminated in the forward- and specular-scatter directions. In the forward-scatter direction, the coefficient, Kf, is a relatively weak function with respect to the imaginary component of grain permittivity. This is due to the fact that the grain head is relatively thin in cross section, and therefore the magnitude of the signal that propagates through any one grain head suffers only minor absorption losses. The forward-scatter coefficient is found to be -forw = 1 + 0.022e-018c'(E - 2 (7.20) K forw 1 + 0.012-0186'(C/)2 (7.2 1) hh 1+ 0.012 -For the specular direction, a "Brewster angle" absorption effect is observed for the vertical polarization. This is related to the two components of the bracketed term in (7.12) summing to zero. For a low frequency dielectric cylinder, the Brewster angle is given by tan0 = /(c + 1)/2 whereas, for a dielectric half space, this angle is tan = AV. The

184 Brewster angle for the wheat grain model is likely to fall at an intermediate value between these two angles. By using the functional form of the dielectric half space however, it is possible to mimic the behavior of the scattered field using the vertical polarization reflection coefficient given by R,, R rew -COS - 1 - sin2 /ebrew /brewcos 0 + /1- sin2 0/brew where Cbrew is the permittivity related to the Brewster angle observed for vertically polarized scattered fields. This angle can be determined for each L-, C-, and X-band frequencies for the range of grain permittivities under study and then related to the Brewster angle permittivity, ebrew. By examining the behavior of Ebrew of the parameter space of c (or equivalently Mv), a relation can be derived between (brew and c. The linear relation that can be used to model this dependence is (brew = bo + bi(E - 1) where the coefficients bo and bl are given by L C X b0 1.86 1.60 1.97 b1 0.10 0.14 0.066 It should be noted that while this model provides an excellent fit for e' between 5 and 25, it does not approach the correct limiting value of (brew = 1 when c = 1. Thus the dielectric portion of this model should not be employed for grain permittivities below 5 (at these permittivities, the grain heads would be uncharacteristically dry and would not be an occurrence found in nature). Thus, for the specular component of the scattered field, Kspec _ Rvv and shec 1. VV hh A final comparison of the semi-empirical model, the dielectric cylinder model and the numerical grain head model is given in Figure 7.6 for typical values of grain head length, diameter and moisture content. Overall, it can be seen that the semi-empirical model agrees very well with the numerically derived model. This agreement can best be illustrated by comparing the results directly between the semi-empirical model and the numerical model over the entire range for which it was constructed (Figure 7.6 illustrates the agreement for one particular permittivity). Figures 7.7 through 7.9 illustrate this for L-, C- and Xbands, for varying real and imaginary permittivity such that 5 < c' < 25,0.0 < c" < 10.0 and observation angle varying from 20~ < 0 < 70~ for a typical grain head length and diameter (L = 7cm andD = l.lcm). For each frequency this gives a total number of 600 points. As can be seen from the set of figures, while the overall performance of the semiempirical model is good, this performance degrades as frequency is increased. This behavior

185 is expected because the volume currents of the grain head can no longer be represented by the the first few terms of a cosine expansion. Increasing the order of this expansion and the order of dependence on the polarizability tensor would undoubtably increase the accuracy at higher frequencies, but the model would become much more complex. 7.4 Extinction Coefficient The optical theorem relates the extinction coefficient (proportional to the amount of energy a signal loses per unit volume) to the imaginary component of the electric field [Tsang et al., 1986] by PP = 47rn0 (Im {Spp}) (7.22) where no is the density of the scatterers whose scattering matrix, S is known. In this case, that scattering matrix is for the grain heads. Because the extinction coefficient is dependent on the imaginary component of the scattered field, rather than its magnitude, the angular dependence modeled previously by (7.14) requires more terms to track both the magnitude and phase. Rather than using the cosine expansion, a first order angular dependence can be determined by deriving the first term in the exact series solution for oblique incidence on a dielectric cylinder given by [Ruck et al., 1970, pg. 266]. Following this method of analysis, the sought after angular dependence of the imaginary component of the scattered field is given by Im{S C} = iao+ai cos52 (7.23) LA Im {hh} = 3 ao (7.24) to The dependence of the coefficients ao and al on the real and imaginary components of grain permittivity may also be determined through a similar analysis. Overall, the dependence of the imaginary component of the scattered field on incidence angle and grain permittivity is not as straightforward a function as the magnitude functions that were laid out in the previous section. Nevertheless, an overall good fit was found by making small adjustments in the derived dependencies mentioned above. A consequence of making these adjustments however is that the functional form of the coefficients ao and al changes depending on the

t 11 4t —4 I I — A n 5 PO t r-)t 11 n 5 0 CD CD Cf2 CD 0 CD x CD 0 0 rC'D 0 -Sn 0 0 0 CDt Cr) CD L-Band C-Band X-Band = -40 a. -50 c~c -60 CA) o -70 Asemi-emperical model......... vv-pol (T-matnx) -hh-poi (T-nJidtrix) -SM I 1.. 0 10 20 30 40 50 60 70 80 90 = -30 &0 un Z -5-s.......... 0 -20 -40 -601 10 5 0 -5 -10 -10 0 7"K 7 —. IX - I 10 20 30 40 50 60 70 80 90 I L 20 0 -20 -40,.x-x~X-)~- )(2 )K-x;.K I \ 7 .XX >K -- - V '.. V i -60L 0 15 -10 - 5 - 0 - 0 - - - X.-.. I'll 1....... X X, X,,, X>S — - -. -- - - ----- i 10 20 30 40 50 60 70 80 90 I L 3 10 20 30 40 50 60 70 80 90 ) 10 20 30 40 50 60 70 80 90 10 20 30 40 50 60 70 80 90 = -20 Cr) *- -60 03 -80 X "K> --- —-. '. xxx4 >">(X.. -10F -----. )............> --- 20 10 0.. x... -20 -30, r, I I I 3 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 80 90 -10U 70 809 0 10 20 30 40 50 60 7 09 n angle from broadside (deg.) CD 0 CD

187 -35 Z -40 o -45 E, -50 ' -55 6 -M),* -60 -65 F -65 1[ II -70 -70 -60 -50 -40 Semi-empirical model (dB) Figure 7.7: Agreement between the semi-emperical model and numerical model for grain head scattering at L-band. Shown are 600 points for back-, forward-, and specularscattering over the full range oof incidence angles (20 < < 70) and grain permittivities (5 < c' < 25, 0 < " < 10). 10r "0 "0 a) O x 7.'; E 0 -10 k:gK.4 " "A -20 -30 -40, I... i -40 -30 -20 -10 0 Semi-empirical model (dB) 10 Figure 7.8: Agreement between the semi-emperical model and numerical model for grain head scattering at C-band. See Figure 7.7 for details.

188 20 2 10 E 0. -20 a- - *30 * -40 ' -40 -20 0 20 Semi-empirical model (dB) Figure 7.9: Agreement between the semi-emperical model and numerical model for grain head scattering at X-band. See Figure 7.7 for details. observing frequency. These coefficients, which are used in (7.23) and (7.24) are L-Band ao = 4.5 x 10-3 + 7.31m {Phh( + 1)} a = 6" 1.7 - 0.27V/) + 0.00144(c')1 5 d ao - 0.14 + 0.Oll0 ' + 10.31m {Phh(6 + 2)} C-Band (7.26) a =- 0.8c' - 2.7 + C" (2.1 - 0.37 ) X-Band ao = 0.47 + 0.04c' + 17.61m {Phh(C + 4)} a = 7.3 log(c') - 8.8 + c"(0.67 - 0.026c') The performance of the semi-empirical extinction model given by (7.23) through (7.27), is demonstrated at the L-, C- and X- band frequencies in Figures 7.10 through 7.12 for typical values of grain head length, diameter and moisture content. In this same set of figures, the equivalent results obtained from the dielectric cylinder model are also shown as dashed lines, where it is seen that the dielectric cylinder grossly overestimates the extinction caused by the grain heads. This effect may possibly be due to the reduced volume that the grain head occupies in relation to a filled cylinder. Further accuracy of the semiempirical model may be obtained by retaining more cosine dependent terms derived from the dielectric cylinder expansion; this approach would give a significantly better match between the numerically derived extinction and the semi-empirical model for the horizontal component of the extinction coefficient.

189 Extinction cross section 0 0 10 20 30 40 50 60 70 Angle from broadside (deg) 80 90 Figure 7.10: Imaginary component of Spp at L-Band to be used in the calculation of extinction coefficient. This figure shows results for the grain head model (solid), cylinder model (dashed) and semi-empirical model (x); vv- and hh-polarizations. In this example, L = 7cm, A = 1.3cm2, c = 20 + i6.25(M, = 14%). Extinction cross section 30 20 10 0 -10 -20 L 0 10 20 30 40 50 60 70 Angle from broadside (deg) 80 90 Figure 7.11: Imaginary component of Spp at C-Band to be used in the calculation of extinction coefficient. This figure shows results for the grain head model (solid), cylinder model (dashed) and semi-empirical model (x); vv- and hh-polarizations. In this example, L = 7cm, A = 1.3cm2, c = 15 + i5.0(Mv = 14%).

190 Extinction cross section 40 30 - 401 ---. --- —---------- 20 10 0 10 20 30 40 50 60 70 80 90 Angle from broadside (deg) Figure 7.12: Imaginary component of Spp at X-Band to be used in the calculation of extinction coefficient. This figure shows results for the grain head model (solid), cylinder model (dashed) and semi-empirical model (x); vv- and hh-polarizations. In this example, L = 7cm, A = 1.3cm2, = 15 + i5.0(Mv = 14%). 7.5 Conclusion Chapters previous to this one have addressed a variety of experimental, theoretical and numerical issues generally associated with the science of remote sensing. In this chapter, the machinery and understanding developed in the previous chapters was put to use to address a specific problem for a first order coherent scattering model of a wheat field. Specifically, that problem was to determine an efficient, accurate method of calculating the principle scattering components (i.e. back-, forward-, and specular-scatter) from a wheat grain head. A general methodology was presented that first began with the physical modeling of grains on the wheat head, followed by a full numerical calculation of the scattering from this dielectric structure. To make the model usable, a semi-empirical model was developed based on the numerical calculations. The semi-empirical model relied on physical knowledge of the interaction of electromagnetic waves with a dielectric cylinder. Thus, the machinery and generality of the dielectric cylinder equations were used to create a broad scope over which the semi-empirical model would be valid. This validity region extends for grain head lengths between 5 and 9 cm, grain head diameters in the nominal range of 1 cm, and volumetric moisture contents ranging from very dry at 8% up to very lush at 19%.

191 vv-forward scatter coeffs freq ao a2 elo = 6.4 e12 = 2.1 L e20 = -0.34 e22 =-0.03 e30 = 6.8(-3) e32 = 0.0 e1o = 6.8 e2 = 2.1 C 20 = -0.30 622 = 0.07 e30 = 5.2(-3) 32 =-2.8(-3) elo = 7.4 e12 = 2.4 X e20 =-0.42 e22 =-0.13 e30 = 8.4(-3) 632 = 2.2(-3) Table 7.3: Forward-scatter coefficients for vv-polarization. 7.6 Semi-Empirical Model Coefficients Coefficients for the model given in Tables 7.3 - 7.8 were obtained by fitting the semiempirical expression to the mean scattered fields evaluated by the T-matrix model for many realizations of the grain head. As a method of saving space, exponential components are given in parenthesis (i.e. 10' 4 (x)).

192 hh-forward scatter coeffs freq a0 a2 elo = 16.5 e12 =-3.6(-3) L e20 = 5.0 e22 = -0.011 e30 = 0.0 32= 0.0 elo = 15.7 e12 = 0.84 C e20 = 7.4 e22 =-2.5 e30 = 0.0 e32 - 0.0 elo = 14.6 el2 =3.5 X e20 = 10.6 e22 = -9.3 e30 = 0.0 e32 = 0.0 Table 7.4: Forward-scatter coefficients for hh-polarization. vv-backscatter coeffs freq ao a; elo = 5.0 el2 = 1.8 ei = 0.073 L 20 = -0.25 e22 = -0.03 e2i = -6.2(-3) e30 5.2(-3) e32 = 0.0 e3 = 1.5(-4) elo = 6.5 e12 = 2.6 eli = 4.3(-4) C e20 =-0.43 e22 = -0.055 e2 = 1.5(-3) e30 = 9.2(-3) e32 = 0.0 e3i = 2.9(-5) elo = 4.3 e12 = 4.3 eli = 4.6(-3) X 20 = -0.39 e22 = -0.40 e2i = 3.0(-3) e30 = 9.4(-3) e32 = 9.6(-3) e3i = 7.0(-5) Table 7.5: Back-scatter coefficients for vv-polarization.

193 hh-backscatter coeffs freq ao a2 ai 61o = 12.9 e12 = 0.45 eli = 0.74 L 620 = 3.4 e22 = 0.49 c2i =-0.91 e30 = 0.0 e32 = 0.0 e3i = 0.37 elo = 15.9 e12 = 1.7 eli =-1.7(-3) C e20 = 6.5 622 =-5.8 e2i = 5.9(-3) ~30 = 0.0 e32 0 =C3i 0 e10 = 17.3 e12 = 1.74 ei = 0.076 X 620 = 4.0 e22 =-7.3 e2i -0.2 e30 - 0.0 e32 = 0.0 3i = 0.14 Table 7.6: Back-scatter coefficients for hh-polarization. vv-specular scatter coeffs freq ao a2 e1o = 28.8 e12 = 19.3 L 620 = -2.1 e22 = -1.4 630 = 0.046 632 = 0.031 e1o = 31.3 el2 = 22.2 C e20 =-1.9 e22 =-1.3 e3o = 0.037 e32 = 0.024 elo = 39.4 e12 = 29.7 X C20 =-3.0 e22 -2.7 e3o = 0.068 e32 = 0.063 Table 7.7: Specular-scatter coefficients for vv-polarization.

194 hh-specular scatter coeffs freq ao a2 elo = 16.5 e2 = -0.01 L 620 = 5.0 e22 - -0.018 e30 = 0.0 e32 = 0.0 eo1 = 15.4 e12 = 0.52 C e20 = 7.2 e22 =-2.7 e30 = 0.0 632 = 0.0 elo = 13.7 e12 = 2.6 X e20 = 10.3 e22 =-9.7 e30 = 0.0 e32 = 0.0 Table 7.8: Specular-scatter coefficients for hh-polarization.

CHAPTER 8 Conclusions and Directions for Future Work The contents of this dissertation summarizes a large volume of work whose aim it was to better understand how electromagnetic waves propagate through and interact with random media at frequencies where the wavelength exists on a similar dimensional scale as the discontinuities that make up the media. To investigate this phenomenon, Chapter 2 described an experimental proceedure to determine effective permittivity of a layer of sand at 35 GHz. The execution of this experiment led to several important theoretical questions, that of how to relate the permittivity of the sand to the physical and electromagnetic characteristics of the sand constituents (i.e. grains of silica). A variety of theoretical methods (the Polder-VanSanten mixing formula, the effective field approximation, the quasi-crystalline approximation, the Born approximation and Radiative transfer) for characterizing electromagnetic wave behavior for this purpose were described and developed in Chapter 3. These two chapters laid the groundwork for the treatment that followed. In Chapter 4, a packing algorithm was developed for simulating particle arrangements in dense, gravity deposited media (such as snow or sand). Results from this algorithm were shown to provide useful unknown functions for characterizing random media such as the correlation function of permittivity fluctuations or the pair distribution function. The packing algorithm may also be used as a first step in a full Monte-Carlo electromagnetic analysis for fields within random media. In two dimensions, this analysis can be conveniently be performed with the aid of the well developed method of moments technique for numerically solving Maxwell's equations. In three dimensions, the treatment requires a more elaborate technique that is amenable to Monte-Carlo analysis of a rondom collection of spherically shaped discontinuities. To this end, Chapter 5 was dedicated to the development of the T-matrix algorithm which relates the spherical wave expansion for an incident field upon a collection of scatterers to the spherical wave expansion of the scattered field. In this chapter, the recursive aggregate T-matrix algorithm (RATMA) was described where it was 195

196 shown in theory and by example the limited ability of RATMA to account for interactions between closely spaced particles. In Chapter 6, the work of the previous chapters is used to numerically address the calculation of effective permittivity in two and three dimensions. This was accomplished by comparing the average scattered fields from bounded sample of a random medium to that of a homogeneous dielectric. In two dimensions, this technique was used to evaluate the performance of the Polder-VanSanten mixing formula, the effective field approximaiton and the quasi-crystalline approximation and to demarcate regions of validity for these techniques. In three dimensions, the direct T-matrix algorithm was used to perform a similar task where estimates of effective permittivity were compared with results derived from the quasi-crystalline approximation with coherent potential, the effective field approximation and an incoherent numerical method for determing the extinction coefficient of a random medium. In the final chapter, Chapter 7, the tools and methodologies that had been developed in the previous six chapters were put to use to address the problem of scattering from the grain head of a wheat plant. In this chapter, the T-matrix algorithm is used to analyze an a grain model of the wheat head. Results from this analysis were used to generate a semi-empirical formulation based upon a low frequency analysis of the scattering from a dielectric cylinder. The purpose of this treatment was to utilize a high fidelity numerical method to calculate the principle scattering components and then in turn develop a simplified polynomial based scattering model that can be used as a tool in computational analysis of the wheat field as a whole. 8.1 Summary of Accomplishments The accomplishments that have been made in this thesis can be summarized as follows 1. Characterization and bistatic measurement of find sand at 35 GHz. Measurements were used to detrimine the real component of effective permittivity for the medium. 2. Development of the two-dimensional version of the quasi-crystalline approximation. 3. Development of the two-dimensional version of the Born approximation and graphic relation of the predicted backscatter to the power spectral density of the permittivity fluctuation correlation function.

197 4. The development of the packing algorithm for efficiently determining particle arrangements in two and three dimensions. 5. Implementation and evaluation of the recursive aggregate T-matrix algorithm [Wang and Chew 1993]. 6. Numerical determination of effective permittivity in both two and three dimensions. 7. Evaluation of the quasi-crystalline approximation, the effective field approximation and the Polder-VanSanten mixing formula in both two and three dimensions. 8. Development of a semi-empirical model for the scattering and extinction due to wheat grain heads. 8.2 Directions for Future Work The breadth of analysis that has been performed for this dissertation leads to a variety of different research directions that could be followed in the future. The directions envisioned by the author are * Analyze bistatic and backscatter controlled measurements performed on materials such as sand to provide the imaginary component of effective permittivity of this medium. The theoretical/numerical tools developed in this dissertation could be used to compare with results obtained from this experiment. * Quantify the effects of material anisotropy for three-dimensional theoretical methods such as the Born approximation and the quasi-crystalline approximation. To this end, alternative functional forms of the correlation function and pair distribution function could be developed. * Expand the three-dimensional version of the packing algorithm to account for other than spheroidal particles. The packing algorithm may be of great use to such disciplines as material science, civil and soil engineering. These avenues need to be further explored. * Implenient alternative numerical techniques (including different formulations of the T-matrix algorithm) to improve or surpass the difficulties that RATMA has with accounting for particle interactions.

198 * Use the Monte-Carlo results obtained in determining effective permittivity as a basis for developing a theoretical method to circumvent much of the numerical burden that the method requires for addressing the three-dimensional problem. * Use the methodology for deriving a semi-empirical scattering model of the wheat grain head to address other remote sensing problems of interest.

BIBLIOGRAPHY [1] Appleton, P.N., P.R. Siqueira, and J.P. Basart, "A Morphological Filter for Removing 'Cirrus-Like' Emission from Far-Infrared Extragalactic IRAS Fields", Astron. J., 106(4), 1664-1678, 1993. [2] Broyles, A.A., S.U. Chung, and H.L. Sahlin, "Comparison of Radial Distribution Functions from Integral Equations and Monte-Carlo," J. Chem. Phys., 37(10), 2462-2469, 1962. [3] Buchalter, B.J. and R.M. Bradley, "Orientational Order in Random Packings of Ellipses", Phys. Rev. A, 46(6), 3046-3056, 1992. [4] Buchalter, B.J. and R.M. Bradley, "Orientational Order in Amorphous Packings of Ellipsoids," Europhys. Lett., 26(3), 159-164, 1994. [5] Chandrasekhar, S., Radiative Transfer, 393 pp., Dover Publications, New York, 1960. [6] Chew, W.C. and C.C. Lu, "The Recursive Aggregate Interaction Matrix Algorithm -for Multiple Scatterers," IEEE Trans. Ant. Prop., AP-43(12), 1483-1486, 1995. [7] Chew, W.C., C.C. Lu, and Y.M. Wang, "Efficient Computation of Three-Dimensional Scattering of Vector Electromagnetic Waves," J. Opt. Soc. Amer. A, 11(4), 1528-1537, 1994. [8] Chew, W.C. and Y.M. Wang, "Efficient Ways to Compute the Vector Addition Theorem," J. Electromagn. Wav. Appl., 7(5), 651-665, 1993. [9] Chew, W.C., Y.M. Wang, and L. Gurel, "Recursive Algorithm for Wave-Scattering Solutions Using Windowed Addition Theorem," J. Electromagn. Waves Appl., 6(11), 1537-1560, 1992. [10] Chew, W.C., "Recurrence Relations for Three-Dimensional Scalar Addition Theorenm," J. Electromagn. Waves Appl., 6(2), 133-142, 1992. [11] Chew, W.C., Waves and Fields in Inhomogeneous Media, 608 pp., IEEE Press, New York, 1990. [12] Chew, W.C. and Y.M. Wang, "A Fast Algorithm for Solution of a Scattering Problem Using a Recursive Aggregate r Matrix Method," Micro. Opt. Tech. Lett, 3(5), 164-169, 1990. [13] Debye, P., H.R. Anderson, and H. Brumberger, "Scattering by an Inhomogeneo-is Solid. II. The Correlation Function and Its Applications," J. Appl. Phys., 28(6), 679 -683, 1957. 199

200 [14] Debye, P. and Bueche, A.M. "Scattering by an Inhomogeneous Solid," J. Appl. Phys., 20, 518-525, 1949. [15] DeLoor, G., "Method of obtaining information on the internal dielectric constant of mixtures," Appl. Sci. Res., 3(B) 479-482, 1953. [16] Ding, K.H., C.E. Mandt, L. Tsang, and J.A. Kong, "Monte-Carlo Simulations of Pair Distribution Functions of Dense Discrete Random Media With Multiple Sizes of Particles", J. Electromagn. Waves Appl., 6(8), 1015-1030, 1992. [17] Ding, K.H. and L. Tsang, "Effective Propagation Constants and Attenuation Rates in Media of Densely Distributed Coated Dielectric Particles with Size Distributions," J. Electromagn. Waves Appl., 5(2), 117-142, 1991. [18] El-Rayes, M.A. and F.T. Ulaby, "Microwave dielectric spectrum of vegetation, Part I: Experimental Observations," IEEE Trans Geosci. Rem. Sens., GE-25(5), 541-549, 1987. [19] Fisher, I.Z., Statistical Theory of Liquids, 335 pp., University of Chicago Press, Chicago, IL 1961. [20] Fung, A.K., "A review of Volume Scatter Theories for Modeling Applications," Radio Sci., 17(5), 1007-1017, 1982. [21] Giardina, C.R. and E.R. Dougherty, Morphological Methods in Signal and Image Processing, 321 pp., Prentice Hall, Englewood Cliffs, NJ, 1988. [22] Gurel, L. and W.C. Chew, "Recursive T-Matrix Algorithm for the Solution of Electromagnetic Scattering from Strip and Patch Geometries," IEEE Trans. Ant. Prop., AP-41(1), 91-99, 1993. [23] Holloway, D., The Physical Properties of Glass, Springer-Verlag, New York, 1973. [24] Lax, M., "Multiple Scattering of Waves. II. The Effective Field in Dense Systems," Phys. Rev., 85(4), 621-629, 1952. [25] Mandt, C.E., Y. Kuga, L. Tsang and A. Ishimaru, "Microwave Propagation and Scattering in a Dense Distribution of Non-tenuous Spheres: Experiment and Theory," Waves in Random Media, 2, 225-234, 1992. [26] McQuarrie, D.A. Statistical Mechanics. Harper & Row, Evanston, IL. 1976. [27] Nashashibi, A.Y., Microwave and Millimeter-wave Propagation and Scattering in Dense Random Media: Modeling and Experiments, 180 pp., Univ. of Michigan Dissertation Thesis, 1995. [28] Percus, J. and G. Yevick, "Analysis of Classical Statistical Mechanics by Means of Collective Coordinates", Phys. Rev., 110, 1-13, 1958. [29] Peterson, B. and S. Strom, "T-Matrix for Electromagnetic Scattering from an Arbitrary Number of Scatterers and Representations of E(3)," Phys. Rev. D, 8(10), 3661-3678, 1973.

201 [30] Polder, D. and J. van Santen,"The Effective Permeability of Mixtures of Solids," Physica, 12(5), 1257-271, 1946. [31] Richmond, J.H., "Scattering by a dielectric cylinder of arbitrary cross-section shape," IEEE Trans. Ant. Prop., AP-13, 334-341, 1965. [32] Richmond, J.H., "TE-Wave Scattering by a Dielectric Cylinder of Arbitrary CrossSection Shape," IEEE Trans. Ant. Prop., AP-14, 460-464, 1966. [33] Ruck, G.T., D.E. Barrick, W.D. Stuart and C.K. Krichbaum, Radar Cross Section Handbook: Volume 1, 472 pp., Plenum Press, New York, 1970. [34] Sarabandi, K. and T.B.A. Senior, "Low-Frequency Scattering from Cylindrical Structures at Oblique Incidence," IEEE Trans. Geosci. Rem. Sens., 28(5), 879-885, 1990. [35] Sarabandi, K., "A Technique for Dielectric Measurement of Cylindrical Objects in a Rectangular Waveguide", IEEE Trans. Instr. Meas., 43(6), 793-798, 1994. [36] Sarabandi, K., F. Ulaby and M. Tassoudji. "Calibration of Polarimetric Radar Systems with Good Polarization Isolation." IEEE Trans. Geosci. Rem. Sens., 28(1), 70-75,1990. [37] Sarabandi, K. Electromagnetic Scattering from Vegetation Canopies, 345 pp., Univ. of Michigan Dissertation Thesis, 1989. [38] Serra, J.P. Image Analysis and Mathematical Morphology, 601 pp., Academic Press, New York, 1982. [39] Siqueira, P.R. and K. Sarabandi, "Method of Moments Evaluation of the TwoDimensional Quasi-Crystalline Aprroximation," IEEE Trans. Ant. Prop., AP-44(8), 1996. [40] Siqueira, P.R., K. Sarabandi and F.T. Ulaby, "Numerical Simulation of Scatterer Positions in a Very Dense Media," IEEE Ant. Prop. Conf. Digest, 1994. [41] Stark, H., and J.W. Woods, Probability, Random Processes, and Estimation Theory for Engineers, 618 pp. Prentice-Hall, Englewood Cliffs, New Jersey, 1986. [42] Stiles, J.M. A Coherent, Polarimetric Microwave Scattering Modelfor Grassland Structures and Canopies, 307 pp., Univ. of Michigan Dissertation Thesis, 1996. [43] Tai, C.T. Dyadic Green Functions in Electromagnetic Theory (second edition), 343 pp., IEEE Press, New York, 1994. Tinga W., W. Voss and D. Blossey. "Generalized Approach to Multiphase Dielectric Mixture Theory." J. Appl. Phys., 44(9), 3897-3902, 1973. [44] Tsang, L., C.E. Mandt, and K.H. Ding, "Monte Carlo simulations of the extinction rate of dense media with randomly distributed dielectric spheres based on solution of Maxwell's equations," Opt. Lett., 17(5), 314-316, 1992. [45] Tsang, L. "Dense Media Radiative Transfer Theory for Dense Discrete Random Media with Particles of Multiple Sizes and Permittivities," Prog. Electromagn. Res., Elsevier, New York, 1992.

202 [46] Tsang, L. and J.A. Kong, "Scattering of Electromagnetic Waves from a Dense Medium Consisting of Correlated Mie Scatterers with Size Distributions and Applications to Dry Snow," J. Electromagn. Waves Appl., 6(3), 265-286, 1992. [47] Tsang, L. and A. Ishimaru, "Radiative Wave Equations for Vector Electromagnetic Propagation in Dense Non-tenuous Media," J. Electromagn. Waves Appl., 1(1), 59-72, 1987. [48] Tsang, L., J. Kong, and R. Shin, Theory of Microwave Remote Sensing, 613 pp., John Wiley & Sons, New York, 1985. [49] Ulaby, F.T. and C. Elachi, Radar Polarimetry for Geoscience Applications, 364 pp., Artech House, Norwood, MA, 1990. [50] Ulaby, F.T., R.K. Moore and A.K. Fung. Microwave Remote Sensing: vols. I, II, and III, 2162 pp., Artech House, Norwood, MA, 1981. [51] Ulaby, F.T., T.F. Haddock and M.E.Coluzzi. "Millimeter-Wave Bistatic Radar Measurements of Sand and Gravel," IGARSS '87 Symp. Proc., Ann Arbor, MI, 1987. [52] Ulaby, F.T., M. Whitt and K. Sarabandi, "AVNA Based Scatterometers," IEEE Ant. Prop. Mag., 32(10), 6-17, 1990. [53] Vallese, F. and J.A. Kong, "Correlation Function Studies for Snow and Ice," J. Appl. Phys., 52(8), 4921-4925, 1981. [54] van Beek, L., Dielectric Behaviour of Heterogeneous Systems, Ed: J. Birks, London Heywood Books, London, 1967. [55] Visscher, W.M. and M. Bolsterli "Random Packing of Equal and Unequal Spheres in Two and Three Dimensions," Nature, 239, 504-507, 1972. [56] von Hippel, A., Dielectric Materials and Applications, John Wiley and Sons, New York, 1954. [57] Wang, Y.M. and W.C. Chew, "A Recursive T-Matrix Approach for the Solution of Electromagnetic Scattering by Many Spheres," IEEE Trans. Ant. Prop., AP-41(12), 1633-1639, 1993. [58] Wang, Y.M. and W.C. Chew, "An Efficient Algorithm for Soluiton of a Scattering Problem," Microwave Opt. Tech. Lett., 3(3), 102-106, 1990. [59] Waterman, P.C., "Matrix Formulation of Electromagnetic Scattering," IEEE Proc., 53, 805-811, 1956. [60] Wen, B., L. Tsang, D. Winebrenner, and A. Ishimaru, "Dense Media Radiative Transfer Theory: Comparison with Experiment and Application to Microwave Remote Sensing Polarimetry," IEEE Trans. Geosci. Rem. Sens., GE-28(1), 46-59, 1990.