1. Introduction

Scattering of electromagnetic (EM) radiation is a ubiquitous phenomenon, which determines how we see the world (literally), plays an important role in regulating the temperature of our atmosphere, and underpins various spectroscopic techniques used for materials characterisation. Classical theoretical description of EM scattering is provided by the vector Helmholtz equation (which follows from Maxwell’s equations under certain physical assumptions), and when supplemented with appropriate numerical methods and approximations it constitutes a framework for computational modelling of optical properties of matter. Within this framework the $T$-matrix formalism introduced by Waterman [1–4] is particularly effective. Unlike purely numerical methods such as finite-element, the $T$-matrix framework shares strong analytical connections with the exact Mie [5] solution (expressed as a series of electric and magnetic multipoles) for spherical scatterers, and extends it to nonspherical scatterers such as spheroids [6]. This approach provides efficient and accurate simulations of orientation-averaged optical properties, as well as near-field calculations, which are of great value in applications such as nano-optics [7,8], and can serve as benchmark results to test other numerical methods. The superposition $T$-matrix method (STM) [9–11] generalises the framework to collections of disjoint scatterers, such as a clump of interstellar dust [12,13], a collection of ice crystals [14] or soot particles [15,16] in air, or colloidal clusters of plasmonic nanoparticles in solution [8,17–19]. The major stated restriction of STM method is the condition that the particles’ smallest circumscribed spheres should not overlap [9,11,20,21], which has limited its application to either sparse clusters or spherical particles. This “no-overlap” condition can be traced back to Peterson and Ström [9], who used it to avert potential problems when translating the underlying series expansions. Additional problems may also be expected from violation of the Rayleigh hypothesis (RH) [22–24], which postulates that the scattered field immediately outside a particle can be expanded into a series of outgoing vector spherical wave functions (VSWFs). This hypothesis is known to be invalid in certain regions [25], and the scattered field derived from the extended boundary condition method (EBCM) [1,4,26] may therefore diverge in parts of a particle’s circumscribed sphere. This divergence would seem to provide another justification for the no-overlap condition, since two adjacent particles affect one another through their scattered field. Most derivations [4,27] and applications [28,29] of the STM method either explicitly start with the no-overlap assumption or abide by it implicitly by considering only sparse configurations [30,31]. While the assumption is irrelevant if all the constituent scatterers are spheres [11,32,33], it becomes increasingly more restrictive for highly anisotropic scatterers, effectively ruling out potential studies of densely packed aggregates. Recently, modifications have been proposed [21,34] to overcome this limitation, but at the expense of added complexity and computational demands.

We here demonstrate on numerical examples that the “no-overlap” condition is in fact not strictly necessary. For this, we systematically examine the effect of overlap on optical properties of two adjacent prolate spheroids and show that the conventional superposition $T$-matrix method [11,28,32,33,35] still produces valid results even for large overlap, provided the calculations are performed with sufficient numerical precision. We implemented the STM method in a custom-made program called terms ($T$-matrix for Electromagnetic Radiation with Multiple Scatterers), and validated the results against finite-element (FE) calculations [36]. We focus on an experimentally-relevant example in plasmonics, interacting silver nanospheroids, but the interpretations and conclusions are general. The apparent but unexpected validity of the STM method is discussed in the context of the Rayleigh hypothesis as well as the translation-addition theorem for vector spherical wave functions [27,37,38]. This work aims to dispel some lingering misconceptions in the literature about the no-overlap assumption and to demonstrate the surprising robustness of the STM method with numerically stable implementations.

2. Theoretical background

We start with a brief summary of the underlying theoretical formalism.

2.1 The $T$-matrix method for a single scatterer

The $T$-matrix formalism introduced by Waterman [1–3] can be viewed as a generalisation of Mie theory to nonspherical scatterers.

At a point coordinate $\mathbf {r}$ exterior to the scatterer (zones 1 and 2 in Fig. 1(a)), the local EM field $\mathbf {f}(\mathbf {r})$ is a superposition of the incident field $\mathbf {f}_\textrm {inc}(\mathbf {r})$ and the scattered field $\mathbf {f}_\textrm {sca}(\mathbf {r}) = \mathbf {f}(\mathbf {r}) - \mathbf {f}_\textrm {inc}(\mathbf {r})$, each represented by a series expansion in terms of VSWFs. From physical considerations, $\mathbf {f}_\textrm {inc}(\mathbf {r})$ and $\mathbf {f}_\textrm {sca}(\mathbf {r})$ are described by the regular and irregular VSWFs, herein denoted by $\widetilde {\mathbf {w}}_{snm}(\mathbf {r})$ and $\mathbf {w}_{snm}(\mathbf {r})$, respectively, with $s$ specifying the mode (electric or magnetic), $n =1,2,3,\dots$ the multipole order, and $m = -n,\dots ,0,\dots, n$ the degree. After mapping $s$, $n$, and $m$ onto a composite index $l(s, n, m)$, the field can be expressed compactly as

 

