1. Introduction

Accurate and efficient computation of light scattering by dielectric particles is important in various fields of research and applications, such as radiative transfer, remote sensing [1,2], bio-optics [3] and nanophotonics [4]. In an event of light scattering by a particle, the incident electromagnetic wave interacts with the particle to stimulate the scattered field. The scattered field propagates, starting from the particle, in all directions. The incident and scattered fields can be expanded into infinite series in terms of vector spherical wave functions (VSWFs). The expansion coefficients of the incident and scattered fields are related by a matrix called the T-matrix. The single-scattering properties, including the extinction and scattering cross sections and the scattering phase matrix, can be analytically derived from the T-matrix [2].

The extended boundary condition method (EBCM) [5,6], based on the electromagnetic surface integrals, is used to compute the T-matrix of a scattering particle. Extensive efforts are devoted to developing and improving EBCM numerical algorithms [2,7–11]. The EBCM has long been regarded as a synonym of the ‘T-matrix method’ in the literature. However, some other methods were also developed to compute a T-matrix such as the superposition T-matrix (STM) method [12,13] and the invariant imbedding T-matrix (IITM) method [14–17]. Other numerical methods to solve the light scattering problem include the Discrete Dipole Approximation [18–20] and Finite Difference Time Domain [21–23] methods. The T-matrix methods have a unique advantage over other methods to calculate single-scattering properties for randomly-oriented particles [7], such as for applications to remote sensing of clouds and aerosols. Because the T-matrix for a specific orientation can be transformed to the T-matrix for an arbitrary orientation with respect to a laboratory coordinate system through known analytical relations [2], this unique property of the T-matrix allows us to analytically average the single-scattering properties over particle random orientations [7,8] that are common for particles in nature.

In contrast to EBCM and STM, the IITM method solves an electromagnetic volume integral equation by converting a boundary-value problem to an initial-value problem [14]. The electromagnetic volume integral equation is a Fredholm equation of the second kind, while the surface integral equation is a Fredholm equation of the first kind [24]. For numerical implementation, the second kind is generally more stable than the first kind. Practically, the EBCM numerical algorithm [2,25,26] is applicable mainly to symmetric particles such as spheroids, cylinders, and cubes [27]. In recent years, numerically stable algorithms based on the IITM method have been developed to compute the single-scattering properties of nonspherical particles with arbitrary geometries and inhomogeneities [28–31].

In theory, the expansion coefficients for the incident and scattered fields are infinite series. The infinite T-matrix in numerical computation must be truncated to a finite dimension by terminating the series at an integer degree N. As discussed in section 2.2, a larger N results in a matrix which better approximates the infinite matrix, but computational time increases rapidly with N. In the IITM algorithm, the T-matrix dimension is set a priori because increasing the truncation number from N to $N + 1$ requires repeating the entire algorithm with larger matrices, not by simply computing another term. Thus, we must determine an optimal VSWF degree to truncate the T-matrix, optimal in the sense that it achieves a balance between accuracy and efficiency. Wiscombe [32] proposes an empirical formula to determine the appropriate VSWF degree in an efficient computational algorithm based on the Lorenz-Mie theory for homogeneous spheres. Because T-matrix methods are based on VSWF expansion, the Lorenz-Mie theory can be thought of a special case of the T-matrix method. T-matrix computations [2,15,33,34] use empirical formulas similar to the one proposed by Wiscombe to determine the VSWF degree. In particular, Zhai et al. [33] and Hu et al. [34] slightly adjust Wiscombe’s formula in the IITM computations for spheroids and hexagonal columns without detailed explanation and validation. In Wiscombe (1980) [32], the predetermined optimal truncation number is a function of only the size parameter x, which is defined as the product of the radius of the scattering sphere ($r$) and the wavenumber ($2\pi /\lambda $). The IITM algorithm involves nonspherical particles and is substantially different from the Lorenz-Mie algorithm [35]. For the IITM computation, the truncation number given by Wiscombe’s formula may be either inefficient (more expansion terms than necessary for desired accuracy) or inaccurate (insufficient expansion terms). Thus, it is imperative to investigate the relationships between the optimal truncation number and the size parameter, refractive index, and particle shape for efficient and accurate IITM computations.

