1.1 A Brief History of the Geometric Optics Method for Light-Scattering Calculations

Knowledge of light scattering by particles is inevitably necessary for numerous disciplines, including astrophysics, planetary sciences, bio-optics, radiative transfer, remote sensing, climate science, and a wide range of subjects covered by applied optics. In the case of spherical particles, the rigorous Lorenz–Mie theory has been well established (Lorenz, Reference Lorenz1890; Mie, Reference Mie1908). Of particular note is that the theory for light scattering by a sphere is usually referred to as the Mie theory, without giving appropriate credit to Ludvig Valentin Lorenz (1829–1891) for his pioneering work in this regard. For fairness from the historical perspective, it has been advocated that the theory of light scattering by a sphere be called the Lorenz–Mie theory (Kragh, Reference Kragh1991). Early studies of the scattering of plane waves by spheres were reviewed by Logan (Reference Logan1965). More recently, a volume edited by Hergert and Wriedt (Reference Hergert and Wriedt2012) is dedicated to a comprehensive review of the history of the Lorenz–Mie theory and its extensions. Although the Lorenz–Mie theory is regarded as a solved problem in classical physics, some fascinating new features associated with light scattering by a sphere were discovered in the twentieth and twenty-first centuries, such as the narrow resonance structure in the Lorenz–Mie scattering characteristics (Chýlek et al., Reference Chýlek, Kiehl and Ko1978) and the zero-backscattering of a sphere when the particle’s relative impedance with respect to the counterpart of the ambient medium is unity (Zhai et al., Reference Zhai2008).

Historically, light-scattering computations based on the Lorenz–Mie theory were challenging because of numerical instability, particularly in the case of an absorptive sphere. Fortunately, numerical stability has been achieved through the utilization of a downward recurrency scheme in light-scattering computation based on the Lorenz–Mie theory, as reported by Kattawar and Plass (Reference Kattawar and Plass1967). Furthermore, we call special attention to a milestone contribution made by Wiscombe (Reference Wiscombe1980) for the computational convergence and efficiency in the numerical implementation of the Lorenz–Mie theory. With various robust open-access computational programs based on the Lorenz–Mie theory, it is now regarded as a trivial effort to compute the optical properties of a sphere of essentially any practical size, ranging from the Rayleigh regime (the particle size is much smaller than the wavelength of the incident light) to the geometric optics regime (the particle size is much larger than the wavelength of the incident light). In particular, special software was already developed for the Lorenz–Mie computation using a mobile phone even in the 2000s (Wriedt, Reference Wriedt2008).

Despite the success in dealing with the scattering of light by a sphere based on the Lorenz–Mie theory, many particles in nature or those generated in laboratories, such as ice crystals in cirrus clouds, airborne dust aerosol particles in the atmosphere, interstellar dust particles, many oceanic particles, and artificial dust aerosol particles, are nonspherical. For example, Figure 1.1 shows the scanning electron microscope (SEM) images of feldspar particles that have various irregular shapes. These particles were generated by crushing large rocks (Volten et al., Reference Volten2005). It is highly unlikely to account for the morphological details of these particles in light-scattering computation. As a first-order approximation, the spheroid model has often been used to account for the nonsphericity effects in downstream applications in remote sensing of dust aerosol properties. The spheroid model represents a significant improvement compared to the sphere model because the spheroidal geometry has two degrees of freedom (i.e., particle’s maximum dimension and the aspect ratio, which is the ratio of particle’s major axis to its minor axis), whereas spherical geometry has only one degree of freedom (i.e., particle size).

Figure 1.1

SEM image of feldspar particles. The particles are highly irregular in shape and show considerable differences in size. The scale of the illustration is indicated by the vertical bar on the right.

Adapted from Volten et al. (2005). Reprinted from Journal of Quantitative Spectroscopy and Radiative Transfer, Vol. 90, H. Volten, O. Muñoz, J. W. Hovenier, J. F. de Haan, W. Vassen, W. J. van der Zande, and L. B. F. M. Waters, WWW scattering matrix database for small mineral particles at 441.6 and 632.8 nm, Pages 191–206, Copyright (2005), with permission from Elsevier.

Figure 1.2 demonstrates the impact of a particle’s nonsphericity on the atmospheric radiation field and shows two solar reflectance spectra at the top of the atmosphere (TOA), which are simulated with the US standard atmosphere for an overcast, dusty sky. The upper and lower spectral curves in Figure 1.2 correspond to using the spheroid and sphere models, respectively. For the simulations, the solar zenith angle is 60° and the viewing angle is 0°, the aerosol optical depth (AOD) of the dust layer is 1.0 at a wavelength of 0.55 µm, and the surface albedo is 0.1. It can be seen from Figure 1.2 that neglecting the nonsphericity of dust particles can significantly underestimate the TOA reflectance. This underestimation has an important implication for assessing the radiative forcing effect of airborne dust, which is an important factor involved in climate studies.

Figure 1.2

Simulated solar reflectance spectra from 0.2 to 2.0 microns at the top of the atmosphere under overcast, dusty sky conditions with spherical and nonspherical (specifically, spheroidal) particles.

Data adapted from Yang et al. (2007). Reprinted from Journal of Aerosol Science, Vol. 38, Ping Yang, Qian Feng, Gang Hong, George W. Kattawar, Warren J. Wiscombe, Michael I. Mishchenko, Oleg Dubovik, Istvan Laszlo, and Irina N. Sokolik, Modeling of the scattering and radiative properties of nonspherical dust-like aerosols, Pages 995–1014, Copyright (2007), with permission from Elsevier.

To demonstrate the effect of a particle’s nonsphericity on remote sensing of dust aerosol properties from spaceborne observations, we select a scene observed by the Moderate Resolution Imaging Spectroradiometer (MODIS) on the Terra satellite (King et al., Reference King1992) over West Africa on March 2, 2003. In particular, Figure 1.3 (Feng et al., Reference Feng2009) shows the RGB (R: band 1 covering from 620 to 670 nm; G: band 4 covering from 545 to 565 nm; and B: band 3 covering from 459 to 479 nm) false-color image based on MODIS observations. Dust transport from land to the ocean can be clearly seen in the image. The area indicated by the rectangular box is selected for the inference of dust AOD.

Figure 1.3

MODIS RGB image on March 2, 2003, showing a dust plume over West Africa. The area indicated by the small box is used to retrieve dust AOD.

Adapted from Feng et al. (2009). Reprinted from Journal of Aerosol Science, Vol. 40, Qian Feng, Ping Yang, George W. Kattawar, Christina N. Hsu, Si-Chee Tsay, and Istvan Laszlo, Effects of particle nonsphericity and radiation polarization on retrieving dust properties from MODIS observations, Pages 776–789, Copyright (2009), with permission from Elsevier.

Figure 1.4a,b shows the dust AOD retrievals for the area indicated by the small box in Figure 1.3, for which nonspherical (spheroidal) and spherical dust models are assumed, respectively. From the comparisons shown in Figure 1.4c,d, it can be seen that the impact of dust particle nonsphericity on the inferred dust AOD is significant, particularly for those pixels corresponding to scattering angles larger than ~140°. Moreover, in an earlier study, Dubovik et al. (Reference Dubovik2006) showed that using the sphere model may result in an erroneous fine mode in the retrieved aerosol size distribution and a spurious spectral variation of the real part of the refractive index in the visible wavelength regime, as demonstrated by remote sensing of aerosol properties from the ground-based radiometric measurements provided by the AERosol Robotic NETwork (AERONET) (Holben et al., Reference Holben1998). These studies (Dubovik et al., Reference Dubovik2006; Feng et al., Reference Feng2009) show that it is necessary to consider particle nonsphericity for spaceborne and ground-based remote-sensing implementation.

Figure 1.4

(a)–(b) The retrieved dust AOD values based on the spheroid (a) and sphere (b) models. (c) Retrieved dust AOD based on the sphere model versus the counterpart based on the spheroid model. (d) The relative differences between the AOD retrievals based on the spheroid and sphere models, with the former as the benchmark.

Adapted from Feng et al. (2009). Reprinted from Journal of Aerosol Science, Vol. 40, Qian Feng, Ping Yang, George W. Kattawar, Christina N. Hsu, Si-Chee Tsay, and Istvan Laszlo, Effects of particle nonsphericity and radiation polarization on retrieving dust properties from MODIS observations, Pages 776–789, Copyright (2009), with permission from Elsevier.

Although the nonsphericity effect of scattering particles is critical to many applications, computing the optical properties of a nonspherical particle poses a significant challenge, particularly from the perspective of numerical computation (Mishchenko et al., Reference Mishchenko, Hovenier and Travis2000). Exact methods include the method of separation of variables, which in practice is applied only to regular geometries, such as spheroidal particles (Oguchi, Reference Oguchi1973; Asano and Yamamoto, Reference Asano and Yamamoto1975; Ding and Yang, Reference Ding and Yang2023), the T-matrix method implemented via the extended boundary condition approach (Waterman, Reference Waterman1965, Reference Waterman1971; Barber and Yeh, Reference Barber and Yeh1975; Tsang et al., Reference Tsang, Kong and Shin1985; Barber and Hill, Reference Barber and Hill1990; Mishchenko and Travis, Reference Mishchenko and Travis1994; Mishchenko et al., Reference Mishchenko, Travis and Mackowski1996; Doicu et al. Reference Doicu, Wriedt and Eremin2006), and the invariant imbedding T-matrix method (IITM) (Johnson, Reference Johnson1988; Bi et al., Reference Bi2013; Bi and Yang, Reference Bi and Yang2014a). Numerically feasible and robust methods, such as the finite-difference time domain (FDTD) method (Yee, Reference Yee1966; Yang and Liou, Reference Yang and Liou1996a; Sun et al., Reference Sun, Fu and Chen1999; Taflove and Hagness, Reference Taflove and Hagness2005) and the discrete dipole approximation (DDA) (Purcell and Pennypacker, Reference Purcell and Pennypacker1973; Draine and Flatau, Reference Draine and Flatau1994; Yurkin and Hoekstra, Reference Yurkin and Hoekstra2007), are usually applied to small-to-moderate-size particles in comparison with the wavelength of the incident electromagnetic wave. Moreover, these methods require considerable computational effort for moderate-size particles. Many nonspherical particles in nature, such as ice crystals and coarse-mode dust particles in the atmosphere and phytoplankton in the ocean, are much larger than visible wavelengths. In addition, it is necessary to consider particle orientations (usually random orientations) for many practical applications, and averaging particle optical properties over orientations substantially increases computational burden. To compute the optical properties of these particles at the aforementioned wavelengths, the only viable option is to employ the principles of geometric optics as an approximate method, with accuracy increasing as the particle size increases; in other words, the geometric optics method represents a case where there is “no alternative” for computing the optical properties of nonspherical particles that are much larger than the wavelength of the incident electromagnetic wave.

