The Computational Optimization of the Invariant Imbedding T Matrix Method for the Particles with N-Fold Symmetry

: The invariant imbedding T-matrix (IIM T-matrix) model is regarded as one of the most promising models for calculating the scattering parameters of non-spherical particles. However, the IIM T-matrix model needs to be iterated along the radial direction when calculating the T-matrix, which involves complex calculations such as matrix inversion and multiplication. Therefore, how to improve its computational efﬁciency is an important problem to be solved. Focused on particles with N-fold symmetric geometry, this paper deduced the symmetry in the calculation process of the IIM T-matrix model, derived the block iteration scheme of the T-matrix, and contracted the IIM T-matrix program for particles with N-fold symmetric geometry. Discrete Dipole Approximation (DDA) and Geometrical Optics Approximation (IGOA) were employed to verify the accuracy of the improved IIM T-matrix model. The results show that the six phase matrix elements (P 11 , P 12 /P 11 , P 22 /P 11 , P 33 /P 11 , P 34 /P 11 and P 44 /P 11 ) calculated by our model are in good agreement with other models. The computational efﬁciency of the improved IIM T-matrix model was further investigated. As demonstrated by the results, the computational efﬁciency for the particles with N-fold symmetry improved by nearly 70% with the improvement of the symmetry of U matrix and T matrix. In conclusion, the improved model can remarkably reduce the calculation time while maintaining high accuracy.


