Stellar Aureole Measurements

Retrieving Cirrus Microphysical Properties from Stellar Aureoles

John G. DeVore11affiliation: Visidyne, Inc., Santa Barbara, California, USA; , Joseph A. Kristl22affiliation: Visidyne, Inc., Burlington, Massachusetts, USA , Saul A. Rappaport33affiliation: 37-647 M.I.T. Department of Physics, 70 Vassar St., Cambridge, MA, 02139
Journal of Geophysical Research: Atmospheres, in press

The aureoles around stars caused by thin cirrus limit nighttime measurement opportunities for ground-based astronomy but can provide information on high-altitude ice crystals for climate research. In this paper we attempt to demonstrate quantitatively how this works. Aureole profiles can be followed out to from stars and from Jupiter. Interpretation of diffracted starlight is similar to that for sunlight, but emphasizes larger particles. Stellar diffraction profiles are very distinctive, typically being approximately flat out to a critical angle followed by gradually steepening power-law falloff with slope less steep than . Using the relationship between the phase function for diffraction and the average Fourier transform of the projected area of complex ice crystals we show that defining particle size in terms of average projected area normal to the propagation direction of the starlight leads to a simple, analytic approximation representing large-particle diffraction that is nearly independent of crystal habit. A similar analytic approximation for the diffraction aureole allows it to be separated from the point spread function and the sky background. Multiple scattering is deconvolved using the Hankel transform leading to the diffraction phase function. Application of constrained numerical inversion to the phase function then yields a solution for the particle size distribution in the range between and . Stellar aureole measurements can provide one of the very few, as well as least expensive, methods for retrieving cirrus microphysical properties from ground-based observations.

slugcomment: blank space

1. Introduction

Accurate retrieval of the properties of particles in the atmosphere, e.g., ice clouds, is of prime importance to understanding their role in climate change and modeling their effects in global simulations (Comstock et al., 2007). Knowledge of the impact of aerosols on climate change has improved so much that in their contribution to the “Fourth Assessment Report of the Intergovernmental Panel on Climate Change” Forster et al. (2007) called AERONET (Holben et al. 1998) a “significant advancement”. However, the impact of cirrus cloud particles is much less certain because they occur high in the atmosphere and are more difficult to monitor.

DeVore et al. (2009) described a new Sun and Aureole Measurement (SAM) instrument that measures the radiance profile of the solar disk and surrounding aureole (from to ) in a narrow spectral band centered on a wavelength of 0.67 m. Applying the ‘diffraction approximation’ (presented in their paper) to SAM measurements they were able to retrieve size distributions of cloud particles in the range from 5 to 50 m. In practice, SAM has been most useful for collecting information on cirrus clouds because of their thinness (i.e., optical depth to ), which is required for the appearance of aureoles.

In this paper we report what we believe to be a new method for making and interpreting aureole measurements at night. Although the Moon provides an alternative to the Sun, we think that stars provide the best sources and that planets may have some utility as well. The angular size of the Moon is basically the same as that of the Sun and therefore does not allow for retrieving information on large particles that comes from measurements at small angles. Moreover, both the Moon and planets limit where in the sky measurements can be made, whereas suitable stars are distributed across the entire sky. We have used inexpensive technology to measure starlight diffraction in the angular range from to , which corresponds to the size range from to . Since the climate impact of cirrus is sensitive to its microphysical properties through both particle scattering asymmetry and emissivity/absorptivity (Stephens et al. 1990), it is important to be able to measure the size distibutions of large ice crystals such as those we discuss in this work.

Our approach is based on aureole measurements with a lens and a medium-quality astronomical CCD camera. This remote sensing method has advantages over in-situ measurements in that (i) it can readily be carried out on virtually any night when thin cirrus clouds are visible, (ii) it is relatively inexpensive to implement, and (iii) the measurements do not disturb the cloud environment. The aureole approach to monitoring cirrus provides one of the few methods available to retrieve cirrus particle microphysical properties from ground-based observations.

§2 describes the various measurements involved in interpreting stellar aureoles, starting with the determination of two important parameters used in aureole profile interpretation, optical depth (§2.1) and stellar exo-atmospheric irradiance (§2.2), followed by a discussion of aureole imagery (§2.3). Examples of star images illustrate the difference between the weak aureoles associated with small particles, e.g. thin water clouds, and the strong ones attributable to the large ice crystals in cirrus.

§3 presents various aspects of the theory and modeling of diffraction by ice crystals used to interpret the measurements. The reader is reminded that despite the wide range of crystal habits that need to be considered, the phase function representing diffraction can be calculated simply using the Fourier transform. Moreover, this solution illuminates a simple wavelength scaling relationship for aureole profiles. Examination of phase functions derived for different habits with the same size measure leads to the selection of a size metric based on average projected area as being particularly useful. Phase functions representing diffraction from distributions of particle sizes are discussed next, followed by discussion of the potentially important role of multiple scattering and an analytic approximation to represent aureole profiles involving multiple scattering. Deconvolution of multiple scattering is presented in the appendix. Next the steps involved in retrieving the particle size distribution from the phase function are presented (§3.4), either by fitting the data with an analytic form or using constrained numerical inversion. Finally, the effects of the point spread function are shown to be taken into account using an appropriate analytic approximation (§3.5).

§4 illustrates the application of theory and modeling to interpret aureole profiles using an example. The measured aureole profile is separated into three components so that the diffraction profile can be distinguished from the point spread function and the sky background (§4.1). Next the phase function is deconvoled from the diffraction profile (§4.2) and then the particle size distribution (PSD) is determined (§4.3) from it. The section concludes with some more examples of stellar aureoles covering a range of optical depths (§4.4) and a brief look at an example of an aureole around the planet Jupiter.

§5 provides a short summary of our work. This is followed by an appendix (§7) discussing the deconvolution of multiple scattering, where it is approximated as a Poisson process (§7.1). The resulting series is solved for single scattering, double scattering, and then for all orders (§7.4).

2. Stellar Aureole Measurement