Classical optics has a long history. Born and Wolf (Reference Born and Wolf1999) and Strong (Reference Strong1958) briefly reviewed the milestones in the development of classical optics. Of the principles of classical optics, which were essentially fully developed before the twentieth century, in this book, we are mainly concerned with the concepts of a geometric optics “ray,” Snell’s law, the reflection and refraction of a ray, and Fraunhofer diffraction. The reason is that scattered light within the framework of classical optics is assumed to be a combination of classical diffraction and geometrical ray reflection and transmission (Hodkinson and Greenleaves, Reference Hodkinson and Greenleaves1963). It is not feasible for us to trace the origin of the concept of a ray. As articulated by Born and Wolf (Reference Born and Wolf1999), the concept of a ray can be well illustrated by thinking of the passage of collimated light through a tiny opening in an opaque screen, which, at first sight, illuminates a region with a sharp boundary. As we will explain in detail in Chapter 2, the validity of the ray concept actually depends on the cross section of the ray beam. Studies of basic ray or beam propagation, as explained by Strong (Reference Strong1958), can be traced to the effort of Euclid (~300 BC) describing ray reflection and of Ptolemy (AD 70-147) describing ray refraction. The law of refraction, which is known as Snell’s law, was attributed to the experimental work of Willebrord Snellius (c. 1580–1626) or Snell. It has been advocated that Snell’s law be called Snel’s law because of the original spelling of his name, Willebrord Snel van Royen (Adams and Kattawar, Reference Adams and Kattawar1997).

Light (and electromagnetic radiation in general) has a very complex nature, and scientific knowledge of the many phenomena associated with light has only gradually evolved over many centuries. Historically, light propagation was regarded as being in the form of minute particles, for example, by Newton. Thomas Young’s double-slit interference experiment conducted at the beginning of the nineteenth century played an instrumental role in the general acceptance of the wave theory of light, for which Robert Hooke (1635–1703) and Christiaan Huygens (1629–1695) were pioneers. In 1816, Augustin-Jean Fresnel (1788–1827) developed his theory of diffraction without considering polarization. In 1823, Fresnel derived formulas, which are now called Fresnel equations, for the amplitude, phase, and polarization of the transmitted and reflected waves. For a more in-depth study of the development of polarization of light, the reader is referred to Brosseau (Reference Brosseau1998).

The earliest quantitative applications of the principles of geometric optics to light scattering were developed for spherical particles. For example, Mecke (Reference Mecke1920) studied the scattering of light incident on a sphere by using the principles of geometric optics and the diffraction associated with a circular disc of the same radius as the sphere. He fully considered the phase of all rays. Bricard (Reference Bricard1946) suggested that, when $2x(m-1)\gg 1$, where m is the index of refraction and x is the size parameter defined as 2$\pi$ times the ratio of particle radius to the wavelength of the incident light, the interference between contributions of the geometric optics rays and diffraction is essentially incoherent, and the two contributing components can be treated separately. It should be pointed out that, in classical geometric optics methods, the geometric optics rays and diffraction equally contribute to the total extinction. The extinction due to diffraction stems from the corrective field in the shadow region of a sphere illuminated by a collimated beam, as illustrated through a rigorous analysis by Brillouin (Reference Brillouin1949), who clarified some inconsistencies and misinterpretations found in several early texts regarding the fact that the extinction cross section is twice the geometric projected area in the case of light scattering by a large particle compared to the wavelength of incident light.

In the history of studying light scattering by large particles, there was intense debate as to whether the extinction efficiency (defined as the ratio of the extinction cross section to the geometric projected area of the particle) is always 2. To clarify the apparently dubious physical process, let us consider the scattering of light by a perfectly conducting sphere. In this case, since there is no absorption involved, the extinction efficiency equals the scattering efficiency (defined as the ratio of the scattering cross section to the geometric projected area of the particle). As shown in Figure 1.5, blocking of the incident light leads to a conical shadow region marked as ABC in the schematic diagram in Figure 1.5. Outside the shadow region, there is a penumbra region. Beyond the penumbra region, the scattered field is essentially given by the external reflection of the incident light by the sphere. With an increase in the ratio of particle size to the wavelength of the incident light, angle $\alpha$ in Figure 1.5 decreases, the shadow extends, and the penumbra shrinks. In the geometric optics limit (i.e., in the case of light scattering by a very large particle), it is assumed that angle $\alpha$ is essentially zero and the blocking of the incident light leads to a cylindrical shadow region behind the sphere. The blocking of the incident light corresponds to electromagnetic radiation flux given by $\pi r^2{I}^i$, where $I^i$ is the incident irradiance. Thus, the contribution to the extinction efficiency due to the blocking effect is 1. Furthermore, near the surface of the sphere within the shadow region, the total electric and magnetic fields are zero, given by

${\overrightarrow{E}}^t={\overrightarrow{E}}^i+{\overrightarrow{E}}^c=0,$(1.1.1a)

${\overrightarrow{H}}^t={\overrightarrow{H}}^i+{\overrightarrow{H}}^c=0,$(1.1.1b)

where superscripts t, i, and c indicate the total, incident, and corrective fields, respectively. Equations (1.1.1a) and (1.1.1b) indicate that the corrective field cancels the incident field. It is the corrective field that leads to an additional attenuation of the incident light with an amount of $\pi r^2{I}^i$. Thus, the total extinction efficiency is 2. The reader is referred to Brillouin (Reference Brillouin1949) for a detailed mathematical derivation of the contribution of the corrective fields to the extinction efficiency. The contribution of the corrective fields to the scattered fields is a physical phenomenon referred to as diffraction, which is fundamentally due to the incomplete wave front of a propagating wave. Section 1.5 provides a more detailed explanation of diffraction, particularly far-field diffraction, known as Fraunhofer diffraction (Born and Wolf, Reference Born and Wolf1999).

Figure 1.5

Schematic diagram illustrating the scattering of light by a perfectly reflective opaque sphere (i.e., a perfectly reflective sphere) within the geometric optics framework. Diagram concept based on Brillouin (Reference Brillouin1949).

Following the studies of Mecke (Reference Mecke1920), Bricard (Reference Bricard1946), and van de Hulst (Reference van de Hulst1957), Ellison and Peetz (Reference Ellison and Peetz1959) reported the forward scattering of light by spheres according to geometric optics. Hodkinson and Greenleaves (Reference Hodkinson and Greenleaves1963) also performed numerical computation of light scattering by a sphere using the geometric optics method for scattering angles up to 40º. They also compared the geometric optics results with the Lorenz–Mie solutions for size parameters up to 30. A comprehensive formalism of the light scattering by a sphere within the framework of geometric optics can be found in van de Hulst (Reference van de Hulst1957). Liou and Hansen (Reference Liou and Hansen1971) applied the geometric optics method to light scattering by polydisperse spheres (an ensemble of spheres with a distinct distribution of sizes) to compute the phase function and the degree of linear polarization. They compared the geometric optics method with the exact Lorenz–Mie theory and found that the two methods were in close agreement for real refractive indices m = 1.33 and 1.50 for an ensemble of polydisperse, homogeneous, and spherical particles when the effective particle size (Hansen and Travis, Reference Hansen and Travis1974) parameter is larger than approximately 100. At a size parameter of 400, the phase function results based on the geometric optics method converge to the Lorenz–Mie results, although the backscattering peak corresponding to an optical phenomenon known as the glory does not appear in the geometric optics solution in the case of m = 1.33, which is approximately the index of refraction of water at a visible wavelength, regardless of the size parameter.

Application of the geometric optics method has been extended to nonspherical particles, such as the scattering of light by spheroids (Lock, Reference Lock1996a,Reference Lockb). In particular, Lock (Reference Lock1996a) obtained a closed-form solution for diffraction and external reflection for an arbitrarily oriented spheroid. In the case of light scattering by a transparent spheroid, the curvature effect of the spheroidal surface poses a significant challenge for performing rigorous ray-tracing calculations. In particular, the equal-phase wavefront associated with the refracted rays is distorted in reference to that of the incident collimated beam (a bundle of rays) after passing through a curved surface. As a result, the cross-section of a beam with a curved wavefront during propagation in a homogeneous medium varies with the ray-bundle propagation distance. Moreover, the phase of a beam with a convergent wavefront undergoes a phase shift of π/2 with each passage of a focal line. The contribution to the scattering amplitude by a beam after passing over a scatterer is determined by three elements: the cross-section variation of its wavefront, the total wave phase, and the refraction coefficients determined by the Fresnel equations. Ren et al. (Reference Ren2011) and Sun et al. (Reference Sun2016) developed a complex vector ray model to overcome the difficulties encountered in considering the wavefront and wave phase in the case of light scattering by a spheroid.

The geometric optics method has been especially applied to light scattering by faceted particles (ice crystals with flat, smooth, polygonal faces), mainly due to the importance of ice clouds in atmospheric remote sensing and climate science. Ice crystals in the atmosphere are almost exclusively nonspherical particles (e.g., Heymsfield and Iaquinta, Reference Heymsfield and Iaquinta2000; Bailey and Hallett, Reference Bailey and Hallett2009), ranging from pristine plates and columns to highly irregular aggregates with small-scale irregularities and surface roughness. The earliest study of the optical properties of hexagonal ice crystals was conducted by Jacobowitz (Reference Jacobowitz1971), as part of the author’s PhD dissertation at Massachusetts Institute of Technology (MIT). In the 1980s and 1990s, the geometric optics method, by means of the ray-tracing technique, was extensively used to investigate the single-scattering properties of nonspherical ice crystals (Wendling et al., Reference Wendling, Wendling and Weickmann1979; Cai and Liou, Reference Cai and Liou1982; Takano and Jayaweera, Reference Takano and Jayaweera1985; Takano and Liou, Reference Takano and Liou1989; Macke, Reference Macke1993; Macke et al., Reference Macke, Mueller and Raschke1996). Using the Monte Carlo method, Wendling et al. (Reference Wendling, Wendling and Weickmann1979) might be the first to quantify the 22° and 46° halo peaks and a scattering maximum at 154° in the phase function of randomly oriented ice crystals. However, in the Wendling et al. (Reference Wendling, Wendling and Weickmann1979) study, diffraction by a hexagonal ice crystal is approximated by that of a rectangular aperture regardless of the orientation of the ice crystal; moreover, polarization is not included; consequently, their results are limited only to the phase function. Cai and Liou (Reference Cai and Liou1982) fully incorporated the polarization configuration and phase interference in the ray-tracing calculation for randomly oriented hexagonal ice crystals. They also rigorously considered the diffraction pattern associated with the projected area of a three-dimensional (3D) hexagonal particle. However, the phase matrix elements P₃₃ and P₄₄ reported by Cai and Liou (Reference Cai and Liou1982) are incorrect, which were later corrected by Takano and Jayaweera (Reference Takano and Jayaweera1985). The theoretical framework for the “conventional” geometric ray-tracing technique and the associated computational algorithm developed by Cai and Liou (Reference Cai and Liou1982) for hexagonal columns and plates were further improved and refined by a number of researchers. In particular, the single-scattering properties associated with more complex particle morphologies (Macke, Reference Macke1993; Macke et al., Reference Macke, Mueller and Raschke1996; Takano and Liou, Reference Takano and Liou1995; Yang and Liou, Reference Yang and Liou1998) have been studied with a similar methodology implemented with the Monte Carlo technique. Moreover, Takano and Liou (Reference Takano and Liou1989) incorporated the effects of birefringence of ice into the ray-tracing computation.