Introduction
Atmospheric radiative transfer models have been widely used in the fields of atmospheric remote sensing [1] and climate numerical simulation. The accuracy of radiative transfer simulation primarily depends on the scattering parameters of atmospheric particles, which is the basic input of atmospheric radiative transfer model [2]. Previously, the particles' scattering properties were generally calculated using the Lorenz-Mie theory. However, the scattering parameters of the Lorenz-Mie theory may seriously deviate from the actual light scattering characteristics due to the irregular shape of atmospheric particles, especially for the polarization of the scattered light [3,4]. In recent years, polarized remote sensing gradually developed, and accurately understanding the influence of particle shape on the polarization characteristics of the atmospheric radiation has become a key step to improving the remote sensing of aerosols and ice clouds [4]. Therefore, the calculation of the scattering parameters of the non-spherical particles has become an international research hotspot [5][6][7][8].
The T-matrix method is one of the important methods to calculate the atmospheric nonspherical particles' light scattering properties. Once the T matrix is obtained in a one-time calculation, the scattering properties of the randomly oriented particles can be analytically calculated. Waterman first proposed using the extended boundary method (EBCM) to calculate the T matrix in 1965 [9]. In the 1990s, Mishchenko improved the extended boundary reduce the matrix calculation and memory consumption. The idea should be generalized and we applied it to accelerate the IIM T-matrix method in this study. In this paper, the symmetry of the U matrix was derived and the iteration of T matrix was also simplified for N-fold symmetric geometrical particles. This work can greatly improve the calculation efficiency of particles with N-fold symmetric geometry and provide an effective calculation scheme to construct the particle scattering parameter database.
This paper is divided into five sections. In Section 2, the basic principle of the invariant imbedding T-matrix method is introduced. In Section 3, the symmetry of the U matrix is systematically derived, and the block iterative method for solving the T matrix is proposed subsequently. In Section 4, in order to verify the accuracy of the model, the calculation results are compared with those of discrete dipole approximation, geometric optics approximation and other models. Moreover, the calculation efficiency of the model is also tested. In Section 5, a discussion of the results is presented. Lastly, a summary of this paper is provided in Section 6.
The invariant imbedding T-matrix model calculates the T matrix by solving the Helmholtz volume integral equation in the frequency domain. After expansion and discretization of the spherical harmonic vector wave function, its form can be converted as follows: E m n (r, θ, ϕ) = Y m n (θ, ϕ)J n (kr) + R 0 dr ∞ ∑ n=1 n ∑ m=−n Y mn (θ, ϕ)H n (kr)T mnm n (2) In which, J n (kr) is a large matrix composed of Bessel functions of the first kind, and H n (kr) is a large matrix composed of Hankel functions of the first kind. Y nm (θ, ϕ) is the angle function matrix Y mn (θ, ϕ) = [A mn (θ, ϕ), B mn (θ, ϕ), C mn (θ, ϕ)]; A mn (θ, ϕ) = ( 2n+1 4πn(n+1) ) 1/2 (−1) m exp(imϕ)[θiπ mn (θ) −φτ mn (θ)]; B mn (θ, ϕ) = ( 2n+1 4πn(n+1) ) 1/2 (−1) m exp(imϕ)[θiτ mn (θ) −φπ mn (θ)]; C mn (θ, ϕ) = ( 2n+1 4πn(n+1) ) 1/2 (−1) m exp(imϕ)[ √ n(n+1) kr d n om (θ)r]; (3) In which, d n 0m represents Wigner-d functions, and π mn (θ) and τ mn (θ) are the angular functions derived from Wigner-d functions. These three can be combined into regular vector spherical harmonics and vector spherical harmonics, and then the specific calculation form of the T matrix can be derived: In the formula above, ω n is Gaussian weight and F mnm n (r j ) is calculated as follows: in the formula above, g n (r i , r j ) = ikH(r n )J T (r n ), i is the unit complex number, the U matrix is defined as By using the invariant imbedding technique, the iterative equation can be derived from Equations (4)-(6), written as [22] T(r n ) = Q 11 (r n ) + (I + Q 12 where, r n is the radius of the nth layer, T(r n ) and T(r n−1 ) are the T matrix of the spheres of n layers and n − 1 layers, Q 11 (r n ), Q 12 (r n ), Q 21 (r n ) and Q 22 (r n ) are the optical matrix of the nth spherical shell, which can be calculated using the Q(r n ) matrix, written as The calculation formula of the Q(r n ) matrix is as follows: In the actual calculation process of the IIM T-matrix method, the non-spherical particles are regarded as an inhomogeneous sphere and discretized into a certain number of inhomogeneous spherical shells from R 0 to R N in the spherical coordinate system (R 1 and R 2 are the radii of the inscribed and circumscribed sphere of the particle), as shown in Figure 1. To improve the modeling efficiency, the T matrix of the inscribed sphere is firstly calculated by employing the Lorenz-Mie theory, and then, the T matrix is updated layer by layer using the invariant imbedding technique (see Equation (7)).
In the T-matrix calculation process, the U matrix is the most important matrix because it contains optical and geometrical information about particles, such as particle shape, particle complex refractive index, etc. Therefore, the calculation of the U matrix is an important step in the calculation. However, its calculation is a highly time-consuming process in the IIM T-matrix method. Therefore, in the third section, the optimization theory for calculating the U matrix is firstly derived. process in the IIM T-matrix method. Therefore, in the third section, the optimization theory for calculating the U matrix is firstly derived.

Optimization Method for Calculating Scattering Parameters of Particles with N-Fold Symmetric Geometry
The particles with N-fold symmetric geometry (After the particle rotates a fixed angle around an axis perpendicular to the center of the particle plane, it completely coincides with the original particle) are generally taken as the approximation of non-spherical particles, such as a hexagonal prism, disk-shaped particles, super ellipsoid, etc. [24,29], as shown in Figure 2. In this section, we focus on deriving the symmetrical properties of the U matrix, and on improving the implementation of the IIM T-matrix method.

Computational Optimization Based on Symmetry Property of U Matrix
The calculation formula of the U matrix can be written as: Figure 1. Cross section diagram of non-spherical particles. In the figure, R1 is the maximum inscribed sphere radius and R2 is the minimum circumscribed sphere radius.

Optimization Method for Calculating Scattering Parameters of Particles with N-Fold Symmetric Geometry
The particles with N-fold symmetric geometry (After the particle rotates a fixed angle around an axis perpendicular to the center of the particle plane, it completely coincides with the original particle) are generally taken as the approximation of non-spherical particles, such as a hexagonal prism, disk-shaped particles, super ellipsoid, etc., [24,29], as shown in Figure 2. In this section, we focus on deriving the symmetrical properties of the U matrix, and on improving the implementation of the IIM T-matrix method. process in the IIM T-matrix method. Therefore, in the third section, the optimization theory for calculating the U matrix is firstly derived.

Optimization Method for Calculating Scattering Parameters of Particles with N-Fold Symmetric Geometry
The particles with N-fold symmetric geometry (After the particle rotates a fixed angle around an axis perpendicular to the center of the particle plane, it completely coincides with the original particle) are generally taken as the approximation of non-spherical particles, such as a hexagonal prism, disk-shaped particles, super ellipsoid, etc. [24,29], as shown in Figure 2. In this section, we focus on deriving the symmetrical properties of the U matrix, and on improving the implementation of the IIM T-matrix method.

Computational Optimization Based on Symmetry Property of U Matrix
The calculation formula of the U matrix can be written as:

Computational Optimization Based on Symmetry Property of U Matrix
The calculation formula of the U matrix can be written as: For the particles with N-fold symmetric geometry, we can extract the integral of ϕ as follows: ϕ 1 and ϕ 2 in the above formula are shown in Figure 3.
( 1) ( 1) 0 0 0 0 ( , , ) n n n n n n d d m m r r For the particles with N-fold symmetric geometry, we can extract the integral of ϕ as follows:  ϕ π ϕ π ϕ π ϕ π ϕ π ϕ ϕ ε θ ϕ ϕ ϕ ε θ ϕ ϕ ϕ ε θ ϕ ϕ ϕ ε θ ϕ ϕ ϕ ε θ ϕ ϕ 1 ϕ and 2 ϕ in the above formula are shown in Figure 3. Through simplification, Equation (14) can be expressed as follows (see Appendix A for detailed derivation): Through simplification, Equation (14) can be expressed as follows (see Appendix A for detailed derivation): From the equation above, it can be found that for the particles with N-fold symmetric geometry, the calculation of the U matrix only needs to be carried out at |m − m |= 0/Nl , which greatly reduces the calculation amount.

Calculation Optimization of T-Matrix Iteration Process
It has been found that when |m − m | = 0/Nl , the U matrix is 0. In order to reduce the computational complexity of matrix multiplication, we rearranged the U matrix in the following way: extracting the sub-U matrix of |m − m |= 0/Nl and arranging it along the diagonal of the U matrix, as shown in Figure 4.
Similarly, the Q matrix can also be transformed into a block diagonal matrix. According to the IIM T-matrix theory, the Q matrix is calculated using the U matrix with the following equation: Because the matrix element of g(r n , r n ) is independent of the azimuthal mode m [20,22,31], we can rearrange the layout of the matrix element of g(r n , r n ) into an N-block diagonal matrix according to the arrangement rule of the U matrix [20,22], as shown in Figure 5.
From the equation above, it can be found that for the particles with N-fold symmetric geometry, the calculation of the U matrix only needs to be carried out at , which greatly reduces the calculation amount.

Calculation Optimization of T-Matrix Iteration Process
It has been found that when | | 0/ m m Nl ′ − ≠ , the U matrix is 0. In order to reduce the computational complexity of matrix multiplication, we rearranged the U matrix in the following way: extracting the sub-U matrix of | | 0/ m m Nl ′ − = and arranging it along the diagonal of the U matrix, as shown in Figure 4. Similarly, the Q matrix can also be transformed into a block diagonal matrix. According to the IIM T-matrix theory, the Q matrix is calculated using the U matrix with the following equation: Because the matrix element of ( , ) n n r r g is independent of the azimuthal mode m [20,22,31], we can rearrange the layout of the matrix element of ( , ) n n r r g into an N-block diagonal matrix according to the arrangement rule of the U matrix [20,22], as shown in Figure 5.  Similarly, the Q matrix can also be transformed into a block diagonal matrix. According to the IIM T-matrix theory, the Q matrix is calculated using the U matrix with the following equation: Because the matrix element of ( , ) n n r r g is independent of the azimuthal mode m [20,22,31], we can rearrange the layout of the matrix element of ( , ) n n r r g into an N-block diagonal matrix according to the arrangement rule of the U matrix [20,22], as shown in Figure 5. It is rearranged as a block diagonal matrix with N blocks referring to the structure of the U matrix. Different colors correspond to different azimuth mode "m", brown indicates m = i, red and green indicate increase and decrease on this basis respectively.
Obviously, after g(r n , r n ) is transformed into a block diagonal matrix, it can be easily found that the Q matrix should also be a block diagonal matrix (with N blocks). The ith sub-matrix block, i.e., Q i , can be calculated as follows: where w represents the wth block of the diagonal matrix. Since the Bessel matrix and the Hankel matrix are also independent of the azimuth mode m, they can also be converted to the form of a block diagonal matrix. For this, the calculation method of each block of Q ij (r n ) can be found as follows: Remote Sens. 2022, 14, 4061 So far, all the matrices in Equation (7) become a block diagonal matrix. So, the whole calculation process can be converted to the calculation of the block diagonal matrix, which is written as (22) After all the block diagonal sub-matrices are obtained, they can be combined into a complete T-matrix. Thus, the scattering parameters can calculated in the same manner as that in the EBCM T-matrix method developed by Mishchenko et al.

Results
To validate the effectiveness of our model, hexagonal prism, octagonal prism and twelve prism particles were selected as examples for verification. The calculation results of their scattering parameters were compared with those of DDA (small particles case) and IGOA (large particle case). Further, the computational efficiency of the optimized IIM T-matrix model was also investigated.

Small Particle Case
In this section, we compare the numerical results computed by the IIM T-matrix method with their counterparts computed using DDA and ADDA, where the particles used for model validation include the hexagonal prism particles, octagonal prism particles and twelve prism particles.
(1) Hexagonal prism case In the comparison of hexagonal prism particles, the wavelength of incident light was set as λ = 500 nm, the refractive index of the particles was taken as m = 1.33 + 0.0008 i, the bottom side length was set as a = 1.0 um and the height was set to be L = 2.0 um. The scattering phase matrix of the randomly oriented particles was simulated by IIM T-matrix and DDA methods, as illustrated in Figure 6. In the DDA model calculation, 2548 directions were selected to obtain the random-orientation averaged results.
From the figure, it can be found that the simulation results of the IIM T-matrix model and the DDA model show a good consistency, where the P 11 curves of the phase matrix elements obtained by the IIM T-matrix model and the DDA model almost coincided with each other. The curves of P 12 /P 11 basically overlapped, and there was only a small error at 140 • . When the scattering angle was 80 •~1 10 • , the P 33 /P 11 and P 34 /P 11 curves of the IIM T-matrix model were slightly lower than those of DDA model. The P 22 /P 11 and P 44 /P 11 curves of the IIM T-matrix model were generally lower than those of DDA.
The integral scattering properties were also compared in this paper, and the results are shown in Table 1. As can be seen, the relative errors of C ext (extinction section) and (scattering cross section) are less than 0.2%. The relative errors of C abs (absorption cross section) are a little larger, reaching 7%.  The integral scattering properties were also compared in this paper, and the results are shown in Table 1. As can be seen, the relative errors of Cext (extinction section) and (scattering cross section) are less than 0.2%. The relative errors of Cabs (absorption cross section) are a little larger, reaching 7%. The improved IIM T-matrix model was also validated for a hexagonal prism with a large refractive index. In this comparison, the refractive index of the particles was taken as m = 1.2762 + 0.4133 i and other parameters were unchanged. In the calculation of the The improved IIM T-matrix model was also validated for a hexagonal prism with a large refractive index. In this comparison, the refractive index of the particles was taken as m = 1.2762 + 0.4133 i and other parameters were unchanged. In the calculation of the ADDA model, 4097 orientation directions were used to obtain the average results (as illustrated in Figure 7). Additionally, the integral scattering properties were also compared, as shown in Table 2. The results show that there is a good consistency between two models. Especially in forward scattering directions, the deviations between the two models are much smaller. For P 12 /P 11 , in the scattering angle of 60 •~9 0 • , the IIM T-matrix model calculation result was slightly higher. The curves of P 44 /P 11 , P 33 /P 11 and P 34 /P 11 are basically consistent. The P 22 /P 11 curve of the IIM T-matrix model was generally higher than that of ADDA. For the integral scattering parameters, the relative errors of C ext were less than 0.1% and the relative errors of C sca and C abs were less than 0.6%. basically consistent. The P22/P11 curve of the IIM T-matrix model was generally higher than that of ADDA. For the integral scattering parameters, the relative errors of Cext were less than 0.1% and the relative errors of Csca and Cabs were less than 0.6%.  In the comparison of octagonal prism particles, the wavelength of incident light was set as λ = 1000 nm, the refractive index of the particles was taken as m = 1.33 + 0.0008 i, the bottom side length as set as a = 1.0 um and the height was set as L = 2.0 um. The scattering phase matrix of the randomly oriented particles was simulated using the IIM T-matrix method and DDA model, as illustrated in Figure 8. In DDA model calculation, 1690 orientation directions were selected to obtain the average results. (2) Octagonal prism case In the comparison of octagonal prism particles, the wavelength of incident light was set as λ = 1000 nm, the refractive index of the particles was taken as m = 1.33 + 0.0008 i, the bottom side length as set as a = 1.0 um and the height was set as L = 2.0 um. The scattering phase matrix of the randomly oriented particles was simulated using the IIM T-matrix method and DDA model, as illustrated in Figure 8. In DDA model calculation, 1690 orientation directions were selected to obtain the average results.
From the graph, it can be found that better agreement was achieved between the scattering phase matrix calculated by DDA and IIM T-matrix method. There was basically no difference between the P 11 curves. For P 12 /P 11 , P 33 /P 11 and P 34 /P 11 , the variation patterns of the deviations between the two models were similar to that of P 11 , and the absolute differences between the two models were within 0.05. The calculated results of small octagonal particles show higher consistency than those of the small hexagonal prism particle. For P 22 /P 11 and P 44 /P 11 , in the scattering angle of 0 •~1 40 • , the absolute errors were within 0.02. The integral scattering properties are given in Table 3. The relative errors of C ext and C abs were less than 0.1% and the relative errors of C sca were less than 0.5%.  From the graph, it can be found that better agreement was achieved between the scattering phase matrix calculated by DDA and IIM T-matrix method. There was basically no difference between the P11 curves. For P12/P11, P33/P11 and P34/P11, the variation patterns of the deviations between the two models were similar to that of P11, and the absolute differences between the two models were within 0.05. The calculated results of small octagonal particles show higher consistency than those of the small hexagonal prism particle. For P22/P11 and P44/P11, in the scattering angle of 0°~140°, the absolute errors were within 0.02. The integral scattering properties are given in Table 3. The relative errors of Cext and Cabs were less than 0.1% and the relative errors of Csca were less than 0.5%. In the comparison of the twelve prism particles, the wavelength of incident light was set as λ = 1000 nm, the refractive index of the particles was taken as m = 1.33 + 0.0008 i, the bottom side length was set as a = 0.5 um and the height was L = 2.0 um. The scattering In the comparison of the twelve prism particles, the wavelength of incident light was set as λ = 1000 nm, the refractive index of the particles was taken as m = 1.33 + 0.0008 i, the bottom side length was set as a = 0.5 um and the height was L = 2.0 um. The scattering phase matrix of the randomly oriented particles was simulated using IIM T-matrix and DDA methods, as illustrated in Figure 9. In DDA model calculation, 1960 orientation directions were selected to obtain the average results.
The results show that good agreement was achieved between the IIM T-matrix model and DDA model, and their P 11 curves completely coincided with each other. Additionally, the simulation errors of P 12 /P 11 , P 33 /P 11 and P 34 /P 11 were less than 0.02. It can be found that the calculation accuracy of twelve prism particles is higher than that of hexagonal and octagonal particles. For P 22 /P 11 , the absolute errors are generally within 0.01. For P 44 /P 11 , in the scattering angle of 0 •~1 50 • , the absolute errors are within 0.02. The integral scattering properties are given in Table 4. The relative errors of C ext and C sca are less than 0.06%. The relative errors of C abs are less than 0.3%.  The results show that good agreement was achieved between the IIM T-matrix model and DDA model, and their P11 curves completely coincided with each other. Additionally, the simulation errors of P12/P11, P33/P11 and P34/P11 were less than 0.02. It can be found that the calculation accuracy of twelve prism particles is higher than that of hexagonal and octagonal particles. For P22/P11, the absolute errors are generally within 0.01. For P44/P11, in the scattering angle of 0°~150°, the absolute errors are within 0.02. The integral scattering properties are given in Table 4. The relative errors of Cext and Csca are less than 0.06%. The relative errors of Cabs are less than 0.3%.

) The comparison of the improved and traditional T-matrix method
The differences in the calculation results of the T-matrix model before and after optimization were also compared. In this comparison, the particles used for model validation were hexagonal prism particles. The wavelength of incident light was set as λ = 500 nm, the refractive index of the particles was taken as m = 1.33 + 0.0008 i, the bottom side length was set as a = 1.0 um and the height was set to be L = 2.0 um. The scattering Figure 9. The phase matrices and their errors calculated by the IIM T−matrix and DDA methods for small twelve prism particles with random orientations. The abscissa is the scattering angle. The solid line is the IIM T−matrix model and the dashed is DDA model. The errors of P 12 /P 11 , P 22 /P 11 , P 33 /P 11 , P 34 /P 11 and P 44 /P 11 are absolute errors.

