Flexible implementation of the particle shape and internal inhomogeneity in the invariant imbedding T-matrix method

Zheng Wang; Lei Bi; Senyi Kong

doi:10.1364/OE.498190

1. Introduction

The invariant imbedding T-matrix (IITM) is one of the powerful methods for computing the optical properties of dielectric particles with complex shapes and compositions. Although the IITM was first proposed by Johnson [1] in 1988, the numerical capabilities of the IITM have only been recognized since 2013 due to the extensive studies of Bi et al. [2] and Bi and Yang [3] for improving the IITM’s theoretical framework and numerical implementation. Analytical algorithms of random-orientation average [4,5] give the IITM higher computational efficiency than other methods that compute the optical properties of particles at all discrete orientations such as discrete dipole approximation (DDA, [6–8]) and finite-difference time domain (FDTD, [9,10]) methods. Versus other T-matrix methods such as the extended boundary condition method (EBCM, [4,5,11]), the IITM is based on a volume integral equation making itself more stable for computing the optical properties of various non-spherical and inhomogeneous particles. In principle, light scattering by arbitrarily shaped non-spherical and inhomogeneous particle can be solved by the IITM, and the maximum size parameter can be up to geometric optical domains [12,13].

The IITM has been used for computing the optical properties of various non-spherical and inhomogeneous particles. For example, Bi and Yang [14] used the IITM to calculate light scattering by ice crystals including droxtals, hexagonal columns and plates, bullet rosettes, and aggregates of columns or plates; ice particles with a rough surface could also be studied. Super-spheroidal models combining IITM are often used for aerosols [15–20], and some have an inhomogeneous structure [15,17,18,20]. The most complex shape using the IITM may be the coccolith and coccolithophore models [21]. However, these works require careful programming of particle shape in the IITM to ensure the numerical correctness. This procedure is usually tricky. To overcome this inconvenience, Bi and Yang [3] proposed a numerical shape implementation algorithm based on the triangular mesh of particle surface. More specifically, the particle surface is represented by triangular facets. One can then determine whether spherical patches are inside or outside the particle. However, the algorithm based on triangular meshes is difficult to apply to arbitrary inhomogeneous particles. It has required tedious efforts for implementing inhomogeneous particles in the IITM.

This paper introduces a new approach for representing the particle shape and composition in the IITM. Inspired by ADDA [8], in which a particle is discretized in cubical subvolumes (dipoles) and the inhomogeneity of the particle is represented by assigning refractive indices at different volume domains, a particle in the IITM is represented through elements of spherical grids and refractive index is associated with each element. Details of the shape file will be presented in Section 2. Based on this new approach, the level of complexity of representing the particle shape and inhomogeneity in the IITM is similar to that in the DDA and FDTD. This new version of IITM is referred to as IITM-discrete hereafter.

The paper is organized as follows. Section 2 introduces the IITM-discrete algorithm. In Section 3, the IITM-discrete is tested, and numerical results are compared with the original IITM and ADDA. Section 4 summarizes this study.

2. Method

A complete algorithm of the IITM is described in previous studies (Bi et al. [2], Bi and Yang [3], Sun et al. [22]). Here, we focus on the U-matrix, which is the only part of the IITM algorithm that depends on particle shape and the refractive indices. The formula of the U-matrix is shown below [3]:

(1)$${U_{mn,m^{\prime}n^{\prime}}}(r) = \frac{{{k^2}{r^2}}}{{4\pi }}\int_0^\pi {d\theta \sin \theta \int_0^{2\pi } {d\phi \textrm{exp} [ - i(m - m^{\prime})\phi ] \times [\varepsilon (r,\theta ,\phi ) - 1]{K_{mn,m^{\prime}n^{\prime}}}(r,\theta ,\phi )} }$$

(2)$${K_{mn,m^{\prime}n^{\prime}}}(r,\theta ,\phi ) = \left( {\begin{array}{{ccc}} {{{\tilde{\pi }}_{mn}}{{\tilde{\pi }}_{m^{\prime}n^{\prime}}} + {{\tilde{\tau }}_{mn}}{{\tilde{\tau }}_{m^{\prime}n^{\prime}}}}&{ - i({{\tilde{\pi }}_{mn}}{{\tilde{\tau }}_{m^{\prime}n^{\prime}}} + {{\tilde{\tau }}_{mn}}{{\tilde{\pi }}_{m^{\prime}n^{\prime}}})}&0\\ {i({{\tilde{\pi }}_{mn}}{{\tilde{\tau }}_{m^{\prime}n^{\prime}}} + {{\tilde{\tau }}_{mn}}{{\tilde{\pi }}_{m^{\prime}n^{\prime}}})}&{{{\tilde{\pi }}_{mn}}{{\tilde{\pi }}_{m^{\prime}n^{\prime}}} + {{\tilde{\tau }}_{mn}}{{\tilde{\tau }}_{m^{\prime}n^{\prime}}}}&0\\ 0&0&{\frac{{\tilde{d}_{0m}^n\tilde{d}_{0m^{\prime}}^{n^{\prime}}}}{{\varepsilon (r,\theta ,\phi )}}} \end{array}} \right)$$