An innovative method was developed by Muinonen (Reference Muinonen1989) to improve the conventional ray-tracing method by employing a modified Kirchhoff approximation, which maps the near field obtained with the ray-tracing technique to the far field. In this form, the spread of a ray due to the phase interference can be accounted for. As a result, unlike the results given by the conventional ray tracing techniques, the scattering peaks associated with halos in the case of the scattering of visible light by a hexagonal ice crystal under the random orientation condition are not independent of the size parameter for the solution given by the modified Kirchhoff approximation approach. However, the modified Kirchhoff approximation approach developed by Muinonen (Reference Muinonen1989) separates the contributions of diffraction and the geometric rays in the same way as the conventional geometric optics method, by assuming that the extinction efficiency is 2 regardless of the size parameter. To overcome this shortcoming, the near field computed with the ray-tracing procedure with a full inclusion of the phase and the field polarization configuration can be mapped to the far field through an exact electromagnetic surface integral equation (Yang and Liou, Reference Yang and Liou1996b) or an exact electromagnetic volume integral equation (Yang and Liou, Reference Yang and Liou1997). The surface- and volume-integral-based methods produce consistent results, particularly for relatively large particles compared to the wavelength of the incident light (Yang et al., Reference Yang2019). To the best of our knowledge, these improvements for the first time capture the variations of the extinction efficiency as a function of the size parameter within the framework of the geometric optics method (Yang and Liou, Reference Yang and Liou1995), and this improvement is referred to as the physical-geometric optics method (PGOM). The singularity known as the delta-transmission (Takano and Liou, Reference Takano and Liou1989), inherent to the conventional ray-tracing calculation in the case of a scattering particle with parallel surface facets, is overcome by the aforementioned PGOM method. Furthermore, PGOM is applicable to a much wider size parameter range compared to the conventional geometric optics method, as demonstrated by Yang et al. (Reference Yang2019). The volume-integral approach can be numerically implemented in terms of a ray-by-ray integration (RBRI) technique (Yang and Liou, Reference Yang and Liou1997) that, unfortunately, is computationally inefficient. Bi et al. (Reference Bi2011) and Sun et al. (Reference Sun, Yang, Kattawar and Zhang2017b) further improved the RBRI technique by implementing simultaneous tracing of a bundle of rays (referred to as a beam) impinging on the same facet. In addition, an improved method referred to as the beam-splitting method (Bi et al., Reference Bi2011) in conjunction with computer graphics techniques (Sun et al., Reference Sun, Yang, Kattawar and Zhang2017b) has been developed for a faceted convex particle that may contain a core or an inclusion with a different refractive index (Sun et al., Reference Sun, Yang, Kattawar and Zhang2017b). Most recently, Yang et al. (Reference Yang2019) extended the previous work to concave particles, particularly aggregates, and a more efficient computational algorithm was developed to map the near field on a cubic surface enclosing the particle to the far field based on an exact electromagnetic surface-integral equation. A synergistically unified method (Yang and Kattawar, Reference Yang and Kattawar2023) based on combining IITM (for small-to-moderate-size parameters) and PGOM (for moderate-to-large-size parameters) provides a practical approach for light-scattering computations covering the entire size-parameter range (Yang et al., Reference Yang2019). This method was recently used to develop databases of the optical properties of dust aerosols (Saito et al., Reference Saito2021) and ice clouds (Coy et al., Reference Coy2024). For computational efficiency a simplified form of the PGOM has been developed, which is referred to as the improved geometric optics method (IGOM) (Yang and Liou, Reference Yang and Liou1996b). The IGOM considers the interference of rays locally through accounting for the ray-spreading effect. The PGOM and IGOM obtain essentially the same results except in backscattering directions, but the latter is much more computationally efficient.

When a scattering particle is absorptive, the localized waves inside the particle may be inhomogeneous; that is, the planes of constant phase may not be in parallel with those of constant amplitude (Born and Wolf, Reference Born and Wolf1999). Early studies of the optical properties of absorptive particles highly simplified the treatment of inhomogeneous waves. In the past decade, significant progress has been made in this specific area of research (Dupertuis et al., Reference Dupertuis, Acklin and Proctor1994; Chang et al., Reference Chang, Walker and Hopcraft2005; Yang and Liou, Reference Yang and Liou2009a,Reference Yang and Lioub; Bi and Yang, Reference Bi and Yang2017). The ray-tracing technique involving inhomogeneous waves is quite complicated, which will be discussed in detail in Chapter 2.

The GOM does not consider the effect of the incident rays beyond the projected area of a scattering particle. This effect is referred to as the tunneling effect in the semiclassical-physics domain and has been extensively investigated in the case of a spherical particle by Nussenzveig and Wiscombe (Reference Nussenzveig and Wiscombe1980), Nussenzveig (Reference Nussenzveig1992), and Grandy (Reference Grandy2000). Their results have been extended to a spheroid (Bi and Yang, Reference Bi and Yang2014b). The tunneling effect is beyond the scope of the classical geometric optics; however, empirical approaches have been employed to incorporate the tunneling effect into the geometric optics solutions to compute the extinction and absorption coefficients for general nonspherical particles (Mitchell et al., Reference Mitchell2001; Yang et al., Reference Yang2013) within a unified and self-consistent framework.

In summary, many improvements to the geometric optics method for light-scattering computation have been made from various perspectives in the past three decades. Thus, it is beneficial to formally document this method by systematically incorporating all improvements.

1.2 Basic Physical Quantities of Quantifying Light Scattering

1.2.1 Maxwell’s Equations, Poynting Vector, and Electromagnetic Boundary Conditions

The most fundamental equations governing electromagnetic scattering are Maxwell’s equations. For simplicity, we consider Maxwell’s equations in the source-free form with Gaussian units for a non-absorbing medium (Jackson, Reference Jackson1975):

$\frac{1}{c}\frac{\partial \overrightarrow{D}}{\partial t}=\nabla \times \overrightarrow{H},$(1.2.1a)

$\frac{1}{c}\frac{\partial \overrightarrow{B}}{\partial t}=-\nabla \times \overrightarrow{E},$(1.2.1b)

$\nabla \cdot \overrightarrow{D}=0,$(1.2.1c)

$\nabla \cdot \overrightarrow{B}=0,$(1.2.1d)

where $\overrightarrow{E}$ is the electric field, $\overrightarrow{H}$ is the magnetic field, $\overrightarrow{D}$ is the electric displacement, and $\overrightarrow{B}$ is the magnetic induction. As standard notations, c indicates the speed of light in a vacuum and t indicates time. Note that $\overrightarrow{H}$ and $\overrightarrow{D}$ are derived quantities (Jackson, Reference Jackson1975) that, in the case of a linear and isotropic medium, are given by the following constitutive relations in the form

$\overrightarrow{D}=\epsilon \overrightarrow{E},$(1.2.2a)

$\overrightarrow{H}=\overrightarrow{B}/\mu .$(1.2.2b)

In Eqs. (1.2.2a) and (1.2.2b), $\epsilon =(1+4\pi \chi )$ and $\mu$ are the electric permittivity and magnetic permeability of the medium, respectively, where $\chi$ is the electric susceptibility. Quantity $\sqrt{\mu /\epsilon }$ is referred to as the impedance, and the backscattering by a sphere is sensitive to the impedance (Zhai et al., Reference Zhai2008) relative to that of the surrounding medium.

For two linear and isotropic media, Eqs. (1.2.1a)–(1.2.1d) lead to electromagnetic boundary conditions at the medium interface as follows:

$\widehat{n}\cdot \left({\overrightarrow{D}}_1-{\overrightarrow{D}}_2\right)=0,$(1.2.3a)

$\widehat{n}\cdot \left({\overrightarrow{B}}_1-{\overrightarrow{B}}_2\right)=0,$(1.2.3b)

$\widehat{n}\times \left({\overrightarrow{E}}_1-{\overrightarrow{E}}_2\right)=0,$(1.2.3c)

$\widehat{n}\times \left({\overrightarrow{H}}_1-{\overrightarrow{H}}_2\right)=0,$(1.2.3d)

where $\widehat{n}$ is a unit vector normal to the interface of the two media. Using Eqs. (1.2.1a) and (1.2.1b), we can obtain the following relation:

$\nabla \cdot \left(\overrightarrow{E}\times \overrightarrow{H}\right)=\left(\nabla \times \overrightarrow{E}\right)\cdot \overrightarrow{H}-\overrightarrow{E}\cdot \left(\nabla \times \overrightarrow{H}\right)$

$=\left(-\frac{1}{c}\frac{\partial \overrightarrow{B}}{\partial t}\overrightarrow{H}\cdot -\frac{1}{c}\frac{\partial \overrightarrow{D}}{\partial t}\cdot \overrightarrow{E}\right)$

$=-\frac{1}{2c}\frac{\partial }{\partial t}\left(\overrightarrow{B}\cdot \overrightarrow{H}+\overrightarrow{D}\cdot \overrightarrow{E}\right),$(1.2.4)

where we implicitly assume that the medium is homogeneous, isotropic, and linear. Based on Eq. (1.2.4), we can introduce the Poynting vector $(\overrightarrow{S})$ specifying electromagnetic energy flow and the electromagnetic energy density (u) as follows:

$\overrightarrow{S}=\frac{c}{4\pi }\left(\overrightarrow{E}\times \overrightarrow{H}\right),$(1.2.5)

$u=\frac{1}{8\pi }\left(\overrightarrow{B}\cdot \overrightarrow{H}+\overrightarrow{D}\cdot \overrightarrow{E}\right).$(1.2.6)

Thus, Eq. (1.2.4) is an expression of electromagnetic energy conservation in the form

$\nabla \cdot \overrightarrow{S}=-\frac{\partial u}{\partial t}.$((1.2.7))

To consider the transport of electromagnetic energy, let us compute ∇ × Eq. (1.2.1b), which leads to

$\nabla \times \nabla \times \overrightarrow{E}=-\frac{\mu }{c}\frac{\partial }{\partial t}\left(\nabla \times \overrightarrow{H}\right).$(1.2.8)

Expanding the left-hand side of Eq. (1.2.8) and substituting (1.2.1a) into the right-hand side of Eq. (1.2.8), we obtain

$\nabla \left(\nabla \cdot \overrightarrow{E}\right)-{\nabla }^2\overrightarrow{E}=-\frac{\mu \epsilon }{c^2}\frac{{\partial }^2\overrightarrow{E}}{\partial t^2}.$(1.2.9)

Because of $\nabla \cdot \overrightarrow{E}=0$ based on Eqs. (1.2.1c) and (1.2.2a), we obtain the wave equation for the electric field as follows:

${\nabla }^2\overrightarrow{E}-\frac{1}{{\left(c/\sqrt{\mu \epsilon }\right)}^2}\frac{{\partial }^2\overrightarrow{E}}{\partial t^2}=0.$(1.2.10)

Similarly, we can derive the wave equation for the magnetic field:

${\nabla }^2\overrightarrow{H}-\frac{1}{{\left(c/\sqrt{\mu \epsilon }\right)}^2}\frac{{\partial }^2\overrightarrow{H}}{\partial t^2}=0.$(1.2.11)

The quantity $\sqrt{\mu \epsilon }$ is referred to as the refractive index. We denote it by m.

Equations (1.2.10) and (1.2.11) together govern the propagation of an electromagnetic wave. In these two equations, $c\text{/}\sqrt{\mu \epsilon }$ is the phase speed of the wave in the medium. In this book, we only consider time harmonic waves and assume that the temporal dependence is in the form of $e^{-i\omega t}$, where $i=\sqrt{-1}$ and $\omega$ is the angular frequency. This choice is consistent with the counterpart in many texts such as Born and Wolf (Reference Born and Wolf1999), Bohren and Huffman (Reference Bohren and Huffman1983), and Mishchenko et al. (Reference Mishchenko, Travis and Lacis2002). With the temporal-dependence factor $e^{-i\omega t}$, the imaginary part of the complex index of refraction in the frequency domain is non-negative (i.e., the quantity can be zero or a positive number) for an absorbing medium. Alternatively, as an equally valid approach, the temporal-dependence factor can be specified as $e^{i\omega t}$, such as in van de Hulst (Reference van de Hulst1957) and Liou (Reference Liou2002). The latter choice leads to a non-positive (i.e., the quantity can be zero or a negative number) imaginary part of the complex index of refraction in the frequency domain for an absorbing medium. Whether we select $e^{-i\omega t}$ or $e^{i\omega t}$, the final analysis gives the same measurable quantities. However, a specific choice of the temporal factor must be used in a consistent form in downstream derivations and formulae as pointed out by Bohren and Huffman (Reference Bohren and Huffman1983) and Mishchenko et al. (Reference Mishchenko, Travis and Lacis2002). Throughout this text, we will use the temporal-dependence factor $e^{-i\omega t}$ to specify a time harmonic wave; therefore, the electric and magnetic fields can be expressed in the form