Fig. 1. (a) Prolate spheroid with major semi-axis $c$, minor semi-axis $a$, and focal distance $f = \sqrt {c^{2}-a^{2}}$. Dashed black and red circumferences represent the circumscribed (CS) and focal sphere (FS) of radius $c$ and $f$, respectively. Digits $1$, $2$, and $3$ demarcate the far, near, and interior zones. (b) Two spheroids and three reference frames with aligned orientation but different origins: $O$, $O_{i}$, and $O_{j}$. (c–f) Clusters of two identical spheroids side by side, with $c$ varied to increase the degree of overlap: (c) CSs intersect but not FSs, (d) FSs intersect, (e) FSs intersect with other spheroid, (f) $r_{ji}<c$. In the latter cases, the solid red circumference represents the transformed singularity of irregular VSWFs when translated by a vector $\mathbf {r}_{ji}$.

Download Full Size | PDF

2.2 The Rayleigh hypothesis

That Eq. (1) provides a valid (i.e. converging to the correct value as the multipole order $n$ increases) representation of $\mathbf {f}(\mathbf {r})$ everywhere outside the scatterer is only a conjecture, often referred to as the RH. For $\mathbf {f}_\textrm {sca}(\mathbf {r})$ arising from an arbitrarily shaped scatterer, the validity of $\mathbf {W}(\mathbf {r})\mathbf {a}$ has never been proven inside the smallest circumscribed sphere. Whether the RH holds true anywhere inside that sphere but still exterior to the scatterer (zone 2 in Fig. 1(a)) is not generally well known, though for spheroids the RH has been shown to break down only inside the focal sphere (red region in Fig. 1(a)) [25]. This breakdown puts limits on where the local field can be evaluated with a series expansion, but it does not affect the far-field quantities in (3), because (2) does not actually require (1) to be valid everywhere in the exterior. The $T$-matrix ansatz only requires that $\widetilde {\mathbf {a}}$ and $\mathbf {T}$ be well-defined, in which case $\mathbf {a}$ can be uniquely determined using (2), thus allowing the calculation of cross-sections using (3), and $\mathbf {f}(\mathbf {r})$ can be evaluated from $\mathbf {a}$ and $\widetilde {\mathbf {a}}$ using (1) anywhere inside the region of validity.

2.3 Superposition $T$-matrix method

For $N$ particles with known individual $T$-matrices $\mathbf {T}_{i}$, where $i=1,\dots ,N$, the collectively scattered field is a linear superposition of $N$ contributions,

2.4 Beyond the no-overlap condition

Let us consider what could happen in the special case of two prolate spheroids placed side-by-side (Fig. 1). For a single spheroid, the region of divergence is not the full circumscribed sphere (CS), but only the focal sphere (FS). Our intuition suggests that there should be no issues for the dimer in Fig. 1(c), where the FSs do not intersect, even if the CSs do. This would already relax the no-overlap condition. We can however expect something interesting and potentially problematic when the FSs overlap (Fig. 1(d)), or worse, extend inside the other spheroid (Fig. 1(e)). An even more extreme situation occurs when the centre-to-centre distance $r_{ji} < c$ (Fig. 1(f)), where $c$ is the major semi-axis, with potentially problematic issues also arising from the translation of the irregular VSWFs.

3. Results

