1. Introduction

1.1. Ice crystals

Microphysical properties of cirrus clouds [1, 2], such as the scattering asymmetry parameter and single scattering albedo of ice crystals, affect climate feedback processes [3]. These micro-physical properties are governed by ice crystal shape, size, refractive index, and the mass of ice crystals per unit volume of air. Simple ice crystals are non-spherical, often hexagonal columns or plates, sometimes with irregularities such as inclusions or voids [4, 5]. Scattering by more complicated crystal shapes can also be calculated [6]. Besides their importance for atmospheric radiative transfer, they are also responsible for intricate optical displays [7–9] and polarization effects [10]. Their single scattering properties play an important role in remote sensing in visible and infrared light as well as in microwaves [11]. Light scattering properties of ice crystals have been typically calculated in the geometric optics (GO) approximation by ray-tracing algorithms [8, 12–17]. When crystal dimensions are comparable to the incident wavelength λ, diffractive effects have been estimated using circular cylinders [18, 19] and combinations of several techniques [6, 20].

Analytic solutions to Maxwell’s equations are known only for special geometries such as spheres, spheroids, or infinite cylinders, so approximate methods are in general required. The Discrete Dipole Approximation (DDA) is one such method. The basic theory of the DDA has been presented elsewhere [21]. Conceptually, the DDA consists of approximating the target of interest by an array of polarizable points. Once the polarizability tensors are specified, Maxwell’s equations can be solved to arbitrary accuracy for the array, and the scattering problem of interest can be solved to arbitrary accuracy by reducing the dipole separation d → 0 (within computational limits). Developed originally to study scattering from isolated, finite structures such as dust grains, the DDA was extended to treat singly or doubly periodic structures [22, 23]; we also generalized the scattering amplitude matrix and the 4×4 Mueller matrix to describe scattering by such periodic targets. The accuracy of DDA calculations using the open-source code DDSCAT was demonstrated by comparison with exact results for infinite cylinders and infinite slabs. A fast method for using the DDA solution to obtain fields within and near the target was also presented [24].

 

Fig. 1 Hexagonal ice crystal with vertex-to-vertex diameter D. There are six side faces and two end (basal) faces. The c-axis is parallel to the sides and perpendicular to the end faces. Prisms are rotated randomly around the c-axis. In the GO approximation, halos discussed in this paper peak at the minimum deviation angle ζ_min for which once-refracted light is propagating parallel to a side face before the second refraction.

Download Full Size | PDF

In this paper we use the DDA for periodic targets to calculate light scattering properties of infinite hexagonal ice crystals with size parameter x ≡ πD/λ ≤ 400, where D is the diameter. Such infinite hexagonal prisms approximate light scattering by atmospheric ice crystals with a large aspect ratio [15]. We show that finite diameter ice crystals have scattering properties that deviate systematically from geometric optics; these deviations can be used to infer the cystal diameter D from measured halos.

1.2. Geometry and relation to atmospheric optics

Water ice crystals form hexagonal columns, with the column symmetry axis parallel to the crystalline c-axis. Long ice crystals fall with their c-axis oriented horizontally [7, 8]. Such ice crystals will be randomly rotated about the c-axis and about a vertical axis perpendicular to the c-axis [8]. In our calculations we include the effects of random rotations around the c-axis, but we consider only two values of the angle Θ between the incident beam and the c-axis. Real atmospheric halos also include a contribution from reflections by hexagonal prism end faces. End faces are not included in the present DDA calculations because we assume infinitely long crystals to allow use of the periodic boundary condition version of the DDA to study very large values of D.

We consider an infinitely long hexagonal prism, with light incident at an angle Θ relative to the c-axis. If the sides are numbered j = 1,...,6 in order, then rays entering through side j and exiting through side j + 2 will be refracted, with the emergent ray scattered by an angle θ_s. We will refer to the scattering associated with light passing through sides j and j + 2 (this is path 35 in [8]) as the “22° feature”. (The term “22° halo” is a standard term for the halo which is generated by light passing through sides j and j + 2 of randomly-oriented crystals. However, here we deal with preferentially oriented crystals.) The photon momentum parallel to the c-axis is unchanged. As the hexagonal prism is rotated around the c-axis, the minimum scattering angle θ_s,min occurs when the ray passes through the prism traveling parallel to one of the sides; for randomly-rotated prisms, the halo peak is at θ_s,min. Let ζ be the deflection of the ray projected on the basal plane. The minima of ζ and θ_s are [25, 26]

2. Discrete dipole approximation method

2.1. Scattering