(4) The comparison of the improved and traditional T-matrix method
The differences in the calculation results of the T-matrix model before and after optimization were also compared. In this comparison, the particles used for model validation were hexagonal prism particles. The wavelength of incident light was set as λ = 500 nm, the refractive index of the particles was taken as m = 1.33 + 0.0008 i, the bottom side length was set as a = 1.0 um and the height was set to be L = 2.0 um. The scattering phase matrix of the randomly oriented particles was simulated by implementing IIM T-matrix methods before and after optimization, as illustrated in Figure 10, and the integral scattering parameters are presented in Table 5. The results show that the consistency between the IIM T-matrix models before and after optimization was high. The relative errors of P 11 were less than 0.1%. Additionally, the absolute errors of P 12 /P 11 , P 33 /P 11 and P 34 /P 11 were less than 0.03. For P 22 /P 11 , in the scattering angle of 0 •~1 40 • , the absolute errors were generally within 0.01. For P 44 /P 11 , in the scattering angle of 0 •~1 40 • , the absolute errors were within 0.05. The integral scattering properties also demonstrated good agreement, and the relative errors of C ext , C sca and C abs were all less than 2%.  The results show that the consistency between the IIM T-matrix models before and after optimization was high. The relative errors of P11 were less than 0.1%. Additionally, the absolute errors of P12/P11, P33/P11 and P34/P11 were less than 0.03. For P22/P11, in the scattering angle of 0°~140°, the absolute errors were generally within 0.01. For P44/P11, in the scattering angle of 0°~140°, the absolute errors were within 0.05. The integral scattering properties also demonstrated good agreement, and the relative errors of Cext, Csca and Cabs were all less than 2%.

