1. Introduction

Aerosols are common in air and can be divided into natural-source and anthropogenic ones. Natural-source aerosols mainly come from sandstorms, volcanic eruptions, evaporation of sea waves, forest burning, and other natural phenomena. Anthropogenic aerosols are mainly from fossil and biomass combustion, transportation, and industrial emissions. From the industrial era to nowadays, content of anthropogenic aerosols such as sulfate and carbonaceous aerosols has kept increasing, which greatly increased the global average of aerosol particles [1]. Atmospheric aerosols have important influences not only on human health, but also on climate and environment. The influence of aerosols on climate system has both direct and indirect aspects [2]. On the one hand, aerosol particles can directly absorb and scatter solar (short-wave) and earth (long-wave) radiation, which directly affects the radiation balance of earth-atmosphere system. On the other hand, aerosol particles can act as condensation nuclei of clouds, thus changing the microphysical structure of clouds such as the size and density of cloud droplets, affecting precipitation and radiation characteristics of clouds, and ultimately indirectly affect the climate. Therefore, the simulation of radiative properties of aerosols is of great significance to climate research.

However, due to the great variability of shape, size and chemical composition of aerosol particles, there are still great uncertainties in the study of aerosol's influence on radiation. For example, black carbon (BC) aerosols have significant positive radiative forcing due to their strong absorption of solar radiation [3]. However, through the aging process, BC tends to be coated by secondary aerosol species (e.g., organics and sulfate), and form complicated mixed structure [4,5], resulting in complex microphysical and optical properties. This leads to one of the largest uncertainties in aerosol radiative forcing. Spherical model is the simplest and most common aerosol model, the optical properties of which can be accurately solved with the Lorenz-Mie theory. However, the shapes of real aerosol particles, like the mentioned coated BC, are always non-spherical and irregular. The use of spherical models necessarily leads to large errors in the inversion of aerosol optical thickness and the estimation of climatic effects [6,7]. Therefore, accurately simulating the optical properties of non-spherical aerosol particles like coated BC is of great importance.

With the rapid development of computers, many numerical methods for simulating particle optical properties have been developed and widely used. According to the operating principles, they can be divided into time-domain and frequency-domain ones, and the latter – into those based on surface and volume discretization. Time-domain based methods include the finite-difference time-domain method (FDTD) [8,9], the pseudo-spectral time domain method [10,11], etc. Surface-discretization methods are commonly semi-analytical (T-matrix), which include the multiple sphere T-matrix (MSTM) [12] and the extended boundary condition method [13]. Fully-numerical boundary equation method has also been recently applied to complicated aggregates [14]. Volume-discretization frequency-domain methods include the discrete-dipole spproximation (DDA) and the invariant imbedding T-matrix (especially, with the latest developments [15]), and the DDA is the one most commonly used for aerosols [16–18]. The DDA is especially popular for calculating the optical properties of particles with complex geometric structure, especially those composed of large number of simple subunits (monomers), such as BC aggregates [19–23]. Coated BC features not only the complex geometric structure, but also a huge size difference between different components. Simply refining discretization to accurately describe the shape of the monomers will lead to a waste of computational time and memory due to the coat which has simple geometric structure, but much larger size.

In this paper, we improve the efficiency of the DDA by combining it with the effective-medium approximation (EMA). Originally, the EMA was based on the assumption that a heterogeneous particle can have a homogeneous counterpart possessing similar scattering and absorption properties. Thus, the EMA has been widely used to model complex heterogeneous substances as being homogeneous [24–30]. Yang and Liou [8] applied the EMA to smooth the voxels on the particle edge in the FDTD simulations to decrease the corresponding discretization errors, i.e., it was used for a single relatively large voxel rather than for a whole particle. A similar approach has been discussed in the DDA framework [21,31–33], which is related to the notion of super-dipoles or meta-atoms, that are used in simulation of metamaterials [34–37]. Moreover, a more complicated “weighted discretization” has been proposed [38,39], but it has never been used in any open DDA code due to the complexity of implementation.