The scattering properties of targets with 1-dimensional translational symmetry are described by the generalized Mueller matrix elements $S_{ij}^{(1\text{d})}$ relating the Stokes parameters of the scattered wave to the Stokes parameters of the incident plane wave [23]. For unpolarized incident light, Stokes parameters I and Q for the scattered wave are proportional to $S_{11}^{(1\text{d})}$ and $S_{21}^{(1\text{d})}$, respectively, and the polarization $P=Q/I=S_{21}^{(1\text{d})}/S_{11}^{(1\text{d})}$.

The $S_{ij}^{(1\text{d})}$ are directly related to differential scattering cross sections per unit length. In particular,

For each allowed θ_s, there are two values of ζ satisfying (4): let these be 0 ≤ ζ₁ ≤ π and ζ₂ = 2π − ζ₁. The usual differential scattering cross section is related to $S_{11}^{(1\text{d})}$ by

2.2. Accuracy

We have performed calculations using DDSCAT [21, 23] for regular infinite hexagonal prisms with size parameters up to x = πD/λ = 400. Birefringence is fully included in the DDA treatment. We use complex refractive index

The DDA calculations were done for 10 different rotations of the target around the c-axis, at intervals of Δβ = 3°; with the symmetry of the target, these 10 orientations are used to average over rotations of the target around the c-axis.

Calculations for larger size parameters required extensive computer resources and several days of computer time. We used the OpenMP option in DDSCAT. Parameters for the dipole arrays are given in Table 1, including the dipole spacing d, and the number N of dipoles in a single hexagonal layer (the “target unit cell”).

Table 1. Parameters for DDA calculations

View Table

We use DDSCAT with conjugate gradient algorithm “GPBICG” [28, 29] to iteratively solve the system of linear equations. Iteration terminates when a specified error tolerance tol is achieved; tol is listed in Table 1, as well as N_iter, the average number of iterations required to reach the specified error tolerance for a single orientation and incident polarization. For a given x, the required CPU time is approximately proportional to N_iter × N ln(N). For the largest problem, x = 400, with N > 3 × 10⁶ dipoles, the error tolerance was raised to tol = 5 × 10⁻⁴ in order to lower N_iter and reduce the required CPU time. Table 1 also gives the total CPU time required for convergence (for 10 orientations and two polarizations) for single-precision calculations using 4 cores running at 2.6 GHz.

Previous studies comparing DDA results with exact solutions for spheres [21] indicate that DDA results will be accurate provided |m|k₀d < 0.5, with the DDA results converging to the exact solution as d → 0. The |m|k₀d < 0.5 condition is satisfied in the calculations presented here, so that the DDA should provide an accurate representation of the infinite prism.

As a check, we have repeated the calculations for x = 200 for smaller interdipole spacing d. Figure 2 shows S₁₁ and S₂₁ for 15° < θ_s < 35° for DDA realizations with |m_e|k₀d = 0.486 and 0.328; the two calculations are in excellent agreement. For a given target, we expect the DDA results (if fully converged) to be exact in the limit d → 0, with the error varying approximately in proportion to d. Better accuracy can be obtained by reducing d, but this becomes computationally expensive for large x.

 

