Light scattering by hexagonal columns in the discrete dipole approximation

Scattering by infinite hexagonal ice prisms is calculated using Maxwell's equations in the discrete dipole approximation for size parameters up to x=400. Birefringence is included in the calculations. Applicability of the geometric optics approximation is investigated. Excellent agreement between wave optics and geometric optics is observed for large size parameter in the outer part of the 22 degree halo feature. For smaller ice crystals halo broadening is predicted, and there is appreciable"spillover"of the halo into shadow scattering angles<22 degrees. Ways to retrieve ice crystal sizes are suggested based on the full width at half-maximum of the halo, the power at<22 deg, and the halo polarization.


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 microphysical 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][8][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][13][14][15][16][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, 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.
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]. 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.

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] where n is the real part of the refractive index. At λ = 550nm and T = −25C, H 2 O ice has n o (E ⊥ĉ) = 1.3112, n e (E ĉ) = 1.3126 for "ordinary" and "extraordinary" polarizations [10,27]. For Θ = 90 • , the halo peak is at θ s,min = 21.930 • and 22.037 • for the ordinary and extraordinary polarizations, respectively.

Scattering
The scattering properties of targets with 1-dimensional translational symmetry are described by the generalized Mueller matrix elements S (1d) i j 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 are directly related to differential scattering cross sections per unit length. In particular, where C sca is the scattering cross section for unpolarized light, k 0 ≡ 2π/λ , and 0 ≤ ζ < 2π is the azimuthal angle around the "scattering cone", with ζ = 0 for forward scattering. The scattering angle θ s ≥ 0 is related to ζ by for given Θ, the largest allowed scattering angle is θ s,max = arccos(1 − 2 sin 2 Θ).

(5)
After averaging over target rotations around the c-axis, the present problem has reflection symmetry, S (1d) and for ζ = π

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 weak absorption [Im(m) = 10 −5 ], introduced to accelerate the conjugate gradient solution, is small enough that the scattering results are not affected appreciably. 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").
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 6 dipoles, the error tolerance was raised to tol = 5 × 10 −4   Table 1) 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 0 d < 0.5, with the DDA results converging to the exact solution as d → 0. The |m|k 0 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 11 and S 21 for 15 • < θ s < 35 • for DDA realizations with |m e |k 0 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.  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.
Also shown in Fig. 3 is S (1d) 11 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 2 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 2 Θ) and 0.5 sin 2 Θ for the ordinary (E ⊥ĉ) and extraordinary (E ĉ) cases. This treatment is correct for Θ = 90 • , but only an approximation when Θ < 90 • .
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).

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" The upper limit 21.9 • adopted for the integral in the numerator in Equation (10) is chosen to be just interior to the inner edge of the halo for the ordinary polarization, for Θ = 90 • . In the GO limit, the only power interior to 21.9 • comes from processes other than refraction through faces j and j + 2, thus Ψ is small: Ψ = 0.0060.    ≈ 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 .
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.

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 We have also smoothed the results using a 0.53 • diameter disk [see Fig. 5(b)] but the effects are minimal for x ≤ 400. 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.

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 = 0.037 for n = 1.312.
directed parallel to the scattering plane (P > 0), i.e., perpendicular to the halo arc. However, birefringence also contributes to the polarization. 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.
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.

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.