1. Introduction

Atmospheric particles, e.g., ice crystals and aerosols, play a significant role in radiative transfer and remote sensing by scattering and absorbing solar radiation and by terrestrial thermal emission. Various numerical methods have been developed to study the single-scattering properties of these particles [1,2]. Among the existing methods, the method of separation of variables (SOV, e.g., the Lorenz-Mie theory) [3] and the T-matrix method [4,5] are analytical or semi-analytical methods. The SOV is applicable to both spheres and spheroids. Although the T-matrix method is, in principle, applicable to arbitrarily shaped particles, its use is normally limited to axially rotational symmetric particles because the simulation is computationally efficient with the particle symmetry. The discrete dipole approximation (DDA) [6–11], the finite-difference time domain (FDTD) method [8,12–15], and the pseudo-spectral time domain (PSTD) method [16,17] share similar areas of applicability and are numerically rigorous methods based on solving Maxwell’s equations for electromagnetic scattering by arbitrarily shaped particles. Specifically, the DDA solves an electromagnetic integral equation in the frequency domain by discretizing a scattering particle in terms of a number of electric dipoles; whereas, both the FDTD and PSTD solve Maxwell’s curl equations in the time domain. The FDTD and PSTD differ in the numerical procedure chosen to perform the spatial derivative of electromagnetic fields. Specifically, the traditional FDTD is based on the second-order finite difference method, and the PSTD employs a Fourier pseudo-spectral method.

The DDA and FDTD have been extensively studied [6–14] and systematically compared for simulating light scattering by spheres for the size parameters x up to 80 (x = 2πr/λ, where r is the radius of the sphere and λ is the incident wavelength) and the real part of refractive index m up to 2 [18]. The numerical performances of the two methods are found to strongly depend on the refractive index of the scattering particles; the DDA is faster for smaller m, and the FDTD is more computationally efficient for larger values of m. The “cross-over” refractive index between the two methods is at approximately 1.4 [18].

Stemming from the traditional FDTD method, the PSTD applies the pseudo-spectral method, instead of the finite difference method, to calculate the spatial derivatives in Maxwell’s equation [19]. The method has been employed to simulate light scattering by non-spherical ice crystals by Chen et al. [16] and has been improved by Liu et al. [17]. In comparison with the finite difference method, the pseudo-spectral method has higher orders of accuracy and smaller numerical dispersion errors. The higher order of accuracy means that a errors smaller than a given tolerance level can be achieved by PSTD with coarser spatial resolution in terms of number of grid points per wavelength, than is needed by FDTD [19]. This is one way in which PSTD methods can save CPU time. In time-stepping methods a larger time-step can be taken if coarser resolution may be used, which is the second way in which PSTD methods can save time. However, the relative strengths of the PSTD in comparison with other numerical methods than FDTD are not clearly known.

Our study compares the use of the PSTD and DDA methods for the numerical simulation of light scattering by dielectric particles. Specifically, we focus on spheres and spheroids because the accuracy of the results can be well quantified by comparison with their counterparts simulated from the Lorenz-Mie theory [3] and the T-matrix method [4,5]. The comparison is performed with the same prescribed accuracy criteria for both methods and covers a broad range of size parameters (up to 100 for spheres and 50 for spheroids) and non-absorbing refractive indices. The specific numerical implementations and techniques used will be detailed in Section 2, and the comparison and results are discussed in Section 3. Section 4 reports the conclusions drawn based on the study.

2. The DDA and PSTD implementations

