#  $\mathcal{T}$-matrix method for computation of second-harmonic generation
upon optical wave scattering from clusters of arbitrary particles

Ivan Sekulic<EMAIL_ADDRESS>Department of Electronic and Electrical
Engineering, University College London, Torrington Place, London WC1E 7JE,
United Kingdom  Jian Wei You<EMAIL_ADDRESS>School of Information Science
and Engineering, Southeast University, Nanjing, Jiangsu, China  Nicolae C.
Panoiu<EMAIL_ADDRESS>Department of Electronic and Electrical
Engineering, University College London, Torrington Place, London WC1E 7JE,
United Kingdom

###### Abstract

We derive the $\mathcal{T}$-matrix formalism tailored for the numerical
analysis of second-harmonic (SH) generation from arbitrarily-shaped particles
made of centrosymmetric optical materials. First, the transfer matrix of a
single particle is computed via the extended boundary condition method, in
which the electromagnetic fields both at fundamental frequency and SH are
expanded in vector spherical wave functions and the integral formulation is
satisfied away from the surface of the scatterer. We allow for the accurate
physical description of the SH sources by taking into account both local
surface and nonlocal bulk polarization contributions to the nonlinear
polarization density source responsible for the generation of the SH signal in
a particle. This single-particle formalism is then extended to arbitrary
distributions of particles by incorporating into the formalism linear and
nonlinear electromagnetic wave scattering from the particles in the cluster.
Importantly from a practical point of view, our method can be applied to
particles of arbitrary shape made of optical materials characterized by
general frequency-dispersion relations, so that it can describe the linear and
nonlinear optical response of clusters of metallic, semiconductor, or
polaritonic particles, as well as mixtures of such particles. The approach
proposed here is faster and more memory-efficient than well-established
numerical techniques, especially in the analysis of spheroidal particles, due
to the favourable symmetries of spherical wave basis functions used in the
wave scattering analysis.

##  I Introduction

Fast and accurate full-wave computational analysis of resonant electromagnetic
response of nanoparticles is of great interest to the scientific community,
chiefly due to their widespread use in many areas of science and technology,
including plasmonic nanoantennas, remote sensing, optical metamaterials,
biology, nonlinear optics, imaging, and colloidal chemistry [6, 7, 8, 9, 3, 4,
5, 1, 2], just to name a few. Modern nanofabrication tools, such as electron
beam lithography, have made possible to produce particles of virtually any
shape, thus the ability to accurately analyze the electromagnetic properties
of nanoparticles with different morphologies is a prerequisite in the design
process of certain modern optical systems. These nanoparticles are often used
as building blocks of metamaterials, so-called meta-atoms, and can be shaped
as crosses, split-rings, cuboids, spheroids, ellipsoids, etc. [10, 11].
Combined with deep-neural networks, fast numerical solvers have dramatically
sped up the computer-aided design phase of the final structures [12, 13, 14].

Second-harmonic generation (SHG) is perhaps the most ubiquitous nonlinear
optical process encountered in modern photonics. It consists of the generation
of waves with frequency twice as large as that of an excitation (pump) optical
wave that interacts with a nonlinear optical medium [1, 15]. The full-wave
analytical solutions of this electromagnetic problem are restricted to a very
few canonical geometries, such as infinite cylinders and spheres [16, 17, 18,
19, 20, 21]. These solutions rely on the expansion of the fundamental-
frequency (FF) and second-harmonic (SH) fields into Fourier series of
cylindrical or spherical harmonic functions and the application of the
tangential field boundary conditions at the surface of the particle. In this
way, the Fourier series coefficients can be efficiently computed.

Numerical methods become indispensable when trying to go beyond these
canonical geometries. The well-known finite-element method (FEM) [22, 23] and
finite-difference time-domain (FDTD) method [24, 25] have been successfully
applied to analyze the optical response at the SH of systems containing a
small number of particles. Their drawbacks pertaining to the SHG analysis are
inherited from the linear regime and can be traced to the discretization of
the computational domain and application of absorbing boundary conditions in
order to artificially truncate it. On the other hand, electromagnetic integral
equations implicitly incorporate the radiation boundary condition at infinity
making them particularly appealing candidates for tackling complex scattering
problems [26, 27].

The Poggio–Miller–Chang–Harrington–Wu–Tsai (PMCHWT) surface integral equation
(SIE) [28, 29] is one of the most popular descriptions of the optical wave
scattering from particles in the linear regime, chiefly due to its accuracy
and resonance-free behaviour. It has recently been extended to the analysis of
SH processes in centrosymmetric media [30, 31]. The PMCHWT integral
formulation is defined on the surface of the scatterer, mapping the unknown
equivalent surface electric and magnetic currents to the tangential field
components. In its classical conforming Galerkin method of moments (MoM)
numerical discretization [32], the lowest order Rao-Wilton-Glisson (RWG)
functions [33] are used as current expansion and field testing set. These are
sub-sectional basis functions defined on the tessellation of the particle
surface (mesh), linked to the internal edges and straddling two adjacent
triangles. This restriction of the RWG expansion functions to the mesh
covering the particle surface implies that a large number of such functions
are needed to adequately interpolate the unknown surface currents and
consequently leads to MoM impedance matrices of particularly large size.
Therefore, intense efforts have been invested into developing new approaches
to accelerate execution of MoM-based codes by combining highly parallel and
performance-oriented implementations with proficient algorithms aimed at
reducing the computational complexity and execution time of the code [34, 35].
However, even if one employs state-of-the-art acceleration tools, an efficient
scattering analysis of clusters containing more than just a few nanoparticles
is seriously impeded by the prohibitively large computational resources that
are required.

In many instances of practical interest the particles under investigation
exhibit axial symmetry, e.g. ellipsoids and cylinders, or the response of a
cluster of particles with different shapes can be well approximated with one
of an assembly of spherical scatterers. In these cases certain symmetries can
be exploited by adopting vector spherical wave functions (VSWFs) as the basis
functions for the expansion of fields and currents [36, 37]. This approach is
usually coupled to the EBCM [38, 39, 40], whereby the surface integral
equations are imposed at some distance away from the boundary of the particle,
in the null-field region. In this process, a $\mathcal{T}$-matrix of the
system of particles relating the Fourier coefficients of the incident and
scattered fields is introduced [41, 42].

A key strength of the TMM consists in the fact that the $\mathcal{T}$-matrices
of particles in a cluster depend solely on the intrinsic physical properties
of the particles. Therefore, they only need to be computed once for a
particular frequency, as they do not depend on the angle and polarization of
the incident field, nor on the arrangement of particles in the cluster. This
feature makes the TMM uniquely suitable for calculations of certain physical
quantities, such as radar cross sections. Furthermore, in the calculation of
the cluster response, the multiple scattering (electromagnetic coupling)
between nanoparticles is conveniently described using the translation-addition
theorem for VSWFs [43, 44, 45], and the scattering matrix of the whole
nanoparticle aggregate is readily obtained from the single-particle
$\mathcal{T}$-matrices and translation-addition matrices. Due to this key
characteristic property of the TMM, it is alternatively called the multiple-
scattering matrix (MSM) method. The TMM has been extensively developed and
used in the linear scattering regime, including for particles of different
shapes, sizes, and electric permittivity (see [42] and references therein). It
has been extended to the analysis of the SHG in centrosymmetric media, but so
far, only the particular cases of spheres and infinite cylinders have been
considered [16, 20, 21].

In this article, we present a rigorous computational analysis based on the TMM
of optical wave scattering from particles made of centrosymmetric materials
and of arbitrary shape. Several examples of axially symmetric particles, where
the TMM performs the best, are also presented. In particular, we compare our
computational predictions with results derived using analytical methods, where
they are applicable, or otherwise with results obtained by employing
commercially available software based on the FEM. Our study is structured as
follows: in the next section, we describe the physics pertaining to the
generation of the optical waves at the SH upon the interaction of optical
fields applied at the FF with systems of nanoparticles made of centrosymmetric
materials. Then, in section III, we introduce the TMM formalism describing the
linear and nonlinear optical wave scattering from a single nanoparticle. We
extend this concept to the multiparticle case in section IV, and illustrate
how our method performs by considering a few relevant examples in section V.
In the last section we summarise the main ideas of this work.

Figure 1: A particle of arbitrary shape made of centrosymmetric optical
material described by permittivity $\epsilon_{i}$ occupying the region $V_{i}$
embedded in a background medium with permittivity $\epsilon_{e}$.

##  II System geometry and the origin of SHG in centrosymmetric media

The central goal of this study is to develop a framework in which to analyze
the SHG upon electromagnetic wave scattering from an arbitrary distribution of
$N$ particles made of centrosymmetric, non-magnetic and isotropic optical
material, occupying region $V_{i}$ and with electrical properties described by
a frequency-dependent complex permittivity function, $\epsilon_{i}$, with
$i=1,\ldots,N$. The particles are embedded in a lossless non-dispersive
background medium $V_{e}$ (usually vacuum) characterized by permittivity
$\epsilon_{e}$ (see Fig. 1). The system is excited by a time-harmonic
electromagnetic plane wave with unity amplitude,
$\\{\mathbf{E}_{inc}^{\omega},\mathbf{H}_{inc}^{\omega}\\}$, oscillating at
FF, $\omega$. In this study, we assume an $e^{-i\omega t}$ time-harmonic
variation of all fields and sources, but throughout the mathematical
derivation of our results this multiplicative harmonic factor is omitted. The
interior ($V_{i}$) and the exterior ($V_{e}$) regions are separated by an
orientable surface $\mathcal{S}$ equipped at each surface point with an
outwardly oriented unit normal vector, $\hat{\mathbf{n}}$.

The nonlinear waves generated at the SH oscillate at frequency,
$\Omega=2\omega$, and originate from two types of polarization sources. The
centrosymmetric crystal lattice of the material forbids the generation of a
local dipole polarization source inside the bulk of constituent medium.
However, due to the inversion symmetry breaking in the vicinity of the surface
of the particle, a dipole surface polarization density source is induced onto
the surface. Although in general this interfacial component of the nonlinear
polarization describes accurately enough many of the properties of the
experimentally measured SHG signal from metallic particles, for an accurate
analysis of nonlinear wave scattering from dielectric particles, a nonlocal
polarization source distribution generated inside the bulk of the particle and
originating from the induced nonlinear electric quadrupoles and nonlinear
magnetic dipoles has to be properly taken into account [46].

The nonlinear surface polarization density is induced within a few atomic
layers just underneath the surface of the particle and is related to the
electric field applied at the FF via a third-rank surface susceptibility
tensor, $\hat{\bm{\chi}}_{s}^{(2)}$ [1, 15]:

$\mathbf{P}_{s}^{\Omega}(\mathbf{r})=\epsilon_{0}\hat{\boldsymbol{\chi}}_{s}^{(2)}:\mathbf{E}^{\omega}(\mathbf{r})\mathbf{E}^{\omega}(\mathbf{r})\delta(\mathbf{r}-\mathbf{r}_{s}),$
(1)

where $\epsilon_{0}$ is the vacuum electric permittivity and the Dirac
function $\delta(\mathbf{r})$ ensures that the applied electric field is
evaluated onto the surface, $\mathbf{r}_{s}\in\mathcal{S}$ (more precisely,
the field must be evaluated at a point located on the inner side of the
particle surface). The susceptibility tensor
$\hat{\boldsymbol{\chi}}_{s}^{(2)}$ is defined at boundary points
$\mathbf{r}_{s}$, and is expressed in a coordinate system defined by a local
right-handed orthonormal triplet of vectors
$(\hat{\mathbf{n}},\hat{\mathbf{t}}_{1},\hat{\mathbf{t}}_{2})$, where
$\hat{\mathbf{t}}_{1}$ and $\hat{\mathbf{t}}_{2}$ are unit vectors tangential
to the surface at $\mathbf{r}_{s}$. Assuming isotropic surface symmetry, only
three tensor components are nonzero and independent:
$\chi_{s,\perp\perp\perp}^{(2)}$, $\chi_{s,\perp\parallel\parallel}^{(2)}$,
and
$\chi_{s,\parallel\perp\parallel}^{(2)}=\chi_{s,\parallel\parallel\perp}^{(2)}$
[15], where the symbols $\perp$ and $\parallel$ refer to the directions
perpendicular and parallel to the boundary, respectively. Taking this into
account, $\mathbf{P}_{s}^{\Omega}(\mathbf{r}_{s})$ can be expressed as [47]:

$\displaystyle\mathbf{P}_{s}^{\Omega}$
$\displaystyle(\mathbf{r}_{s})=\epsilon_{0}\left\\{\hat{\mathbf{n}}\left[\chi_{s,\perp\perp\perp}^{(2)}E_{n}^{\omega}(\mathbf{r}_{s})E_{n}^{\omega}(\mathbf{r}_{s})\right.\right.$
(2)
$\displaystyle+\left.\left.\chi_{s,\perp\parallel\parallel}^{(2)}E_{t}^{\omega}(\mathbf{r}_{s})E_{t}^{\omega}(\mathbf{r}_{s})\right]+2\chi_{s,\parallel\perp\parallel}^{(2)}E_{n}^{\omega}(\mathbf{r}_{s})\mathbf{E}_{t}^{\omega}(\mathbf{r}_{s})\right\\},$

where $\mathbf{E}_{n}^{\omega}(\mathbf{r}_{s})$ and
$\mathbf{E}_{t}^{\omega}(\mathbf{r}_{s})$ are normal and tangential components
of the linear field $\mathbf{E}^{\omega}(\mathbf{r})$ taken at the inner side
of the surface $\mathcal{S}$, at point $\mathbf{r}_{s}$, i.e.

$\displaystyle\mathbf{E}_{n}^{\omega}(\mathbf{r}_{s})$
$\displaystyle=\left[\hat{\mathbf{n}}\cdot\mathbf{E}^{\omega}(\mathbf{r})|_{\mathbf{r}=\mathbf{r}_{s}}\right]\hat{\mathbf{n}},$
(3a) $\displaystyle\mathbf{E}_{t}^{\omega}(\mathbf{r}_{s})$
$\displaystyle=-\hat{\mathbf{n}}\times\hat{\mathbf{n}}\times\mathbf{E}^{\omega}(\mathbf{r})|_{\mathbf{r}=\mathbf{r}_{s}}.$
(3b)

The nonlocal bulk polarization is induced in the interior of the particle by
the fundamental field and its spatial derivatives, and is described by the
following expression:

$\mathbf{P}_{b}^{\Omega}(\mathbf{r})=\epsilon_{0}\hat{\boldsymbol{\chi}}_{b}^{(2)}\vdots\mathbf{E}^{\omega}(\mathbf{r})\nabla\mathbf{E}^{\omega}(\mathbf{r}),\quad\mathbf{r}\in
V_{i},$ (4)

where $\hat{\boldsymbol{\chi}}_{b}^{(2)}$ is a fourth-rank bulk susceptibility
tensor. This polarization can be alternatively expressed as [48]:

$\mathbf{P}_{b}^{\Omega}=\epsilon_{0}\left[\beta\mathbf{E}^{\omega}\nabla\cdot\mathbf{E}^{\omega}+\gamma\nabla(\mathbf{E}^{\omega}\cdot\mathbf{E}^{\omega})+\delta^{\prime}(\mathbf{E}^{\omega}\cdot\nabla)\mathbf{E}^{\omega}\right],$
(5)

with $\beta,\gamma$, and $\delta^{\prime}$ being material material parameters.
The first term in this equation vanishes because we assume that for all the
particles there are no free charges, so that
$\nabla\cdot\mathbf{E}^{\omega}=0$. The relative magnitude of the third term
as compared to the second one is a question that is still debated. Thus, the
contribution of the third term increases when the field at the FF inside the
particle becomes strongly inhomogeneous, as it is the case around steep
variations of the morphology of the scatterer or nearby sharp corners where
so-called hot spots can form. Since here we are mostly concerned with axially
symmetric and smooth particles, rapid variations of the linear field are not
expected and therefore we set $\delta^{\prime}=0$ as in fact most authors do
[49]. It should be noted, however, that the term proportional to
$\delta^{\prime}$ can be accounted for in a rigorous description of generation
of SH upon optical wave scattering from particles, but the mathematical
descrition can become quite intricate even in the simple case of a single
particle [19]. Finally, the total nonlinear polarization density source is
equal to the sum of the surface and bulk components:

$\mathbf{P}^{\Omega}(\mathbf{r})=\mathbf{P}_{b}^{\Omega}(\mathbf{r})+\mathbf{P}_{s}^{\Omega}(\mathbf{r}),\quad\mathbf{r}\in
V_{i}.$ (6)