Large Particle Case
In this section, we compare the numerical results computed from the IIM T-matrix Figure 10. The phase matrices and their errors calculated by the IIM T−matrix models before and after optimization of small hexagonal prism particles with random orientations. The abscissa is the scattering angle. The solid line is the IIM T−matrix model before optimization and the dashed is the IIM T−matrix model after optimization. The errors of P 12 /P 11 , P 22 /P 11 , P 33 /P 11 , P 34 /P 11 and P 44 /P 11 are absolute errors.

Large Particle Case
In this section, we compare the numerical results computed from the IIM T-matrix model and their counterparts computed using IGOA, where the particles used for model validation include the hexagonal prism particles and twelve prism particles.
In the comparison of hexagonal prism particles, the refractive index of the particles was taken as m = 1.3078 + 1.66 × 10 −8 i, the size parameters were, respectively, set at 80 and 120 (2πH/λ = 80, 120), and the ratio of the maximum distance between any two points on the base to the height of the prism was 1:1. (as illustrated in Figure 11). The halos emerged when the size parameter was 80, when the size parameter was 120, and the ice halo peaks become sharp. This proves that the IIM T-matrix model can well reflect the diffraction properties of large hexagonal prism particles. For P 12 /P 11 , P 33 /P 11 , P 22 /P 11 , P 44 /P 11 and P 34 /P 11 , good agreement was also achieved between the results of the IIM T-matrix method and IGOA. The P 12 /P 11 curves were basically coincident, but when the scattering angle was 60 •~1 20 • , for P 33 /P 11 and P 34 /P 11 curves, the calculation result of the IIM T-matrix model was slightly lower than that of IGOA. halo peaks become sharp. This proves that the IIM T-matrix model can well reflect the diffraction properties of large hexagonal prism particles. For P12/P11, P33/P11, P22/P11, P44/P11 and P34/P11, good agreement was also achieved between the results of the IIM T-matrix method and IGOA. The P12/P11 curves were basically coincident, but when the scattering angle was 60°~120°, for P33/P11 and P34/P11 curves, the calculation result of the IIM T-matrix model was slightly lower than that of IGOA. In the comparison of twelve prism particles, the refractive index of the particles was taken as m = 1.3078 + 1.66 × 10 −8 i, the size parameters were, respectively .set at 80 and 120 (2 / 80,120) H π λ = , and the ratio of the maximum distance between any two points on the base to the height of the prism was 1:1 (as illustrated in Figure 12). From the results, it can be found that a good consistency was achieved between the results of IIM T-matrix method and IGOA. For P12/P11 and P22/P11, the two models displayed good consistency. For P33/P11, P44/P11 and P34/P11, in forward scattering directions, the deviations between the two models were small, but with the increased scattering angles, their differences increased gradually. Figure 11. Comparison of four phase matrix elements computed from the IGOA and the II−TM for hexagonal prism particles. The particle size parameters for curve pairs from lower to upper are 80, 120. For clarity, the phase functions of 120 are multiplied by 1000, P 12 /P 11 , P 22 /P 11 , P 33 /P 11 , P 44 /P 11 and P 34 /P 11 of 120 are shifted by 1. The solid line is the IIM T-matrix model and the dashed is IGOA model.
In the comparison of twelve prism particles, the refractive index of the particles was taken as m = 1.3078 + 1.66 × 10 −8 i, the size parameters were, respectively .set at 80 and 120 (2πH/λ = 80, 120), and the ratio of the maximum distance between any two points on the base to the height of the prism was 1:1 (as illustrated in Figure 12). From the results, it can be found that a good consistency was achieved between the results of IIM T-matrix method and IGOA. For P 12 /P 11 and P 22 /P 11 , the two models displayed good consistency. For P 33 /P 11 , P 44 /P 11 and P 34 /P 11 , in forward scattering directions, the deviations between the two models were small, but with the increased scattering angles, their differences increased gradually.
halo peaks become sharp. This proves that the IIM T-matrix model can well reflect the diffraction properties of large hexagonal prism particles. For P12/P11, P33/P11, P22/P11, P44/P11 and P34/P11, good agreement was also achieved between the results of the IIM T-matrix method and IGOA. The P12/P11 curves were basically coincident, but when the scattering angle was 60°~120°, for P33/P11 and P34/P11 curves, the calculation result of the IIM T-matrix model was slightly lower than that of IGOA. In the comparison of twelve prism particles, the refractive index of the particles was taken as m = 1.3078 + 1.66 × 10 −8 i, the size parameters were, respectively .set at 80 and 120 (2 / 80,120) H π λ = , and the ratio of the maximum distance between any two points on the base to the height of the prism was 1:1 (as illustrated in Figure 12). From the results, it can be found that a good consistency was achieved between the results of IIM T-matrix method and IGOA. For P12/P11 and P22/P11, the two models displayed good consistency. For P33/P11, P44/P11 and P34/P11, in forward scattering directions, the deviations between the two models were small, but with the increased scattering angles, their differences increased gradually. Figure 12. Comparison of four phase matrix elements computed from the IGOA and the II-TM for twelve prism particles. The particle size parameters for curve pairs from lower to upper are 80, 120. For clarity, the phase functions of 120 are multiplied by 1000, P 12 /P 11 , P 22 /P 11 , P 33 /P 11 , P 44 /P 11 and P 34 /P 11 of 120 are shifted by 1. The solid line is the IIM T-matrix model and the dashed is IGOA model.

