1. Introduction

The phenomenon of resonance in dielectric spheres was first predicted by Debye [1], and describes the phenomena of optical-and-particle interaction with certain modes of electric and magnetic vibrations. Ashkin and Dziedzic [2] observed sharp resonances of a small suspended particle in a given position using optical levitation. The resonances attract more attention after Chýlek et al. [3] performed accurate simulations and explanations on the levitation experiment. These resonances are now generally referred to as morphology-dependent resonances (MDR). The MDR of spheres can be explained analytically by the theory of quantum-mechanical shape resonances [4]. A physical interpretation is that the electromagnetic wave is trapped by total internal reflection as it propagates around the inside surface of the sphere, and, after circumnavigating the sphere, the wave returns to its starting point in phase [4].

An easy way to identify resonance on scattering properties is the occurrence of the sharp peaks on scattering property curves as a function of particle size parameter x = 2πa/λ [5–11], where a and λ are the sphere radius and the incident wavelength, respectively. Further, we define the full width of a MDR at half maximum as Δx (as shown in Fig. 1), which specifies the (very small) difference of x for two spheres to have significantly different scattering properties. The relative width of a MDR is defined as Δx/x. Chýlek et al. [3] classified these resonances into three kinds based on the relative width of the MDR. The first-order resonance has relative width in the order of 10⁻⁷~10⁻⁸, which is extremely narrow and is hard to detect in experiments [3]. The second- and third-order resonances have relative width in the order of 10⁻⁵ and 10⁻², respectively. Lam et al. [12] gives simple approximated formulae showing that resonance width rapidly decreases with the increases of x and/or the refractive index m.

 

Fig. 1 Extinction efficiency (Q_ext) and asymmetry factor (g) of spheres with m = 1.6 given by the Mie simulations with different range of x.

Download Full Size | PDF

The MDRs of spheres or spheroids can be numerically obtained and exhibited using the Mie theory [2–12] or the T-matrix method [13,14], which both solve the Maxwell’s equations semi-analytically. However, this may be not true for discretization-based numerical methods, or at least the MDRs can be significantly distorted. For example, Liu et al. [15] evaluated the discrete dipole approximation (DDA) and pseudo-spectral time domain (PSTD) method for light scattering simulations by comparing with the Mie results, and obtained excellent agreements (relative errors on the extinction cross sections less than 1%) for most considered spheres. There was, however, a single outlier with x = 20 and m = 1.6 (relative errors over 5%, see Table 1 in [15]), which was neither specially mentioned nor explained. The retrospective analysis of that artifact motivated current study. In short, this case was a coincidental hit on a MDR, while both the DDA and PSTD missed this MDR producing the scattering quantities close to the baseline (see Section 2 for details). A positive hint was given by [16], where the DDA was shown to qualitatively reproduce the wider (third-order) MDRs (as a function of m), but at a shifted position. By contrast, the study of high sensitivity of MDRs to slight variation of the particle shape [13] suggested that reproducing narrow MDRs with discretization-based methods, including the DDA, is “essentially impossible”.

Table 1. The absolute errors of the extrapolated values of both peak and baseline parameters for the MDR with m = 1.2 and x = 50.5369 (see Figs. 5 and 6) using three different grid ranges of DDA simulations in comparison with the Mie result. The values for n = 768 are from the direct fit of the DDA results.

View Table | View all tables in this article

To avoid misleading results related to the MDR and to better understand the DDA capabilities, this study investigates whether the DDA can reproduce the MDRs of spherical particles, and how fine the dipoles should be set to obtain the MDRs. This study is organized as follows. The MDRs from the Mie simulation is discussed in Section 2, and Section 3 presents the performance of the DDA in reproducing a MDR. In the latter we start with tedious brute-force simulations for very fine discretization levels, but then employ Lorentzian fits and extrapolation technique to further improve the accuracy. Section 4 is devoted to another MDR, for which the DDA simulations result in peak splitting. Section 5 gives the guidelines for simulating MDRs with the DDA, while Section 6 concludes this study.

2. The morphology-dependent resonances

The MDRs of spheres can be numerically presented through Mie simulations with fine x resolution, or one may stumble upon a MDR without expecting it. Multiple studies have been performed to study the influence of MDR on scattering properties and here we only present a simple introduction. Figure 1 shows the extinction efficiency and asymmetry factors of spheres with an m = 1.6. Here and further we used the Mie code by Bohren and Huffman [17]. A step of 10⁻⁴ is used for x, and clear spikes appear. The peaks appear quasi randomly with respect to the main oscillations and have relative amplitudes of ~10% of the baseline value. Moreover, they have the Lorentzian shape [17,18]. The upper right panel details a peak with x around 20, i.e., the specific case Liu et al. [15] encountered. More specifically, Q_ext peaks at 2.725 when x = 20.00010 with width Δx = 2.4 × 10⁻⁴, which classifies it as a second-order MDR. At exactly x = 20 we have Q_ext = 2.608, which is almost half-way to the maximum (different by 6% from the baseline value of 2.450). We illustrate the asymmetry factor g in the bottom panels. Obviously, the spikes on g occur at almost the same x as those for the Q_ext but are directed downwards.