The numerical and analytical methods used in the treatment of the SH signal
are almost exclusively developed within the so-calledundepleted-pump
approximation. In this framework it is assumed that no back-coupling between
the nonlinear and linear fields occurs, i.e. the down-conversion process is
disregarded. This assumption is valid in most situations of practical interest
due to the very low energy conversion efficiency of the SHG process.

##  III The transfer-matrix method for single particle scattering

In this section, we introduce the mathematical apparatus characteristic to the
EBCM, both at the FF and SH. More specifically, we discretize the associated
integral operators in a basis consisting of VSWFs and establish a relation
between the Fourier coefficients of the series expansions in this basis of the
scattered and excitation fields. Importantly, this relation defines the
$\mathcal{T}$-matrix of the particle, which is at the heart of the TMM.

###  III.1 Linear scattering regime

The linear part of our TMM, which is used to determine the FF distribution
around an arbitrarily shaped particle, is constructed in a well-known manner
described in detail in the literature (see, for example, Ref. [50]). For the
sake of completeness and to facilitate the presentation of the nonlinear part
of the algorithm, here we will briefly outline the main steps of the
derivation, whereas the explicit formulae for the $\mathcal{T}$-matrix
elements are given in Appendix A.

We start by invoking the surface equivalence principle, which divides the
original physical problem into two problems equivalent to the initial one. The
first problem is concerned with the exterior domain $V_{e}$ and states that
the fields in this region are uniquely determined by the rotated tangential
traces of the total fields located at the boundary $\mathcal{S}=\partial
V_{i}$, the equivalent electric and magnetic currents defined as
$\mathbf{J}_{e}^{\omega}=\hat{\mathbf{n}}^{\prime}\times\mathbf{H}_{e}^{\omega}$
and
$\mathbf{M}_{e}^{\omega}=-\hat{\mathbf{n}}^{\prime}\times\mathbf{E}_{e}^{\omega}$,
and the incident fields
$\left\\{\mathbf{E}_{inc}^{\omega}(\mathbf{r}),\mathbf{H}_{inc}^{\omega}(\mathbf{r})\right\\}$.
These (non-physical) equivalent currents radiate outwardly in $V_{e}$ the
scattered fields
$\left\\{\mathbf{E}_{s}^{\omega}(\mathbf{r}),\mathbf{H}_{s}^{\omega}(\mathbf{r})\right\\}$
and exactly cancel the incident fields inside $V_{i}$, so as the equalities
$\mathbf{E}_{e}^{\omega}(\mathbf{r})=\mathbf{E}_{s}^{\omega}(\mathbf{r})+\mathbf{E}_{inc}^{\omega}(\mathbf{r})$
and
$\mathbf{H}_{e}^{\omega}(\mathbf{r})=\mathbf{H}_{s}^{\omega}(\mathbf{r})+\mathbf{H}_{inc}^{\omega}(\mathbf{r})$
hold throughout the domain $V_{e}$. These conditions are mathematically
expressed by the Stratton-Chu integral field representations [51]:

$\displaystyle\oiint_{\mathcal{S}}$
$\displaystyle\left\\{i\omega\mu_{0}\hat{\mathbf{G}}_{e}^{\omega}(\mathbf{r},\mathbf{r}^{\prime})\cdot\left[\hat{\mathbf{n}}^{\prime}\times\mathbf{H}_{e}^{\omega}(\mathbf{r}^{\prime})\right]+\left[\nabla\times\hat{\mathbf{G}}_{e}^{\omega}(\mathbf{r},\mathbf{r}^{\prime})\right]\right.$
$\displaystyle\left.\cdot\left[\hat{\mathbf{n}}^{\prime}\times\mathbf{E}_{e}^{\omega}(\mathbf{r}^{\prime})\right]\right\\}dS+\mathbf{E}_{inc}^{\omega}(\mathbf{r})=0,\quad\mathbf{r}\in
V_{i},$ (7a) $\displaystyle\oiint_{\mathcal{S}}$
$\displaystyle\left\\{i\omega\mu_{0}\hat{\mathbf{G}}_{e}^{\omega}(\mathbf{r},\mathbf{r}^{\prime})\cdot\left[\hat{\mathbf{n}}^{\prime}\times\mathbf{H}_{e}^{\omega}(\mathbf{r}^{\prime})\right]+\left[\nabla\times\hat{\mathbf{G}}_{e}^{\omega}(\mathbf{r},\mathbf{r}^{\prime})\right]\right.$
$\displaystyle\left.\cdot\left[\hat{\mathbf{n}}^{\prime}\times\mathbf{E}_{e}^{\omega}(\mathbf{r}^{\prime})\right]\right\\}dS=\mathbf{E}_{s}^{\omega}(\mathbf{r}),\quad\mathbf{r}\in
V_{e},$ (7b)

where the primed coordinates refer to the integration points on $\mathcal{S}$.
According to the first equivalent problem, the electromagnetic medium filling
$V_{i}$ is irrelevant for the field solution in $V_{e}$, so that we can in
principle replace it, to our convenience, with the same background material
occupying $V_{e}$. This way we can use the dyadic Green’s function of the
infinite space,
$\hat{\mathbf{G}}_{e}^{\omega}(\mathbf{r},\mathbf{r}^{\prime})$, defined as:

$\hat{\mathbf{G}}_{e}^{\omega}(\mathbf{r},\mathbf{r}^{\prime})=\left[\hat{\mathbf{I}}+\frac{\nabla\nabla}{(k_{e}^{\omega})^{2}}\right]\frac{e^{ik_{e}^{\omega}|\mathbf{r}-\mathbf{r}^{\prime}|}}{4\pi|\mathbf{r}-\mathbf{r}^{\prime}|},$
(8)

where $\hat{\mathbf{I}}$ is the unit dyad and
$k_{e}^{\omega}=\omega\sqrt{\epsilon_{e}\mu_{0}}$ is the wave number in the
background medium at the FF. The integral field representations (7) form the
basis of EBCM. Since the equivalent currents radiate in $V_{i}$ fields that
exactly cancel the incident fields, these equations are sometimes called the
null-field or extinction equations. A similar set of electric field integral
equations can be derived for the interior region, by exploiting the second
equivalent problem. Then, by invoking the duality principle, the magnetic
field integral counterparts of Eqs. (7) can be obtained, too [52].

The electric fields incident, scattered, and internal to the particle are
expanded in VSWFs [36, 37]:

$\displaystyle\mathbf{E}_{inc}^{\omega}(\mathbf{r})=\sum_{\nu\geq
1}\left[q_{\nu}^{\omega}\mathbf{M}_{\nu}^{(1)}(k_{e}^{\omega}\mathbf{r})+p_{\nu}^{\omega}\mathbf{N}_{\nu}^{(1)}(k_{e}^{\omega}\mathbf{r})\right],\quad\mathbf{r}\in
V_{e},$ (9a) $\displaystyle\mathbf{E}_{s}^{\omega}(\mathbf{r})=\sum_{\nu\geq
1}\left[b_{\nu}^{\omega}\mathbf{M}_{\nu}^{(3)}(k_{e}^{\omega}\mathbf{r})+a_{\nu}^{\omega}\mathbf{N}_{\nu}^{(3)}(k_{e}^{\omega}\mathbf{r})\right],\quad\mathbf{r}\in
V_{e},$ (9b) $\displaystyle\mathbf{E}_{i}^{\omega}(\mathbf{r})=\sum_{\nu\geq
1}\left[c_{\nu}^{\omega}\mathbf{M}_{\nu}^{(1)}(k_{i}^{\omega}\mathbf{r})+d_{\nu}^{\omega}\mathbf{N}_{\nu}^{(1)}(k_{i}^{\omega}\mathbf{r})\right],\quad\mathbf{r}\in
V_{i},$ (9c)

where $k_{i}^{\omega}=\omega\sqrt{\epsilon_{i}^{\omega}\mu_{0}}$ is the wave
number at the FF of the medium filling the particle. The incident and the
internal fields are expanded in a basis of regular spherical functions
$\\{\mathbf{M}_{\nu}^{(1)},\mathbf{N}_{\nu}^{(1)}\\}$, which are finite at the
origin, the corresponding expansion coefficients being
$\\{q_{\nu}^{\omega},p_{\nu}^{\omega}\\}$ and
$\\{c_{\nu}^{\omega},d_{\nu}^{\omega}\\}$, respectively. The scattered field
is expanded using radiating (outgoing) VSWFs,
$\\{\mathbf{M}_{\nu}^{(3)},\mathbf{N}_{\nu}^{(3)}\\}$, which satisfy the
Silver-Müller radiation boundary condition, the Fourier coefficients in this
case being $\\{b_{\nu}^{\omega},a_{\nu}^{\omega}\\}$. Since the incident field
is known, the corresponding Fourier coefficients are known physical
quantities. In particular, they are given in Appendix B for the specific case
of an incident linearly polarized plane wave. On the other hand, the expansion
coefficients of the internal and scattered fields are unknown and their
calculation is the main goal of the TMM.

The expanded incident, scattered, and internal magnetic fields oscillating at
frequency $\omega$ can be obtained directly from the electric fields (9), by
dint of Maxwell-Faraday equation for harmonic fields,
$\nabla\times\mathbf{E}=i\omega\mathbf{B}$:

$\displaystyle\mathbf{H}_{inc}^{\omega}(\mathbf{r})=\frac{1}{iZ_{e}}\sum_{\nu\geq
1}\left[p_{\nu}^{\omega}\mathbf{M}_{\nu}^{(1)}(k_{e}^{\omega}\mathbf{r})+q_{\nu}^{\omega}\mathbf{N}_{\nu}^{(1)}(k_{e}^{\omega}\mathbf{r})\right],\;\mathbf{r}\in
V_{e},$ (10a)
$\displaystyle\mathbf{H}_{s}^{\omega}(\mathbf{r})=\frac{1}{iZ_{e}}\sum_{\nu\geq
1}\left[a_{\nu}^{\omega}\mathbf{M}_{\nu}^{(3)}(k_{e}^{\omega}\mathbf{r})+b_{\nu}^{\omega}\mathbf{N}_{\nu}^{(3)}(k_{e}^{\omega}\mathbf{r})\right],\;\mathbf{r}\in
V_{e},$ (10b)
$\displaystyle\mathbf{H}_{i}^{\omega}(\mathbf{r})=\frac{1}{iZ_{i}^{\omega}}\sum_{\nu\geq
1}\left[d_{\nu}^{\omega}\mathbf{M}_{\nu}^{(1)}(k_{i}^{\omega}\mathbf{r})+c_{\nu}^{\omega}\mathbf{N}_{\nu}^{(1)}(k_{i}^{\omega}\mathbf{r})\right].\;\mathbf{r}\in
V_{i},$ (10c)

In these equations, $Z_{e}=\sqrt{\nicefrac{{\mu_{0}}}{{\epsilon_{e}}}}$ and
$Z_{i}^{\omega}=\sqrt{\nicefrac{{\mu_{0}}}{{\epsilon_{i}^{\omega}}}}$ are the
impedances of the background and the medium inside the particle, respectively.
In the derivation of Eqs. (10), we have relied on the curl identities that
relate the VSWFs, given in the Appendix B. Furthermore, the summation over the
multi-index $\nu$ combines summations over the orbital and azimuthal indices
$l$ and $m$, respectively, so that the summation over this index should be
understood as $\sum_{\nu\geq 1}\equiv\sum_{l\geq 1}\sum_{m=-l}^{l}$.

The dyadic Green function $\hat{\mathbf{G}}_{e}^{\omega}$ is the fundamental
solution of the vector Helmholtz equation in 3D and can be expanded in a sum
of dyadic products of VSWFs [53]:

$\displaystyle\hat{\mathbf{G}}_{e}^{\omega}(\mathbf{r},\mathbf{r}^{\prime})=$
$\displaystyle ik_{e}^{\omega}\sum_{\nu\geq
1}\left[\mathbf{M}_{\bar{\nu}}^{(3)}(k_{e}^{\omega}\mathbf{r}^{\prime})\mathbf{M}_{\nu}^{(1)}(k_{e}^{\omega}\mathbf{r})\right.$
$\displaystyle\left.+\mathbf{N}_{\bar{\nu}}^{(3)}(k_{e}^{\omega}\mathbf{r}^{\prime})\mathbf{N}_{\nu}^{(1)}(k_{e}^{\omega}\mathbf{r})\right],\quad|\mathbf{r}|<|\mathbf{r}^{\prime}|,$
(11a)
$\displaystyle\hat{\mathbf{G}}_{e}^{\omega}(\mathbf{r},\mathbf{r}^{\prime})=$
$\displaystyle ik_{e}^{\omega}\sum_{\nu\geq
1}\left[\mathbf{M}_{\bar{\nu}}^{(1)}(k_{e}^{\omega}\mathbf{r}^{\prime})\mathbf{M}_{\nu}^{(3)}(k_{e}^{\omega}\mathbf{r})\right.$
$\displaystyle\left.+\mathbf{N}_{\bar{\nu}}^{(1)}(k_{e}^{\omega}\mathbf{r}^{\prime})\mathbf{N}_{\nu}^{(3)}(k_{e}^{\omega}\mathbf{r})\right],\quad|\mathbf{r}|>|\mathbf{r}^{\prime}|,$
(11b)

where the multi-index $\bar{\nu}$ combines the indices $(l,-m)$.

In the absence of nonlinear optical interactions, the tangent components of
the electric and magnetic fields are continuous functions across the surface
of the particle, a property that is expressed as:

$\displaystyle\hat{\mathbf{n}}\times\left[\mathbf{E}_{e}^{\omega}(\mathbf{r})-\mathbf{E}_{i}^{\omega}(\mathbf{r})\right]=0,\quad\mathbf{r}\in\mathcal{S},$
(12a)
$\displaystyle\hat{\mathbf{n}}\times\left[\mathbf{H}_{e}^{\omega}(\mathbf{r})-\mathbf{H}_{i}^{\omega}(\mathbf{r})\right]=0,\quad\mathbf{r}\in\mathcal{S}.$
(12b)

Now we turn our attention to Eq. (7a), and plug in this integral equation the
VSWF-expansions of the internal electric field (9c), internal magnetic field
(10c), the incident electric field (9a), as well as the dyadic Green function
(11a) defined at points inside the particle. Equating the coefficients
multiplying regular spherical harmonics $\mathbf{M}_{\nu}^{(1)}$ and
$\mathbf{N}_{\nu}^{(1)}$ with the same multi-index $\nu$, we obtain a matrix
equation relating the Fourier coefficients of the internal and incident fields
[41]:

$\mathbf{Q}^{\omega(1,3)}\mathbf{h}^{\omega}=-\mathbf{g}^{\omega}.$ (13)

In this equation, the (unknown) internal field coefficients
$\\{c_{\nu}^{\omega},d_{\nu}^{\omega}\\}$ are assembled in the vector
$\mathbf{h}^{\omega}$, whereas the (known) incident field coefficients
$\\{q_{\nu}^{\omega},p_{\nu}^{\omega}\\}$ are stacked in the vector
$\mathbf{g}^{\omega}$. The superscript $(1,3)$ is used as a reminder that the
matrix elements of $\mathbf{Q}^{\omega(1,3)}$ are expressed as combinations of
regular and outgoing VSWFs. The specific mathematical expressions of these
matrix elements are provided in the Appendix A.

We now consider the second of Eqs. (7), from which one can establish a
relation between the expansion coefficients of the incident and scattered
fields. Analogously to the process just described, we first use the field
boundary conditions (12) to substitute in this equation the external fields
with the internal ones, and then insert the expansions of the internal and
scattered fields as well as the expansion of the dyadic Green function (11b).
Finally, in the resulting expression we equate the coefficients of the like-
terms $\mathbf{M}_{\nu}^{(3)}$ and $\mathbf{N}_{\nu}^{(3)}$ that correspond to
the same multi-index $\nu$. As a result of these mathematical manipulations we
are left with the matrix equation connecting the Fourier expansion
coefficients of the internal and scattering fields:

$\mathbf{Q}^{\omega(1,1)}\mathbf{h}^{\omega}=\mathbf{f}^{\omega}.$ (14)

In this relation, the (unknown) scattering Fourier coefficients at the FF,
$\\{b_{\nu}^{\omega},a_{\nu}^{\omega}\\}$, are collected in the vector
$\mathbf{f}^{\omega}$, and the superscript $(1,1)$ in Eq. (14) indicates that
the matrix elements of $\mathbf{Q}^{\omega(1,1)}$ contain only regular VSWFs.