Calculation Efficiency by IIM T-Matrix Model before and after Optimization
To validate the effectiveness of our improvements, the computational efficiency of the IIM T-matrix method was also evaluated. In this simulation, the refractive index was set to 1.33 + 0.0008 i, and the scatterers were set as hexagonal prism particles (with different size parameters). All the scattering processes were simulated on the same computer. Their computational time for the traditional and improved IIM T-matrix code is presented in Table 1. A histogram of their calculation time is shown in Figure 13.
Remote Sens. 2021, 13, x FOR PEER REVIEW 16 of 20 Figure 13. Comparison of calculation efficiency before and after improving, the abscissa is the size parameter, and the ordinate is the calculation time.

Discussion
The modeling accuracy of the IIM T-matrix method was validated by comparing its results with those of DDA, ADDA, IGOA and the IIM T-matrix method before optimization. In these comparisons, the relative errors of P11 were basically less than 5% where the forward scattering errors of the twelve small prism particles were less than 2%. The absolute errors of P12/P11, P33/P11, P22/P11, P44/P11 and P34/P11 were mostly less than 0.05, and some of them were less than 0.01. Relative errors of extinction cross section were basically less than 1% and the errors of scattering cross section and absorption cross section were less than 3%. Whether regarding the calculation of large particles or small particles, the twelve prisms have better accuracy. The reason is that with the increase in the side edges of prisms, the side surface becomes gradually smooth and closer to the cylindrical shape, and the computational accuracy of rotational symmetric particles is always better than other types of particles in most calculation models.
The computational time for the particles with different sizes was compared by using the traditional and improved IIM T-matrix model. In this comparison, the computational efficiency of the IIM T-matrix model improved notably, especially for particles with a large scale size. When the size parameter reaches 50, its calculation efficiency improves by more than 70%.