The peaks on the scattering properties caused by the MDRs can be mathematically explained by the Mie equations. The Mie theory gives the extinction of a sphere by [19,20]:

Besides the integral scattering properties, the MDRs also significantly affect the scattering matrix elements [16,21,22]. Figure 2 compares the non-zero scattering matrix elements of three spheres at and around the MDR with x = 20.0001, 19.9981, and 20.0021, respectively, i.e., the green star for the MDR and the red and blue ones for the baseline cases marked in Fig. 1. The elements at its two sides are almost identical, whereas the one between them at the peak position shows more significant oscillations compared to its two neighbors. In particular, the phase functions with scattering angles larger than 40° differ significantly, and a stronger backward scattering is noticed for the MDR case. This explains the smaller asymmetry factor at the MDR.

 

Fig. 2 Comparison of the scattering matrix elements of spheres with m = 1.6 and x around 20. The green lines characterize the phase matrix elements of the sphere with x = 20.0001, and the red and blue lines characterize the ones with x = 19.9981 and 20.0021, respectively.

Download Full Size | PDF

Now we are ready to explain the results of volume-discretization methods for this case (x = 20), i.e. the outliers in [15]. Both the DDA and PSTD discretize the scatterer using dipoles or voxels, and a few tens of voxels within a wavelength are normally considered. Up to 40 dipoles per wavelength were used in the DDA simulations of [15] (and up to 20 in the PSTD), corresponding to a dipole size parameter kd = 0.16 (d is dipole size, k is the wave number), which is much larger than the width of this (and any other) MDR. This necessarily implies that in both DDA and PSTD simulations the resonance peak is significantly distorted (shifted and widened, see Section 3). We do not reproduce this DDA peak in details due to large computational requirements, but the PSTD and DDA values for x = 20 are 2.460 and 2.467, respectively. They are within 0.7% of each other and of the baseline Mie value (2.450), but more than 5% different from the true Q_ext at x = 20. The latter difference could even be twice larger if we were lucky to hit the top of the peak instead its side lobe. The conclusion of the above analysis is that MDRs should be explicitly avoided when using the spheres for benchmarking other codes. Fortunately, that can be easily done by calculating Q_ext(x) in a small range of x around the chosen value.

To facilitate the choice of MDRs for further analysis, we show in Fig. 3 the extinction efficiency of spheres with different values of m. Obviously, the MDRs occur at smaller x for spheres with larger m, and become narrower as x increases (see Fig. 3). There are no MDRs for spheres with m = 1.2 and x< 30. This guides us the choice of MDR cases for the DDA simulations. First, it is easier to consider spheres with smaller x with respect to the required computational resources. Second, the computing time of the DDA simulations for same-sized spheres also increases with m due to the slower convergence of the iterative solver [15]. As a result, in the next section we consider a MDR with relatively small Δx (to make it non-trivial) but with a relatively small x or m to make it feasible for the DDA.

 

Fig. 3 Extinction efficiency calculated by the Mie theory versus x with m varying from 1.2 to 1.6 with a step of 0.2.

Download Full Size | PDF

3. DDA simulations for the MDR

The DDA is a popular and robust method for numerical simulation of light scattering by dielectric particles [23–27]. The dipoles (discretization voxels) should be much smaller than both the incident wavelength and any characteristic dimension of the particle. For a sphere the evident scale is its radius a or diameter. But close to the MDR, some other scale should appear, corresponding to the narrow width of the MDR. Thus, we can expect that the dipole size parameter kd should be comparable to (or even smaller than) Δx for the DDA to accurately reproduce the MDR. In the following we use the number of dipoles along the sphere diameter, named as ‘grid number’ n, to specify the degree of the discretization. Fixing n is much more convenient for simulation of resonances, than the standard number of dipoles per lambda (“dpl”), since we do not need to worry about potential changes of the discretization geometry (set of cubes approximating a sphere) when varying x in a narrow range. This recipe is used in all the DDA simulations below.

We used code ADDA v. 1.3b4 [25], which can parallelize a single DDA simulation on a compute cluster, allowing very large n (we used up to 1024). We have used the default simulation parameters throughout the simulations. Studying more advanced formulations implemented in ADDA (e.g. filtered coupled dipoles [28] and integration of Green’s tensor [29]) for such simulations is an interesting research direction, but we haven’t pursued it in this paper. One issue associated with the default settings is the volume correction – ADDA automatically adjusts the size of the dipole, so that the volume of a discretized shape exactly equals that of the original sphere [25]. In other words, the DDA simulation with a given d and n corresponds to two spheres with slightly different x, corresponding to the diameter (x_D = nkd/2) and volume of the discretized shape (x_V = kd[3N_d/(4π)]^1/3), where N_d is the total number of dipoles. Fortunately, x_V/x_D = f(n), i.e. a coefficient depending only on n. It is close to 1, but the difference may be comparable to Δx/x, and thus be not negligible. Importantly, the results of simulations using the volume correction (i.e. based on x_V) can be trivially transformed into that based on x_D, and vice versa. The same transformation can also be directly applied to processed quantities, such as the peak position. Our simulations showed that, the results based on x_D converge more smoothly to the Mie reference. Therefore, we use it for data presentation throughout the paper, corresponding to the default setting of ADDA.

