1. Introduction

In the electromagnetic scattering problem one considers an incident field ${\textbf E}_\textrm {inc}$, which, after introducing a scattering target, is modified to a field ${\textbf E}_\textrm {tot}$. The scattered field in the region outside the target is the difference in the fields, ${\textbf E}_\textrm {sca}={\textbf E}_\textrm {tot}-{\textbf E}_\textrm {inc}$. The field inside the target is referred to as the internal field ${\textbf E}_\textrm {int}$. A common approach to scattering problems is to expand the (complex) fields in a complete set of wave functions, and to determine the unknown complex expansion coefficients of the scattered and internal fields in terms of the known expansion coefficients of the incident field. Since the boundary conditions on the surface of the target are linear, one obtains linear relations among the expansion coefficients. One is specifically interested in the linear relation that gives the expansion coefficients of the scattered field in terms of those of the incident field. This relation is expressed by the T-matrix.

The traditional approach to computing the T-matrix is Waterman’s null-field method [1]. Other methods have been developed, based on, e.g., the separation-of-variables method [2], the point-matching method [3], and the invariant-imbedding T-matrix method [4,5]. For particles with a basic regular geometry and an impressed regular or irregular small-scale surface structure an approach has been proposed [6,7] that is based on solving an iterative equation (similar to a Lippmann-Schwinger equation) for computing the T-matrix. A short review of the method is given in the next section. A detailed explanation is given in [8].

The convergence properties of this method are, as yet, poorly understood. In previous applications [7] one has simply taken a pragmatic approach and tested the convergence by numerical experiments. However, this can be quite impractical, because it requires us to go through a potentially long iterative process before discovering whether or not the algorithm converges. It would be a significant improvement of the method if we had a convergence prediction before starting the iteration. If the prediction is negative, one can choose a different method, such as a group theoretical method based on the irreducible representations of the particle’s symmetry group [9–11]. It is also mathematically unsatisfactory to develop a numerical method without formulating any criteria for its convergence. The aim of this note is to fill this gap. We will show that a convergence criterion can not only be used for making predictions, but also for systematically improving the chances of convergence of the iterative method.

2. Iterative T-matrix approach

Several methods for numerically solving electromagnetic scattering problems are based on expanding the incident, scattered, and internal fields according to

The linearity of the boundary conditions entails linear relations among the expansion coefficients, which in matrix-vector notation take the form

A potential disadvantage of the last equation is that it relies on inversion of the matrix ${\mathbf {Q}}$, which can become numerically ill-conditioned. Several methods have been contrived to alleviate this problem, such as the use of expansion functions defined in spheroidal coordinates [2], the expansion of the fields about several expansion points [12], as in the method of discrete sources [13], or the use of group-theoretical methods [9,10]. Here we will consider the iterative T-matrix method [6–8], which proceeds as follows. Consider a reference geometry with Q-matrix ${\mathbf {Q}}_0$, as well as a geometry that can be seen as a small perturbation of the reference case with Q-matrix ${\mathbf {Q}}$. We formally introduce the difference of the two Q-matrices,

Equation (10) has one major drawback. It is only an implicit equation for computing the T-matrix, while Eq. (7) is an explicit equation. Equation (10) can be solved by the following iteration scheme.

Can we make a prediction whether or not this iteration scheme will converge? In the following section we derive a sufficient condition for the existence of a unique solution of the method. The main tool for the derivation will be the Banach fixed-point theorem, which we briefly review in the following section.

3. Sufficient condition for the convergence of the iterative T-matrix approach

3.1 Derivation based on the Banach fixed point theorem

For the following derivation we need to borrow some basic concepts from functional analysis.

Contraction mapping. Let $(X,g)$ denote a metric space, where $X$ is a vector space and $g:\, X\times X\rightarrow \mathbb {R}$ denotes a metric defined on $X$. A mapping $\hat {C}:\, X\rightarrow X$ is called a contraction mapping, if there exists a real number $\alpha <1$ such that

Contraction theorem or Banach fixed point theorem. Let $(X,g)$ denote a complete, non-empty set $X$ with metric $g$, and let $\hat {C}:\, X\rightarrow X$ denote a contraction on $X$. Then there exists one and only one point $x\in X$ for which $\hat {C}x=x$. (Such a point is called a fixed point of the mapping $\hat {C}$.)

The proof of this theorem can be found in the standard literature on functional analysis (e.g. [14,15]).