This study aims to explore the factors determining the optimal T-matrix truncation number in the IITM algorithm. We first briefly introduce the IITM method and raise the major issues concerned in its numerical implementation in Section 2. In Section 3, we investigate the impact of the size parameter and the index of refraction on the selection of an optimal truncation number in a case of a homogeneous hexagonal column particle. In Section 4, we discuss particle shape effects on the optimal T-matrix truncation number. Section 5 summarizes the findings.

2. Invariant imbedding T-matrix method

2.1 Theory and algorithm

For a light scattering problem, the incident and scattered electric fields can be expanded in terms of VSWFs [2]:

The T-matrix transforms the expansion coefficients of the incident electric vector to the scattered wave counterparts,

The IITM method is based on a recursive algorithm [36], given by [15] as

Figure 1 illustrates the IITM algorithm. The computation starts from the inscribed homogeneous sphere of a particle and ends with the circumscribed sphere surrounding the particle. In Eq. (4), ${\mathbf {T}}({{r_l}} )$ is the T-matrix of the portion of the particle enclosed by a spherical shell of radius ${r_l}$. Equation (4) computes ${\mathbf {T}}({{r_l}} )$ given ${\mathbf {T}}({{r_{l - 1}}} )$ (where ${r_l} > {r_{l - 1}}$) and particle shape information at radius ${r_l}$. The initial T-matrix in the recursive algorithm is either a zero matrix for ${r_1} = 0$ or is computed by the Lorenz-Mie theory for the inscribed sphere. After the recursive computation is completed, ${\mathbf {T}}({{r_{\textrm{max}}}} )$ is the T-matrix of the whole particle.

 

Fig. 1. Illustration of IITM algorithm. ${r_{\textrm{max}}}$ is the circumscribed sphere radius, ${r_1}$ is the inscribed sphere radius, and ${r_l}$ is the radius of a sphere enclosing a portion of the particle in recursion step l.

Download Full Size | PDF

The particle size measured with respect to the incident wavelength can also be specified in terms of the size parameter defined as $x = k{r_c}$. For a spherical particle, ${r_c}$ is the radius. For a nonspherical particle, ${r_c}$ is the radius of the smallest sphere circumscribing the particle.

In Eq. (4), ${{\mathbf {Q}}_{jh}}({j,h = 1,2} )$ are matrices involving Bessel functions of the first kind (matrix ${\mathbf {J}}$) and Hankel functions of the first kind (matrix ${\mathbf {H}}$):

The matrix ${\mathbf {U}}({{r_l}} )$ in Eq. (6) specifies the particle shape and the index of refraction. Its elements are angular integrals over the spherical shell at radius ${r_l}$. The dimension of ${\mathbf {U}}$ is $N({N + 2} )\times N({N + 2} )$, and each element of ${\mathbf {U}}$ is a $3 \times 3$ matrix [15,37]:

After the T-matrix of a particle is obtained, we can analytically average the particle single-scattering properties, namely, the scattering phase matrix (${\mathbf {P}}$, including the phase function ${\textrm{P}_{11}}$), the extinction cross section (${C_{\textrm{ext}}}$), and the scattering cross section (${C_{\textrm{sca}}}$), over either random orientations or preferred orientations. Detailed discussions of the relationships between the particle single-scattering properties and the T-matrix can be found in Chapter 5.5 and 5.6 of Mishchenko et al. [2].

2.2 Numerical computational accuracy and efficiency

The accuracy of IITM computation is determined by the dimensions of the T-matrix or the maximum VSWF degree, the radial resolution (number of thin spherical shells) in the recursion, and the precision of $\mathbf{U}$ elements that describe the particle shape and refractive index.

The computation of the ${\mathbf {U}}$ matrix elements involves integration on a spherical surface. Practically, the integral along the $\varphi $ direction can be evaluated analytically by finding the intersections of the spherical surface with the particle. A geometric approach is developed by Bi and Yang [16] to find the intersections. The analytical algorithm developed in [38] is used to compute the ${\mathbf {U}}$ matrix of axially symmetric particles. The integral along $\theta $ contains the Wigner-d function and its derivative and must be evaluated numerically. Zhai et al. [33] discusses the accuracy and efficiency of different quadrature schemes in computing any ${\mathbf {U}}$ matrix element. Hu et al. [34] proposes an empirical formula to determine the number of Gauss-Legendre (GL) quadrature points in evaluating the integral along $\theta $ in ${\mathbf {U}}$ matrix computation, based on the numerical results for spheroidal particles.