(3)$${\tilde{\pi }_{mn}}(\theta ) = {\left[ {\frac{{2n + 1}}{{n(n + 1)}}} \right]^{1/2}}\frac{m}{{\sin \theta }}d_{0m}^n(\theta )$$

(4)$${\tilde{\tau }_{mn}}(\theta ) = {\left[ {\frac{{2n + 1}}{{n(n + 1)}}} \right]^{1/2}}\frac{d}{{d\theta }}d_{0m}^n(\theta )$$

(5)$$\tilde{d}_{0m}^n(\theta ) = \sqrt {2n + 1} d_{0m}^n(\theta )$$

where k is the wave number, $(r,\theta ,\phi )$ is the spherical coordinate, m is the projected angular momentum, n is the angular momentum, ɛ is the permittivity, and $d_{0m}^n$ is the Wigner-d function [3]. The U matrix can be analytically obtained for spherical particles. For axially symmetric particles (e.g., spheroids), a semi-analytical algorithm is proposed in Bi and Yang [23]. For an n-fold symmetry particle (e.g., hexagonal ice crystals), the integral can still be simplified for efficient computation. Here, we only consider the universal algorithm that is applicable to arbitrary shapes. In Bi et al. [2], the integral process of Eq. (1) is to integrate along the azimuth angle (ϕ) first and then integrate along the zenith angle (θ) where numerical methods are applied. The procedure is illustrated as follows:

(6)$$\begin{aligned} {U_{mn,m^{\prime}n^{\prime}}}(r) &= \int_0^\pi {d\theta \sin \theta \int_0^{2\pi } {F(\varepsilon ,r,\theta ,\phi )d\phi } } \\ &= \int_0^\pi {G(\varepsilon ,r,\theta )d\theta } = \sum\limits_\theta {G(\varepsilon ,r,\theta )\Delta \theta } \end{aligned}$$

in which, $F(\varepsilon ,r,\theta ,\phi )$ is the integrand of Eq. (1), and $G(\varepsilon ,r,\theta )$ is the integral of $F(\varepsilon ,r,\theta ,\phi )$ along the azimuth angle. Note that $G(\varepsilon ,r,\theta )$ is always the exact solution of the integration in the original code. However, this requirement makes the calculation of $G(\varepsilon ,r,\theta )$ complicated—especially for inhomogeneous particles. To provide a universal solution, we rewrite the integration as:

(7)$${U_{mn,m^{\prime}n^{\prime}}}(r) = \int_\Omega {F(\varepsilon ,r,\theta ,\phi )d\Omega } = \sum\limits_\Omega {F(\varepsilon ,r,\theta ,\phi )\Delta \Omega }$$

where $d\Omega $/$\Delta \Omega $ are the solid angle. Here, the entire integration domain is discretized—not just in the zenith angle dimension. It is clear that Eq. (7) may have more computational costs than Eq. (6), and it usually has a lower accuracy if $\Delta \Omega $ is not small enough. The only advantage of Eq. (7) is that it is quite easy to implement even for very complicated shapes with internal inhomogeneity.

In principle, any spherical grid is adaptable. However, the choice of the discrete spherical grid is important to ensure the numerical accuracy and meanwhile minimize the computational cost. For example, if azimuth and zenith angles are evenly discretized, then the grids are much finer at the poles than in the equator. To achieve the best performance, we tend to adopt a relatively uniform spherical grid provided no additional shape information is known. General circulation models (GCMs) encounter similar problems, and their solutions for the discrete grid are insightful. For example, the cubed sphere grid is adopted by the Finite-Volume Cubed-Sphere Dynamical Core (FV3), which projects the inner cubic grid onto the outer bounding sphere [24]. A twisted icosahedral grid is used for solving shallow-water equations in GCM [25]. Its use is also found in radiative transfer simulations [26]. The Model for Prediction Across Scales (MPAS, [27]) uses centroidal Voronoi tessellation (CVT) and provides a quasi-uniform grid. CVT is a special Voronoi tessellation whose generating points coincide with cell centroids. CVT is useful for many applications such as data compression, clustering, and especially the optimal quadrature rule [28]. Du and Ju [29] illustrated quadratic error when solving the convection-diffusion problem with the CVT grid. The advantages of the CVT grid allow it to solve Eq. (7). We selected Hüttig’s method to create the CVT grid, which is based on the spherical spiral line [30]. The formula of the spherical spiral is shown as follows:

(8)$$\begin{array}{l} x = \cos (\alpha )\cos ( - \frac{\pi }{2} + \frac{\alpha }{{{\alpha _{max}}}}\pi )\\ y = \sin (\alpha )\cos ( - \frac{\pi }{2} + \frac{\alpha }{{{\alpha _{max}}}}\pi )\\ z ={-} \sin ( - \frac{\pi }{2} + \frac{\alpha }{{{\alpha _{max}}}}\pi )\\ {\alpha _{max}} = \frac{{3{\pi ^2}}}{{2 \cdot Resolution}} \end{array}$$

in which, ${\alpha _{max}}$ determines the number of revolutions from north pole to south pole and depends on the grid resolution. Term $\alpha$ ranges from 0 to ${\alpha _{max}}$. The spherical spiral is illustrated in Fig. 1(a). In Hüttig’s method, the generating points are first generated by evenly dividing the spherical spiral line and then adjusted to the centroids iteratively based on Lloyd’s method [31]. The grid resolution is defined as the length of the evenly-divided spherical spiral segments. Figure 1 shows examples of the CVT grid with resolutions of 0.1, 0.05, and 0.03.

Fig. 1. CVT grids with resolutions of 0.1, 0.05, and 0.03 based on Hüttig’s method [30]. The blue line in (a) is a spherical spiral.

Download Full Size | PDF

The IITM-specified shape file could be easily defined based on the CVT grid. According to the IITM formalism, the entire particle can be viewed as an inhomogeneous sphere (the refractive index in the outer medium is different from that in the particle). The inhomogeneous sphere is radially split into spherical layers ranging from the inscribed sphere to the circumscribed sphere. Medium information (an index representing refractive index) on the CVT grid should be provided for each layer. For example, Fig. 2 shows how to label medium information for a super-spheroid in which grids inside the super-spheroid are rendered in red while the others are rendered in white.

Fig. 2. A super-spheroid and its inhomogeneous spherical layers discretized by CVT grid with a resolution of 0.1. Grids inside the super-spheroid are labelled as red, while the others are labelled as white.

Download Full Size | PDF

Cell coordinates (radius, θ and ϕ) of the CVT grid should be given in addition to the medium number that is associated with refractive index. The solid angles of cells are also necessary. The particle geometry and composition are thus approximately encoded in the shape file. The accuracy of the shape file can be adjusted by changing the CVT grid resolution. Importantly, generating the IITM-specified shape file is not much harder than generating the ADDA-specified shape file. It is also possible to convert the ADDA-specified shapefiles to the IITM-specified shapefiles. However, any potential shape distortions resulting from the conversions should be carefully examined. The procedures are summarized as follows:

1. Split the particle into spherical layers ranging from the inscribed sphere to the circumscribed sphere.
2. Choose a discrete spherical grid. The spiral CVT grid is used here, although the resolution is critical (similar to the number of dipoles in the DDA).
3. For each layer, label each cell with a number indicating its composition. An index of zero is reserved for the outer medium and is usually the vacuum.
4. Finally, describe the medium and grid information and record them in a shape file.

Both the CVT grid resolution and the radial step (related to the number of spherical layers) are important for resolving the particle shape and inhomogeneity. In defining shapes, the radial step (Δr_shape) should be sufficiently small to guarantee a fine representation of the particle shape and inhomogeneity. However, in the IITM computational program, the radial step (Δr_code) is a tunable parameter. When one spherical layer in the IITM is involved for calculating the U matrix, the nearest layer in the shape file is identified to feed the medium and grid information. Therefore, the only requirement is to make Δr_shape smaller than Δr_code. Empirically, kΔr_code (k is wavenumber) is fixed to 0.1∼0.2 and kΔr_shape is much smaller than 0.1, which is found to be sufficient for present calculations. Practically, optimal values of kΔr_code and kΔr_shape could be dependent on specific shape and inhomogeneity.

3. Result

The accuracy and performance of IITM-discrete are tested in this section.

3.1 Compared with original IITM