$\overrightarrow{E}\left(\overrightarrow{r},t\right)=\overrightarrow{E}\left(\overrightarrow{r}\right)e^{-i\omega t},\overrightarrow{H}\left(\overrightarrow{r},t\right)=\overrightarrow{H}\left(\overrightarrow{r}\right)e^{-i\omega t},$(1.2.12)

where $\overrightarrow{r}$ indicates the position vector. The light-scattering processes involved in this book are constrained to the case of elastic scattering; that is, the electromagnetic wave frequency is not changed in the scattering process.

With the expression in Eq. (1.2.12), we can simplify Maxwell’s equations in their time-independent form as follows:

$-ik\sqrt{\frac{\epsilon }{\mu }}\overrightarrow{E}\left(\overrightarrow{r}\right)=\nabla \times \overrightarrow{H}\left(\overrightarrow{r}\right),$(1.2.13a)

$ik\overrightarrow{H}\left(\overrightarrow{r}\right)=\sqrt{\frac{\epsilon }{\mu }}\nabla \times \overrightarrow{E}\left(\overrightarrow{r}\right),$(1.2.13b)

$\nabla \cdot \overrightarrow{E}\left(\overrightarrow{r}\right)=0,$(1.2.13c)

$\nabla \cdot \overrightarrow{H}\left(\overrightarrow{r}\right)=0,$(1.2.13d)

where $k=\omega \sqrt{\epsilon \mu }/c$ is referred to as the modified wavenumber in a medium with the index of refraction $m=\sqrt{\epsilon \mu }$.

For a time-harmonic wave in optical analysis, we normally consider the average of a measurable quantity within a temporal duration T much longer than $1/\omega$. Thus, the time-averaged Poynting vector is given by

$〈\overrightarrow{S}〉=\frac{1}{T}{\displaystyle {\int }_0^T\frac{c}{4\pi }}\left\{\mathrm{Re}\left[\overrightarrow{E}\left(\overrightarrow{r}\right)e^{-\omega t}\right]\times \mathrm{Re}\left[\overrightarrow{H}\left(\overrightarrow{r}\right)e^{-\omega t}\right]\right\}dt$

$=\frac{c}{16\pi T}{\displaystyle {\int }_0^T\left\{\left[\overrightarrow{E}\left(\overrightarrow{r}\right)e^{-i\omega t}+{\overrightarrow{E}}^{*}\left(\overrightarrow{r}\right)e^{i\omega t}\right]\times \left[\overrightarrow{H}\left(\overrightarrow{r}\right)e^{-i\omega t}+{\overrightarrow{H}}^{*}\left(\overrightarrow{r}\right)e^{i\omega t}\right]\right\}}dt,$(1.2.14)

where the asterisk indicates complex conjugation. Note that because $\frac{1}{T}{\displaystyle {\int }_0^T{e}^{\pm i2\omega t}dt\to 0}$ if $T\gg 1\text{/}\omega$, Eq. (1.2.14) can be simplified in the form

$〈\overrightarrow{S}〉=\frac{c}{16\pi }\left[\overrightarrow{E}\left(\overrightarrow{r}\right)\times {\overrightarrow{H}}^{*}\left(\overrightarrow{r}\right)+{\overrightarrow{E}}^{*}\left(\overrightarrow{r}\right)\times \overrightarrow{H}\left(\overrightarrow{r}\right)\right]$

$=\frac{c}{8\pi }\mathrm{Re}\left[\overrightarrow{E}\left(\overrightarrow{r}\right)\times {\overrightarrow{H}}^{*}\left(\overrightarrow{r}\right)\right].$(1.2.15)

Thus, we can now define the complex Poynting vector, where only the real part is physically meaningful, as:

$\overrightarrow{S}=\frac{c}{8\pi }\overrightarrow{E}\left(\overrightarrow{r}\right)\times {\overrightarrow{H}}^{*}\left(\overrightarrow{r}\right).$(1.2.16)

Similarly, the harmonic electromagnetic energy density is given by

$u=\frac{1}{16\pi }\left(\overrightarrow{E}\cdot {\overrightarrow{D}}^{*}+\overrightarrow{H}\cdot {\overrightarrow{B}}^{*}\right).$(1.2.17)

Furthermore, for a homogeneous (i.e., the planes of constant phase are parallel to the planes of constant amplitude) plane electromagnetic wave propagating along a direction indicated by a unit vector $\widehat{k}$ with the electric field in the form of

$\overrightarrow{E}\left(\overrightarrow{r}\right)={\overrightarrow{E}}_0{e}^{i\overrightarrow{k}\cdot \overrightarrow{r}},$(1.2.18a)

$\overrightarrow{k}=\frac{\omega }{c}\sqrt{\mu \epsilon }\widehat{k},$(1.2.18b)

where ${\overrightarrow{E}}_0$ is the complex amplitude of the wave, and $\widehat{k}$ is a unit vector pointing along vector $\overrightarrow{k}$, that is, $\widehat{k}=\frac{\overrightarrow{k}}{\left|\overrightarrow{k}\right|}$.

If we select $e^{-i\omega t}$ as the time-dependent factor for a plane electromagnetic wave in an absorbing (with a nonzero effective electric conductivity $\sigma \ne 0$) nonferromagnetic, medium ($\mu =1$), the complex refractive index in the frequency domain is given by $m=m_r+i{m}_i={({\epsilon }_r+i{\epsilon }_i)}^{1/2}={\left[(1+4\pi \chi )+i4\pi \sigma /\omega \right]}^{1/2}$, or ${\epsilon }_r=1+4\pi \chi =m_r^2-m_i^2$, and ${\epsilon }_i=4\pi \sigma /\omega =2m_r{m}_i$. In this case, the time-averaged Poynting vector is given by

$〈\overrightarrow{S}\left(\overrightarrow{r}\right)〉=\frac{c}{8\pi }\mathrm{Re}\left[\overrightarrow{E}\left(\overrightarrow{r}\right)\times {\overrightarrow{H}}^{*}\left(\overrightarrow{r}\right)\right]$

$=\frac{c}{8\pi }{m}_r{\left|{\overrightarrow{E}}_0\right|}^2\text{exp}\left[-2\left(\omega /c\right)m_{i}k̂\cdot \overrightarrow{r}\right]k̂.$(1.2.19)

It is clear from the above equation Eq. (1.2.19) that the magnitude of the Poynting vector decreases with the electromagnetic wave propagation if $m_i$ is nonzero, that is, $\sigma \ne 0$.

1.2.2 Stokes Parameters

An electromagnetic wave is a vector wave; that is, its electric and magnetic fields must be specified in terms of both their magnitudes and directions. The vibration configuration of the electric field vector is referred to as the polarization state of the wave. Perrin (Reference Perrin1942) presented a thorough review of the polarization state. The polarization state is referred to as linear, circular, or elliptical if the projection of the end point of the electric vector on a plane transverse to the propagation direction traces out a straight line, a circle, or an ellipse, respectively, as schematically illustrated in Figure 1.6. In the case of circular or elliptical polarization, the polarization state is said to be right-handed or left-handed for an observer facing the incoming light (i.e., the observer looks toward the light source) if the end point of the electric field vector describes a helix in a clockwise or anticlockwise form, respectively (Born and Wolf, Reference Born and Wolf1999). The opposite convention has also been used in the literature by defining the observation to be looking along the light propagation direction (van de Hulst, Reference van de Hulst1957; Hovenier et al., Reference Hovenier, van der Mee and Domke2004). In this text, the convention of right-/left-handed polarization follows the traditional definition described in Born and Wolf (Reference Born and Wolf1999) by looking toward the light source, which is also consistent with Bohren and Huffman (Reference Bohren and Huffman1983) and Mishchenko et al. (Reference Mishchenko, Travis and Lacis2002).

Figure 1.6

Illustration of the polarization states of linearly, circularly, and elliptically polarized waves. The waves propagate along the direction specified by unit vector $\widehat{k}.$

Most optical instruments do not directly measure the electric and magnetic fields, mainly because the fields are highly oscillatory complex quantities. To overcome this difficulty, four real-valued optically measurable quantities referred to as the Stokes parameters, introduced originally by Stokes (Reference Stokes1862), can be used. To define the Stokes parameters, let us consider a plane electromagnetic wave in a vacuum as follows:

$\overrightarrow{E}\left(\overrightarrow{r},t\right)=\left(E_{\parallel }e{̂}_{\parallel }+E_{\perp }e{̂}_{\perp }\right)e^{i\left(kk̂\cdot \overrightarrow{r}-\omega t\right)},$(1.2.20)

where ${\widehat{e}}_{\parallel }$ and ${\widehat{e}}_{\perp }$ are two orthogonal unit vectors, and unit vector $\widehat{k}={\widehat{e}}_{\perp }\times {\widehat{e}}_{\parallel }$ points along the propagation direction of the wave, as illustrated in Figure 1.7. The quantities $E_{\parallel }$ and $E_{\perp }$ are the projections of the electric field vector onto the directions specified by ${\widehat{e}}_{\parallel }$ and ${\widehat{e}}_{\perp }$, respectively. Note that $E_{\parallel }$ and $E_{\perp }$ can be complex quantities and are often referred to as the parallel and perpendicular electric field components with respect to the meridional plane specified in terms of zenith angle $\theta$ and azimuthal angle $\phi$, as shown in Figure 1.7. With such a decomposition of the electric field vector associated with a wave in a nonabsorptive medium with relative permittivity $\epsilon$ and relative permeability $\mu$, the Stokes parameters can be defined in the form

$I=\frac{c}{8\pi }\sqrt{\frac{\epsilon }{\mu }}\left(E_{\parallel }{E}_{\parallel }^{*}+E_{\perp }{E}_{\perp }^{*}\right),$(1.2.21a)

$Q=\frac{c}{8\pi }\sqrt{\frac{\epsilon }{\mu }}\left(E_{\parallel }{E}_{\parallel }^{*}-E_{\perp }{E}_{\perp }^{*}\right),$(1.2.21b)

$U=\frac{c}{8\pi }\sqrt{\frac{\epsilon }{\mu }}\left(E_{\parallel }{E}_{\perp }^{*}+E_{\parallel }^{*}{E}_{\perp }\right),$(1.2.21c)

$V=i\frac{c}{8\pi }\sqrt{\frac{\epsilon }{\mu }}\left(E_{\parallel }{E}_{\perp }^{*}-E_{\parallel }^{*}{E}_{\perp }\right),$(1.2.21d)