Code ADDA, developed by Yurkin and Hoekstra [20], is a widely used DDA implementation for light scattering simulations. Using a cluster of processors, the ADDA can simulate light scattering by particles much larger than the incident wavelength, and the reported maximum size parameters for spheres with refractive indices of 1.05 and 1.2 were 320 [20] and 130 [21], respectively. We used ADDA v.0.79 with the default settings for dipole polarizability (lattice dispersion relation) and iterative method (quasi minimal residual method). The convergence criterion of the iterative solver was set to be 10⁻³; larger than the default value (10⁻⁵) but sufficient to reach the accuracy required by this study. These code settings are identical to those used in [18] and correspond to the mainstream DDA. In particular, the code settings are similar to the default settings of the DDSCAT [22], another widely used implementation of the DDA. Thus, we focus on the practical performance of the ADDA (with default settings), and not on the best theoretically possible. We also believe the conclusions will be valid for the DDA method in general. However, we briefly discuss possible consequences of using the latest (1.1) version of the ADDA at the end of Section 3. We use spheres and spheroids, both of which have symmetry utilized by the ADDA to halve the computational time compared to nonsymmetric shapes [20].

The PSTD Fortran90 implementation, based on the original algorithm described by Chen et al. [16], was improved and applied to large sized particles by Liu et al. [17]. The resultant PSTD implementation was parallelized with the OpenMP Application Program Interface (API) that supports shared-memory to model parallel computation. Unlike the FDTD, which defines the electric and magnetic field components at the edges and face centers of Yee cells [12], the PSTD employs a centered grid scheme that specifies all field components at the grid cell centers. The Gibbs phenomena, encountered in previous FDTD algorithms, was eliminated by truncating the high wavenumber terms in the pseudo-spectral simulations. Consequently, the applicability and accuracy of the PSTD was substantially improved and allowed the PSTD implementation to be applied to spheres with size parameters up to 200 at three ice refractive indices with real parts close to 1.3 [17]. The advanced package “fastest Fourier transform in the west” (FFTW) [23] was used to perform the fast Fourier transform (FFT) and inverse FFT for the Fourier pseudo-spectral method. A Gaussian pulse multiplied by a cosine term (with a center frequency equal to the frequency of interest) was used as the incident wave. The uniaxial perfectly matched layer (UPML) absorbing boundary condition was applied to truncate the spatial domain [24, 25]. The electromagnetic field values in the time domain were transformed to the frequency domain through a discrete Fourier transformation. The Kirchhoff surface integral equation was used to transform the near-field values into their far-field counterparts [26]. The discrete Fourier transformation involves a numerical integration in time, and the time interval of integration sufficient for convergence must be determined by experimentation because of its dependence on the particle size and the index of refraction. For the results in this study, we found an integration time interval between 4T_p and 5T_p to be sufficient, where T_p is the time required for the pulse to cross a distance equal to the largest diameter of the scattering particle.

We simulated the single-scattering properties of spheres and spheroids with different sets of x and m by using the PSTD and DDA, and compared the results with the exact solutions to quantify the accuracies of the two numerical methods. Considering the axially rotational symmetry of the scattering particles, one simulation of a linearly polarized incident wave was sufficient to yield the 4 by 4 phase matrix P. The phase matrix in one scattering plane was calculated with the scattering angle θ varying from 0° to 180° in steps of 0.25°. The extinction efficiency Q_ext and the normalized phase function P₁₁(θ) are the two major quantities in the estimation of the accuracy of the two methods, but the asymmetry factor g and phase matrix element P₁₂(θ) will also be compared with the analytic solutions. With a specified accuracy criteria for the Q_ext and P₁₁(θ), the computational time required to achieve the criteria became the most meaningful parameter to describe the overall performance of the methods.

The PSTD and DDA discretize the scattering particles with grid points (PSTD) and dipoles (DDA), but their computational times are dependent on the number of grid points or dipoles in the computational domain and on the number of time steps for the PSTD or iterations by the iterative solver of a large linear system for the DDA. For a particle with a fixed size, the computational domain scales cubically with the spatial resolution, i.e., the number of grid intervals or dipoles per wavelength (λ/Δx). The accuracy of each method increases with an increase in the spatial resolution. We increase the spatial resolution until the required accuracy criteria are achieved, namely that the relative errors (RE) of Q_ext are less than 1%, and the root mean square relative errors (RMSRE) of P₁₁(θ) are less than 25% (same as in [15]). This procedure should not be considered as a one-fits-all solution, because it is not suitable for certain applications. In particular, it results in over a 50% relative error in backscattering intensity (see Section 3). However, both methods can produce smaller errors at the expense of extra computational resources, and the chosen procedure describes the general trends.