Super-spheroid is one of the typical geometries used for modelling aerosols [15,16]. Here, the single-scattering properties of a super-spheroid (the aspect ratio is 1.0 and the roundness parameter is 2.5) are calculated with the IITM-discrete and the original IITM. The particle size parameter is 10, and the refractive index is 1.5 + 0.001i. Two shape files with grid resolutions of 0.1 (Discrete-0.1) and 0.05 (Discrete-0.05) are tested with the IITM-discrete. The comparison results are shown in Fig. 3. The results of the original IITM are used as a “benchmark”. In general, the IITM-discrete results agree with the original IITM especially for the phase function (P₁₁). However, the Discrete-0.1 test yields noticeable errors in both the extinction/scattering cross section and the phase matrix except P₁₁. However, these errors can be reduced by increasing the grid resolution as shown in the Discrete-0.05 test. For example, the relative error of extinction cross section decreases from 1.4% to 0.4% by comparing a Discrete-0.05 test with a Discrete-0.1 test. Errors of the IITM-discrete majorly come from incomplete sampling of the particle boundaries. In this case, the sharp tips of the super-spheroid are often ignored when using a shape file, while the original IITM would always consider them with careful programming. The price of the improved accuracy in the Discrete-0.05 test is increased computation time. In practice, the accuracy and performance should be balanced for need as discussed in Section 3.3. Note that the application superiority of the IITM-discrete is not a homogeneous particle for which traditional algorithms are more accurate.

Fig. 3. Single-scattering properties of a super-spheroid using the original IITM and the IITM-discrete with grid resolutions of 0.1 and 0.05. The size parameter is 10, and the refractive index is 1.5 + 0.001i.

Download Full Size | PDF

The next example is an inhomogeneous particle, and the shape is a super-spheroid containing several spheres. This was first studied in Zong et al. [17]. The refractive indices of the super-spheroids and spheres are 1.5 + 0.001i and 2.5 + 0.9i. Comparison results are shown in Fig. 4. The IITM-discrete agrees very well with the original IITM. The relative error of the extinction cross section is 0.02% versus the original IITM. This is much lower than the errors in Fig. 3 even if the same grid is used. This is because the IITM-discrete does not take any advantage of particle symmetry as it is designed for complex particles. Symmetric particles, such as super-spheroids having four-fold symmetry along the z-axis, are more sensitive to symmetry-breaking operations than asymmetric particles. This is similar to the case in [32] that the DDA errors for spheres are usually larger than other shapes. Therefore, the original IITM and IITM-discrete have more similar results when calculating particles with no symmetry.

Fig. 4. Single-scattering properties of an inhomogeneous super-spheroid containing spheres using the original IITM and the IITM-discrete with a grid resolution of 0.05. The size parameter is 10, and the refractive index of super-spheroid (spheres) is 1.5 + 0.001i (2.5 + 0.9i).

Download Full Size | PDF

3.2 Compared with ADDA

To illustrate the generality of IITM-discrete, another inhomogeneous super-spheroid was used for testing. The super-spheroid is divided by three planes into seven regions (Fig. 5). The refractive indices of the seven regions are selected causally: 1.3 + 0.0001i, 1.5 + 0.001i, 1.35 + 0.0005i, 1.4 + 0.003i, 1.55 + 0.01i, 1.45 + 0.005i, and 1.6 + 0.01i. For the original IITM, dealing with this shape would be a challenging task. In contrast, the easily generated shape file makes the IITM-discrete quite convenient. In fact, the IITM-specified shape file is similar with the ADDA-specified shape file. They both use 3D pixels to represent the particle shape while the coordinates of pixels are different. The ADDA is based on the Cartesian cubic grid, while the IITM-discrete is based on the spherical grid. Comparison results of the IITM-discrete and ADDA of the inhomogeneous super-spheroid are shown in Fig. 5. The phase matrix of IITM-discrete nicely fits the ADDA.

Fig. 5. Single-scattering properties of an inhomogeneous super-spheroid using the ADDA and the IITM-discrete with a grid resolution of 0.05. The super-spheroid is divided by three planes into seven regions. The size parameter is 10, and the refractive indices of seven regions are 1.3 + 0.0001i, 1.5 + 0.001i, 1.35 + 0.0005i, 1.4 + 0.003i, 1.55 + 0.01i, 1.45 + 0.005i, and 1.6 + 0.01i.

Download Full Size | PDF

The next example is a coated fractal model that is often used for black carbon (BC) aerosols mixed with other aerosols [33,34] as shown in Fig. 6. The monomers of the fractal are overlapped, and the overlap parameter ($1 - d/D$) is 0.2 in which d is the distance of two adjacent monomers, and D is the diameter of monomers. The refractive index of the fractal is 1.95 + 0.79i as indicated by Bond and Bergstrom [35] for BC aerosols. The refractive index of the coating is 1.5 + 0.0001i. The size parameter of the volume-equivalent sphere of the BC core is 5. The relative difference of extinction cross section is 0.4% between the IITM-discrete and the ADDA. The phase matrix of IITM-discrete is also consistent with the ADDA.

Fig. 6. Single-scattering properties of a coated BC model using the ADDA and the IITM-discrete with a grid resolution of 0.05. The size parameter of the volume-equivalent sphere of BC is 5, and the refractive index of BC (coating) is 1.95 + 0.79i (1.5 + 0.0001i).