Given an equation of the form $\hat {C}x=x$ (such as the implicit T-matrix equation), we are now in a position to make a prediction about the existence and uniqueness of a solution to this equation. We have to check whether or not $\hat {C}$ is a contraction. This is the main idea for deriving a convergence criterion for the implicit T-matrix equation.

The T-matrix has elements $T_{n,m,q,n',m',q'}$, where $n,n'=1,\ldots ,N_{\mathrm {cut}}$, $m=-n,\ldots ,n$, $m'=-n',\ldots ,n'$, and $q,q'=1,2$. The T-matrix is an element of a vector space ${\cal M}_N$ of $(N\times N)$ matrices, where $N=2N_{\mathrm {cut}}(N_{\mathrm {cut}}+2)$. In this finite dimensional space, we can simply define the matrix norm by

Now we consider the mapping

3.2 Derivation based on a direct comparison test

Consider the matrix equation

4. Applications of the convergence criterion

4.1 Testing convergence by use of the contraction criterion

As an illustration of the method, let us consider a 3D-Chebyshev particle characterised by the surface parameterisation in Eq. (11). We choose $\ell =160$, a size parameter $x=2\pi r_0/\lambda =40$ (where $\lambda$ is the wavelength of light), a refractive index of $m=3+0.1$i (typical for hematite at visible wavelengths), and five different deformation parameters, namely, $\epsilon =$0.0175, 0.02, 0.025, 0.03, and 0.04. This is a well-studied test case that has been considered in [11], where the reciprocity condition has been used [16] to show that the iterative T-matrix approach converges for $\epsilon \le$0.03, but not for $\epsilon =$0.04. In fact, for $\epsilon =$0.0175 it was found that the reciprocity error has dropped below 0.1 % after only three iterations. For $\epsilon =$0.03, the reciprocity error was less than 1.5 % after 39 iteration, while for $\epsilon =$0.04 the reciprocity error increased with the number of iterations, thus indicating divergence of the algorithm. As an illustration, Fig. 1 shows optical properties computed for a deformation parameter of $\epsilon =$0.03. The panels show the linear polarisation $-F_{12}/F_{11}$ (top), as well as $F_{22}/F_{11}$ (bottom), where $F_{ij}$, $i,j=1,\ldots ,4$ denote the elements of the Stokes scattering matrix. The insert in the bottom panel shows the linear depolarisation ratio $\delta _L$ as a function of scattering angle in the backscattering hemisphere, where

 

Fig. 1. Degree of linear polarisation $-F_{12}/F_{11}$ (top), scattering matrix element $F_{22}/F_{11}$ (bottom), and linear depolarisation ratio $\delta _L$ (bottom insert) of randomly oriented 3D-Chebyshev particles with size parameter $x=40$, Chebyshev order $\ell$=160, refractive index $m=3+0.1$i, and deformation parameter $\epsilon =0.03$, computed by performing 6, 15, 27, and 39 iterations of Eq. (13).

Download Full Size | PDF

The left hand side of the inequality in (26) has been implemented into Tsym and applied to the five test cases. Table 1 shows that the matrix norm is smaller than 1 for $\epsilon =$0.0175, 0.2, 0.25 and 0.03. Thus the algorithm is guaranteed to converge. The table also shows the number of iterations $n_{\mathrm {max}}$ that are required to reach convergence in Eq. (13). This can be taken as a measure for the speed of convergence. Evidently, the magnitude of $\left \lVert {\Delta {\mathbf {Q}}\cdot {\mathbf {Q}}_0^{-1}}\right \rVert$ provides us with a clear indication of the speed of convergence.

Table 1. $\left \lVert {\Delta{\mathbf Q}\cdot{\mathbf Q}_0^{-1}}\right \rVert$ from Eq. (26) and number of iterations $n_{\mathrm{max}}$ for Chebyshev particles with size parameter $x=40$, Chebyshev order $\ell$=160, refractive index $m=3+0.1$i, and with different deformation parameters $\epsilon$.

View Table | View all tables in this article

For $\epsilon =$0.04 we have $\left \lVert {\Delta {\mathbf {Q}}\cdot {\mathbf {Q}}_0^{-1}}\right \rVert >1$. Since the contraction criterion is not necessary for convergence, we cannot make any predictions about the convergence of the iterative T matrix equation. The numerical results show that the iteration diverges.

We can conclude that the numerical tests of the iterative T-matrix approach are consistent with the predictions of the contraction criterion.

4.2 Nested iteration to compute the inverse Q-matrix

We can take the iterative solution of the T-matrix problem one step further. In the previous section, we had a reference geometry with Q-matrix $\mathbf {Q}_0$, and a target geometry with Q-matrices $\mathbf {Q}$ and $Rg\mathbf {Q}$. From these three matrices we obtained $\mathbf {T}$ by an iterative approach. Now we want to add an additional intermediate step.

Suppose we have a reference geometry with Q-matrix $\mathbf {Q}_0$, a second reference geometry with Q-matrix $\mathbf {Q}_1$, and a target geometry with Q-matrix $\mathbf {Q}$. To take a specific example, the first reference geometry may be a sphere, the second one a Chebyshev particle of order $\ell$ with a small deformation parameter $\epsilon _1$, and the target geometry a Chebyshev particle of the same order $\ell$ and a somewhat larger deformation parameter $\epsilon > \epsilon _1$. We can determine the T-matrix of the target geometry by the following steps.

We apply this method to the Chebyshev particle considered in the previous section. The first reference geometry is a sphere, the second one a Chebychev particle with deformation parameter $\epsilon _1$, and the target geometry a Chebyshev particle with deformation parameter $\epsilon$. First we take $\epsilon _1$=0.02 and $\epsilon$=0.025. As indicated in Table 2, the convergence criterion yields 0.53 for the first iteration (which is identical with the corrsponding value for $\epsilon$=0.02 in Table 1), and 0.40 for the second iteration. The iterations given in Eqs. (37)–(38) and (39)–(40) converge after $n_{\mathrm {max}}^{(1)}$=5 and $n_{\mathrm {max}}^{(2)}$=4 iterations, respectively. This should be compared to the corresponding values in Table 1 for a single application of the iteration scheme without a second reference geometry, which was $\left \lVert {\Delta {\mathbf {Q}}\cdot {\mathbf {Q}}_0^{-1}}\right \rVert$=0.73 and $n_{\mathrm {max}}$=15. Thus convergence is reached with relative ease in each of the two steps. The price we pay for this is that we need to compute an extra Q-matrix $\mathbf {Q}_1$, and we need to iteratively compute $\mathbf {Q}_1^{-1}$.

Table 2. Convergence criterion in Eq. (35) for Chebyshev particles with size parameter $x=40$, Chebyshev order $\ell$=160, refractive index $m=3+0.1$i, and with different deformation parameters $\epsilon$. The second reference geometry in the iteration has deformation parameter $\epsilon _1$. The number of iterations in the first iteration (37)–(38) is denoted by $n_{\mathrm {max}}^{(1)}$, that in the second iteration (39)–(40) is denoted by $n_{\mathrm {max}}^{(2)}$.

View Table | View all tables in this article

Similarly, for $\epsilon$=0.03 and $\epsilon _1$=0.025 we obtain $\left \lVert {{\mathbf {Q}}_0^{-1}\cdot \Delta {\mathbf {Q}}_1}\right \rVert$=0.73 and $n_{\mathrm {max}}^{(1)}$=15 (which, again, is identical with the corresponding result in Table 1 for $\epsilon$=0.25), and we get $\left \lVert {{\mathbf {Q}}_1^{-1}\cdot \Delta {\mathbf {Q}}}\right \rVert$=0.57, $n_{\mathrm {max}}^{(2)}$=4. Thus, in total, we perform $n_{\mathrm {max}}^{(1)}+n_{\mathrm {max}}^{(2)}$=19 iteration steps, instead of $n_{\mathrm {max}}$=39 (see Table 1).

Finally, for $\epsilon$=0.04 and $\epsilon _1$=0.03, the convergence criterion yields a norm larger than 1, and the iterative method does not converge. To tackle this problem, it would be possible to generalise the method in Eqs. (37)–(40). To this end, one could introduce $m$ reference geometries with Q-matrices $\mathbf {Q}_0,\, \mathbf {Q}_1,\ldots ,\mathbf {Q}_{m-1}$. If we label the target Q-matrix $\mathbf {Q}$ by $\mathbf {Q}_m$, then we can compute $\mathbf {Q}^{-1}$ by the nested iteration scheme

We have made no attempt to test the nested iteration method beyond a second-order outer iteration. Instead, we found a different approach for extending the applicability of the iterative T-matrix method, which is discussed in the following sections.

4.3 Optimising the convergence of the iterative approach by use of the contraction criterion

In the preceding sections, we have applied the contraction criterion to assess whether or not the iterative T-matrix scheme will converge. Now we turn this around and ask: Can we use the contraction criterion to set up the iterative method in a clever way such as to ensure its convergence? To get there, it is important to understand that our choice of the reference matrix $\mathbf {Q}_0$ was in no way unique. In fact, this matrix does not even have to be based on an actual reference geometry; it can be literally any regular matrix. Merely the matrices $\mathbf {Q}$ and $Rg\mathbf {Q}$ have to be computed for the actual target particle. Thus we will now drop the rather constraining assumption that $\mathbf {Q}_0$ is based on a reference geometry. Instead, we will tailor the matrix $\mathbf {Q}_0$ so that the iterative T-matrix scheme (13) converges. The guiding principle in the construction of $\mathbf {Q}_0$ is the contraction criterion (26).

Since we need the inverse of the matrix $\mathbf {Q}_0$, the matrix has to be regular. Further, it greatly simplifies the inversion if we assume that $\mathbf {Q}_0$ be diagonal. In that case the matrix product that occurs in the contraction criterion (26) has components

Let

 

Fig. 2. Generic plot of the functions $f({q_0}_j) = \bar {q}_j/{q_0}_j$ and $g({q_0}_j) = \lvert q_{j,j}/{q_0}_j - 1\rvert$.

Download Full Size | PDF

In summary, we construct the matrix $\mathbf {Q}_0$ by the following recipe.

4.4 Application of the optimisation of $\mathbf {Q}_0$

We revisit the Chebyshev particles encountered in Sect. 4.1. Before, we took $\mathbf {Q}_0$ to be the Q-matrix of the unperturbed sphere. Now we employ the procedure described in Eqs. (57)–(61). We obtain perfect agreement with the phase matrix elements computed earlier (not shown), but with significantly faster convergence rates.

Table 3 shows that the matrix norm in the contraction criterion is substantially reduced in all cases as compared to Table 1. For instance, for a Chebyshev deformation parameter of $\epsilon$=0.03, the norm has been reduced from 0.90 to 0.11. Correspondingly, the number of iteration $n_{\mathrm {max}}$ required for obtaining convergent results decreases from 39 to only 4. For $\epsilon$=0.04 the contraction criterion was previously violated, and we were not able to obtain convergent results. Now the contraction criterion is fulfilled, and convergence is achieved after only 10 iterations.

Table 3. As Table 1, but computed by construction $\mathbf {Q}_0$ according to Eqs. (57)–(61)

View Table | View all tables in this article

Next, we consider a case in which the method described in Eqs. (57)–(61) fails. We consider an oblate spheroid with size parameter $x=20$ and aspect ratio $ar=a/b$=1.5 (where $b$ is the extend along the spheroid’s rotational symmetry axis, and $a$ is the maximum extent in the direction perpendicular to that axis). The refractive index is $m$=1.6 + 0.01i. The surface $r_s(\theta )$ of the spheroid is perturbed with a Chebyshev polynomial analogous to Eq. (11), i.e.

4.5 Construction of a non-diagonal reference matrix $\mathbf {Q}_0$

In Sect. 4.3 we tried something clever to construct $\mathbf {Q}_0$. Now let us try something simple. But as before, we use the contraction criterion as a compass. The goal is to alleviate ill-conditioning problems in the inversion of the Q-matrix. Ill conditioning problems often occur whenever the matrix contains elements of largely different magnitude. Thus, let us introduce a number $P$, and partition the Q-matrix as follows.

We implemented this method with an automated procedure for determining $P$. If $P$ is chosen too large, then $\mathbf {Q}_0$ may end up with row and/or column vectors that only contain zeros, thus making the matrix singular. If $P$ is chosen too small, then $\mathbf {Q}_0$ may be too similar to $\mathbf {Q}$, so that no significant improvement in the conditioning of the matrix is achieved. Thus, our algorithm starts with a high value of $P$ and successively lowers that value until $\mathbf {Q}_0$ becomes regular. Then the norm in the contraction criterion is computed, and $P$ is decreased further until the norm becomes as small as possible, but before $\mathbf {Q}_0$ may become ill-conditioned. In practice, we found that the iterative method converges very fast if the matrix norm is smaller than 0.3. So, we can usually exit the search algorithm whenever $\left \lVert {\Delta {\mathbf {Q}}\cdot {\mathbf {Q}}_0^{-1}}\right \rVert < 0.3$. The inversion of $\mathbf {Q}_0$ is performed by using standard Lapack routines for LU decomposition. Also, we combined the iterative method with group theoretical methods. The use of irreducible representations, as described in [8,10], block-diagonalises the Q-matrix prior to applying the iterative method. This way we can determine a partition number $P$ and a corresponding matrix $\mathbf {Q}_0$ for each block-matrix in $\mathbf {Q}$.

We now return to the Chebyshev spheroid mentioned in Sect. 4.4. The method of constructing a diagonal matrix $\mathbf {Q}_0$, developed in Sect. 4.3, failed to give convergent results for this problem. The simple construction of a non-diagonal matrix $\mathbf {Q}_0$ proposed here converges after only 3 iterations. The matrix norm in the contraction criterion is smaller than 0.3.

Figure 3 shows the elements $F_{11}$ (top), $-F_{12}/F_{11}$ (centre), and $F_{22}/F_{11}$ (bottom) of the scattering matrix for randomly oriented Chebyshev spheroids (solid line) and for randomly oriented unperturbed spheroids (dashed line). While the phase function (top) is hardly impacted by the surface perturbation, the other two elements display noticeable differences between smooth and rough spheroids. The differences are less dramatic than in Fig. 1. The main reason for this is that the spheroids considered here are optically softer and less strongly absorbing, which reduces the significance of surface roughness — see the discussion in [18].

 

Fig. 3. Phase function $F_{11}$ (top), degree of linear polarisation $-F_{12}/F_{11}$ (centre), as well as $F_{22}/F_{11}$ (bottom) for smooth oblate spheroid (dashed line) and Chebyshev spheroid (solid line). The particles have a size parameter of $x$=20, an aspect ratio of 1.5, and a refractive index of $m$=1.6+0.001i. The Chebyshev order and deformation parameter are $\ell$=160 and $\epsilon$=0.03, respectively.

Download Full Size | PDF

5. Summary

The main goal of this study was to derive a convergence criterion for the iterative T-matrix method. From the contraction theorem, we obtained a sufficient condition for the convergence of the method, which is given in Eq. (13). Specifically, if there exists an $\alpha <1$ so that

The applications we showed were meant to illustrate the potential usefulness of the contraction criterion. We employed Eq. (65) to assess the chances of convergence of the iterative T-matrix method. First, we iteratively computed the T-matrix by use of a single reference geometry. Second, we performed a two-step iterative computation of the inverse Q-matrix by use of two reference geometries. In either case, the numerical tests confirm that the condition correctly predicts convergence of the algorithm. Also, if the criterion is fulfilled, then the magnitude of the matrix norm $\left \lVert {\Delta {\mathbf {Q}}\cdot {\mathbf {Q}}_0^{-1}}\right \rVert$ provides us with an indication of the prospective speed of convergence. Such a criterion is of great practical value, since it allows us to decide beforehand whether or not the iterative approach is a promising method for computing the T-matrix.

Next, we employed the contraction criterion in optimising the convergence of the method. To this end, we used the contraction criterion as a compass in constructing the reference matrix $\mathrm {Q}_0$. In a first attempt, we assumed $\mathrm {Q}_0$ to be diagonal and determined its elements by minimising the matrix norm in the contraction criterion. In a second attempt, we dropped the assumption that $\mathrm {Q}_0$ be diagonal. The obvious advantage is that one has more degrees of freedom in the construction of $\mathrm {Q}_0$. A major disadvantage is that it is not straightforward for a general matrix to impose the constraint that $\mathrm {Q}_0$ be regular, and even well-conditioned. We approached this problem by simply partitioning the Q-matrix into a matrix $\mathrm {Q}_0$ containing all the large elements of $\mathrm {Q}$, and a matrix $\Delta \mathrm {Q}$ containing all the small elements of $\mathrm {Q}$. The demarcation between small and large elements was defined by a parameter $P$. Although this is only an ad hoc approach, it gave us the flexibility of adjusting $P$ by use of the contraction criterion, with the constraint that the matrix $\mathrm {Q}_0$ must be well-conditioned. Numerical examples confirmed that the use of the contraction criterion in the construction of $\mathrm {Q}_0$ can, indeed, improve the range of applicability of the iterative T-matrix method, as well as speed up its convergence. The partitioning method, in conjunction with group theoretical methods, seemed to be a particularly promising candidate for further studies.

The matrix $\mathrm {Q}_0$ determines the starting value of the iteration — see Eq. (12) — which is critical for iterative methods. Other methods for constructing $\mathrm {Q}_0$ than the ones considered here are conceivable. For instance, it may be worth to investigate whether the optimisation procedure discussed in sect. 4.3 can be mapped to a corresponding eigenvalue analysis of the matrix $\Delta \mathrm {Q} \cdot \mathrm {Q}_0^{-1}$. The salient point in the context of this work is that the contraction criterion provides us with a powerful tool for optimising the choice of the reference matrix $\mathrm {Q}_0$ in the iterative T-matrix method.