Combining Eqs. (13) and (14), we arrive to the $\mathcal{T}$-matrix that maps
the Fourier coefficients of the incident field to those of the scattering
field,

$\mathbf{T}^{\omega}\mathbf{g}^{\omega}=\mathbf{f}^{\omega},$ (15)

with the $\mathcal{T}$-matrix of the particle given by the product of
Q-matrices [41]:

$\mathbf{T}^{\omega}=-\mathbf{Q}^{\omega(1,1)}\left[\mathbf{Q}^{\omega(1,3)}\right]^{-1}.$
(16)

Although in principle this formalism can be applied to scatterers of
practically any shape, because of the spherical nature of the underlying
expansion functions, the TMM is especially suitable for the analysis of
optical wave scattering from particles with axial symmetry, such as ellipsoids
and cylinders. In these cases, the numerical method converges relatively fast.
More precisely, both the near- and far-field physical quantities can be
accurately calculated using a fairly small number of VSWFs in the field
expansions. Moreover, the VSWFs can be viewed as entire-domain basis
functions, i.e. they are defined irrespective of the surface discretization
mesh, so that the constraints pertaining to the mesh generation process can be
more relaxed when compared to the MoM. Specifically, the RWG basis functions
used in the MoM are edge-based sub-sectional based, so that the underlying
surface discretization mesh has to be conformal with all pairs of adjacent
facets sharing a certain edge. On the other hand, the surface discretization
mesh needed for the evaluation of the elements of the $\mathcal{T}$-matrix can
be nonconformal as well, where adjacent triangles may not share single edges.
This is a truly significant alleviation on the mesh generation constraints.
However, if the shape of the particle is highly irregular, with sharp corners,
tips, and deep corrugations, strongly inhomogeneous near-fields are induced.
In such cases, MoM could be preferable to TMM since the flexibility provided
by the RWG basis functions allows one to accurately interpolate the rapidly
varying surface current densities [54].

###  III.2 Wave scattering at the second-harmonic

Unlike the linear electromagnetic field, whose sources are not contained in
the physical domain of interest, the nonlinear fields are generated by sources
(nonlinear polarization) distributed inside the particle or at its surface.
Hence, the electromagnetic fields oscillating at the SH frequency, $\Omega$,
must satisfy the following inhomogeneous Maxwell equations valid inside the
particle:

$\displaystyle\nabla\times\mathbf{H}_{i}^{\Omega}(\mathbf{r})+i\Omega\epsilon_{i}^{\Omega}\mathbf{E}_{i}^{\Omega}(\mathbf{r})=-i\Omega\mathbf{P}_{b}^{\Omega}(\mathbf{r}),\quad\mathbf{r}\in
V_{i},$ (17a)
$\displaystyle\nabla\times\mathbf{E}_{i}^{\Omega}(\mathbf{r})-i\Omega\mu_{0}\mathbf{H}_{i}^{\Omega}(\mathbf{r})=0,\quad\mathbf{r}\in
V_{i},$ (17b)

and the homogeneous system of equations defined in the surrounding medium:

$\displaystyle\nabla\times\mathbf{H}_{e}^{\Omega}(\mathbf{r})+i\Omega\epsilon_{e}\mathbf{E}_{e}^{\Omega}(\mathbf{r})=0,\quad\mathbf{r}\in
V_{e},$ (18a)
$\displaystyle\nabla\times\mathbf{E}_{e}^{\Omega}(\mathbf{r})-i\Omega\mu_{0}\mathbf{H}_{e}^{\Omega}(\mathbf{r})=0,\quad\mathbf{r}\in
V_{e}.$ (18b)

The external field solutions
$\left\\{\mathbf{E}_{e}^{\Omega},\mathbf{H}_{e}^{\Omega}\right\\}$ at the SH
must satisfy in the far-field region the Silver-Müller radiation condition,
whereas on the boundary of the particle the system of Maxwell equations is
supplemented with a set of nonlinear tangential field boundary conditions:

$\displaystyle\hat{\mathbf{n}}\times\left[\mathbf{E}_{e}^{\Omega}(\mathbf{r})-\mathbf{E}_{i}^{\Omega}(\mathbf{r})\right]=-\bm{\mathcal{M}}^{\Omega},\quad\mathbf{r}\in\mathcal{S},$
(19a)
$\displaystyle\hat{\mathbf{n}}\times\left[\mathbf{H}_{e}^{\Omega}(\mathbf{r})-\mathbf{H}_{i}^{\Omega}(\mathbf{r})\right]=\bm{\mathcal{J}}^{\Omega},\quad\mathbf{r}\in\mathcal{S}.$
(19b)

It can be seen from these equations that, unlike the linear case described by
Eqs. (12), the tangent fields at the SH are discontinuous. These field
discontinuities are induced by the surface nonlinear polarization, and can be
described by introducing nonlinear surface electric and magnetic currents,
$\bm{\mathcal{J}}^{\Omega}$ and $\bm{\mathcal{M}}^{\Omega}$, respectively.
Specifically, these nonlinear surface currents are related to the surface
nonlinear polarization density as [55]:

$\displaystyle\bm{\mathcal{J}}^{\Omega}(\mathbf{r}_{s})=i\Omega\hat{\mathbf{n}}\times\left[\hat{\mathbf{n}}\times\bm{\mathbf{P}}_{s}^{\Omega}(\mathbf{r}_{s})\right],$
(20a)
$\displaystyle\bm{\mathcal{M}}^{\Omega}(\mathbf{r}_{s})=\frac{1}{\epsilon_{0}}\hat{\mathbf{n}}\times\nabla_{S}\left[\hat{\mathbf{n}}\cdot\bm{\mathbf{P}}_{s}^{\Omega}(\mathbf{r}_{s})\right],$
(20b)

where the operator $\nabla_{S}$ is defined on a plane tangent to the surface
at the field evaluation point.

We follow the well-known technique of solving inhomogeneous systems of partial
differential equations by seeking the solution as the sum of the general
solution of the associated homogeneous system satisfying the homogeneous
boundary conditions (12) and a particular solution of the inhomogeneous system
obeying the nonlinear boundary conditions (19). Applying this procedure to the
Maxwell equations at the the SH (17), we look for the electromagnetic fields
$\left\\{\bar{\mathbf{E}}_{i}^{\Omega},\bar{\mathbf{H}}_{i}^{\Omega}\right\\}$
satisfying the homogeneous system of equations obtained from Eqs. (17) by
setting $\mathbf{P}_{b}^{\Omega}=0$, and a particular solution
$\left\\{\mathbf{E}_{i,p}^{\Omega},\mathbf{H}_{i,p}^{\Omega}\right\\}$ that
satisfies the inhomogeneous system (17).

As can be easily verified by direct substitution into Eqs. (17), one
particular solution to Eqs. (17) is:

$\displaystyle\mathbf{H}_{i,p}^{\Omega}(\mathbf{r})=0,\quad\mathbf{r}\in
V_{i},$ (21a)
$\displaystyle\mathbf{E}_{i,p}^{\Omega}(\mathbf{r})=-\frac{1}{\epsilon_{i}^{\Omega}}\mathbf{P}_{b}^{\Omega}(\mathbf{r}),\quad\mathbf{r}\in
V_{i},$ (21b)

so that the full solution to the system (17) is:

$\displaystyle\mathbf{H}_{i}^{\Omega}(\mathbf{r})=\bar{\mathbf{H}}_{i}^{\Omega}(\mathbf{r})+\mathbf{H}_{i,p}^{\Omega}(\mathbf{r}),\quad\mathbf{r}\in
V_{i},$ (22a)
$\displaystyle\mathbf{E}_{i}^{\Omega}(\mathbf{r})=\bar{\mathbf{E}}_{i}^{\Omega}(\mathbf{r})+\mathbf{E}_{i,p}^{\Omega}(\mathbf{r}),\quad\mathbf{r}\in
V_{i}.$ (22b)

Solutions
$\left\\{\bar{\mathbf{E}}_{i}^{\Omega},\bar{\mathbf{H}}_{i}^{\Omega}\right\\}$
and $\left\\{\mathbf{E}_{e}^{\Omega},\mathbf{H}_{e}^{\Omega}\right\\}$ of the
corresponding homogeneous systems of equations valid inside and outside the
particle, respectively, must obey the following boundary conditions on the
particle surface, $\mathcal{S}$:

$\displaystyle\hat{\mathbf{n}}\times[\mathbf{E}_{e}^{\Omega}(\mathbf{r})-\bar{\mathbf{E}}_{i}^{\Omega}(\mathbf{r})]=-\bm{\mathcal{M}}^{\Omega}(\mathbf{r})+\hat{\mathbf{n}}\times\mathbf{E}_{i,p}^{\Omega}(\mathbf{r}),\quad\mathbf{r}\in\mathcal{S},$
(23a)
$\displaystyle\hat{\mathbf{n}}\times[\mathbf{H}_{e}^{\Omega}(\mathbf{r})-\bar{\mathbf{H}}_{i}^{\Omega}(\mathbf{r})]=\bm{\mathcal{J}}^{\Omega}(\mathbf{r}),\quad\mathbf{r}\in\mathcal{S}.$
(23b)

Analogously to the integral field representations (7) derived for the fields
at the FF, following the equivalence principle, we introduce the Stratton-Chu
integral equations for the SH fields:

$\displaystyle\oiint_{\mathcal{S}}$
$\displaystyle\left\\{i\Omega\mu_{0}\hat{\mathbf{G}}_{e}^{\Omega}(\mathbf{r},\mathbf{r}^{\prime})\cdot\left[\hat{\mathbf{n}}^{\prime}\times\mathbf{H}_{e}^{\Omega}(\mathbf{r}^{\prime})\right]+\left[\nabla\times\hat{\mathbf{G}}_{e}^{\Omega}(\mathbf{r},\mathbf{r}^{\prime})\right]\right.$
$\displaystyle\left.\cdot\left[\hat{\mathbf{n}}^{\prime}\times\mathbf{E}_{e}^{\Omega}(\mathbf{r}^{\prime})\right]\right\\}dS=0,\quad\mathbf{r}\in
V_{i},$ (24a) $\displaystyle\oiint_{\mathcal{S}}$
$\displaystyle\left\\{i\Omega\mu_{0}\hat{\mathbf{G}}_{e}^{\Omega}(\mathbf{r},\mathbf{r}^{\prime})\cdot\left[\hat{\mathbf{n}}^{\prime}\times\mathbf{H}_{e}^{\Omega}(\mathbf{r}^{\prime})\right]+\left[\nabla\times\hat{\mathbf{G}}_{e}^{\Omega}(\mathbf{r},\mathbf{r}^{\prime})\right]\right.$
$\displaystyle\left.\cdot\left[\hat{\mathbf{n}}^{\prime}\times\mathbf{E}_{e}^{\Omega}(\mathbf{r}^{\prime})\right]\right\\}dS=\mathbf{E}_{e}^{\Omega}(\mathbf{r}),\quad\mathbf{r}\in
V_{e},$ (24b)

where $\hat{\mathbf{G}}_{e}^{\Omega}(\mathbf{r},\mathbf{r}^{\prime})$ is the
unbounded dyadic Green function at the SH. It is given by Eq. (8), in which we
replace $k_{e}^{\omega}\rightarrow k_{e}^{\Omega}$. In a similar way to that
described by Eqs. (11), this dyadic Green function can be expanded in VSWFs:

$\displaystyle\hat{\mathbf{G}}_{e}^{\Omega}(\mathbf{r},\mathbf{r}^{\prime})=$
$\displaystyle ik_{e}^{\Omega}\sum_{\nu\geq
1}\left[\mathbf{M}_{\bar{\nu}}^{(3)}(k_{e}^{\Omega}\mathbf{r}^{\prime})\mathbf{M}_{\nu}^{(1)}(k_{e}^{\Omega}\mathbf{r})\right.$
$\displaystyle\left.+\mathbf{N}_{\bar{\nu}}^{(3)}(k_{e}^{\Omega}\mathbf{r}^{\prime})\mathbf{N}_{\nu}^{(1)}(k_{e}^{\Omega}\mathbf{r})\right],\quad|\mathbf{r}|<|\mathbf{r}^{\prime}|,$
(25a)
$\displaystyle\hat{\mathbf{G}}_{e}^{\Omega}(\mathbf{r},\mathbf{r}^{\prime})=$
$\displaystyle ik_{e}^{\Omega}\sum_{\nu\geq
1}\left[\mathbf{M}_{\bar{\nu}}^{(1)}(k_{e}^{\Omega}\mathbf{r}^{\prime})\mathbf{M}_{\nu}^{(3)}(k_{e}^{\Omega}\mathbf{r})\right.$
$\displaystyle\left.+\mathbf{N}_{\bar{\nu}}^{(1)}(k_{e}^{\Omega}\mathbf{r}^{\prime})\mathbf{N}_{\nu}^{(3)}(k_{e}^{\Omega}\mathbf{r})\right],\quad|\mathbf{r}|>|\mathbf{r}^{\prime}|.$
(25b)

The relations expressed in Eqs. (24) are conceptually similar to the
corresponding integral equations at the FF, with the operators mapping the
rotated tangential traces of the total external SH fields
$\\{\hat{\mathbf{n}}^{\prime}\times\mathbf{E}_{e}^{\Omega},\hat{\mathbf{n}}^{\prime}\times\mathbf{H}_{e}^{\Omega}\\}$
to the SH electric field, $\mathbf{E}_{e}^{\Omega}(\mathbf{r})$, outside the
particle [as per Eq. (24b)], and to the zero field inside it – see Eq. (24a).
It is important to notice that there are no SH excitation field sources in the
exterior region, thus the total fields outside the particle are equal to the
scattered fields. As we will see later on, the SH excitation sources are
introduced in the integral formalism via boundary conditions and the
particular solution given by Eqs. (21).

We now proceed with the expansion of SH electric fields, which are solutions
of the homogeneous Maxwell equations at the SH, into Fourier series of VSWFs:

$\displaystyle\mathbf{E}_{e}^{\Omega}(\mathbf{r})=\sum_{\nu\geq
1}\Big{[}b_{\nu}^{\Omega}\mathbf{M}_{\nu}^{(3)}(k_{e}^{\Omega}\mathbf{r})+a_{\nu}^{\Omega}\mathbf{N}_{\nu}^{(3)}(k_{e}^{\Omega}\mathbf{r})\Big{]},\quad\mathbf{r}\in
V_{e},$ (26a)
$\displaystyle\bar{\mathbf{E}}_{i}^{\Omega}(\mathbf{r})=\sum_{\nu\geq
1}\Big{[}c_{\nu}^{\Omega}\mathbf{M}_{\nu}^{(1)}(k_{i}^{\Omega}\mathbf{r})+d_{\nu}^{\Omega}\mathbf{N}_{\nu}^{(1)}(k_{i}^{\Omega}\mathbf{r})\Big{]},\quad\mathbf{r}\in
V_{i},$ (26b)

with $k_{e}^{\Omega}$ and $k_{i}^{\Omega}$ representing wave numbers
associated to the exterior and the interior regions calculated at the SH
frequency, respectively. The corresponding Fourier series expansions of the SH
magnetic fields are derived from Eqs. (26) combined with Eqs. (17b) and (18b),
and can be written as:

$\displaystyle\mathbf{H}_{e}^{\Omega}(\mathbf{r})=\frac{1}{iZ_{e}}\sum\limits_{\nu\geq
1}\Big{[}a_{\nu}^{\Omega}\mathbf{M}_{\nu}^{(3)}(k_{e}^{\Omega}\mathbf{r})+b_{\nu}^{\Omega}\mathbf{N}_{\nu}^{(3)}(k_{e}^{\Omega}\mathbf{r})\Big{]},~{}\mathbf{r}\in
V_{e},$ (27a)
$\displaystyle\bar{\mathbf{H}}_{i}^{\Omega}(\mathbf{r})=\frac{1}{iZ_{i}^{\Omega}}\sum\limits_{\nu\geq
1}\Big{[}d_{\nu}^{\Omega}\mathbf{M}_{\nu}^{(1)}(k_{i}^{\Omega}\mathbf{r})+c_{\nu}^{\Omega}\mathbf{N}_{\nu}^{(1)}(k_{i}^{\Omega}\mathbf{r})\Big{]},~{}\mathbf{r}\in
V_{i}.$ (27b)