Download Full Size | PDF

The two examples show the generality of IITM-discrete and guarantees its numerical accuracy with respect to other computational methods.

3.3 Performance and accuracy

The key factor in generating the shape file is the grid resolution. Finer grids usually lead to higher accuracy but with a higher computational cost in the IITM-discrete. In this section, we compare the numerical accuracy of the extinction cross section and computation cost with different grid resolutions as shown in Fig. 7. The particle is the same as that in Fig. 4, and the result of the original IITM works as the “benchmark” and is not explicitly shown. Particles of size parameters of 10 and 20 are considered. When the size parameter is 10, the relative error with a grid resolution of 0.2 is very large—about 18%, which means that the grid is too coarse. The relationship between extinction cross section errors and resolutions satisfies a power law approximately when the resolution is larger than 0.05, but the order is larger than two, indicating that C_ext converges fast. When the grid resolution is sufficiently small, the relative error is nearly unchanged. At the same grid resolution, the relative errors of C_ext at the size parameter 20 is larger than their counterparts at the size parameter of 10. However, in general, the relationship between the relative errors of C_ext and the grid resolution is similar.

Fig. 7. Relative errors of extinction cross section and relative total computation time of the IITM-discrete using different grid resolutions compared to the original IITM. The refractive index of super-spheroid (spheres) is 1.5 + 0.001i (2.5 + 0.9i).

Download Full Size | PDF

The relative computation time with respect to the original IITM increases as the grid resolution decreases. When the grid resolution is larger than 0.05, the total computation time grows relatively slowly as the resolution decreases. However, the total computation time grows quickly in nearly second order with respect to resolutions smaller than 0.05, indicating that the computational time of calculating U-matrix becomes dominant in the IITM.

Figure 7 also shows that the IITM-discrete is faster than the original IITM when the grid resolution is larger than 0.03. One of the reasons is that the grids are same for every layer. This implies that the angular function depending on (θ, ϕ) needed by the U matrix can be calculated only once. The (θ, ϕ) pairs vary in different layers in the original IITM; thus, the angular function must be calculated separately for each layer.

3.4 Importance of inhomogeneity for large dust aerosols

Dust aerosols are composed of different minerals, and are inhomogeneous in reality [36]. Although homogenous non-spherical dust models are extensively studied [16,37,38], several studies [39,40] show the importance of the internal inhomogeneity on dust optical properties. The inhomogeneous super-spheroid model shown in Fig. 5 is used here and represents inhomogeneous dust aerosols. The refractive indices of the seven regions are 1.406 + 0.0015i, 1.502 + 2 × 10⁻⁷i, 1.572 + 1 × 10⁻⁵i, 1.587 + 3 × 10⁻⁵i, 1.431 + 0.002i, 1.509 + 4.6 × 10⁻⁴i, and 2.5 + 0.928i—these are associated with illite, quartz, calcite, feldspar, chlorite, kaolinite, and hematite, respectively. The inhomogeneous dust model is then compared with the homogeneous super-spheroid model. The effective refractive index of the homogeneous model is 1.486 + 0.0034i based on the Bruggeman mixing rule. The phase matrices of the two models are compared in Fig. 8 at size parameters of 30, 50, and 80. The two phase matrices are quite similar when the size parameters are 30 and 50, thus indicating that the inhomogeneity has only minor effects when the particles are relatively small. However, noticeable differences appear at a size parameter of 80 for the P₁₁, P₂₂, P₃₃, P₄₄, and P₄₃ elements. Backscattering intensity is particularly underestimated when the inhomogeneity is not considered. The inhomogeneity may cause larger biases for size parameters larger than 80. This example shows the limitation of effective medium approximation at large size parameters, consistent with previous studies [39,40]. Inhomogeneity should not be ignored. Coarse and giant dust aerosols are commonly found over the Saharan region [41], and the inhomogeneity may be an important factor for lidar observations such as the lidar ratio. Extensive investigation of the impact of internal inhomogeneity on the dust optical properties in a broad range of size parameters will be reported in future.

Fig. 8. Comparison of the phase matrix of inhomogeneous dust (using the IITM-discrete with a resolution of 0.03) and homogeneous dust (using the original IITM) at size parameters of 30, 50, and 80.

Download Full Size | PDF

4. Summary