The DDA is the preferred method for optically soft particles (particles with refractive indices near 1) [18,21], and may outperform specialized methods, like the discrete sources method, for axisymmetric particles [27]. Therefore, our comparison focused on refractive indices larger than 1.2, and for spheres, we used a real m ranging from 1.2 to 2.0 in steps of 0.2. The minimum size parameter for the comparison was 10. To keep the computational time manageable and achieve the accuracy criteria, especially for the DDA simulations, the upper limit of the x we considered decreased from 100 to 40 as m increased from 1.2 to 2.0. The exact sets of x and m involved in the computation are shown in Tables 1 and 2 . Moreover, the comparison was extended to spheroids with realistic refractive indices of ice (m = 1.312 + 1.489 × 10⁻⁹i at 0.532 μm [28]) and mineral dust (approximate m = 1.55 + 0.001i at the visible wavelength [29]). The size parameter of a spheroid was specified in terms of its equivalent-volume sphere. Aspect ratio values a/b of 0.5 (corresponding to oblate spheroids) and 2 (prolate spheroids) were used, where a was the equatorial radius and b the semi-length of the symmetric axis. Spheroid size parameters ranging from 10 to 50 in steps of 10 were chosen for the simulation. The propagation direction of the incident field coincided with the symmetry axis. All simulations were carried out using a single node containing 8 64-bit 2.8 GHz processors (a cluster at the Texas A&M Supercomputing Facility). It should be noted that for such a shared-memory configuration, parallelization scheme, the ADDA (MPI) is less efficient than the OpenMP used in the PSTD, because the MPI was originally designed for distributed-memory (multi-node) hardware. However, we estimate that the effect due to the difference in the parallelization scheme should not exceed 20% of the computational times and, hence, does not influence the final conclusions.

Table 1. The Comparison of Numerical Performances of the PSTD and ADDA Methods for Spheres with Different x and m.

View Table | View all tables in this article

Table 2. Same as Table 1 but for Some Accuracy Results

View Table | View all tables in this article

3. Results

Table 1 lists both the computational parameters and the simulation results and illustrates the numerical performance of the PSTD and DDA. In addition to m and x, Table 1 includes the spatial resolution, computational time, REs of Q_ext, and RMSREs of P₁₁. Indicated within parentheses are the results of cases in which the PSTD or DDA failed to reach the prescribed accuracy even with a very fine spatial resolution. Computations too time-consuming (taking more than 4 days) to reach the prescribed accuracy for the DDA are marked as “N” in the table. The PSTD simulations covered all the sets of x and m chosen for the study and achieved the prescribed accuracy in 24 of the 28 total pairs. To achieve the prescribed accuracy criteria, the spatial resolutions used by the PSTD varied from 10 to 30 without systematic dependence on x or m. The accuracy values for the DDA do show significant sensitivity to m. The DDA used spatial resolutions smaller than 10 for a refractive index of 1.2 and increased monotonically to 40 for m = 2.0. As an extra verification of the DDA results, we note that they agree well with those of [18], where an earlier version of the ADDA code (version 0.76) was used. In particular, we obtained almost identical values of spatial resolution and simulation error for the two ADDA versions when using the same x and m (x ≤ 60, 40, and 10 for m = 1.2, 1.4, and 2.0, respectively).