The recursion equation Eq. (4) is indeed a discretization of a differential matrix equation, as proved by Johnson [14] and Sun et al. [17]. With a finer radial resolution defined as $\mathrm{\Delta }r = {r_l} - {r_{l - 1}}$ or equivalently, more recursive computations, the T-matrix result becomes more accurate. The numerical experiments conducted by Hu et al. [34] show that in the case of spheroidal shapes, using Gaussian points with the same number of recursive computations to determine ${r_l}$ can be more accurate than specifying ${r_l}$ with a constant spacing. Based on numerical results of spheroidal particles, Hu et al. [34] reports an empirical formula to determine the number of Gaussian radial points in recursive computations. We have nevertheless chosen to use uniform spacing, as discussed below.

The T-matrix dimensions, or maximum VSWF degree N, determine the multipoles used to represent the scattered field, which provide clear physical insights such as the edge effect and the transfer of angular momentum [39]. When we truncate the T-matrix up to the maximum VSWF degree N, N indeed determines the VSWF expansion of the incident and scattered fields, in which only the terms with degrees from $1$ to N are kept. The information carried by higher-degree terms (from $N + 1$ to infinity) is lost. Thus, N must be large enough to guarantee that the included terms can account for the physical process of scattering. The uncertainties attributed to omitting the higher-degree terms must be negligible. The optimal N that meets the above criterion is named ${N_{\textrm{ideal}}}$ here. ${N_{\textrm{ideal}}}\; $ should be independent of the algorithms used to compute a T-matrix. Because the algorithms such as IITM and EBCM have different numerical computational errors, an optimal N chosen in numerical computations (here named as ${N_{\textrm{num}}}$) may be different among various numerical algorithms. ${N_{\textrm{num}}}$ is expected to be larger than ${N_{\textrm{ideal}}}$, because more expansion terms are needed to reduce the uncertainties in numerical computations. In the remaining text, the optimal N refers to ${N_{\textrm{num}}}$ in the IITM computation if not specifically mentioned.

The CPU time and memory required by the IITM method are mainly spent on matrix operations including addition, multiplication, and inversion. Figure 2 shows the CPU time and memory requirements for IITM computations on the Intel x86-64 Linux cluster at the Texas A&M High Performance Research Computing (HPRC) Facilities. As shown in Fig. 2(a) and Fig. 2(b), both the memory requirement and CPU time increase nearly exponentially with increasing N for the various particle shapes shown in Fig. 2(c). We applied the refractive index of ice at 660 nm as the default refractive index setting [16]: $\tilde{m} = 1.3078 + 1.67 \times {10^{ - 8}}i\; $. Since the imaginary part is small, we simply set the imaginary part to $0$ in our computations shown in Figs. 2 and 3. The aspect ratio of a hexagonal column is defined as the ratio of the maximum cross-section width ($2a$) to the height $(h )$, as shown in Fig. 2. Although particle symmetry can be exploited to reduce the computational burden such as in the case of a hexagonal column, the memory requirement and total CPU time are substantial when N is large. Therefore, it is important to identify an optimal maximum VSWF degree or truncation number to economize on computational resources, while we retain the computational accuracy.

 

Fig. 2. (a) The memory requirement and (b) CPU time cost for particles with different shapes in IITM computation with increasing N. (c) Particle shapes. All particles have the same size parameter $x = 25$ and refractive index $\tilde{m} = 1.3078 + 0i$. The hexagonal column has an aspect ratio $2a/h = 0.8$.

Download Full Size | PDF

 

Fig. 3. Radial resolution test for a single hexagonal column with an aspect ratio of $2a/h = 0.8$ and 8-column aggregate. In these computations, $x = 25$ and $\tilde{m} = 1.3078 + 0i$. The wavelength is assumed to be $2\pi \; \mu m$ in the computation. (a) Relation between computed ${C_{\textrm{sca}}}$ and truncation number N for different radial resolutions $\Delta x$ in the single hexagonal column case. (b) Relative change in ${C_{\textrm{sca}}}$, or ${\textrm{R}_{\textrm{Mi}02}}[{C(N )} ]$, from truncation number $N - 1$ to N for different $\Delta x$ in the single hexagonal column case. (c-d) Same as (a-b), except that the particle shape is an 8-column aggregate.

Download Full Size | PDF