Fig. 2 Discrete dipole approximation results for incidence angle Θ = 90° and x = 200, averaged over rotations around the c-axis (sampled at Δβ = 3° intervals), for 2 values of d (see Table 1). Vertical dashed lines show the geometric optics halo angles θ_s,min for ordinary and extraordinary polarizations. (a) Muller matrix element $S_{11}^{(1d)}$ (total halo intensity $I\propto S_{11}^{(1d)}$ (b) Fractional polarization P. The short dashed line indicates the polarization (P = 0.037) of the 22° halo peak in the GO limit.

Download Full Size | PDF

3. Results

3.1. Halo features

Figure 3 shows d²C_sca/d cosθ_sdL calculated for Θ = 90°, for selected values of the scattering parameter x ≡ πD/λ. In the GO calculation, the inner edge of the 22° halo is at 21.93° and 22.04° for ordinary and extraordinary polarizations. In the wave optics calculation for x = 25, the peak occurs at ∼18°, well interior to the inner edge of the GO halo. As x is increased, the halo power becomes more concentrated, and closer to the expected halo angle 22°. For x = 50 we find a double-peaked halo, centered near 22°. For x ≥ 100 the peak becomes increasingly narrow and pronounced as x is increased.

 

Fig. 3 (1/D)d²C_sca/d cosθ_sdL as a function of scattering angle θ_s for hexagonal prism with light incoming perpendicular to c-axis, averaged over prism rotations around the c-axis. (a) Full range of θ_s; (b) zoom on 15° < θ_s < 30°. DDA results are shown for selected values of x ≡ πD/λ, as well as GO results (see text). Vertical dashed lines show the expected position θ_s,min of the inner edge of the 22° halo for ordinary and extraordinary polarizations in the ray-tracing limit.

Download Full Size | PDF

The peaks occuring every 6° from 33° to 177° in Fig. 3(a) are produced by specular reflection. Each orientation β contributes specular reflection peaks at specific scattering angles. Because we sample rotations around the c-axis at Δβ = 3° intervals, the peaks from different orientations are separated by Δζ = 2Δβ = 6° (for Θ = 90°, Δθ_s = Δζ). These peaks, with diffractive broadening (δθ_s)_diff ≈ λ/D = 180°/x, become narrower as x is increased, and hence are most conspicuous for the x = 400 case. For finer rotational sampling Δβ → 0, these specular reflection peaks would blend into a continuum.

Also shown in Fig. 3 is $S_{11}^{(1\text{d})}$ obtained from a GO calculation using a ray-tracing code written by A. Macke [13], with the scattering binned with Δθ_s = 0.02°. The conspicuous spike in the GO results at θ_s = 120° comes from rays that enter through side j, undergo internal reflections from sides j + 1, j + 2, j + 4, and j + 1, and exit through side j + 5 with a deflection of exactly 120°. With only ∼0.1% of the scattered power in this feature, diffractive broadening washes out the corresponding scattering in the DDA results. Overall, the DDA results appear to be clearly converging to the GO limit as x → ∞.

To demonstrate that these effects are not confined to the special case of Θ = 90°, Fig. 4 shows results of both DDA and GO calculations for incidence angle Θ = 60°. The DDA calculations take the birefringence into account exactly, but the ray-tracing code [13] is for isotropic material. To approximate the results for weakly anisotropic H₂O ice, we take a weighted sum of separate GO calculations for refractive indices n_o and n_e, with weights 0.5(1 + cos² Θ) and 0.5sin² Θ for the ordinary (E ⊥ ĉ) and extraordinary (E ‖ ĉ) cases. This treatment is correct for Θ = 90°, but only an approximation when Θ < 90°.

 

Fig. 4 Similar to Fig. 3, but for oblique incidence Θ = 60°. (a) Full range of allowed scattering angles θ_s (0 – 120°); (b) zoom on 18° < θ_s < 33°. For Θ = 60° the GO cusps are at 24.896° and 25.017° for ordinary and extraordinary rays. Note the good agreement between GO and the x = 400 results for θ_s > 26°.

Download Full Size | PDF

For Θ = 60° incidence, the GO caustics are at θ_s,min = 24.90° and 25.02° for n_o and n_e (see Eq. 2). The DDA results again show convergence to the GO solution as x → ∞. As for Θ = 90°, the halo is clearly visible for x = 100, although appreciably broadened [18], with the halo increasing in peak intensity and in sharpness as the size parameter increases.

When the sun is not at the zenith, the angle between the incident light and the c-axis of the horizontally-oriented crystals will range from Θ = 90° to 90° − θ_z, where θ_z is the zenith angle. Thus the present results, for only Θ = 90° and Θ = 60°, can only be compared to the 22° halo for the sun at the zenith (θ_z = 0).

3.2. Particle sizing: power spillover index

For finite size crystals, ray optics ceases to be a good approximation, as is evident from comparison of the GO and wave optics results in Figs. 3 and 4. For x ≤ 400, an appreciable amount of scattered light spills into the “forbidden” (so-called “shadow”) region with θ_s < θ_s,min. The power scattered into this forbidden zone can be used to constrain the diameter of the ice crystals doing the scattering.

We define a quantitative “power spillover index”

Figure 5(a) shows the power spillover index Ψ as a function of x ≡ πD/λ for Θ = 90°. For x ≈ 100 (D ≈ 18μm), we see that Ψ ≈ 0.4 – much larger than the value for GO. It is evident that measurement of the scattered power interior to the GO caustic at θ_s,min = 22° provides information about the particle size. The diagnostic is most straightforward if the light source is at the zenith, so that all of the scattering is for Θ = 90°. If the light source is at zenith angle θ_z > 0, then one must take an appropriate average over Θ values between 90° and 90° − θ_z.

 

Fig. 5 (a) The “power spillover index” Ψ as a function of x for Θ = 90°. Ψ is a measure of scattered power interior to 21.9° (see Eq. 10). Line connecting calculated points is only to guide the eye. (b) Full width at half-maximum (FWHM) for the scattering peak near 22°, as a function of 1/x = λ/πD. Line connecting calculated points is only to guide the eye.

Download Full Size | PDF

If the sun (or full moon) is the source of illumination, then the differential scattering cross section must be convolved with the 0.53° diameter uniform brightness disk. This solar convolution increases the GO value of Ψ from 0.006 to 0.058. For x ≤ 400 the solar convolution has only a minimal effect on the wave optics results for Ψ [see Fig. 5(a)].

Of course, interpretation of atmospheric scattering must recognize that particles other than long ice columns may also be contributing to scattering at θ_s < θ_s,min, so that a metric like Ψ provides only a lower limit on the diameter of the halo-producing hex columns.

3.3. Particle sizing: FWHM

Another particle-size diagnostic is provided by the FWHM of the scattering peak near 22°. In the GO limit, this peak is a caustic, with dC_sca/dθ_s → ∞ and FWHM = 0. For finite D, the peak dC_sca/dθ_s is finite, and FWHM > 0.

The 22° halo is produced by radiation that enters through a side face j (see Fig. 1) with width perpendicular to the beam w = (D/2)cos(30° + ζ_min/2) = 0.377D for ζ_min = 22°. Diffraction through a slit w has

It should be stressed that the FWHM in Fig. 5(b) is for the special case of Θ = 90° (i.e., sun at the zenith, and horizontally-oriented hex columns). For the more general case of solar zenith angle θ_z > 0, the FWHM must be recalculated with appropriate averaging over Θ values between 90° and 90° − θ_z. Integration over a range of Θ values will increase the FWHM, with the effects most pronounced for large values of D.

3.4. Polarization

For x → ∞, application [30] of the Fresnel equations shows that the halo light is linearly polarized. If birefringence is neglected, for Θ = 90° light passing through faces j and j + 2, the fractional polarization at θ_s = θ_s,min is

Können and co-workers [31, 32] have discussed the expected birefringence peak in the polarization for geometric optics. For Θ = 90°: (1) For θ_s < 21.93°, the 22° halo is not present, and P ≈ −0.25 from weak small-angle scattering by oblique reflections off the edge faces (see Fig. 3). (2) For 21.930° < θ_s < 22.037° the 22° halo is ∼100% polarized with E ⊥ ĉ, i.e., parallel to the scattering plane (P > 0). (3) Just beyond 22.037°, both polarizations contribute to the 22° halo, but the contribution with E ‖ ĉ will dominate at and just outside the 22.037° caustic, and the net polarization will be perpendicular to the scattering plane (P < 0). (4) At larger angles, θ_s > 24°, the polarization has a value close to P_h ≈ 0.037 from eq. (13). Figure 6 shows the polarization over the 20–25° range for x → ∞ and Θ = 90°, calculated for x → ∞ using ray-tracing for random rotations of the target around the c-axis.

 

Fig. 6 Polarization P = S₂₁/S₁₁ as a function of scattering angle θ_s for hexagonal prism with unpolarized light incident perpendicular to the c-axis (Θ = 90°), for x = 100, 200, 300, 400. (a) 19° < θ_s < 25°; (b) zoom into 21° < θ_s < 24°. Also shown (labelled x = ∞) are GO results, including birefringence. Vertical dashed lines show location of halo edge for ordinary and extraordinary polarizations. Horizontal dashed line shows P_h = 0.037 expected for GO neglecting birefringence (see Eq. 13). See text for discussion.

Download Full Size | PDF

Können and Tinbergen [31] discussed the diffractive broadening of the 22° polarization peak produced by birefringence, and showed how multiwavelength polarization observations could be used for particle sizing.

Figure 6 shows the polarization obtained from the full solutions to Maxwell’s equations, including birefringence. Our exact results for x = 100 – 400 show pronounced oscillations with angle, but have generally positive P in the 20–22° region where the scattered light is dominated by “spillover” of the 22° halo into this forbidden region. For 21.930° < θ_s < 22.037° where GO predicts P ≈ 1, the x = 400 solution (corresponding to diameter D ≈ 70μm) has P ≈ 0.10 at 550 nm, comparable to P = 0.10 at 615 nm and 0.16 at 435 nm observed by Können, Wessels, and Tinbergen [33]. From the progression of P values at x = 200, 300, 400, it appears that D ≈ 90μm (x ≈ 500 for λ = 0.55μm) would likely reproduce the polarizations reported by Können et al.

4. Summary

The public-domain code DDSCAT [21] was used to calculate light scattering by infinite hexagonal ice prisms for size parameter of up to x = 400 in the discrete dipole approximation. Birefringence is fully included in the discrete dipole approximation calculations. The 22° halo like feature is predicted for size parameter as low as x = 100, but with an angular structure very different from the halo calculated using geometric optics. The halo peak becomes more pronounced and narrower for increasing size parameters.

The results show that accurate solutions of Maxwell’s equations can be used to relate observed properties of atmospheric optics displays, including some aspects of halos [8], to the particle size. The power spillover index Ψ [Fig. 5(a)] and the FWHM [Fig. 5(b)] are both measurable and directly related to the diameter D of the hexagonal columns. Future research may employ the methodology demonstrated in this paper to study the effects on halo properties of various factors, such as porosity [5] or non-evenness of side edges [32]; such target geometries can be treated by DDSCAT.