From our Visidyne measurement site in West Acton, MA (Fig. 1), we collected an extensive stellar aureole dataset comprising approximately 42 hours of imagery from 17 nights between 19 Jun 2011 and 26 Feb 2012 in a narrow 5 nm band around 672 nm. We used multiple exposures to extend the dynamic range. Targets ranged from bright planets such as Jupiter down to magnitude 3 stars. Sky conditions varied from clear to sufficiently cloudy that stars were practically indistinguishable from their aureoles, analogous to the situation with a “fuzzy Sun” (Linskens and Bohren, 1994). The clouds consisted primarily of cirrus although water clouds were present at times. The two types of clouds are readily distinguishable from their visual appearance on all-sky imagery as well as from the strength of their aureoles. As with the SAM instrument for solar aureoles, two cameras were used, one to measure particulate optical depth while the second measures the aureole radiance profile at the same time. Supporting data include all-sky imagery, laser cloud altimetry, and both local and synoptic weather data. Fig.  2 shows a typical image from the all-sky camera with targets identified through the cirrus cloud.

Figure 1.— West Acton Observatory where the stellar aureole data were collected for Phase I.
Figure 2.— Example all-sky image from the evening of 16 Dec 2011 showing the locations of the star Capella and planet Jupiter.

2.1. Optical Depth

The central star pixels in the stellar aureole images (discussed in §2.3) were saturated in order to set the aureole camera exposure times to yield the best signal for the relatively faint aureole. To measure optical depth at the aureole star, a second, co-aligned camera collected short-exposure, unsaturated images of the aureole star. A filter was used to allow short exposures (0.05 to 1 second), even during high optical depth conditions. The line-of-sight transmittance was measured by referencing the measured flux to a clear-sky measurement of the same star. The data were scaled by the frame exposure time, and then the optical depth along the line of sight, , was calculated directly as follows:


where is the measured irradiance and 2.2) is a reference measurement on a clear night. Optical depth was measured once per second, and up to 5 depth values were averaged to reduce scatter due to atmospheric scintillation. Determination of optical depth using Sun or star photometry needs to take into account forward scattering. Since the correction is a strong function of the angular width of the source irradiance measurement (e.g., Shiobara and Asano,1994; DeVore et al, 2009), the smaller angular width of stars as compared with the Sun significantly reduces the size of the corrections required. Also, to the extent that the data analysis can distinguish between direct and scattered source radiance, the error attributable to forward scattering can be reduced further.

Some typical optical depth data are shown in Fig. 3. The line-of-sight optical depth is plotted in blue, while the start of each 30-second aureole measurement is indicated with a red diamond. In general, the cirrus clouds are not spatially uniform so that rapid changes in occur as the clouds drift through the line-of-sight. The high sample rate is used to identify and then discard aureole measurements that contain significant variation (as indicated by the fluctuations in the blue curve) during the longer aureole integration times (as indicated by the separations between the red diamonds). The primary source of error is atmospheric scintillation, which is seen in this type of optical depth measurement even on clear nights. Analysis of clear night data over 30-second time intervals (the duration of aureole image exposures) returns RMS errors of 0.1 to 0.25 in line-of sight optical depth. We recommend use of the conservative error value for data presented here. Temporal variations caused by cloud motion also limit the total exposure time practical for any single star aureole measurement.

Figure 3.— Illustrative measurement data for the star Capella on 17 Dec 2011. The red diamonds indicate the start of 30-second aureole radiance measurements.

2.2. Stellar Exo-Atmospheric Irradiance

As with the determination of , absolute irradiance values are not necessary for specifying . It suffices to define irradiance for a specific camera in terms of the number of photons counted divided by the product of the aperture area of the lens and the exposure time , giving units of . For radiance one divides the irradiance by the pixel field of view , giving units of .

We examined 16 unsaturated images containing Capella from 26 Nov 2011, when the sky was very clear. For each image we determined the maximum, , and average, , number of counts and the corresponding standard deviation, , in the pixel square surrounding Capella. This size square was large enough so that the minimum number of counts at the edges was down by 2 orders of magnitude from the peak. In order to identify those counts associated with the directly transmitted stellar photons, we separated the individual counts into “peak” and “background” bins using the following process. First, we selected a somewhat arbitrary threshold value, , to perform an initial separation. Next we calculated the average, , and standard deviation, , of the counts in the background bin. Guided by the three-sigma rule (e.g., Kreyszig, 1979), we reset and redivided the counts into the two bins. Then we calculated the sum of the counts, , and the number of pixels, , in the peak bin, and recalculated the average number of counts, , in the background bin. Then we calculated as:


For the 16 images analyzed, the average with a standard deviation of . We selected the maximum of the set, , to use under the assumption that it represents the clearest sky conditions.

2.3. Aureole Profiles

The stellar aureole profiles were collected with a QSI-583wg astronomical camera, interfaced to a Canon 70-200 mm f/# 2.8 lens set to 200 mm focal length. The camera and lens were mounted to a German equatorial mount that located stars and corrected for sidereal motion. The filter was an Astrodon 672 nm filter with 5 nm bandwidth, chosen to minimize light pollution and provide a narrow band for model calculations. The camera was operated for most collects using pixel binning to enhance the relatively faint aureoles, resulting in single pixel angular sizes of approximately . In later collects the binning was reset to , providing higher angular resolution (). Each aureole image was a 30 second camera exposure. After the collection, frames were grouped by optical depth (as measured by the optical depth camera) and stacked to improve signal-to-noise. Images that contained significant cloud spatial gradients were discarded prior to stacking.

Fig. 4 shows a montage of images of aureoles around the star Capella taken through a thin cirrus layer that varied over the course of the measurements. The images are annotated with measured by a second camera. Note the expected correlation between the apparent size and intensity of the aureole and . The dashed yellow circle shows the size of the Moon for reference. The overall brightening of the background represents a combination of the “wings” of the aureoles, scattering of cityshine from the greater Boston area, and detector noise. Fig. 5 compares the images of the star Capella as seen through cirrus, clear sky, and a water cloud on a different night. Note that despite the fact that in this case is close to the value that should maximize aureole radiance (see below), the water cloud does not produce a discernable aureole. This is due to the fact that small droplets of water diffract light out to angles larger than our field of regard. (Since the color scalings applied in the two sets of figures differ slightly, one should not place any particular interpretation on the light blue versus black background shading.)