This study focuses on finding an optimal maximum VSWF expansion degree in IITM computation, for various particle sizes, refractive indices, and particle shapes. In the following numerical experiments, the number of GL quadrature points in the ${\mathbf {U}}$ matrix computation must be sufficient to have converged results. As shown in Zhai et al. [33], with more than 200 GL quadrature points, the relative error of the extinction efficiency can be smaller than $0.01\%$ for hexagonal columns.

As mentioned above, the radial points in our recursion computation are selected with equal spacing. Taking radial resolution to be defined as $\Delta x = {x_l} - {x_{l - 1}}$, we investigate how the IITM calculated scattering cross section depends on $\Delta x$ over a range of truncation numbers in the case of two different particles. Results are shown in Figs. 3(a) and 3(c), from which it is apparent that $\Delta x = 0.1$ is sufficient in the cases of the homogeneous hexagonal column and 8-column aggregate to obtain converged scattering cross section results.

In our study, we consider convergence criteria that are formulated in terms of the relative change in computed value of some optical quantity or quantities as N increases by 1. We first consider the measure ${\textrm{R}_{\textrm{Mi}02}}[{C(N )} ]$ proposed by Mishchenko et al. [2]

3. Convergence and truncation of the T-matrix in IITM method

3.1 Motivation for T-matrix truncation

The truncation of VSWFs and convergence in Lorenz-Mie and T-matrix computations has been discussed in many previous studies [32,40–45]. Wiscombe [32] proposes an empirical formula (hereafter indicated with subscript ‘Wis80’) to determine the maximum VSWF degree as a function of size parameter x,

T-matrix methods are also based on VSWF expansion. For spherical particles, the Lorenz-Mie series coefficients ${a_n}$ and ${b_n}$ are equal to T-matrix elements $- {\delta _{mm^{\prime}}}{\delta _{nn^{\prime}}}T_{mnm^{\prime}n^{\prime}}^{22}$ and $- {\delta _{mm^{\prime}}}{\delta _{nn^{\prime}}}T_{mnm^{\prime}n^{\prime}}^{11}$, respectively.

In the EBCM, Mishchenko et al. [2] uses a formula (hereafter shown with subscript ‘Mi02’) that is similar to Eq. (12),

For the IITM method, Zhai et al. [33] (Zh19) states that the constant term in Eq. (13) can be smaller for spheroids and hexagonal columns, given by

Hu et al. [34] (Hu20) determine optimal N by multiplying ${N_{\textrm{Mi}02}}$ by a constant $1.1$, or

In the IITM recursive algorithm, the radius of an intermediate sphere ${r_l}$ is smaller than the circumscribed sphere radius ${r_{\textrm{max}}}$. We can specify different N for individual recursive steps depending on the size parameters of the intermediate spheres, ${x_l} = k{r_l}$. This can save more than half of CPU time in the recursive computation, but uneven results of smaller N values in intermediate spheres may be the reason for non-monotonic decreases of a convergence measure as the final truncation number increases, such as in Fig. 3(b).

Following Mishchenko et al. [2] we could judge the computation (here by the IITM) to have “converged” when

That is, the computation has converged if neither ${C_{\textrm{ext}}}(N )$ nor ${C_{\textrm{sca}}}(N )$ changes by $< 0.1\%$ from truncation number $N - 1$ to N. We note that, while one might hope that the same is true for larger values of N, such monotonic decrease in ${\textrm{R}_{\textrm{Mi}02}}[{C(N )} ]$ does not occur for all particles (for example, see Fig. 4).

 

Fig. 4. Convergence criteria $\Delta [{{\textrm{P}_{11}}(N )} ]$ and ${\textrm{R}_{\textrm{Mi}02}}[{C(N )} ]$ with increasing N in IITM computations. The particle is a hexagonal column with $2a/h = 0.8,x = 50,\; \tilde{m} = 1.3078 + 0i$.

Download Full Size | PDF

Using this tentative measure of relative improvement with increasing N, we now examine how useful the analytic determination of N provided by Eq. (13) is, a determination that involves only the size parameter. The panels in the upper row of Figs. 5 (a-c) show how ${\textrm{R}_{\textrm{Mi}02}}[{C(N )} ]$ varies for homogeneous hexagonal columns of varying size parameters and indices of refraction. As can be seen, ${\textrm{R}_{\textrm{Mi}02}}[{C(N )} ]$ varies significantly when the particle size parameter and refractive index vary. For particles with a larger real part of the refractive index, the convergence criterion is not satisfied. Evidently Eq. (13) does not correctly predict the relation between the truncation number needed in the IITM computations and the size parameter, and the equation must be modified to consider its refractive index dependence. Before considering such a modification, we turn to an improved measure of convergence.

 

