## 2 The Principle of the Invariant Imbedding (IIM) T-Matrix Method
### 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\omega t)\), then the equation in frequency domain can be expressed as [40]:  
\[\mathbf{E}(\mathbf{r})=\mathbf{E}_{inc}(\mathbf{r})+\iiint\limits_{V}\mathbf{G }_{0}(\mathbf{r},\mathbf{r}^{\prime})u(\mathbf{r}^{\prime})\mathbf{Z}(\mathbf{r })\mathbf{E}(\mathbf{r}^{\prime})d^{3}\mathbf{r}^{\prime} \tag{6}\]in which, E denotes the electric vector, \(u(\mathbf{r})=k^{2}(\varepsilon_{r}(\mathbf{r})-1)\), \(k\) is the wavenumber, \(\mathbf{r}\) is the position vector,\(\varepsilon_{r}(\mathbf{r})\) is the relative permittivity of the medium, \(\mathbf{E}_{inc}(\mathbf{r})\) is the electric vector of incident wave, \(\mathbf{G}_{0}(\mathbf{r},\mathbf{r}^{\prime})\) denotes the free space dyadic Green's function, and \(\mathbf{Z}(\mathbf{r})\) is a transformation matrix determined by the particle's permittivity, the expressions of \(\mathbf{G}_{0}(\mathbf{r},\mathbf{r}^{\prime})\) and \(\mathbf{Z}(\mathbf{r})\) are shown in Equations (7) and (8). \[\mathbf{G}_{0}(\mathbf{r},\mathbf{r}^{\prime})=(\mathbf{I}+\frac{1}{k^{2}} \
abla\otimes\
abla)\frac{\exp(ik|\mathbf{r}-\mathbf{r}^{\prime}|)}{4\pi| \mathbf{r}-\mathbf{r}^{\prime}|} \tag{7}\]  
\[\mathbf{Z}(\mathbf{r})=\frac{1}{\varepsilon_{r}(\mathbf{r})}\mathbf{\hat{r}} \otimes\mathbf{\hat{r}}+\mathbf{\hat{\theta}}\otimes\mathbf{\hat{\theta}}+ \mathbf{\hat{\varphi}}\otimes\mathbf{\hat{\varphi}} \tag{8}\]  
In the equations above, the operator \"\(\otimes\)\" denotes the dyad product, i.e., if \(\mathbf{a}\) and \(\mathbf{b}\) are two column vectors, then \(\mathbf{a}\otimes\mathbf{b}=\mathbf{a}\cdot\mathbf{b}^{T}\). In order to realize the discretization of Equation (6), both field components and dyadic Green's function \(\mathbf{G}_{0}(\mathbf{r},\mathbf{r}^{\prime})\) 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 \(\mathbf{G}_{0}(\mathbf{r},\mathbf{r}^{\prime})\). The expansion expression of the incident field has been given by Equation (1), which can be written into matrix form, given as:  
\[\mathbf{E}_{inc}(\mathbf{r})=\sum_{n=1}^{\infty}\sum_{m=-n}^{n}\left(Rg \mathbf{M}_{mn},Rg\mathbf{N}_{mn}\right)\left( \begin{array}{c}a_{mn}\\ b_{mn}\end{array} \right)=\sum_{n=1}^{\infty}\sum_{m=-n}^{n}\overline{\mathbf{ \widetilde{Y}}}_{mn}(\theta,\varphi)\overline{\mathbf{\widetilde{J}}}_{n}(kr) \left( \begin{array}{c}a_{mn}\\ b_{mn}\end{array} \right) \tag{9}\]  
where, \(\overline{\mathbf{\widetilde{Y}}}_{mn}(\theta,\varphi)\) is an angular function matrix, and \(\overline{\mathbf{\widetilde{J}}}_{n}(kr)\) is the radial Bessel function matrix, which can be written as:  
\[ \begin{array}{l}\overline{\mathbf{\widetilde{Y}}}_{mn}(\theta,\varphi)= (-1)^{m}\exp(imp\varphi)\bigg{[}\frac{2n+1}{4\pi in(n+1)}\bigg{]}^{1/2}\\ \left(\begin{array}{ccc}0&0&\sqrt{n(n+1)}d_{0m}^{n}(\theta)\\ i\pi_{mn}(\theta)&\tau_{mn}(\theta)&0\end{array}\right)\,\ \overline{\mathbf{ \widetilde{J}}}_{n}(kr)=\left(\begin{array}{ccc}J_{n}(kr)&0\\ 0&\frac{1}{kr}\frac{d}{d(kr)}(krJ_{n}(kr))\\ 0&\frac{\sqrt{n(n+1)}}{kr}J_{n}(kr)\end{array}\right)\end{array}