Importantly, previously EMA-DDA combination has only been applied to aggregates or even simpler-shaped homogeneous particles. Here we apply it to coated aggregates in order to solve the conflict between the complex internal structure and large size. In doing so we investigate the DDA + EMA for particles of different shapes, sizes, and refractive indices in terms of both accuracy and computational resources. This study is organized as follows. The methods including the DDA, the EMA, and the reference ones (for accuracy evaluations) are introduced in Section 2. Section 3 presents the performance of the DDA + EMA on spheres, including the sensitivity tests for different size parameters and refractive indices. Section 4 applies the DDA + EMA to coated BC particles, which model is based on transmission electron microscope observations, and gives the guidelines for using the DDA + EMA method. Section 5 concludes the study.

2. Methods

The DDA discretizes the scatterer using dipoles (or voxels), and, for a given particle, its efficiency and accuracy are largely determined by the so-called dpl (dipoles per lambda/wavelength). The dipoles should be much smaller than both the incident wavelength and the characteristic particle length scales. Thus, the dpl from 10 to a few tens are normally considered, but much larger values may be required for particles smaller than the wavelength. The ADDA implementation (v1.3b4) [40] is used in this study, which is parallelized with the MPI to be run on computer clusters. The code is used with the default parameters, including point-dipole interaction, lattice-dispersion relation for polarizability, and ${10^{ - 5}}$ threshold for the convergence of the iterative solver. Finally, all simulated particles are considered in random orientation, using the default orientation-averaging parameters in ADDA. This may seem redundant for spheres, but removes the dependence on the incidence direction that is significant for crude discretizations.

Whenever the DDA accuracy is questionable, the straightforward solution is to increase the dpl, in other words, to decrease the dipole size. This leads to huge simulation time and required memory [41,42]. The same applies to the coated BC aggregates [20,43,44], for which discretization of the monomers significantly smaller than the wavelength (leading to dpl up to a few hundreds) is combined with a large coating around the aggregate. A redundantly large number of such dipoles to fill the coating leads to waste of computational resources. Importantly, using adaptive discretization (i.e., dipoles of different sizes) is not an option, since it will break the FFT acceleration inside the DDA [17]. We further use EMA to enable larger dipoles with satisfactory accuracy.

Figure 1 illustrates the application of the EMA to boundary dipoles during discretization of a sphere (the two-dimensional cross section is shown). The blue grids indicate the dipoles (cubical voxels), while the red circles represent spherical boundaries. The left panel shows the original DDA discretization (without the EMA), and the right one is for the DDA + EMA discretization. During the standard discretization, dipoles are considered as either particle (dark gray) or vacuum (white) according to whether their centers are inside or outside of the particle. This leads to the well-known staircasing effect, i.e., the difference in geometry between the model and the real particle (so-called, shape errors). By contrast, the EMA smoothes the boundary, since the boundary dipoles are assigned the values of refractive index in between that for the particle and vacuum, depending on their volume fraction. The latter is represented by the fill color in the right panel of Fig. 1; the darker the dipole is – the larger is the volume fraction, and the larger is the effective refractive index. This smooth transition (from dark to light) between the BC material and the vacuum (or a coating) is expected to reduce the shape errors.

 

Fig. 1. Illustration of the standard (left) and EMA-weighted (right) DDA discretization for a cross section of a sphere. Fill color represents the refractive index of the dipoles.

Download Full Size | PDF

The important side effect of the EMA discretization is the increase of the total number of non-void dipoles. It has only minor effect on the computational cost, since the circumscribing rectangular box (relevant for FFT-accelerated operations in the DDA [18]) usually stays the same. However, it notably increases the total volume of the discretized particle, especially for aggregates. This volume is used by ADDA to calculate the volume-equivalent radius ${r_{\textrm{eq}}}$ and corresponding geometric cross section ${C_{\textrm{eq}}} = \pi r_{\textrm{eq}}^2$ of the particle; hence, these values are not accurate. Thus, any optical cross section ${C_\textrm{X}}$ (where X is ext, abs, or sca) is directly usable from ADDA output, while the efficiencies ${Q_\textrm{X}} = {C_\textrm{X}}/{C_{\textrm{eq}}}$ are not. Therefore, to obtain accurate ${Q_\textrm{X}}$ we divide ${C_\textrm{X}}$ by independently computed ${C_{\textrm{eq}}}$.