Figure 4.— Images of aureoles around the star Capella in a narrow 5 nm band around 672 nm during the evening of 16 Dec 2011 annotated with the cirrus optical depth along the line of sight, , as measured by a second camera.
Figure 5.— Images of the star Capella in a narrow 5 nm band around 672 nm through different cloud conditions during the evening of 11 Jan 2012 annotated with the particulate optical depth along the line of sight, , as measured by a second camera.

Fig. 6 shows more quantitatively some illustrative aureole profiles about the star Capella as a function of angle. The aureoles span a range of from 0.1 to 2.3. The colored symbols indicate measured radiances and the solid curves show fits to analytic functions representing the three physical phenomena, the point spread function, stellar diffraction, and sky background, responsible for the shapes of aureole profiles (see § 4 for details). The point spread function results from instrumental scattering and atmospheric turbulence and dominates the aureole profiles at small angles; the resulting radiance is anticorrelated with . The background primarily represents diffusely scattered cityshine and dominates the aureole profiles at large angles; the resulting radiance is positively correlated with . Stellar diffraction gives the aureole its shape at intermediate angles and is the main subject of this paper; the radiance from diffraction peaks for between about 1 and 2.

Figure 6.— Examples of aureole profiles and physical model fits for the star Capella from the night of 27 Nov 2011, annotated by (“OD” in the figure legend).

3. Aureole Theory and Modeling

Various aspects of the theory and modeling of aureoles are covered in this section to provide a basis for their application in the next section. Despite the wide range of crystal habits that need to be considered, the phase function describing diffraction from individual ice crystals can be calculated simply using the Fourier transform (§3.1). The solution exhibits a potentially useful wavelength scaling relationship. Airy’s analytic solution for a sphere is presented both to provide validation for the implementation of the numerical solution and as a guide for a simpler analytic approximation. Examination of phase functions derived for different habits with the same size measure leads to the selection of a size metric based on the average projected area as being particularly useful. Distributions of particle sizes are discussed next (§3.2) considering both power-law and exponential analytic forms, which leads to a useful analytic approximation for the phase function for a distribution of particles. Next the diffraction radiance profile is considered (§3.3) as well as the potentially important role of multiple scattering. Example calculations suggest a useful analytic approximation to represent multiply-scattered diffraction profiles, which can be deconvolved using the equations presented in the appendix. Next the steps involved in retrieving the particle size distribution from the phase function are presented (§3.4), either by fitting the data with an analytic form or using constrained numerical inversion. Finally the effects of the point spread function are taken into account using an appropriate analytic function (§3.5).

3.1. Diffraction From Individual Ice Crystals

Diffraction is the dominant mechanism involved in forming aureoles. The angles involved in stellar aureole measurements tend to be smaller than the size of the Sun, and hence smaller than the angles measured in solar aureoles. This difference in angular measurement range has important implications for modeling the particles most affecting the aureole shape. Whereas solar aureole measurements range between and , and correspond to particle dimensions between and according to the diffraction approximation (DeVore et al., 2009), stellar aureole measurements range between and and correspond to particle dimensions an order of magnitude larger. Ice crystals in the sensitive size range for solar aureoles tend to be compact, e.g., droxtals, and hence are well approximated as spheres. However, the ice crystals relevant to stellar aureole modeling require consideration of more complex shapes (e.g., Baum et al. 2005).

3.1.1 Diffraction Calculation

Consider diffraction from the 8 crystal habits shown in Fig. 7. Triangular facets are used to describe the outer surfaces of each crystal. Although not representative of a naturally occurring ice crystal we include the spherical shape both as an approximation for small ice crystals and to check our numerical calculations of the diffraction patterns against the well-known solution of Airy (1835) for a sphere. We approximate a sphere (panel a) using a surface with 512 triangular facets, generated by starting with an octahedron and sequentially subdividing each triangle into 3 smaller triangles 3 times, each time moving the vertices to the surface of the circumscribing sphere. Panel (b) shows a droxtal with maximum volume relative to the circumscribing sphere used to represent small, compact ice crystals (Yang et al., 2003). While columns, plates, and bullets are observed, aggregations of bullets are more commonly used to represent the larger ice crystals in cirrus. Panels (c) and (d) show a solid hexagonal column and plate with ratios of the width to length or height based on the scaling relations for crystal habits C1e and P1a provided by Heymsfield and Platt (1984). The bullet and bullet rosettes in panels (e) to (h) use the tip angle and width-to-length scaling relations given by Um and McFarquhar (2007).

Figure 7.— Triangular facetized surfaces used to represent 8 different crystal habits. The sphere is included for comparison with the Airy function solution for diffraction from a circular aperture.

Bi et al. (2011) discuss several ways of calculating the scattering of light by particles larger than the wavelength of the incident radiation. Born and Wolf (1959) relate Fraunhofer diffraction to the 2-dimensional Fourier transform of the projection of the particle on the plane normal to the incident direction of the radiation. Although the line integral formulation of Gordon (1975) may be slightly faster, we have found numerical calculation using the fast Fourier transform (Press et al., 1992) adequate and we can calculate the projected area, , which we have found useful, in the process. Using the small angle approximation (Lenoble 1985), a phase function representing the angular distribution of diffracted radiation can be written as:


where is the scattering angle, is an azimuthal angle, and are cartesian coordinates in the projection plane with their origin within the particle projection, is the 2-dimensional Fourier transform operator, and is the “aperture” function describing the projection of the particle:


We represent the aperture function by a two-dimensional, uniformly spaced grid of points. We set the grid dimensions to be 1024 1024 in the examples presented here unless specified otherwise. We set the physical size of the grid in each dimension to the product of a characteristic crystal dimension (e.g., the diameter of an hexagonal plate or the length of a solid column) and the square root of the number of grid points. We project the crystal facets onto a gridded plane oriented normal to the direction of the incident radiation. Each grid point is initialized to zero and set to one if any projected facet covers it.

