Design and Validation of the Invariant Imbedded T-Matrix Scattering Model for Atmospheric Particles with Arbitrary Shapes

: Light scattering by non-spherical particles is an important factor inﬂuencing atmospheric radiative transfer. To accurately simulate the scattering properties of non-spherical particles, the Invariant Imbedded T-matrix method (IIM T-Matrix) is developed by combining the Lorenz–Mie theory and invariant imbedding technique. In this model, the non-spherical particle is regarded as an inhomogeneous sphere and discretized into multiple spherical layers in the spherical coordinate system. The T-matrix of the inscribed sphere is ﬁrstly calculated by the Lorenz–Mie theory, and then taking it as the initial value, the T-matrix is updated layer by layer by using the invariant imbedding technique. To improve the computational e ﬃ ciency, the model is further parallelized by the OpenMP technique. To verify the simulation accuracy of the IIM T-Matrix method, the results of the model are compared with those of the EBCM (Extended Boundary Condition Method) T-Matrix method, DDA (Discrete Dipole Approximation) and MRTD (Multi-Resolution Time Domain). The results show that the scattering phase matrix simulated by the IIM T-Matrix method closely agrees with that of the well-tested models, indicating that the IIM T-Matrix method is a powerful tool for the light scattering simulation of non-spherical particles. Since the IIM T-Matrix method is derived from the volume integral equation, compared to the T-Matrix method which is based on surface integral principles (i.e., “EBCM” or the “null ﬁeld method”), it can be applied to the scattering calculations of particle with arbitrary shapes and inhomogeneous compositions, which can greatly expand the application scope of the T-Matrix method. IIM T-Matrix method was compared with that of EBCM, DDASCAT, and MRTD. The results show that good agreements were achieved between the results of the IIM T-Matrix method and those well-tested scattering models, indicating that the IIM T-Matrix method can accurately simulate the scattering properties of non-spherical particles with di ﬀ erent shapes and inhomogeneous compositions. After being improved by the iterative acceleration scheme based on the Lorenz–Mie theory, the computational e ﬃ ciency of the IIM T-Matrix method was improved notably. Owning to its excellent performance, the IIM T-Matrix model might become a powerful tool for the light scattering simulation of non-spherical particles in the atmospheric radiation ﬁeld.


Introduction
The light scattering and absorption by non-spherical particles (such as ice crystals, dust aerosols, etc.) plays an important role in the atmospheric radiative transfer [1][2][3][4]. However, because of the irregular shape of these particles, there is still considerable uncertainty in their scattering characteristics [5][6][7], which has become an important factor restricting the accuracy of atmospheric remote sensing [8,9]. For polarized remote sensing, the effect of the "non-spherical effect" is much more remarkable since the polarization of diffuse light is very sensitive to the particles' shape [10,11]. Therefore, how to obtain the scattering properties of these non-spherical particles accurately has become a hot issue in the field of atmospheric radiation [9,[12][13][14].
In order to simulate the light scattering by non-spherical particles, many scattering computational models were established [15,16]. Typical models include the T-matrix method [17,18], Discrete Dipole introduced in Section 2.3. Further, the acceleration scheme based on the Lorenz-Mie theory is proposed in Section 2.4.

A Brief Introduction of the T-Matrix
The T-matrix is a linear transformation matrix of the expansion coefficients of incident and scattered fields. In order to obtain the T-matrix, the incident and scattered fields should first be expanded by vector spherical harmonic functions [16,29], shown as: n m=−n a mn RgM mn (kr, θ, ϕ) + b mn RgN mn (kr, θ, ϕ) (1) n m=−n p mn M mn (kr, θ, ϕ) + q mn N mn (kr, θ, ϕ), (r > R) (2) where, R is the radius of the smallest circumscribed sphere of the scatterer centered at the origin of the laboratory coordinate system. RgM mn (kr, θ, ϕ) and RgN mn (kr, θ, ϕ) are the regular vector spherical wave functions, M mn (kr, θ, ϕ) and N mn (kr, θ, ϕ) are the vector spherical wave functions,p mn and q mn are the expansion coefficients of the scattering electric field, and a mn and b mn are the expansion coefficients of the incident field. The definition of the vector spherical wave functions can be found in the Appendix A of this paper.
Since the Maxwell equations are linear, the expansion coefficients of incident light (a mn and b mn ) and scattered light (p mn and q mn ) can also be converted by a linear transformation. The linear transformation matrix is just the T-matrix, which can be expressed as: In order to express the T-matrix more intuitively and compactly, the two formulas above can be written into matrix form, given as: in which, the T-matrix is a 2l max × 2l max super-matrix, where, l max = n max (n max + 2) (n max is the highest expansion order of the field components).