Fig. 5. (a-c) ${\textrm{R}_{\textrm{Mi}02}}[{C(N )} ]$ defined in Eq. (11) as a function of size parameter and refractive index from IITM computational results of a hexagonal column with aspect ratio $2a/h = 0.8$; (d-f) Similar plots as (a-c), but for $\Delta [{{\textrm{P}_{11}}(N )} ]$ defined in Eq. (17a). In all panels $N = {N_{\textrm{Mi}02}}\; $ and are obtained by Eq. (13). In (a) and (d), all particles have $\tilde{m} = 1.3078 + 0i$; in (b) and (e), all particles have $x = 50,\; {m_\textrm{i}} = 0$; and in (c) and (f), all particles have $x = 50,\; {m_\textrm{r}} = 1.3078$.

Download Full Size | PDF

3.2 Convergence criterion for phase matrix

The convergence criterion phrased in terms of the measure in Eq. (11) considers only the extinction and scattering cross sections. As reported in [44], other single-scattering properties, including the scattering phase matrix, may have different convergence characteristics. The accurate computation of the scattering phase matrix ${\mathbf {P}}$ of a particle, and in particular the scattering phase function or matrix element ${\textrm{P}_{11}}$, is of great importance in remote sensing applications. We therefore consider an alternative measure of convergence that focuses on ${\textrm{P}_{11}}$:

The lower row of panels in Figs. 5 (d-f) shows how $\Delta [{{\textrm{P}_{11}}(N )} ]$ varies, again with $N\; $ chosen using Eq. (13). Figure 4 shows the variation of $\Delta [{{\textrm{P}_{11}}(N )} ]$ and ${\textrm{R}_{\textrm{Mi}02}}[{C(N )} ]$ with increasing N. $\Delta [{{\textrm{P}_{11}}(N )} ]$ tends to be a stricter criterion than ${\textrm{R}_{\textrm{Mi}02}}[{C(N )} ]$. In most cases, when $\Delta [{{\textrm{P}_{11}}(N )} ]\le \varepsilon $ is satisfied, ${\textrm{R}_{\textrm{Mi}02}}[{C(N )} ]\le \varepsilon $ is also satisfied. This may not be surprising because the extinction and scattering cross sections are integrated quantities, and errors in integrands often cancel out in the integration.

According to our numerical experiments using hexagonal column particles, the convergence criterion $\Delta [{{\textrm{P}_{11}}(N )} ]\le 0.001$ is sufficient for determining the convergence of all the single-scattering properties. As an example, Fig. 6 shows the nonzero scattering phase matrix elements computed with different N for a hexagonal column. The convergence criterion $\Delta [{{\textrm{P}_{11}}(N )} ]\le 0.001$ is satisfied at $N \ge 138$. The scattering phase matrix elements with $N = 158$ are almost identical to those with $N = 138$. We also compare the results with cases $N = 120\; $ predicted by ${N_{\textrm{Wis}80}}$ and $N = 126$ predicted by ${N_{\textrm{Mi}02}}$, but both N values produce large deviations in the backscattering angle range and nondiagonal scattering phase matrix elements. In the following experiments, we use the convergence criterion $\Delta [{{\textrm{P}_{11}}(N )} ]\le 0.001$ to find the optimal truncation number because it is more stringent than ${\textrm{R}_{\textrm{Mi}02}}[{C(N )} ]$.

 

Fig. 6. Scattering phase matrix elements computed by IITM with different truncation numbers N. The convergence criterion $\Delta [{{\textrm{P}_{11}}(N )} ]\le 0.001$ is satisfied at $N = 138$ (red dashed curves). The particle is a hexagonal column with $2a/h = 0.8,x = 100,$ and $\tilde{m} = 1.3078 + 0i$.

Download Full Size | PDF