Unlike an ordinary phase function, is defined only in the forward hemisphere and is meant to represent the diffraction component of the angular distribution of radiation in the near-forward direction. For cases of interest here, i.e., ice crystals large compared with and , diffraction dominates the internal or body scattering component, which can be ignored. However, in order to relate correctly to extinction and optical depth , both the diffraction and internal components are included in the extinction cross section (from the optical theorem, (e.g., Liou, 2002).

3.1.2 Wavelength Scaling

From Eqn. 3 we see that the amplitude of diffraction phase functions scales inversely with and that the scattering angle scales with . These observations allow for scaling diffraction phase functions in the near forward direction from one wavelength to a second wavelength :


for and . The interpretation of Eqn. 5 is simple. For a given size crystal as the wavelength increases (decreases) the projected area of the crystal becomes smaller (larger) relative to the wavelength and the scattered radiation spreads to larger (smaller) angles. However, since the total diffracted power remains the same, the radiance at any angle decreases (increases) by the square of the wavelength to conserve energy. Although we shall retain in a number of equations to emphasize its role, we shall take in the calculations and examples in this work.

3.1.3 Airy Solution for a Spherical Particle

It is useful to consider the analytic solution for a sphere (Airy 1835). Expressed as a diffraction phase function, Airy’s solution for a sphere, , is given by:


where is a non-dimensional scaling parameter, is the diameter of the sphere, is the Bessel function of the first kind of order 1, and the small angle approximation () has been applied to find the expression on the right.

Consider the case of light of wavelength incident on a spherical particle [Fig. 7 (a)] with diameter . The dark blue line in Fig. 8 shows the Airy solution (Eqn. 6) for the diffraction phase function. We approximated the sphere using a facetted surface as described in §3.1.1. However, in this case we subdivided the initial octahedon surface 4 times so that the resulting spherical surface had 2048 facets. Each facet was projected onto a plane of points spaced apart. The 1,596 points determined to be inside of projected facets were used to generate an aperture function. They represent a projected area of , which is approximately 1% less than the projected cross section of a diameter true sphere. The red line in the figure shows the phase function generated from application of the Fast Fourier transform (Press et al. 1992) to the aperture function using Eqn. 3. The two curves are remarkably close but can be distinguished because the minima of the Airy solution are somewhat deeper than those of the numerical one. Although we expect that the limited domain of the numerical Fourier transform is responsible for these differences, they are small compared with the approximations introduced in the following sections.

Figure 8.— Diffraction phase functions calculated for a sphere or facetized sphere of diameter and . Also shown are two analytic approximations.

3.1.4 Analytic Approximation

It will prove useful to introduce a simple, analytic approximation, , to the Airy solution for a sphere in the small angle limit by considering its asymptotic behavior. For small angles, , and for large angles, . A simple analytic function with these asymptotic behaviors is:


where parameter is found from the normalization of using the small angle approximation and taking the upper limit of the integral to infinity:


The green line in Fig. 8 shows this analytic approximation. By design, the approximations’s plateau at small angles agrees with the Airy and Fourier transform numerical solutions, while at the larger angles the approximation cuts through their oscillations near the peaks.

3.1.5 Area Diameter

Fig. 9 (a) compares calculations of the diffraction phase functions for the 8 crystal habits shown in Fig. 7, all with the same maximum size, and averaged over particle orientation. The diffraction phase functions are similar, with plateaus at small scattering angles followed by power-law decreases with slopes of . The oscillations in the power-law regions tend to average out with distributions of particle sizes and/or shapes. However, the overall amplitudes for the same maximum size vary by a factor of depending upon particle shape. We found that diffraction phase functions for particles with the same volumes exhibited a factor of 3 variability versus crystal habit. For particles with the same ratios of their volumes to their average projected areas, the variability is over an order of magnitude. In other words, none of these measures of particle size proves self-consistent or useful for aureole work. However, the diffraction patterns for particles with the same average projected areas are very similar as can be seen in Fig. 9 (b), prompting us to define the “area diameter” as the diameter of the circle with the same area as that of the average over orientations of the particle projections. The area diameter therefore appears to be a useful measure of ice crystal size. This should not be surprising since diffraction by a particle is closely related to its projected area.

Figure 9.— Diffraction phase functions calculated for the 8 crystal habits shown in Fig. 7 with (a) maximum dimension = and (b) area diameter = . The latter also shows the analytic approximation (Eqn. 10) for comparison.

Given that phase functions are similar for different crystal habits with the same , including spherical particles, and that the phase functions for the latter can be approximated with a simple analytic function, we are led to replace in the approximation with . In terms of the scattering angle and diameter , the analytic approximation in Eqn. 7 becomes:


The black line in Fig. 9 (b) illustrates how well this approximation fits the diffraction calculations. Therefore, for the remainder of this work we adopt the form of Eqn. 10 to represent a “universal” phase function for any shape particle with area-diameter . The extinction cross section is simply twice the projected area:


3.2. Particle Size Distributions

The simple, analytic approximation is useful for individual ice crystal habits assuming that they are part of a distribution so that the fine structures of the individual phase functions average out. Next we consider different size distributions of ice crystals and look for ones which are both realistic and have convenient parameterizations.

Using the total diffraction phase function representing the average over all sizes and habits is given by the following integral over the particle size distribution, PSD:


where the extinction (rather than the scattering) cross section has been used (since absorption is negligible at visible wavelengths), is the differential number density of particles along the path through the atmosphere per unit area diameter, and is related to as follows:


3.2.1 Power-Law Distribution

Consider a PSD with the power-law form (e.g., Heymsfield and Platt, 1984):


where is a normalization constant, is the power-law slope, and the PSD is truncated below and above . Substitute Eqn. 14 into Eqn. 13 and use Eqn. 11 to find:


when . Fig. 10 (a) shows numerical integrations of Eqn. 12 for with extending from to using . The asymptotic behaviors of , constant for small and proportional to for large , are similar to those of the individual analytic phase function except that for the steeper PSD slopes there is a transition region with an intermediate power-law slope between and .

Figure 10.— (a) Total diffraction phase functions calculated for a range of from to and (b) comparisons with least-squares fits using the analytic expression in Eqn. 27.

Although an analytic expression can be found for , it involves hypergeometric functions and is not very insightful. However, it is illuminating to apply a different approximation to the single particle phase function by replacing with what we call the broken-line power-law approximation :




where is determined from the normalization requirement for to be:


Fig. 8 compares this approximation (the purple curve) with the continuous version (the green curve) from Eqn. 10.

With the limits of the integral in Eqn. 12 again set to and the analytic result for the total phase function naturally divides into three regions depending upon as follows:






where, importantly, . Substituting this expression into Eqn. 21 gives an equation with the following form for the transition region:


where , , and are constants that depend upon , , and . Eqns. 19 to 22 show a total phase function which is constant for small and proportional to for large . If the intermediate transition region of the phase function is modeled with a power-law form, i.e., , then . This result is consistent with the relationship between the phase function (aureole) and PSD power-law slopes found using the diffraction approximation (DeVore et al, 2009).

3.2.2 Exponential Distribution

For later use we note that an exponential form is also sometimes used to model ice crystal PSDs (e.g., Field and Heymsfield, 2003):


where is a normalization constant, is a characteristic size, and the PSD is assumed to be truncated below and above . Substituting Eqn. 23 into Eqn. 13 and using Eqn. 11 gives:






Analytic forms for are even more complicated than those using the power-law form and are not investigated in this work. However, it is useful to note that since the exponential form tends to fall off very rapidly for , frequently little is lost by taking . As a result, the number of free parameters tends to be 1 fewer than for the power-law PSD.

Figure 11.— The phase functions used to determine the importance of multiple scattering.

3.2.3 Total Phase Function Approximation

The example phase functions calculated for power-law PSDs presented above [Fig. 10 (a)] and the analytic solution using (Eqns. 19 and 22) suggest that an analytic function representing in the first two regions, i.e., for scattering angles before the behavior sets in, might be approximated by replacing the term with in the denominator of the approximate individual particle phase function :


where , , and are constants to be determined from fitting either calculations or measurement data. In fact, based on our results for the power-law PSD, we might expect . Fig. 10 (b) compares four of the calculations of with fits, . The fits were carried out using the Levenberg-Marquardt method (Press et al. 1992) for values of below where the curves appear to start bending over to the region. The closeness of the fits confirms the potential utility of this functional form.

3.3. Radiance Profiles

The aureole radiance resulting from the single scattering of starlight by ice crystals uniformly distributed in a plane-parallel layer depends upon their total phase function as follows (e.g., pg. 302, Liou, 2002):


where the single scattering albedo (the ratio of the scattering to the extinction cross section) has been taken as 1, is the exo-atmospheric irradiance of the star, and the small-angle approximation has been applied. When Eqn. 28 applies, it is a relatively simple matter to solve for given measurements of the aureole radiance and , and using a previously determined value of . However, numerous authors (e.g., Dave, 1964; Hovenier, 1971; Korkin et al., 2012) have commented on the limited applicability of the single scattering approximation, e.g., to . These studies, however, typically considered scattering from small particles at all angles. The situation is distinctly different for scattering from large particles at small angles as discussed next.

3.3.1 Multiple Scattering

To investigate the limits of the applicability of the single scattering approximation for use in modeling aureole radiance profiles, we used four phase functions representing the full range of scattering angles. Starting with the least forward scatterering particles, we looked at Rayleigh scattering, which has a simple, analytic phase function. As an example of aerosols we used the phase function retrieved by AERONET (Holben et al., 1998) at NASA Goddard Space Flight Center at 18.636 UT on 8 January 2010. We calculated a phase function for altostratus using the PSD shown by Liou (2002) and the Mie code of Bohren and Huffman (1983). For an example of cirrus we used the phase function calculated by Baum et al. (2005) for a distribution of ice crystal habits and sizes with effective size of . Fig. 11 plots these phase functions on a log-log scale to emphasize the forward scattering direction. This clearly demonstrates the dominance of forward scattering for large particles.

Figure 12.— Monte Carlo calculations of the aureole radiance profiles for cirrus with (“OD = 2.0” in the figure legend) including 1, 2, 3, 4, and 30 scatterings and a fit using Eqn. 29.
Figure 13.— Multiple scattering scattering correction factor for (a) Rayleigh scatterers, (b) aerosols, (c) altostratus, and (d) cirrus for 5 values of (“OD” in the figure legend) from to .

We digress briefly to describe the radiative transfer (RT) method we used to calculate the scattering of starlight through a uniform, plane-parallel particulate layer. The strongly forward peaked phase functions characteristic of ice crystals in cirrus (Fig. 11) tend to cause problems for many RT algorithms. Methods relying on truncating the peak of the phase function at small angles (for example, Joseph et al., 1976; Wiscombe, 1977) work well for hemispheric flux calculations, but are not concerned with details of the angular distribution of radiation and therefore are not appropriate for calculating the profiles of stellar aureole radiance scattered at small angles. We used the successive orders of scattering method (Evans and Marshak, 2005), where the integrals were calculated using a Monte Carlo method (Sobol, 1994). In order to make the calculations more efficient we solved the adjoint problem, tracing photons from the sensor back to the source. Fig. 12 shows an example of calculations for cirrus with . The exo-atmospheric stellar irradiance was taken as that of the Sun () for these calculations. The colored lines labelled “Max Scatterings =” show calculations where the maximum number of scatterings was limited to 1, 2, 3, 4, and 30. (The black curve labelled “Approximation Fit” is discussed later.) The dark blue curve shows the single scattering solution. The fine structure for is numerical noise from the Monte Carlo calculations. This noise decreases somewhat as the maximum number of scatterings increases. The radiance in this region falls off as . As expected, the aureole broadens somewhat as the maximum number of scatterings increases. For the range of angles shown, the solution appears to have converged using a maximum of three or four scatterings.

We define a multiple scattering correction factor as the ratio of the aureole radiance calculated using a maximum of 30 scatterings to that calculated for single scattering. Fig. 13 shows calculations of the correction factors using the phase functions shown in Fig. 11 for = 0.1, 0.2, 0.5, 1.0, and 2.0. The Rayleigh scattering case (panel a) shows a significant increase in radiance as the maximum number of scatterings increases, and has not converged by four scatterings (not shown). The correction factor is nearly constant with , reflecting the shape of the phase function. The aerosol case (panel b) is interesting in that although the phase function is much flatter than that of either altostratus or cirrus, the single scattering solution appears adequate. Unlike the other three cases, which represent conservative scattering, this aerosol case is highly absorptive. The single scattering albedo is only 0.12, which inhibits multiple scattering significantly. The correction factors for the altostatus (panel c) and cirrus (panel d) are close to 1 or nearly constant for less than and , respectively, meaning that the shapes of the aureoles differ only slightly from those of the underlying phase functions.

Notwithstanding the fact that there is typically relatively little difference between single and multiple scattering in the stellar aureoles from cirrus clouds, it is worthwhile to explore correction algorithms, especially for the higher optical depth cases. And it appears appropriate to try to correct for multiple scattering in these cases, especially at the larger angles. The appendix presents a technique for deconvolving the effects of multiple scattering.

Figure 14.— Comparison of the phase function (Fig. 11) used to calculate the aureole radiance for cirrus with (Fig. 12) with the phase function retrieved from the diffraction radiance profile using Eqn. 30 and divided by 2 for clarity.
Figure 15.— (a) Comparison of PSDs retrieved by fitting power-law and exponential PSDs and by application of numerical inversion using first and second difference constraints to the phase function generated from an input model power-law PSD. (b) Comparison of phase functions calculated numerically from the retrieved PSDs with the input model phase function. To separate the curves four of them have been reduced by a factor as indicated in the figure legends.

3.3.2 Aureole Approximation

To the extent that the shapes of aureoles differ only slightly from those of the underlying phase functions, it is reasonable to assume that an analytic approximation similar to the one found for phase functions could be useful in the near forward direction. Consider the aureole profile approximation based on Eqn. 27:


where , , and are constants to be determined from fitting either model calculations or measurement data.

As an example, the black line in Fig. 12 shows an illustrative fit to the red curve representing the multiple scattering aureole radiance profile for a cirrus cloud with . The fit for is obviously very good.

3.3.3 Phase Function Retrieval

The phase function can be retrieved or deconvolved from a measured diffraction radiance profile numerically using the Hankel transform solution described in §7.4, specifically Eqns. 74 and  47:


where is the exo-atmospheric stellar irradiance and is the Hankel transform (Eqn. 54).

Fig. 14 compares the phase function (the blue curve) used to calculate the aureole through a cirrus cloud with (shown in Fig. 12) with a retrieval (the red curve) using Eqn. 30. The two curves overlap so well that we have divided the latter curve by a factor of 2 so that it can be distinguished in the plot. The deconvolution result exhibits numerical noise for , which some simple smoothing could mitigate. To the best of our knowledge this deconvolution technique is new, and the effects of measurement errors and noise have yet to be worked out. The result shown in Fig. 14 lacks such noise and serves simply to confirm that this technique can work under ideal conditions. The reasonableness of the results presented later in the paper for actual stellar aureole data serve to add confidence in the technique under more realistic conditions. Although RT calculations could be used to check whether a retrieved phase function reproduces the measured aureole profiles, we prefer to perform such checks starting with the PSDs derived from the phase functions and seeing how well the phase functions are reproduced.

3.4. PSD Retrieval

Given the phase function and the optical depth there are several ways of attempting to retrieve the PSD. We consider two complementary techniques.

3.4.1 Fitting an Analytic Form

In §3.2.1 and §3.2.2 two parametric approximations were shown for representing PSDs. Given , the power-law PSD has three free parameters (, , and ) and the exponential PSD two ( and , with set to ). A error metric can be defined as:


where is a vector of the free parameters, are the values of the phase function derived (deconvolved) from the aureole measurements, are suitable weights, e.g., setting is equivalent to assuming errors, and are the values of the phase function calculated using the model PSD with parameters and Eqns. 10 to 13. The Levenberg-Marquardt method (e.g., Press et al. 1992) can be applied to minimize with the required derivatives calculated numerically.

As an example fit consider again a model power-law PSD extending from to with . The plot in Fig. 15 (a) compares the input model PSD (the black curve) used to generate the synthetic aureole with retrievals using analytic models involving the power-law and exponential forms. (The constrained retrievals are discussed below.) Not surprisingly the retrieval using the power-law form (the red curve) does quite well, while the retrieval using the exponential form (the orange curve) is only close over part of the size range. As with the multiple scattering deconvolution (§3.3.2) these calculations lack noise and serve simply to confirm that the technique can work under ideal conditions and assuming that the correct underlying form of the PSD is used. Fig. 15 (b) compares the phase functions calculated numerically using the retrieved PSDs shown in panel (a). The fact that the phase function calculated from the retrieved PSD when an exponential form is used does not match the input PSD in this case reassures us that any arbitrary form cannot be assumed. In principle then a variety of PSD models could be used and the one producing the best fit to the input phase function selected. However, the next section presents an alternative approach that has worked well for a variety of other physical retrieval problems.

3.4.2 Constrained Numerical Inversion

To avoid assuming a specific functional form for the PSD as in the previous section we consider the more direct numerical solution of Eqn. 12. The numerical solution of this Fredholm integral equation for as a function of can be made easier (King et al., 1978) by replacing the independent variable with the product of a slowly varying function, , and a more rapidly varying one, e.g., :


where was suggested by the large-particle power-law slopes found by Heymsfield and Platt (1984) and seems to work reasonably well for the current problem. Substituting Eqn. 32 into Eqn. 12 gives:


Following the usual inversion approach (for example, Twomey, 1977), we discretize this equation as follows:


where the form a discrete sequence of values spanning the range of area diameters from to , is suitably defined, e.g., if the values are linearly spaced, and is the average particle density in the size interval about . It is conventional to express Eqn. 34 using matix notation:






Substituting from Eqns. 10 and 11 we find:


where is calculated approximately and analytically by assuming a power-law dependence of the integrand between the end points of each particle size interval.

Figure 16.— (a) Analysis of the components comprising the aureole profile around Capella at 01:10 UT on 17 Dec 2011 and (b) comparison of the diffraction profile fit with the data minus the PSF and background fits.

In general, is not square and hence Eqn. 35 cannot be inverted directly. The least squares solution to this equation, , is problematic because it is under constrained (Twomey, 1977). Frequently, a practical solution makes use of the addition of an appropriate constraint, e.g., that the solution be smooth in some sense. The constraint takes the form of a second equation that is fit simultaneously. A Lagrange multiplier is selected to govern the balance between the least squares errors of the two equations:


To look for a smooth solution, one can constrain (minimize) the first differences by taking as:


or the second differences by taking as:


Press et al. (1992) recommend the following choice for , which balances the two components of the minimization:


where is the trace of the matrix.

Fig. 15 also shows the results of constrained inversions. Both constraints work reasonably well for the power-law model case, although the second difference did marginally better than the first difference constraint.

3.5. Point Spread Function

We model stellar aureoles as the sum of the instrumental/atmospheric point spread function (PSF), starlight diffraction, and background, e.g., diffusely illuminated sky. The PSF describes the spreading of the radiation from a point source, such as a star, and results from a combination of physical processes, e.g., diffraction by optical elements, saturation of detectors, bleeding of charge on the focal plane, and spreading by atmospheric turbulence. Since the latter process tends to be variable with time and space, we approach modeling the PSF using an empirical parameterization. The Gaussian function is frequently used to model PSFs (Racine, 1996):


where and are fitting parameters. The effects of atmospheric turbulence appear to be modeled somewhat better using the Moffat function (e.g., Moffat, 1969; Trujillo et al. 2001):


where , , and are fitting parameters. However, in this work we have found the Gaussian function preferable because its functional form is more distinct from that of the aureole approximation (Eqn. 29). This difference allows the two functional forms to be fit simultaneously without their roles being confused in the process. Therefore we have adopted the Gaussian functional form for modeling PSFs.

Figure 17.— The green curve shows the multiple scattering probability per unit solid angle, , calculated from the diffraction profile fit shown in Fig. 16 (b). The red curve shows the single scattering probability per steradian, , deconvolved from using Eqn. 73.

4. Aureole Profile Interpretation

The process of interpreting aureole profiles begins with separating the three major physical components, the PSF, the stellar diffraction profile, and the background. Next the stellar diffraction profile is deconvolved (for the effects of multiple scattering) to give the phase function, from which the PSD is then derived. The process is illustrated in this section using aureole data for Capella from about 01:10 UT on 17 Dec 2011.

4.1. Diffraction Profile Extraction

The center of Capella on the image was identified, first by locating the brightest pixel and then refining the location using the maximum of a quadratic surface fit to the pixels within of the brightest point. A set of concentric annuli was defined, and the pixel values within each annulus were averaged to produce the aureole radiance values shown as the black diamonds in Fig. 16 (a).

Using Eqn. 44 for the PSF, Eqn. 29 for the stellar aureole, and a constant, , for the sky background, the sum of these components, , is modeled as:


We find the 6 parameters, , , , , , and using the Levenberg-Marquardt method to fit the data points. We tried using the formal RMS flucutations of the pixel values in each annulus to specify the weights used in the fitting process, but this produced error bars that were manifestly smaller than the actual error bars. Consequently, we somewhat arbitrarily adopted error bars for all of the data points to take into account other, systematic errors. We found empirically that the specific choice of weighting functions had relatively little effect on the retrieved diffraction pattern of the aureole profile. Fig. 16 (b) compares the fit for the stellar aureole with data points after subtracting off the fitted PSF and background values. The plot suggests that the aureole approximation (Eqn. 29) has done a reasonable job of representing the diffraction of starlight by the cirrus cloud particles.

Figure 18.— The black curve in panel (a) shows the retrieved phase function calculated by multiplying the single scattering probability shown in Fig. 17 (the green curve) by . The red and green curves in panel (b) show the PSDs calculated from this phase function using numerical inversion with first and second difference constraints, respectively. To check the retrieved PSDs, the red and green curves in panel (a) are the phase functions calculated numerically from the corresponding PSDs in panel (b).

4.2. Phase Function Deconvolution

Next the phase function of the scatterers is derived from the diffraction radiance profile fit, , and the exo-atmospheric stellar irradiance, , using the deconvolution expression of Eqn. 30. Although they do not need to be calculated explicitly, it may be helpful to show two of the intermediate products in this process. The green curve in Fig. 17 shows the multiple scattering probability per unit solid angle, (Eqn. 47). The red curve shows the single scattering probability per steradian, , deconvolved from using Eqn. 73. Note that is both slightly lower and narrower than . This is expected since the effect of multiple scattering is both to increase the total number of photons scattered and to widen their angular distribution. The phase function, the end product of the process, is simply times (Eqn. 58).

4.3. PSD Solution

The last step is to retrieve the PSD from the phase function. We use the numerical inversion technique (§ 3.4.2), trying both first and second difference constraints, and selecting the one that fits better. Fig. 18 (a) compares the phase functions derived from the fits (the colored lines) with the phase function (the black curve) scaled from (the red curve in Fig. 17). The numerical inversion using the second difference constraint (the green curve) produces a slightly better fit than does the one using the first difference constraint (the red curve).

Fig. 18 (b) compares the PSDs corresponding to the fits in Fig. 18 (a). The matrix used in the inversion solution has dimensions . Examination of the 12 eigenvalues for the unconstrained problem shows magnitudes ranging from a maximum of to a minimum of , indicative of ill-conditioning. By contrast, when the second difference constraint is added in using the Lagrange multiplier as indicated in Eqn. 43 the eigenvalues range from to , a considerable improvement that enables a solution to be found. Consideration of the first few eigenvalues of the unconstrained problem, , , , , and , would suggest that, practically speaking, there are perhaps 2 or maybe 3 free parameters. Note that although these PSDs correspond well to the phase function deconvolved from the aureole radiance measurement we cannot say that they are unique, merely reasonable.

Figure 19.— Additional examples of aureole profiles and least-squares fits for the star Capella from the night of 16-17 Dec 2011, annotated by (“OD” in the figure legend).
Figure 20.— (a) Phase functions deconvolved from the diffraction profile fits in Fig. 19 and (b) PSD solutions derived from the phase functions using constrained numerical inversion.
Figure 21.— (a) Example aureole profile data and least-squares fits for the planet Jupiter at 0041 UT on 27 Nov 2011 when the line-of-sight optical depth was 5.1 and (b) PSD calculated using the second difference constraint.

4.4. More Examples

From the datasets we have collected, some 30% of which have been analyzed in detail, we show a few illustrative examples in this subsection. Fig. 19 shows more measurement data from the night of 16-17 Dec 2011. The PSF and background magnitudes are anti-correlated as one might expect, the PSF being brightest when is least, while the background (resulting largely from backscattered “cityshine”) is brightest when is greatest. The diffraction profile radiance is maximized for or so, as is consistent from consideration of the single scattering equation (Eqn. 28).

Fig. 20 (a) shows the phase functions deconvolved from the diffraction profile fits in Fig. 19. Interestingly, the magnitude of the phase function forward scattering peak seems to be loosely correlated with . Fig. 20 (b) shows the PSD solutions derived from the phase functions using constrained numerical inversion. The correlation noted above carries through to the PSDs. We are encouraged by the consistency of the retrievals.

4.4.1 Planetary Aureole Example

Fig. 21 (a) shows an example of aureole measurement data for the planet Jupiter from the early morning of 27 Nov 2011. Note that the aureole profile does not drop down to the background level until about in contrast with the case for the star Capella in Fig. 16 (a), offering a potentially useful increase in particle measurement range. Nevertheless the PSD retrieved using the second difference constraint [Fig. 21 (b)] covers roughly the same size range as the retrievals for the star Capella.

Examination of the first few eigenvalues of the unconstrainted problem, , , and , suggests that there are perhaps 2 free parameters, i.e., no more than in the Capella example shown in §4.3. This may not be surprising given that the least squares fitting procedure used to extract the diffraction profile employs an analytic approximation (Eqn. 29) that has 3 free parameters and yields satisfactory fits.

5. Summary

In summary, we have demonstrated retrieving and analyzing stellar aureole profiles for the purpose of determining the particle size distributions in ice clouds. (1) We have collected a useful stellar aureole database covering different atmospheric conditions and cloud optical depths. (2) Stellar aureole profiles have clearly been detected and measured through a wide range of cloud conditions. (3) The aureole profiles can be followed out to from stars ( from Jupiter). (4) The stellar aureoles from cirrus that we have examined have very distinctive profiles, typically being approximately flat out to a critical angle, followed by a steepening power-law decline with a slope less steep than . (5) The retrieved phase functions cover the range of sizes from to .

In the process, we noted that (6) the relation between the diffraction phase functions of complicated ice crystal habits and the Fourier transforms of their projected areas results in a simple wavelength scaling relationship. We compared the phase functions for a variety of crystal habits with the same size metric. We defined the “area diameter” and found that (7) the phase functions for different habits are very similar for the same area diameter and (8) could be approximated using a simple analytic formula. We found a similar analytic approximation to represent possibly mutltiply scattered diffraction aureole profiles. (9) We utilized this simple aureole approximation along with a Gaussian function for the PSF and a constant for the background to separate the aureole profile data into these components. (10) We developed a technique for deconvolving multiply scattered diffraction profiles to find the phase functions from which (11) we were able to retrieve reasonable PSDs using constrained numerical inversion.

The next steps are (a) to develop an instrument for the routine, automatic measurement of thin cirrus microphysical properties, (b) to deploy it at instrumented sites, such the Atmospheric Radiation Measurement facilities operated by the Department of Energy, and (c) to compare its cirrus microphysical property retrievals with those of other instruments. In this latter category we include specifically (d) independent validation of the retrieved size distributions of cirrus ice crystals by comparison with coincident in situ measurements. In the overall process, a number of instrument improvements are anticipated. For example, we expect that the angular dynamic range of the aureole measurements, which is closely related to the range of particle sizes measured, can be increased by a factor of 2 to 3. Several ways of doing this include (i) reducing the background level by locating the instrument to a darker area, (ii) reducing the PSF width by improving the camera optics, and (iii) using multiple color filters to take advantage of the dependence of diffraction on wavelength (Eqn. 5). In connection with instrument development it will be important to investigate the sensitivity of the phase function retrieval to noise in the radiance measurements.

6. Acknowledgements

The authors acknowledge and thank Rickey Petty and the Atmospheric System Research Program in the Climate and Environmental Sciences Division of the Department of Energy for support through Phase I SBIR grant DE-SC0006325. The authors also thank the anonymous reviewers for their recommendations, particularly their encouragement to take into account the effects of multiple scattering. Finally, the authors thank Brent Holben and his staff for their work in establishing and maintaining AERONET and especially the site at Goddard Space Flight Center.

7. Appendix - Multiple Scattering Deconvolution

Small-angle, multiple scattering has been studied extensively for well over half a century, particularly in the context of particle or electron scattering from beams (e.g., and references therein, Bethe, 1953). Mathematically the problem involves convolution in the forward direction, or deconvolution in the case of determining the single scattering pattern from measurements (e.g., Misell and Burge, 1969; Jones and Misell, 1967). Assuming cylindrical symmetry in the aureole profiles, we adopt an approach that takes advantage of the Hankel transform (Bracewell, 2000).

Consider plane-wave radiation from a distant point source, e.g., a star, of irradiance incident normally on a uniform, plane-parallel, particulate layer, e.g., cirrus. Let be the total radiance (luminance) scattered by the particulate layer at scattering angle . The total or multiple scattering probability per unit solid angle into an annular ring at scattering angle is:


It is convenient to discretize the scattering angle domain, i.e., replacing continuous functions with vectors. Let be the vector element representing at scattering angle , where and covers the range from . Then Eqn. 47 can be written as:


Consider the equation for single scattering in the forward direction from radiaion incident normally on a uniform, plane-parallel, particulate layer (e.g., pg. 302, Liou 2002). The singly scattered aureole radiance is:


where is the optical depth measured normal to the layer and is the phase function normalized to integrate to , is the phase funnction normalized to integrate to , and the single scattering albedo is omitted assuming conservative scattering. Or, dividing by we find that:


7.1. Poisson Model

Approximate multiple scattering as a Poisson process (e.g., Ning et al. 1995). Then is the sum over the probabilities of all orders of scattering:


where represents the 2-dimensional convolution operator and is the Poisson probability for scatterings:


and is the average number of scatterings, which is equal to the optical depth . For the particles of interest the forward scattering radiance forming the aureole is much more intense than scattering at other angles. Consequently the integral in the convolution operator can be approximated by an integral that goes to infinity:


where is an azimuthal angle and the small angle approximation, , has been employed.

It is convenient to apply the Hankel transform to Eqn. 51 so that the convolution operation is replaced by multiplication in the transformed domain using the convolution theorem for Hankel transforms. Let be the Hankel transform operator:


where is a function to be transformed and is the Bessel function of order 0. The inverse Hankel transform is also . Let represent the Hankel transform of so that Eqn. 51 becomes: