Open Access
Add to BibSonomy
Abstract
The discrete dipole approximation (DDA) is capable of treating scatterers with arbitrary shape and composition. However, large values of the refractive index require additional considerations. DDA calculations are performed for small spheres with a stronger absorbing inclusion and compared to T-matrix results. A natural phenomenon with strong optical contrasts is the melting of ice hydrometeors at microwave frequencies, and hence corresponding refractive indices are chosen. The obtained extinction and absorption efficiencies are found to depend mainly on the dipole size, whereas the phase function closely follows the T-matrix results by choosing a smaller stopping criterion and changing the polarizability formulation.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
Introduction. The discrete dipole approximation (DDA) possesses a broad user base thanks to its ability to handle scatterers of arbitrary shape and composition [1]. Traditionally, the DDA has been considered to be limited to cases in which the scatterer’s complex refractive index satisfies $|m-1|<2$ [2]; however, it has been shown that the DDA is capable of treating problems with extreme values of the refractive index [3].
One example of scatterers consisting of two materials with a high optical contrast is the melting of ice particles in the microwave range. Current and future satellite missions carrying microwave instruments, as well as providing guidance for radar applications, give rise to obtaining the single scattering properties of ice hydrometeors covering the frequency range from 1 to 886 GHz, as done in [4]. For this entire spectral range, the complex refractive index of ice [5] satisfies $|m-1|<2$, whereas liquid water only satisfies this conditions for frequencies above $\sim 150$ GHz [6].
The DDA has previously been employed for studying melting ice particles with a high optical contrast (see, for example, [7–11]). Based on insights from prior studies, especially those gained from Refs. [3,11], we attempt to provide some guidance when considering scatterers with high optical contrast, such as melting ice particles in the microwave region. Specifically, we aim to assess what dipole sizes are necessary to obtain accurate results using the DDA at these frequencies and what computational burden would be incurred in doing so with random orientation. We represent the melting ice particles as concentric spheres with radii smaller than the wavelengths.
Single scattering calculations. Single scattering calculations were performed at three frequencies, 10 GHz, 88.8 GHz, and 166.9 GHz. The corresponding wavelength in mm and the respective refractive indices of ice and liquid water are shown in Table 1, assuming a temperature of $0^{\circ }$C. While for ice the refractive index fulfills the condition of $|m-1|<2$, which should be considered as a practical guideline [3], at all three frequencies, for liquid water, the condition holds only at 166.9 GHz.
Table 1. Frequencies, Corresponding Wavelengths, and Refractive Indices Used in This Study
Reference calculations were performed using the Multiple Spheres T-Matrix (MSTM, version 3.0) code [12]. While MSTM was constructed for collections of spherical scatterers, in each of our calculations, we apply it to a single concentric sphere. The code is based on the superposition T-matrix method. The method and hence the code itself are applicable for scatterers composed of spheres with non-intersecting surfaces, i.e., the surfaces may, at most, share a single point. The code solves the interaction equation iteratively. Orientation-averaged optical properties are computed from the T-matrix by using an analytical solution to the integration over Euler angles. MSTM directly calculates the scattering matrix, reducing post-processing. We further validated the scattering efficiencies using the Mie code BHCOAT [13].
In the case of the DDA, the scatterer is divided into small fully polarizable volume-elements or dipoles, which are much smaller than the wavelength. The interaction of the dipoles with each other and the incident field result in a set of linear equations. These equations can be solved using standard numerical methods. The DDA’s ability to treat scatterers with arbitrary shape and composition stems from dividing the scatterer into dipoles. Here the DDA code ADDA (version 1.3b4) was used [2].
Two sets of DDA calculations were performed: the first set combines the lattice dispersion relation (LDR) formulation for polarizability and a stopping criterion of the iterative solver of $\epsilon = 10^{-5}$. The LDR and this stopping criterion are ADDA default settings (see, e.g., [2]). These settings were used for pure ice particles in Ref. [4]. The second set combined a stopping criterion of $\epsilon =10^{-10}$ with using the filtered coupled dipoles (FCD) formulation. The FCD formulation is considered well suited for large values of refractive index. The latter choice of the stopping criterion and the polarizability formulation follows the results reported in Ref. [3] and was also used in Ref. [11]. The dipole size was changed by increasing the number of dipoles per wavelength (dpl) from 15 to 600 (in the case of $f=10$ GHz from 35 to 15,000) when executing ADDA. The input values of dpl were subject to changes, due to ADDA performing volume correction. The dipole size was then evaluated according to the resulting ADDA log files.
For the purpose of providing a lower estimate of the computational burden that would be incurred during calculations of orientation-averaged scattering by non-spherically symmetric scatterers, we perform orientation averaging, even though it does not affect the scattering results for our concentric spheres. In ADDA, particle orientation is specified by three Euler angles. Orientationally averaged results are calculated internally by averaging over the Euler angles. The averaging is performed by using Romberg integration, which leads to uniformly spaced angles. The number of orientations is controlled by the minimum and maximum number of subdivisions, $J_\mathrm {min}$ and $J_\mathrm {max}$, respectively. (The maximum number of angles can the be calculated by $2^{J_\mathrm {max}}+1$.) Orientation averaging is performed until the maximum number of subdivisions is reached or the error of $C_\mathrm {ext}$ is below a predefined stopping criterion for the orientation averaging $\mathrm {eps}$ [2]. Here we used orientation averaging settings of $J_\mathrm {max}=4$ and $\mathrm {eps}=10^{-6}$ for the first set (again in accordance with Ref. [4]) and $J_\mathrm {max}=5$ and $\mathrm {eps}=10^{-10}$ for the second set. For both sets, we chose $J_\mathrm {min}=2$.
In every case, the $F_{11}$ element of the scattering matrix was calculated for a scattering angle resolution of $\Delta \vartheta =0.25^{\circ }$. For the ADDA calculations, the enhanced, stabilized bi-conjugate gradient solver [Bi-CGStab(2) or bcgs2] was used.
MSTM calculations were performed on an eight-core desktop computer with Intel i7-4790 processors at 3.6 GHz. Due to their higher computational demand, ADDA calculations were performed on the cluster Vera at the Chalmers Centre for Computational Science and Engineering. The cluster Vera is equipped with 64-core (of which 32 cores are virtual cores due to hyperthreading) Intel Xeon Gold 6130 processors at 2.1 GHz.
The required dipole spacing for a scatterer is frequently considered to be sufficient if $|m|kd<0.5$ [4,10,11]. Here, $k$ is the wavenumber of the incident electro-magnetic wave and $d$ the length of an individual dipole. This criterion can be related to the number of dipoles per wavelength $\mathrm {dpl}$, which was used to control the dipole size, by $|m|kd=2|m|\pi /\mathrm {dpl}$ [14].
When considering the physically more realistic case of an outer layer of liquid water, small melting fractions are difficult to achieve using the DDA, as the dipoles have to be rather fine, as already pointed out in Ref. [11]. To overcome this limitation, we used an inverted geometry with a liquid water core within an ice hull. While this inverted geometry is not representative for melting ice particles, it allows us to test the effects of dipole sizes at different melting fractions more accurately. This is the main reason for limiting ourselves to concentric spheres, as more complex morphologies, even the only moderately more complex case of a dimer, i.e., two spheres in point-contact, would require a finer dipole spacing to accurately discretize the melting water. By using concentric spheres, we can perform reference calculations using publicly available code.
In this study, melting fractions were expressed as mass fractions, for which the density of ice and the density of liquid water at $0^{\circ }C$ were assumed to be $\rho _\mathrm {ice}=0.9167\,\mathrm {g}\,\mathrm {cm}^{-3}$ and $\rho _\mathrm {liq. water}=0.9998\,\mathrm {g}\,\mathrm {cm}^{-3}$ [15]. Here three different melting mass fractions $f_\mathrm {water}$ were considered: $f_\mathrm {water}=0.05$; $f_\mathrm {water}=0.1$; and $f_\mathrm {water}=0.3$. The radius of the corresponding pure ice spheres was set to $r=100\,\mathrm{\mu}$m. Assuming preservation of mass, adding the liquid water core requires volume correction of the sphere. This leads to a radius of $r=99.2\,\mathrm{\mu}$m for $f_\mathrm {water}=0.3$. This small change in radius can, in addition, increase the robustness of our results, since spheres are known to be difficult to treat with the DDA due to size-sensitive resonances [16]. For the ADDA calculations, the shapes are created by using the predefined shape “coated" (cf. [7]).
Results and discussion. The absolute percentage error between DDA results and MSTM results for the extinction [Figs. 1(a)–1(c)] and absorption efficiencies [Figs. 1(d)–1(f)] as function of dipole size are shown.
Fig. 1. Absolute percentage error between MSTM reference results and DDA calculations as a function of dipole size in (a)–(c) $Q_{ext}$ and (d)–(f) $Q_{abs}$ at (a),(d) 10 GHz; (b),(e) 88.8 GHz; and (c),(f) 166.9 GHz for three different melting mass fractions: 0.05 (violet); 0.1 (green); and 0.3 (magenta). The dotted gray line indicates 10% and the gray-shaded area indicates 5% deviation.
The accuracy of the phase function at different dipole sizes is gauged by using the mean absolute percentage error (MAPE), which can be calculated as
Here $A_t$, the actual value, refers to the MSTM result at a scattering angle with index $t$ and the forecast $F_t$ refers to the DDA result at index $t$ for the respective dipole size. Additionally, $n$ denotes the number of angles. The MAPE as a function of dipole size in micrometers is shown in Figs. 2(a)–2(c).
Fig. 2. (a)–(c) Mean absolute percentage error between the MSTM reference result of the phase function and the DDA calculations as a function of dipole size. (d)–(f) Discretization error of the melting mass fraction. Columns and colors are as in Fig. 1.
To disentangle the effects of the optical contrast and effects of discretization, we calculated the relative discretization error of the melting mass fraction:
Here, $N_\mathrm {water}$ refers to the number of dipoles associated with liquid water and $N_\mathrm {total}$ refers to the total number of dipoles. The prescribed melting mass fraction is denoted as $f_\mathrm {water, analytical}$. The discretization error is shown in Figs. 2(d)–2(f). As expected, the discretization error is reduced with decreasing dipole size, as it becomes easier to accurately represent the liquid water inclusion.
Figures 1 and 2 cover the first set of DDA calculations ($\epsilon =10^{-5}$ and LDR). The results for the second set are basically identical, with the exception of the MAPE at $f=10\,$GHz. Hence, we performed additional calculations at $f = 10$ GHz using $\epsilon = 10^{-10}$ and FCD. The resulting MAPE values, which are now well below 0.02%, are displayed in Fig. 3. For reference, the differences stemming from the different sets of DDA calculations are shown in Fig. S1.
Fig. 3. Mean absolute percentage error between the MSTM reference result of the phase function and the DDA calculations with $\epsilon =10^{-10}$ and using FCD as a function of dipole size at $f=10$ GHz.
For the smallest dipole sizes, i.e., $d\lesssim 10\,\mathrm{\mu}$m, the discretization error of the water inclusion is close to zero. In contrast, larger dipole sizes, resulting in large values of $\delta f_\mathrm {water}$, coincide with large values of $\delta Q_\mathrm {ext}$ and $\delta Q_\mathrm {abs}$ at $f=88.8$ GHz and $f=166.9$ GHz. This indicates discretization errors to impact the absolute percentage error of the efficiencies at these frequencies. In contrast, at $f=10$ GHz, both $\delta Q_\mathrm {ext}$ and $\delta Q_\mathrm {abs}$ show high values for dipole sizes between $10\mathrm {-}30\,\mathrm{\mu}$m, which coincide with small values of the discretization error, pointing toward the high optical contrast as an error source. For the MAPE of set 1 (see Fig. 2), we see the largest values for dipole sizes that are lower than the sizes at which the largest discretization error occurs. To investigate this further, we performed additional calculations at $f=10$ GHz and $f=88.8$ GHz with $\epsilon =10^{-6}$ and using LDR, but using the refractive index at $f=166.9$ GHz. In other words, we took the refractive indices for these additional calculations to be $m=1.7856+i0.0022$ for the ice shell and $m=2.5036+i1.0209$ for the liquid water inclusion, which fulfill the recommendation of $|m-1|<2$. The results of these calculations are shown in Fig. S2. The MAPE at $f=10$ GHz, obtained from the additional calculations, shows a similar behavior, i.e., it reaches its maximum, when the discretization error’s absolute value is reducing. Therefore, this pattern is consistent at both high and low optical contrasts, and hence, it is unlikely to be caused by the strong optical contrast. However, $\delta Q_\mathrm {ext}$ and $\delta Q_\mathrm {abs}$ at $f=10$ GHz show different behavior in the additional simulations, pointing toward effects of the high optical contrast as an error source. This means that both the discretization error of the water inclusion and the high optical contrast contribute to errors in the DDA calculations.
Figure 4 shows the CPU time for the ADDA calculations in hours as a function of dipole size for both sets of DDA calculations, as denoted by full and dotted lines. In general, the required computational burden increases with decreasing dipole size. This is entirely expected, as smaller dipoles increase the number of dipoles and hence the complexity of the calculations. The computational burden of DDA calculations is known to increase with increasing refractive index and with increasing size parameter [17]. For $f=10$ GHz, the impact of the high refractive index of liquid water outweighs the impact of the small size parameter. Following the recommendations for high refractive indices laid out in Ref. [3] (reducing the stopping criterion $\epsilon$ and using FCD) decreases the computational burden compared to the default settings ($\epsilon =10^{-5}$ and LDR). This reduction in computation time comes despite extending the maximum number subdivisions used in orientation averaging and decreasing the stopping criterion used within the orientation averaging scheme from $\mathrm {eps}=10^{-5}$ to $\mathrm {eps}=10^{-10}$ between set 1 and 2 of the DDA calculations. For moderate refractive indices at $f=88.8$ GHz and $f=166.9$ GHz, this no longer holds. Please note that the computation time might also be affected by the choice of iterative solver.
Fig. 4. CPU time in hours required by ADDA calculations as a function of dipole size. Colors are as in Fig. 2. Full lines indicate DDA calculations performed with $J_\mathrm {max}=4$, $\mathrm {eps}=10^{-6}$, $\epsilon =10^{-5}$, and using LDR. Calculations performed with $J_\mathrm {max}=5$, $\mathrm {eps}=10^{-10}$, $\epsilon =10^{-10}$, and using FCD are represented by dotted lines.
The comparison of DDA and T-matrix results for small spheres with high optical contrast allows us to arrive at the following conclusions.
- • The dipole size is a key parameter controlling the accuracy of the DDA results. The common guideline of $|m|kd<0.5$ is insufficient when treating particles with strong optical contrast, as already discussed in Refs. [3,11]. This criterion is already fulfilled for dipoles with $d\leq 625\,\mathrm{\mu}$m at $f=10$ GHz, $d\leq 160\,\mathrm{\mu}$m at $f=88.8$ GHz and $d\leq 105\,\mathrm{\mu}$m at $f=166.9$ GHz. We show that uncertainties in $Q_\mathrm {ext}$ and $Q_\mathrm {abs}$ below $10\%$ for $f=10$ GHz and below $5\%$ for $f=88.8$ GHz and $f=166.9$ GHz, require dipole sizes of $d\leq 2.5\,\mathrm{\mu}$m, $d\leq 8\,\mathrm{\mu}$m, $d\leq 10\,\mathrm{\mu}$m, respectively. An additional source of uncertainty is insufficient representation of the melting mass fraction on the dipole grid. Thus reducing the dipole size reduces two different errors affecting the optical properties.
- • When investigating extinction and absorption efficiencies alone, the ADDA default stopping criterion ($\epsilon =10^{-5}$) in conjunction with the LDR are sufficient to reduce the uncertainties solely by decreasing the dipole size. However, when also modeling the phase function (i.e., the scattering matrix element $F_{11}$), decreasing the stopping criterion in accordance with Ref. [3] to $\epsilon =10^{-10}$, and employing the FCD formulation drastically reduces the MAPE of the DDA results compared to the MSTM results.
- • Considering melting ice as an application of this study, the dipole size required for low uncertainties ($<10\%$ or even $<5\%$) is far below the dipole size used for calculating the single scattering properties of the ice particles included in the database presented in Ref. [4]. That is, an extension of the database to include melting particles would require drastically smaller dipole sizes applied and thus long calculation times are expected.
- • While exceeding the traditionally recommended applicability guidance of $|m-1|<2$, we see a lower computational demand for calculating the optical properties with $\epsilon =10^{-10}$ using the FCD, an extended number of orientations at $f=10$ GHz compared to standard ADDA settings (LDR and $\epsilon =10^{-5}$), and fewer orientations considered. This lower computational demand can be attributed to a faster calculation of the internal fields. At $f=88.8$ GHz, at which the refractive indices still exceed the applicability guidance, using $\epsilon =10^{-10}$, FCD, and an increased number of orientations increase the computational burden without discernible benefits in accuracy. However, our reported CPU times should only be considered as a lower estimate. When considering more complex shapes, such as the melting snowflakes presented in Refs. [9–11,18], the computational burden is expected to increase.
Funding
Swedish National Space Agency (166/18).
Acknowledgments
We are grateful to Maxim Yurkin and Alfons Hoekstra for making their ADDA code publicly available and to Dan Mackowski for making his MSTM code publicly available . The calculations were partly performed by resources provided by the Swedish National Infrastructure for Computing (SNIC) at Chalmers Centre for Computational Science and Engineering (C3SE).
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are available in Dataset 1, Ref. [19].
Supplemental document
See Supplement 1 for supporting content.
REFERENCES
1. M. A. Yurkin, D. A. Smunev, S. A. Glukhova, A. A. Kichigin, A. E. Moskalensky, and K. G. Inzhevatkin, in Radiation and Scattering of Electromagnetic Waves (RSEMW), (IEEE, 2021), pp. 111.
2. M. A. Yurkin and A. G. Hoekstra, J. Quant. Spectrosc. Radiat. Transfer 112, 2234 (2011). [CrossRef]
3. M. A. Yurkin, M. Min, and A. G. Hoekstra, Phys. Rev. E 82, 036703 (2010). [CrossRef]
4. P. Eriksson, R. Ekelund, J. Mendrok, M. Brath, O. Lemke, and S. A. Buehler, Earth Syst. Sci. Data 10, 1301 (2018). [CrossRef]
5. C. Mätzler, W. Ellison, B. Thomas, A. Sihvola, and M. Schwank, “Di electric properties of natural media,” (2006).
6. W. J. Ellison, J. Phys. Chem. Ref. Data 36, 1 (2007). [CrossRef]
7. J. Tyynelä, T. Nousiainen, S. Göke, and K. Muinonen, J. Quant. Spectrosc. Radiat. Transfer 110, 1654 (2009). XI Conference on Electromagnetic and Light Scattering by Non-Spherical Particles: 2008. [CrossRef]
8. J. Tyynelä, J. Leinonen, D. Moisseev, T. Nousiainen, and A. von Lerber, J. Quant. Spectrosc. Radiat. Transfer 133, 504 (2014). [CrossRef]
9. D. Ori, T. Maestri, R. Rizzi, D. Cimini, M. Montopoli, and F. S. Marzano, J. Geophys. Res.: Atmos. 119, 9931 (2014). [CrossRef]
10. B. T. Johnson, W. S. Olson, and G. Skofronick-Jackson, Atmos Meas. Tech. 9, 9 (2016). [CrossRef]
11. D. Ori and S. Kneifel, J. Quant. Spectrosc. Radiat. Transfer 217, 396 (2018). [CrossRef]
12. D. Mackowski and M. Mishchenko, J. Quant. Spectrosc. Radiat. Transfer 112, 2182 (2011). [CrossRef]
13. C. F. Bohren and D. R. Huffman, Absorption and scattering of light by small particles (Wiley, New York, 1983).
14. M. A. Yurkin, V. P. Maltsev, and A. G. Hoekstra, J. Opt. Soc. Am. A 23, 2592 (2006). [CrossRef]
15. W. M. Haynes, D. R. Lide, and T. J. Bruno, CRC Handbook of Chemistry and Physics (CRC Press, 2016).
16. Y. Zhu, C. Liu, and M. A. Yurkin, Opt. Express 27, 22827 (2019). [CrossRef]
17. M. Yurkin and A. Hoekstra, “User manual for the discrete dipole approximation code adda 1.3b4,” https://github.com/adda-team/adda/releases/tag/v1.3b4 (2014).
18. J. Leinonen and A. von Lerber, J. Geophys. Res.: Atmos. 123, 1811 (2018). [CrossRef]
19. F. Kanngießer and P. Eriksson, “Optical properties of small scatterers with high optical contrast,” zenodo, 6726136, (2022), https://zenodo.org/record/6500523#.Yntfz1zP3Jk.