Discretization and Vectorization of the Helmholtz Volume Integral Equation
The physical basis of the IIM T-matrix method is the Helmholtz volume integral equation. If we take the time dependence of the electromagnetic wave as exp(−jωt), then the equation in frequency domain can be expressed as [40]: Appl. Sci. 2019, 9, 4423 4 of 22 in which, E denotes the electric vector, u(r) = k 2 (ε r (r) − 1), k is the wavenumber, r is the position vector,ε r (r) is the relative permittivity of the medium, E inc (r) is the electric vector of incident wave, G 0 (r, r ) denotes the free space dyadic Green's function, and Z(r) is a transformation matrix determined by the particle's permittivity, the expressions of G 0 (r, r ) and Z(r) are shown in Equations (7) and (8).
In the equations above, the operator "⊗" denotes the dyad product, i.e., if a and b are two column vectors, then a ⊗ b = a · b T . In order to realize the discretization of Equation (6), both field components and dyadic Green's function G 0 (r, r ) need to be expanded by vector spherical harmonic functions. Therefore, the expansion and vectorization of the field components and dyad Green's function are firstly introduced, and then the discretization scheme of volume integral equation is further derived.
a. Expansion and vectorization of the field components and dyad Green's function G 0 (r, r ). The expansion expression of the incident field has been given by Equation (1), which can be written into matrix form, given as: where, Y mn (θ, ϕ) is an angular function matrix, and J n (kr) is the radial Bessel function matrix, which can be written as: where, d n 0m (θ) is the Wigner-d function, while π mn (θ) and τ mn (θ) are the angular functions, given by: Similarly, the total field can also be expanded by the vector spherical harmonic functions and written into matrix form, expressed as: where, E mn is the expansion coefficients of the total field, and H n (kr) is the radial Hankel function matrix, whose form is similar to J n (kr), expressed as: Appl. Sci. 2019, 9, 4423 5 of 22 Similar to the field components, the dyadic Green's function G 0 (r, r ) should also be expanded by the vector spherical harmonic functions, but in the case that r > r' and r < r', their expansion forms are slightly different, written as: Because the expansion of the dyadic Green's function is similar for r > r' and r < r' on the whole, here we only introduce the matrix expansion process of G 0 (r, r ) for r > r'. From the symmetry of the dyadic Green's function, the following equation can be easily derived, given as: By using the symmetric relations of the vector spherical harmonic function [16], Equation (13) can be simplified and rewritten it into the matrix form, expressed as: where, g n (r, r ) is a matrix constructed by the Bessel and Hankel functions, written as: Similarly, for r < r , the dyad Green's function G 0 (r, r ) can also be written into the matrix form, similar to Equation (15), written as: On the condition that r = r , the expansion coefficients of G 0 (r, r ) are taken as the average value of those for r > r and r < r . In this case, the matrix form of the dyad Green's function is similar to that of Equation (15), but the matrix g n (r, r ) should be rewritten as: g n (r, r ) = ik J n (kr)H n (kr ) + H n (kr)J n (kr ) /2 (19) b. The discretization of the volume integral equation. Substitute the expansion equations of the field components and the dyadic Green's function into Equation (6), and use the orthogonality of the vector spherical harmonic function for simplification, then, the following integral formulas can be obtained for each (m, n), written as: E m n (r, θ, ϕ) = Y m n (θ, ϕ)J n (kr) Appl. Sci. 2019, 9, 4423 6 of 22 For further simplification, we define the vector amplitude density function F mnm n (r) as: T u(r)Z(r)E m n (r)dΩ (21) This parameter is related to the electric current density and can be thought of as the source of the scattered electric wave. Then, by substituting the equation above into Equation (20), we can get: Further, we insert Equation (22) into Equation (21), then the following relation can be obtained after some simplification, given as: n m=− n U mnm n (r)g n (r, r )F m nm n (r ) (23) where, U mnm n (r) is the surface integral matrix of the spherical shell, written as: From the formula above, we can find that the U-matrix contains the spatial distribution information of the scatterer, that is, the irregular shape of particles and its material distribution can be manifested by this matrix, so its calculation accuracy is very important in the IIM T-Matrix model.
Since the scattering field usually refers to the electromagnetic field at infinite distance, in this case, matrix g n has the form of g n (r, r ) = ikH(kr)J n T n (kr ), since r is larger than r' in the far field. Then, substitute matrix g n into Equation (26), and we can obtain: where, T mnm n is the T-matrix introduced in Section 2.1, which can be calculated by: Equations (23), (25), and (26) constitute the basic equation set for the calculation of the T-matrix. By solving these equations, the T-matrix can be obtained directly, and the scattering parameters can also be calculated based on the T-matrix elements.
In order to solve the equations above, the radial integrals in Equations (25) and (26) should be discretized by Gaussian quadrature, which yields: n m=− n U mnm n (r i )g n (r i , r j )F m nm n (r j ) (27) T mnm n = ik N j=1 ω j J T n (kr j )F mnm n (r j ) (28) where, r j ( j = 1, 2 . . . , N) are the Gaussian quadrature points, and ω j is the weight factor for each Gaussian point. Further, the above two formulas can be written into the matrix form, expressed as: where, J(kr i ) and g(r i , r j ) are the super-matrices with J n (kr i ) and g n (r, r ) as the element, and = F(r i ) and U(r i ) are the super-matrices comprised by F mnm n and U mnm n . The expression of T matrix is shown in Equation (5).