The only remaining design choice is that of a particular EMA rule. There are numerous approximations for the effective refractive index m of mixtures; we further test two most popular ones, namely the Maxwell Garnett (MG) [45] and the Bruggeman (BR) theories [46]. The corresponding refractive indices are given by:

Here, ${m_1}$ and ${m_2}$ are the refractive indices of two components of the mixture (BC and coat/vacuum in this study), and their corresponding volume fractions in each dipole are ${f_1}$ and ${f_2}$, respectively (${f_1} + {f_2} = 1$). Note that the two materials with m₁ and m₂ in the MG are asymmetric, and they are understood as the inclusion (particle) and host medium, respectively. While the two materials in the BR are symmetric, i.e., they can be freely interchanged.

In addition, there are known Wiener bounds of the effective relative permittivities [47]:

However, they were originally derived for the case of non-absorbing constituents, for which ${m_{\textrm{W}1}} \,>\, {m_{\textrm{W}2}}$ and the bounds are also called maximum and minimum, respectively. Their generalization to complex input permittivities is known as the Bergman-Milton bounds [44]. The latter describe a region in the complex plane of permittivity, bound by two circular arcs, passing through ${\varepsilon _{\textrm{W}1}}$, ${\varepsilon _{\textrm{MG}}}$, ${\varepsilon _{\textrm{W}2}}$ and ${\varepsilon _{\textrm{W}1}}$, $\varepsilon _{\textrm{MG}}^{\textrm{com}}$, ${\varepsilon _{\textrm{W}2}}$, respectively, where $\varepsilon _{\textrm{MG}}^{\textrm{com}}$ is complementary Maxwell Garnett value obtained by Eq. (1) with indices 1 and 2 interchanged. Finding the bounds for real and imaginary parts of m, corresponding to this region of permittivity is straightforward but cumbersome. Fortunately, for the two cases considered in this paper, these bounds are always realized by either ${m_{\textrm{W}1}}$ or ${m_{\textrm{W}2}}$, i.e., the Bergman–Milton bounds coincide with the Wiener bounds in terms of real and imaginary parts of the refractive index.

These bounds (W1 and W2) are shown in Fig. 2 as a function of ${f_1}$, together with MG and BR functions. The BC refractive index ${m_1}$ is assumed to be $1.8 + 0.6\textrm{i}$ [48], while ${m_2}$ is 1 and $1.53 + 0.006\textrm{i}$ [49] for vacuum and coating (left and right columns), respectively. As expected, the real and imaginary parts of refractive indices increase monotonically from that of vacuum/coating to the BC value (with the only exception of $\textrm{Re}\,{m_{\textrm{W}2}}$ for ${f_1}$ close to 1 in the right column). However, the bounds W1 and W2 are interchanged relative to each other when mixing BC with the coating. This highlights the complexity of Bergman–Milton bounds – even when they coincide with the Wiener bounds, the latter are not necessarily maximum and minimum ones in the same order as for real ${m_1}$ and ${m_2}$. The difference between the MG and BR values is always much smaller than the attainable range according to Bergman–Milton bounds. This difference can be considered significant only for the imaginary part when mixing with vacuum and (to a lesser extent) for the real part when mixing with the coating. Anyway, both the MG and BR are within the range given by the Bergman–Milton bounds, and are further tested in simulations.

 

Fig. 2. Comparison of refractive indices given by the two effective medium approximations and Bergman–Milton bounds, which for these cases coincide with the Wiener bounds. The left and right rows correspond to the cases of BC mixing with vacuum and the coating, respectively.

Download Full Size | PDF

Straightforward implementation of the EMA into the ADDA code will require significant changes in its internal structure to enable arbitrary refractive index for each of the dipoles. However, the current code (with minor adjustments during the compilation) supports inhomogeneous particle comprised of up to 255 different materials. Moreover, approximate computation of ${f_1}$ can easily be done using a subgrid with s times smaller dipoles, accounting for occupancy of each sub-voxel as in the standard discretization. The number of potential discrete values of ${f_1}$ is then ${s^3} + 1$, and so is the number of different values of m (for each pair of materials). We created a separate script, which prepares shape files and a list of refractive-index values for ADDA, similarly to the one used previously in [21]. We used $s = 2$ in simulations, resulting in 8 and 16 different values of m for homogeneous particles and coated BC aggregates, respectively, excluding the surrounding vacuum. Moderately larger values of s (e.g., 3 and 4) can be used with almost the same computational time, but it was tested to have negligible effect on the final accuracy. Also, the script can be easily adapted to other shapes, like prism or ellipsoid, which is left for future research.

For spherical particles, the simulation accuracy is normally validated by comparing with results of the exact Lorenz–Mie (LM) theory. We use the code by Bohren and Huffman [50] in this study. As for coated BC particles, an efficient numerical solution for idealized aggregates of spheres, MSTM 3.0 [51], is used as a reference. However, the MSTM requires constituent spheres not to overlap with their boundaries. Therefore, when studying the realistic coated BC aggregates (for moderate size of coating), we manually remove the aggregate monomers that intersect the coating boundary. This highlights another advantage of this study and the DDA in general, as it allows consideration of arbitrary complicated (realistic) particle shapes (although this is also true for other volume-discretization methods [52]).

All simulations are carried out on a single node with 28 2.5 GHz processors and 128 GB memory.

3. Performance for spheres

We first discuss the performance of the DDA + EMA for homogeneous spheres. We consider two refractive indices of the particle: $1.6 + 0.0\textrm{i}$ and $1.2 + 0.6\textrm{i}$. Figure 3 shows the comparison of the extinction efficiency ${Q_{\textrm{ext}}}$ computed with the original DDA (“Ori”), the DDA with the EMAs (MG or BR), and the Lorenz–Mie theory for spheres ($m = 1.6 + 0.0\textrm{i}$) with size parameter x varying from 1 to 15. Note that the original spherical shape is also specified by an input file instead of being generated inside ADDA to avoid any artefacts due to the so-called volume correction. The latter slightly adjusts the size of each dipole to ensure that the volume of the dipole of the dipole representation of the particle is exactly correct [40]. As expected, the accuracy of the DDA improves with refining discretization (increasing dpl). Three widely spaced values of dpl (4, 8, and 16) are used to illustrate the general trend. The dpl = 4 results are sensitive to the use of orientation averaging (in contrast to a single default orientation), as well as to the DDA formulation. For instance, the default LDR polarizability formulation depends on the incident direction, but may be inadequate for such small dpl [40]. Still, the employed simulation parameters are directly generalizable to complex shapes, while the detailed study of various DDA formulations is outside the scope of this paper.

 

Fig. 3. Comparison of the DDA and Lorenz–Mie results (extinction efficiency) for spheres ($m = 1.6 + 0.0\textrm{i}$) with size parameter x varying from 1 to 15. The DDA results are with (solid lines for the MG results and dash-dotted lines for the BR results) and without (colored dashed lines) the EMAs and with dpl values of 4 (blue), 8 (green), and 16 (red).

Download Full Size | PDF

Moreover, the dpl is a fully meaningful parameter for particles larger than the wavelength (i.e., $x\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel> \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }\pi $) [18], including the aggregates discussed in the next section. For significantly smaller particles, a more relevant parameter is the number of dipoles along the sphere diameter ${n_D} = x \cdot \textrm{dpl}/\pi $, since its small value necessarily implies poor description of the particle shape (even for large dpl). Specifically, the benefit of using the EMA is significant for ${n_D}\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }4$, which corresponds to $x\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }3$, $x\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }1.5$, and none of the shown data points for dpl = 4, 8, and 16, respectively. In other words, for sufficiently large x the results with or without the EMA differ much less than the error of each of them (or difference of the DDA values with different dpl). By contrast, the improvement can be tremendous for small x and dpl. For instance, the relative error is reduced from 445% to –15% for $x = 1.5$ and dpl = 4. The MG results are slightly better than the BR results at dpl = 4, but the difference between the two is generally much smaller than their difference to the original DDA. The first conclusion is, thus, that EMA helps for ${n_D}\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }4$, while the choice of specific EMA formula is not important. Also, the use of EMA for the cases, where no systematic accuracy improvement is obtained, is not recommended since it incurs certain increase of computational time, as explained above.

Figure 4 shows similar results for spheres with $m = 1.2 + 0.6\textrm{i}$. Here, results of extinction, absorption, and scattering efficiency (${Q_{\textrm{ext}}}$, ${Q_{\textrm{sca}}}$, ${Q_{\textrm{abs}}}$), as well as single scattering albedo (SSA) are shown in the four panels. Three different dpl values and the line colors and styles are the same as in Fig. 3. As expected, all DDA variants converge to the exact solution with increasing dpl. For extinction efficiency, the results are similar to Fig. 3. For dpl = 8 and 16 the EMA leads to significant accuracy improvement for sufficiently small particles ($x \le 3$ and $x \le 2$, respectively) corresponding to ${n_D}\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }10$. For the smallest ${n_D} \approx 5$ with dpl = 16 the improvement is about 60 times – compare relative error of 15% for $x = 1$ with the original DDA to 0.23% after applying the MG. At these smaller values of x, results of EMA with dpl = 8 (green solid and dash-dotted lines) are close to that of original DDA with dpl = 16 (red dashed line), but requires about 8 times smaller computational time and memory. For dpl = 4 the EMA also shows the largest accuracy improvement for $x \le 2$. However, in contrast to Fig. 3, the improvement is present over the whole range of x for this dpl, although it is only 2-3 times for $x \ge 7$. The results for ${Q_{\textrm{sca}}}$, ${Q_{\textrm{abs}}}$, and SSA show similar behavior, but the errors of the original and EMA DDA are comparable (but still of opposite signs) for $x \ge 5$ with dpl = 4. Moreover, in the same range of x for both dpl = 8 and 16, the original DDA outperforms the EMA, which may be related to the non-monotonous convergence of the former (other examples of such convergence are discussed, e.g., in [16]). In particular, the errors of the original DDA for dpl = 8 is smaller than that for dpl = 16 for these scattering quantities, which is an artefact of the specific shape and refractive index. In addition, different from the Fig. 3 ($m = 1.6 + 0.0\textrm{i}$), the BR method performs slightly better than the MG one in this case.

 

Fig. 4. Similar to Fig. 3, but for $m = 1.2 + 0.6\textrm{i}$. The extinction, scattering, and absorption (${Q_{\textrm{ext}}}$, ${Q_{\textrm{sca}}}$, ${Q_{\textrm{abs}}}$), and single scattering albedo (SSA) are shown.

Download Full Size | PDF

The performance of the DDA + EMA for spheres with different refractive indices is illustrated in Fig. 5. Here we consider only the MG formulation of the EMA and fix $x = 2$, which is typical for or a bit larger than BC monomer size parameters at ultraviolet, visible, and shortwave infrared wavelengths. EMA results with dpl = 8 (middle column) are compared with that without EMA for the same and twice larger dpl (left and right columns, respectively). Application of the EMA at dpl = 8 systematically reduces the errors by at least 2 times for almost all values of the refractive indices. The notable exception is ${Q_{\textrm{abs}}}$ for weakly absorbing but strongly refracting particles (i.e., large $\textrm{Re}\, m$ and small $\textrm{Im}\,m$). Moreover, the errors for the DDA + EMA are comparable to that of the original DDA with twice larger dpl.

 

Fig. 5. Relative errors of extinction, scattering, and absorption efficiencies given by the DDA simulations with (middle column, dpl = 8) and without (left and right columns – dpl = 8 and 16, respectively) the EMA (MG) treatment as functions of the complex refractive index.

Download Full Size | PDF

To conclude this section, the EMA (both MG and BR) significantly improves the DDA accuracy for single homogeneous spheres when the discretization is crude; the maximum value of ${n_D}$ is up to 10, depending on the refractive index. In these cases, the application of the EMA is often comparable to doubling the value of the dpl in the original DDA, thus reducing the simulation time by about an order of magnitude for the same accuracy. The EMA effectively smoothes the surface of a sphere, reducing large-scale roughness that is the most important factor for small particles with relatively small refractive index. More specifically, the EMA (especially, the MG that is based on Clausius–Mossotti polarizability) is expected to accurately reproduce the total polarizability of a small sphere. On the one hand, the EMA cannot reproduce the intricate details of the electric field near the boundary present in the Lorenz–Mie theory (which are, e.g., the reason for oscillations of ${Q_{\textrm{ext}}}$ with x in Fig. 3). Thus, the EMA is not beneficial for sufficiently large x, m, and dpl, where weighted discretization [38] may be more successful (left for future research). On the other hand, the simple spherical shape allows other ways to optimize the discretization, such as ensuring the overall grid dimension to be equal to sphere diameter (in contrast to Fig. 1) or the abovementioned volume correction. We do not test these options, since they are not easily generalizable to complex shapes.

4. Performance for coated BC

In this part, we discuss the application of the DDA + EMA to randomly oriented coated BC. For BC particles, typical diameter of a spherical monomer ranges between 10 and 100 nm. A model was built based on transmission electron microscope observations in Fig. 6(a) [53]. The fractal dimension (D_f) of BC aggregates indicates whether the particle is lacy or compact. The D_f of fresh BC aggregates is commonly smaller than 2, but it can be close to 3 for aged BC. Referring to the microscope images, the model D_f is assumed to be 2.8. Moreover, we assume the fractal prefactor of 1.2 and used 200 monomers with diameter of 100 nm, and refractive index of ${m_1} = 1.8 + 0.6\textrm{i}$. The incident wavelength is fixed at 500 nm. With all these parameters set, the fractal aggregates are generated by a diffusion-limited aggregation algorithm. Then it is partly coated by a sphere with diameter of 800 nm and refractive index of ${m_2} = 1.53 + 0.006\textrm{i}$. The values of refractive indices are the same as in Section 2. Moreover, as also explained there, the monomers intersecting the coating boundary are removed to ensure the robust operation of the MSTM; only 137 monomers remain.

 

Fig. 6. Coated BC model (b) constructed based on the transmission electron microscope observations (a) [53] and the approximated geometry by the DDA with (d and f) and without (c and e) the EMA with a dpl = 4 (c and d) and 8 (e and f). In the discretized shapes, the yellow and brown parts denote BC and the coating, while the shades of light yellow and pink are used for the mixture between BC and either vacuum or the coating, respectively.

Download Full Size | PDF

The added simplification of such exclusion is that we almost avoid the mixture of three components in a single cubical voxel, thus making it sufficient to use Eqs. (1) and (2), although their generalization to larger number of components is readily available. Moreover, we do not use the EMA for the interface between the coating and vacuum since no significant accuracy improvement is expected for such large sphere size, while the computational time slightly increase (see Section 3). This conclusion is also supported by a limited number of tests (data not shown). Thus, the EMA discretization is performed by first setting ${m_2}$ to the value of either coating or the vacuum if the dipole center is inside or outside the coating, respectively. Then this refractive index is used for the dipole, if it does not intersect any of BC monomers, and effective refractive index is computed otherwise.

The resulting coated BC model is shown in Fig. 6(b), while their standard discretizations using dpl = 4 and 8 are shown in Figs. 6(c) and (e), respectively. The corresponding EMA discretizations are given in Figs. 6(d) and (f), respectively.

Figure 7 shows the comparison of the DDA results (extinction, scattering, and absorption efficiencies) with and without the EMAs and MSTM results of the case shown in Fig. 6. The black dashed lines are for MSTM results. The colored solid and dash-dotted lines are for the DDA results with the EMAs (MG and BR), and the colored dashed lines are for those without the EMAs. Compared to the previous section, the curves are smoother and the accuracy improvement due to the EMA is present over the wider range of dpl, up to 15. The latter, however, corresponds to ${n_D} = 3$, which is comparable to that in Fig. 1. This again highlights the relevance of ${n_D}$ for estimating the performance of EMA-DDA. A smaller value (3 instead of 4) can be explained by averaging of shape errors between multiple monomers, as well as addition of constant contribution from the coating. Coming back to the potential use of EMA for coating/vacuum interface, it may help for similar small number of dipoles (3 or 4) per the coating diameter. But this leads to dpl value of about 2, which in turn implies unacceptably large errors for any variant of the DDA, making such improvement meaningless in real applications.

 

Fig. 7. The comparison of the DDA (with the MG, BR, and without any EMA) and MSTM results for the cases shown in Fig. 6. Red, green, and blue colors correspond to absorption, scattering, and extinction efficiencies, while line style denotes the DDA variant.

Download Full Size | PDF

For ${Q_{\textrm{ext}}}$ results with the MG-DDA, the relative error is within 1% for dpl > 5 (and within 0.5% for most of this dpl range). While for the original DDA, the maximum relative error is 5% (for dpl = 6) and the errors fall within 1% only starting at dpl around 17. If such accuracy level is desired, EMA leads to much smaller dpl (3–4 times). For ${Q_{\textrm{sca}}}$ and ${Q_{\textrm{abs}}}$, the conclusions are the same. Specifically, for dpl = 5 (i.e., ${n_D} = 1$) the relative error of ${Q_{\textrm{ext}}}$, ${Q_{\textrm{sca}}}$, and ${Q_{\textrm{abs}}}$ are reduced from 4.7%, 13.5%, and 5.9% to 0.3%, 0.9%, and 1.9%, respectively, by using the MG. This is not surprising, since the volume of the cubical dipole is then almost twice the volume of the spherical BC monomer (although this error is then averaged over multiple monomers). It is also worth mentioning that the MG method is superior to the BR one, especially in absorption, where the BR results are comparable to the original DDA. Thus, we strongly recommend MG-DDA for coated BC aggregates, which agrees with previous results of bare BC aggregates [21].

The normalized phase functions ${P_{11}}$, computed with the MG and original DDA using dpl = 4 and 8, are shown in Fig. 8 in comparison with the MSTM. For dpl = 4 (red lines) the MG-DDA has good accuracy (and significant improvement with respect to the original DDA) for the forward direction (up to 30°), which is expected from the results for integral scattering quantities (Fig. 7). MG-DDA is also more accurate than the original one for scattering angles above 60°, but that is hardly practically relevant since the errors are still huge. By contrast, for dpl = 8 (green lines) the MG-DDA is accurate (almost within the line width) up to 160° and much better than the original DDA. The only problematic region is that of (near-)backscattering. On the one hand, it is known that accurate backscattering calculation requires larger dpl and number of orientations in the DDA [54,55] (and larger computational resources in many other methods), so it is not surprising that this applies to EMA-DDA as well. On the other hand, the behavior of phase function at backscattering is an interference phenomenon, related to the coherent backscattering for discrete random medium [56,57]. This interference is sensitive to the detailed structure of the internal electric fields, which is not necessarily made more accurate by a simple EMA (as discussed in Section 3).

 

Fig. 8. Comparison of the normalized phase function, computed with the DDA for two dpl values (4 and 8 – red and green, respectively) with and without the MG EMA (solid or dashed lines), to the reference MSTM results (black dashed line).

Download Full Size | PDF

Although the previous results provide some hints on the potential acceleration of computations by using the EMA, it is meaningful to compare the computational time only for the same accuracy, which is not-trivial when comparing different methods or variants of the same method [41,58,59]. To address this ambiguity, we present part of the results of Fig. 7 as the computational time versus the relative error of ${Q_{\textrm{ext}}}$ in Fig. 9 for the original DDA and that with the MG. The dpl is the implicit parameter in this plot, which scales roughly as cubic root of the computational time. The exception to this scaling and the difference in computational time between the original and MG-DDA result for the same dpl is due to the adaptive orientation averaging in ADDA, which in some cases uses up to 3 times smaller number of independent evaluations of the internal fields than in other cases. Overall, this plot clearly shows the steep increase of the simulation time required to improve the accuracy – a common feature of the DDA and other discretization methods. The MG-DDA is more efficient for the vast majority of error values, but that is hardly meaningful when the errors are larger than 10%. In the range of the MG-DDA errors from 1 to 10% original DDA may seem to be superior but that is related to change of error sign which happens at these values of dpl. Moreover, the convergence is inherently oscillating due to the random interaction between a dipole grid (determined by dpl value) and real particle geometry. The most important factor is probably the position of the coating boundary with respect to the grid (cf. Figure 1).

 

Fig. 9. The computational time of the DDA (original and with MG) versus the absolute value of the relative error of ${Q_{\textrm{ext}}}$ (with regards to the MSTM results) in log-log scale. The simulation data is the same as in Fig. 7.

Download Full Size | PDF