We now turn to the numerical calculation of optical properties using terms. Our implementation supports double and quad precision calculations, using Amos’ library to calculate the Bessel functions [42], and BLAS/LAPACK routines to solve (5). It can also import accurate $T$-matrices calculated for a spheroid using our improved algorithms, smarties [43]. Note that our discussion of the underlying formalism was based on infinite series expansions, but in practice all the series are truncated at a specified multipole order $n_\textrm {c}$ — an input parameter for our calculations. After the linear system in (5) is set up for a given $n_\textrm {c}$, preconditioned using the analytic matrix-balancing scheme of Stout et al. [28], and then numerically solved, the resulting coefficients $\mathbf {c}_{1}^{(1)},\dots ,\mathbf {c}_{N}^{(N)}$ are “contracted” to the common origin using the regular offset matrices $\widetilde {\mathbf {O}}^{(0, j)}$ to produce $\mathbf {a} = \sum _{j=1}^{N} \widetilde {\mathbf {O}}^{(0, j)}\mathbf {c}_{j}^{(j)}$, which is plugged into (3) to calculate optical cross-sections. Near-fields are calculated from Eqs. (1,4) and presented in terms of the electric field intensity enhancement factor $I/I_0$, where $I=|\mathbf {f}|^{2}$ is the total field intensity, and $I_0=|\mathbf {f}_\textrm {inc}|^{2}$ is the incident field’s intensity.

In order to cover the overlap scenarios illustrated in Figs. 1(c)–1(f), we consider a sequence of silver prolate spheroids dimers with $a=20$ nm, $r_{ji} = 50$ nm ($10$ nm gap), and $c=20$ (2 spheres), $28$ (case Fig. 1(c)), $36$ (Fig. 1(d)), $44$ (Fig. 1(e)), $52$ (Fig. 1(f)), and $60$ nm (more extreme than Fig. 1(f)). We use Yang et al.’s [44] dielectric function for silver, choose water ($\sqrt {\epsilon _\textrm {med}} = 1.33$) as the surrounding medium, and focus on an incident plane wave $\mathbf {f}_\textrm {inc}(\mathbf {r}) = \mathbf {E}_{0}\exp (\mathbf {k}\cdot \mathbf {r})$ satisfying $\mathbf {k}\cdot \mathbf {r}_{ji} = 0$ and $\mathbf {E}_{0}\cdot \mathbf {r}_{ij} = E_{0}r_{ij}$, i.e. travel direction along the $c$-axis to maintain a constant geometric cross-section $\sigma _\textrm {geo} = 2\pi a^{2}$, and linear polarisation along the dimer axis $\mathbf {r}_{ji}$ to strongly excite modes located around the gap. Spectra of the scattering ($\sigma _\textrm {sca}$) and absorption ($\sigma _\textrm {abs} = \sigma _\textrm {ext} - \sigma _\textrm {sca}$) cross-sections are shown in Fig. 2, where the $T$-matrix results are compared with FE calculations (using COMSOL following Ref. [36]), and with a reference $T$-matrix calculation for two equivalent spheroids in isolation, to show the effect of inter-particle coupling. This particular configuration was chosen for illustration as it is relevant to plasmonics studies, but our conclusions are more general. The increased interaction between spheroids evidently produces a complicated change in the spectrum, while for non-interacting spheroids the dipolar resonance around 400 nm is simply amplified and blue-shifted slightly as $c$ increases. Leaving the specifics of the interaction for further studies, here we only highlight the good agreement between the $T$-matrix and FE results for $c \leq 44$ nm, which is unexpected considering the high degree of overlap ($r_{ij}/c \approx 1.14$), already well beyond the no-overlap condition ($r_{ij}/c > 2$). The agreement deteriorates for the more extreme case of $c=52$ nm, where the $T$-matrix spectra remain smooth but do not quite fit the FE data; and the agreement strongly deteriorates for $c=60$, with wild fluctuations in the spectrum indicating numerical problems.

 

Fig. 2. Spectra of the scattering and absorption cross-sections (normalised by the geometrical cross-section $\sigma _\textrm {geo} = 2\pi a^{2}$) for the dimers in Figs. 1(c)–1(f). Six $c$-values are considered with fixed $a = 20$ nm and gap of $10$ nm. Finite-element calculations (black circles) are compared with $T$-matrix data for $n_\textrm {c} = 20$ (solid red line). The reference $T$-matrix data (dashed blue line) is for two spheroids isolated from each other (non-interacting).

Download Full Size | PDF