With the same accuracy criteria achieved by the PSTD and DDA simulations, the behavior of both with respect to the computational time show substantial variations for different x and m. With size parameters up to 100, the PSTD simulations were finished within 9.0 × 10⁴ seconds (i.e., 25 hours), and the most time-consuming simulation was for a sphere with x = 100 and m = 1.4. Furthermore, neither the efficiency nor the accuracy of the PSTD was significantly influenced by an increase of m. However, the computational time used by the DDA simulations increases dramatically with both particle size and refractive index due to the simultaneous increases of the spatial resolution, computational domain, and iteration number. For example, for m = 1.2, only a few seconds were required for spheres with x = 10 and 20; whereas, a sphere with x = 80 took 7.3 × 10⁴ seconds (over 20 hours). When m is larger than 1.4, the DDA encounters difficulties with respect to both efficiency and accuracy. A sphere with a size parameter of 30 and m = 2.0 took 5.1 × 10⁵ seconds (almost 6 days) and obtained Q_ext with RE of 2.0% and P₁₁ with RMSRE of 55%. The DDA achieved the prescribed accuracy criteria for spheres at a size parameter of 30 for m = 1.6 and only a size parameter of 10 for m = 1.8 and 2.0. The DDA did not achieve convergence for most large (spheres with x > 60 and m = 1.4 or 1.6 and x > 40 and m = 1.8 or 2.0) cases (7 spheres out of 28).

On the contrary, the DDA was very efficient for spheres with small m and x and was one to two orders of magnitude faster than the PSTD. However, a critical size parameter existed above which the PSTD outperformed the DDA for small refractive indices (1.2 or 1.4) and the value of x decreased from 80 to 30 as the refractive index increased from 1.2 to 1.4. For m larger than 1.4, the PSTD became more efficient for all size parameters in the range from 10 to 60 and was almost two orders of magnitude faster than the DDA for spheres with x larger than 30.

In the (x, m) domain of this study, Fig. 1 clearly illustrates the strengths of both the PSTD and DDA. In the figure, the dark region corresponds to the parameters of x and m for which the PSTD simulations outperformed the DDA, and the light region indicates the area where the DDA was more efficient. The PSTD area of high performance is at the top right of the (x, m) plane, and the DDA counterpart is at the lower left. The region in Fig. 1 is made up of only the (x, m) sets we used for comparison, but both methods have been applied to larger sizes without significant loss of accuracy [17,21]. Comparing Fig. 1 with results of [18], we can conclude that the PSTD is similar to the FDTD when compared with the DDA, except for an increase in the relative performance of the PSTD with an increase in size parameter, even for e.g., m = 1.2.

 

Fig. 1 The relative performances of the PSTD and DDA on the (x, m) plane. The light area corresponds to the DDA-preferred region and the dark area to the PSTD-preferred region.

Download Full Size | PDF

Table 2 lists other numerical optical property errors with which to compare the performances of the two methods, but no separate accuracy criterion was prescribed for these quantities. The table includes the REs of the asymmetry factor g, the maximum REs of P₁₁, the REs of P₁₁ at 180° (i.e., backscatter), and the root mean square absolute errors (RMSAE) of P₁₂/P₁₁. Overall, the differences between the four errors given by the PSTD and DDA are relatively small and were dependent on x or m in the following manner:

A further comparison of the numerical accuracy of the two methods, the P₁₁ of spheres with the same size parameter of 40 and different refractive indices, is illustrated in the left panels of Fig. 2 . From the top to the lower panel, the refractive index is increased from 1.2 to 2.0 in steps of 0.2. Shown in the figure are the exact solutions given by Mie theory (dashed lines) and the results of the PSTD (dark lines) and DDA (light lines) simulations. The relative errors of P₁₁ are shown in the right panels of Fig. 2. The RMSREs of P₁₁ for the spheres, as given by the PSTD and DDA, range from 15% (m = 1.4 for DDA) to 33% (m = 1.6 for DDA). For spheres with x = 40, the DDA simulation achieves the prescribed accuracy only for refractive indices of 1.2 and 1.4. When the refractive indices reach to 1.8 or 2.0, numerical convergence is not achieved for the DDA simulations at x = 40, hence no DDA results are shown. However, the PSTD results exceed the 25% criterion only for the sphere with m = 2.0, and the RMSRE of P₁₁ is 26%. The relative errors of the PSTD and DDA are smaller than 30% at most scattering angles, but became significant, even as large as 100%, near the angles where sharp troughs or peaks occurred in the phase function. In comparison, the REs of the phase functions given by the PSTD and DDA were of the same order for spheres with refractive indices of 1.2 and 1.4. At a refractive index of 1.6, the REs, simulated by the DDA, of the backward scattering at scattering angles larger than 140° became 50% or larger. The spheres with refractive indices of 1.8 and 2.0 could only be simulated by the PSTD. The REs are then comparable to those with small m and indicate the weak influence of m on the PSTD simulation accuracy. When particle size distributions or different orientations of non-spherical particles are taken into consideration in practical applications, the strong oscillations in the phase function are smoothed. Thus, both methods will provide much more accurate and reliable phase matrix elements. The comparison shown in Fig. 2 indicate the results for forward scattering are apparently more accurate than those of backward scattering.

 