Figures 3 and 4 show that neither of the convergence measures strictly monotonically decrease with the truncation number N. It is possible that although the convergence criterion is met, the computational result is indeed biased from the truth if the first N is chosen that makes the convergence criterion satisfied. We find, from the numerical experiments, that the oscillation of $\Delta [{{\textrm{P}_{11}}(N )} ]$ is much weaker than ${\textrm{R}_{\textrm{Mi}02}}[{C(N )} ]$, such as in Fig. 4. If we choose the first N that satisfies $\Delta [{{\textrm{P}_{11}}(N )} ]\le 0.001$, it is less likely that the chosen N is not sufficient than using ${\textrm{R}_{\textrm{Mi}02}}[{C(N )} ]\le 0.001$ to choose N. Furthermore, to avoid the case where an insufficient N is chosen, we compare the scattering phase matrix result with the chosen N to the result with a larger N. If the two results are consistent with negligible difference, we assume the N is properly selected. For example, in Fig. 6, the chosen N is $138$. The scattering phase matrix results with $N = 138$ and a larger value $N = 158$ are consistent.

3.3 Empirical formula to determine optimal truncation number

According to the results shown in Figs. 4–6, the truncation number predicted by the empirical formulas Eqs. (12)-(15) may not be sufficient for IITM computations. The single-scattering properties of nonspherical particles are determined by three parameters, namely the size parameter, the index of refraction, and particle shape [48]. The optimal truncation number should also be related to these three factors. In particular, the optimal truncation number is highly dependent on the refractive index, which is not considered in previous studies. We thus propose an empirical formula of the form

We first consider only one particle shape, namely a hexagonal column with aspect ratio $0.8$. IITM numerical experiments are performed with size parameters from $5$ to $100$ (with an interval of $5$) and real part of refractive indices $1.01$ and $1.05$ to $2.0$ with an interval of $0.05$, to determine the optimal truncation number ${N_{\textrm{opt}}}$ according to the ${\textrm{P}_{11}}$-based convergence criterion expressed in Eq. (17b). The imaginary part of the index of refraction is assumed to be $0$. The optimal truncation numbers are plotted in Fig. 7(a). Figure 7(a) clearly shows the dependence of N on both the size parameter and real refractive index. ${N_{\textrm{Wis}80}}$ overestimates the optimal truncation number for the cases with refractive indices less than $1.1$ and underestimates the optimal truncation number for the cases with refractive indices larger than $1.2$.

 

Fig. 7. (a) Black contours with labeled numbers show the optimal truncation number ${N_{\textrm{opt}}}$ for hexagonal columns with an aspect ratio $2a/h = 0.8$ as a function of the size parameter and the real part of the refractive index. Colored contours with the colorbar represent the difference ${N_{\textrm{opt}}} - {N_{\textrm{Wis}80}}\; .$ (b) Relative difference $1 - {N_{\textrm{new},\textrm{hc}}}/{N_{\textrm{opt}}}$ for truncation number computation using the new estimation formula. (c) Absolute difference ${N_{\textrm{new},\textrm{hc}}} - {N_{\textrm{opt}}}$ for truncation number computation using the new estimation formula. The dashed contour lines show $|{N_{\textrm{new},\textrm{hc}}} - {N_{\textrm{opt}}}|= 2$.

Download Full Size | PDF

Then, we determine the three constants in Eq. (18) by fitting the optimal truncation numbers. Because N is an integer, small changes of ${c_1}$ and ${c_2}$ may not affect the results in the ranges of x and ${m_\textrm{r}}$ where IITM is applicable. In the fitting process, ${c_1}$ varies from $0.60$ to $1.20$ with step $0.01$, ${c_2}$ varies from $1.0$ to $8.0$ with step $0.1$, and ${c_3}$ is an integer varying from $- 3$ to $5$. We search numerically among combinations of ${c_1}$, ${c_2}$ and ${c_3}$ to minimize the root mean square error (RMSE) of Eq. (18) compared to the optimal truncation number obtained in the numerical experiments. The resulting formula is

Figure 7 shows comparisons between IITM-calculated values of optimal truncation ${N_{\textrm{opt}}}$ using criterion expressed in Eq. (17b) and the two analytical formulas ${N_{\textrm{Wis}80}}$ and ${N_{\textrm{new},\,\textrm{hc}}}$ . Figure 7(a) shows contours of constant IITM-calculated ${N_{\textrm{opt}}}$ and uses color to show the difference ${N_{\textrm{opt}}} - {N_{\textrm{Wis}80}}$. In Fig. 7(c), the RMSE of the difference ${N_{\textrm{new},\textrm{hc}}} - {N_{\textrm{opt}}}$ is $2.01$, and we can see from the figure that the absolute error of the estimate is mostly within two integers. Figure 7(b) shows the relative error between truncation numbers predicted by Eq. (19) and those obtained from numerical experiments, while Fig. 7(c) shows the difference of truncation numbers. For small particles with a large refractive index, Eq. (19) evidently underestimates the required truncation number. We suggest increasing the constant term to ensure convergence for these cases.