First, we consider an easier case of $b_{55}^{(1)}$ resonance for m = 1.2, positioned at x = 50.5369 with the width as large as Δx = 0.096. Thus, the corresponding relative width is ~0.002, which is slightly narrower than the typical values for the third-order resonances. For this small refractive index, large x is required to have a MDR width of ~0.1, while larger Δx are not that interesting for the subject of this paper. Figure 4 illustrates the DDA results for spheres with x between 50.4 and 50.7 and m = 1.2 in comparison with the reference Mie results. Simulations are performed in sets of fixed n (from 128 to 768) and are a pinnacle of brute-force approach – total simulation time is about 9 core-years on a compute cluster. The spacing of x values varies between the sets since it evolves during the development of the data processing algorithms. In the end, we decided to use all the data we have.

 

Fig. 4 DDA simulations of the extinction efficiency versus x, corresponding to $b_{55}^{(1)}$ MDR (with m = 1.2). Markers indicate simulated data (with n from 128 to 768) and the solid lines show the Lorentzian fits with a cubic baseline (see text). Mie results are also shown together with a fit; the corresponding two lines overlap on the whole range.

Download Full Size | PDF

The immediate conclusion of these brute-force simulations is that the DDA is capable of directly reproducing this MDR, converging to the exact solution with increasing n, but only at a great computational cost. This agrees with the qualitative discussion in the beginning of this section. The simulated peak is significantly different for smaller n, including that for n = 192, which corresponds to the rule of thumb for using 10 dipoles per wavelength in the medium [25], i.e. kd = π/5m. By contrast, the most accurate result (n = 768) corresponds to kd = 0.13 – close to the MDR width Δx. The errors for smaller n become even more pronounced if one considers the simulated Q_ext at a specific x. For instance, when n = 128 the errors of Q_ext values at x = 50.605 and 50.537 (maxima of this DDA simulation and of the Mie reference) are both over 0.03 (i.e., relative error over 1.6%). That is a lot, since the whole peak amplitude (compared to the baseline) is about 0.08. However, convergence of Q_ext for fixed x (although being a standard approach [30]) may be misleadingly bad. Instead, Fig. 4 suggests tracking the evolution of the whole peaks (curves).

To quantify the evolution of the peaks, we parameterize them by fitting with a Lorentzian with a cubic baseline:

The Lorentzian fits are almost perfect, which is well known for the Mie theory [17,18], but was not a priori evident for the DDA. Still, that is not completely unexpected since the DDA result for a fixed grid mesh is the exact solution for the corresponding set of point dipoles. The latter also has resonances (poles in the complex plane of x), leading to the same Lorentzian peak shape. In particular, x_c–iw/2 is approximately the location of this pole (see Section 2). Hence, the convergence of x_c and w with n directly corresponds to the convergence of the pole position in complex plane of x. Moreover, these poles may be related to the spectrum of the DDA interaction matrix [28,31], but we leave it for future research. The only visible exception is the result for n = 128, which is effected by peak splitting, discussed in Section 4.

By performing those fits we compressed a lot of simulated data (from 50 to 300 points for each n) into several parameters. In the following we study the convergence of these parameters with increasing n. A natural parameter of discretization fineness is kd [30,32], but here it slightly varies over the data set with fixed n. Hence, we use 1/n = kd/(2x_D), almost linearly proportional to kd, that was previously used in the case of very small particles (when kd is irrelevant) [28]. Figures 5 and 6 illustrate the convergence of peak and baseline parameters, respectively, to the Mie results, where the error bars correspond to standard deviations of parameters obtained from the non-linear fits. Those errors are significant for smaller n, which correlates to the quality of fit in Fig. 4. Overall, the DDA results definitely converge to the true value, but significant oscillations are present.

 

Fig. 5 Convergence of the peak parameters (x_c, w, and A) from the DDA simulations, shown in Fig. 4 (for the MDR with m = 1.2 and x = 50.5369), with refining discretization. Solid lines correspond to the fits of c₀ + c₂n⁻² + c₃n⁻³ in three different ranges of n, extrapolated to 1/n→ 0 (see text for details). The insets show the results for 6 largest grids.

Download Full Size | PDF

 

Fig. 6 Same as Fig. 5 but for the baseline parameters (y₀, y₁, y₂, and y₃).