Fig. 2 Comparison of the PSTD and DDA results with the exact (Lorenz-Mie) solutions for simulation of P₁₁ (left panel) for spheres with x = 40 and a refractive index increase from 1.2 to 2.0 from the top to lower panel. The right panels show the relative errors.

Download Full Size | PDF

For the same spheres with size parameters of 40 and refractive indices ranging from 1.2 to 2.0, the left panels of Fig. 3 show the ratios of P₁₂/P₁₁ as functions of scattering angles and the right panel the absolute errors. The RMSAEs of the ratios given by both the PSTD and DDA are between 0.13 and 0.25 (from Table 2), and the absolute errors at most scattering angles are less than 0.2. Again, the results given by the PSTD and DDA for m = 1.2, 1.4, and 1.6 are comparable. For refractive indices up to 1.8 and 2.0, the PSTD simulations give results with similar accuracy to those with small m.

 

Fig. 3 Same as Fig. 2 but for P₁₂/P₁₁ (left panels) and the absolute errors (right panels).

Download Full Size | PDF

Similar to Table 1, Table 3 shows the comparison between the PSTD and DDA for spheroids, and the results compared with the exact solutions given by the T-matrix method [4,5]. With the refractive index of ice, both the PSTD and DDA fail only at a spheroid with x = 50 and a/b = 2. When the spheroids had the refractive index of mineral dust, the PSTD achieved the accuracy criteria for all sizes and aspect ratios except the one with x = 50 and a/b = 0.5. However, the DDA simulations could only be carried out for x smaller than 30 and achieved the criteria for sizes less than or equal to 20. As expected, the PSTD outperforms the DDA for large spheroids with x = 50 when m = 1.3117 + 1.489 × 10⁻⁹i, and was the preferable method for all spheroids with m = 1.55 + 0.001i, except the one with a/b = 2 and x = 10. The relative performance of the two methods shows no dependence on the spheroid aspect ratio. Generally, the refractive indices of ice at different wavelengths have a real part of approximately 1.3, and those of aerosol particles, i.e., dust and back carbon, are 1.5 or larger. Thus, our comparison suggests the DDA to be suitable for numerical simulations of ice crystals with size parameters smaller than 50, whereas the PSTD is more efficient and more accurate for ice crystals with size parameters larger than 50 and aerosol particles of all sizes. The PSTD and DDA results with respect to the REs of g, maximum REs of P₁₁, REs of P₁₁(180°), and RMSAEs of P₁₂/P₁₁ for spheroids are similar to those of spheres and will not be included here.

Table 3. Same as Table 1 but for Spheroids with Size Parameters from 10 to 50, Aspect Ratios of 0.5 and 2.0, and Refractive Indices of 1.313 + 1.489 × 10⁻⁹i and 1.55 + 0.001i.

View Table | View all tables in this article

The left panels of Fig. 4 show P₁₁ of the spheroids with x = 30 and the right panels the relative errors of the PSTD and DDA results compared with the T-matrix solutions. The aspect ratios and refractive indices are labeled in the figure. In general, the PSTD and DDA results both had excellent agreement with the T-matrix results, although the relative errors became significant at a few scattering angles around the troughs or peaks in P₁₁. For the spheroid with x = 1.55 + 0.001i and a/b = 2.0, the DDA gave the RMSREs of P₁₁ to be 31%, which is larger than the criterion, but Fig. 4 shows the relative errors to be larger than 50% only at the scattering angles around 80° and the ones larger than 130°.

 