where the asterisk indicates complex conjugation and $i=\sqrt{-1}$. The Stokes parameters have the dimension of irradiance, namely power per unit area. Note that the convention of the Stokes parameters in Eqs. (1.2.21a)–(1.2.21d) is consistent with many texts such as van de Hulst (Reference van de Hulst1957), Bohren and Huffman (Reference Bohren and Huffman1983), and Mishchenko et al. (Reference Mishchenko, Travis and Lacis2002). However, Liou (Reference Liou2002) introduces a negative sign to component V. Moreover, the sign difference of the V component may cause confusion in the definition of left-handed and right-handed circular polarization. Note that the I, Q, and U components are not affected by the choice of the time dependence, but the sign of the V component depends on both the sign choice in Eq. (1.2.21d) and the representation of harmonic waves in the form of either $e^{-i\omega t}$ or $e^{i\omega t}$. Because the Stokes parameters are usually considered in a relative sense, the constant factor $(c/8\pi )\sqrt{\epsilon /\mu }$ in Eqs. (1.2.21a)–(1.2.21d) can be omitted in many cases without affecting the relevant analysis. Thus, following Bohren and Huffman (Reference Bohren and Huffman1983) and Mishchenko et al. (Reference Mishchenko, Travis and Lacis2002), we will not include this constant factor in the following discussions.

Figure 1.7

Reference plane for defining the Stokes parameters.

For a monochromatic harmonic wave, it will be shown later that the four Stokes parameters are not independent, and they satisfy the following relation:

$I^2=Q^2+U^2+V^2.$(1.2.22)

The Stokes parameters fully specify the characteristics of a light beam; in particular, any optical analysis is not capable of distinguishing two light beams that have the same Stokes parameters (Chandrasekhar, Reference Chandrasekhar1960). Among the Stokes parameters, the quantity I specifies the intensity of light, the quantities U and V specify the linear polarization, and V specifies circular polarization, as illustrated by Figure 1.8. Thus, it is natural to introduce the degree of polarization (DoP):

$DoP={(Q^2+U^2+V^2)}^{1/2}/I,$(1.2.23)

the degree of circular polarization (DoCP) as

$DoCP=V/I,$(1.2.24)

and the degree of linear polarization (DoLP) as

$DoLP={(Q^2+U^2)}^{1/2}/I.$(1.2.25)

In the case of $U=0$, DoLP is usually defined as follows (Hovenier et al., Reference Hovenier, van der Mee and Domke2004):

$DoLP=-Q/I.$(1.2.26)

Figure 1.8

Polarization configurations for linearly polarized light (panels a, b, d, and e) and circularly polarized light (panels c and f). Note that the light propagation direction points out of the paper (the observer looks toward the light source). Panel (c) illustrates right circular polarization and panel (f) illustrates left circular polarization.

Let us assume the phases of the parallel and perpendicular electric field components to be ${\alpha }_{\parallel }$ and ${\alpha }_{\perp }$, and the real amplitudes of the two electric field components to be $A_{\parallel }$ and $A_{\perp }$. Then, we can rewrite Eq. (1.2.20) in the form

$\overrightarrow{E}\left(\overrightarrow{r},t\right)=\left(A_{\parallel }{e}^{-i{\alpha }_{\parallel }}e{̂}_{\parallel }+A_{\perp }{e}^{-i{\alpha }_{\perp }}e{̂}_{\perp }\right)e^{i\left(kk̂\cdot \overrightarrow{r}-\omega t\right)}.$(1.2.27)

With the representation of the electric field in Eq. (1.2.27), along with omitting the constant factor $(c/8\pi )\sqrt{\epsilon /\mu }$, the Stokes parameters can be given in the form

$I=A_{\parallel }^2+A_{\perp }^2,$(1.2.28a)

$Q=A_{\parallel }^2-A_{\perp }^2,$(1.2.28b)

$U=2A_{\parallel }{A}_{\perp }\cos \left({\alpha }_{\parallel }-{\alpha }_{\perp }\right),$(1.2.28c)

$V=2A_{\parallel }{A}_{\perp }\sin \left({\alpha }_{\parallel }-{\alpha }_{\perp }\right).$(1.2.28d)

The most general polarization state is elliptical polarization, where the projection of the end point of the electric vector on a plane transverse to the propagation direction traces out an ellipse. If we define two unit vectors $\widehat{u}$ and $\widehat{v}$ to be along the semi-major and semi-minor axes of the ellipse as shown in Figure 1.9, the real part of the electric vector can be expressed in the form

$\mathrm{Re}\left[\overrightarrow{E}\left(r,t\right)\right]=ûA_u\cos \left(kr-\omega t+\eta \right)+v̂A_v\sin \left(kr-\omega t+\eta \right),$(1.2.29)

where $\eta$ denotes a constant phase difference. Because the intensity of the light beam is independent of the coordinate system employed to specify the polarization state, the following relation must hold:

$A_u^2+A_v^2=A_{\perp }^2+A_{\parallel }^2=A^2.$(1.2.30)

Figure 1.9

Geometric configuration of elliptical polarization. Note that this diagram depicts quantities in the $\widehat{u}-\widehat{v}$ plane. Angles $\beta$ and $\gamma$ are defined in the text.

Because of the relation $A_u^2+A_v^2=A^2$, we introduce an angle $\beta \in [-\pi /4,\pi /4]$ such that

$A_u=A\cos \beta ,$(1.2.31a)

$A_v=A\sin \beta .$(1.2.31b)

The sign of $\beta$ is defined positive for right-handed polarization and negative for left-handed polarization. The geometric meaning of angle $\beta$ is illustrated in Figure 1.9, that is, $\tan \beta$ indicates the ellipticity of the ellipse, namely, $A_v/A_u$

To project the real part of the electric vector in Eq. (1.2.29) onto the two directions specified by ${\widehat{e}}_{\parallel }$ and ${\widehat{e}}_{\perp }$, we obtain

$\begin{array}{l}{A}_{\parallel }\cos \left(kr-\omega t-{\alpha }_{\parallel }\right)=A\cos \beta \cos (kr-\omega t+\eta )\cos \gamma \\ +A\sin \beta \sin (kr-\omega t+\eta )\sin \gamma ,\end{array}$(1.2.32a)

$\begin{array}{l}{A}_{\perp }\cos \left(kr-\omega t-{\alpha }_{\perp }\right)=A\cos \beta \cos (kr-\omega t+\eta )\sin \gamma \\ -A\sin \beta \sin (kr-\omega t+\eta )\cos \gamma .\end{array}$(1.2.32b)

The above equations Eqs. (1.2.32a) and (1.2.32b) can be further expanded and rearranged as follows:

$\begin{array}{l}\cos \left(kr-\omega t\right)\left(A_{\parallel }\cos {\alpha }_{\parallel }-A\cos \beta \cos \gamma \cos \eta -A\sin \beta \sin \gamma \sin \eta \right)=\hfill \\ \sin \left(kr-\omega t\right)\left(-A_{\parallel }\sin {\alpha }_{\parallel }-A\cos \beta \cos \gamma \sin \eta +A\sin \beta \sin \gamma \cos \eta \right),\hfill \end{array}$(1.2.33a)

$\begin{array}{l}\cos \left(kr-\omega t\right)\left(A_{\perp }\cos {\alpha }_{\perp }-A\cos \beta \sin \gamma \cos \eta +A\sin \beta \cos \gamma \sin \eta \right)=\hfill \\ \sin \left(kr-\omega t\right)\left(-A_{\perp }\sin {\alpha }_{\perp }-A\cos \beta \sin \gamma \sin \eta -A\sin \beta \cos \gamma \cos \eta \right).\hfill \end{array}$(1.2.33b)

Because $\cos (kr-\omega t)$ and $\sin (kr-\omega t)$ are linearly independent, Eqs. (1.2.33a) and (1.2.33b) lead to

$A_{\parallel }\cos {\alpha }_{\parallel }=A\cos \beta \cos \gamma \cos \eta +A\sin \beta \sin \gamma \sin \eta ,$(1.2.34a)

$A_{\parallel }\sin {\alpha }_{\parallel }=-A\cos \beta \cos \gamma \sin \eta +A\sin \beta \sin \gamma \cos \eta ,$(1.2.34b)

$A_{\perp }\cos {\alpha }_{\perp }=A\cos \beta \sin \gamma \cos \eta -A\sin \beta \cos \gamma \sin \eta ,$(1.2.34c)

$A_{\perp }\sin {\alpha }_{\perp }=-A\cos \beta \sin \gamma \sin \eta -A\sin \beta \cos \gamma \cos \eta .$(1.2.34d)

Based on Eqs. (1.2.34a)–(1.2.34d), the Stokes parameters can be written with respect to the amplitude and the orientation angle as depicted in Figure 1.9. On squaring and adding Eqs. (1.2.34a) and (1.2.34b), we obtain the intensity along the ${\widehat{e}}_{\parallel }$ direction

$A_{\parallel }^2=A^2\left({\cos }^2\beta {\cos }^2\gamma +{\sin }^2\beta {\sin }^2\gamma \right).$(1.2.35)

Similarly, Eqs. (1.2.34c) and (1.2.34d) lead to

$A_{\perp }^2=A^2\left({\cos }^2\beta {\sin }^2\gamma +{\sin }^2\beta {\cos }^2\gamma \right).$(1.2.36)

Adding the product of Eqs. (1.2.34a) and (1.2.34c) to that of (1.2.34b) and (1.2.34d), we obtain

$2A_{\parallel }{A}_{\perp }\cos \alpha =A^2\sin 2\gamma \cos 2\beta ,$(1.2.37)

where $\alpha ={\alpha }_{\parallel }-{\alpha }_{\perp }$. Subtracting the product of Eqs. (1.2.34b) and (1.2.34c) from that of (1.2.34a) and (1.2.34d) yields

$2A_{\parallel }{A}_{\perp }\sin \alpha =A^2\sin 2\beta .$(1.2.38)

Taking the quotient of Eqs. (1.2.34b) and (1.2.34a), we have

$\tan \left({\alpha }_{\parallel }+\eta \right)=\tan \beta \tan \gamma .$(1.2.39)

Similarly, the quotient of Eqs. (1.2.34d) and (1.2.34c) gives

$\tan \left({\alpha }_{\perp }+\eta \right)=-\tan \beta \cot \gamma ,$(1.2.40)

therefore

$\tan \alpha =\tan \left[\left({\alpha }_{\parallel }+\eta \right)-\left({\alpha }_{\perp }+\eta \right)\right]=\frac{\tan 2\beta }{\sin 2\gamma }.$(1.2.41)

Note that sines and cosines of the angles $\beta$, $\gamma$, and $\eta$ can be simultaneously switched in Eqs. (1.2.32)–(1.2.34), to represent the geometry of the polarization. In that case, Eqs. (1.2.35)–(1.2.38) and (1.2.41) will not be changed, but Eqs. (1.2.39) and (1.2.40) will be given by

$\tan \left(\eta -{\alpha }_{\parallel }\right)=\cot \beta \cot \gamma ,$(1.2.42)

$\tan \left(\eta -{\alpha }_{\perp }\right)=-\cot \beta \tan \gamma .$(1.2.43)

With Eqs. (1.2.28), (1.2.35)–(1.2.38), it is straightforward to show that the Stokes parameters in geometric notation can be expressed in the form (van de Hulst, Reference van de Hulst1957)

$I=A^2,$((1.2.44a))

$Q=A^2\cos 2\beta \cos 2\gamma ,$(1.2.44b)

$U=A^2\cos 2\beta \sin 2\gamma ,$(1.2.44c)

$V=A^2\sin 2\beta .$(1.2.44d)

The four Stokes parameters associated with the electric field specified by Eq. (1.2.27) are not completely independent, as they depend on only three quantities, namely, A, $\gamma$, and $\beta$ in Eqs. (1.2.44a)–(1.2.44d). Based on Eq. (1.2.44), it is straightforward to show that $I^2=Q^2+U^2+V^2$. From Eq. (1.2.44), the orientation of the polarization ellipse and the ellipticity can be specified by