Download Full Size | PDF

Next, we employ the Richardson-style extrapolation, which was successfully used previously [28,32] based on the approximation of the DDA errors by a quadratic function of kd or 1/n [30]. In particular, we select a range of n (specified below) and perform a polynomial fit through corresponding data points, minimizing the sum of residual squares. Standard prescription consists in using quadratic function with weighting of residual squares by n⁶, i.e. assuming that expected fit residuals are proportional to n⁻³ – the first omitted term in the fitted function [32]. We will use this prescription in Section 4. However, for the data in this section we noticed that the linear in 1/n term is negligible. Thus, we use the cubic function without the linear term (c₀ + c₂n⁻² + c₃n⁻³) and weigh residual squares by n⁸ (expecting the residuals to scale as n⁻⁴). Overall, three is the optimal number of terms in approximating function given the amount of data and oscillations around the polynomial trend. The reason for disappearance of the linear term is not clear, but it is related to the large size of the sphere, since the same effect was seen for x = 15 sphere in [33]. Anyway, the main argument for using the modified approximating function is its empirical success in accurate reproduction of the Mie result, as discussed below.

Importantly, this is the first time the extrapolation is applied to the processed parameters, which have considerable errors apart from the standard DDA simulation errors. In particular, the error bars in the figures are determined by the number and spacing of x values in the original DDA data. But using the error weighting described above we ignore these processing errors. This may cause issues both for the smaller and larger n. For the former, the processing errors are definitely large enough to influence extrapolation and, although they increase with 1/n, the particular scaling may be different from the used exponent in error weighting. For the larger n the processing errors only weakly depend on n, which definitely contradicts the used residual weighting. Thus, if these errors cannot be neglected in comparison with the DDA simulation errors, our procedure will assign too large weights to these data points. We leave the detailed analysis of these issues for future research, as well as attempts to rigorously account for processing errors, combining them with n⁻³ or n⁻⁴ scaling.

The results of fits with various ranges of n are shown in Figs. 5 and 6. The value of a quadratic function at 1/n→ 0, i.e. its coefficient c₀, is the estimated (extrapolated) value of the corresponding parameter. Their errors (in comparison with the Mie values) are given in Table 1. In principle, the extrapolation also provides the standard error of the estimated value (internal error estimate) as a part of the standard polynomial fit [30]. In most cases they are larger than the actual errors (data not shown). However, it may be not adequate, when applied to our data sets with small number of points and large oscillations around the polynomial trend [30], and due to the unaccounted effect of processing errors. Thus, we do not further discuss these standard errors.

We consider three extrapolation ranges of n: the full range [128,768], the range [256,768] aiming at the best accuracy, and the range [128,320] aiming at satisfactory accuracy at much smaller computational cost. The ratios of total computational times for these ranges and for a single set of n = 768 are approximately 11:11:1:6, respectively (if using the same x values). The extrapolation works perfectly for x_c (see the inset in Fig. 5(a) and Table 1): in the range [256,768] it decreases the resulting error 42 times compared to already small error for n = 768 (from 2.8 × 10⁻³ to 6.6 × 10⁻⁵ in Table 1). In the ranges [128,768] and [128,320] the resulting errors are (8.0 × 10⁻⁴ and 7.4 × 10⁻⁴ in Table 1), respectively, 3 and 4 times smaller than that for n = 768. The latter is especially surprising, considering the 6 times smaller computational cost. Extrapolation in the range [256,768] is also perfect for w and A, decreasing the error 21 and 53 times, respectively, compared to that for n = 768 (from 2.6 × 10⁻⁴ and −4.1 × 10⁻⁵ to −1.2 × 10⁻⁵ and 7.8 × 10⁻⁷), see the insets in Figs. 5(b) and 5(c). However, in the full range the extrapolation error for w and A is slightly larger than that for n = 768, but describes the overall trend qualitatively well. When using the range [128,320] the error is significantly larger, comparable to that for n = 320, so in this range the extrapolation does not add value apart from qualitative description. The DDA simulation with n = 320 already results in 6% accuracy of w and A, which may be considered satisfactory, considering the relatively weak effect of these parameters on Q_ext at any specific x.

The extrapolation quality for baseline parameters (Fig. 6 and Table 1) is similar to that for w and A. The results of using the ranges [128,768] and [256,768] are generally comparable, the former leads to better accuracy for y₀ and y₁, while the latter – for y₂ and y₃. The best extrapolation error for y₀ and y₃ is, respectively, 5 and 2.3 times smaller than that for n = 768, see the insets in Figs. 6(a) and 6(d), but no improvement is observed for y₁ and y₂ (the errors are almost the same as that for n = 768). Still, the approximating function correctly captures the qualitative trends. Using the range [128,320] the extrapolation can be considered satisfactory only for y₀ and y₂, but that is not very important since the baseline is already very accurate for n = 320 (relative error of y₀ is 0.08%).

