We are interested in solving direct (as opposed to inverse) electromagnetic (EM) scattering problems with multiple spheroidal scatterers. Denoting the EM filed by \(\mathbf{E}(\mathbf{r},t) \in \mathbb{C}^{3}\), where \(\mathbf{r}\) is a point in \(\mathbb{R}^{3}\) and \(t\) is time, we assume time dependence so that \(e^{-i\omega t}\) factors out and the time dimension can be treated completely independently. We also make certain assumptions about the medium and the scatterer material properties, for which the relevant equation to be satisfied is the homogeneous Helmholtz equation, as well as the appropriate boundary conditions (BCs).
All contributions to the total EM field (“radiation”), including the incident and the scattered fields, are represented by solutions of the vector Helmholtz equation \[\begin{equation}\label{eqn:vectorHelmholtz} \nabla^{2}\mathbf{E}+ k^{2}\mathbf{E}= 0, \end{equation}\] where \(\nabla^{2}\) is the vector Laplacian, \(\mathbf{E}(\mathbf{r}):\mathbb{R}^{3} \mapsto \mathbb{C}^{3}\) must be divergence-free (i.e.~\(\nabla\cdot\mathbf{E}= 0\)), and \(k \in \mathbb{C}\) is the or , related to the \(\mathbf{k} \in \mathbb{C}^3\) \[\begin{equation} k^{2}=\mathbf{k}\cdot\mathbf{k} = \frac{\omega^2}{c^2}\epsilon(\omega) = \frac{4\pi^2}{\lambda^2}\epsilon(\lambda), \end{equation}\] where \(\omega = 2\pi c / \lambda\) is the angular frequency, \(c\) is the speed of light, \(\lambda\) is the wavelength, and \(\epsilon\) is the frequency-/wavelength-dependent dielectric function. It is often convenient to write \(\mathbf{k} = k\mathbf{n}\), where \(\mathbf{n} \in \mathbb{C}^3\) satisfies \(\mathbf{n} \cdot \mathbf{n} = 1\), but is not necessarily a unit vector. In Cartesian coordinates \((x,y,z)\) it is relatively straightforward to verify that a stationary of the form \(\mathbf{E}_{0}\exp{(i\mathbf{\mathbf{k}\cdot\mathbf{r}})}\) is a solution of (\(\ref{eqn:vectorHelmholtz}\)) for any constant \(\mathbf{E}_{0} \in \mathbb{C}^{3}\) satisfying \(\mathbf{E}_{0} \cdot \mathbf{k} = 0\), since EM waves are transverse.
The solution must also satisfy the appropriate BCs at interfaces and/or infinity. Given the geometry of problems we address, the BCs are better expressed in spherical polar coordinates \((r,\theta,\phi)\), with the general solution to (\(\ref{eqn:vectorHelmholtz}\)) being a linear combination of waves. To build up the relevant formulae, it is useful to first consider the Helmholtz equation \[\begin{equation}\label{eqn:scalarHelmholtz} \nabla^{2}E + k^{2}E = 0, \end{equation}\] where \(E(\mathbf{r}):\mathbb{R}^{3} \mapsto \mathbb{C}\) and \(\nabla^{2}\) is the scalar Laplacian.
In spherical polar coordinates, (\(\ref{eqn:scalarHelmholtz}\)) is satisfied by linear combinations of \[\begin{equation} \psi_{nm}^{(\zeta)}(k\mathbf{r}) = z_{n}^{(\zeta)}(kr) Y_{nm}(\theta,\phi), \end{equation}\] where \(z_{n}^{(\zeta)}\) are Bessel functions of kind \(\zeta\in\{1,2,3,4\}\) and order \(n=0\dots\infty\), and \(Y_{nm}\) are the spherical harmonics with \(|m| \leq n\), defined using the Condon-Shortley phase. Hence, a general solution to (\(\ref{eqn:scalarHelmholtz}\)) can be expressed as \[\begin{equation} \label{eqn:ScalGenSol} E^{(\zeta)}(\mathbf{r}; k) = \sum_{p} a_{p} \psi_{p}^{(\zeta)}(k\mathbf{r}) \equiv \mathbf{a}\cdot \pmb{\psi}^{(\zeta)}, \end{equation}\] where \(\pmb{\psi}^{(\zeta)}\) is a vector, \(\mathbf{a}\in \mathbb{C}\) is a vector of coefficients (\(a_{p} \in \mathbb{C}\)), and \[\begin{equation} p(n,m)=n(n+1)+m \end{equation}\] is a composite index introduced for convenience, such that \[\begin{eqnarray} n & = & \mathrm{Int}(\sqrt{p}) \\ m & = & p-n(n+1). \end{eqnarray}\] Note that, while in principle \(p\) can run from \(0\) up to \(\infty\), in practice the series in (\(\ref{eqn:ScalGenSol}\)) is made finite by truncating at some maximum multipole order \(n_{\rm max}\), so that the column vectors \(\mathbf{a} \in \mathbb{C}^{p_{\rm max}}\) and \(\pmb{\psi}^{(\zeta)} \in \mathbb{C}^{p_{\rm max}}\) are of finite length \(p_{\rm max} = n_{\rm max}^{2} + 2n_{\rm max}\).
For our purposes we only require \(\zeta = 1\) (\(z_{n}^{(1)} = j_{n}\), spherical Bessel functions of the first kind) and \(\zeta = 3\) (\(z_{n}^{(3)} = h_{n}\), spherical Hankel functions of the first kind), referred to as the and the functions, respectively, which are linearly independent. Hence, for brevity and notational convenience we may refer to \(\psi_{nm}^{(3)}\) as simply \(\psi_{nm}\), and \(\psi_{nm}^{(1)}\) as \(\widetilde{\psi}_{nm}\). Furthermore, the tilde will also be placed over the coefficients (e.g.~\(\widetilde{\mathbf{a}}\)) to explicitly indicate the regular basis-set flavour. Note that the existence of two linearly independent basis-sets stems from the fact that the underlying (Helmholtz) partial differential equation is of second order and is homogeneous.
Now, solving (\(\ref{eqn:vectorHelmholtz}\)) in spherical polar coordinates leads to two sets of vector-valued functions \[\begin{eqnarray} \mathbf{M}_{nm}^{(\zeta)}(k\mathbf{r}) & = & \frac{1}{\sqrt{n(n+1)}} \nabla \times (\psi_{nm}^{(\zeta)} (k\mathbf{r})\mathbf{r}), \\ \mathbf{N}_{nm}^{(\zeta)}(k\mathbf{r}) & = & \frac{1}{k} \nabla \times \mathbf{M}_{nm}^{(\zeta)}(k\mathbf{r}). \end{eqnarray}\] \(\mathbf{M}_{nm}^{(\zeta)}\) and \(\mathbf{N}_{nm}^{(\zeta)}\) for \(n=0\dots\infty\) and \(|m| \leq n\) form a complete set of mutually orthogonal angular-dependent functions, independent of \(\zeta\). That is, \[\begin{equation} \label{eqn:orthogonality} \int_{\phi=0}^{2\pi} \int_{\theta=0}^{\pi} \mathbf{A}_{nm}^{(\zeta)*}\mathbf{B}_{n'm'}^{(\zeta')} \sin\theta \mathrm{d}\theta \mathrm{d}\phi = \delta_{nn'} \delta_{mm'} \delta_{\mathbf{AB}} f_{n}^{(\mathbf{A},\zeta,\zeta')} (kr), \end{equation}\] where \(\mathbf{A}\) and \(\mathbf{B}\) are placeholders for \(\mathbf{M}\) and \(\mathbf{N}\), \(\mathbf{A}_{nm}^{(\zeta)*}\) denotes the complex conjugate of \(\mathbf{A}_{nm}^{(\zeta)}\), and \(\delta_{ij}\) is the Kronecker delta. The radial function \(f_{n}^{(\mathbf{A}, \zeta,\zeta')}(kr)\) is actually a functional of \(j_{n}(kr)\) and/or \(h_{n}(kr)\), depending on \(\zeta\) and \(\zeta'\), and its exact form also depends on whether \(\mathbf{A}\) is \(\mathbf{M}\) or \(\mathbf{N}\).
A general solution to (\(\ref{eqn:vectorHelmholtz}\)), up to a maximum multipole index \(n_{\rm max}\), can be expressed as \[\begin{eqnarray} \widetilde{\mathbf{E}}({\bf r}; k) & = & \sum\limits_{n=1}^{n_{\rm max}} \sum\limits_{m=-n}^{n} [\widetilde{c}_{1,nm}\widetilde{\bf M}_{nm}(k{\bf r}) + \widetilde{c}_{2,nm}\widetilde{\bf N}_{nm}(k{\bf r})] \nonumber \\ &=& \sum_{q=1}^{2}\sum_{p=1}^{p_{\rm max}} \widetilde{\mathbf{w}}_{qp}(k\mathbf{r}) \widetilde{c}_{qp} \nonumber \\ &=& \sum_{l=1}^{l_{\rm max}} \widetilde{\mathbf{w}}_{l}(k\mathbf{r}) \widetilde{c}_{l} ~ \equiv ~ \widetilde{\mathbf{W}}(k\mathbf{r}) \widetilde{\mathbf{c}}, \label{eqn:VectorGenSol} \end{eqnarray}\] where \(\widetilde{\mathbf{c}}\in \mathbb{C}^{l_{\rm max}}\) is a column vector of coefficients, \(\widetilde{\mathbf{W}}\) is a of dimension \(3\times l_{\rm max}\) (i.e. a row vector composed of column vectors \(\widetilde{\mathbf{w}}_{l(q=1,n,m)} = \widetilde{\mathbf{M}}_{nm}\) and \(\widetilde{\mathbf{w}}_{l(q=2,n,m)} = \widetilde{\mathbf{N}}_{nm}\)), and \(l_{\rm max}\) is the maximal value of another, more convenient, composite index \[\begin{eqnarray} l(q,n,m; n_{\rm max}) & = & (q-1)p_{\rm max} + p(n,m) \nonumber \\ & = & (q-1)n_{\rm max}(n_{\rm max} + 2) + n(n+1) + m \label{eqn:defLindex} \\ & \leq & l_{\rm max} ~=~2n_{\rm max}(n_{\rm max} + 2), \label{eqn:defLmax} \end{eqnarray}\] where \(q \in \{1,2\}\) is sometimes referred to as the or index, which is used to distinguish between the \(\mathbf{M}\) and \(\mathbf{N}\) functions.
The irregular variant of (\(\ref{eqn:VectorGenSol}\)) is obtained by simply removing the overhead tildes (e.g.~\(\widetilde{\mathbf{W}}\to \mathbf{W}\)), which corresponds to switching the radial dependence from \(j_{n}(kr)\) to \(h_{n}(kr)\) throughout. Also, the most general solution to (\(\ref{eqn:vectorHelmholtz}\)) is actually given by the superposition of regular and irregular solutions: \[\begin{equation} \mathbf{E}+ \widetilde{\mathbf{E}}= \mathbf{W}\mathbf{c}+ \widetilde{\mathbf{W}}\widetilde{\mathbf{c}}. \end{equation}\]
Outside a given scatterer, the total field \(\mathbf{E}_{\rm tot}(\bf r) = \widetilde{\mathbf{E}}_{\rm inc}({\bf r}) + \mathbf{E}_{\rm sca}({\bf r})\) is partitioned into a known incident contribution \(\widetilde{\mathbf{E}}_{\rm inc}({\bf r})\) and unknown scattered contribution \(\mathbf{E}_{\rm sca}({\bf r})\). Both partitions are expanded in terms of VSWs up to some multipole order \(n_{\rm max}\), \[\begin{eqnarray} \widetilde{\mathbf{E}}_{\rm inc}({\bf r}) & = & E \widetilde{\mathbf{W}}(\mathbf{r}) \widetilde{\mathbf{c}}_{\rm inc}, \\ \mathbf{E}_{\rm sca}({\bf r}) & = & E \mathbf{W}(\mathbf{r}) \mathbf{c}_{\rm sca}, \end{eqnarray}\] where \(E\) is a parameter determining the field amplitude (usually \(E = |\widetilde{\mathbf{E}}_{\rm inc} |\)), and \(\widetilde{\mathbf{c}}_{\rm inc}\in \mathbb{C}^{l_{\rm max}}\) with \(\mathbf{c}_{\rm sca} \in \mathbb{C}^{l_{\rm max}}\) are the incident and scattered coefficients, respectively. The association of \(\widetilde{\mathbf{E}}_{\rm inc}\) with regular (or ) and \(\mathbf{E}_{\rm sca}\) with irregular (or ) VSWs is a choice motivated by physical reasoning: (i) \(\widetilde{\mathbf{E}}_{\rm inc}(\mathbf{r})\) ought to be well defined everywhere within a finite distance from the origin, which rules out irregular VSWs due to their singular behaviour at \(\mathbf{r}= 0\); and (ii) \(\mathbf{E}_{\rm sca}(\mathbf{r})\) ought to satisfy the outgoing (Sommerfeld?) BCs, requiring that \(|\mathbf{E}_{\rm sca}(\mathbf{r})| \to 0\) at a particular (?) rate as \(|\mathbf{r}| \to \infty\), which rules out regular VSWs due to their divergence in the far field (?). Given the linearity of the governing (Helmholtz) equations, imposition of particular BCs at the scatterer’s surface introduces a linear dependence among \(\widetilde{\mathbf{c}}_{\rm inc}\) and \(\mathbf{c}_{\rm sca}\), causing them to satisfy a relation of the form \[\begin{equation} \label{eqn:genTmatrix} \mathbf{c}_{\rm sca} = \mathbf{T}\, \widetilde{\mathbf{c}}_{\rm inc}, \end{equation}\] where \(\mathbf{T}\) is the so-called “transition” or “transfer” matrix (\(T\)-matrix for short), which depends on the scatterer characteristics and is independent of \(\widetilde{\mathbf{c}}_{\rm inc}\). The \(T\)-matrix elements are (?) determined by matching the BCs over the scatterer’s surface in a particular way, which is generally a non-trivial task. However, for a spherically symmetric, homogeneous scatterer centred at the origin, \(\mathbf{T}\) is a diagonal matrix with the diagonal elements determined analytically by Mie theory. \(T\)-matrices can be calculated numerically for spheroids, as well as more general, axially symmetric scatterers, and… For multiple scatterers, the collective \(T\)-matrix can be built up iteratively from the individual one-body T-matrices.
In multiple scattering problems it is often either necessary or simply more convenient to express the same (partial) field in terms of different basis-sets, mutually related by a translation and/or rotation of the coordinate system. This section is concerned with relations between VSW basis-sets or coefficients defined in different reference frames.
Let \((\theta_{1},\phi_{1})\) and \((\theta_{2},\phi_{2})\) be the spherical angles of the same vector \(\mathbf{r}\) in coordinate systems \(1\) and \(2\), respectively, sharing the same origin. If coordinate system \(2\) is obtained by rotating coordinate system \(1\) through Euler angles \((\alpha,\beta,\gamma)\), here defined in the “\(zyz\)” convention with \(0 \leq \alpha < 2\pi\), \(0 \leq \beta \leq \pi\), and \(0 \leq \gamma < 2\pi\), then \[\begin{eqnarray} \label{eqn:rotTrans} \psi_{nm}(kr,\theta_{2},\phi_{2}) & = & \sum_{\mu=-n}^{n} \psi_{n\mu}(kr,\theta_{1},\phi_{1}) D_{\mu m}^{n}(\alpha,\beta,\gamma), \\ \psi_{nm}(kr,\theta_{1},\phi_{1}) & = & \sum_{\mu=-n}^{n} \psi_{n\mu}(kr,\theta_{2},\phi_{2}) D_{\mu m}^{n}(-\gamma,-\beta,-\alpha), \end{eqnarray}\] where \(D_{\mu m}^{n} = e^{-i\mu\alpha}d_{\mu m}^{n}(\beta) e^{-im\gamma}\) and \(d_{\mu m}^{n}\) are the Wigner \(D\)- and \(d\)-functions. Conveniently, \(\widetilde{\psi}_{nm}\), \(\widetilde{\mathbf{M}}_{nm}\), \(\widetilde{\mathbf{N}}_{nm}\), \(\mathbf{M}_{nm}\) and \(\mathbf{N}_{nm}\) transform in exactly the same manner under rotation, so substituting \(\psi_{nm}\) by a desired basis function in (\(\ref{eqn:rotTrans}\)) will give the appropriate expression. See equations (5.23)-(5.24) of Ref. for details. In our nomenclature, a basis set \(\mathbf{W}^{(2)} \equiv \mathbf{W}(k\mathbf{r}_{2})\) in coordinate system \(2\) is related to \(\mathbf{W}^{(1)} \equiv \mathbf{W}(k\mathbf{r}_{1})\) in coordinate system \(1\) via \[\begin{equation} \label{eqn:basisRot} \mathbf{W}^{(2)} = \mathbf{W}^{(1)} \mathbf{R}(\alpha,\beta,\gamma) \end{equation}\] where \(\mathbf{R}(\alpha,\beta,\gamma)\) is a unitary block-diagonal matrix (of size \(l_{\rm max} \times l_{\rm max}\)), satisfying \[\begin{equation} \mathbf{R}^{-1}(\alpha,\beta,\gamma) = \mathbf{R}^{\dagger}(\alpha,\beta,\gamma) = \mathbf{R}(-\gamma,-\beta,-\alpha) \end{equation}\] and the matrix elements given by \[\begin{equation} R_{ll'}(\alpha,\beta,\gamma) = \delta_{qq'}\delta_{nn'} D_{m'm}^{n}(\alpha,\beta,\gamma), \end{equation}\] where the index \(l(q,n,m)\) is defined in (\(\ref{eqn:defLindex}\)). Note that (\(\ref{eqn:basisRot}\)) also applies to regular waves \(\widetilde{\mathbf{W}}\). Now, if a (regular or irregular) spherical wave expansion is described by a vector of coefficients \(\mathbf{c}^{(1)}\) in coordinate system \(1\), as well as \(\mathbf{c}^{(2)}\) in coordinate system \(2\), then \[\begin{equation} \label{eqn:rotMat} \mathbf{c}^{(2)} = \mathbf{R}^{\dagger}(\alpha,\beta,\gamma) \mathbf{c}^{(1)}, \end{equation}\] which follows from equating the field expansions and using (\(\ref{eqn:basisRot}\)), i.e.
Let us re-label coordinate system \(1\) as \(G\) to indicate a global, space-fixed reference frame, and coordinate system \(2\) as \(L\) for local, scatterer-fixed frame. A \(T\)-matrix \(\mathbf{T}^{(L)}\) expressed in the local frame is transformed into \(\mathbf{T}^{(G)} = \mathbf{R}\mathbf{T}^{(L)} \mathbf{R}^{\dagger}\) in the global frame, where \(\mathbf{R}(\alpha,\beta,\gamma)\) depends on the Euler angles \((\alpha,\beta,\gamma)\) that rotate frame \(G\) onto frame \(L\) (as opposed to \(L\) onto \(G\)). For explanation, consider \[\begin{eqnarray*} \mathbf{c}_{\rm sca}^{(L)} & = & \mathbf{T}^{(L)}\widetilde{\mathbf{c}}_{\rm inc}^{(L)} \\ \mathbf{R}^{\dagger} \mathbf{c}_{\rm sca}^{(G)} & = & \mathbf{T}^{(L)}\mathbf{R}^{\dagger} \widetilde{\mathbf{c}}_{\rm inc}^{(G)} \\ \mathbf{c}_{\rm sca}^{(G)} & = & \underbrace{\mathbf{R}\mathbf{T}^{(L)}\mathbf{R}^{\dagger}}_{\mathbf{T}^{(G)}} \widetilde{\mathbf{c}}_{\rm inc}^{(G)} \implies \mathbf{T}^{(G)} = \mathbf{R}\mathbf{T}^{(L)} \mathbf{R}^{\dagger}. \end{eqnarray*}\]
If the scatterer is rotationally symmetric about the local \(z\)-axis, which is tilted by spherical polar angles \((\theta,\phi)\) relative to the global \(z\)-axis, then \(\alpha = \phi\), \(\beta = \theta\), and the value of \(\gamma\) is irrelevant due to axial symmetry, so we can choose \(\gamma = 0\) to have \(\mathbf{T}^{(G)} = \mathbf{R}(\phi,\theta,0) \mathbf{T}^{(L)} \mathbf{R}(0,-\theta,-\phi)\). See Sec.,5.2 of Ref.~ for more details.
Consider a point \(P\) with coordinates \(\mathbf{r}_1\) in coordinate systems \(1\), with the origin at \(O_{1}\). If we choose another origin \(O_{2}\) with coordinates \(\mathbf{r}_{12}\) relative to \(O_{1}\), then the new coordinates of \(P\) relative to \(O_{2}\) will be \(\mathbf{r}_{2} = \mathbf{r}_{1} - \mathbf{r}_{12}\), as illustrated in Fig.~\(\ref{fig:transadd}\). The translation-addition theorem for vector spherical waves states that, in the limit \(n_{\rm max} \to \infty\), \[\begin{eqnarray} \mathbf{W}(k\mathbf{r}_1) & = & \left\{ \begin{array}{lcl} \mathbf{W}(k\mathbf{r}_2) \widetilde{\mathbf{O}}(k\mathbf{r}_{12}), & {\rm if} & r_{2} > r_{12}, \\ \widetilde{\mathbf{W}}(k\mathbf{r}_2) \mathbf{O}(k\mathbf{r}_{12}), & {\rm if} & r_2 < r_{12}, \end{array} \right. \label{eqn:iVTACS1} \\ \widetilde{\mathbf{W}}(k\mathbf{r}_1) & = & \quad ~\widetilde{\mathbf{W}}(k\mathbf{r}_2) \widetilde{\mathbf{O}}(k\mathbf{r}_{12}), \label{eqn:rVTACS1} \end{eqnarray}\] where \(\widetilde{\mathbf{O}}(k\mathbf{r}_{12})\) and \(\mathbf{O}(k\mathbf{r}_{12})\) are (\(l_{\rm max} \times l_{\rm max}\)) matrices of regular and irregular translation-addition coefficients (TACs), respectively. Note the conditional statement for irregular waves: the transformation depends on the relative length of \(\mathbf{r}_{2}\) and \(\mathbf{r}_{12}\). In Fig.~\(\ref{fig:transadd}\), an irregular basis centred at \(O_{1}\) is mapped onto a regular basis centred at \(O_{2}\) via the irregular TACs. However, if \(O_{2}\) were to the left of the bisector, so that \(r_{2} > r_{12}\), then the irregular basis centred at \(O_{1}\) would be mapped onto an irregular basis centred at \(O_{2}\) via the regular TACs. [What exactly happens when \(r_{12} = r_{2}\), near the bisector, or if one simply uses the wrong expression?] Perhaps a more helpful interpretation is that, if the origin is translated by \(\mathbf{r}_{12}\) (from \(O_{1}\) to \(O_{2}\)), then at points with the new coordinates \(\mathbf{r}_{2}\) satisfying \(r_{2} < r_{12}\) (i.e.~within a ball of radius \(r_{12}\)) the irregular waves are transformed into regular waves, and at all other points the irregular waves remain irregular. Why? To better handle the singularity in the irregular waves, I suppose…
Note that \(O_{1}\) and \(O_{2}\) are themselves points with coordinates \(\mathbf{x}_{1}\) and \(\mathbf{x}_{2}\), respectively, in some designated global “laboratory” frame with a fixed origin \(\mathbf{x}_{0}\). If we denote the lab-frame coordinates of \(P\) by \(\mathbf{r}\), then \(\mathbf{r}_{1} = \mathbf{r}- \mathbf{x}_{1}\), \(\mathbf{r}_{2} = \mathbf{r}- \mathbf{x}_{2}\), and \(\mathbf{r}_{12} = \mathbf{r}_{1} - \mathbf{r}_{2} = \mathbf{x}_{2} - \mathbf{x}_{1} \equiv \mathbf{x}_{21}\). Hence, we follow Stout_et al_ and adopt the shorthand notation \(\widetilde{\mathbf{O}}^{(i,j)} \equiv \widetilde{\mathbf{O}}(k\mathbf{x}_{ij}) = \widetilde{\mathbf{O}}(k\mathbf{r}_{ji})\), and likewise for \(\mathbf{O}^{(i,j)}\), yielding \[\begin{eqnarray} \mathbf{W}^{(i)} & = & \left\{ \begin{array}{lcl} \mathbf{W}^{(j)} \widetilde{\mathbf{O}}^{(j,i)}, & {\rm if} & r_{j} > r_{ij}, \\ \widetilde{\mathbf{W}}^{(j)} \mathbf{O}^{(j,i)}, & {\rm if} & r_j < r_{ij}, \end{array} \right. \label{eqn:iVTACS2} \\ \widetilde{\mathbf{W}}^{(i)} & = & \quad ~\widetilde{\mathbf{W}}^{(j)} \widetilde{\mathbf{O}}^{(j,i)}. \label{eqn:rVTACS2} \end{eqnarray}\] Note the reversal of indices \(i\) and \(j\) in \(\mathbf{x}_{ij} = \mathbf{r}_{ji}\), which has caused some confusion in the past.
To express the translation-addition theorem in terms of the coefficients (\(\mathbf{c}\)) of a VSW expansion, multiply (from the right) both sides of equations (\(\ref{eqn:iVTACS2}\)) and (\(\ref{eqn:rVTACS2}\)) by column vector \(\mathbf{c}^{(i)}\), where the superscript (\(i\)) indicates where the VSW expansion is centred. From inspection of the right-hand side we find that \[\begin{eqnarray} \mathbf{c}^{(j)} & = & \widetilde{\mathbf{O}}^{(j,i)} \mathbf{c}^{(i)}, \quad {\rm if} \quad r_{j} > r_{ij}, \label{eqn:iVTACS3a} \\ \widetilde{\mathbf{c}}^{(j)} & = & \mathbf{O}^{(j,i)} \mathbf{c}^{(i)}, \quad {\rm if} \quad r_{j} < r_{ij}, \label{eqn:iVTACS3b} \\ \widetilde{\mathbf{c}}^{(j)} & = & \widetilde{\mathbf{O}}^{(j,i)}\widetilde{\mathbf{c}}^{(i)}. \label{eqn:rVTACS3} \end{eqnarray}\] Note that (\(\ref{eqn:iVTACS3a}\)), (\(\ref{eqn:iVTACS3b}\)) and (\(\ref{eqn:rVTACS3}\)) are in a similar form to the rotation equation (\(\ref{eqn:rotMat}\)). Also, it is worth reiterating that the above equations hold exactly only in the limit of \(n_{\rm max} \to \infty\).
A general translation, from \(\mathbf{x}_{i}\) to \(\mathbf{x}_{j}\) by \(\mathbf{x}_{ji} = \mathbf{r}_{ij} = (r_{ij},\theta_{ij},\phi_{ij})\), can be separated into three steps:
This factorisation can be expressed in matrix form as \[\begin{equation}\label{eqn:rot-tranz-rot} \mathbf{O}^{(j,i)} = \mathbf{R}(\phi_{ij},\theta_{ij},0) \mathbf{O}^{(z)}(r_{ij}) \mathbf{R}(0,-\theta_{ij}, - \phi_{ij}) \end{equation}\] for the irregular case, where \(\mathbf{O}^{(z)}(r_{ij})\) represents the matrix of axial translation coefficients, many (how many?) of which are zero due to the special case of axial translation along \(z\). Note that the three aforementioned steps correspond to stepwise movement of the local axis, from the perspective of the initial point \(i\), and the transformation corresponds to reading the matrix multiplication in (\(\ref{eqn:rot-tranz-rot}\)) from left to right; but the sequence of steps and the direction of movement is actually reversed from the perspective of the scatterer at the destination point. More importnantly, since all three matrix-factors on the right-hand side of (\(\ref{eqn:rot-tranz-rot}\)) will contain many (how many?) zeroes, operating on a vector of VSW coefficients in a sequence of three steps can actually reduce the scaling of the net computational cost from \(\sim n_{\rm max}^{4}\) to \(\sim n_{\rm max}^{3}\), when the na"ive matrix multiplication on each step is replaced by a custom operation that sums just over the relevant (non-zero) components. See Appendix~\(\ref{app:transadd}\) for more details.
Another potential advantage of using (\(\ref{eqn:rot-tranz-rot}\)) is that, after obtaining the three factors for \(\mathbf{O}^{(j,i)}\), they can be recycled when computing the reverse translation \(\mathbf{O}^{(i,j)}\) to reduce the overall computational cost. First, beware that \(\mathbf{O}^{(z)}(-r_{ij})\) is the inverse of \(\mathbf{O}^{(z)}(r_{ij})\), and note that \(\mathbf{O}^{z}(r_{ij})\) is invariant to interchanging \(i\) and \(j\) (\(r_{ij} = r_{ji} \geq 0\)). Actually, \([\mathbf{O}^{(z)}(r_{ij})]^{-1} = \mathbf{R}(0,\pi,0) \mathbf{O}^{(z)}(r_{ij}) \mathbf{R}(0,-\pi,0)\), where \(\mathbf{R}(0,\pi,0)\) is block-diagonal and each block is -diagonal. Second, since \(\phi_{ij} = \pi + \phi_{ji}\) and \(\theta_{ij} = \pi - \theta_{ji}\), \(\mathbf{R}(\phi_{ji}, \theta_{ji},0)\) can be calculated from \(\mathbf{R}(\phi_{ij}, \theta_{ij},0)\) using the symmetry relation \(d_{\mu m}^{n}(\pi - \theta) = (-1)^{n-m}d_{-\mu m}^{n}(\theta)\) (see eqn.~B.7 in ref.), so that \(D_{\mu m}^{n}(\phi_{ji},\theta_{ji},0) = (-1)^{n-m+1}D_{-\mu m}^{n}(\phi_{ij},\theta_{ij},0)\), where the extra prefactor of \(-1\) comes from multiplying by \(e^{-i\pi} = -1\).
To conclude this section, we note that computing \(\mathbf{O}\) in general is actually more costly than the combined effort of computing \(\mathbf{O}^{(z)}\), \(\mathbf{R}\), \(\mathbf{R}^{-1}\). However, performing matrix multiplication of the three factors in the right sequence (\(\mathbf{R}\mathbf{O}^{(z)}\mathbf{R}^{-1}\)) does not seem to produce the expected result (\(\mathbf{O}\)), perhaps symptomatic of the numerical subtleties involved in the calculation of \(\mathbf{O}\).
This section was written after first reading two papers by Stout_et al_~ and then three earlier papers by Mackowski_et al_~.
For a cluster of \(N\) individual scatterers (each of arbitrary shape), the collectively scattered field \(\mathbf{E}_{\rm sca}^{\rm cl}\) can be separated into additive contributions from the individuals, i.e.: \[\begin{equation} \mathbf{E}_{\rm sca}^{\rm cl} (\mathbf{r}) = \sum\limits_{j=1}^{N} \mathbf{E}_{\rm sca}^{(j)}(\mathbf{r}) = E \sum\limits_{j=1}^{N} \mathbf{W}({\bf r}_{j}) \mathbf{c}_{\rm sca}^{(j)} \label{eqn:scaCoeffsJ} \end{equation}\] where \(\mathbf{r}_{j} = \mathbf{r}- \mathbf{x}_{j}\), \(\mathbf{x}_{j}\) is the position of the \(j\)th scatterer relative to a fixed origin \(\mathbf{x}_{0}\). It is important to stress that the VSW expansion for each \(\mathbf{E}_{\rm sca}^{(j)}(\mathbf{r})\) in (\(\ref{eqn:scaCoeffsJ}\)) is developed in terms of irregular waves centred at \(\mathbf{x}_{j}\), as indicated in the subscript of \(\mathbf{r}_{j}\) and the superscript of the corresponding coefficients \(\mathbf{c}_{\rm sca}^{(j)}\). It is even more important to note that (\(\ref{eqn:scaCoeffsJ}\)) does not actually prescribe how exactly the collective field is partitioned among the individuals.
It is useful to define the excitation field \(\widetilde{\mathbf{E}}_{\rm exc}^{(j)}(\mathbf{r})\) for each scatterer \(j\), and then develop it in terms of regular VSWs centred at \(\mathbf{x}_{j}\), i.e. \[\begin{eqnarray} \widetilde{\mathbf{E}}_{\rm exc}^{(j)} (\mathbf{r}) & \equiv & \widetilde{\mathbf{E}}_\mathrm{inc}(\mathbf{r}) + \sum\limits_{\substack{l=1\\l \neq j}}^{N} \mathbf{E}_{\rm sca}^{(l)} (\mathbf{r}) \nonumber \\ & = & E \widetilde{\mathbf{W}}(\mathbf{r}) \widetilde{\mathbf{c}}_{\rm inc} + E \sum\limits_{\substack{l=1\\l\neq j}}^{N} \mathbf{W}(\mathbf{r}_{l}) \mathbf{c}_{\rm sca}^{(l)} \nonumber \\ & = & E \widetilde{\mathbf{W}}(\mathbf{r}_{j}) \left( \widetilde{\mathbf{c}}_{\rm inc}^{(j)} + \sum\limits_{\substack{l=1\\l\neq j}}^{N} \underbrace{\mathbf{O}^{(j,l)} \mathbf{c}_{\rm sca}^{(l)}}_{\widetilde{\mathbf{c}}^{(j\leftarrow l)}} \right) \quad \mathrm{for~} \|\mathbf{r}_{j}\| < \min\|\mathbf{x}_{j} - \mathbf{x}_{l}\|\label{eqn:excE} \\ & \equiv & E \widetilde{\mathbf{W}}(\mathbf{r}_{j}) \widetilde{\mathbf{c}}_{\rm exc}^{(j)} \label{eqn:excCoeffsJ} \end{eqnarray}\] where \(\widetilde{\mathbf{c}}_{\rm inc}^{(j)} = \widetilde{\mathbf{O}}^{(j,0)} \widetilde{\mathbf{c}}_{\rm inc}\) and \(\widetilde{\mathbf{O}}^{(j,0)} \equiv \widetilde{\mathbf{O}}(\mathbf{x}_{j} - \mathbf{x}_{0})\). In the last equality of (\(\ref{eqn:excE}\)), \(\mathbf{W}(\mathbf{r}_{l})\) is transformed into \(\widetilde{\mathbf{W}}(\mathbf{r}_{j})\), with the correponding coefficients given by \(\widetilde{\mathbf{c}}^{(j\leftarrow l)} = \mathbf{O}^{(j,l)} \mathbf{c}_{\rm sca}^{(l)}\), where the superscript \((j\leftarrow l)\) specifies the scattered field partition of scatterer \(l\) developed (as a regular VSW expansion) about particle \(j\). Note the application of translation-addition theorem clause that applies only inside the ball of radius \(r_{jl}\) (for each \(l\neq j\)) centred at \(\mathbf{x}_{j}\)). This approach is valid only if the translation ball fully contains the target scatterer’s surface, where the BCs are supposed to be matched. While non-overlapping spherical scatters are always guaranteed to satisfy this condition, Fig.~\(\ref{fig:offsetSpheroids}\) shows that spheroids with aspect ratio greater than \(2\) can be problematic, because the t-ball can intersect the target scatterer’s surface if it is sufficiently close (yet still not overlapping). Note that Mackowski explicitly mentions this issue when discussing equations 5-8 in ref.~, but only in the context of spherical scatterers.
From linearity and imposition of appropriate boundary conditions, again, the excitation coefficients \(\widetilde{\mathbf{c}}^{(j)}_{\rm exc}\) defined in (\(\ref{eqn:excCoeffsJ}\)) must be related to the corresponding scattering coefficients \(\mathbf{c}^{(j)}_{\rm sca}\) via \[\begin{equation} \mathbf{c}_{\rm sca}^{(j)} \equiv \mathbf{T}_{1}^{(j)} \widetilde{\mathbf{c}}^{(j)}_{\rm exc}, \label{eqn:Tmat1} \end{equation}\] where the \(\mathbf{T}_{1}^{(j)}\) are the one-body transfer matrices. Note that \(\widetilde{\mathbf{c}}^{(j)}_{\rm exc} = \widetilde{\mathbf{c}}_{\rm inc}\) for an isolated (\(j=N=1\)) scatterer at the origin; and if the scatterer is spherically symmetric then \(\mathbf{T}_{1}^{(1)}\) is a diagonal matrix with Mie coefficients running along the diagonal. For non-spherical scatterers, it may also be helpful to rotate (as well as translate) the coordinates to achieve a favourable orientation. Anyway, from (\(\ref{eqn:excE}\)), (\(\ref{eqn:excCoeffsJ}\)), and (\(\ref{eqn:Tmat1}\)) we obtain an equation expressed just in terms of the field coefficients: \[\begin{equation} \label{eqn:FunMultiScat} \widetilde{\mathbf{c}}^{(j)}_{\rm exc} = \widetilde{\mathbf{c}}^{(j)}_{\rm inc} + \sum\limits_{\substack{i=1\\i\neq j}}^{N} \underbrace{\mathbf{O}^{(j,i)} \mathbf{T}_{1}^{(i)}}_{\mathbf{T}_{1}^{(j \leftarrow i)}} \widetilde{\mathbf{c}}^{(i)}_{\rm exc}, \end{equation}\] which is what Stout_et al_ regard as “the fundamental multiple scattering equation” (eqn.~9 in ref.~). It is helpful to rewrite the linear system (\(\ref{eqn:FunMultiScat}\)) in block-matrix form: \[\begin{equation} \label{eqn:FunMultiScatMat} \begin{bmatrix} \mathbf{I}& -\mathbf{O}^{(1,2)} \mathbf{T}_{1}^{(2)} & \cdots & -\mathbf{O}^{(1,N)}\mathbf{T}_{1}^{(N)} \\ -\mathbf{O}^{(2,1)} \mathbf{T}_{1}^{(1)} & \mathbf{I}& \cdots & -\mathbf{O}^{(2,N)} \mathbf{T}_{1}^{(N)} \\ \vdots & \vdots & \ddots & \vdots \\ -\mathbf{O}^{(N,1)} \mathbf{T}_{1}^{(1)} & -\mathbf{O}^{(N,2)} \mathbf{T}_{1}^{(2)} & \cdots & \mathbf{I} \end{bmatrix} \begin{pmatrix} \widetilde{\mathbf{c}}^{(1)}_{\rm exc} \\ \widetilde{\mathbf{c}}^{(2)}_{\rm exc} \\ \vdots \\ \widetilde{\mathbf{c}}^{(N)}_{\rm exc} \end{pmatrix} = \begin{pmatrix} \widetilde{\mathbf{c}}^{(1)}_{\rm inc} \\ \widetilde{\mathbf{c}}^{(2)}_{\rm inc} \\ \vdots \\ \widetilde{\mathbf{c}}^{(N)}_{\rm inc} \\ \end{pmatrix}. \end{equation}\] We can use (\(\ref{eqn:Tmat1}\)) to substitute the excitation field coefficients for the scattered field coefficients and, assuming the one-body T-matrices are invertible, obtain \[\begin{equation} \label{eqn:FunMultiScat2} \left[\mathbf{T}_{1}^{(j)}\right] ^{-1}\mathbf{c}^{(j)}_{\rm sca} - \sum\limits_{\substack{i=1\\i\neq j}}^{N} \mathbf{O}^{(j,i)} \mathbf{c}^{(i)}_{\rm sca} = \widetilde{\mathbf{c}}^{(j)}_{\rm inc}, \end{equation}\] \[\begin{equation} \label{eqn:FunMultiScatMat2} \begin{bmatrix} \left[ \mathbf{T}_{1}^{(1)} \right]^{-1} & -\mathbf{O}^{(1,2)} & \cdots & -\mathbf{O}^{(1,N)} \\ -\mathbf{O}^{(2,1)} & \left[ \mathbf{T}_{1}^{(2)}\right]^{-1} & \cdots & -\mathbf{O}^{(2,N)} \\ \vdots & \vdots & \ddots & \vdots \\ -\mathbf{O}^{(N,1)} & -\mathbf{O}^{(N,2)} & \cdots & \left[ \mathbf{T}_{1}^{(N)}\right]^{-1} \end{bmatrix} \begin{pmatrix} \mathbf{c}^{(1)}_{\rm sca} \\ \mathbf{c}^{(2)}_{\rm sca} \\ \vdots \\ \mathbf{c}^{(N)}_{\rm sca} \end{pmatrix} = \begin{pmatrix} \widetilde{\mathbf{c}}^{(1)}_{\rm inc} \\ \widetilde{\mathbf{c}}^{(2)}_{\rm inc} \\ \vdots \\ \widetilde{\mathbf{c}}^{(N)}_{\rm inc} \\ \end{pmatrix}. \end{equation}\] We can also rearrange the linear system into \[\begin{equation} \label{eqn:FunMultiScat3} \mathbf{c}^{(j)}_{\rm sca} -\mathbf{T}_{1}^{(j)} \sum\limits_{\substack{i=1\\i\neq j}}^{N} \mathbf{O}^{(j,i)} \mathbf{c}^{(i)}_{\rm sca} = \mathbf{T}_{1}^{(j)} \widetilde{\mathbf{c}}^{(j)}_{\rm inc}, \end{equation}\] \[\begin{equation} \label{eqn:FunMultiScatMat3} \begin{bmatrix} \mathbf{I}& -\mathbf{T}_{1}^{(1)} \mathbf{O}^{(1,2)} & \cdots & -\mathbf{T}_{N}^{(1)}\mathbf{O}^{(1,N)} \\ -\mathbf{T}_{1}^{(2)} \mathbf{O}^{(2,1)} & \mathbf{I}& \cdots & -\mathbf{T}_{1}^{(2)}\mathbf{O}^{(2,N)} \\ \vdots & \vdots & \ddots & \vdots \\ - \mathbf{T}_{1}^{(N)} \mathbf{O}^{(N,1)} & -\mathbf{T}_{1}^{(N)} \mathbf{O}^{(N,2)} & \cdots & \mathbf{I} \end{bmatrix} \begin{pmatrix} \mathbf{c}^{(1)}_{\rm sca} \\ \mathbf{c}^{(2)}_{\rm sca} \\ \vdots \\ \mathbf{c}^{(N)}_{\rm sca} \end{pmatrix} = \begin{pmatrix} \mathbf{T}_{1}^{(1)}\widetilde{\mathbf{c}}^{(1)}_{\rm inc} \\ \mathbf{T}_{1}^{(2)}\widetilde{\mathbf{c}}^{(2)}_{\rm inc} \\ \vdots \\ \mathbf{T}_{1}^{(3)}\widetilde{\mathbf{c}}^{(N)}_{\rm inc} \\ \end{pmatrix}, \end{equation}\] which is equivalent to Mackowski’s eqn.~4 in ref.~. In Mackowski’s earlier papers, however, beware that the same “fundamental” equation (eqn.~13 in ref.~ and eqn.~3 in ref.~) has a plus sign instead of the minus, which may be entirely due to a dubious minus sign featuring in the incident field expansion (see eqn.~4 in ref.~). The dubious minus sign is absent in eqn.~2 of the more recent ref.~, and the subsequent eqn.~4 matches our equation (\(\ref{eqn:FunMultiScat3}\)). Phew!
Note that (\(\ref{eqn:FunMultiScatMat}\)), (\(\ref{eqn:FunMultiScatMat2}\)), and (\(\ref{eqn:FunMultiScatMat3}\)) are all in the form of a general matrix equation \(\mathbf{A}\mathbf{x}= \mathbf{y}\), which can be solved for the column vector \(\mathbf{x}\) without necessarily inverting the matrix \(\mathbf{A}\). However, formal inversion facilitates definition of auxiliary \(T\)-matrix constructions for multiple scatterers.
The system of coupled matrix equations in (\(\ref{eqn:FunMultiScatMat2}\)) can be solved for \(\mathbf{f}^{(j)}\) by inverting the matrix to obtain
\[\begin{equation}\label{eqn:scatCentTMat} \mathbf{c}^{(j)}_{\rm sca} = \sum_{i=1}^{N} \mathbf{T}_{N}^{(j,i)} \widetilde{\mathbf{c}}^{(i)}_{\rm inc} = \underbrace{\left(\sum_{i=1}^{N} \mathbf{T}_{N}^{(j,i)} \widetilde{\mathbf{O}}^{(i,0)}\right)}_{\mathrm{Mackowski's}~ \mathbf{T}_{N}^{(j)}} \widetilde{\mathbf{c}}_{\rm inc} = \underbrace{\left(\sum_{i=1}^{N} \mathbf{T}_{N}^{(j,i)} \widetilde{\mathbf{O}}^{(i,j)}\right)}_{\mathrm{Stout's}~ \mathbf{T}_{N}^{(j)}} \widetilde{\mathbf{c}}^{(j)}_{\rm inc}, \end{equation}\] where \(\mathbf{T}_{N}^{(j,i)}\) represents pairwise \(N\)-body T-matrices, arranged and defined as follows \[\begin{eqnarray} \left[\mathbf{T}_{N}^{(j,i)}\right] & = & \begin{bmatrix} \mathbf{T}_{N}^{(1,1)} & \mathbf{T}_{N}^{(1,2)} & \cdots & \mathbf{T}_{N}^{(1,N)} \\ \mathbf{T}_{N}^{(2,1)} & \mathbf{T}_{N}^{(2,2)} & \cdots & \mathbf{T}_{N}^{(2,N)} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{T}_{N}^{(N,1)} & \mathbf{T}_{N}^{(N,2)} & \cdots & \mathbf{T}_{N}^{(N,N)} \end{bmatrix} = \begin{bmatrix} \left[\mathbf{T}_{1}^{(1)}\right]^{-1} & -\mathbf{O}^{(1,2)} & \cdots & -\mathbf{O}^{(1,N)} \\ -\mathbf{O}^{(2,1)} & \left[\mathbf{T}_{1}^{(2)}\right]^{-1} & \cdots & -\mathbf{O}^{(2,N)} \\ \vdots & \vdots & \ddots & \vdots \\ -\mathbf{O}^{(N,1)} & -\mathbf{O}^{(N,2)} & \cdots & \left[\mathbf{T}_{1}^{(N)}\right]^{-1} \end{bmatrix}^{-1} \label{eqn:TjiMatrixDefinition} \\ & = & \begin{bmatrix} \mathbf{T}_{1}^{(1)} & 0 & \cdots & 0 \\ 0 & \mathbf{T}_{1}^{(2)} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \mathbf{T}_{1}^{(N)} \end{bmatrix} \begin{bmatrix} \mathbf{I}& -\mathbf{O}^{(1,2)} \mathbf{T}_{1}^{(2)} & \cdots & -\mathbf{O}^{(1,N)} \mathbf{T}_{1}^{(N)} \\ -\mathbf{O}^{(2,1)}\mathbf{T}_{1}^{(1)} & \mathbf{I}& \cdots & -\mathbf{O}^{(2,N)} \mathbf{T}_{1}^{(N)} \\ \vdots & \vdots & \ddots & \vdots \\ -\mathbf{O}^{(N,1)}\mathbf{T}_{1}^{(1)} & -\mathbf{O}^{(N,2)} \mathbf{T}_{1}^{(2)} & \cdots & \mathbf{I} \end{bmatrix}^{-1} \\ & = & \begin{bmatrix} \mathbf{I}& - \mathbf{T}_{1}^{(1)} \mathbf{O}^{(1,2)} & \cdots & -\mathbf{T}_{1}^{(1)} \mathbf{O}^{(1,N)} \\ -\mathbf{T}_{1}^{(2)} \mathbf{O}^{(2,1)} & \mathbf{I}& \cdots & -\mathbf{T}_{1}^{(2)} \mathbf{O}^{(2,N)} \\ \vdots & \vdots & \ddots & \vdots \\ -\mathbf{T}_{1}^{(N)} \mathbf{O}^{(N,1)} & -\mathbf{T}_{1}^{(N)} \mathbf{O}^{(N,2)} & \cdots & \mathbf{I} \end{bmatrix}^{-1} \begin{bmatrix} \mathbf{T}_{1}^{(1)} & 0 & \cdots & 0 \\ 0 & \mathbf{T}_{1}^{(2)} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \mathbf{T}_{1}^{(N)} \label{eqn:TjiMatrixFactored2} \end{bmatrix} \end{eqnarray}\]
with the last two lines merely showing how the one-body T-matrices can be factored out in two different ways (left or right). [Q: How to interpret \(\mathbf{O}^{(i,j)}\mathbf{T}_{1}^{(j)}\) and \(\mathbf{T}_{1}^{(i)} \mathbf{O}^{(i,j)}\)?]
The \(\mathbf{T}_{N}^{(j,i)}\) matrices form a complete multiple-scattering solution of the \(N\)-particle system. Stout et al calls these matrices “scatterer-centred tranfer matrices” (see eqn.~12 in ref.~), while Mackowski refers to them as “sphere-centred” (see eqn.~16 in ref.~ and eqn.~4 in ref.~). However, both terms seem equally applicable to \(\mathbf{T}_{N}^{(j)}\), which Mackowski defines one way (see eqn.~61 in ref.~) but does not give it a name, while Stout et al define \(\mathbf{T}_{N}^{(j)}\) differently and call it “individual \(N\)-body transfer matrices” (see eqns.~14 and 27 in ref.~). To avoid any possible confusion, we shall refer to \(\mathbf{T}_{N}^{(j,i)}\) and Mackowski’s \(\mathbf{T}_{N}^{(j)}\) as “pairwise” and “individual” \(N\)-body T-matrices, respectively, and we stress that \(\mathbf{T}_{N}^{(j)}\) differs from \(\mathbf{T}_{1}^{(j)}\) in equation (\(\ref{eqn:Tmat1}\)).
Mackowski’s “contraction” of the index \(i\) in \(\mathbf{T}_{N}^{(j,i)}\) to obtain \(\mathbf{T}_{N}^{(j)}\) essentially links the individually-centred scattering coefficients, \(\mathbf{f}^{(j)}\), with the incident field coefficients, \(\widetilde{\mathbf{a}}\), for the regular VSW expansion centred at the common origin \(\mathbf{x}_{0}\). Mackowski also considers the collective scattering coefficients \(\mathbf{f}^{(0)}\) for the irregular VSW expansion about the common origin, i.e.
\[\begin{equation} \mathbf{E}_{\rm sca}^{\rm cl} (\mathbf{r}) = E \mathbf{W}(\mathbf{r}) \mathbf{c}^{(0)}_{\rm sca} = E \mathbf{W}(\mathbf{r}) \left(\sum\limits_{j=1}^{N} \widetilde{\mathbf{O}}^{(0,j)}\mathbf{c}^{(j)}_{\rm sca}\right),\quad \mathrm{for} \quad \|\mathbf{r}\| > \max\|\mathbf{x}_{j}\|\label{eqn:fcl}, \end{equation}\]
where the second equality relies on a particular clause of the translation-addition theorem, which is valid only outside of the smallest circumscribing sphere (encompassing all \(N\) scatterers) centred at \(\mathbf{x}_{0}\). From (\(\ref{eqn:fcl}\)) and (\(\ref{eqn:scatCentTMat}\)) we have
\[\begin{equation} \label{eqn:Tcl} \mathbf{c}^{(0)}_{\rm sca} = \sum_{j=1}^{N} \widetilde{\mathbf{O}}^{(0,j)} \sum_{i=1}^{N} \mathbf{T}_{N}^{(j,i)} \widetilde{\mathbf{O}}^{(i,0)} \widetilde{\mathbf{c}}_{\rm inc} = \sum_{j=1}^{N} \widetilde{\mathbf{O}}^{(0,j)} \mathbf{T}_{N}^{(j)} \widetilde{\mathbf{c}}_{\rm inc} \equiv \mathbf{T}^{(0)}_{N} \widetilde{\mathbf{c}}_{\rm inc}, \end{equation}\]
which is analogous to (\(\ref{eqn:genTmatrix}\)) and provides an explicit expression for the collective \(\mathbf{T}^{(0)}_{N}\) in terms of the pairwise and (Mackowski’s) individual \(N\)-body T-matrices:
\[\begin{equation} \mathbf{T}^{(0)}_{N} = \sum_{j=1}^{N}\sum_{i=1}^{N} \widetilde{\mathbf{O}}^{(0,j)} \mathbf{T}_{N}^{(j,i)} \widetilde{\mathbf{O}}^{(i,0)} = \sum_{j=1}^{N} \widetilde{\mathbf{O}}^{(0,j)} \mathbf{T}_{N}^{(j)} \end{equation}\]
which is equivalent to Mackowski’s “cluster-centred” T-matrix (see eqn.~19 in ref.~, eqn.~64 in ref.~, and eqn.~29 in ref.~). Again, it is worth reiterating that (\(\ref{eqn:Tcl}\)) is valid only outside of the cluster’s smallest circumscribing sphere centred at the common origin, and Mackowski explicitly acknowledges this fact in refs.~! I suspect this limitation is present in most other studies that utilise a collective T-matrix, such as ref.~
The total energy transfer from an incident plane wave \(\mathbf{E}_{\rm inc}\) into \(\mathbf{E}_{\rm sca}\) and/or \(\mathbf{E}_{\rm int}\) as a result of scattering can be quantified by the cross-section, \(\sigma_{\rm ext}\), which in the current formalism can be calculated from the scatterer-centred coefficients, here for convenience denoted by \(\mathbf{f}^{(j)} = \mathbf{c}_{\rm sca}^{(j)}\) and \(\mathbf{a}^{(j)} = \widetilde{\mathbf{c}}_{\rm inc}^{(j)}\), using
\[\begin{equation} \label{eqn:sigExt} \sigma_{\rm ext} = -\frac{1}{k_{\rm med}^{2}} \Re \left\{ \sum_{j=1}^{N} \mathbf{a}^{(j)\dagger}\mathbf{f}^{(j)}\right\} = -\frac{1}{k_{\rm med}^{2}} \sum_{j=1}^{N}\sum_{q=1}^{2}\sum_{n=1}^{n_{\rm max}} \sum_{m=-n}^{n} \Re \left\{ a^{(j)*}_{qnm} f^{(j)}_{qnm} \right \} \equiv \sum_{j=1}^{N} \sigma_{\rm ext}^{(j)}, \end{equation}\]
where \(\Re\{ \dots \}\) indicates taking the real part of the quantity inside the braces, \(\mathbf{a}^{(j)\dagger}\) is conjugate transpose of the column vector \(\mathbf{a}^{(j)}\), \(a^{(j)*}_{qnm}\) is the conjugate of the vector component \(a^{(j)}_{qnm}\), and \(k_{\rm med}\) is the (real-valued) wave number for the (non-absorbing) medium. Note that (\(\ref{eqn:sigExt}\)) is from eqn.~29 of ref.~, which follows from simplifying eqn.~43 of ref.~ and substituting into eqn.~42 (of same reference); but it apparently differs from eqn.~H.65 of ref.~, where the conjugation is swapped, nonetheless producing the same result since \(\Re\{ab^{*}\} = \Re\{ a^{*}b\}\) for any complex numbers \(a\) and \(b\). Also note that \(\sigma_{\rm ext}\) is separable into additive contributions from individual scatterers, i.e. \(\sigma_{\rm ext} = \sum_{j}\sigma_{\rm ext}^{(j)}\). Mackowski cites the optical theorem when stating analogous expressions for \(\sigma_{\rm ext}\) in eqns.~23-26 of ref.~, missing the minus sign (because of an extra minus sign in Mackowski’s definition of incident field coefficients) and featuring an extra factor of \(4\pi\) (due to different normalisation condition for the angular part of VSWs, which can be seen by comparing the \(E_{mn}\) in eqn.~21 of ref.~ with the \(\gamma_{mn}\) in eqn.~C.22 of ref.~). To spice things up, in Mackowski’s earlier paper the \(4\pi\) is \(2\pi\) in the extinction cross-section formula (eqn.~31 of ref.~), where the minus sign is also missing yet there is no extra minus sign in the definition of incident field coefficients to explain the discrepancy. In Mackowski’s latest paper, the \(4\pi\) again becomes \(2\pi\) in the extinction cross-section formula (eqn.~24 in ref.), which also includes the minus sign, consistent with the absence of the dubious minus sign in the definition of incident field coefficients. Note that Suryadharma and Rockstuhl seem to have the same “\(4\pi\) and no minus sign” convention as in ref.~, and the comparison between eqns.~64-65 of ref.~ may offer an alternative explanation for this pervasive discrepancy.
The scattering cross-section for a given excitation, \(\sigma_{\rm sca}\), can be computed using \[\begin{equation} \label{eqn:sigSca} \sigma_{\rm sca} = \frac{1}{k_{\rm med}^2} \Re\left\{ \sum_{j=1}^{N} \sum_{i=1}^{N} \mathbf{f}^{(j)\dagger} \widetilde{\mathbf{O}}^{(j,i)} \mathbf{f}^{(i)} \right\} = \frac{1}{k_{\rm med}^2} \sum_{j=1}^{N} \sum_{i=1}^{N} \sum_{l=1}^{L} \sum_{l'=1}^{L} \Re \left\{ f_{l}^{(j)*} \widetilde{O}_{ll'}^{(j,i)} f_{l'}^{(i)} \right\} \equiv \sum_{j=1}^{N} \sigma_{\rm sca}^{(j)}, \end{equation}\]
where the composite index \(l\) defined in (\(\ref{eqn:defLmax}\)) has been used for brevity. Note that (\(\ref{eqn:sigSca}\)) is also from eqn.~29 of ref.~, which follows from substituting eqn.~45 of ref.~ into eqn.~42 (of same reference). It seems to differ by a factor of \(4\pi\) from eqn.~35 of ref.~, but for single scatterer (\(N=1\)) it does reduce to eqn.~H.64 of ref.~. My understanding is that (\(\ref{eqn:sigSca}\)) is derived by considering the outward radial component of the Poynting vector for the scatterer field and integrating it over the surface of a sphere enclosing all the scatterers.
The absorption cross-section, \(\sigma_{\rm abs}\) can be deduced using conservation of energy: \[\begin{equation} \sigma_{\rm abs} = \sigma_{\rm ext} - \sigma_{\rm sca}, \end{equation}\] or
Recall that \(\sigma_{\rm ext}\), \(\sigma_{\rm sca}\), and \(\sigma_{\rm abs}\) are defined for a particular incident field direction. The corresponding orientationally averaged quantities are obtained by averaging over all possible orientations, producing \(\langle \sigma_{\rm ext}\rangle\), \(\langle \sigma_{\rm sca}\rangle\) and \(\langle \sigma_{\rm abs}\rangle\).
The \(T\)-matrix formalism is founded on the premise that everything converges below a finite multipole order \(n_{\rm max} < \infty\), i.e.~the field coefficients and \(T\)-matrix elements become negligible for \(n > n_{\rm max}\). Hence, while the formalism (including the translation-addition theorem for VSWs) is exact only in the impractical limit of \(n\to\infty\), sufficiently accurate numerical results can be obtained with finite vectors and matrices truncated above \(n_{\rm max}\). The question is: what is the appropriate choice of \(n_{\rm max}\) for a given system and given numerical scheme? Should one choose just one value of \(n_{\rm max}\), or multiple values, each used at a different stage of the scheme? How exactly does the truncation error propagate from stage to stage? Answering these questions requires some numerical analysis and convergence testing…
To “solve” a multiple scattering problem can mean different things. Usually it means to determine the scattering coefficients \(\mathbf{c}^{(j)}_{\rm sca}\) for a given individual scatterer properties (described by \(\mathbf{T}_{1}^{(j)}\)) and an excitation field (described by \(\widetilde{\mathbf{c}}^{(j)}_{\rm inc}\)) impinging from a particular direction. This kind of “solution” can be achieved by solving the linear system of equations in (\(\ref{eqn:FunMultiScatMat}\)) for \(\mathbf{a}^{(j)}_{\rm sca}\) performing matrix inversion, thus producing a complete description of the scattered field for the given excitation. The linear system would have to be solved again (entirely from scratch) for a different excitation. Performing the matrix inversion to determine the \(N\)-body T-matrices becomes worthwhile only when many impinging directions are to be considered, or when rotationally-averaged quantities are of primary interest.
A brute force approach to solving the multiple scattering problem would be to construct the \(NL \times NL\) matrix in equation (\(\ref{eqn:FunMultiScatMat}\)) and then invert it to obtain the pairwise matrices \(\mathbf{T}^{(i,j)}\). However, this approach is computationally demanding. To help alleviate the cost of one large matrix inversion, Stout_et al_ proposed a recursive scheme where a smaller (\(L \times L\)) matrix is inverted \(N-1\) times.
In the recursive algorithm described by Stout_et al_, the elements of \(\mathbf{T}_{N}^{(j,i)}\) are accumulated recursively from auxiliary subsystems, incrementally built up from one to \(N\) particles. The recursive system is prescribed by the following four equation: \[\begin{eqnarray} \mathbf{T}_{N}^{(N,N)} & = & \left[ \left[\mathbf{T}_{1}^{(N)}\right]^{-1} - \sum_{j=1}^{N-1} \mathbf{O}^{(N,j)} \left( \underline{ \sum_{i=1}^{N-1} \mathbf{T}_{N-1}^{(j,i)} \mathbf{O}^{(i,N)} } \right) \right]^{-1} \equiv \mathbf{S}^{-1}, \label{eqn:TNNN} \\ \mathbf{T}_{N}^{(N,i)} & = & \mathbf{T}_{N}^{(N,N)} \left( \sum_{j=1}^{N-1} \mathbf{O}^{(N,j)} \mathbf{T}_{N-1}^{(j,i)} \right), \quad i\neq N,\\ \mathbf{T}_{N}^{(j,N)} & = & \left( \underline{ \sum_{i=1}^{N-1} \mathbf{T}_{N-1}^{(j,i)} \mathbf{O}^{(i,N)} } \right) \mathbf{T}_{N}^{(N,N)}, \quad j\neq N,\\ \mathbf{T}_{N}^{(j,i)} & = & \mathbf{T}_{N-1}^{(j,i)} + \left( \underline{ \sum_{l=1}^{N-1} \mathbf{T}_{N-1}^{(j,l)} \mathbf{O}^{(l,N)} } \right) \mathbf{T}_{N}^{(N,i)}, \quad j,i\neq N, \end{eqnarray}\] where a common matrix sum has been underlined. Note that only one \(L\times L\) matrix is inverted on each of the \(N-1\) iterations. The inverted matrix becomes ill-conditioned for large \(n_{\rm max}\), but the associated problems can be (at least partly) circumvented by applying the recursive scheme to appropriately “balanced” matrices and coefficients: \[\begin{eqnarray*} \left[ \widehat{\mathbf{T}}_{N}^{(j,i)} \right]_{qp,q'p'} & \equiv & \left[ \mathbf{D}^{(j)} \mathbf{T}_{N}^{(j,i)} [\widetilde{\mathbf{D}}^{(i)}]^{-1} \right]_{qp,q'p'} = \frac{\xi_{n(p)}(k_{\rm M}R_{j})}{\psi_{n'(p')}(k_{\rm M}R_{i})} \left[ \mathbf{T}_{N}^{(j,i)} \right]_{qp,q'p'} \\ \left[ \widehat{\mathbf{T}}_{1}^{(j)} \right]_{qp,q'p'} & \equiv & \left[ \mathbf{D}^{(j)} \mathbf{T}_{1}^{(j)} [\widetilde{\mathbf{D}}^{(j)}]^{-1} \right]_{qp,q'p'} = \frac{\xi_{n(p)}(k_{\rm M}R_{j})}{\psi_{n'(p')}(k_{\rm M}R_{j})} \left[ \mathbf{T}_{1}^{(j)} \right]_{qp,q'p'} \\ \left[ \widehat{\mathbf{O}}^{(j,i)} \right]_{qp,q'p'} & \equiv & \left[ \widetilde{\mathbf{D}}^{(j)} \mathbf{O}^{(j,i)} [\mathbf{D}^{(i)}]^{-1} \right]_{qp,q'p'} = \frac{\psi_{n(p)}(k_{\rm M}R_{j})}{\xi_{n'(p')}(k_{\rm M}R_{i})} \left[ \mathbf{O}^{(j,i)} \right]_{qp,q'p'}, \end{eqnarray*}\] where \(\widetilde{\mathbf{D}}^{(j)}\) and \(\mathbf{D}^{(j)}\) are regular and irregular diagonal matrices with the Riccati-Bessel functions \(\psi_{n}(x) = x j_{n}(x)\) or \(\xi_{n}(x) = x h_{n}(x)\) on the diagonal. (Here, \(j_{n}(x)\) and \(h_{n}(x)\) are, respectively, the spherical Bessel and Hankel functions of the first kind). Note that, what Stout_et al_~ call “matrix balancing” could also be described as “left and right preconditioning”, because a matrix to be inverted is essentially left- and right-multiplied (i.e.~preconditioned) by two different matrices to improve its condition number, which is supposed to make numerical inversion more robust and accurate. Instead of balancing throughout, as Stout_et al_~ seem to do, it may be better to localise the balancing act just at the inversion step in (\(\ref{eqn:TNNN}\)), i.e. \[\begin{eqnarray*} \mathbf{T}_{N}^{(N,N)} = \mathbf{S}^{-1} & = & \left[ \mathbf{D}^{(N)} \right]^{-1} \mathbf{D}^{(N)} \mathbf{S}^{-1} \left[ \widetilde{\mathbf{D}}^{(N)} \right] ^{-1}\widetilde{\mathbf{D}}^{(N)} \\ & = & \left[ \mathbf{D}^{(N)} \right]^{-1} \left[ \widetilde{\mathbf{D}}^{(N)} \mathbf{S} \left[ \mathbf{D}^{(N)} \right]^{-1} \right]^{-1}\widetilde{\mathbf{D}}^{(N)} \equiv \left[ \mathbf{D}^{(N)} \right]^{-1} \widehat{\mathbf{S}}^{-1} \widetilde{\mathbf{D}}^{(N)} \end{eqnarray*}\] where \(\widehat{\mathbf{S}}\) is obtained by balancing \(\mathbf{S}\) analogously to \(\widehat{\mathbf{H}}^{(j,i)}\).
Note that equation (\(\ref{eqn:TNNN}\)) can be factored in two ways: \[\begin{eqnarray} \label{eqn:TNNN2} \mathbf{T}_{N}^{(N,N)} & = & \mathbf{T}_{1}^{(N)} \left[ \mathbf{I}- \sum_{j=1}^{N-1} \mathbf{O}^{(N,j)} \left( \underline{ \sum_{i=1}^{N-1} \mathbf{T}_{N-1}^{(j,i)} \mathbf{O}^{(i,N)} } \right) \mathbf{T}_{1}^{(N)} \right]^{-1} \quad \equiv \quad \mathbf{T}_{1}^{(N)} \mathbf{S}^{-1}_{R}\\ & = & \left[ \mathbf{I}- \mathbf{T}_{1}^{(N)} \sum_{j=1}^{N-1} \mathbf{O}^{(N,j)} \left( \underline{ \sum_{i=1}^{N-1} \mathbf{T}_{N-1}^{(j,i)} \mathbf{O}^{(i,N)} } \right) \right]^{-1} \mathbf{T}_{1}^{(N)} \quad \equiv \quad \mathbf{S}^{-1}_{L} \mathbf{T}_{1}^{(N)} , \end{eqnarray}\] either of which may be preferred if \(\mathbf{T}_{1}^{(N)}\) is difficult to invert. To facilitate the inversion of \(\mathbf{S}_{L}\) and \(\mathbf{S}_{R}\), slightly different balancing (and subsequent unbalancing) should be used: \[\begin{eqnarray} \mathbf{S}_{L}^{-1} = \left[\mathbf{D}_{N}^{(N)}\right]^{-1} \left[ \mathbf{D}_{N}^{(N)} \mathbf{S}_{L} \left[\mathbf{D}_{N}^{(N)}\right]^{-1} \right]^{-1} \mathbf{D}_{N}^{(N)} \equiv \left[\mathbf{D}_{N}^{(N)}\right]^{-1} \widehat{\mathbf{S}}_{L}^{-1} \mathbf{D}_{N}^{(N)}, \\ \mathbf{S}_{R}^{-1} = \widetilde{\mathbf{D}}_{N}^{(N)} \left[ \left[\widetilde{\mathbf{D}}_{N}^{(N)}\right]^{-1} \mathbf{S}_{R} \widetilde{\mathbf{D}}_{N}^{(N)} \right]^{-1} \left[\widetilde{\mathbf{D}}_{N}^{(N)}\right]^{-1} \equiv \widetilde{\mathbf{D}}_{N}^{(N)} \widehat{\mathbf{S}}_{R}^{-1} \left[\widetilde{\mathbf{D}}_{N}^{(N)}\right]^{-1}. \end{eqnarray}\] Here \(\mathbf{S}_{L}\) is balanced like irregular VTACs (or the inverse of a \(T\)-matrix) using only irregular weights, while \(\mathbf{S}_{R}\) is balanced like a \(T\)-matrix using only regular weights. In my experience, \(\widehat{\mathbf{S}}_{R}\) is much better conditioned for inversion than \(\widehat{\mathbf{S}}\), and I haven’t compared with \(\widehat{\mathbf{S}}_{L}\).
From the equations (\(\ref{eqn:TjiMatrixDefinition}\)) and (\(\ref{eqn:TjiMatrixFactored2}\)) it follows that \[\begin{equation} \begin{bmatrix} \mathbf{I}& \cdots & - \mathbf{T}_{1}^{(1)} \mathbf{O}^{(1,N)} \\ \vdots & \ddots & \vdots \\ - \mathbf{T}_{1}^{(N)} \mathbf{O}^{(N,1)} & \cdots & \mathbf{I} \end{bmatrix} \begin{bmatrix} \mathbf{T}_{N}^{(1,1)} & \cdots & \mathbf{T}_{N}^{(1,N)} \\ \vdots & \ddots & \vdots \\ \mathbf{T}_{N}^{(N,1)} & \cdots & \mathbf{T}_{N}^{(N,N)} \end{bmatrix} = \begin{bmatrix} \mathbf{T}_{1}^{(1)} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \mathbf{T}_{1}^{(N)} \end{bmatrix}, \end{equation}\] or, equivalently, \[\begin{equation} \mathbf{T}_{N}^{(i,i)} - \sum_{j\neq i}\mathbf{T}_{1}^{(i)} \mathbf{O}^{(i,j)}\mathbf{T}_{N}^{(j,i)} = \mathbf{T}_{1}^{(i)}, \end{equation}\] which is a linear system of the general form \(\mathbf{A}\mathbf{X} = \mathbf{B}\), where we want to find matrix \(\mathbf{X}\) for a given \(\mathbf{A}\) and \(\mathbf{B}\), which contains multiple right hand sides. That is, each column of \(\mathbf{X}\) and \(\mathbf{B}\) can be treated independently, so we have to solve many linear systems of the form \(\mathbf{A}\mathbf{x}_{\nu} = \mathbf{b}_{\nu}\), where \(\mathbf{x}_{\nu}\) and \(\mathbf{b}_{\nu}\) are the \(\nu\)’th column of \(\mathbf{X}\) and \(\mathbf{B}\), respectively. This approach suggests that matrix inversion should be more difficult than solving a linear system, because the inversion actually requires solving many linear systems…]
Mackowski and Mishchenko (M{&}M) “contract” the second particle index of \(\mathbf{T}_{N}^{(j,i)}\), using \(\mathbf{T}_{N}^{(j)} = \sum_{i} \mathbf{T}_{N}^{(j,i)} \widetilde{\mathbf{O}}^{(i,0)}\), to rewrite the linear system in terms of (as opposed to pairwise) scatterer-centred T-matrices \(\mathbf{T}_{N}^{(j)}\), i.e. \[\begin{equation} \label{eqn:Mackowski} \begin{bmatrix} \mathbf{I}& \cdots & - \mathbf{T}_{1}^{(1)} \mathbf{O}^{(1,N)} \\ \vdots & \ddots & \vdots \\ - \mathbf{T}_{1}^{(N)} \mathbf{O}^{(N,1)} & \cdots & \mathbf{I} \end{bmatrix} \begin{pmatrix} \mathbf{T}_{N}^{(1)} \\ \vdots \\ \mathbf{T}_{N}^{(N)} \end{pmatrix} = \begin{pmatrix} \mathbf{T}_{1}^{(1)} \widetilde{\mathbf{O}}^{(1,0)} \\ \vdots \\ \mathbf{T}_{1}^{(N)} \widetilde{\mathbf{O}}^{(N,0)} \end{pmatrix}, \end{equation}\] or, equivalently, \[\begin{equation} \mathbf{T}_{N}^{(i)} - \sum_{j\neq i}\mathbf{T}_{1}^{(i)} \mathbf{O}^{(i,j)}\mathbf{T}_{N}^{(j)} = \mathbf{T}_{1}^{(i)} \widetilde{\mathbf{O}}^{(i,0)}, \end{equation}\]
where the number of independent linear systems of the form \(\mathbf{A}\mathbf{x}_{\nu} = \mathbf{b}_{\nu}\) is now reduced. [Does this index contraction affect the “optimal” choice of \(n_{\rm max}\)?] M{&}M claim to use (bi)conjugate gradient method to solve \(\mathbf{A}\mathbf{x}_{\nu} = \mathbf{b}_{\nu}\) for each \(\nu\), where the row order (i.e.~length of each column vector) is predetermined using Lorentz/Mie theory result for each (spherical) scatterer in isolation, and the truncation limit for the column order (i.e.~the maximum value of \(\nu\)) is determined on-the-fly from the convergence of each scatterer’s contribution to the collective rotationally-veraged extinction cross-section (see eqn.~66 in ref.~).