Conclusions
The IIM T-matrix model is regarded as one of the most promising models to simulate light scattering by atmospheric non-spherical particles. Compared with other computational models, e.g., DDA and PSTD, this model can not only compute the scattering parameters of randomly oriented particles analytically but also simulate them in much larger sizes effectively. However, because the IIM T-matrix model needs to be iterated layer by layer along the radial direction, the larger the particle size is, the more iterations are required. Therefore, ways to reduce the matrix computation represent an important research topic for the invariant imbedding T-matrix method. In this paper, an efficient realization scheme was proposed to improve the calculation efficiency of the IIM T-matrix method. Additionally, then the IIM T-matrix method's modeling accuracy was Figure 13. Comparison of calculation efficiency before and after improving, the abscissa is the size parameter, and the ordinate is the calculation time.
As can be seen from Table 6, after being improved by the N-fold symmetric geometry, the computational efficiency of the IIM T-matrix model improved notably, especially for particles of a large size. When the size parameter reached 50, its calculation efficiency improved by more than 70%. From Figure 13, it can be seen that as the particle size increases, the improvement in the modeling efficiency becomes much more remarkable. This phenomenon can be explained by the increase in the computational amount of the U-matrix inversion processes. As the particle size becomes large, the dimensions of the U matrix increase dramatically, so the ratio of the computational amount of U-matrix inversion in the total computational task improves as well. From the results, it can be seen that the improved IIM T-matrix model has significant computational advantages for the calculation of large particles.