In summary, the proposed extrapolation procedure works remarkably well. Starting from the large-grid DDA simulations it reaches unprecedented accuracy, especially for peak parameters (see insets in Fig. 5). We stress that the Mie result is not anyhow used in the extrapolation procedure; thus, hitting it with such precision at large distance from the data points cannot be explained by pure coincidence. This accuracy is the main argument for the validity of this empirical technique. The extrapolation also allows accurately estimating x_c and y₀ based on the grid range [128,320], which potentially opens the way to simulation of narrower MDRs. At the same time, the choice of this coarser range remains empirical warranting further study before systematic use of this technique.

To finalize this section, we take a quick look at the angle-resolved scattering functions. Figure 7 shows the P₁₁ and P₁₂ elements of the MDR from the Mie theory and the DDA simulations with n = 192 and 768. The other two non-zero elements show the same behavior. The DDA result with n = 768 agrees well with the exact solution. For the default discretization case n = 192, there are clear differences from the Mie results, and the stronger backscattering is not obtained. Actually the coarse-grid DDA result looks similar to the off-resonance values, which is not shown here. It should be possible to apply the extrapolation procedure to scattering matrix as well, but we leave it for future research.

 

Fig. 7 Comparison of the P₁₁ and P₁₂ elements at the MDR with m = 1.2 and x = 50.5369. Shown are the Mie results and DDA simulations with n = 768 and 192.

Download Full Size | PDF

4. Peak splitting in the DDA simulations

A somewhat narrower resonance is considered in the following. It is $b_{27}^{(1)}$ for m = 1.6, positioned at x = 13.24417 with width Δx = 0.012. Thus, the relative width is ~9 × 10⁻⁴. We performed the DDA simulations for 11 values of n between 128 and 1024. In contrast to the previous MDR, this one is well isolated from the neighboring ones; hence, we performed simulations in a range of x from 13.15 to 13.35 – 20 times wider than its width to get a clear view of the baseline. The results are shown in Fig. 8, but only for three values of n (128,384, and 1024) and for the Mie theory to avoid clutter. The results for all simulated grids are shown in Fig. 9 in a narrower range of x (from 13.23 to 13.29) highlighting the peak itself. Note that we used varying spacing of data points from 0.001 or smaller for the peak itself (for most grids from 13.235 to 13.28) to 0.01 for the tails. At this x the best discretization of n = 1024 corresponds to kd = 0.026, which is about twice larger than the peak width.

 

Fig. 8 DDA simulations of the extinction efficiency versus x, corresponding to $b_{17}^{(1)}$ MDR with m = 1.6. Markers indicate simulated data (three representative n from 128 to 1024) and the solid lines show the fits by two Lorentzians with a cubic baseline (see text). Mie results are also shown together with a single-peak fit (solid line for simulated data and dashed line for fit); the corresponding two lines overlap on the whole range.

Download Full Size | PDF

 

Fig. 9 Same as Fig. 8, but showing the results for all levels of discretization and in the reduced range of x.

Download Full Size | PDF

The most striking feature of Figs. 8 and 9 is the peak splitting, clearly present for n≤ 256. In short, this can be explained by the symmetry breaking. For a sphere, the scattering quantities, including MDRs, are independent of the sphere orientation or incidence direction and polarization. Thus, the MDR of a sphere is a combination of several identical MDRs. When we discretize a sphere we reduce the particle symmetry from O(3) to O_h [34], then the MDR splits into several ones, and the separation decreases with increasing n. Similar effects appear when a sphere is morphed into a spheroid/ellipsoid [13]. Predicting the number of resonances may be possible through the group-theoretical analysis, accounting for the symmetry of the incident field (which changes with incidence direction), but we leave it for future research.

It is definitely a bad idea to fit these DDA results by a single Lorentzian (discussed below); hence we modify Eq. (3) to account for two peaks:

Another interesting feature of Fig. 8 is very different levels of DDA errors on the left and right tails (far from the peak). In particular, the difference between the DDA result and the Mie theory is 6 (n = 128) and 10 (n = 1024) times larger for x = 13.15 than for x = 13.35. We do not have a clear explanation for this fact, but a similar effect is also visible in Fig. 4 although distorted by the effect of the other neighboring MDRs.