The IITM is a stable and efficient light scattering method for non-spherical inhomogeneous particles. Although a flexible algorithm based on the triangular mesh has already been introduced in the IITM (Bi and Yang, 2014 [3]), it is only applicable for homogeneous particles. Implementing the inhomogeneous particles in the IITM has required tedious or substantial efforts on programming. The adaptability of IITM is not yet fully explored in practice. In this paper, we developed a volume discretization approach that is specifically designed for an arbitrary and inhomogeneous particle. The shape file stores a collection of discrete spherical layers ranging from the inscribed sphere to the circumscribed sphere, thus describing the refractive index information of grids. The numerical integrals in the U-matrix are greatly simplified with the aid of a shape file, and no additional programming is required to modify the IITM program. Appropriate spherical grids should be selected in consideration of the numerical accuracy and computational efficiency. Practically, we adopted the spherical CVT grid because the grids are nearly uniform on the spherical surface and have relatively low truncation errors.

Numerical tests of the IITM-discrete were compared with other light scattering methods including the original IITM and the ADDA. The results of a few examples, i.e., whether the particle is homogeneous or inhomogeneous, show that the IITM-discrete is reliable. The grid resolution in the IITM-discrete should be appropriately tuned. A coarse grid may lead to inaccurate results. The computational cost increases as the grid resolution is refined. This requires users to select a suitable grid resolution balancing accuracy and efficiency. An empirical rule on grid resolution could be obtained, but it should depend on the level of complexity of particle shape and internal inhomogeneity. This could be the limitations of the IITM-discrete. For example, if particles have very tiny inclusions or extreme aspect ratios, then they require a very small grid resolution to ensure the correct shape. A minimum number of layers in the IITM algorithm may also need to be increased manually. In such cases, programming the particle shape manually in the IITM would be a better choice. Therefore, the IITM-discrete provides a universal solution for the IITM, but it may be not the most computationally efficient option. Further investigation could be necessary to identify an optimized grid scheme, which costs less. Recent studies [42,43] show that spherical design and Lebedev quadrature schemes are promising for optimal spherical integration.

Funding

National Natural Science Foundation of China (42022038); National Key Research and Development Program of China (2022YFC3004004).

Acknowledgments

We acknowledge Ms. Rui Liu at the Training Center of Atmospheric Sciences of Zhejiang University for her help with managing computing resources. We also used the National Key Scientific and Technological Infrastructure project“Earth System Numerical Simulation Facility” (EarthLab), and the computing facilities at China HPC Cloud Computing Center. We also thank M. A. Yurkin and A. G. Hoestra for using their ADDA code.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. B. R. Johnson, “Invariant imbedding T matrix approach to electromagnetic scattering,” Appl. Opt. 27(23), 4861–4873 (1988). [CrossRef]

2. L. Bi, P. Yang, G. W. Kattawar, and M. I. Mishchenko, “Efficient implementation of the invariant imbedding T-matrix method and the separation of variables method applied to large nonspherical inhomogeneous particles,” J. Quant. Spectrosc. Radiat. Transfer 116, 169–183 (2013). [CrossRef]

3. L. Bi and P. Yang, “Accurate simulation of the optical properties of atmospheric ice crystals with the invariant imbedding T-matrix method,” J. Quant. Spectrosc. Radiat. Transfer 138, 17–35 (2014). [CrossRef]

4. M. I. Mishchenko, “Light scattering by randomly oriented axially symmetric particles,” J. Opt. Soc. Am. A 8(6), 871–882 (1991). [CrossRef]

5. M. I. Mishchenko, L. D. Travis, and A. A. Lacis, Scattering, Absorption, and Emission of Light by Small Particles (Cambridge University, 2002).

6. B. T. Draine and P. J. Flatau, “Discrete-Dipole Approximation For Scattering Calculations,” J. Opt. Soc. Am. A 11(4), 1491–1499 (1994). [CrossRef]

7. M. A. Yurkin, V. P. Maltsev, and A. G. Hoekstra, “The discrete dipole approximation for simulation of light scattering by particles much larger than the wavelength,” J. Quant. Spectrosc. Radiat. Transfer 106(1-3), 546–557 (2007). [CrossRef]

8. M. A. Yurkin and A. G. Hoekstra, “The discrete-dipole-approximation code ADDA: Capabilities and known limitations,” J. Quant. Spectrosc. Radiat. Transfer 112(13), 2234–2247 (2011). [CrossRef]

9. K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. 14(3), 302–307 (1966). [CrossRef]

10. P. Yang and K. N. Liou, “Finite-difference time domain method for light scattering by small ice crystals in three-dimensional space,” J. Opt. Soc. Am. A 13(10), 2072–2085 (1996). [CrossRef]

11. M. I. Mishchenko and L. D. Travis, “Capabilities and limitations of a current FORTRAN implementation of the T-matrix method for randomly oriented, rotationally symmetric scatterers,” J. Quant. Spectrosc. Radiat. Transfer 60(3), 309–324 (1998). [CrossRef]