Regular spherical harmonics are used in the expansion of the internal SH
fields
$\left\\{\bar{\mathbf{E}}_{i}^{\Omega},\bar{\mathbf{H}}_{i}^{\Omega}\right\\}$,
as they must be finite at the origin, while the external SH fields
$\left\\{\mathbf{E}_{e}^{\Omega},\mathbf{H}_{e}^{\Omega}\right\\}$ are
expanded in the basis of outgoing vector spherical harmonics. The Fourier
expansion coefficients $\\{c_{\nu}^{\Omega},d_{\nu}^{\Omega}\\}$ for the
internal fields and $\\{b_{\nu}^{\Omega},a_{\nu}^{\Omega}\\}$ for the external
ones are unknown physical quantities and will be determined by solving a
certain linear system of equations that we will construct in what follows.

In the construction of the $\mathcal{T}$-matrix at the SH, which relates the
expansion coefficients of the excitation and external (scattering) fields, we
employ the SH integral relation (24a) valid in the interior of the particle in
which we replace the external field traces
$\left\\{\hat{\mathbf{n}}^{\prime}\times\mathbf{E}_{e}^{\Omega},\hat{\mathbf{n}}^{\prime}\times\mathbf{H}_{e}^{\Omega}\right\\}$
with the homogeneous internal field traces
$\left\\{\hat{\mathbf{n}}^{\prime}\times\bar{\mathbf{E}}_{i}^{\Omega},\hat{\mathbf{n}}^{\prime}\times\bar{\mathbf{H}}_{i}^{\Omega}\right\\}$.
This latter step is performed by taking into account the boundary conditions
(23). After these mathematical manipulations, we arrive to the following
integral equation:

$\displaystyle\oiint_{\mathcal{S}}\left\\{i\Omega\mu_{0}\hat{\mathbf{G}}_{e}^{\Omega}(\mathbf{r},\mathbf{r}^{\prime})\cdot\left[\hat{\mathbf{n}}^{\prime}\times\bar{\mathbf{H}}_{i}^{\Omega}(\mathbf{r}^{\prime})\right]\right.$
$\displaystyle+\left.\left[\nabla\times\hat{\mathbf{G}}_{e}^{\Omega}(\mathbf{r},\mathbf{r}^{\prime})\right]\cdot\left[\hat{\mathbf{n}}^{\prime}\times\bar{\mathbf{E}}_{i}^{\Omega}(\mathbf{r}^{\prime})\right]\right\\}dS=-i\Omega\mu_{0}\oiint_{\mathcal{S}}\hat{\mathbf{G}}_{e}^{\Omega}(\mathbf{r},\mathbf{r}^{\prime})\cdot\bm{\mathcal{J}}^{\Omega}(\mathbf{r}^{\prime})dS$
$\displaystyle+\oiint_{\mathcal{S}}\left[\nabla\times\hat{\mathbf{G}}_{e}^{\Omega}(\mathbf{r},\mathbf{r}^{\prime})\right]\cdot\left[\bm{\mathcal{M}}^{\Omega}(\mathbf{r}^{\prime})-\hat{\mathbf{n}}^{\prime}\times\mathbf{E}_{i,p}^{\Omega}(\mathbf{r}^{\prime})\right]dS,\quad\mathbf{r}\in
V_{i},$ (28)

where the r.h.s. of this equation represents the contribution of the SH
polarization sources. Substituting in the equation above the VSWF-expansions
of the internal homogeneous fields at the SH, described by Eqs. (26b) and
(27b), as well as the expansion in VSWFs of the dyadic Green function at the
SH and valid for the points inside the target, given by Eq. (25a), and
equating the coefficients multiplying regular spherical harmonics with the
same multi-index, we are left with a linear system of equations represented in
the matrix form as:

$\mathbf{Q}^{\Omega(1,3)}\mathbf{h}^{\Omega}=\mathbf{g}^{\Omega(3)}.$ (29)

The matrix elements $\mathbf{Q}^{\Omega(1,3)}$ are of a form similar to that
of the FF matrix elements $\mathbf{Q}^{\omega(1,3)}$, but computed at the
frequency $\Omega$. Their specific expressions are given in Appendix A.
Moreover, the unknown Fourier coefficients of the internal homogeneous fields,
$\\{c_{\nu}^{\Omega},d_{\nu}^{\Omega}\\}$, are stacked in the vector
$\mathbf{h}^{\Omega}$. The nonlinear excitation source components are
collected in the vector $\mathbf{g}^{\Omega(3)}$ and are defined by surface
integrals of products of tangential components of the nonlinear current
sources and outgoing VSWFs. The specific formulae of the components of
$\mathbf{g}^{\Omega(3)}$ are provided in Appendix A, too. Note that
$\mathbf{g}^{\Omega(3)}$ can be computed once the fields at the FF are
determined, so that as far as the calculation of the nonlinear fields is
concerned the vector $\mathbf{g}^{\Omega(3)}$ should be viewed as a known
physical quantity.

In order to evaluate the Fourier coefficients
$\\{b_{\nu}^{\Omega},a_{\nu}^{\Omega}\\}$ of the external fields at the SH, we
use Eq. (24b) from which again the external fields are eliminated with the
help of the boundary conditions (23). The resulting equation can be cast in
the form:

$\displaystyle\oiint_{\mathcal{S}}\left\\{i\Omega\mu_{0}\hat{\mathbf{G}}_{e}^{\Omega}(\mathbf{r},\mathbf{r}^{\prime})\cdot\left[\hat{\mathbf{n}}^{\prime}\times\bar{\mathbf{H}}_{i}^{\Omega}(\mathbf{r}^{\prime})\right]\right.$
$\displaystyle+\left.\left[\nabla\times\hat{\mathbf{G}}_{e}^{\Omega}(\mathbf{r},\mathbf{r}^{\prime})\right]\cdot\left[\hat{\mathbf{n}}^{\prime}\times\bar{\mathbf{E}}_{i}^{\Omega}(\mathbf{r}^{\prime})\right]\right\\}dS=\mathbf{E}_{e}^{\Omega}(\mathbf{r})-i\Omega\mu_{0}\oiint_{\mathcal{S}}\hat{\mathbf{G}}_{e}^{\Omega}(\mathbf{r},\mathbf{r}^{\prime})\cdot\bm{\mathcal{J}}^{\Omega}(\mathbf{r}^{\prime})dS$
$\displaystyle+\oiint_{\mathcal{S}}\left[\nabla\times\hat{\mathbf{G}}_{e}^{\Omega}(\mathbf{r},\mathbf{r}^{\prime})\right]\cdot\left[\bm{\mathcal{M}}^{\Omega}(\mathbf{r}^{\prime})-\hat{\mathbf{n}}^{\prime}\times\mathbf{E}_{i}^{\Omega,p}(\mathbf{r}^{\prime})\right]dS,\quad\mathbf{r}\in
V_{e}.$ (30)

After replacing the external and the internal fields with their corresponding
spherical harmonic expansions, then substituting Eq. (25b) for the dyadic
Green function $\hat{\mathbf{G}}_{e}^{\Omega}$ in the outer region, and
finally collecting the Fourier coefficients multiplying radiating spherical
harmonics $\mathbf{M}_{\nu}^{(3)}$ and $\mathbf{N}_{\nu}^{(3)}$ of the same
multi-index $\nu$, we are left with the following matrix equation:

$\mathbf{f}^{\Omega}+\mathbf{g}^{\Omega(1)}=\mathbf{Q}^{\Omega(1,1)}\mathbf{h}^{\Omega}.$
(31)

Here, the vector $\mathbf{f}^{\Omega}$ contains the unknown coefficients of
the SH external fields, $\\{b_{\nu}^{\Omega},a_{\nu}^{\Omega}\\}$, and
$\mathbf{g}^{\Omega(1)}$ is the source vector with components expressed as
surface integrals of products involving tangential components of the nonlinear
current sources and regular VSWFs.

Finally, we eliminate the vector $\mathbf{h}^{\Omega}$ between Eqs. (29) and
(31), and so arrive to the following matrix equation relating the vector
containing the expansion coefficients of the external (scattering) field,
$\mathbf{f}^{\Omega}$, and the vectors $\mathbf{g}^{\Omega(3)}$ and
$\mathbf{g}^{\Omega(1)}$ holding the expansion coefficients of the nonlinear
excitation sources:

$\mathbf{f}^{\Omega}=-\mathbf{g}^{\Omega(1)}+\mathbf{Q}^{\Omega(1,1)}[\mathbf{Q}^{\Omega(1,3)}]^{-1}\mathbf{g}^{\Omega(3)}.$
(32)

This formula can be used to define the $\mathcal{T}$-matrix of the particle
that characterizes the nonlinear scattering wave interaction at the SH:

$\mathbf{T}^{\Omega}=\mathbf{Q}^{\Omega(1,1)}[\mathbf{Q}^{\Omega(1,3)}]^{-1}.$
(33)

Similarly to the TMM formulation at the FF, the SH TMM developed here can in
principle be used to analyze nonlinear wave scattering from particles with
arbitrary morphologies. One important idea related to the nonlinear field
calculations, which we want to highlight here, is that in this two-step
procedure the accurate evaluation of the linear near-field on the inner side
of the particle surface is of paramount importance, as these near-fields fully
determine the polarization sources at the SH. This aspect of the TMM should be
considered carefully especially when studying metallic particles or particles
with surface that deviates significantly from a sphere, as it is well known
that in these cases the TMM calculation of the linear near-fields converges
slowly.

##  IV Extension to the multiparticle systems

In this section, we generalize the formalism introduced in the previous
section, valid for a single particle, to an arbitrary distribution of
particles. This extension is conveniently achieved by dint of translation-
addition matrices [43, 44, 45], which are used to capture the multiple wave
scattering nature of the physical process investigated.

###  IV.1 Linear wave scattering

The physical system under investigation is schematically depicted in Fig. 2.
It consists of a general distribution of $N$ particles whose shape and
material parameters are given but arbitrary otherwise. We associate to each
particle a (local) coordinate system with origin at
$\mathcal{O}_{n},n=1,\ldots,N$, and describe the entire cluster of particles
in a (global) coordinate system with origin at $\mathcal{O}$, chosen in such a
manner that the particles are arranged around it in some loosely defined
balanced way.

Figure 2: Schematic illustration of a distribution of particles, the global
coordinate system with origin at $\mathcal{O}$, and local coordinate systems
attached to each of the particles.

We can express the field exciting the nth particle, $n=1,\ldots,N$, as the sum
of the incident field at the location of the particle and the fields scattered
from all the other particles (after possibly being scattered multiple times
among the particles):

$\displaystyle\mathbf{E}_{n,ex}^{\omega}$
$\displaystyle(\mathbf{r})=\left[\mathbf{M}_{\nu}^{(1)}(k_{e}^{\omega}\mathbf{r}),\mathbf{N}_{\nu}^{(1)}(k_{e}^{\omega}\mathbf{r})\right]\mathbf{g}^{\omega}$
(34) $\displaystyle+\sum\limits_{\begin{subarray}{c}k=1\\\ k\neq
n\end{subarray}}^{N}\left[\mathbf{M}_{\nu}^{(3)}(k_{e}^{\omega}\mathbf{r}_{k}),\mathbf{N}_{\nu}^{(3)}(k_{e}^{\omega}\mathbf{r}_{k})\right]\mathbf{f}_{k}^{\omega},\quad
n=1,\ldots,N.$

To simplify the notation, in this equation and in what follows the summation
over the multi-index $\nu$ is omitted. This field can be expanded in the local
coordinate system associated to the nth particle as:

$\displaystyle\mathbf{E}$
${}_{n,ex}^{\omega}(\mathbf{r}_{n})=\left[\mathbf{M}_{\nu}^{(1)}(k_{e}^{\omega}\mathbf{r}_{n}),\mathbf{N}_{\nu}^{(1)}(k_{e}^{\omega}\mathbf{r}_{n})\right]\bm{\beta}_{n,0}^{\omega}\mathbf{g}^{\omega}$
(35) $\displaystyle+\sum\limits_{\begin{subarray}{c}k=1\\\ k\neq
n\end{subarray}}^{N}\left[\mathbf{M}_{\nu}^{(1)}(k_{e}^{\omega}\mathbf{r}_{n}),\mathbf{N}_{\nu}^{(1)}(k_{e}^{\omega}\mathbf{r}_{n})\right]\bm{\alpha}_{n,k}^{\omega}\mathbf{f}_{k}^{\omega},~{}n=1,\ldots,N.$

Here, the matrix $\bm{\alpha}_{n,k}^{\omega}$ is the FF irregular translation-
addition operator mapping the radiative VSWFs defined in the local coordinate
system with origin at $\mathcal{O}_{k}$ to the regular VSWFs expressed in the
local coordinate system with origin at $\mathcal{O}_{n}$ (see Fig. 2).
Similarly, $\bm{\beta}_{n,0}^{\omega}$ is the regular translation-addition
matrix computed at the FF $\omega$, which maps the regular VSWFs defined in
the global coordinate system with origin at $\mathcal{O}$ to their expansion
in the local coordinate system with origin at $\mathcal{O}_{n}$, associated to
the $n$th particle. The total field incident on the $n$th particle,
$\mathbf{E}_{n,ex}^{\omega}(\mathbf{r}_{n})$, can itself be expanded into a
series of regular VSWFs expressed in the local coordinate system with origin
at $\mathcal{O}_{n}$ and weighted with coefficients stacked in a vector
denoted with $\mathbf{e}_{n}^{\omega}$. Therefore, we are left with the
following equation containing only Fourier expansion coefficients:

$\mathbf{e}_{n}^{\omega}=\bm{\beta}_{n,0}^{\omega}\mathbf{g}^{\omega}+\sum\limits_{\begin{subarray}{c}k=1\\\
k\neq
n\end{subarray}}^{N}\bm{\alpha}_{n,k}^{\omega}\mathbf{f}_{k}^{\omega},\quad
n=1,\ldots,N.$ (36)

Note that $\mathbf{e}_{n}^{\omega}$ describes the field incident onto the
$n$th particle of the cluster, so that it plays the rôle of
$\mathbf{g}_{n}^{\omega}$ for an isolated individual particle labeled by $n$,
in the absence of other particles. The key difference between these two
vectors is that $\mathbf{e}_{n}^{\omega}$ consists of both the field directly
incident onto the cluster, which reaches the particle $n$ without being
scattered by any other particle in the cluster, and the field that reaches the
particle $n$ after being scattered (possibly multiple times) by the other
particles in the cluster; $\mathbf{g}_{n}^{\omega}$, on the other hand,
consists simply of the incident field exciting the cluster. With this in mind,
it should be clear that the relation
$\mathbf{f}_{n}^{\omega}=\mathbf{T}_{n}^{\omega}\mathbf{e}_{n}^{\omega}$
holds, where $\mathbf{T}_{n}^{\omega}$ is the $\mathcal{T}$-matrix of the
$n$th particle.

Multiplying Eqs. (36) with the $\mathcal{T}$-matrices of individual scatterers
$\left\\{\mathbf{T}_{1}^{\omega},\ldots,\mathbf{T}_{N}^{\omega}\right\\}$,
respectively, which are precomputed and stored beforehand, we are left with a
linear system of equations whose solution gives us the scattering field
Fourier coefficients:

$\mathbf{f}_{n}^{\omega}=\mathbf{T}_{n}^{\omega}\bm{\beta}_{n,0}^{\omega}\mathbf{g}^{\omega}+\sum\limits_{\begin{subarray}{c}k=1\\\
k\neq
n\end{subarray}}^{N}\mathbf{T}_{n}^{\omega}\bm{\alpha}_{n,k}^{\omega}\mathbf{f}_{k}^{\omega},\quad
n=1,\ldots,N.$ (37)

This linear system of equations can be cast in the following matrix form [41]:

$\displaystyle\begin{bmatrix}\textbf{I}&-\textbf{T}_{1}^{\omega}\bm{\alpha}_{1,2}^{\omega}&\cdots&-\textbf{T}_{1}^{\omega}\bm{\alpha}_{1,N}^{\omega}\\\
-\textbf{T}_{2}^{\omega}\bm{\alpha}_{2,1}^{\omega}&\textbf{I}&\cdots&-\textbf{T}_{2}^{\omega}\bm{\alpha}_{2,N}^{\omega}\\\
\vdots&\vdots&\ddots&\vdots\\\
-\textbf{T}_{N}^{\omega}\bm{\alpha}_{N,1}^{\omega}&-\textbf{T}_{N}^{\omega}\bm{\alpha}_{N,2}^{\omega}&\cdots&\textbf{I}\end{bmatrix}\begin{bmatrix}\textbf{f}_{1}^{\omega}\\\
\textbf{f}_{2}^{\omega}\\\ \vdots\\\ \textbf{f}_{N}^{\omega}\end{bmatrix}=$
$\displaystyle\begin{bmatrix}\textbf{T}_{1}^{\omega}\bm{\beta}_{1,0}^{\omega}\textbf{g}^{\omega}\\\
\textbf{T}_{2}^{\omega}\bm{\beta}_{2,0}^{\omega}\textbf{g}^{\omega}\\\
\vdots\\\
\textbf{T}_{N}^{\omega}\bm{\beta}_{N,0}^{\omega}\textbf{g}^{\omega}\\\
\end{bmatrix},$ (38)

or in the more condensed form,

$\mathbf{S}^{\omega}\mathbf{F}^{\omega}=\mathbf{G}^{\omega}.$ (39)

The matrix $\mathbf{S}^{\omega}$ is the scattering matrix of the system of
particles. It contains all the information regarding the geometry and material
parameters of the cluster. The unknown scattering field coefficients
corresponding to all of the particles are collected in the vector
$\mathbf{F}^{\omega}$ and are retrieved by solving the linear system (39). The
corresponding expansion coefficients associated to the fields internal to each
of the particles are obtained with the help of Eq. (14). After the scattering
and internal field expansion coefficients are computed for all of the
particles the corresponding fields can be directly obtained from Eqs. (9b) and
(9c) (electric fields) and (10b) and (10c) (magnetic fields). Finally, the
scattering cross section spectrum of the cluster of particles,
$\sigma_{sc}^{FF}(\omega)$, is easily obtained using the scattering field
coefficients of all of the particles:

$\sigma_{sc}^{FF}(\omega)=\frac{1}{(k_{e}^{\omega})^{2}}\sum\limits_{n,k=1}^{N}\mathfrak{Re}\\{{\textbf{f}_{n}^{\omega}}^{{\dagger}}\bm{\beta}_{n,k}^{\omega}\textbf{f}_{k}^{\omega}\\},$
(40)

where the symbol “†” denotes Hermitian conjugation and $\mathfrak{Re}\\{z\\}$
is the real part of the complex number $z$. To obtain this expression one
integrates over the scattering angles the differential scattering cross-
section of each particle then sums the results over the particles (a detailed
derivation of this formula can be found in Ref. [56]).

###  IV.2 Wave scattering at the second-harmonic

The analysis of the wave scattering at the SH follows the same main steps as
in the case of the linear wave scattering we just discussed. The main
difference between the linear and nonlinear wave scattering amounts to the
fact that the sources of the SH field are not located at an infinite distance
from the cluster of particles but inside the particles and at their surface.
To account for this key difference, we add to the right hand side of the
integral equation (24a) the excitation field at the SH, representing the
source field for the $n$th particle. This excitation field consists of the sum
of the fields scattered from all the other particles.

This excitation field at the SH can be expanded into series of regular VSWFs
around the origin $\mathcal{O}_{n}$ associated to the $n$th particle, so that
we can establish a relation between the expansion coefficients of the SH
internal homogeneous fields corresponding to the $n$th particle and the
expansion coefficients of the total nonlinear field impinging onto the
particle:

$\mathbf{Q}_{n}^{\Omega(1,3)}\mathbf{h}_{n}^{\Omega}=\mathbf{g}_{n}^{\Omega(3)}-\sum\limits_{\begin{subarray}{c}k=1\\\
k\neq
n\end{subarray}}^{N}\bm{\alpha}_{n,k}^{\Omega}\mathbf{f}_{k}^{\Omega},\quad
n=1,\ldots,N,$ (41)

In these relations, the translation-addition matrices
$\bm{\alpha}_{n,k}^{\Omega}$ are evaluated at the SH frequency, $\Omega$.
Substituting the internal SH field expansion coefficients
$\mathbf{h}_{n}^{\Omega}$ to Eq. (31) governing the Fourier coefficients of
the scattered field at the SH, and iterating through all of the particles, we
obtain:

$\mathbf{f}^{\Omega}_{n}=-\mathbf{g}^{\Omega(1)}_{n}+\mathbf{T}^{\Omega}_{n}\mathbf{g}^{\Omega(3)}_{n}-\sum\limits_{\begin{subarray}{c}k=1\\\
k\neq
n\end{subarray}}^{N}\mathbf{T}_{n}^{\Omega}\bm{\alpha}_{n,k}^{\Omega}\mathbf{f}_{k}^{\Omega},\quad
n=1,\ldots,N,$ (42)

These relations can be expressed in matrix form as:

$\displaystyle\begin{bmatrix}\textbf{I}&\textbf{T}_{1}^{\Omega}\bm{\alpha}_{1,2}^{\Omega}&\cdots&\textbf{T}_{1}^{\Omega}\bm{\alpha}_{1,N}^{\Omega}\\\
\textbf{T}_{2}^{\Omega}\bm{\alpha}_{2,1}^{\Omega}&\textbf{I}&\cdots&\textbf{T}_{2}^{\Omega}\bm{\alpha}_{2,N}^{\Omega}\\\
\vdots&\vdots&\ddots&\vdots\\\
\textbf{T}_{N}^{\Omega}\bm{\alpha}_{N,1}^{\Omega}&\textbf{T}_{N}^{\Omega}\bm{\alpha}_{N,2}^{\Omega}&\cdots&\textbf{I}\end{bmatrix}\begin{bmatrix}\textbf{f}_{1}^{\Omega}\\\
\textbf{f}_{2}^{\Omega}\\\ \vdots\\\ \textbf{f}_{N}^{\Omega}\end{bmatrix}=$
$\displaystyle\begin{bmatrix}-\textbf{g}_{1}^{\Omega(1)}+\textbf{T}_{1}^{\Omega}\textbf{g}_{1}^{\Omega(3)}\\\
-\textbf{g}_{2}^{\Omega(1)}+\textbf{T}_{2}^{\Omega}\textbf{g}_{2}^{\Omega(3)}\\\
\vdots\\\
-\textbf{g}_{N}^{\Omega(1)}+\textbf{T}_{N}^{\Omega}\textbf{g}_{N}^{\Omega(3)}\\\
\end{bmatrix},$ (43)

or, more concisely,

$\mathbf{S}^{\Omega}\mathbf{F}^{\Omega}=\mathbf{G}^{\Omega}.$ (44)

Similarly to the wave scattering at the FF, $\mathbf{S}^{\Omega}$ represents
the system scattering matrix at the SH, the vector $\mathbf{F}^{\Omega}$ holds
the Fourier coefficients of the scattering field corresponding to all of the
particles, and the vector $\mathbf{G}^{\Omega}$ contains the information on
the nonlinear polarization sources responsible for the generation of the
optical field at the SH. After the SH scattering field expansion coefficients
corresponding to all of the particles are computed by solving Eq. (44), for
each particle one can compute via Eq. (31) the expansion coefficients
corresponding to the internal homogeneous SH field. The associated scattering
and internal homogeneous fields are computed, respectively, from Eqs. (26a)
and (26b) (electric fields) and Eqs. (27a) and (27b) (magnetic fields).
Finally, in order to compute the total fields inside each particle, the
particular (inhomogeneous) solution defined by Eqs. (21) must be added to the
homogeneous one. This is conveniently done by expanding the product of linear
electric fields inside the particle into Fourier series of VSWFs [36]. The SH
scattering cross-section spectra are then calculated in an analogous fashion
to the procedure used at the FF, as per (40), with all the physical quantities
being evaluated at SH frequency.

##  V Applications and discussion

In this section we will illustrate how the numerical method developed in the
preceding section can be used to compute the scattering cross-section spectra
of of systems of particles with different morphology, the calculations being
performed both at the fundamental and second-harmonic frequencies. In the case
in which the particles are spherical, we compare our results with analytical
solutions [19] (for a single sphere), or predictions of the TMM with
analytically computed $\mathcal{T}$-matrices (for a multiple sphere system).
Because in this latter case the analytic calculation of the matrix elements is
based on the Mie theory [57], we call it Mie-TMM [21]. In the case of
particles with non-spherical shape, we validate our results against results
obtained using CST Studio [58], a commercial software based on the FEM. We
consider particles made of two centrosymmetric optical materials widely used
in nanophotonics, namely silicon (Si) and gold (Au). Our choice was guided by
the fact that both these optical materials present challenging testing
conditions for any numerical method, namely silicon particles have large
electric permittivity and hence scatter significantly the incident field,
whereas (plasmonic) particles made of Au possess plasmon resonances that
induce strongly inhomogeneous, enhanced near-fields.

The cluster of particles is assumed to be embedded in vacuum, whereas the
excitation field is a $\hat{\bm{\theta}}$ polarized plane wave whose
propagation direction is defined by the incident angles $\theta_{inc}=\pi/4$
and $\phi_{inc}=\pi/2$. The frequency-dependent permittivity of silicon and
gold is interpolated from the experimental data provided in Refs. [59] and
[60], respectively. The tensor components of the nonlinear surface
susceptibility $\bm{\hat{\chi}_{s}^{(2)}}$ of gold, as well as its bulk
nonlinear susceptibility $\gamma$, are calculated using the free-electron
hydrodynamic model [49]:

$\displaystyle\chi_{s,\perp\perp\perp}^{(2)}=-\frac{a}{4}\left[\epsilon_{r}(\omega)-1\right]\frac{e}{m\omega^{2}},$
(45a)
$\displaystyle\chi_{s,\parallel\perp\parallel}^{(2)}=\chi_{s,\parallel\parallel\perp}^{(2)}=-\frac{b}{2}\left[\epsilon_{r}(\omega)-1\right]\frac{e}{m\omega^{2}},$
(45b)
$\displaystyle\gamma=-\frac{d}{8}\left[\epsilon_{r}(\omega)-1\right]\frac{e}{m\omega^{2}},$
(45c)

where the Rudnick-Stern coefficients [61] are $a=1$, $b=-1$, and $c=1$ [49],
whereas $e$ and $m$ denote the charge and the mass of the electron,
respectively. The corresponding nonlinear susceptibilities of silicon are
assumed to be constant in the whole frequency range of interest, with values
equal to
$\chi_{s,\perp\perp\perp}^{(2)}=$65\text{\times}{10}^{-19}\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{V}}^{-1}$$,
$\chi_{s,\parallel\perp\parallel}^{(2)}=\chi_{s,\parallel\parallel\perp}^{(2)}=$3.5\text{\times}{10}^{-19}\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{V}}^{-1}$$,
and
$\gamma=$1.3\text{\times}{10}^{-19}\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{V}}^{-1}$$
[62]. For both optical materials, gold and silicon, we set
$\chi_{s,\perp\parallel\parallel}^{(2)}=0$, which is a commonly used
assumption [63],[64].

###  V.1 Convergence tests

Figure 3: Relative scattering cross-section error of two Si spheres at FF (a)
and SH (b).

We begin our investigation of the properties of the numerical method developed
in this study with an analysis of its convergence characteristics. For this
purpose, we choose two representative numerical examples and monitor the
evolution of the relative scattering cross-section error with the maximum
expansion order, $l_{max}$. This analysis was performed for several values of
the frequency of the incident wave. In the first example we considered two
silicon nanospheres with radii $R_{1}=$250\text{\,}\mathrm{nm}$$ and
$R_{2}=$200\text{\,}\mathrm{nm}$$ centered at
$O_{1}($0\text{\,}\mathrm{,}$$0\text{\,}\mathrm{,}$$0\text{\,}\mathrm{)}$$ and
$O_{2}($0\text{\,}\mathrm{,}$$0\text{\,}\mathrm{,}$$550\text{\,}\mathrm{nm}$)$.
In the second example the particle system consists of two spheroids, a prolate
spheroid with equatorial and polar radii $a_{1}=$200\text{\,}\mathrm{nm}$$ and
$c_{1}=$300\text{\,}\mathrm{nm}$$, respectively, and an oblate spheroid with
radii $a_{2}=$320\text{\,}\mathrm{nm}$$ and $c_{2}=$230\text{\,}\mathrm{nm}$$.
The two spheroids are made of gold and centered at
$O_{1}($0\text{\,}\mathrm{,}$$-350\text{\,}\mathrm{nm}$,$0\text{\,}\mathrm{)}$$
and
$O_{2}($0\text{\,}\mathrm{,}$$350\text{\,}\mathrm{nm}$,$0\text{\,}\mathrm{)}$$,
respectively.

Figure 4: Relative scattering cross-section error of two Au spheroids at FF
(a) and SH (b).

We define the relative scattering cross-section error as:

$\displaystyle\mathcal{E}^{\omega,\Omega}=\left|\frac{\sigma_{l_{max}}^{\omega,\Omega}-\sigma_{REF}^{\omega,\Omega}}{\sigma_{REF}^{\omega,\Omega}}\right|.$
(46)

The reference scattering cross-section $\sigma_{REF}^{\omega,\Omega}$ is
calculated within the Mie-TMM formalism [21] for the case of nanospheres, and,
for the case of spheroids, with the TMM method described in this work and
expansion order $l_{max}=10$. According to the data plotted in Figs. 3 and 4,
a reduction in the relative scattering cross-section error, both in FF and SH
regimes, is observed as the number of spherical harmonics included in the
Fourier expansion is increased. This is the characteristic trend observed when
a numerical method converges. We also note a general trend, more pronounced in
the case of scattering cross-section error evaluated at the FF, namely that
the error is smaller for larger values of the wavelength of the incident wave.
This is explained by the fact that as the wavelength increases the wave
scattering process approaches the quasi-static regime described by the first
order of the series expansion. This conclusion is not necessarily valid at the
SH, as in this case the nonlinear source distribution can be markedly affected
by the excitation of optical resonances of the particles, a phenomenon that is
strongly dependent on the wavelength.

###  V.2 Wave scattering from nanospheres

In this numerical example, we choose a nanosphere made of either gold or
silicon and of radius equal to $R=$50\text{\,}\mathrm{nm}$$ in both cases.
Nanospheres of this characteristic size are commonly used in different
applications, including plasmonic sensing and colloidal chemistry. Although it
is a simple physical system, it has the advantage that analytic solutions to
the linear and nonlinear scattering problems exist, so that it provides an
effective validation test of our numerical method. Thus, for this system we
calculated the spectra of the scattering cross-section, at both the FF and SH,
in a spectral domain of practical interest. We meshed the surface of each
nanosphere into 2000 triangles, which was enough to reach convergence. The
scattering cross-section spectra were obtained with the method introduced in
this study, where the corresponding surface integrals were computed
numerically over the boundary tessellation, and the results were compared to
the conclusions derived from the analytical Mie theory [19]. In both
approaches, numerical and analytical, we used a maximum expansion order of
$l_{max}=6$, which corresponds to a total of 96 VSWFs per sphere at the FF and
SH.

Figure 5: Scattering cross-sections at FF (left panels) and SH (right panels)
for a gold (top panels) and silicon (bottom panels) sphere with radius
$R=$50\text{\,}\mathrm{nm}$$.

The spectra computed using this procedure are plotted in Fig. 5. It is
interesting to notice that the FF and SH scattering spectra of the gold
particle exhibit a narrow spectral peak at the incident wavelength
$\lambda_{inc}=$520\text{\,}\mathrm{nm}$$, Figs. 5(a) and 5(b). These peaks
can be attributed to the plasmonic resonances of a gold nanosphere. On the
other hand, Mie resonances, characteristic to dielectric particles, are not
observed for the silicon nanosphere in the frequency range considered here, as
for particles with size considered in this example they exist for shorter
wavelengths. Furthermore, according to the data presented in Fig. 5, there is
a good agreement, both at fundamental and second-harmonic frequencies, between
scattering cross-section spectra obtained with the analytical and numerical
procedures.

Figure 6: Electric field amplitude in the $y=0$ median plane of a Au (top
panels) and Si (bottom panels) sphere with radius of $50\text{\,}\mathrm{nm}$,
at FF (left panels) and SH (right panels).

In Fig. 6 we show the spatial distribution of the amplitude of the electric
near-field for the silicon and gold nanospheres, computed at the FF and SH. We
choose the incident wavelength $\lambda_{inc}=$520\text{\,}\mathrm{nm}$$ for
the gold particle, so as to excite the plasmonic resonance seen in the spectra
of the scattering cross-section corresponding to the FF and SH. In the case of
the silicon nanosphere we pick an off-resonance incident wavelength
$\lambda_{inc}=$1100\text{\,}\mathrm{nm}$$ for our field plots, since Mie
resonances are not present in the frequency range considered here.