The next step is to apply Richardson extrapolation to the fit parameters, however, some of the latter are determined with large errors from the non-linear fits (see Fig. 10). This issue is especially prominent for the largest n, due to fitting two strongly overlapping peaks; then variation of parameters of one peak can be almost compensated by that of another peak. This would hamper the extrapolation, as discussed in Section 3. Since both peaks converge to the Mie peak with increasing n, x_c1,2→x_c (its Mie value) and w_1,2→w and the same applies to the complex poles: x_c1,2−iw_1,2/2→x_c−iw/2 (the pole of $b_{17}^{(1)}$). The latter analogy also implies that the residues of two poles should add up to the residue of the Mie pole, i.e. A₁ + A₂→A. Therefore, we combine the parameters of the two DDA poles to have smoother convergence towards the Mie one. Specifically, we set the definition of A₀ = A₁ + A₂, x₀ as defined above, and w₀ = (A₁w₁ + A₂w₂)/(A₁ + A₂). By this, we average positions of two poles with weights A₁ and A₂, which guarantees that the averaged pole will have the same first two terms of asymptotic expansion at large distances as the sum of two original poles. To obtain standard errors of x₀, w₀, and A₀ we transformed Eq. (4) to have those as fit parameters instead of x_c1, w₁, and A₁ and performed the second fit of the same data (this affects neither the fit function nor the baseline). The results of this averaging are presented in Fig. 10. The large standard errors for parameters of individual peaks are greatly reduced, although noticeable errors remain for w₀ at the largest values of n.

 

Fig. 10 Comparison of the original parameters of two peaks derived from the DDA simulations, shown in Figs. 8 and 9 (for the MDR with m = 1.6 and x = 13.24417), to that of the averaged peak (see text for details). Parts (a), (b), and (c) show the peak positions, widths, and integrals, respectively.

Download Full Size | PDF

The averaged peak parameters are further used in the extrapolation procedure; the results are given in Fig. 11. In contrast to the previous section, we resort to the standard quadratic approximating function with n⁶ weights for residual squares, since here the linear term is definitely significant. The latter is caused by smaller size parameter of the particle and, hence, of the dipoles (for the same n). This extrapolation also ignores standard errors of fit parameters, as discussed in Section 3. We consider three different ranges of n: the full one [128,1024], the intermediate one [160,448] aiming at the best accuracy, and the fastest one [128,256]. The specification of the intermediate range is affected by the remaining uncertainties of determination of average peak parameters for the three largest grids (512–1024), which is especially visible for w₀, see Fig. 10(b). The ratios of total computational times for these ranges and for a single set of n = 1024 are approximately 45:5:1:26. The errors of the extrapolated values, i.e. the differences between c₀ and the Mie values, are shown in Table 2.

 

Fig. 11 Convergence of the average peak parameters (x₀, w₀, and A₀) derived from the DDA simulations for the MDR with m = 1.6 and x = 13.24417, with refining discretization (same as in Fig. 10). Solid lines correspond to the fits of c₀ + c₁n⁻¹ + c₂n⁻² in three different ranges of n, extrapolated to 1/n→ 0 (see text for details).

Download Full Size | PDF

Table 2. The absolute errors of the extrapolated values of average peak parameters for the MDR with m = 1.6 and x = 13.24417 (see Fig. 11) using different grid ranges of DDA simulations in comparison with the Mie result. The values for n = 1024 are from the direct fit of the DDA results.

View Table | View all tables in this article

The extrapolation works perfectly and regularly for x₀ [see Fig. 11(a) and Table 2], similar to the previous MDR, thanks to the averaging procedure. The resulting error compared to the Mie results is 22 (range [128,1024]), 18 (range [160,448]), and 6.6 (range [128,256]) times smaller than that for the largest n of the corresponding range. Moreover, the extrapolation error for the fastest range is still smaller than that for n = 1024 (1.0 × 10⁻³). For w₀ [see Fig. 11(b) and Table 2] the extrapolation leads to the similar error in the ranges [128,1024] and [160,448], about 3 times smaller than the smallest errors of single DDA simulations (that is similar for n = 448 and 1024), but the approximation is much more reliable (closer to the data points) for the range [160,448]. For A₀ the extrapolation is definitely not useful – it neither decreases the error nor describes the general trend, see Fig. 11(c). This can be explained by already small variation of A₀ due to proper averaging, in comparison to Fig. 5(c). In other words, most of the DDA variation (that is well described by a low-order polynomial) is removed during peak averaging, and the remaining part has a more complicated dependence. The relative error of A₀ for n = 320 is only 0.3% – good enough for any purposes. The same argument partly applies to w₀ – its relative accuracy is already 1% for n = 384.

The extrapolation quality for baseline parameters (Fig. 12 and Table 3) is similar to that for w₀ and A₀. The results are marginally useful for the full range, but are not robust for other ranges. We exemplify the latter by a single range [192,512], which produces small error for y₀ [Fig. 12(a)] but does not lead to reliable extrapolation for other baseline parameters. Similar to the previous MDR (Section 3), the baseline is well reproduced by a single DDA simulation (relative error of y₀ is within 0.04% for n≥ 160) due to the properly designed two-peak fitting of the simulated data.

 

Fig. 12 Same as Fig. 11 but for the baseline parameters (y₀, y₁, y₂, and y₃) and using other ranges of n.

Download Full Size | PDF

Table 3. Same as Table 2 but for the baseline parameters (y₀, y₁, y₂, and y₃) and using other ranges of n (the data is in Fig. 12).