$\tan 2\gamma =U/Q.$(1.2.45a)

$\sin 2\beta =V/{(Q^2+U^2+V^2)}^{1/2}.$(1.2.45b)

Furthermore, it is obvious that the relations in Eqs. (1.2.44a)–(1.2.44d) imply that Q, U, and V can be used as spherical coordinates. Thus, the Stokes parameters can be represented using a spherical coordinate system, the so-called Poincaré sphere as illustrated in Figure 1.10, following Born and Wolf (Reference Born and Wolf1999).

Figure 1.10

Representation of the Stokes parameters on the Poincaré sphere.

The Stokes parameters are specified with respect to a specific reference plane. As shown in Figure 1.11, the decomposition of the electric field vector associated with the rotation of a coordinate system is

${E^{\prime }}_{\parallel }=E_{\parallel }\cos \varphi -E_{\perp }\sin \varphi ,$(1.2.46a)

${E^{\prime }}_{\perp }=E_{\parallel }\sin \varphi +E_{\perp }\cos \varphi ,$(1.2.46b)

where $\varphi$ is defined as the clockwise rotation angle of the ${\widehat{e}}_{\perp }$ axis to the ${{\widehat{e}}^{\prime }}_{\perp }$ axis for an observer looking along the propagation direction specified by unit vector $\widehat{k}$. The transformation represented by Eqs. (1.2.46a) and (1.2.46b) can be expressed in a matrix form as

$\left[\begin{array}{c}{E}_{\parallel }^{\prime }\\ E_{\perp }^{\prime }\end{array}\right]=\left[\begin{array}{cc}\cos \varphi & -\sin \varphi \\ \sin \varphi & \cos \varphi \end{array}\right]\left[\begin{array}{c}{E}_{\parallel }\\ E_{\perp }\end{array}\right].$(1.2.47)

Figure 1.11

Rotation of a reference coordinate system for the electric field decomposition.

As a result of the field transform in Eqs. (1.2.46a) and (1.2.46b), the rotational transform for the Stokes parameters is given by

$\left[\begin{array}{c}I\prime \\ Q\prime \\ U\prime \\ V\prime \end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos 2\varphi & \sin 2\varphi & 0\\ 0& -\sin 2\varphi & \cos 2\varphi & 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}I\\ Q\\ U\\ V\end{array}\right].$(1.2.48)

It can be shown that, with the rotational transform given by Eq. (1.2.48), $I$, $Q^2+U^2$, and $V$ are invariant. In addition, the degree of polarization defined by Eq. (1.2.23), the degree of circular polarization defined by Eq. (1.2.24), the degree of linear polarization defined by Eq. (1.2.25), and the ellipticity specified in terms of $\tan \beta$ (see Figure 1.9) are also invariant.

Later we recapture the experimental methodology described by Chandrasekhar (Reference Chandrasekhar1960) and Bohren and Huffman (Reference Bohren and Huffman1983) for measuring the Stokes parameters. Specifically, Figure 1.12 schematically illustrates an experimental setup for measuring the Stokes parameters. Specifically, a wave retarder introduces a certain amount of retardation in the phase of one of the two components of the electric field decomposed with respect to the ${\widehat{e}}_{\parallel }$ and ${\widehat{e}}_{\perp }$ directions. In the following discussion, we assume the phase difference between ${\widehat{e}}_{\parallel }$ and ${\widehat{e}}_{\perp }$ components of the electric field vector to be $\alpha$. The polarizer constrains the measurement of the electric field vector component vibrating in only one direction and deviating from the ${\widehat{e}}_{\parallel }$ axis by an angle $\phi$, which is called the direction of the transmission axis of the linear polarizer.

Figure 1.12

Schematic diagram illustrating the measurement of the Stokes parameters.

Following Chandrasekhar (Reference Chandrasekhar1960), the measured radiance can be expressed as a function of angles $\alpha$ and $\phi$ in the form of

$I\left(\alpha ,\phi \right)=I_{\parallel }{\cos }^2\phi +I_{\perp }{\sin }^2\phi +\frac{1}{2}\left(U\cos \alpha -V\sin \alpha \right)\sin 2\phi ,$(1.2.49)

where $I_{\parallel }=|E_{\parallel }{|}^2$ and $I_{\perp }=|E_{\perp }{|}^2$. Using Eq. (1.2.49), let us consider the following four polarization cases:

Case 1: Natural light, that is, unpolarized with $I_{\parallel }=I_{\perp }=I/2$, $U=0$, and $V=0$. Thus, we have

$I=2I(\alpha ,\phi ).$(1.2.50)

Case 2: Horizontal and vertical polarization configurations without a retarder for a polarized incident beam ($\alpha =0$ and two measurements with $\phi =0°$ and $\phi =90°$). In this scenario, we have

$I\left(\phi =0°,\alpha =0\right)=I_{\parallel },$(1.2.51a)

$I\left(\phi =90°,\alpha =0\right)=I_{\perp }.$(1.2.51b)

Thus, the first two Stokes parameters can be given by

$I=I_{\parallel }+I_{\perp }=I\left(\phi =0°,\alpha =0\right)+I\left(\phi =90°,\alpha =0\right),$(1.2.52a)

$Q=I_{\parallel }-I_{\perp }=I\left(\phi =0°,\alpha =0\right)-I\left(\phi =90°,\alpha =0\right).$(1.2.52b)

Case 3: +45° and −45° polarization configurations without a compensator (for a polarized incident beam), $\alpha =0$ and two measurements with $\phi =45°$ and $\phi =-45°$. In this scenario, we have

$I\left(\phi =45°,\alpha =0\right)=\frac{1}{2}\left(I_{\parallel }+I_{\perp }\right)+\frac{1}{2}U,$((1.2.53a))

$I\left(\phi =-45°,\alpha =0\right)=\frac{1}{2}\left(I_{\parallel }+I_{\perp }\right)-\frac{1}{2}U.$(1.2.53b)

From Eq. (1.2.53), we obtain

$U=I\left(\phi =45°,\alpha =0\right)-I\left(\phi =-45°,\alpha =0\right).$(1.2.54)

Case 4: Measurement of circular polarization using +45° and −45° polarization configurations with a quarter-wave plate ($\alpha =\pi /2$ and two measurements with $\phi =45°$ and $\phi =-45°$). In this scenario, we have

$I\left(\phi =45°,\alpha =\pi /2\right)=\frac{1}{2}\left(I_{\parallel }+I_{\perp }\right)-\frac{1}{2}V,$(1.2.55a)

$I\left(\phi =-45°,\alpha =\pi /2\right)=\frac{1}{2}\left(I_{\parallel }+I_{\perp }\right)+\frac{1}{2}V.$(1.2.55b)

From Eq. (1.2.55), we obtain

$V=I\left(\phi =-45°,\alpha =\pi /2\right)-I\left(\phi =45°,\alpha =\pi /2\right).$(1.2.56)

For an arbitrarily polarized light beam, the measurements associated with Eqs. (1.2.51), (1.2.52), (1.2.54), and (1.2.56) can fully quantify the Stokes parameters of the light beam.

1.2.3 Amplitude Scattering Matrix and the Optical Theorem

Figure 1.13 illustrates the general configuration of the scattering of light by a particle. The observer is far away from the scattering particle, and thus, the observed scattered field is in the far-field zone (i.e., the distance between the observer and the particle is much larger than the wavelength of the incident light). The direction pointing from the particle to the observer is referred to as the scattering direction, and the scattering angle is the angle between the incident and scattering directions. The plane containing the incident and scattering directions is called the scattering plane. If the scattering particle has axial symmetry around the z-axis, the scattering properties are independent of the rotation of the scattering plane around the symmetric axis of the particle, that is, the z-axis in the present case. In the most general case, the scattering properties depend on the rotation of the scattering plane around the z-axis.

Figure 1.13

Schematic diagram illustrating the configuration of the scattering of light by a particle, where superscripts i and s refer to the incident and scattering directions, respectively.

As shown in Figure 1.13, both incident and scattered electric vectors can be decomposed into parallel and perpendicular components with respect to the scattering plane. With the decomposition, we have the following relation

$\left(\begin{array}{c}{E}_{\parallel }^s\\ E_{\perp }^s\end{array}\right)=\frac{\exp \left[ik\left(r-z\right)\right]}{-ikr}\left(\begin{array}{cc}{S}_2& S_3\\ S_4& S_1\end{array}\right)\left(\begin{array}{c}{E}_{\parallel }^i\\ E_{\perp }^i\end{array}\right),$(1.2.57)

where k is referred to as the modified wavenumber in a vacuum, given by $k=2\pi /\lambda$ where λ is the wavelength. r (>> λ) is the distance between the observer and the scattering particle, and $S_{i,i=1-4}$ are the elements of the amplitude scattering matrix. Here we use the popular notation introduced originally by van de Hulst (Reference van de Hulst1957) for the amplitude scattering matrix, namely,

$S=\left[\begin{array}{cc}{S}_2& S_3\\ S_4& S_1\end{array}\right].$(1.2.58)

The amplitude scattering matrix contains essentially all the optical information about the elastic scattering and absorption properties of the scattering particle. In particular, the amplitude scattering matrix elements in the forward direction are directly related to the extinction cross section of the particle. To illustrate this feature, we consider the far-field region near the forward direction $\left[(x^2+y^2)/z^2\ll 1\right]$, where the following approximation is valid:

$R={(x^2+y^2+z^2)}^{1/2}\approx z+(x^2+y^2)/(2z)+O(1/z^2).$(1.2.59)

Thus, the scattered field in the forward direction can be given by

$\left(\begin{array}{c}{E}_{\parallel }^s\\ E_{\perp }^s\end{array}\right)\approx \frac{\exp \left[ik\left(x^2+y^2\right)/\left(2z\right)\right]}{-ikz}\left(\begin{array}{c}{S}_2\left(0°\right)E_{\parallel }^i+S_3\left(0°\right)E_{\perp }^i\\ S_4\left(0°\right)E_{\parallel }^i+S_1\left(0°\right)E_{\perp }^i\end{array}\right).$(1.2.60)

The first Stokes parameter associated with the total electromagnetic field (the incident field plus the scattered field) can be given in the form

$I^t=\frac{c}{8\pi }\left[\left(E_{\parallel }^i+E_{\parallel }^s\right)\left(E_{\parallel }^{i*}+E_{\parallel }^{s*}\right)+\left(E_{\perp }^i+E_{\perp }^s\right)\left(E_{\perp }^{i*}+E_{\perp }^{s*}\right)\right],$(1.2.61)

where superscript t indicates the total irradiance. The earlier irradiance can be decomposed into two components:

$I^t=I^i-I^e+O(1/z^2),$(1.2.62)

where $I^i$ is the incident irradiance and $I^e$ is associated with the extinction process. These components are given by

$I^i=\frac{c}{8\pi }\left(E_{\parallel }^i{E}_{\parallel }^{i*}+E_{\perp }^i{E}_{\perp }^{i*}\right),$(1.2.63)

$\begin{array}{l}{I}^e=\frac{c}{4\pi }\mathrm{Re}\left\{\frac{\exp \left[ik\left(x^2+y^2\right)/\left(2z\right)\right]}{-ikz}\right.\\ \left.\left[S_2\left(0°\right)E_{\parallel }^i{E}_{\parallel }^{i*}+S_3\left(0°\right)E_{\perp }^i{E}_{\parallel }^{i*}+S_4\left(0°\right)E_{\parallel }^i{E}_{\perp }^{i*}+S_1\left(0°\right)E_{\perp }^i{E}_{\perp }^{i*}\right]\right\}\end{array}$