12. S. Zhai, R. L. Panetta, and P. Yang, “Improvements in the computational efficiency and convergence of the Invariant Imbedding T-matrix method for spheroids and hexagonal prisms,” Opt. Express 27(20), A1441–A1457 (2019). [CrossRef]

13. P. Yang, J. Ding, R. L. Panetta, K.-N. Liou, G. W. Kattawar, and M. Mishchenko, “On the Convergence of Numerical Computations for Both Exact and Approximate Solutions for Electromagnetic Scattering by Nonspherical Dielectric Particles,” Prog. Electromagn. Res. 164, 27–61 (2019). [CrossRef]

14. L. Bi and P. Yang, “Improved ice particle optical property simulations in the ultraviolet to far-infrared regime,” J. Quant. Spectrosc. Radiat. Transfer 189, 228–237 (2017). [CrossRef]

15. L. Bi, W. Lin, Z. Wang, X. Tang, X. Zhang, and B. Yi, “Optical Modeling of Sea Salt Aerosols: The Effects of Nonsphericity and Inhomogeneity,” J. Geophys. Res.: Atmos. 123(1), 543–558 (2018). [CrossRef]

16. W. Lin, L. Bi, and O. Dubovik, “Assessing Superspheroids in Modeling the Scattering Matrices of Dust Aerosols,” J. Geophys. Res.: Atmos. 123(24), 917–943 (2018). [CrossRef]

17. R. Zong, F. Weng, L. Bi, X. Lin, C. Rao, and W. Li, “Impact of hematite on dust absorption at wavelengths ranging from 0.2 to 1.0 µm: an evaluation of literature data using the T-matrix method,” Opt. Express 29(11), 17405–17427 (2021). [CrossRef]

18. L. Bi, Z. Wang, W. Han, W. Li, and X. Zhang, “Computation of Optical Properties of Core-Shell Super-Spheroids Using a GPU Implementation of the Invariant Imbedding T-Matrix Method,” Front. Remote Sens. 3, 903312 (2022). [CrossRef]

19. W. Lin, L. Bi, F. Weng, Z. Li, and O. Dubovik, “Capability of Superspheroids for Modeling PARASOL Observations Under Dusty-Sky Conditions,” J. Geophys. Res.: Atmos. 126(1), e2020JD033310 (2021). [CrossRef]

20. Z. Wang, L. Bi, H. Wang, Y. Wang, W. Han, and X. Zhang, “Evaluation of a new internally-mixed aerosol optics scheme in the weather research and forecasting model,” J. Quant. Spectrosc. Radiat. Transfer 283, 108147 (2022). [CrossRef]

21. L. Bi and P. Yang, “Impact of calcification state on the inherent optical properties of Emiliania huxleyi coccoliths and coccolithophores,” J. Quant. Spectrosc. Radiat. Transfer 155, 10–21 (2015). [CrossRef]

22. B. Sun, L. Bi, P. Yang, M. Kahnert, and G. Kattawar, Invariant Imbedding T-Matrix Method for Light Scattering by Nonspherical and Inhomogeneous Particles (Elsevier, 2020).

23. L. Bi and P. Yang, “Tunneling effects in electromagnetic wave scattering by nonspherical particles: A comparison of the Debye series and physical-geometric optics approximations,” J. Quant. Spectrosc. Radiat. Transfer 178, 93–107 (2016). [CrossRef]

24. W. M. Putman and S.-J. Lin, “Finite-volume transport on various cubed-sphere grids,” J. Comput. Phys. 227(1), 55–78 (2007). [CrossRef]

25. R. Heikes and D. A. Randall, “Numerical Integration of the Shallow-Water Equations on a Twisted Icosahedral Grid. Part I: Basic Design and Results of Tests,” Mon. Weather Rev. 123(6), 1862–1880 (1995). [CrossRef]

26. C. Wang, P. Yang, S. L. Nasiri, S. Platnick, B. A. Baum, A. K. Heidinger, and X. Liu, “A fast radiative transfer model for visible through shortwave infrared spectral reflectances in clear and cloudy atmospheres,” J. Quant. Spectrosc. Radiat. Transfer 116, 122–131 (2013). [CrossRef]

27. W. C. Skamarock, J. B. Klemp, M. G. Duda, L. D. Fowler, S.-H. Park, and T. D. Ringler, “A Multiscale Nonhydrostatic Atmospheric Model Using Centroidal Voronoi Tesselations and C-Grid Staggering,” Mon. Weather Rev. 140(9), 3090–3105 (2012). [CrossRef]

28. Q. Du, V. Faber, and M. Gunzburger, “Centroidal Voronoi Tessellations: Applications and Algorithms,” SIAM Rev. 41(4), 637–676 (1999). [CrossRef]