Fig. 4 Same as Fig. 2 but for spheroids with x = 30, aspect ratios of 0.5 and 2.0, and refractive indices of 1.312 + 1.489 × 10⁻⁹i and 1.55 + 0.001i.

Download Full Size | PDF

The ratios of P₁₂/P₁₁ for the spheroids with the same size parameters and the absolute errors are illustrated in Fig. 5 . The absolute errors of P₁₂/P₁₁ are no more than 0.2 at most scattering angles. For the spheroid with a/b = 2.0 and m = 1.55 + 0.001i simulated by the DDA (lower panels), the errors became larger than 0.5 at the scattering angles that had relative errors of P₁₁ larger than 50%. From Figs. 4 and 5, we notice that, for spheroids, the PSTD results of P₁₁ and P₁₂/P₁₁ are slightly more accurate than those of the DDA.

 

Fig. 5 Same as Fig. 4 but for P₁₂/P₁₁ (left panels) and the absolute errors (right panels).

Download Full Size | PDF

Before drawing our final conclusions, we reflect once more on the parameters of the DDA simulations. As noted before, the paper is based on an older version of the ADDA and (almost) default settings, but we performed a limited set of simulations (for two spheres) with the current development version of ADDA (1.1b6 as of May 1, 2012), trying different DDA formulations and iterative solvers to maximize performance.

For the sphere with x = 10 and m = 2.0 the best result was obtained with the filtered coupled dipoles (FCD) formulation of the DDA [30,31]. The version decreased the number of iterations by 25% and largely improved the accuracy such that a spatial resolution of 20 was sufficient for a prescribed accuracy threshold. The resulting computational time was 140 s – 14 times faster than the DDA with the default settings, but still 3 times slower than the PSTD. For the sphere with x = 80 and m = 1.2, the FCD formulation halves the number of iterations but has little effect on the accuracy. The best performance, however, was achieved by using a CSYM iterative solver [32], which resulted in an almost four times smaller number of iterations, and a computational time of 2.0 × 10⁴ s, only two times larger than that of PSTD.

4. Conclusions

A systematic comparison between the PSTD and DDA for ice crystal and atmospheric aerosol light scattering computations was made by using the parallelized implementations of the two methods on the same multi-processor hardware, although we have no reason to believe that the relative performance is significantly affected by the hardware used. For spheres, size parameters up to 100 and refractive indices up to 2.0 were used, and for spheroids, two aspect ratios and two realistic refractive indices of ice and dust were used with equivalent-volume size parameters up to 50. The same prescribed accuracy criteria were required for the extinction efficiency and the phase function, and the computational time was used as the key parameter to evaluate and compare the two methods. The DDA was more economical for numerical simulations of spheres with small refractive indices and small size parameters; whereas, the PSTD was more economical for large x and m. The critical size parameter, above which the PSTD outperformed the DDA, decreased from 80 to 30 as the refractive index increased from 1.2 to 1.4. The PSTD was more CPU-efficient and applicable to a wider range of x when the refractive index was larger than 1.4. Similar conclusions were obtained for the spheroids. Furthermore, the overall accuracy of the asymmetry factor, backscatter, and linear polarization given by the PSTD and DDA were in agreement.

The implementation of each of the two compared methods has been substantially enhanced through our efforts. For instance, we found that the latest formulations of the DDA can decrease the required computational time by an order of magnitude. However, this is expected only to shift the boundary between the methods in the (x, m) plane but not to principally affect the conclusion of this comparison. Moreover, potential users are advised to test (and fine-tune) these and other light-scattering methods for their particular applications before performing large-scale simulations. Finally, we note that the comparison was performed only for real refractive indices or those with negligible imaginary parts. Significant absorption is known to largely improve the convergence of the iterative solver in the DDA [33]. Therefore, a comparison between the PSTD and DDA for absorbing, including metallic, refractive indices is an interesting topic for future research.