Thus, it is more reasonable to consider not an error for a single simulation, but rather the maximum error for this and all more computationally intense simulations (see, e.g., [58]). In particular, let us consider the typical requirement of 1% relative accuracy illustrated by the vertical line in Fig. 9. The MG-DDA gets into these bounds at about 20 s computational time and further stays there, while the original DDA does the same at only 600 s (30 times slower). Similar improvement is obtained, if 2% or 3% accuracy is required. For smaller errors the MG EMA still provides significant improvement (e.g., about 10 times acceleration for a threshold of 0.5% error), but further numerical studies are required to substantiate this claim. These results support the observations above that the EMA is most useful and practically relevant for BC aggregates when ${n_D}$ (for a monomer) is from 1 to 4. We do not show the corresponding plots for ${Q_{\textrm{sca}}}$ and ${Q_{\textrm{abs}}}$ since they are very similar (as can be deduced from Fig. 7). The only notable difference is that the error of the original DDA changes sign at different values of dpl (and, hence, computational time) highlighting the random nature of this crossing.

5. Conclusion

We investigated the smoothing of the DDA discretization of a particle using two EMA theories (MG and BR) to compute the effective refractive index for boundary dipoles. Using the Lorenz–Mie theory and MSTM as references for a number of test cases (spheres and coated BC particle), we contrasted it with the traditional DDA, where the representations of the particle shapes can be improved (to reduce simulation errors) only by increasing the number of dipoles.

For spheres, we varied both x and dpl (in addition to $m$), but the most relevant parameter with regards to the EMA performance is ${n_D}$, since it describes how fine the cubical discretization reproduces the spherical shape. In particular, for the refractive indices of 1.6 + 0.0i and 1.2 + 0.6i, both MG and BR significantly improve the DDA accuracy of integral scattering quantities (up to 60 times) when ${n_D} \,< \,4$ and 10, respectively. This necessarily implies that x should also be sufficiently small ($x$ < 3), since otherwise very small dpl leads to poor simulation accuracy (even with the EMA). The same conclusion holds for MG-DDA in a wide range of m (for $x$ = 2). In these cases, the application of the EMA is comparable to doubling the value of the dpl in the original DDA. Therefore, good accuracy (a few percent errors in ${Q_{\textrm{ext}}}$) can be obtained faster by about an order of magnitude.

Fortunately, such small size parameters are typical for BC monomers at ultraviolet, visible, and shortwave infrared wavelengths. Thus, we further tested the EMA smoothing for a coated BC model, based on transmission electron microscope observations. The coating diameter is 8 times larger than that of a BC monomer (which is 1/5 of the wavelength), thus the shape errors for the coating boundary are relatively small. For this coated aggregate, the MG is more accurate than the BR, and is most useful when ${n_D}$ (for a monomer) is from 1 to 4. In particular, for ${n_D}$ = 1 the relative error of ${Q_{\textrm{ext}}}$ is reduced from 4.7% to 0.3%. Moreover, the MG-DDA needs 3 times smaller dpl to achieve the accuracy of 1% for extinction, absorption, and scattering efficiencies than the original DDA, corresponding to about 30 times faster simulations. The EMA also improves the accuracy of the phase function, except for near-backscattering direction.

Overall, the combination of the DDA + EMA is expected to largely improve the simulation accuracy at the same dpl or to save computational time for the same accuracy for a wide class of particles that are similar to coated aggregates, i.e., combine simply-shaped parts larger than wavelength with large number of particles much smaller than the wavelength. This ensures that using a few dipoles per a smaller constituent (where the largest efficiency of the EMA is expected) results in sufficiently fine discretization of a larger one. However, any practical application will necessarily require preliminary tests to tune the method for the specific particle, wavelength, and refractive indices. This includes choosing the most efficient EMA formulation and optimal discretization level (for the required accuracy).

Despite the wide potential applicability, the presented approach does have its limitations. First, it may be an overkill for small homogeneous particles with simple shapes, where the volume correction employed in the standard DDA codes has almost the same effect as the EMA smoothing. Second, the EMA is a relatively crude approximation, which is not beneficial for sufficiently large x and m of smaller constituents, as well as when larger dpl is required to achieve the desired accuracy. In this case, weighted discretization [38] seems to be a promising approach. Moreover, we plan to study the combination of the EMA-DDA with modern DDA formulations, such as the integration of Green’s tensor [60] and filtered coupled dipoles [61], as well as Richardson extrapolation [21,62].