All the $T$-matrix calculations in Fig. 2 were done in double precision with $n_\textrm {c} = 20$. However, as the spheroids become more elongated, higher multipolar orders are required to accurately describe the scattered field; the resulting VSWFs series, in turn, require increasing numerical precision. This is clearly illustrated in Fig. 3, where convergence of the extinction spectrum with respect to increasing $n_\textrm {c}$ for $c=52$ nm displays wild fluctuations for $n_\textrm {c} \gtrsim 30$. To check that these fluctuations occur because of limited numerical precision, rather than more “intrinsic” problems, we calculated $\sigma _\textrm {ext}$ in double- and quadruple-precision for increasing $n_\textrm {c} \leq 50$ at two resonant wavelengths, $388$ and $484$ nm. The results (Fig. 3) demonstrate numerical instability of the double-precision values from about $n_\textrm {c} \gtrsim 25$, while the quadruple-precision values continue converging to the FE result (until also showing signs of numerical instability for $n_\textrm {c} \gtrsim 50$). Hence, the underlying problems are clearly alleviated by resorting to higher precision, and they could potentially be circumvented entirely with better numerical algorithms. The difficulties in convergence are not therefore a fundamental limitation, and the no-overlap condition appears unwarranted in these cases. In fact, converging results are obtained even for $r_{ji} < c$, revealing a counter-intuitive robustness of the STM method in situations where the aforementioned singularity spheres cross the target scatterer boundary. One would have reasonably expected that $r_{ij}>c$ was a necessary condition to ensure the validity of the translation relations on the surface of the other scatterer.

 

Fig. 3. Top: Spectra of normalised extinction cross-section for the $c=52$ nm case in Fig. 2, calculated using terms with different $n_\textrm {c}$ values and compared with finite-element (FE) calculations (black dots). Bottom: The FE value (thick line) at specified wavelength compared to terms calculations in double (dashed line) and quad (solid line) precision for increasing $n_\textrm {c}$.

Download Full Size | PDF

Another counter-intuitive result comes from consideration of the near field calculated using (4) for the dimer with $c=36$ at the resonant wavelength $\lambda =436$ nm. The local field intensity and the error relative to FE results are mapped in Figs. 4(a)–4(b), showing slices along the three planes of symmetry (note that we do not evaluate the fields inside spheroidal particles, which presents no additional difficulties but is not the focus of this study). The exterior field intensity converges outside the individual focal spheres, as for a single spheroid, but surprisingly it also converges (and agrees with FE calculations) inside the focal spheres near the centre of the overlap region. The interaction between the spheroids appears to somehow increase the region of validity of the RH in parts of the overlap region, where (4) converges to the correct value. Convergence with increasing $n_\textrm {c}$ is demonstrated in Fig. 4(c), where the intensity profile is scanned in the direction of the incident wave vector $\mathbf {k}$ (parallel to the spheroids’ $c$-axis) passing through the origin (the point where all three symmetry planes intersect). Convergence and good agreement with the FE results are reached for $n_\textrm {c} < 30$ without any numerical instabilities in double precision. However, much slower and more oscillatory convergence happens for the dimer with $c=44$ nm, as shown in Fig. 4(d), where the profile from $T$-matrix calculations oscillates with increasing frequency and diminishing amplitude as $n_\textrm {c}$ rises. Although these calculations do not fully converge for $n_{c} \leq 60$ in quadruple-precision and become numerically unstable for $n_\textrm {c} \geq 70$, the trend nonetheless clearly shows that (1) is convergent in the overlap region, and getting closer to the limit of convergence will require going beyond quadruple-precision with the current approach.

 

Fig. 4. (a) Local-field enhancement factor $I/I_{0}$ and (b) magnitude of the error relative to FE results mapped along three symmetry planes for the $c=36$ nm case in Fig. 2 at $\lambda = 436$ nm, calculated using terms in double precision with $n_\textrm {c} = 40$. White regions outside the spheroids but inside their focal spheres indicate divergence. (c) Line scan of the enhancement factor along $\mathbf {k}$ and through the point of intersection of all three symmetry planes. Calculations with different values of $n_\textrm {c}$ show convergence to the FE results, all plotted on linear and logarithmic scales for clarity. (d) Analogous line scan for the dimer with $c=44$ nm, calculated using quad precision and showing slower convergence in the overlap region.

Download Full Size | PDF