View Table | View all tables in this article

We also applied two-peaks fitting to the DDA data for the m = 1.2 and x = 50.5369 from Section 3, but present only results for average peak parameters here (Fig. 13 and Table 4). It definitely helps for n = 128 and decreases the overall spread of the data points for n≤ 256, when compared to the single-peak fit. However, for n> 256 the fit errors of peak parameters significantly increase (as in Fig. 10 but more pronounced), which is not fully resolved by peak averaging (see Fig. 13). This is also accompanied by increased standard errors of baseline parameters (data not shown). Applying the same cubic extrapolation as in Section 3 (Fig. 5) in the full range of n (Fig. 13), we seem to capture the general trend but get systematically biased results for c₀ (see Table 4). This is probably caused by biases (and fit errors) for the larger n due to the overfitting. However, applying extrapolation in the range [128,320] (where the two-peaks fit works fine) we get more reliable and accurate results than that for the single-peak fit (compare Table 4 to Table 1). The extrapolation in the range [256,768] is unsatisfactory, similar to that for the full range (data not shown).

 

Fig. 13 Convergence of the average peak parameters (x₀, w₀, and A₀) derived by two-peaks fit from the DDA simulations, shown in Fig. 4 (for the MDR with m = 1.2 and x = 50.5369), with refining discretization. Solid lines correspond to the fits of c₀ + c₂n⁻² + c₃n⁻³ in two different ranges of n, extrapolated to 1/n→ 0 (the same procedure as in Fig. 5).

Download Full Size | PDF

Table 4. Same as Table 3, but for the MDR with m = 1.2 and x = 50.5369, based on the data in Fig. 13.

View Table | View all tables in this article

Similar to the discussion for the range [128,320] above, the single-peak fit for the MDR with m = 1.6 and x = 13.24417 is less accurate than or comparable to the corresponding two-peak fit, in terms of both the fit parameters and further extrapolation errors.

5. Guidelines for simulating MDRs with the DDA

The previous two sections are devoted to development of the processing algorithms for the DDA, with a focus on rationale, i.e. why particular design choices were made. Those design choices are interweaved with the data for two specific MDRs, and the data was obtained in batches in parallel with the development of data processing. Our goal in this section is to generalize this complicated experience and produce guidelines for efficient simulation of the MDRs with the DDA. Most of them are relevant to other simulation methods as well.

The algorithms consist of data simulating (sampling a range of x), peak fitting, and extrapolation. First, there is no need to sample the peak maximum in details, since the peak is further fitted as a whole. Instead, about 30 uniformly spaced values of x inside the main peak body (e.g., down to 10% of its value) can be sufficient. The amount of the peak tails to cover (to the left and right of the main body) is determined by the presence the neighboring peaks. If the simulated MDR is well isolated (relative to its width), up to 10w can be sampled on each side; then the sampling of the tails should be much coarser to have about 10 points for each tail (as in Section 4). By this, one can hope for good accuracy for the baseline parameters, which will in turn improve the accuracy of peak parameters. Otherwise, the tails can be sampled up to the minima between the peaks – the goal is to be able to describe the baseline by a cubic polynomial. In this case, the range of x should also cover the peak details for all used n in addition to the Mie peak (see Section 3).

It is also important to have all simulation parameters fixed, when varying x, including those that affect the accuracy. Fortunately, for the DDA this boils down to fixing n, while other parameters, such as the convergence threshold of the iterative solver are typically small enough. Actually, our DDA simulations show quasi-random oscillations with amplitude of the order of this threshold (10⁻⁵), but it can be safely neglected in processing. A related issue is the volume correction (i.e. using x_V or x_D), discussed in Section 3. While the DDA results can be trivially transformed from one convention to the other (and we showed that x_V is more robust), most importantly is to use a single convention during the whole processing. The presence of a single parameter determining the simulation accuracy is also critical for further extrapolation. This is not necessarily easy to achieve with other simulation methods.

When fitting the simulated data, the number of peaks (resonances) should be determined empirically, since both under- and over-fitting worsens the accuracy of parameters. It is easiest to fix the number of peaks for all the data, but it can also be varied with n. Some statistical test can be used for the latter, but we have not tried it. If two peaks are fitted, the averaging of corresponding poles should be performed: pole residues (peak integrals) are summed and pole positions are averaged with weights given by the residues (see Section 4). This procedure can be easily generalized to any number of poles. The order of a polynomial approximating the baseline should be determined empirically, by increasing it and monitoring the standard errors of fit parameters. Third order seems to be optimal for many cases, but that is also related to the simulated range of x (length of tails), discussed above. By using the Lorentzian peak shape we assumed that pole residue for the Mie coefficient b_n is purely imaginary, since the peak of Q_ext is proportional to Re(b_n) [Eq. (1)]. The fit can be improved by introducing the real part of this residue as an additional fit parameter (adding asymmetry to the peak), but we have not systematically tested this option yet.