We also calculate the ${N_{\textrm{opt}}}$ for absorptive particles with different imaginary parts of the refractive index. The values ${N_{\textrm{opt}}} - {N_{\textrm{new},\textrm{hc}}}$ for various size parameters, the real and imaginary parts of the refractive index are given in Table 1. For $0 < {m_\textrm{i}} \le 0.2$, most of the values $|{N_{\textrm{opt}}} - {N_{\textrm{new},\textrm{hc}}}|$ are less than $5$. According to the results listed in Table 1, although Eq. (19) does not consider the imaginary part of the refractive index, and is fitted based on a single value ($0$) of the imaginary part, it still can be applied to IITM computations of particles with imaginary part of refractive index less than $0.2$. For particles with small refractive index and strong absorption (e.g., $\tilde{m} = 1.2 + 0.5i$), the measure of convergence Eq. (17b) may not be quite as useful. This is because the values of ${\textrm{P}_{11}}$ in backscattering angles are very small, which leads to large relative differences. In these cases, as shown in Fig. 5(c) and Fig. 5(f), $\Delta [{{\textrm{P}_{11}}(N )} ]$ values become larger than nonabsorptive cases with the same truncation number, while ${\textrm{R}_{\textrm{Mi}02}}[{C(N )} ]$ values do not.

Table 1. The values ${N_{\textrm{opt}}} - {N_{\textrm{new},\textrm{hc}}}$ for absorptive hexagonal column with $2a/h = 0.8$

View Table | View all tables in this article

4. Effect of particle shape

Section 3 only considers the case of a single shape, a hexagonal column with aspect ratio $0.8$. In the following numerical experiments, we consider two other shapes, namely a hexagonal plate with aspect ratio $2.0$ (i.e., a hexagonal plate) and an eight-hexagonal column aggregate (see Fig. 2(c)) defined in Yang et al. [49]. IITM numerical experiments are conducted with size parameters ranging from $5$ to $60$ with an interval of $5$, and real refractive indices $1.01$ and $1.05$ to $2.0$ with an interval of $0.05$, and determine the optimal truncation number according to the convergence criterion Eq. (17) for the hexagonal plate, while for the aggregate shape, size parameters vary from $5$ to $50$ with an interval of $5$ and real refractive indices $1.01$ and $1.1$ to $2.0$ with an interval of $0.1$.

We find ${N_{\textrm{opt}}}$ values for the hexagonal plate and aggregate according to the convergence criterion Eq. (17) at individual size parameters and indices of refraction. For the two shapes, we assume that their ${N_{\textrm{opt}}}$ values are dependent on x and ${m_\textrm{r}}$ in the same manner as Eq. (18). Then, we obtain the coefficients ${c_1}$, ${c_2}$ and ${c_3}$ for the two shapes by minimizing the RMSE over the suite of calculations. The resulting empirical formulas for the hexagonal plate (${N_{\textrm{new},\textrm{hp}}}$) and aggregate (${N_{\textrm{new},\textrm{agg}}}$) are

The differences ${N_{\textrm{new},\textrm{hp}}} - {N_{\textrm{opt}}}$ and ${N_{\textrm{new},\textrm{agg}}} - {N_{\textrm{opt}}}$ are shown in Fig. 8. The minimal RMSEs obtained with ${N_{\textrm{new},\textrm{hp}}}$ and ${N_{\textrm{new},\textrm{agg}}}$ are $1.91$ and $2.94,$ respectively and are close to the RMSE value of $2.01$ obtained with ${N_{\textrm{new},\textrm{hc}}}$ in the hexagonal column cases considered in Fig. 7. This suggests that the ${N_{\textrm{opt}}}$ for different shapes may have similar dependences on x and ${m_\textrm{r}}$. Because the term ${c_2}{x^{1/3}}$ is understood to be related to the edge effect, differences in values of ${c_2}$ for different shapes indicates the plausible result that the edge effect depends on particle shape.

 