29. Q. Du and L. Ju, “Finite Volume Methods on Spheres and Spherical Centroidal Voronoi Meshes,” SIAM J. Numer. Anal. 43(4), 1673–1692 (2005). [CrossRef]

30. C. Hüttig and K. Stemmer, “The spiral grid: A new approach to discretize the sphere and its application to mantle convection,” Geochem., Geophys., Geosyst. 9(2), 1 (2008). [CrossRef]

31. S. Lloyd, “Least squares quantization in PCM,” IEEE Trans. Inf. Theory 28(2), 129–137 (1982). [CrossRef]

32. M. A. Yurkin, M. Min, and A. G. Hoekstra, “Application of the discrete dipole approximation to very large refractive indices: Filtered coupled dipoles revived,” Phys. Rev. E 82(3), 036703 (2010). [CrossRef]

33. J. Luo, Q. Zhang, J. Luo, J. Liu, Y. Huo, and Y. Zhang, “Optical Modeling of Black Carbon With Different Coating Materials: The Effect of Coating Configurations,” J. Geophys. Res.: Atmos. 124(23), 13230–13253 (2019). [CrossRef]

34. Y. Wang, Y. Pang, J. Huang, L. Bi, H. Che, X. Zhang, and W. Li, “Constructing Shapes and Mixing Structures of Black Carbon Particles With Applications to Optical Calculations,” J. Geophys. Res.: Atmos. 126(10), e2021JD034620 (2021). [CrossRef]

35. T. C. Bond and R. W. Bergstrom, “Light Absorption by Carbonaceous Particles: An Investigative Review,” Aerosol Sci. Technol. 40(1), 27–67 (2006). [CrossRef]

36. G. Y. Jeong and T. Nousiainen, “TEM analysis of the internal structures and mineralogy of Asian dust particles and the implications for optical modeling,” Atmos. Chem. Phys. 14(14), 7233–7254 (2014). [CrossRef]

37. M. Saito, P. Yang, J. Ding, and X. Liu, “A Comprehensive Database of the Optical Properties of Irregular Aerosol Particles for Radiative Transfer Simulations,” J. Atmos. Sci. 78(7), 2089–2111 (2021). [CrossRef]

38. O. Dubovik, A. Sinyuk, T. Lapyonok, B. N. Holben, M. Mishchenko, P. Yang, T. F. Eck, H. Volten, O. Muñoz, B. Veihelmann, W. J. van der Zande, J.-F. Leon, M. Sorokin, and I. Slutsker, “Application of spheroid models to account for aerosol particle nonsphericity in remote sensing of desert dust,” J. Geophys. Res.: Atmos. 111(D11), D11208 (2006). [CrossRef]

39. M. Kahnert, “Modelling radiometric properties of inhomogeneous mineral dust particles: Applicability and limitations of effective medium theories,” J. Quant. Spectrosc. Radiat. Transfer 152, 16–27 (2015). [CrossRef]

40. J.-M. Perrin and P. L. Lamy, “On the Validity of Effective-Medium Theories in the Case of Light Extinction by Inhomogeneous Dust Particles,” ApJ 364, 146 (1990). [CrossRef]

41. C. L. Ryder, E. J. Highwood, A. Walser, P. Seibert, A. Philipp, and B. Weinzierl, “Coarse and giant particles are ubiquitous in Saharan dust export regions and are radiatively significant over the Sahara,” Atmos. Chem. Phys. 19(24), 15353–15376 (2019). [CrossRef]

42. M. A. Yurkin, “Fair Evaluation of Orientation-Averaging Techniques in Light-Scattering Simulations: Comment on “Evaluation of Higher-Order Quadrature Schemes in Improving Computational Efficiency for Orientation-Averaged Single-Scattering Properties of Nonspherical Ice Particles” by Fenni et al,” JGR Atmospheres 128(2), e2021JD036088 (2023). [CrossRef]

43. I. Fenni, K.-S. Kuo, M. S. Haynes, Z. S. Haddad, and H. Roussel, “Evaluation of Higher-Order Quadrature Schemes in Improving Computational Efficiency for Orientation-Averaged Single-Scattering Properties of Nonspherical Ice Particles,” J. Geophys. Res.: Atmos. 126(11), e2020JD034172 (2021). [CrossRef]

Flexible implementation of the particle shape and internal inhomogeneity in the invariant imbedding T-matrix method

Abstract

1. Introduction

2. Method

3. Result

3.1 Compared with original IITM

3.2 Compared with ADDA

3.3 Performance and accuracy

3.4 Importance of inhomogeneity for large dust aerosols

4. Summary

Funding

Acknowledgments

Disclosures

Data availability

References

Data availability

Cited By

Figures (8)

Equations (8)

Optics Express