As a final and even more provocative example we consider a silver sphere of radius $r=15$ nm fully inside the focal sphere of a silver spheroid with $a=20$ and $c=60$ nm. This means that the scattered field expansion for the spheroid by itself diverges on the entire region where we add the sphere. A schematic of this dimer configuration is shown in Fig. 5, together with the extinction spectrum and a map of the local field intensity for the same incidence as in previous examples. Again, without delving into a modal interpretation of the multiple resonances, we only stress the fact that the conventional STM method produces valid spectra for far-field optical properties, albeit at a high computational cost of quadruple precision and high multipole orders. And again, the presence of the second scatterer enlarges the region of convergence of the near-field solution. Remarkably, we find that the field inside the sphere can be calculated, even though it is fully contained within the region where the Rayleigh Hypothesis is invalid for the neighbouring spheroid. The internal field inside the sphere is computed using Mie theory from the sphere-centred excitation field coefficients (a similar procedure can be used to compute the internal field of the spheroids, but the practical implementation is more cumbersome and outside the scope of this study).

 

Fig. 5. Top: normalised extinction spectra for a silver sphere of radius $r=15$ nm near a silver spheroid of $a=20$ and $c=60$ nm, with $\sigma _\textrm {geo} = \pi (a^{2}+r^{2})$. The geometry has two symmetry planes and $5$ nm gap ($r_{ij}=40$ nm). Bottom: maps of the local-field intensity enhancement (calculated using terms in quadruple precision with $n_\textrm {c} = 50$) and the error relative to FE calculations for $\lambda =430$ nm.

Download Full Size | PDF

4. Discussion and conclusion

When considering multiple scattering by nonspherical particles in close proximity, the STM method is often assumed to be restricted to configurations where the smallest circumscribed spheres of each scatterer do not intersect. Our results demonstrate that such no-overlap condition is not necessary and can be relaxed. Does it mean that the STM method has no such restriction and that scatterers may be placed arbitrarily close to each other, or do we need to consider another, less-restrictive condition? Beside the Rayleigh Hypothesis, the other natural condition that arises in the derivation of the results would be $r_{ij}>c$. If violated, the translation relation used to express the scattered field as a series centred on particle $i$ as a function of that centred on particle $j$ would not be valid everywhere on particle $i$’s surface. Such a situation should invalidate the standard derivation of the STM method. However, our numerical calculations suggest here again that the STM method remains valid even beyond this more restrictive condition and is much more robust than generally assumed.

Beyond these practical conclusions, this work raises the important question: why, or indeed how does the STM method produce correct results when scatterers intrude on the no-overlap condition? While a full theoretical study is outside the scope of this work, we will conclude by discussing a possible explanation for these counter-intuitive results. We believe the considerations that have been put forward to justify why the RH is valid/invalid for a single scatterer [25,45,46] can also explain, at least in part, the “unreasonable applicability” of the STM method in the context of multiple particles. The solution to the electromagnetic problem in the exterior region is analytic and therefore defines a unique analytic continuation in the entire space, except possibly at singular points [45]. For a spheroid excited by a plane wave, the singular region is the segment between the focal points [25,45]. When expanding such a solution as a series, the resulting sums will necessarily diverge at the singular points. Since the series expansion involves VSWFs centred about the spheroid, the region of divergence becomes the entire focal sphere, namely the smallest sphere containing all the singularities of the analytic continuation of the solution. It should however be clear that this enlarged region of divergence is merely a consequence of the basis chosen for the expansion, which in this case is not well-suited to the nature of the singular region (a segment). The fact that the scattered field series diverges in a region outside the spheroid does not invalidate the solution elsewhere; further, we submit that the solution remains formally well-defined everywhere outside the spheroid as the series expansion, even when divergent, uniquely defines a finite analytic continuation. We believe a similar interpretation can be extended to multiple scatterers: what matters is the location of the singularities, not the “artificial” region of divergence resulting from the series expansion employed. As long as the singularities remain inside the scatterers, the exterior solution is always well-defined through analytic continuation, even where diverging series are present. The STM framework is therefore not affected by the divergence. Moreover, the addition of a scatterer can change the configuration of singularities, which is likely to affect the region of convergence/divergence, as observed numerically in our field maps. These considerations help justify some of the counter-intuitive results obtained numerically, but clearly further work will be needed to support these arguments with a rigorous proof. More generally perhaps, similar considerations around the divergence of series expansions and their analytic continuation should be applicable to a much wider class of differential equations in theoretical physics.

We hope these results will stimulate renewed theoretical efforts to improve the derivations of the STM framework and carefully pin down its exact range of applicability. From a practical point of view, this contribution highlights that the STM method can be applied to many relevant problems comprising densely-packed, strongly-interacting nonspherical particles, which were previously shunned because of the prevailing no-overlap condition.