Fig. 8. (a) Absolute difference ${N_{\textrm{new},\textrm{hp}}} - {N_{\textrm{opt}}}$. The hexagonal plate has aspect ratio $\sigma = 2a/h = 2.0$. (b) Absolute difference ${N_{\textrm{new},\textrm{agg}}} - {N_{\textrm{opt}}}$. The 8 columns aggregate shape is shown in Fig. 2(c). The dashed contour lines in the two panels denote $|{N_{\textrm{new},\textrm{hp}}} - {N_{\textrm{opt}}}|= 2$ and $|{N_{\textrm{new},\textrm{agg}}} - {N_{\textrm{opt}}}|= 2$, respectively.

Download Full Size | PDF

The number of potential particle shapes is unlimited, and we cannot validate the dependence of ${N_{\textrm{opt}}}$ on x and ${m_\textrm{r}}$ for all shapes. We conduct IITM experiments to determine ${N_{\textrm{opt}}}$ for another eight particle shapes for a collection of size parameters and refractive indices. Then, we use Eq. (19) (${N_{\textrm{new},\textrm{hc}}}$) to estimate ${N_{\textrm{opt}}}$ for the eight particle shapes. The differences ${N_{\textrm{opt}}} - {N_{\textrm{new},\textrm{hc}}}$ for the tested cases are listed in Table 2. The shapes include hexagonal column/plate with different aspect ratios, droxtals, hollow columns, 5-plate aggregates, solid/hollow bullets, and hexahedrons. The shapes are shown in Fig. 9. The hexahedron shape has been used to model single-scattering properties of nonspherical dust aerosols [50]. The other shapes have been used to model single-scattering properties of ice crystals [49].

 

Fig. 9. Additional particle shapes considered in the IITM experiments.

Download Full Size | PDF

Table 2. ${N_{\textrm{opt}}} - {N_{\textrm{new},\textrm{hc}}}$ for various particles with other shapes shown in Fig. 9

View Table | View all tables in this article

As shown in Table 2, ${N_{\textrm{opt}}} - {N_{\textrm{new},\textrm{hc}}}$ varies significantly for different shapes. This indicates that the relation between ${N_{\textrm{opt}}}$ and x and ${m_\textrm{r}}$ is highly dependent on particle shape. Although we cannot give a single empirical formula to estimate ${N_{\textrm{opt}}}$ for all the particle shapes, based on our numerical experiments for the three shapes, namely hexagonal column with aspect ratio $0.8$, hexagonal plate with aspect ratio $2.0$, and 8 columns aggregate, the relation between ${N_{\textrm{opt}}}$ and x and ${m_\textrm{r}}$ can largely be described by Eq. (18). The coefficients ${c_1}$, ${c_2}$ and ${c_3}$ are shape-dependent.

5. Conclusion

For a given particle, the dimension of the T-matrix or the maximum VSWF degree, the radial resolution of spherical shells in the IITM recursive computation, and the precision of the ${\mathbf {U}}$ matrix calculation all contribute to the determination of accuracy in IITM computations. By specifying sufficient radial resolution and quadrature points in the ${\mathbf {U}}$ matrix calculation, we focus on determining the T-matrix truncation number. For computation of three single-scattering properties (scattering phase matrix, and extinction and scattering cross sections), we have argued that a criterion of convergence should be based on the scattering phase function, because it will in most cases ensure convergence of other single-scattering properties, and we have introduced a criterion for convergence based on the scattering phase function.

IITM calculations of the sort we have conducted are highly CPU-intensive, and so we were only able to examine a few relatively simple particle shapes that have been commonly studied in the single-scattering literature, and a modest range of choices of the size parameter, the index of refraction, and particle shape. Our results may therefore be better regarded as indications of directions for further research than definitive findings valid for the entire range of particles considered in single scattering studies. One simple conclusion is that a larger T-matrix dimension is needed for nonabsorptive particles with larger refractive index. Beyond that observation, we have investigated how well an empirical formula that has just three adjustable coefficients is able to describe the relation between optimal truncation number (which determines the T-matrix dimension), the size parameter, and the refractive index. Our numerical experiments indicate that all three coefficients are shape-dependent. The empirical formula has been presented as a generalization of earlier formulae for optimal truncation number that ignored the index of refraction. While our results do not suggest that there is a universal empirical formula of so simple a form as we have investigated, the results do clearly illustrate significant sensitivity of optimal truncation number in IITM computations to factors other than the size parameter, in particular to the refractive index and particle shape.