$\begin{array}{l}=\frac{c}{4\pi }\mathrm{Re}\left\{\frac{\exp \left[ik\left(x^2+y^2\right)/\left(2z\right)\right]}{-ikz}\left(I^i\left[S_2\left(0°\right)+S_1\left(0°\right)\right]+Q^i\left[S_2\left(0°\right)-S_1\left(0°\right)\right]\right)\right.\\ \left.+U^i\left[S_3\left(0°\right)+S_4\left(0°\right)\right]+i{V}^i\left[S_3\left(0°\right)-S_4\left(0°\right)\right]\right\}.\end{array}$(1.2.64)

The extinction cross section of the particle can be calculated via

$C_{ext}=\frac{1}{I^i}{\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }{I}^{e}dxdy}}.$(1.2.65)

From Eq. (1.2.65) and the double Fresnel integral given by

${\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }{e}^{i\left(x^2+y^2\right)}dxdy=i\pi }},$(1.2.66)

it is straightforward to formulate the extinction cross section as follows:

$\begin{array}{l}{C}_{ext}=\frac{2\pi }{k^2}\mathrm{Re}\left\{\left[S_2\left(0°\right)+S_1\left(0°\right)\left]+\right[S_2\left(0°\right)-S_1\left(0°\right)\right]\left(Q^i/I^i\right)\right.\\ \left.+\left[S_3\left(0°\right)+S_4\left(0°\right)\left](U^i/I^i)+i\right[S_3\left(0°\right)-S_4\left(0°\right)\right](V^i/I^i)\right\}.\end{array}$(1.2.67)

It is evident from Eq. (1.2.67) that the extinction cross section is determined by the amplitude scattering matrix in the forward direction and by the polarization state of the incident electromagnetic radiation, although scattering occurs in all directions and the extinction cross section is the sum of the scattering and absorption cross sections. If the incident radiation is unpolarized $(Q=U=V=0)$, the extinction cross section given by Eq. (1.2.67) reduces to

$C_{ext}=\frac{2\pi }{k^2}\mathrm{Re}\left[S_2\left(0°\right)+S_1\left(0°\right)\right].$(1.2.68)

If $S_2(0°)$ and $S_2(0°)$ are the same [hereafter, indicated by $S(0°)$], such as in the case of the scattering of light by a sphere, Eq. (1.2.68) reduces to the form of the optical theorem given in most texts:

$C_{ext}=\frac{4\pi }{k^2}\mathrm{Re}\left[S\left(0°\right)\right].$(1.2.69)

Note that the approximation ${(x^2+y^2+z^2)}^{1/2}\approx z+(x^2+y^2)/(2z)$ assumed in Eq. (1.2.59) can be avoided by using an alternate approach based on the asymptotic behavior (Jones and Kline, Reference Jones and Kline1958) of integral $\int \int \exp [ikzf(x,y)]g(x,y)dxdy$ (Bohren and Huffman, Reference Bohren and Huffman1983).

It is apparent that the extinction cross section of a particle may depend on the polarization state of the incident light. This optical phenomenon is analogous to the dichroism observed in conjunction with anisotropic materials (Mishchenko et al., Reference Mishchenko, Travis and Lacis2002; Mishchenko, Reference Mishchenko2008). Shefer and Popov (Reference Shefer and Popov2010) investigated the elements of the extinction matrix based on physical optics and reported that the off-diagonal elements differ from the first element by an order of magnitude.

To quantify the effect of the polarization state on the extinction cross section, we introduce a factor:

$f_{ext}=\frac{C_{ext}(1,Q/I,U/I,V/I)}{C_{ext}(1,0,0,0)}.$(1.2.70)

To illustrate the dependence of the earlier factor on the polarization state of the incident light, let us consider the scattering of polarized light by a hexagonal plate. As shown in Figure 1.14a, OX′Y′Z′ is a particle coordinate system fixed to the particle, and OZ is the z-axis of the incident-light coordinate system, pointing along the incident direction. The orientation of the particle system relative to the incident-light system can be specified in terms of three Euler angles $\alpha$, $\beta$, and $\gamma$. Figure 1.14 illustrates the definition of angles $\beta$ and $\gamma$. Angle $\alpha$ (not shown in Figure 1.14) indicates the rotation of the particle around the Z-axis. The extinction cross section of the particle depends on $\beta$ and $\gamma$ but independent of $\alpha$.

Figure 1.14

(a) Illustration of Euler angles $\beta$ and $\gamma$ in a three-dimensional plot. OX′Y′Z′ is the particle coordinate system. OXYZ is the incident-light coordinate system (only the Z axis is shown). The incident direction is along OZ. $\beta$ is the angle between OZ and OZ′. $\gamma$ is the angle between the projection of OZ on the OX′Y′ plane and OY′. (b) A top view of the particle along OZ′. The dashed line is the projection of OZ on the OX′Y′ plane. The definition of $\gamma$ is shown.

Figures 1.15 and 1.16 show $f_{ext}$ as a function of the orientation of a hexagonal column and plate for the incident light with a Stokes vector of (1,1,0,0) based on the IITM (Johnson, Reference Johnson1988; Bi and Yang, Reference Bi and Yang2014a). The aspect ratio of the hexagonal column is 1.5 and the counterpart of the hexagonal plate is 0.5. The refractive index is 1.3 + i0.0, the size parameter defined in terms of the height of the hexagonal column and the thickness of the hexagonal plate is 30, and angle $\beta$ indicates the tilting of the symmetry axis of the particle from the zenith. As demonstrated in Figures 1.15 and 1.16, $f_{ext}$ can be either larger or less than unity, depending on the particle orientation as well as particle shape (i.e., the particle aspect ratio in the present case).

Figure 1.15

Effect of the polarization state and the orientation of a hexagonal column with a refractive index of m = 1.3 + i0.0 on the extinction of light. Three orientations are considered and specified by $\gamma$ angles equal to 10°, 20°, and 30°.

Figure 1.16

The same as Figure 1.15, but for a hexagonal plate.

In summary, the extinction cross section depends on particle orientation and the polarization state of the incident light. Specifically, the extinction cross section for a linearly polarized radiation beam can be larger or smaller than that for an unpolarized radiation beam, depending on the orientation of the scattering particle with respect to the incident direction. Furthermore, it should be pointed out that the complete extinction matrix rather than the scalar extinction cross section is required for describing the transfer of polarized light (Mishchenko et al., Reference Mishchenko, Travis and Lacis2006) in a scattering medium consisting of oriented particles.

1.2.4 Scattering Phase Matrix

In conjunction with the transformation of the incident electric field to the scattered electric field in Eq. (1.2.57), the corresponding transformation for the Stokes parameters is in the form

$\left[\begin{array}{c}{I}^s\\ Q^s\\ U^s\\ V^s\end{array}\right]=\frac{F}{k^2{r}^2}\left[\begin{array}{c}{I}^i\\ Q^i\\ U^i\\ V^i\end{array}\right],$(1.2.71)

where superscripts s and i indicate the scattered and incident light, respectively. F is a 4 × 4 matrix in the form

$F=\left[\begin{array}{cc}\begin{array}{cc}{F}_{11}& F_{12}\\ F_{21}& F_{22}\end{array}& \begin{array}{cc}{F}_{13}& F_{14}\\ F_{23}& F_{24}\end{array}\\ \begin{array}{cc}{F}_{31}& F_{32}\\ F_{41}& F_{42}\end{array}& \begin{array}{cc}{F}_{33}& F_{34}\\ F_{43}& F_{44}\end{array}\end{array}\right].$((1.2.72))

The F matrix elements are not independent in some scenarios. The relationships of the elements are discussed in many studies (Perrin, Reference Perrin1942; van de Hulst, Reference van de Hulst1957; Takano and Jayaweera, Reference Takano and Jayaweera1985; Hu et al., Reference Hu1987; Hovenier and van de Mee, Reference Hovenier and van der Mee1996; Hovenier and Mackowski, Reference Hovenier and Mackowski1998). Figure 1.17 illustrates the incident-scattering configurations for an arbitrarily shaped particle at its reciprocal and mirror positions and can be used to explain the relationships of the F matrix elements. In the laboratory coordinate system, the incident and scattering directions are fixed. The position of the particle varies in Figure 1.17(b)–(d). In the particle coordinate system, the particle’s position is fixed and the incident and scattering directions vary. In reference to the configuration depicted in Figure 1.17a, the line that is confined on the scattering plane and bisects the angle (π − $\theta$) between the incident and scattering directions is referred to as the bisectrix, and the plane through the bisectrix and normal to the scattering plane is referred to as the bisectrix plane (van de Hulst, Reference van de Hulst1957). Figure 1.17b depicts the reciprocal configurations. In particular, the original particle position is rotated 180° about the bisectrix to result in the reciprocal position. The reader is referred to Hovenier et al. (Reference Hovenier, van der Mee and Domke2004) for a clear explanation of the reciprocity principle in the context of light scattering by an individual particle. The configurations shown in Figure 1.17c and 1.17d correspond to the mirror configurations of those in Figure 1.17a and 1.17b, respectively. Note that the configurations illustrated in Figure 1.17(a) and the left panels of Figure 1.17(b)–(d) are, in principle, the same as van de Hulst’s illustrations presented with respect to the incident-light system, while the illustrations in the right panels of Figure 1.17(b)–(d) are presented with respect to the particle coordinate system in the same form as Sun et al. (Reference Sun2020).

Figure 1.17

Illustrations of the incident and scattering directions in laboratory and particle coordinate systems. The illustration on the left of each panel is in principle the same as that in van de Hulst (Reference van de Hulst1957) and shows the incident $(e{̂}_i)$ and scattering $(e{̂}_s)$ directions and associated polarization directions, in reference to the laboratory coordinate system. The illustration on the right of each panel is in reference to the particle coordinate system. The scattering plane formed by the incident and scattering directions is coplanar with the paper. (a) An arbitrarily shaped particle in its original position. The incident-scattering configurations are identical in laboratory and particle coordinate systems. $\theta$ is the scattering angle. A bisectrix bisects the angle π − $\theta$; (b) the particle is at the reciprocal position (van de Hulst, Reference van de Hulst1957) with respect to (a), and is obtained by rotating it 180° about the bisectrix; (c) the particle is at the mirror position with respect to (a); (d) the particle is at the mirror position with respect to (b). In (c) and (d), the mirror is the scattering plane.

When the amplitude scattering matrix and the phase matrix corresponding to configuration 1.17a are given by Eqs. (1.2.58) and (1.2.72), respectively, the amplitude scattering and phase matrices corresponding to Figures 1.17b, 1.17c, and 1.17d are given by

$S_b=\left[\begin{array}{cc}{S}_2& -S_4\\ -S_3& S_1\end{array}\right],$(1.2.73)

$F_b=\left[\begin{array}{cccc}{F}_{11}& F_{21}& -F_{31}& F_{41}\\ F_{12}& F_{22}& -F_{32}& F_{42}\\ -F_{13}& -F_{23}& F_{33}& -F_{43}\\ F_{14}& F_{24}& -F_{34}& F_{44}\end{array}\right].$(1.2.74)

$S_c=\left[\begin{array}{cc}{S}_2& -S_3\\ -S_4& S_1\end{array}\right],$(1.2.75)