The extrapolation technique is closely related to the used values of n. It is desirable to reach kd of ~2w, i.e. n≈x/w, then the result for the largest n is already a satisfactory approximation for the MDR. It can be further improved by performing extrapolation. The reasonable choice for spacing of values of n is approximately uniform in the logarithmic scale, in particular, for any n₀ = 4q the suitable values of n in the range from n₀ to 2n₀ are {4q,5q,6q,7q,8q} in agreement with previous studies [28,30,32], but we have not tried other options. The best extrapolation is expected using the range from n₁ to n_max, where n_max is the largest achievable n and n₁ is the smallest value that still falls on the same trend. The principal choice is that of an approximating polynomial of 1/n. If the linear trend is noticeable in the considered range of n, then the standard quadratic function should be used (as in Section 3), otherwise a cubic one without a linear term is more suitable (as in Section 4). There are three free parameters and the expected residuals are of the order of the first omitted terms (n⁻³ and n⁻⁴, respectively), i.e. the weights of the residual squares are n⁶ and n⁸, respectively.

Those guidelines combined with a large compute cluster are rather powerful, as exemplified by the results of Sections 3 and 4. But even they have practical limitations. The narrower the resonance is – the larger n is needed. Assuming that the current computational resources can support multiple simulations with n≈10³ one can try to attack MDRs with relative widths down to 10⁻⁴, such as deemed impossible for the DDA in [13]. This hope is also based on the (not yet proved) hypothesis that the extrapolation of the peak parameters will work better for narrower MDRs due to the reduced influence of neighboring MDRs (baseline). Still, this is far from the MDR with x = 20.0001 and m = 1.6 discussed in Section 2, which has relative width around 10⁻⁵. The first-order resonances classified by Chýlek et al. [3] seem completely unreachable for the DDA.

We are not aware of any practical reason to simulate such ultra-narrow resonances except for extreme testing of the method. However, the DDA can be used to simulate MDRs for arbitrary non-spherical shapes. Generally, those MDRs are expected to be relatively wide and, thus, trivial to simulate. But if the particle is only slightly different from a sphere or a spheroid supporting ultra-narrow resonances [14], its MDR may again become narrow making the above guidelines highly relevant. In particular, they can enable the DDA to study the evolution of the MDR under non-axisymmetric deformations of a sphere or a spheroid to complement existing results for axisymmetric deformations [3,13,14].

6. Conclusion

This study investigates the performance of the DDA in reproducing the MDRs. The narrower the resonance is – the smaller the dipole size is required for accurate simulations. Thus, the brute-force application of the DDA reaches a satisfactory approximation for the MDRs only if the dipole size is smaller than twice the MDR width, requiring enormous computing resources and time.

We proposed a more sophisticated approach, consisting of DDA simulations in a range of x, Lorentzian fitting, and extrapolation. Both Mie reference and the DDA simulations of Q_ext can be well fitted with a Lorentzian combined with a cubic baseline, but for coarser discretization grids two peaks have to be included into the fitted model. The latter is caused by the symmetry breaking due to the discretization, and its prominence significantly differs between the two considered MDRs. The corresponding peak and baseline parameters enable robust quantification of the DDA accuracy. Moreover, in case of two peaks the weighted peak averaging (based on averaging of the underlying complex poles) helps to reduce the standard errors of individual peaks parameters, especially when the peaks largely overlap. Relatively smooth convergence of fit parameters is further exploited by the Richardson extrapolation (to 1/n→ 0) to further improve the accuracy. We have summarized this approach into step-by-step guidelines, which partly apply to other simulation methods.

We tested this workflow on two MDRs with relative widths 2 × 10⁻³ and 9 × 10⁻⁴ using massive simulations with n up to 768 and 1024, respectively. After extrapolation the Mie reference can be approached with extreme accuracy: peak positions are reproduced with errors less than 0.07% and 0.4% of their widths, respectively. Alternatively, if only the simulations with n< 320 and 256 are considered (simulation times are 11 and 45 times faster than that for full ranges of n), the peak positions are still reproduced satisfactory: within 0.8% and 7%, respectively. Thus, it is feasible to satisfactory simulate MDRs with relative width of 10⁻⁴ if the computational hardware supports n = 2048 (requires about 6 TB of memory). Addressing MDRs with relative widths below 10⁻⁵ would require more than a PB, making it not feasible at this time. More practical, however, is not to reproduce the well-known Mie results but rather to use this DDA workflow to study the evolution of the MDR under non-axisymmetric deformations of a sphere or a spheroid.

The DDA with standard discretization level is inadequate for most MDRs, the same is expected for any volume- or surface-discretization methods. Thus the last conclusion is that MDRs should be explicitly avoided in systematic sphere-based benchmarks for any simulation methods, unless they are aimed at MDRs themselves. More specifically, x = 20 and m = 1.6 is a combination of round numbers (easy to stumble upon) that should not be used.