T-Matrix Computation Based on the Invariant Embedding Technique
In principle, Equation (29) is essentially a linear equation system, and can be solved by standard numerical methods, such as the Gaussian elimination method, the conjugate gradient method, etc., after F(r i ) is obtained, and then, the T-matrix can be obtained easily by substituting F(r i ) into Equation (30). However, because the matrices involved in this model are all super-matrices, the implementation of this scheme is very difficult, especially for large particles. Therefore, in order to improve the computational efficiency of the T-matrix, Johnson [33] first introduced the invariant embedding technique into the light scattering simulation (this technique was originally developed to solve quantum mechanical scattering problems). This method not only has the advantage of high computational efficiency, but can also calculate the T-matrix directly without calculating the matrix F(r i ), so it can greatly reduce computational time and memory consumption.
According to Equations (29) and (30), the non-spherical particle should be discretized along the radial direction, which is equivalent to dividing the particle into multiple inhomogeneous spherical layers in the spherical coordinate system, as shown in Figure 1. In this way, the non-spherical particle can be viewed as an inhomogeneous sphere, i.e., a portion of the sphere has the dielectric properties of the scattering particle and the rest is regarded as a vacuum. , which satisfies the following equations: n T Figure 1. The discretization of the non-spherical particle in the spherical coordinate system. Assuming n is an integer in (1, N) (N is the number of the Gaussian integral points), then according to Equations (29) and (30), we define the matrix function F(n|r i ) and T(r n ), which satisfies the following equations: where, F(n|r i ) and T(r n ) can be visually regarded as the vector amplitude density function and T-matrix of the sphere with nth spherical shell. From the Equations (31) and (32), it can be found that except for the fact that N is replaced by n, these equations are identical to Equations (29) and (30). Therefore, if we let N = n, we can obtain F(r i ) = F(N|r i ) and T = T(r N ).
Then, separate the nth term in the summation of Equation (31), and define Q(r n ) and q as follows: Then, Equation (31) can be simplified as: where, H(kr) is a diagonal super-matrix with its element defined by Equation (12), ω n is the weight factor of the nth Gaussian point. Because r n is larger than r j in Equation (35), g n (r, r ) has the expression of g n (r, r ) = ikH(kr)J T (kr ). Equation (32) can be reshaped in a similar way. The nth term in the summation is firstly separated, together with Equation (35), then it can be simplified as: where, matrix Q 11 (r n ) and Q 12 (r n ) are two auxiliary matrices, defined as: From the equations above, the relationship between matrix F(n|r n ), T(r n ) and q is established.
However, q is still a super-matrix, therefore, we should further try to eliminate the matrix F(n|r n ) and q to establish the iterative equation only containing the T(r n ) matrix. First, from Equation (32), we can know that: Assuming that F(n|r j ) = F(n − 1|r j )(I + p) (p is a matrix), then from Equation (34), it can rewritten as: At the same time, we insert F(n|r j ) = F(n − 1|r j )(I + p) and g n (r i , into Equation (31), a simplified expression between p and q can be obtained, written as: where, Q 21 (r n ) and Q 22 (r n ) are another two auxiliary matrices, which are defined as follows: Then, by inserting Equation (41) into Equation (40), we can obtain: The equation above can be regarded as a linear matrix equation of q, through which the matrix q can be solved, written as: At last, by substituting Equation (45) into Equation (36), the invariant imbedding iterative equation for the T-matrix can be obtained, given as: From this equation, It can be seen that if the T-matrix T(r n−1 ) of the sphere with (n − 1) layers is known, then the T-matrix T(r n ) of the sphere of nth spherical layers can be modified on the basis of T(r n−1 ), and the correction terms mainly depend on the optical properties of the nth spherical layer (through the matrix Q ij (i, j = 1,2)).

Iterative Acceleration Based on the Lorenz-Mie Theory
According to Equation (46), during the implementation of the invariant embedding T-matrix method, the non-spherical particles are needed to be regarded as an inhomogeneous sphere and divided into a series of spherical shells. In this process, the optical matrices of spherical shells (including U, Q and Q ij (i, j = 1, 2)) should be calculated first, then, the T-matrix of the sphere with n layers can be calculated based on the T-matrix of the smaller one with n − 1 layers.
In the calculation of the T-matrix, if we do the iteration directly from the origin of the sphere (r = 0) to the circumscribed spherical shell (r = R), the computational amount will be very large, which will limit the calculation efficiency of the model. In order to solve this problem, the scattering calculation region of particles is divided into two parts (as shown in Figure 2): one is the spherical core region with a radius of R 0 (R 0 is the radius of the particle's inscribed sphere), the other is the spherical mantle region in R 0 < r < R (the region between the inscribed and circumscribed sphere). The spherical core region is a homogeneous sphere, and its T-matrix can be calculated by the Lorenz-Mie theory. While in the spherical mantle region, the invariant embedding technique is applied to solve the T-matrix, and its initial iteration value is set as the T-matrix of the spherical core. In this way, the unnecessary iteration process can be reduced, and the computational efficiency can be improved as well.
spherical mantle region in R0 < r < R (the region between the inscribed and circumscribed sphere). The spherical core region is a homogeneous sphere, and its T-matrix can be calculated by the Lorenz-Mie theory. While in the spherical mantle region, the invariant embedding technique is applied to solve the T-matrix, and its initial iteration value is set as the T-matrix of the spherical core. In this way, the unnecessary iteration process can be reduced, and the computational efficiency can be improved as well. A specific interface is needed to convert the Mie coefficients into the T-matrix of the inscribed sphere. For spheres, its T-matrix is a diagonal matrix, which satisfies A specific interface is needed to convert the Mie coefficients into the T-matrix of the inscribed sphere. For spheres, its T-matrix is a diagonal matrix, which satisfies T ij mnm n = δ m,m δ n,n · T ij mnm n , so the relationship between the T and Mie coefficients (i.e., a n and b n ) can be expressed as: in which, a n and b n can be calculated by: where, x is the size parameter of the inscribed sphere, and ψ n (x) and ξ n (x) are the Riccati-Bessel functions. Based on the principle introduced above, the implementation of the IIM T-Matrix method can be concluded as follows: Step 1: Determine the inscribed and circumscribed sphere of the non-spherical particle, and calculate the T-matrix of the inscribed sphere with the Lorenz-Mie Theory by using Equations (47)-(50).
Step 2: Taking the T-matrix of the inscribed sphere as the initial value, then the T-matrix of non-spherical particles is updated layer by layer by using Equation (46). The computational flowchart of the optical matrix of the spherical shell is presented in Figure 3. Firstly, the U-matrix is calculated by Equation (24), and then, the Q-matrix is directly calculated from the U-matrix by Equation (33). After the Q-matrix is obtained, Q 11 , Q 12 , Q 21 , and Q 22 can be computed by Equations (37), (38), (42) and (43).
Step 3: After the T-matrix of the non-spherical particles is obtained, then the cross-sections can be directly linked to the elements of the T-matrix, and the phase matrix elements can also be obtained by combining the generalized spherical functions and T-matrix elements, readers can see Reference [16] for the details.
Step 3: After the T-matrix of the non-spherical particles is obtained, then the cross-sections can be directly linked to the elements of the T-matrix, and the phase matrix elements can also be obtained by combining the generalized spherical functions and T-matrix elements, readers can see Reference [16] for the details.

Model Validation and Results Analysis
The invariant imbedded T-matrix model is implemented in Fortran95, and to improve the modeling efficiency, the model is further parallelized by the OpenMP technique (OpenMP is a multithreading implementation technique that allows the compiler to generate code for task and data parallelism). To validate the simulation accuracy of this model, the scattering parameters calculated by the IIM T-Matrix method are compared with those calculated by EBCM and DDASCAT (which have been well tested and are usually taken as a benchmark). The non-spherical particles simulated in this paper include spheroid with different sizes, cylinder, spheroid with a spherical core (inhomogeneous particle), and hexagonal prism (particles with non-rotationally geometry).

Small Spheroidal Particle Case
In this simulation, the incident light wavelength is taken as λ = 0.5 µm (the typical wavelength in the visible band), the refractive index is set as m = 1.60 − 0.008i (the typical value of mineral dust aerosol), the half-length of the horizontal and rotational axis are taken as a = 1.0 µ m and b = 0.5 µ m,

Model Validation and Results Analysis
The invariant imbedded T-matrix model is implemented in Fortran95, and to improve the modeling efficiency, the model is further parallelized by the OpenMP technique (OpenMP is a multi-threading implementation technique that allows the compiler to generate code for task and data parallelism). To validate the simulation accuracy of this model, the scattering parameters calculated by the IIM T-Matrix method are compared with those calculated by EBCM and DDASCAT (which have been well tested and are usually taken as a benchmark). The non-spherical particles simulated in this paper include spheroid with different sizes, cylinder, spheroid with a spherical core (inhomogeneous particle), and hexagonal prism (particles with non-rotationally geometry).

Small Spheroidal Particle Case
In this simulation, the incident light wavelength is taken as λ = 0.5 µm (the typical wavelength in the visible band), the refractive index is set as m = 1.60 − 0.008i (the typical value of mineral dust aerosol), the half-length of the horizontal and rotational axis are taken as a = 1.0 µm and b = 0.5 µm, and the scattering phase matrix of the spheroid is simulated by the IIM T-Matrix method and the EBCM T-Matrix method, respectively. The results are shown in Figure 4.
As shown in the figure, the results obtained by the two models are in good agreement, which verifies the simulation accuracy of the IIM T-Matrix model for the non-spherical particles with small size. For the phase function F 11 , the relative errors are less than 5% in forward scattering directions. Though the calculation errors are relatively large near 180 • , they are still smaller than 12%. For F 12 /F 11 , the calculation errors of the IIM T-Matrix model are all within 0.1, while the simulation errors of F 34 /F 11 and F 44 /F 11 are less than 0.1 in most scattering angles. From the spatial distribution of simulation errors, it can be found that the simulation accuracy of the IIM T-Matrix model in the forward scattering direction is obviously higher than that in large scattering angles. The reason is that, in the IIM T-Matrix model, the particle is discretized in a spherical coordinate system, and similar to PSTD and MRTD, there are also stepped approximation errors in shape construction. Further, because the scattering phase matrix is more sensitive to the particle's shape in backscattering directions, larger simulation errors will be caused in large scattering angles correspondingly.   To further validate the simulation accuracy of the IIM T-Matrix method, the integral scattering parameters (including extinction cross-section C ext , absorption cross-section C abs , scattering cross-section C sca , and single scattering albedo ω) are also compared with those obtained by EBCM, the results are shown in Table 1. It can be found that the relative simulation errors of the integral scattering parameters are generally smaller than 1%, indicating that the integral scattering parameters can be calculated by the IIM T-Matrix method with a high accuracy.

Large Spheroidal Particle Case
To validate the calculation precision of the IIM T-Matrix method for particles with large sizes, light scattering by a large spheroid particle is simulated by the IIM T-Matrix method and EBCM, respectively, and the result is shown in Figure 5. In this simulation, the light wavelength is set as λ = 0.6 µm, the refractive index of the particle is taken to be m = 1.60 − 0.0008i, and the length of the horizontal and rotational axes are set as a = 4.5 µm and b = 6.0 µm, respectively. The region between the inscribed and circumscribed spheres is divided into 50 layers along the radial direction.
From the figures, it can be found that the scattering phase matrix obtained by the IIM T-Matrix model shows a high consistency with that obtained by EBCM, indicating that the IIM T-Matrix method can simulate the large non-spherical particle effectively. For F 11 , the relative simulation errors are all less than 10% in the scattering angles ranging from 0 • to 120 • . Though the simulation errors are slightly larger in the scattering angles in large scattering angles, its maximum error is smaller than 15%. For F 12 /F 11 and F 44 /F 11 , the calculation errors of the IIM-T-matrix model are less than 0.05 in most scattering directions, especially in scattering angles ranging from 0 • to 60 • , the modeling accuracy is much higher than other scattering directions. The simulation accuracy of F 34 /F 11 is slightly lower than that of F 12 /F 11 and F 44 /F 11 , but their simulation errors are still less than 0.1 at most scattering angles. Similar to the small spheroidal particle, the simulation accuracy is lower in the large scattering angles, while in the forward scattering directions, the simulation errors are much smaller.

Cylindrical Particle Case
The modeling capability of the IIM T-Matrix model is investigated for cylindrical particles. In this test, the light wavelength is taken as 0.5 µm, the refractive index of the particle is set as m = 1.60 − 0.0008i, and the diameter (D) and length (L) of the cylinder are set as 2 µm. The scattering phase matrices are calculated by the EBCM and IIM T-Matrix methods, and the results are shown in Figure  6. Similar to Section 3.1, the integral scattering parameters of the large spheroid of our model are also compared with those of EBCM, the results are presented in Table 2. As can be seen, a good consistency is achieved between the results obtained by different models. The relative differences of these parameters are all less than 0.5%, which indicates that the IIM T-Matrix method can simulate the scattering process of large particles accurately.

Cylindrical Particle Case
The modeling capability of the IIM T-Matrix model is investigated for cylindrical particles. In this test, the light wavelength is taken as 0.5 µm, the refractive index of the particle is set as m = 1.60 − 0.0008i, and the diameter (D) and length (L) of the cylinder are set as 2 µm. The scattering phase matrices are calculated by the EBCM and IIM T-Matrix methods, and the results are shown in Figure 6.  It can be seen that the results calculated by the IIM T-Matrix model show a good agreement with those of EBCM, indicating that the IIM T-Matrix method can simulate the light scattering process of a cylinder with a high accuracy. For F 11 , its relative errors are less than 10% in most scattering directions, and for F 12 /F 11 , F 34 /F 11 , and F 44 /F 11 , their absolute errors are generally smaller than 0.1. Similar to the results of spheroidal particles, the spatial distribution of simulation errors is similar to the variation pattern of the phase matrix elements; namely, in the scattering directions where the curves of the phase matrix elements change dramatically, their simulation errors will increase as well.
The integral scattering parameters of different models are also calculated, as shown in Table 3. It can be found that the results of the IIM T-Matrix method show a good agreement with those of EBCM, which validates the modeling accuracy of the IIM T-Matrix method.

Inhomogeneous Particle Case
To verify the modeling capability of the IIM T-Matrix method for particles with inhomogeneous compositions, the results obtained by the IIM T-matrix model are compared with the DDASCAT for a spheroid with a spherical core. In this case, the wavelength of incident light is taken as 0.5 µm, the horizontal and rotational axis of the outer spheroid are set as a = 1.0 µm and b = 1.5 µm respectively, and its complex refractive index is set to be m 1 = 1.44 − 0.000i. The inner part is a spherical core with a radius of 1.0 µm, its refractive index is taken as m 2 = 1.20 − 0.000i. The region between the inscribed and circumscribed spheres is divided into 30 layers along the radial direction. The simulation results are shown in Figure 7.
As can be seen from the figure, the IIM-T-matrix model achieves high calculation accuracy, and the curves of the scattering phase matrix are almost coincided with those of DDASCAT. For F 11 , the relative simulation errors are all less than 10% in scattering angles ranging from 0 • to 30 • . With the increase of scattering angle, its simulation accuracy decreases, but the relative errors are within 26%. The spatial distributions of the simulation errors of F 12 /F 11 , F 34 /F 11 , and F 44 /F 11 are similar to F 11 . In the forward scattering direction, their simulation accuracy is higher and the absolute simulation errors tend to be 0, while as the scattering angle becomes larger, their simulation accuracy is slightly reduced, but it should be noted that their absolute simulation errors are still within 0.1 in most scattering directions. From the discussion, it can be found that the IIM T-Matrix method can simulate the light scattering by inhomogeneous particles with a high accuracy. a spheroid with a spherical core. In this case, the wavelength of incident light is taken as 0.5 µm, the horizontal and rotational axis of the outer spheroid are set as a = 1.0 µm and b = 1.5 µm respectively, and its complex refractive index is set to be m1 = 1.44 − 0.000i. The inner part is a spherical core with a radius of 1.0 µm, its refractive index is taken as m2 = 1.20 − 0.000i. The region between the inscribed and circumscribed spheres is divided into 30 layers along the radial direction. The simulation results are shown in Figure 7. As can be seen from the figure, the IIM-T-matrix model achieves high calculation accuracy, and the curves of the scattering phase matrix are almost coincided with those of DDASCAT. For F11, the relative simulation errors are all less than 10% in scattering angles ranging from 0° to 30°. With the increase of scattering angle, its simulation accuracy decreases, but the relative errors are within 26%. The spatial distributions of the simulation errors of F12/F11, F34/F11, and F44/F11 are similar to F11. In the forward scattering direction, their simulation accuracy is higher and the absolute simulation errors tend to be 0, while as the scattering angle becomes larger, their simulation accuracy is slightly reduced, but it should be noted that their absolute simulation errors are still within 0.1 in most

Hexagonal Prism Case
To validate the simulation accuracy of the IIM T-Matrix method for non-rotational symmetric particles (which cannot be effectively simulated by the EBCM T-Matrix method), the results of the IIM T-Matrix method are further compared with that of the MRTD model developed by our team [27,41].
The MRTD scattering model is established based on the Multi-Resolution Time Domain technique. By using this model, light scattering by arbitrarily shaped particles or ice crystals can be effectively simulated. In this model, the plane wave is introduced into the computational space by the TF/SF (Total Field/Scattering Field) technique, the Convolution Perfectly Matched Layer (CPML) is applied to truncate the computational domain, and the volume integral method is employed to transform the near electric field to the far field.
In this simulation, the light wavelength and refractive index are set as λ = 0.5 µm and m = 1.20 − 0.0008i, and the length and bottom edge length of hexagonal prism are set as L = 2.0 µm and a = 1.0 µm, respectively. The results are shown in Figure 8. As can be seen, a high agreement is achieved between the results of the T-matrix method and those of the MRTD model, which indicates that both MRTD and IIM-T-matrix models can effectively simulate the scattering characteristics of a hexagonal prism. Similar to the conclusion drawn above, the consistency of the two models in forward scattering directions is higher than those at large scattering angles, where for the scattering phase functions F11, the relative simulation errors tend to be zero in scattering angles smaller than 30°, while in scattering directions near 170°, the simulation accuracy is relatively lower, and the maximum simulation errors can reach 36.7%.

Analysis of Modeling Efficiency
To further validate the effectiveness of the iterative acceleration scheme based on the Lorenz-Mie theory, the computational efficiency of the improved model is investigated as well. In those simulations, the light wavelength is set as 0.5 µm, and the scatterers are set as spheroidal particles with different shapes. All the scattering processes are simulated on the same computer (32 bit 3.1 GHz), and the computational time of the traditional and improved IIM T-Matrix code is presented in Table 4.  As can be seen, a high agreement is achieved between the results of the T-matrix method and those of the MRTD model, which indicates that both MRTD and IIM-T-matrix models can effectively simulate the scattering characteristics of a hexagonal prism. Similar to the conclusion drawn above, the consistency of the two models in forward scattering directions is higher than those at large scattering angles, where for the scattering phase functions F 11 , the relative simulation errors tend to be zero in scattering angles smaller than 30 • , while in scattering directions near 170 • , the simulation accuracy is relatively lower, and the maximum simulation errors can reach 36.7%.

Analysis of Modeling Efficiency
To further validate the effectiveness of the iterative acceleration scheme based on the Lorenz-Mie theory, the computational efficiency of the improved model is investigated as well. In those simulations, the light wavelength is set as 0.5 µm, and the scatterers are set as spheroidal particles with different shapes. All the scattering processes are simulated on the same computer (32 bit 3.1 GHz), and the computational time of the traditional and improved IIM T-Matrix code is presented in Table 4. As can been seen, after being improved by the iterative acceleration scheme, the computational efficiency of the IIM T-Matrix method is improved notably, where for the particle with (a, b) = (1.5, 2.0), the computational time is cut down by 75%. With the increase of the aspect ratio (a/b), the improvement of the computational efficiency is much more remarkable.

Conclusions
The T-matrix is a complete scattering information set that only depends on the particle size, refractive index, and shape; therefore, once the T-matrix is obtained, the scattering properties of the particles with arbitrary orientations can be calculated analytically, which is a tremendous advantage over DDA, FDTD, and other scattering models. In this paper, a T-matrix scattering model was developed by combining the invariant imbedding iterative technique and the Lorenz-Mie theory. Compared with the T-Matrix model developed by Mishchenko, it can be applied to the light scattering simulation of particles with arbitrary shape and inhomogeneous compositions. In this paper, the basic principle of the IIM T-Matrix method was firstly derived, and then, an iterative acceleration scheme was proposed based on the Lorenz-Mie theory and the invariant imbedding technique. To improve the computational efficiency, the model was further parallelized by the OpenMP technique. To validate the calculation accuracy of the model, the scattering phase matrix computed by the IIM T-Matrix method was compared with that of EBCM, DDASCAT, and MRTD. The results show that good agreements were achieved between the results of the IIM T-Matrix method and those well-tested scattering models, indicating that the IIM T-Matrix method can accurately simulate the scattering properties of non-spherical particles with different shapes and inhomogeneous compositions. After being improved by the iterative acceleration scheme based on the Lorenz-Mie theory, the computational efficiency of the IIM T-Matrix method was improved notably. Owning to its excellent performance, the IIM T-Matrix model might become a powerful tool for the light scattering simulation of non-spherical particles in the atmospheric radiation field.

Conflicts of Interest:
The authors declare no conflict of interest.