Discussion
The modeling accuracy of the IIM T-matrix method was validated by comparing its results with those of DDA, ADDA, IGOA and the IIM T-matrix method before optimization. In these comparisons, the relative errors of P 11 were basically less than 5% where the forward scattering errors of the twelve small prism particles were less than 2%. The absolute errors of P 12 /P 11 , P 33 /P 11 , P 22 /P 11 , P 44 /P 11 and P 34 /P 11 were mostly less than 0.05, and some of them were less than 0.01. Relative errors of extinction cross section were basically less than 1% and the errors of scattering cross section and absorption cross section were less than 3%. Whether regarding the calculation of large particles or small particles, the twelve prisms have better accuracy. The reason is that with the increase in the side edges of prisms, the side surface becomes gradually smooth and closer to the cylindrical shape, and the computational accuracy of rotational symmetric particles is always better than other types of particles in most calculation models.
The computational time for the particles with different sizes was compared by using the traditional and improved IIM T-matrix model. In this comparison, the computational efficiency of the IIM T-matrix model improved notably, especially for particles with a large scale size. When the size parameter reaches 50, its calculation efficiency improves by more than 70%.

Conclusions
The IIM T-matrix model is regarded as one of the most promising models to simulate light scattering by atmospheric non-spherical particles. Compared with other computational models, e.g., DDA and PSTD, this model can not only compute the scattering parameters of randomly oriented particles analytically but also simulate them in much larger sizes effectively. However, because the IIM T-matrix model needs to be iterated layer by layer along the radial direction, the larger the particle size is, the more iterations are required. Therefore, ways to reduce the matrix computation represent an important research topic for the invariant imbedding T-matrix method. In this paper, an efficient realization scheme was proposed to improve the calculation efficiency of the IIM T-matrix method. Additionally, then the IIM T-matrix method's modeling accuracy was validated by comparing its results with those of the DDA, ADDA and IGOA methods. The computing time of particles with different sizes was compared by using the traditional and improved IIM T-matrix model. Several conclusions were drawn and are as follows: (1) The improved IIM T-matrix method can accurately simulate the light scattering parameters of the particles with N-fold symmetry. The scattering properties obtained by the IIM T-matrix method demonstrated excellent agreement with those of the well-test scattering models, especially for the polarization characteristics of which absolute differences between our model and DDA were within 0.02. (2) The computational time of the model was notably shortened by nearly 50%. As the particle size increased, the improvement in the modeling efficiency was much more remarkable, reaching approximately 75% for large particles.