$F_c=\left[\begin{array}{cccc}{F}_{11}& F_{12}& -F_{13}& -F_{14}\\ F_{21}& F_{22}& -F_{23}& -F_{24}\\ -F_{31}& -F_{32}& F_{33}& F_{34}\\ -F_{41}& -F_{42}& F_{43}& F_{44}\end{array}\right].$(1.2.76)

$S_d=\left[\begin{array}{cc}{S}_2& S_4\\ S_3& S_1\end{array}\right],$(1.2.77)

$F_d=\left[\begin{array}{cccc}{F}_{11}& F_{21}& F_{31}& -F_{41}\\ F_{12}& F_{22}& F_{32}& -F_{42}\\ F_{13}& F_{23}& F_{33}& -F_{43}\\ -F_{14}& -F_{24}& -F_{34}& F_{44}\end{array}\right].$((1.2.78))

For many applications, such as in atmospheric remote sensing implementations and plane-parallel radiative transfer simulations, we need to consider a scattering particle system that can be reasonably assumed to be randomly oriented with an equal number of mirror-image positions with respect to both the scattering plane and the reciprocal position plane. This of course implies that there are equal numbers of particles in positions (a)–(d) in Figure 1.17. In this text, we assume this random orientation condition for the particles unless specifically stated. Furthermore, we assume that the scattering by an ensemble of particles or a particle system is incoherent, implying that the scattering phase matrix of the particle system is given by the sum of those of the individual particles in the particle ensemble. Under such assumptions, the mean matrix F of the particle system can be simplified based on Eqs. (1.2.72) and (1.2.76) as

$F=\left[\begin{array}{cccc}{F}_{11}& F_{12}& 0& 0\\ F_{12}& F_{22}& 0& 0\\ 0& 0& F_{33}& -F_{43}\\ 0& 0& F_{43}& F_{44}\end{array}\right]$

$=\left[\begin{array}{cccc}(M_2+M_3+M_4+M_1)/2& (M_2-M_1)/2& 0& 0\\ (M_2-M_1)/2& (M_2-M_3-M_4+M_1)/2& 0& 0\\ 0& 0& S_{12}+S_{34}& -D_{21}\\ 0& 0& D_{21}& S_{12}-S_{34}\end{array}\right],$(1.2.79)

where

$M_j=|S_j{|}^2,$(1.2.80a)

$S_{kj}=S_{jk}=\left(S_j{S}_k^{*}+S_k{S}_j^{*}\right)/2,$(1.2.80b)

$-D_{kj}=D_{jk}=\left(S_j{S}_k^{\text{*}}-S_k{S}_j^{\text{*}}\right)i/2,j,k=1,2,3,4.$(1.2.80c)

Takano and Jayaweera (Reference Takano and Jayaweera1985), Hu et al. (Reference Hu1987), and Hovenier and Mackowski (Reference Hovenier and Mackowski1998) discussed the symmetry relations in the forward and backward directions for light scattering by nonspherical particles. Special attention is called for a particle with a plane of symmetry under the random orientation condition. In this case, the nonzero elements of the F matrix in the forward direction reduce to

$F_{11}=(M_2+M_3+M_4+M_1)/2,$(1.2.81a)

$F_{22}=F_{33}=(M_2-M_3-M_4+M_1)/4+(S_{12}+S_{34})/2,$(1.2.81b)

$F_{44}=S_{12}-S_{34}.$(1.2.81c)

If the particle is rotationally symmetric, the following relation holds in the 0° forward scattering direction (Hovenier and Mackowski, Reference Hovenier and Mackowski1998):

$F_{44}=2F_{22}-F_{11}.$(1.2.82)

For a particle under the random orientation condition, the nonzero elements of the F matrix in the backward direction reduce to (Takano and Jayaweera, Reference Takano and Jayaweera1985; Hovenier and Mackowski, Reference Hovenier and Mackowski1998)

$F_{11}=(M_2+2M_3+M_1)/2,$(1.2.83a)

$F_{22}=-F_{33}=(M_2+M_1)/4-S_{12}/2,$(1.2.83b)

$F_{44}=M_3+S_{12}.$(1.2.83c)

If the particle has a plane of symmetry, the following relation holds in the 180° backward scattering direction:

$F_{44}=F_{11}-F_{22}.$(1.2.84)

In practice, the relations in Eqs. (1.2.72)–(1.2.77) can largely simplify the numerical computation or they can be used to validate numerical results. Furthermore, the F matrix in Eq. (1.2.79) can be expressed in a normalized form as follows:

$F=\left[\begin{array}{cccc}{F}_{11}& F_{12}& 0& 0\\ F_{12}& F_{22}& 0& 0\\ 0& 0& F_{33}& -F_{43}\\ 0& 0& F_{43}& F_{44}\end{array}\right]=\frac{C_{sca}}{4\pi }\left[\begin{array}{cccc}{P}_{11}& P_{12}& 0& 0\\ P_{12}& P_{22}& 0& 0\\ 0& 0& P_{33}& -P_{43}\\ 0& 0& P_{43}& P_{44}\end{array}\right],$(1.2.85)

where $C_{sca}$ is the scattering cross section given by

$C_{sca}={\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{\pi }{F}_{11}(\theta )\sin \theta d\theta d\phi }}.$(1.2.86)

It can be shown that $P_{11}(\theta )$ (usually, called the phase function) is normalized in the form:

$\frac{1}{4\pi }{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{\pi }{P}_{11}(\theta )\sin \theta d\theta d\phi =1}}.$(1.2.87)

The matrix P in Eq. (1.2.85) is referred to as the normalized scattering phase matrix. In the case of the scattering of light by a sphere, we have $P_{43}=0$, $P_{33}=P_{44}$, and $P_{11}=P_{22}$. Since $P_{22}/P_{11}=1$ for a sphere, the deviation of the ratio $P_{22}/P_{11}$ from unity is an indicator of the degree of nonsphericity (Bohren and Huffman, Reference Bohren and Huffman1983; Mishchenko et al., Reference Mishchenko, Travis and Lacis2006).

For practical applications, the scattering phase matrix elements can be expanded in terms of the generalized spherical function. In this way, a large look-up table for tabulated phase matrix element data can be avoided. In particular, following Mishchenko et al. (Reference Mishchenko, Travis and Lacis2006), we have

$P_{11}\left(\text{cos}\theta \right)=\sum _{l=0}^L{C}_{11,l}{p}_{00}^l\left(\text{cos}\theta \right),$(1.2.88a)

$P_{12}\left(\text{cos}\theta \right)=\sum _{l=2}^L{C}_{12,l}{p}_{02}^l\left(\text{cos}\theta \right),$(1.2.88b)

$P_{22}\left(\text{cos}\theta \right)=\frac{1}{2}\sum _{l=2}^L\left[\left(C_{22,l}+C_{33,l}\right)p_{22}^l\left(\text{cos}\theta \right)+\left(C_{22,l}-C_{33,l}\right)p_{2,-2}^l\left(\text{cos}\theta \right)\right],$(1.2.88c)

$P_{33}\left(\text{cos}\theta \right)=\frac{1}{2}\sum _{l=2}^L\left[\left(C_{22,l}+C_{33,l}\right)p_{22}^l\left(\text{cos}\theta \right)-\left(C_{22,l}-C_{33,l}\right)p_{2,-2}^l\left(\text{cos}\theta \right)\right],$(1.2.88d)

$P_{34}\left(\text{cos}\theta \right)=\sum _{l=2}^L{C}_{34,l}{p}_{02}^l\left(\text{cos}\theta \right),$(1.2.88e)

$P_{44}\left(\text{cos}\theta \right)=\sum _{l=0}^L{C}_{44,l}{p}_{00}^l\left(\text{cos}\theta \right),$(1.2.88f)

where $C_{11,l}$, $C_{12,l}$, $C_{22,l}$, $C_{33,l}$, $C_{34,l}$, and $C_{44,l}$ are the expansion coefficients. $p_{mn}^l$ are the generalized spherical functions that are complex functions given by the Wigner d-functions. The reader is referred to Appendix F of Mishchenko et al. (Reference Mishchenko, Travis and Lacis2006) for detailed discussions of the generalized spherical functions and the Wigner d-functions.

In many practical applications where computational feasibility and efficiency are critical, we do not need the detailed angular variation of the phase function, such as in the radiative transfer models involved in climate models due to consideration of computational feasibility and efficiency. Instead, we introduce a parameter called the asymmetry factor g to quantify the overall degree of the deviation of the scattered light from the incident direction, which is given by

$g=〈\cos \theta 〉=\frac{1}{4\pi }{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{\pi }{P}_{11}\left(\theta ,\phi \right)\cos \theta \sin \theta d\theta d\phi }}.$(1.2.89a)

If the phase function is independent of the azimuthal angle of the scattering plane, Eq. (1.2.89a) is simplified to the following form:

$g=〈\cos \theta 〉=\frac{1}{2}{\displaystyle {\int }_0^{\pi }{P}_{11}\left(\theta \right)\cos \theta \sin \theta d\theta .}$(1.2.89b)

An analytical phase function specified in terms of the asymmetry factor, known as the Henyey–Greenstein (H-G) phase function (Henyey and Greenstein, Reference Henyey and Greenstein1941), has been widely used in radiative transfer simulation, particularly for irradiance computation, and is given by

$P(\theta )=\frac{1-g^2}{{(1+g^2-2g\cos \theta )}^{3/2}}.$(1.2.90)

It is worth noting that the (H-G) phase function can be analytically expanded in terms of the Legendre polynomials in the form:

$P(\cos \theta )={\displaystyle \sum _{n=1}^{\infty }(2n+1)g^n{p}_n(\cos }\theta ),$(1.2.91)

where $P_n(\cos \theta )$ is the nth-order Legendre polynomials. The feature of the H-G phase function shown in Eq. (1.2.91) is convenient for the truncation of the phase function in multiple scattering (i.e., radiative transfer) simulation (e.g., Wiscombe, Reference Wiscombe1977), particularly for the two-stream radiative transfer scheme (Meador and Weaver, Reference Meador and Weaver1980) that has been extensively used for radiative flux calculations in climate models.

However, the H-G phase function can be quite different from a phase function with the same asymmetry factor. For example, Figure 1.18 shows the phase functions of a water cloud and an ice cloud in comparison with their corresponding Henyey–Greenstein phase function approximations. The water and ice cloud phase functions in Figure 1.18 are taken from those used in the operational MODIS cloud property retrievals (Platnick et al., Reference Platnick2017). For remote sensing applications, the detailed phase function is necessary for radiance simulation involved in generating look-up tables for implementing a retrieval algorithm. In this case, the H-G phase function is usually not applicable.

Figure 1.18

Comparison between rigorously computed liquid cloud and ice cloud phase functions and the corresponding Henyey–Greenstein (H-G) phase functions based on the same values of the asymmetry factor. In each plot, the H-G phase function has the same asymmetry factor as the cloud phase function.

(a)

Liquid cloud with effective radius 20 μm, and the asymmetry factor (g) is 0.871;

(b)

Ice cloud with effective radius 30 μm, and g is 0.756. The wavelength is 0.87 μm. Note that the effective radius is defined as $R_e=\frac{3}{4}\frac{\int V(D)N(D)dD}{\int A(D)N(D)dD}$ following Foot (Reference Foot1988) as a generalization of the definition of the effective radius for a polydisperse ensemble of spherical particles (Hansen and Travis, Reference Hansen and Travis1974), where D, V, and A indicate the maximum dimension, volume, and projected area of an ice crystal, respectively. N(D) indicates the particle size distribution.