Figure 7: Scattering cross-section spectra at the FF (left panels) and SH
(right panels) for three nanospheres made of Au (top panels) and Si (bottom
panels).

There are several ideas illustrated by these field profiles. First, the near-
field distribution at the FF shows the dipolar character of the scattering
process at the wavelength considered. At the SH, on the other hand, the near-
fields in the two cases are of a different nature, namely the field in the
case of th eAu nanosphere is quadrupolar, whereas for the Si nanosphere it is
almost radially symmetric. Second, in both cases the SH near-field is confined
primarily at the surface of the nanosphere, suggesting that the contribution
of the nonlinear surface polarization dominates when compared to that of the
bulk polarization. Moreover, if only the contributions of the bulk SHG are
compared, it appears that it is significantly larger in the case of the Au
nanosphere. This conclusion can be inferred from the spatial distribution of
the electric field inside the particle, as the gradients of its components
determine the bulk polarization. Thus, since the wavelength is much larger
than the size of the sphere, the quasistatic regime is valid and hence one can
view the particle as if it were placed in a uniform electric field.
Consequently, the electric field inside is practically uniform and the bulk
nonlinear polarization vanishes.

As part of the second validation test, we analyze the linear and nonlinear
wave scattering from three nanospheres made of gold and silicon, with centers
located at $O_{1}$($-100\text{\,}\mathrm{nm}$, $-100\text{\,}\mathrm{nm}$,
$100\text{\,}\mathrm{nm}$), $O_{2}$($-150\text{\,}\mathrm{nm}$,
$150\text{\,}\mathrm{nm}$, $100\text{\,}\mathrm{nm}$), and
$O_{3}$($200\text{\,}\mathrm{nm}$, $200\text{\,}\mathrm{nm}$,
$100\text{\,}\mathrm{nm}$), and radii equal to
$R_{1}=$50\text{\,}\mathrm{nm}$$, $R_{2}=$100\text{\,}\mathrm{nm}$$ and
$R_{3}=$150\text{\,}\mathrm{nm}$$. Following the numerical scheme presented in
this work, we meshed each sphere with 2000 triangles and expanded the fields
associated to each sphere into 160 VSWFs. As a result, we have used a total of
480 harmonics for both the FF and SH frequency. With this surface
discretization and field expansions we then computed the surface integrals
defining the $\mathcal{T}$-matrices at the FF and SH and subsequently
evaluated the spectra of the corresponding scattering cross-sections.
Alternatively, we calculated analytically the elements of the
$\mathcal{T}$-matrices, which is possible in the case of spherical particles,
and then the spectra of the scattering cross-sections, namely employed what we
call the Mie-TMM approach [21]. The results obtained using the two methods are
compared in Fig. 7, where an excellent agreement between the predictions of
the method developed here and analytical calculations can be observed both at
the FF and SH.

###  V.3 Wave scattering from prolate spheroids with different aspect ratios

We now consider more challenging examples, which are not amenable to an
analytical treatment. The first such numerical example consists of linear and
nonlinear electromagnetic wave scattering from prolate spheroids made of
either gold or silicon and with different aspect ratio, $\xi=c/a$, where $c$
and $a$ are the polar and equatorial radii of a spheroid, respectively. In
this numerical experiment, we keep the equatorial radius constant,
$a=$50\text{\,}\mathrm{nm}$$, and vary the polar radius from
$c_{min}=$50\text{\,}\mathrm{nm}$$ to $c_{max}=$100\text{\,}\mathrm{nm}$$. The
numerical results for high aspect ratios, $\xi>2$, are not presented because
of the slow convergence, especially in the SH frequency case. This problematic
behavior arises from the fact that the electric fields at the FF evaluated on
the surface of the particles converge very slowly as the number of spherical
harmonics is increased [65]. All the spheroids under consideration are meshed
in approximately 2800 triangles and analyzed with 70 vector spherical
harmonics each, both at the FF and SH frequency. To give a better sense vis-à-
vis the resources required by our method, the PMCHWT integral formulation
employing discretization with RWG-MoM on the same computational mesh would
require just for solving the linear scattering problem around 4200 basis
functions per equivalent current densities (electric and magnetic) totalling
8400 unknowns.

Figure 8: Scattering cross-sections at FF (left panels) and SH (right panels)
for prolate spheroids made of Au (top panels) and Si (bottom panels), computed
for $a=$50\text{\,}\mathrm{nm}$$ and different values of the aspect ratio
$\xi=c/a$.

The spectra of the scattering cross-sections corresponding to these spheroids
are presented in Fig. 8. It can be seen in this figure that, both at the FF
and SH, the scattering cross-section spectra of the spheroids converge to that
of a sphere, depicted with red lines in Fig. 8 and computed analytically, as
the aspect ratio $\xi\to 1$. This is the expected evolution because as $\xi\to
1$ the spheroids become less elongated and approach a spherical shape.
Moreover, the scattering cross-sections increase as the radius $c$ increases,
which is explained by the fact that the geometrical area of the spheroids
increases, too. Importantly, it can be seen in Figs. 8(a) and 8(b) that the
wavelengths of the corresponding plasmonic resonances are redshifted as $\xi$
increases, which is a well-known size-dispersion property of plasmonic
resonance wavelength of nanoparticles.

Figure 9: Spatial distribution of the electric field amplitude evaluated in
the $y=0$ plane passing through the middle of a prolate spheroid with radii
$a=$50\text{\,}\mathrm{nm}$$ and $c=$80\text{\,}\mathrm{nm}$$, at the FF (left
panels) and SH (right panels) for Au (top panels) and Si (bottom panels).

In Fig. 9 we plot the spatial profile of the electric near-field amplitude,
$\lvert\mathbf{E}\rvert$, computed in the $y=0$ plane passing through the
middle of a prolate spheroid with radii $a=$50\text{\,}\mathrm{nm}$$ and
$c=$80\text{\,}\mathrm{nm}$$. The spheroid is centered at the origin of the
coordinate system and the plots are evaluated at the incident wavelengths
$\lambda_{inc}=$580\text{\,}\mathrm{nm}$$ and
$\lambda_{inc}=$1020\text{\,}\mathrm{nm}$$ for the gold and silicon spheroids,
respectively. These wavelengths correspond to the peaks in the SH scattering
cross-section spectra indicating the excitation of plasmonic and Mie
resonances, respectively. The most important conclusion revealed by these
field profiles is that, as in the case of a sphere, most of the SH is
generated at the surface of the nanoparticles.

###  V.4 Wave scattering from a dimer of spheroids

The next particle system we considered consists of two identical prolate
spheroids with radii $a=$50\text{\,}\mathrm{nm}$$ and
$c=$70\text{\,}\mathrm{nm}$$ located at
$O_{1}($0\text{\,}\mathrm{,}$$0\text{\,}\mathrm{,}$$0\text{\,}\mathrm{)}$$ and
$O_{2}($0\text{\,}\mathrm{,}$$0\text{\,}\mathrm{,}$$200\text{\,}\mathrm{nm}$)$.
Thus, the spheroids are separated by a gap with size
$g=$60\text{\,}\mathrm{nm}$$, the center-to-center distance being
$d=$200\text{\,}\mathrm{nm}$$. We considered only the case of spheroids made
of gold, because in this case the presence of the narrow gap between the
spheroids can lead to interesting phenomena, such as the formation in the gap
of hot spots. As in the previous cases, for this configuration we computed the
scattering cross-section spectra and spatial field profiles both at the FF and
SH.

Figure 10: Scattering cross-section at the FF (a) and SH (b) for two prolate
spheroids made of Au.

The scattering cross-section spectra are plotted in Fig. 10, for the
fundamental and second-harmonic frequencies. We calculated these spectra using
the TMM numerical technique presented in this study and compared the results
to the spectra obtained using CST Studio, a commercially available software
suit based on the finite-element method. We meshed each spheroid in 2800
triangles and both in the FF and SH regimes we decomposed the electromagnetic
fields by employing a total of 140 spherical harmonics (i.e., 70 harmonics per
spheroid). This is a substantial reduction of unknown quantities as compared
to an approach based on the FEM, in which case about 67,000 elements were
needed. This figure shows a good agreement between the results obtained by the
two methods, both at the FF and SH and in the entire frequency range of
interest. We stress that this large difference between the number of unknowns
needed to represent the electromagnetic fields in the two methods, makes the
TMM the clear choice when one is interested in the optical response of
clusters containing a large number of particles. In particular, whereas the
use of the FEM becomes problematic even for a few particles, the TMM can be
used to investigate clusters with hundreds or even thousands of particles.

In Figs. 11(a) and 11(b) we present the spatial profile of the electric field
magnitude at the FF and SH, respectively, computed in the $y=0$ plane passing
through the middle of the dimer of spheroids. These field profiles correspond
to an incident wavelength $\lambda_{inc}=$550\text{\,}\mathrm{nm}$$, for which
plasmonic resonances can be observed in the spectra of the scattering cross-
sections at the FF and SH. Since the dimer is not symmetrically located with
respect to the origin of the coordinate system, the spatial profiles of both
the linear and nonlinear fields exhibit a certain asymmetry.

Figure 11: Electric field amplitude at FF (a) and SH (b) computed in the $y=0$
plane passing through the middle of two identical prolate Au spheroids with
radii $a=$50\text{\,}\mathrm{nm}$$ and $c=$70\text{\,}\mathrm{nm}$$ located on
the z-axis.

###  V.5 Wave scattering analysis of three particles with different shapes

As the last example that illustrates powerful features of our numerical
method, we consider the linear and nonlinear wave scattering from a cluster of
three particles with different morphologies. We assume that the particle
cluster comprises of an oblate spheroid with radii
$a_{1}=$50\text{\,}\mathrm{nm}$$ and $c_{1}=$30\text{\,}\mathrm{nm}$$ centered
at
$O_{1}($0\text{\,}\mathrm{,}$-$100\text{\,}\mathrm{nm}$,$100\text{\,}\mathrm{nm}$)$,
a prolate spheroid with radii $a_{2}=$50\text{\,}\mathrm{nm}$$ and
$c_{2}=$70\text{\,}\mathrm{nm}$$ and center at
$O_{2}($0\text{\,}\mathrm{,}$$100\text{\,}\mathrm{nm}$,$100\text{\,}\mathrm{nm}$)$,
and a sphere with radius $r=$50\text{\,}\mathrm{nm}$$ positioned at
$O_{3}($150\text{\,}\mathrm{nm}$,$100\text{\,}\mathrm{nm}$,$100\text{\,}\mathrm{nm}$)$.
The spheroids are oriented such that the z-axis is parallel to their polar
axis. All the particles are assumed to be made of gold, so that the resonant
nature of the scattering process would be more apparent.

The scattering cross-section results used to benchmark our method are computed
with CST Studio using a finite-element mesh containing 42400 tetrahedra. In
our calculations based on the TMM developed here, each nanoparticle is meshed
using about 2500 triangles and analyzed with 96 spherical harmonics, hence
leading to a total of 288 degrees of freedom for the FF and SH computations.
This again represents a substantial reduction of memory requirements as
compared to the resources needed in the FEM computations.

In view of Fig. 12, where the data obtained using the two numerical methods
are summarized, there is a good agreement both in the linear and nonlinear
regimes between the FEM and TMM predictions, in the entire frequency range. As
it transpires from these plots, the differences between the results obtained
with the two methods are more noticeable near the plasmonic resonances. This
is explained by the slow convergence of the near-field, which is the main
physical quantity that determined the nonlinear field. Whereas the accuracy of
the FEM results could have been improved by further increasing the number of
finite-elements used in the computations, we were precluded to do this given
the extent of the computational resources available to us.

Figure 12: Scattering cross sections at FF (a) and SH (b) for three Au
particles with different shapes.

##  VI Conclusion

To summarize, in this manuscript we have extended the T-matrix numerical
formalism to the second-harmonic scattering analysis of arbitrarily shaped
targets made of centrosymmetric optical materials. The mathematical apparatus
developed in this text is general in purpose and can be in principle applied
to the SH treatment of targets with arbitrary shapes. However, in traditional
TMM fashion, we restrict ourselves to the analysis of particles with axial
symmetry, where TMM performs the best due to the underlying spherical
expansion functions employed. The numerical technique presented in this text
is developed from the integral representations of the fundamental and second-
harmonic fields and the imposition of the null-field conditions via EBCM. The
integral operators are then discretized with vector spherical harmonic basis
functions and the matrix equation is established relating the unknown field
and the known source expansion coefficients. The SH polarization sources,
induced in the bulk and on the surface of the target under analysis, are
conveniently introduced into TMM formalism via nonlinear tangential boundary
conditions thus allowing the accurate physical description of the second-
harmonic signal.

The method presented in this study is particularly useful for the analysis of
different kinds of ellipsoids, where certain symmetries can be exploited
through the spherical expansion of the fields and sources generating
substantially lower memory imprint when compared to the standard methods based
on finite element discretizations and subsectional bases. Furthermore, this
method is convenient for the analysis of particle clusters since the
individual scatterer T-matrices can be precomputed and connected with
translation-addition operators to form a scattering matrix of a whole particle
assembly, a process which is easily parallelized and independent of the
underlying excitation field. Thus, we expect this method will be useful in the
design process of finite parts of metamaterials or photonic crystals where the
edge effects can not be disregarded.

## Acknowledgement

This work was funded under the embedded CSE programme of the ARCHER2 UK
National Supercomputing Service and used the ARCHER2 UK National
Supercomputing Service (http://www.archer2.ac.uk). We acknowledge the use of
the UCL High Performance Computing Facility and associated support services in
the completion of this work.

##  Appendix A Formulae for $\mathbf{Q}$-matrix elements

In this Appendix we provide the mathematical expressions for the elements of
the matrix $\mathbf{Q}^{\rho(1,3)}$, $\rho\in\\{\omega,\Omega\\}$, at FF and
SH frequency, as well as the components of the excitation vector
$\mathbf{g}^{\Omega(3)}$ at the SH. The corresponding elements of the matrix
$\mathbf{Q}^{\rho(1,1)}$ and components of the vector $\mathbf{g}^{\Omega(1)}$
are obtained by replacing in these expressions the outgoing vector spherical
harmonics, $(3)$, with the regular ones, $(1)$. Furthermore, we assume that a
boundary mesh of the particle at hand is available, and defined as the union
of $N_{t}$ non-overlapping flat triangular facets
$\tilde{S}=\bigcup\limits_{t=1}^{N_{t}}S_{t}$.

The matrix $\mathbf{Q}^{\rho(1,3)}$ can be written as a combination of four
submatrices $\mathbf{Q}^{\rho,11}$, $\mathbf{Q}^{\rho,12}$,
$\mathbf{Q}^{\rho,21}$, and $\mathbf{Q}^{\rho,22}$, i.e.:

$\mathbf{Q}^{\rho(1,3)}=\begin{pmatrix}\mathbf{Q}^{\rho,11}&\mathbf{Q}^{\rho,12}\\\
\mathbf{Q}^{\rho,21}&\mathbf{Q}^{\rho,22}\end{pmatrix},$ (A.1)

with the elements given by:

$\displaystyle Q^{\rho,11}_{\mu\nu}$
$\displaystyle=i(k_{e}^{\rho})^{2}\oiint_{\tilde{S}}\hat{\mathbf{n}}\cdot\left[\mathbf{M}_{\nu}^{(1)}(k_{i}^{\rho}\mathbf{r})\times\mathbf{N}_{\bar{\mu}}^{(3)}(k_{e}^{\rho}\mathbf{r})\right]dS$
$\displaystyle+$ $\displaystyle
ik_{e}^{\rho}k_{i}^{\rho}\oiint_{\tilde{S}}\hat{\mathbf{n}}\cdot\left[\mathbf{N}_{\nu}^{(1)}(k_{i}^{\rho}\mathbf{r})\times\mathbf{M}_{\bar{\mu}}^{(3)}(k_{e}^{\rho}\mathbf{r})\right]dS,$
(A.2a) $\displaystyle Q^{\rho,12}_{\mu\nu}$
$\displaystyle=i(k_{e}^{\rho})^{2}\oiint_{\tilde{S}}\hat{\mathbf{n}}\cdot\left[\mathbf{N}_{\nu}^{(1)}(k_{i}^{\rho}\mathbf{r})\times\mathbf{N}_{\bar{\mu}}^{(3)}(k_{e}^{\rho}\mathbf{r})\right]dS$
$\displaystyle+$ $\displaystyle
ik_{e}^{\rho}k_{i}^{\rho}\oiint_{\tilde{S}}\hat{\mathbf{n}}\cdot\left[\mathbf{M}_{\nu}^{(1)}(k_{i}^{\rho}\mathbf{r})\times\mathbf{M}_{\bar{\mu}}^{(3)}(k_{e}^{\rho}\mathbf{r})\right]dS,$
(A.2b) $\displaystyle Q^{\rho,21}_{\mu\nu}$
$\displaystyle=i(k_{e}^{\rho})^{2}\oiint_{\tilde{S}}\hat{\mathbf{n}}\cdot\left[\mathbf{M}_{\nu}^{(1)}(k_{i}^{\rho}\mathbf{r})\times\mathbf{M}_{\bar{\mu}}^{(3)}(k_{e}^{\rho}\mathbf{r})\right]dS$
$\displaystyle+$ $\displaystyle
ik_{e}^{\rho}k_{i}^{\rho}\oiint_{\tilde{S}}\hat{\mathbf{n}}\cdot\left[\mathbf{N}_{\nu}^{(1)}(k_{i}^{\rho}\mathbf{r})\times\mathbf{N}_{\bar{\mu}}^{(3)}(k_{e}^{\rho}\mathbf{r})\right]dS,$
(A.2c) $\displaystyle Q^{\rho,22}_{\mu\nu}$
$\displaystyle=i(k_{e}^{\rho})^{2}\oiint_{\tilde{S}}\hat{\mathbf{n}}\cdot\left[\mathbf{N}_{\nu}^{(1)}(k_{i}^{\rho}\mathbf{r})\times\mathbf{M}_{\bar{\mu}}^{(3)}(k_{e}^{\rho}\mathbf{r})\right]dS$
$\displaystyle+$ $\displaystyle
ik_{e}^{\rho}k_{i}^{\rho}\oiint_{\tilde{S}}\hat{\mathbf{n}}\cdot\left[\mathbf{M}_{\nu}^{(1)}(k_{i}^{\rho}\mathbf{r})\times\mathbf{N}_{\bar{\mu}}^{(3)}(k_{e}^{\rho}\mathbf{r})\right]dS,$
(A.2d)

where $\nu$ and $\bar{\mu}$ are the combined indices $(l_{1},m_{1})$ and
$(l_{2},-m_{2})$, respectively.

The components of the excitation vector $\mathbf{g}^{\Omega(3)}$ are defined
by surface integrals of products between tangential source currents and VSWFs.
More specifically, the vector $\mathbf{g}^{\Omega(3)}$ can be decomposed into
a stack of vectors $\mathbf{g}_{\mu}^{\Omega}$ defined as:

$\mathbf{g}_{\mu}^{\Omega}=\oiint_{\tilde{S}}\bm{\mathcal{C}}^{\Omega}(\mathbf{r})\cdot\begin{pmatrix}\mathbf{M}_{\bar{\mu}}^{(3)}(k_{e}^{\Omega}\mathbf{r})\\\
\mathbf{N}_{\bar{\mu}}^{(3)}(k_{e}^{\Omega}\mathbf{r})\end{pmatrix}dS,$ (A.3)

where the $2\times 2$ matrix $\bm{\mathcal{C}}^{\Omega}(\mathbf{r})$ is
defined in terms of the electric and magnetic currents,
$\bm{\mathcal{J}}^{\Omega}(\mathbf{r})$ and
$\bm{\mathcal{M}}^{\Omega}(\mathbf{r})$, respectively, and the internal field
$\mathbf{E}_{i,p}^{\Omega}(\mathbf{r})$:

$\displaystyle\bm{\mathcal{C}}^{\Omega}_{11}=\bm{\mathcal{C}}^{\Omega}_{22}=\mu_{0}\Omega
k_{e}^{\Omega}\bm{\mathcal{J}}^{\Omega},$ (A.4a)
$\displaystyle\bm{\mathcal{C}}^{\Omega}_{12}=\bm{\mathcal{C}}^{\Omega}_{21}=i(k_{e}^{\Omega})^{2}\left(\bm{\mathcal{M}}^{\Omega}-\hat{\mathbf{n}}\times\mathbf{E}_{i,p}^{\Omega}\right).$
(A.4b)

The integrals in Eq. (A.3) are of three different types, namely

$\displaystyle I_{1}=\mu_{0}\Omega
k_{e}^{\Omega}\oiint_{\tilde{S}}\bm{\mathcal{J}}^{\Omega}(\mathbf{r})\cdot\mathbf{M}_{\bar{\mu}}^{(3)}(k_{e}^{\Omega}\mathbf{r})\,dS,$
(A.5a) $\displaystyle I_{1}^{\prime}=\mu_{0}\Omega
k_{e}^{\Omega}\oiint_{\tilde{S}}\bm{\mathcal{J}}^{\Omega}(\mathbf{r})\cdot\mathbf{N}_{\bar{\mu}}^{(3)}(k_{e}^{\Omega}\mathbf{r})\,dS,$
(A.5b) $\displaystyle
I_{2}=i(k_{e}^{\Omega})^{2}\oiint_{\tilde{S}}\bm{\mathcal{M}}^{\Omega}(\mathbf{r})\cdot\mathbf{N}_{\bar{\mu}}^{(3)}(k_{e}^{\Omega}\mathbf{r})\,dS,$
(A.5c) $\displaystyle
I_{2}^{\prime}=i(k_{e}^{\Omega})^{2}\oiint_{\tilde{S}}\bm{\mathcal{M}}^{\Omega}(\mathbf{r})\cdot\mathbf{M}_{\bar{\mu}}^{(3)}(k_{e}^{\Omega}\mathbf{r})\,dS,$
(A.5d) $\displaystyle
I_{3}=i(k_{e}^{\Omega})^{2}\oiint_{\tilde{S}}\left[\hat{\mathbf{n}}\times\mathbf{E}_{i,p}^{\Omega}(\mathbf{r})\right]\cdot\mathbf{N}_{\bar{\mu}}^{(3)}(k_{e}^{\Omega}\mathbf{r})\,dS,$
(A.5e) $\displaystyle
I_{3}^{\prime}=i(k_{e}^{\Omega})^{2}\oiint_{\tilde{S}}\left[\hat{\mathbf{n}}\times\mathbf{E}_{i,p}^{\Omega}(\mathbf{r})\right]\cdot\mathbf{M}_{\bar{\mu}}^{(3)}(k_{e}^{\Omega}\mathbf{r})\,dS.$
(A.5f)

First, let us consider the evaluation of the integral $I_{1}$. With the help
of Eqs. (20a) and (2), we can write:

$\displaystyle I_{1}=$ $\displaystyle
2i\chi_{s,\parallel\perp\parallel}^{(2)}k_{e}^{\Omega}(k_{0}^{\Omega})^{2}\oiint_{\tilde{S}}E_{n}^{\omega}(\mathbf{r})$
(A.6)
$\displaystyle\cdot\left\\{\left[\hat{\mathbf{n}}\times\hat{\mathbf{n}}\times\mathbf{E}_{t}^{\omega}(\mathbf{r})\right]\cdot\mathbf{M}_{\bar{\mu}}^{(3)}(k_{e}^{\Omega}\mathbf{r})\right\\}dS,$

where the fields are evaluated on the inner side of the surface of the
particle and $k_{0}^{\Omega}$ is the wave number in vacuum at frequency
$\Omega$. Written in this form, the integral $I_{1}$ is suitable for the
numerical evaluation.

Now we turn our attention to the integral $I_{2}$ and use Eq. (20b) defining
the nonlinear magnetic surface current in terms of the normal component of the
surface polarization density source,
$P_{n,s}^{\Omega}(\mathbf{r})=\hat{\mathbf{n}}\cdot\bm{\mathbf{P}}_{s}^{\Omega}(\mathbf{r})$:

$I_{2}=\frac{i(k_{e}^{\Omega})^{2}}{\epsilon_{0}}\oiint_{\tilde{S}}\left[\hat{\mathbf{n}}\times\nabla_{S}P_{s,n}^{\Omega}(\mathbf{r})\right]\cdot\mathbf{N}_{\bar{\mu}}^{(3)}(k_{e}^{\Omega}\mathbf{r})\,dS.$
(A.7)

We break this integral over the triangular patches as:

$\displaystyle
I_{2}=\frac{i(k_{e}^{\Omega})^{2}}{\epsilon_{0}}\sum\limits_{p=1}^{N_{t}}\left\\{\oint_{C_{p}}P_{s,n}^{\Omega}(\mathbf{r})\mathbf{N}_{\bar{\mu}}^{(3)}(k_{e}^{\Omega}\mathbf{r})\cdot
d\mathbf{l}\right.$
$\displaystyle\left.-\iint_{S_{p}}P_{s,n}^{\Omega}(\mathbf{r})\left[\nabla_{S}\times\mathbf{N}_{\bar{\mu}}^{(3)}(k_{e}^{\Omega}\mathbf{r})\right]\cdot\hat{\mathbf{n}}\,dS\right\\}.$
(A.8)

In this equation, $C_{p}=\partial S_{p}$ is the boundary of the $p$th triangle
in the mesh and the direction of integration along the path $C_{p}$ is chosen
in such a way that the right-hand-screw rule gives the direction of the
outwardly oriented normal, $\hat{\mathbf{n}}$. To derive Eq. (A) from Eq.
(A.7) one first uses the vector calculus relation
$\left[\hat{\mathbf{n}}\times\nabla_{S}\psi(\mathbf{r})\right]\cdot\mathbf{A}(\mathbf{r})=\left\\{\nabla_{S}\times[\psi(\mathbf{r})\mathbf{A}(\mathbf{r})]-\psi(\mathbf{r})\nabla_{S}\times\mathbf{A}(\mathbf{r})\right\\}\cdot\hat{\mathbf{n}}$,
where $\psi(\mathbf{r})$ and $\mathbf{A}(\mathbf{r})$ are arbitrary smooth
scalar and vector functions, respectively, and then employs the Stokes
theorem.

The numerical evaluation of the integral $I_{2}$ is greatly simplified by the
realization that the contour integral on the right hand side of Eq. (A) is
equal to zero when it is evaluated on a closed mesh without cracks and slits.
Under these circumstances, the final expression for the integral $I_{2}$ is:

$I_{2}=-i(k_{e}^{\Omega})^{3}\chi_{s,\perp\perp\perp}^{(2)}\oiint_{\tilde{S}}[E_{n}^{\omega}(\mathbf{r})]^{2}\mathbf{M}_{\bar{\mu}}^{(3)}(k_{e}^{\Omega}\mathbf{r})\cdot\hat{\mathbf{n}}\,dS.$
(A.9)

In the derivation of this expression we took advantage of the curl identities
in Eqs. (B.2) satisfied by the VSWFs, as well as the fact that the surface
tensor component $\chi_{s,\perp\parallel\parallel}^{(2)}$ is vanishingly small
and thus its contribution to the SHG can be disregarded. As such, the normal
component of the surface polarization density is:
$P_{s,n}^{\Omega}(\mathbf{r})=\epsilon_{0}\chi_{s,\perp\perp\perp}^{(2)}[E_{n}^{\omega}(\mathbf{r}_{s})]^{2}$.

The last integral, $I_{3}$, is written as:

$I_{3}=\frac{i(k_{e}^{\Omega})^{2}\epsilon_{0}\gamma}{\epsilon_{i}^{\Omega}}\oiint_{\tilde{S}}\left\\{\hat{\mathbf{n}}\times\nabla\left[\mathbf{E}_{i}^{\omega}(\mathbf{r})\cdot\mathbf{E}_{i}^{\omega}(\mathbf{r})\right]\right\\}\cdot\mathbf{N}_{\bar{\mu}}^{(3)}(k_{e}^{\Omega}\mathbf{r})\,dS,$
(A.10)

where we have used the particular electric field solution given by Eq. (21b)
of the inhomogeneous system of Maxwell equations at the SH. This integral can
be calculated using a similar method to that used to calculate $I_{2}$, the
end result being:

$I_{3}=-\frac{i(k_{e}^{\Omega})^{3}\epsilon_{0}\gamma}{\epsilon_{i}^{\Omega}}\oiint_{\tilde{S}}\left[\mathbf{E}_{i}^{\omega}(\mathbf{r})\cdot\mathbf{E}_{i}^{\omega}(\mathbf{r})\right]\left[\mathbf{M}_{\bar{\mu}}^{(3)}(k_{e}^{\Omega}\mathbf{r})\cdot\hat{\mathbf{n}}\right]\,dS.$
(A.11)

Finally, the primed integrals, $I_{1}^{\prime}$, $I_{2}^{\prime}$, and
$I_{3}^{\prime}$, can be obtained from $I_{1}$, $I_{2}$, and $I_{3}$,
respectively, by interchanging the spherical harmonics
$\mathbf{M}_{\bar{\mu}}^{(3)}(k_{e}^{\Omega}\mathbf{r})\leftrightarrow\mathbf{N}_{\bar{\mu}}^{(3)}(k_{e}^{\Omega}\mathbf{r})$
in Eqs. (A.6), (A.9), and (A.11). All of the integrals can be computed
numerically using Gauss-Legendre quadrature in six points.

##  Appendix B Definition and formulae for vector spherical harmonics

In this Appendix we introduce the definition and main properties of vector
spherical harmonics (VSHs) used in our study to represent the electromagnetic
fields, as well as to expand the plane waves.

### Definition of vector spherical wave functions

Given the unit vector $\hat{\mathbf{r}}=\mathbf{r}/r$, where $\mathbf{r}$ is
the position vector of a point, $P$, we define three sets of mutually
orthogonal VSHs: longitudinal VSHs, $\mathbf{Y}_{lm}$, which are oriented
along $\hat{\mathbf{r}}$, and two sets of transverse VSHs, $\mathbf{X}_{lm}$
and $\mathbf{Z}_{lm}$, oriented perpendicular to $\hat{\mathbf{r}}$ [36, 56].
The triplet $\\{\mathbf{X}_{lm},\mathbf{Y}_{lm},\mathbf{Z}_{lm}\\}$ of VSHs
form a right-handed system of vectors and are defined as:

$\displaystyle\mathbf{Y}_{lm}(\theta,\phi)=Y_{lm}(\theta,\phi)\hat{\mathbf{r}},$
(B.1a)
$\displaystyle\mathbf{X}_{lm}(\theta,\phi)=\frac{1}{\sqrt{l(l+1)}}\hat{\mathbf{r}}\times\nabla_{\hat{\mathbf{r}}}Y_{lm}(\theta,\phi),$
(B.1b)
$\displaystyle\mathbf{Z}_{lm}(\theta,\phi)=\frac{1}{\sqrt{l(l+1)}}\nabla_{\hat{\mathbf{r}}}Y_{lm}(\theta,\phi).$
(B.1c)

Furthermore, we define the regular VSWFs $\mathbf{M}_{lm}^{(1)}$ and
$\mathbf{N}_{lm}^{(1)}$ as follows:

$\displaystyle\mathbf{M}_{lm}^{(1)}(kr,\theta,\phi)=-j_{l}(kr)\mathbf{X}_{lm},$
(B.2a)
$\displaystyle\mathbf{N}_{lm}^{(1)}(kr,\theta,\phi)=\frac{1}{k}\nabla\times\mathbf{M}_{lm}^{(1)},$
(B.2b)
$\displaystyle\mathbf{M}_{lm}^{(1)}(kr,\theta,\phi)=\frac{1}{k}\nabla\times\mathbf{N}_{lm}^{(1)},$
(B.2c)

where $k$ is the wave number calculated at the particular frequency in the
medium of choice and $j_{l}(kr)$ are spherical Bessel functions of the first
kind and order $l$. The radiative (outgoing) VSWFs $\mathbf{M}_{lm}^{(3)}$ and
$\mathbf{N}_{lm}^{(3)}$ are defined by substituting the spherical Bessel
functions in (B.2) with (outgoing) spherical Hankel functions,
$h_{l}^{(1)}(kr)$, of the first kind and order $l$.

Scalar spherical harmonics can be defined in terms of associated Legendre
functions $P_{l}^{m}$ as [56]:

$Y_{lm}(\theta,\phi)=\gamma_{lm}\sqrt{l(l+1)}P_{l}^{m}(\cos\theta)e^{im\phi},$
(B.3)

where $\gamma_{lm}$ is the normalization coefficient given by:

$\gamma_{lm}=\sqrt{\frac{(2l+1)(l-m)!}{4\pi l(l+1)(l+m)!}}.$ (B.4)

In these definitions the scalar spherical harmonic functions
$Y_{lm}(\hat{\mathbf{r}})$ and the three kinds of VSHs are normalized to unity
over a sphere:

$\displaystyle\iint_{\Omega_{\hat{\mathbf{r}}}}Y_{lm}^{*}(\hat{\mathbf{r}})Y_{l^{\prime}m^{\prime}}(\hat{\mathbf{r}})d\Omega_{\hat{\mathbf{r}}}=\delta_{ll^{\prime}}\delta_{mm^{\prime}},$
(B.5a)
$\displaystyle\iint_{\Omega_{\hat{\mathbf{r}}}}\mathbf{A}_{lm}^{*}(\hat{\mathbf{r}})\cdot\mathbf{B}_{l^{\prime}m^{\prime}}(\hat{\mathbf{r}})d\Omega_{\hat{\mathbf{r}}}=\delta_{ll^{\prime}}\delta_{mm^{\prime}}\delta_{AB},$
(B.5b)

where $\mathbf{A}$ and $\mathbf{B}$ can be any of the functions $\mathbf{X}$,
$\mathbf{Y}$, and $\mathbf{Z}$.

Finally, the normalized VSHs can be defined in terms of associated Legendre
functions as:

$\displaystyle\mathbf{Y}_{lm}(\theta,\phi)=\,$
$\displaystyle\gamma_{lm}\sqrt{l(l+1)}P_{l}^{m}(\cos\theta)e^{im\phi}\hat{\mathbf{r}},$
(B.6a) $\displaystyle\mathbf{X}_{lm}(\theta,\phi)=\,$
$\displaystyle\gamma_{lm}\left[-\frac{im}{\sin\theta}P_{l}^{m}(\cos\theta)e^{im\phi}\hat{\bm{\theta}}\right.$
$\displaystyle+\left.\frac{d}{d\theta}P_{l}^{m}(\cos\theta)e^{im\phi}\hat{\bm{\phi}}\right],$
(B.6b) $\displaystyle\mathbf{Z}_{lm}(\theta,\phi)=\,$
$\displaystyle\gamma_{lm}\left[\frac{d}{d\theta}P_{l}^{m}(\cos\theta)e^{im\phi}\hat{\bm{\theta}}\right.$
$\displaystyle+\left.\frac{im}{\sin\theta}P_{l}^{m}(\cos\theta)e^{im\phi}\hat{\bm{\phi}}\right].$
(B.6c)

### Expansion of plane waves in series of vector spherical wave functions

The series expansion of a plane wave in terms of VSWFs is given in (9a) and
(10a), with the expansion coefficients $q_{lm}^{\omega}$ and $p_{lm}^{\omega}$
calculated as follows [56]:

$\displaystyle q_{lm}^{\omega}$ $\displaystyle=-4\pi
i^{l}\mathbf{X}_{lm}^{*}(\hat{\mathbf{k}}_{inc})\cdot\hat{\mathbf{e}}_{inc},$
(B.7a) $\displaystyle p_{lm}^{\omega}$ $\displaystyle=-4\pi
i^{l+1}\mathbf{Z}_{lm}^{*}(\hat{\mathbf{k}}_{inc})\cdot\hat{\mathbf{e}}_{inc},$
(B.7b)

where the symbol “$*$” denotes complex conjugation, $\hat{\mathbf{k}}_{inc}$
is the unit vector in the direction of the wave propagation, and
$\hat{\mathbf{e}}_{inc}$ is the polarization direction of the plane wave.

## References

  * Boyd [1992] R. W. Boyd, Nonlinear Optics (Academic Press, 1992).
  * Ray [2003] P. C. Ray, Size and Shape Dependent Second Order Nonlinear Optical Properties of Nanomaterials and Their Application in Biological and Chemical Sensing, Chem. Rev. 110, 5332 (2010).
  * Hatab [2010] N. A. Hatab, C. H. Hsueh, A. L. Gaddis, S. T. Retterer, J. H. Li, G. Eres, Z. Zhang, and B. Gu, Free-Standing Optical Gold Bowtie Nanoantenna with Variable Gap Size for Enhanced Raman Spectroscopy, Nano Lett. 10, 4952 (2010).
  * Butet [2015] J. Butet, P. F. Brevet, and O. J. F. Martin, Optical second harmonic generation in plasmonic nanostructures: from fundamental principles to advanced applications. ACS Nano 9, 10545 (2015).
  * Panoiu [2018] N. C. Panoiu, W. E. I. Sha, D. Y. Lei, and G. C. Li, Nonlinear optics in plasmonic nanostructures, J. Opt. 20, 083001 (2018).
  * Homola [2008] J. Homola, Surface plasmon resonance sensors for detection of chemical and biological species, Chem. Rev. 108, 462 (2008).
  * Pendry [1999] J. B. Pendry, A. J. Holden, D. J. Robbins, and W. J. Stewart, Magnetism from conductors and enhanced nonlinear phenomena, IEEE Trans. Microw. Theory Tech. 47, 2075 (1999).
  * Liu [2008] N. Liu, H. Guo, L. Fu, S. Kaiser, H. Schweizer, and H. Giessen, Three-dimensional photonic metamaterials at optical frequencies, Nat. Mater. 7, 31 (2008).
  * Ziolkowski [2004] R. W. Ziolkowski, Propagation in and scattering from a matched metamaterial having a zero index of refraction, Phys. Rev. E 70, 046608 (2004).
  * Canfield [2006] B. K. Canfield, S. Kujala, K. Laiho, K. Jefimovs, J. Turunen, and M. Kauranen, Chirality arising from small defects in gold nanoparticle arrays, Opt. Express 14, 950 (2006).
  * Valev [2011] V. K. Valev, X. Zheng, C. G. Biris, A. V. Silhanek, V. Volskiy, B. De Clercq, O. A. Aktsipetrov, M. Ameloot, N. C. Panoiu, G. A. E. Vandenbosch, and V. V. Moshchalkov, The origin of second harmonic generation hotspots in chiral optical metamaterials, Opt. Mater. Express 1, 36 (2011).
  * Malkiel [2018] I. Malkiel, M. Mrejen, A. Nagler, U. Arieli, L. Wolf, and H. Suchowski, Plasmonic nanostructure design and characterization via Deep Learning, Light: Science & Applications 7, 60 (2018).
  * Molesky [2018] S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vuckovic, and A. W. Rodriguez, Inverse design in nanophotonics, Nat. Photon. 12, 659 (2018).
  * Kudyshev [2020] Z. A. Kudyshev, A. V. Kildishev, V. M. Shalaev, and A. Boltasseva, Machine-learning-assisted metasurface design for high efficiency thermal emitter optimization, Appl. Phys. Rev. 7, 021407 (2020).
  * Shen [2003] Y. R. Shen, The Principles of Nonlinear Optics (Wiley, 2003).
  * Biris [2010] C. G. Biris and N. C. Panoiu, Second harmonic generation in metamaterials based on homogeneous centrosymmetric nanowires, Phys. Rev. B 81, 195102 (2010).
  * Pavlyukh [2004] Y. Pavlyukh and W. Hubner, Nonlinear Mie scattering from spherical particles, Phys. Rev. B 70, 245434 (2004).
  * Beer1 [2007] A. G. F. de Beer and S. Roke, Nonlinear Mie theory for second-harmonic and sum-frequency scattering, Phys. Rev. B 79, 155420 (2007).
  * Forestiere1 [2014] A. Capretti, C. Forestiere, L. Dal Negro, and G. Miano, Full-wave analytical solution of second-harmonic generation in metal nanospheres, Plasmonics 9, 151 (2014).
  * Xu [2012] J. Xu and X. Zhang, Second harmonic generation in three-dimensional structures based on homogeneous centrosymmetric metallic spheres, Opt. Express 20, 1668 (2012).
  * Sekulic [2021] I. Sekulic, J. W. You, and N. C. Panoiu, T-matrix method for calculation of second-harmonic generation in clusters of spherical particles, J. Quant. Spectrosc. Radiat. Transfer 268, 107643 (2021).
  * Silvester [1996] P. P. Silvester and R. L. Ferrari, Finite Elements for Electrical Engineers, (Cambridge University Press, 1996).
  * Jin [2002] J. M. Jin, The Finite Element Method in Electromagnetics, (John Wiley & Sons, 2002).
  * Yee [1966] K. Yee, Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media, IEEE Trans. Antennas Propag. AP-14, 302 (1966).
  * Taflove [2005] A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, (Artech House, 2005).
  * Mayergoyz [2005] I. D. Mayergoyz, D. R. Fredkin, and Z. Zhang, Electrostatic (plasmon) resonances in nanoparticles, Phys. Rev. B 72, 155412 (2005).
  * Hohenester1 [2005] U. Hohenester and J. Krenn, Surface plasmon resonances of single and coupled metallic nanoparticles: A boundary integral method approach, Phys. Rev. B 72, 195429 (2005).
  * Poggio [1973] A. J. Poggio and E. K. Miller, Integral equation solutions of three-dimensional scattering problems, in Computer Techniques for Electromagnetics, edited by R. Mittra (Pergamon Press, Oxford, 1973).
  * Wu [1977] T. K. Wu and L. L. Tsai, Scattering from arbitrarily-shaped lossy dielectric bodies of revolution, Radio Sci. 12, 709 (1977).
  * Makitalo [2011] J. Makitalo, S. Suuriniemi, and M. Kauranen, Boundary element method for surface nonlinear optics of nanoparticles, Opt. Express 19, 23386 (2011).
  * Forestiere2 [2013] C. Forestiere, A. Capretti, and G. Miano, Surface integral method for second harmonic generation in metal nanoparticles including both local-surface and nonlocal-bulk sources, J. Opt. Soc. Am. B 30, 2355 (2013).
  * Harrington [1968] R. F. Harrington, Field Computation by Moment Method, (Macmillan 1968).
  * Rao [1982] S. M. Rao, D. R. Wilton, and A. W. Glisson, Electromagnetic scattering by surfaces of arbitrary shape, IEEE Trans. Antennas Propag. 30, 409 (1982).
  * Song [1997] J. Song, C. C. Lu, and W. C. Chew, Multilevel fast multipole algorithm for electromagnetic scattering by large complex objects, IEEE Trans. Antennas Propag. 45, 1488 (1997).
  * Bebendorf [2000] M. Bebendorf, Approximation of boundary element matrices, Numer. Math. 86, 565 (2000).
  * Varshalovich [1988] D. A. Varshalovich, A. N. Moskalev, and V. K. Khersonskii, Quantum Theory of Angular Momentum, (World Scientific Publishing 1988).
  * Jackson [1999] J. D. Jackson, Classical Electrodynamics, (John Wiley & Sons, 1999).
  * Waterman1 [1965] P. C. Waterman, Matrix formulation of electromagnetic scattering, Proc. IEEE 53, 805 (1965).
  * Waterman2 [1969] P. C. Waterman, Scattering by dielectric obstacles, Alta Freq. 38, 348 (1969).
  * Barber [1975] P. Barber and C. Yeh, Scattering of electromagnetic waves by arbitrarily shaped dielectric bodies, Appl. Opt. 14, 2864 (1975).
  * Mishchenko1 [1996] M. I. Mishchenko, L. D. Travis, and D. W. Mackowski, T-matrix computations of light scattering by nonspherical particles: a review, J. Quant. Spectrosc. Radiat. Transfer 55, 535 (1996).
  * Mishchenko1 [2007] M. I. Mishchenko, G. Videen, V. A. Babenko, N. G. Khlebtsov, and T. Wriedt, Comprehensive T-matrix reference database: A 2004-06 update, J. Quant. Spectrosc. Radiat. Transfer 106, 304 (2007).
  * Cruzan [1962] O. R. Cruzan, Translational addition theorems for spherical vector wave functions, Quart. Appl. Math. 20, 33 (1962).
  * Mackowski [1996] D. W. Mackowski and M. I. Mishchenko, Calculation of the T-matrix and the scattering matrix for ensembles of spheres, J. Opt. Soc. Am. A 13, 2266 (1996).
  * Stout1 [2002] B. Stout, J. C. Auger, and J. Lafait, A transfer matrix approach to local field calculations in multiple-scattering problems, J. Mod. Opt. 49, 2129 (2002).
  * Timbrell [2018] D. Timbrell, J. W. You, Y. S. Kivshar, and N. C Panoiu, A comparative analysis of surface and bulk contributions to second-harmonic generation in centrosymmetric nanoparticles, Sci. Rep. 8, 3586 (2018).
  * Dadap [2004] J. I. Dadap, J. Shan, and T. F. Heinz, Theory of optical second-harmonic generation from a sphere of centrosymmetric material: small-particle limit, J. Opt. Soc. Am. B 21, 1328 (2004).
  * Bloembergen [1968] N. Bloembergen, R. K. Chang, S. S. Jha, and C. H. Lee, Optical second-harmonic generation in reflection from media with inversion symmetry, Phys. Rev. 174, 813 (1968).
  * Sipe1 [1980] J. E. Sipe, V. C. Y. So, M. Fukui, and G. I. Stegeman, Analysis of second-harmonic generation at metal surfaces, Phys. Rev. B 21, 4389 (1980).
  * Tsang [1985] L. Tsang, J. A. Kong, and R. T. Shin, Theory of Microwave Remote Sensing, (Wiley, 1985).
  * Stratton [1939] J. A. Stratton and L. J. Chu, Diffraction theory of electromagnetic waves, Phys. Rev. 56, 99 (1939).
  * Chew [2009] W. C. Chew, M. S. Tong, and B. Hu, Integral Equation Methods for Electromagnetic and Elastic Waves, (Morgan & Claypool, 2009).
  * Morse [1953] P. Morse and H. Feshbach, Methods of Theoretical Physics, (McGraw-Hill, 1953).
  * Sekulic [2019] I. Sekulic, D. C. Tzarouchis, P. Yla-Oijala, E. Ubeda, and J. M. Rius, Enhanced discretization of surface integral equations for resonant scattering analysis of sharp-edged plasmonic nanoparticles, Phys. Rev. B 99, 165417 (2019).
  * Heinz [1991] T. F. Heinz, in Nonlinear Surface Electromagnetic Phenomena, edited by H. E. Ponath and G. I. Stegeman (Elsevier, Amsterdam, 1991), p. 353.
  * Stout2 [2001] B. Stout, J. C. Auger, and J. Lafait, Individual and aggregate scattering matrices and cross-sections: conservation laws and reciprocity, J. Mod. Opt. 48, 2105 (2001).
  * Mie1 [1908] G. Mie, Beitrage zur optik trüber medien, speziell kolloidaler metallösungen, Ann. Phys. (Berlin) 330, 377 (1908).
  * [58] CST Studio®, www.cst.com.
  * Schinke [2015] C. Schinke, P. C. Peest, J. Schmidt, R. Brendel, K. Bothe, M. R. Vogt, I. Kroger, S. Winter, A. Schirmacher, S. Lim, H. T. Nguyen, and D. MacDonald, Uncertainty analysis for the coefficient of band-to-band absorption of crystalline silicon, AIP Advances 5, 067168 (2015).
  * Johnson [1972] P. B. Johnson and R. W. Christy, Optical constants of the noble metals, Phys. Rev. B 6, 4370 (1972).
  * Rudnick [1971] J. Rudnick and E. A. Stern, Second-harmonic radiation from metal surfaces, Phys. Rev. B 4, 4274 (1971).
  * Falasconi [2001] M. Falasconi, L. C. Andreani, A. M. Malvezzi, M. Patrini, V. Mulloni, and L. Pavesic, Bulk and surface contributions to second-order susceptibility in crystalline and porous silicon by second-harmonic generation, Surf. Sci. 481, 105 (2011).
  * Corvi [1986] M. Corvi and W. Schaich, Hydrodynamic-model calculation of second-harmonic generation at a metal surface, Phys. Rev. B 33, 3688 (1986).
  * Guyot [1986] P. Guyot-Sionnest, W. Chen, and Y. R. Shen, General considerations on optical second-harmonic generation from surfaces and interfaces, Phys. Rev. B 33, 8254 (1986).
  * Somerville [2015] W. R. C. Somerville, B. Auguie, and E. C. Le Ru, Accurate and convergent T-matrix calculations of light scattering by spheroids, J. Quant. Spectrosc. Radiat. Transfer 160, 29 (2015).