# Inverse design of focused vector beams
for mode excitation in optical nanoantennas

Xiaorun Zang<EMAIL_ADDRESS>Currently at Department of Applied
Physics, Aalto University, P.O Box 13500, FI-00076 Aalto, Finland  Institute
of Photonics, University of Eastern Finland, P.O. Box 111, FI-80101 Joensuu,
Finland  Ari T. Friberg  Institute of Photonics, University of Eastern
Finland, P.O. Box 111, FI-80101 Joensuu, Finland  Tero Setälä  Institute of
Photonics, University of Eastern Finland, P.O. Box 111, FI-80101 Joensuu,
Finland  Jari Turunen  Institute of Photonics, University of Eastern Finland,
P.O. Box 111, FI-80101 Joensuu, Finland

###### Abstract

We propose a free-space, inverse design of nanostructure’s effective mode-
matching fields via a backward propagation of tightly focused vector beams to
the pupil plane of an aplanatic system of high numerical aperture. First, we
study the nanostructure’s eigenmodes without considering any excitation fields
and then extract the modal near fields in the focal plane. Each modal field is
then taken as the desired focal field, the band-limited waves of which are
backward propagated to the pupil plane via a reversal of the Richards–Wolf
vector diffraction formula. The pupil fields can be designed to be genuinely
paraxial by associating the longitudinal electric/magnetic field component
with the radial one on the reference sphere. The inversely designed pupil
field in turn is propagated forwardly into the focal region to generate the
designed focal field, whose distribution over the nanostructure’s surface is
used to evaluate the overlap between the designed focal field and the modal
fields, i.e., the modal expansion coefficients. Studies for a silicon nanodisk
monomer, dimer, and tetramer demonstrate the ability of our inverse approach
to design the necessary tightly focused vector field that can effectively and
exclusively match a certain eigenmode of interest. Compared with the forward
beam-shaping method, the inverse design approach tends to yield quantitatively
more precise mode-matching field profiles. This work can have a significant
impact on optical applications that rely on controllable and tunable mode
excitation and light scattering.

##  I Introduction

Strong local field in plasmonic or high-index dielectric nanostructures is
beneficial for boosting a variety of optical effects, via an enhanced
interaction between the strong light field and nanostructures. In contrast to
plasmonic systems, high-index dielectric nanoparticles, such as silicon
nanodisks [1, 2, 3, 4], do not suffer from strong optical losses but still
sustain deep sub-wavelength field confinements [5]. Moreover, all-dielectric
nanostructures have been demonstrated to be versatile semiconductor-compatible
nanophotonic devices for directional scattering [1], near-field mapping of
optical modes [6], tailoring second-harmonic [7, 8, 9] and third-harmonic [10]
generations, anapole-enhanced Raman scattering [3], etc. With recent advances,
all-dielectric nanostructures are promising to form basic building blocks that
complement or even replace the plasmonic counterparts in modern nanophotonics
[5].

At certain frequencies the local field enhancement achieves maximum values
(i.e., peaks), which is a phenomenon described as the optical resonance [11,
5]. It is known that the resonance response can be tuned by varying size,
shape, and dielectric environment of the nanostructures [12]. The resonance is
commonly analyzed with respect to a certain incident light field, such as a
plane wave used in previous studies [1, 13, 3], which, however, may not probe
all resonances [14]. Recently, it has been shown that vector beams [15, 16]
with spatially varying polarization distributions can reveal otherwise hidden
resonances, such as the dark resonance detected by the radially or azimuthally
polarized cylindrical vector beams [17, 18, 19], or the hybrid resonances by
higher-order vector beams possessing helical phase distributions [20, 21].

The optical resonance response in nanostructures is associated with the
excitation of one or more mode(s) [14] which, strictly speaking, are
independent of any external excitation field and hence are intrinsically
determined by the optical properties of the nanostructure and surrounding
media. The resonance frequencies are generally isolated poles in the complex
frequency plane, a characteristic of quasi-normal modes [22]. At a fixed real-
valued frequency, however, a set of discrete eigenmodes [14, 23] can be
defined. Modes in nanostructures have been investigated by various approaches
including the electrostatic method [24, 25, 26], and several rigorous full-
wave analyses based on the Fourier modal method [27, 28, 29], transfer matrix
method [30], finite element method [31, 22], Green’s tensor method [32], and
boundary element method (BEM) based on the surface integral equations (SIEs)
[14, 23, 33, 34], etc. Among them, BEM has been proven to be efficient
particularly in modelling the interaction between focused vector field and
nanostructures [35, 36, 21], because the numerical discretization as well as
the corresponding focal fields evaluation are reduced to two dimension (2D)
for a general three dimensional (3D) problem.

The excitation efficiency of an eigenmode scales with the overlap between the
modal and the applied incident field profiles, as well as the corresponding
eigenvalue [14]. An eigenmode’s near field in optical nanoantennas is
vectorial and localized such that it may vary considerably within a sub-
wavelength scale, which is a fact that poses a challenge to the choice of an
excitation field with an effective mode-matching profile. A free-space
approach to obtain such a mode-matching vector field is to shape a focal field
that is to match with the modal field, by focusing a paraxial vector beam in
an aplanatic system with high numerical aperture (NA). A forward design
approach has been shown to yield spatially varying field distributions for
efficient hybrid mode excitations [21], but a forwardly designed focal field
is limited to a qualitative mode matching due to energy exchanges between the
radial and azimuthal field components in the focusing process, and then a
deviation of its field profile from the designed pupil field’s. Alternatively,
an inverse design is versatile in shaping focal fields that potentially bring
more refined and effective mode-matching fields.

In the previous work of inverse problem in high-NA aplanatic system, the main
interest has been focused on the focal field total intensity or a single
component from which the complex pupil field is retrieved. For instance, the
depth of focus can be optimized after determining the abberration and
amplitude functions in high-NA imaging system [37], or even retrieving the
birefringence from four focal field intensity distributions generated by pupil
fields of different polarization states [38]. A focal field with a null
longitudinal component is generated by azimuthally polarized light [39].
Needle-, tube-, and bubble-shaped focal field intensity distributions [40],
perfect polarization vortex focal field [41], as well as focal field of
arbitrary homogeneous polarization [42] are also investigated. However, a
mode-matching focal field generally involves more than one or all field
component(s) to be effectively matched with both the amplitude and phase
distributions of the modal field, which has not been considered in previous
work.

In this work, we present an inverse design of the desired focal field for
effective mode excitation in optical nanoantennas by taking into account all
modal field components and the inherent dependencies of all electric and
magnetic field components in the designed beam-like pupil, reference, and
focal fields. Independent of any excitation field, the eigenmodes are fully
determined by the nanostructure’s optical properties themselves and thus they
are studied via a BEM mode solver implemented through the Müller formulation
of the SIEs [43]. The modal fields in the focal plane are evaluated via the
Stratton—Chu surface integral [44], and the extracted focal fields are
propagated back to the pupil plane by reversing the tight focusing process.
The best possible incident paraxial beams before focusing are then obtained
from solving an inverse problem governed by the Richards–Wolf vector
diffraction formula.

This work is organized as follows. In Sec. II, we study the eigenmodes in
nanostructures, where we review the SIEs, describe the eigenvalue problem
based on the Müller formulation, and analyze the eigenmodes in silicon
nanodisk monomer, dimer, and tetramer. In Sec. III, we describe the forward
and backward tight focusing processes which are governed by the Richards–Wolf
vector diffraction formula. The degrees of freedom and strategies in our
inverse design are then discussed through a systematical study of the mode-
matching field design in a silicon nanodisk monomer, dimer, and tetramer. In
Sec. IV, we analyze the quality of the designed focal field in terms of modal
expansion coefficients, i.e., the overlaps between the designed focal field
and all the first twelve eigenmodes considered in each nanodisk oligomer. At
last, conclusions and perspectives are given in Sec. V.

##  II Eigenmodes in nanostructures

###  II.1 Surface integral equations

We first describe the eigenmode problem in optical nanostructures using the
SIE formulation. Without loss of generality, the SIE formulation is reviewed
for the case of a single arbitrarily shaped 3D nanostructure, since it is
readily extended to more general cases of nano-object ensembles. As shown in
Fig. 1, the nano-object occupying domain $\Omega_{1}$ contains a homogeneous
and isotropic medium with permittivity $\epsilon_{1}$ and permeability
$\mu_{1}$. The free space in domain $\Omega_{0}$ has optical properties
$\epsilon_{0}$ and $\mu_{0}$. The monochromatic electric and magnetic fields
of harmonic time dependence $e^{-i\omega t}$ in the scatterer and free space
are denoted by $\\{\mathbf{E}_{1},\mathbf{H}_{1}\\}$ and
$\\{\mathbf{E}_{0},\mathbf{H}_{0}\\}$, respectively, where $\omega$ is the
angular frequency.

Figure 1: An arbitrarily-shaped 3D scatterer in free space. The scatterer that
occupies domain $\Omega_{1}$ has permittivity $\epsilon_{1}$ and permeability
$\mu_{1}$, whereas domain $\Omega_{0}$ is free space with optical properties
are $\epsilon_{0}$ and $\mu_{0}$. The electric and magnetic fields are
$\\{\mathbf{E}_{0},\mathbf{H}_{0}\\}$ in the free space and
$\\{\mathbf{E}_{1},\mathbf{H}_{1}\\}$ in the scatterer. Two domains’ surfaces
are correspondingly denoted by $\partial\Omega_{0}$ and $\partial\Omega_{1}$,
and their directions are determined by normal vectors $\hat{n}_{0}=\hat{n}$
and $\hat{n}_{1}=-\hat{n}$. Equivalent electric and magnetic surface current
densities $\\{\mathbf{J}_{0},\mathbf{M}_{0}\\}$ and
$\\{\mathbf{J}_{1},\mathbf{M}_{1}\\}$ are mathematically introduced on
$\partial\Omega_{0}$ and $\partial\Omega_{1}$, respectively.

Applying the field equivalence principle [45, 44, 46, 47, 48] to the
scatterer, we are interested in the fields
$\\{\mathbf{E}_{1},\mathbf{H}_{1}\\}$ in the domain $\Omega_{1}$ and thus the
fields elsewhere are not of interest. In this regard, the fields
$\\{\mathbf{E}_{0},\mathbf{H}_{0}\\}$ in the domain $\Omega_{0}$ can take any
values while the fields in $\Omega_{1}$ are kept fixed, which in turn implies
that the tangential fields over the surface $\partial\Omega_{1}$ of the domain
$\Omega_{1}$ should be unchanged, according to the uniqueness theorem [49, 50,
48]. In the equivalent problem for the scatterer, making the fields in
$\Omega_{0}$ trivial $\\{\mathbf{E}_{0}=0,\mathbf{H}_{0}=0\\}$ leads to an
arbitrariness in the $\epsilon_{0}$ and $\mu_{0}$ values, and a convenient
choice is to take the same values as in the scatterer with
$\epsilon_{0}=\epsilon_{1}$ and $\mu_{0}=\mu_{1}$, i.e., the full space medium
is homogeneous and isotropic. However, to maintain the tangential fields over
$\partial\Omega_{1}$, equivalent electric and magnetic surface current
densities must be introduced mathematically

$\displaystyle\mathbf{J}_{1}$ $\displaystyle=\hat{n}_{1}\times\mathbf{H}_{1},$
(1) $\displaystyle\mathbf{M}_{1}$
$\displaystyle=-\hat{n}_{1}\times\mathbf{E}_{1}.$ (2)

Applying the Stratton–Chu equation [44] in the equivalent problem for the
scatterer, the fields of interest at a point in the scatterer
($\mathbf{r}\in\Omega_{1}$) can be evaluated by

$\displaystyle\mathbf{E}_{1}(\mathbf{r})$
$\displaystyle=i\omega\mu_{1}\iint_{\partial\Omega_{1}}\left[\mathbf{J}_{1}\left(\mathbf{r}^{\prime}\right)G_{1}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\right]\mathrm{d}S^{\prime}$
$\displaystyle-\frac{1}{i\omega\epsilon_{1}}\iint_{\partial\Omega_{1}}\left[\nabla^{\prime}\cdot\mathbf{J}_{1}\left(\mathbf{r}^{\prime}\right)\nabla
G_{1}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\right]\mathrm{d}S^{\prime}$
$\displaystyle+\iint_{\partial\Omega_{1}}\left[\mathbf{M}_{1}\left(\mathbf{r}^{\prime}\right)\times\nabla
G_{1}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\right]\mathrm{d}S^{\prime},$
(3) $\displaystyle\mathbf{H}_{1}(\mathbf{r})$
$\displaystyle=-\iint_{\partial\Omega_{1}}\left[\mathbf{J}_{1}\left(\mathbf{r}^{\prime}\right)\times\nabla
G_{1}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\right]\mathrm{d}S^{\prime}$
$\displaystyle+i\omega\epsilon_{1}\iint_{\partial\Omega_{1}}\left[\mathbf{M}_{1}\left(\mathbf{r}^{\prime}\right)G_{1}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\right]\mathrm{d}S^{\prime}$
$\displaystyle-\frac{1}{i\omega\mu_{1}}\iint_{\partial\Omega_{1}}\left[\nabla^{\prime}\cdot\mathbf{M}_{1}\left(\mathbf{r}^{\prime}\right)\nabla
G_{1}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\right]\mathrm{d}S^{\prime},$
(4)

where
$G_{1}(\mathbf{r},\mathbf{r}^{\prime})=\exp(ik_{1}|\mathbf{r}-\mathbf{r}^{\prime}|)/(4\pi|\mathbf{r}-\mathbf{r}^{\prime}|)$
is Green’s function in the whole space of the equivalent problem, and $k_{1}$
is the wave number with $k_{1}^{2}=\epsilon_{1}\mu_{1}\omega^{2}$.

In the equivalent problem for the free space, the fields at a point
$\mathbf{r}\in\Omega_{0}$ can be evaluated similarly by

$\displaystyle\mathbf{E}_{0}(\mathbf{r})$
$\displaystyle=\mathbf{E}_{0}^{(i)}(\mathbf{r})+i\omega\mu_{0}\iint_{\partial\Omega_{0}}\left[\mathbf{J}_{0}\left(\mathbf{r}^{\prime}\right)G_{0}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\right]\mathrm{d}S^{\prime}$
$\displaystyle-\frac{1}{i\omega\epsilon_{0}}\iint_{\partial\Omega_{0}}\left[\nabla^{\prime}\cdot\mathbf{J}_{0}\left(\mathbf{r}^{\prime}\right)\nabla
G_{0}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\right]\mathrm{d}S^{\prime}$
$\displaystyle+\iint_{\partial\Omega_{0}}\left[\mathbf{M}_{0}\left(\mathbf{r}^{\prime}\right)\times\nabla
G_{0}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\right]\mathrm{d}S^{\prime},$
(5) $\displaystyle\mathbf{H}_{0}(\mathbf{r})$
$\displaystyle=\mathbf{H}_{0}^{(i)}(\mathbf{r})-\iint_{\partial\Omega_{0}}\left[\mathbf{J}_{0}\left(\mathbf{r}^{\prime}\right)\times\nabla
G_{0}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\right]\mathrm{d}S^{\prime}$
$\displaystyle+i\omega\epsilon_{0}\iint_{\partial\Omega_{0}}\left[\mathbf{M}_{0}\left(\mathbf{r}^{\prime}\right)G_{0}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\right]\mathrm{d}S^{\prime}$
$\displaystyle-\frac{1}{i\omega\mu_{0}}\iint_{\partial\Omega_{0}}\left[\nabla^{\prime}\cdot\mathbf{M}_{0}\left(\mathbf{r}^{\prime}\right)\nabla
G_{0}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\right]\mathrm{d}S^{\prime},$
(6)

where the extra terms $\mathbf{E}_{0}^{(i)}$ and $\mathbf{H}_{0}^{(i)}$
account for the incident fields in $\Omega_{0}$,
$\mathbf{J}_{0}=\hat{n}_{0}\times\mathbf{H}_{0}$,
$\mathbf{M}_{0}=-\hat{n}_{0}\times\mathbf{E}_{0}$ are the equivalent electric
and magnetic surface current densities which are introduced mathematically to
maintain the tangential fields on $\partial\Omega_{0}$, and $G_{0}$ is the
full-space Green function where the corresponding wave number is $k_{0}$ with
$k_{0}^{2}=\epsilon_{0}\mu_{0}\omega^{2}$.

To summarize, the fields in each region are compactly written as

$\displaystyle\mathbf{E}_{l}$
$\displaystyle=\delta_{l0}\mathbf{E}_{l}^{(i)}+\mathcal{D}_{l}\mathbf{J}_{l}-\mathcal{K}_{l}\mathbf{M}_{l},$
(7) $\displaystyle\mathbf{H}_{l}$
$\displaystyle=\delta_{l0}\mathbf{H}_{l}^{(i)}+\mathcal{K}_{l}\mathbf{J}_{l}+\left(\frac{\epsilon_{l}}{\mu_{l}}\right)\mathcal{D}_{l}\mathbf{M}_{l}.$
(8)

Above, $\delta_{l0}$ is the Kronecker delta and the integro-differential
operators are defined as follows,

$\displaystyle\left\\{\mathcal{D}_{l}\mathbf{f}\left(\mathbf{r}^{\prime}\right)\right\\}\left(\mathbf{r}\right)$
$\displaystyle=i\omega\mu_{l}\iint_{\partial\Omega_{l}}\left[\mathbf{f}\left(\mathbf{r}^{\prime}\right)G_{l}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\right]\mathrm{d}S^{\prime}$
$\displaystyle-\frac{1}{i\omega\epsilon_{l}}\iint_{\partial\Omega_{l}}\left[\nabla^{\prime}\cdot\mathbf{f}\left(\mathbf{r}^{\prime}\right)\nabla
G_{l}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\right]\mathrm{d}S^{\prime},$
(9)
$\displaystyle\left\\{\mathcal{K}_{l}\mathbf{f}\left(\mathbf{r}^{\prime}\right)\right\\}\left(\mathbf{r}\right)$
$\displaystyle=-\iint_{\partial\Omega_{l}}\left[\mathbf{f}\left(\mathbf{r}^{\prime}\right)\times\nabla
G_{l}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\right]\mathrm{d}S^{\prime},$
(10)

where $l\in\\{0,1\\}$ is the domain index. Therefore, for a given point on the
scatterer’s surface we have, after taking the operation $\hat{n}_{l}\times$ to
Eqs. (7) and (8),

$\displaystyle\mathbf{J}_{l}$
$\displaystyle=\delta_{l0}\mathbf{J}_{l}^{(i)}+\hat{n_{l}}\times\left[\mathcal{K}_{l}\mathbf{J}_{l}+\left(\frac{\epsilon_{l}}{\mu_{l}}\right)\mathcal{D}_{l}\mathbf{M}_{l}\right],$
(11) $\displaystyle\mathbf{M}_{l}$
$\displaystyle=\delta_{l0}\mathbf{M}_{l}^{(i)}-\hat{n_{l}}\times\left[\mathcal{D}_{l}\mathbf{J}_{l}-\mathcal{K}_{l}\mathbf{M}_{l}\right],$
(12)

where the current densities
$\mathbf{M}_{l}^{(i)}=-\hat{n}_{l}\times\mathbf{E}_{l}^{(i)}$ and
$\mathbf{J}_{l}^{(i)}=\hat{n}_{l}\times\mathbf{H}_{l}^{(i)}$ are due to the
incident fields in the free space.

The N-Müller formulation [14, 43] yields stable solutions even at low-
frequency limit and it can be obtained by combining the above two equations
with weights $\mu_{l}$ and $\epsilon_{l}$, as well as taking into account the
tangential field continuity at the interface (i.e.,
$\mathbf{J}=\mathbf{J}_{0}=-\mathbf{J}_{1}$ and
$\mathbf{M}=\mathbf{M}_{0}=-\mathbf{M}_{1}$),

$\displaystyle\left(\mu_{0}+\mu_{1}\right)\mathbf{J}$
$\displaystyle=\left(\mu_{0}\hat{n}_{0}\times\mathcal{K}_{0}+\mu_{1}\hat{n}_{1}\times\mathcal{K}_{1}\right)\mathbf{J}+\mu_{0}\mathbf{J}_{0}^{(i)}$
$\displaystyle+\left(\epsilon_{0}\hat{n}_{0}\times\mathcal{D}_{0}+\epsilon_{1}\hat{n}_{1}\times\mathcal{D}_{1}\right)\mathbf{M},$
(13) $\displaystyle\left(\epsilon_{0}+\epsilon_{1}\right)\mathbf{M}$
$\displaystyle=-\left(\epsilon_{0}\hat{n}_{0}\times\mathcal{D}_{0}+\epsilon_{1}\hat{n}_{1}\times\mathcal{D}_{1}\right)\mathbf{J}+\epsilon_{0}\mathbf{M}_{0}^{(i)}$
$\displaystyle+\left(\mu_{0}\hat{n}_{0}\times\mathcal{K}_{0}+\mu_{1}\hat{n}_{1}\times\mathcal{K}_{1}\right)\mathbf{M}.$
(14)

This can be further written in the compact form

$\left(\mathcal{Z}_{0}+\mathcal{Z}_{1}\right)\mathbf{x}=\mathbf{b},$ (15)

with the vector $\mathbf{x}=[\mathbf{J},\mathbf{M}]$ containing the equivalent
surface current densities to be solved with the knowledge of the tangential
components of the incident fields
$\mathbf{b}=[\mu_{0}\mathbf{J}_{0}^{(i)},\epsilon_{0}\mathbf{M}_{0}^{(i)}]$.
We also invoked the matrix notation

$\mathcal{Z}_{l}=\begin{bmatrix}\mu_{l}&0\\\
0&\epsilon_{l}\end{bmatrix}+\hat{n}_{l}\times\begin{bmatrix}-\mu_{l}\mathcal{K}_{l}&-\epsilon_{l}\mathcal{D}_{l}\\\
\epsilon_{l}\mathcal{D}_{l}&-\mu_{l}\mathcal{K}_{l}\end{bmatrix}.$ (16)

###  II.2 Eigenmodes

In the absence of an incident field $\mathbf{b}=0$ holds, and the system may
have nontrivial solutions at certain frequencies (usually complex-valued),
each of them being associated with a resonance mode of the optical system [14,
23, 33, 22]. The resonance frequencies are generally isolated poles in the
complex frequency plane. At a fixed real-valued frequency, however, we can
find the scatterer’s eigenmodes [14, 33] by solving

$\left(\mathcal{Z}_{0}+\mathcal{Z}_{1}\right)\mathbf{x}^{(m)}=\lambda^{(m)}\mathbf{x}^{(m)}.$
(17)

Rewriting Eq. (17) in the form of Eq. (15), we have

$\left(\mathcal{Z}_{0}+\mathcal{Z}_{1}\right)\frac{\mathbf{x}^{(m)}}{\lambda^{(m)}}=\mathbf{x}^{(m)}.$
(18)

Comparing Eqs. (15) and (18), it is seen that an incident field
$\mathbf{x}^{(m)}$, which is the $m$th eigenmode, yields a scattered field
$\mathbf{x}^{(m)}/\lambda^{(m)}$. In addition, an eigenmode associated with a
smaller eigenvalue $\lambda^{(m)}$ yields a stronger scattered field
$\mathbf{x}^{(m)}/\lambda^{(m)}$ when the scatterer is excited by an incident
field of $\mathbf{x}^{(m)}$. Therefore, $1/\lambda^{(m)}$ can be viewed as the
scattering strength associated with the $m$th eigenmode.

For a given eigenmode $\mathbf{x}^{(m)}=[\mathbf{J}^{(m)},\mathbf{M}^{(m)}]$,
we can calculate its electric and magnetic fields $\mathbf{E}^{(m)}$ and
$\mathbf{H}^{(m)}$ everywhere, by invoking Eqs. (7) and (8),

$\displaystyle\mathbf{E}_{l}^{(m)}(\mathbf{r})$
$\displaystyle=\mathcal{D}_{l}\mathbf{J}_{l}^{(m)}-\mathcal{K}_{l}\mathbf{M}_{l}^{(m)},$
(19) $\displaystyle\mathbf{H}_{l}^{(m)}(\mathbf{r})$
$\displaystyle=\mathcal{K}_{l}\mathbf{J}_{l}^{(m)}+\left(\frac{\epsilon_{l}}{\mu_{l}}\right)\mathcal{D}_{l}\mathbf{M}_{l}^{(m)},$
(20)

where the subscript $l$ denotes the domain where the field at a given point
$\mathbf{r}\in\Omega_{l}$ is evaluated.

###  II.3 Silicon nanodisks

#### Monomer

Figure 2: The reciprocal of absolute eigenvalues $1/|\lambda^{(m)}|$ of the
first twelve leading eigenmodes in a silicon nanodisk in air as a function of
wavelength. The nanodisk has a thickness of $h=220$ nm and a diameter of
$d=480$ nm and it is situated at the origin of the Cartesian coordinate system
(shown in the inset). Each eigenmode is labelled by M1 to M12, respectively.
In case of degeneracy, one of the two degenerate eigenmodes is labelled in
brackets. These twelve eigenmodes are ranked by their $1/|\lambda^{(m)}|$
values at wavelength $1300$ nm (the dashed line). In the filled grey area
below the zero horizontal line, the scattering cross-section
($\sigma_{\mathrm{sca}}$) in arbitrary units is shown as the black curve when
the nanodisk is illuminated by an $x$-polarized plane wave.

$\mathrm{MD}\downarrow$ |  $\mathrm{ED}\downarrow$ |  $\mathrm{MD}\varodot$ |  $\mathrm{MQ}_{\varotimes}^{\varodot}$ |  $\mathrm{EQ}_{+-}^{-+}$ |  $\mathrm{ED}\varotimes$ |  $\mathrm{MQ}_{+-}^{-+}$
---|---|---|---|---|---|---

Figure 3: Near fields including the instantaneous electric
$\widetilde{e}_{n}^{(m)}$ (first row) and magnetic $\widetilde{h}_{n}^{(m)}$
(second row) fields on the nanodisk’s surface, and far fields (third row) of a
few leading eigenmodes at the wavelength $1300$ nm for the same silicon
nanodisk as in Fig. 2. Each column is for one eigenmode with its
$1/|\lambda^{(m)}|$ value given on the top. Figure 4: Near and far fields of
eigenmodes M1–12 at 1600 nm in a dimer consisting of two identical silicon
nanodisks (same as in Fig. 2) with a center-to-center gap of $480$ nm. Figure
5: Near and far fields of a few leading eigenmodes M1–12 at 1600 nm in a
tetramer consisting of four identical silicon nanodisks (same as in Fig. 2)
with the center-to-center gap of two opposite neighbouring nanodisks being
$480$ nm.

For the purpose of demonstrating the eigenmodes in all-dielectric
nanostructures, we study a silicon nanodisk in air. The nanodisk has a
thickness of $h=220$ nm and a diameter of $d=480$ nm and it is situated at the
origin of the Cartesian coordinate system as shown in the inset of Fig. 2. The
eigenvalues are solved in the wavelength range of $[850,1750]$ nm and the
twelve dominant eigenmodes are shown in Fig. 2. The eigenmodes sharing the
same $1/|\lambda^{(m)}|$ values are degenerate modes, whose fields are of
identical distribution but distinctively oriented. Treating the degenerate
eigenmodes together, we observe seven unique eigenmodes showing different
resonance responses. Mode M1 (M2) has a resonance wavelength around 1330 nm.
Mode M3 (M4) does not show clear resonance as its scattering strength value is
flattened over a relatively wide spectral range. Mode M5 is a non-degenerate
mode showing resonance scattering at 1675 nm. Mode M6 (M7) has a relatively
weak resonance near 1200 nm. Resonances at two shorter wavelengths 970 nm and
1070 nm are associated with the modes M8 (M9) and M11 (M12), respectively. At
last, mode M10 is weakly resonant across the entire spectral range of
interest. At this point, we emphasize that the resonance wavelengths
determined from the scattering strength values are the genuine resonance
wavelengths that do not depend on any external incident field. The calculated
scattering cross section in arbitrary units from the nanodisk when illuminated
by an $x$-polarized plane wave propagating along the $z$-axis is shown as the
black curve in the filled grey area. The scattering cross section under the
illumination of a linearly polarized plane wave shows a resonance peak near
1300 nm (the dashed line), at which we have ranked the aforementioned
eigenmodes by their scattering strength values $1/|\lambda^{(m)}|$.

The corresponding near and far fields of the eigenmodes at the wavelength of
$1300$ nm are inspected in Fig. 3, and only one of degenerate eigenmodes, if
any, is considered. For visualizing both the electric and magnetic near fields
with balanced amplitudes [51, 52], we introduce the scaled instantaneous
electric and magnetic fields as

$\widetilde{\mathbf{e}}^{(m)}\left(\mathbf{r},t\right)=\sqrt{\epsilon_{0}}\mathrm{Re}\left\\{\mathbf{E}_{0}^{(m)}(\mathbf{r})e^{-i\omega
t}\right\\}$ (21)

and

$\widetilde{\mathbf{h}}^{(m)}\left(\mathbf{r},t\right)=\sqrt{\mu_{0}}\mathrm{Re}\left\\{\mathbf{H}_{0}^{(m)}(\mathbf{r})e^{-i\omega
t}\right\\},$ (22)

respectively, where $t$ is a time instant, $r\in\partial\Omega_{0}$, and
$\mathrm{Re}\\{\cdot\\}$ denotes the real part. It is noted that the
instantaneous electric and magnetic energy densities will be written as halves
the square of the scaled instantaneous electric and magnetic fields, i.e.,
$w_{e}=[\widetilde{\mathbf{e}}^{(m)}]^{2}/2$ and
$w_{h}=[\widetilde{\mathbf{h}}^{(m)}]^{2}/2$, respectively [53]. On nanodisk’s
surface, the normal component of the scaled field
$\widetilde{e}_{n}^{(m)}=\hat{n}\cdot\widetilde{\mathbf{e}}^{(m)}$ at a time
instant $t_{0}$ is shown for the eigenmode’s electric field, and the magnetic
counterpart is shown by
$\widetilde{h}_{n}^{(m)}=\hat{n}\cdot\widetilde{\mathbf{h}}^{(m)}$ at a
$\pi/2$-phase-delayed time instant $t_{0}^{\prime}$ with
$\omega(t_{0}^{\prime}-t_{0})=\pi/2$, at which time the magnetic field
generally reaches its maximum.

For a nanodisk monomer, the surface normal near fields
$\widetilde{e}_{n}^{(m)}$ and $\widetilde{h}_{n}^{(m)}$ are shown in the first
and second rows, respectively, in Fig. 3, and the eigenmodes’ corresponding
far fields are shown in the third row. Each column displays the near and far
fields for one eigenmode with its $1/|\lambda^{(m)}|$ value given on the top
of each column. As can be seen from the radiation patterns, the eigenmode M1
is a dipolar mode showing a typical doughnut shape. Furthermore, the near
fields show that the mode M1 is of magnetic dipole with a strong transverse
magnetic field, which is denoted by the symbol MD$\downarrow$ with the
downarrow representing a field vector in the transverse plane. In the second
column, M3 has a smaller $1/|\lambda^{(3)}|$ value, and the radiation pattern
and near fields indicate that it is an electric dipolar mode, which is denoted
by ED$\downarrow$. Under a similar observation, it turns out that M5 is a
magnetic dipolar mode that differs from M1 by a dominating $z$ magnetic field
component, denoted by MD$\varodot$ with the symbol “$\varodot$” implying a
field vector along the positive $z$ axis. Mode M6 emerges from a magnetic
quadrupole lying in a plane parallel to the $z$ axis, which is symbolized by
MQ${}_{\varotimes}^{\varodot}$ where the symbol “$\varotimes$” implies a field
vector along the negative $z$ axis that is $\pi$-out-of-phase with respect to
the “$\varodot$” field vector. The mode M6 is distinct from another magnetic
quadrupolar mode M11 where the quadrupole resides in the $xy$ plane, denoted
by MQ${}_{+-}^{-+}$ where the symbols “$-$” and “$+$” represent negative and
positive field extremes, respectively. An electric quadrupolar mode, denoted
by EQ${}_{+-}^{-+}$, is observed in mode M8. At last, mode M10 is an electric
dipolar mode, denoted by ED$\varotimes$, with a dominating electric $z$ field
component.

#### Oligomers

In silicon nanodisks ensembles, i.e., nanodisk oligomers consisting of several
nanodisks, further local field enhancements and mode hybridization happen via
optical near-field coupling. Here, we also investigate the near fields of a
few leading eigenmodes in nanostructure assemblies that consist of two and
four identical nanodisks, so-called nanodisk dimer and tetramer, respectively.

The eigenmodes’ near and far fields of a dimer are shown in Fig. 4. Modes M1
and M2 are hybridized from two longitudinal magnetic dipolar modes, where the
former arises from two parallel longitudinal MD modes, and the latter is due
to two anti-parallel longitudinal MD modes. Transversal magnetic modes include
M7-10, that are due to the hybridization of, subsequently, two $y$-oriented
anti-parallel transversal MD modes, two $x$-oriented anti-parallel transversal
MD modes, two $x$-oriented parallel transversal MD modes, and two $y$-oriented
parallel transversal MD modes. The electric counterparts of the magnetic modes
M7–10 are observed in transversal electric modes M3–6. In detail, M3 is
hybridized from two $x$-oriented anti-parallel transversal ED modes, M4 is due
to two $y$-oriented anti-parallel transversal ED modes, M5 arises from two
$y$-oriented parallel transversal ED modes, and two $x$-oriented parallel
transversal ED modes yields the hybrid mode M6. Additionally, two transveral
EQ modes, when arranged differently, hybridize into two distinct modes M11 and
M12.

In a tetramer, modes hybridization leads to more types of eigenmodes as a
result of more near-field coupling channels. Near and far fields of twelve
eigenmodes with leading scattering strengths $1/|\lambda^{(m)}|$ are shown in
Fig. 5. It is clear that mode M1 arises from four parallel longitudinal
magnetic modes, whereas the longitudinal magnetic fields in every two
neighbouring nanodisks in M4 are $\pi$-out-of-phase with respect to each
other. Modes hybridized from four transversal ED modes contain M7 where the ED
moments are $\pi$-out-of-phase along the radial direction between every two
neighbouring elements, M8 where the ED moments are $\pi$-out-of-phase between
every two neighbouring elements along the azimuthal direction, M9 with
radially in-phase ED moment in each element, and M12 with in-phase ED moment
in each nanodisk along the azimuthal direction. Hybridization of transversal
MDs leads to mode M10 with the MD moments being along the radial direction and
$\pi$-out-of-phase between the neighbouring nanodisks and M11 with azimuthally
in-phase MD moments. In addition, M10 and M11 are the magnetic counterparts of
modes M7 and M12. It is not straightforward to gain physical insights of the
origins of hybridized modes M2, M3, M5, and M6 without future investigation,
and determining the near-field coupling channels for these modes is less
important and not further discussed in this work.

As it can be seen from Figs. 3, 4, and 5, an eigenmode is associated with its
unique near-field distribution and far-field radiation pattern. Tuning light
scattered from nano-objects via selective excitation(s) of certain
eigenmode(s) could serve as a basic building block for engineering optical
near and/or far field(s). It is, therefore, vital to design an excitation
field which highly overlaps with the field of a specific eigenmode for an
effective mode excitation, i.e., mode-matching field.

##  III Inverse design via tight focusing

As shown in previous Sec. II, the eigenmode’s near fields in optical
nanoantennas are genuinely localized and vectorial at the nanoscale. A
potential free-space approach to design such a 3D mode-matching vector field
is to shape a focal field (that matches with the eigenmode’s near field in
optical nanoantennas that are to be placed in the vicinity of the focal
point), by focusing a paraxial vector beam in an aplanatic system with high
NA. In the previous work, tightly focused cylindrical vector beams of
fundamental order, i.e., radially and azimuthally polarized beams [18, 17],
and of higher-order polarization states [21] have been used for optical
excitation of dark and hybrid modes in optical nanoantennas. In these forward
design approaches, the focused field has a predefined field profile and only
qualitative match of polarization distribution between the hybrid modes and
the excitation fields was considered. To design the best possible focal field
that quantitatively matches with the eigenmode’s field, it is necessary to
solve an inverse problem where the needed paraxial vector beam to form a
desired focal field after tight focusing is to be sought.

###  III.1 Tight focusing: forward and backward propagation

Figure 6: The inverse design scheme of focused vector beams with effective
mode-matching profiles. (a) The forward and backward tight focusing processes.
The nanostructure’s surface is denoted by $\partial\Omega$, where the $m$th
eigenmode’s equivalent surface current density
$\\{\mathbf{J}^{(m)},\mathbf{M}^{(m)}\\}$ resides. The middle symmetry plane
$\Gamma$ coincides with the focal plane ($z=0$ plane), where the modal field
$\\{\mathbf{E}^{(m)},\mathbf{H}^{(m)}\\}$ is evaluated. The reference field
$\mathbf{E}_{\infty}$ resides on the reference spherical cap with a radius
being equal to the focal length $f$. The incident beam is denoted by
$\mathbf{E}_{\mathrm{inc}}$. The nanostructure occupies domain $\Omega_{1}$,
whereas the free space is denoted by $\Omega_{0}$ with a refractive index of
$n_{f}$. The index before the reference sphere is denoted by $n_{i}$. (b) The
inverse design flowchart. Solid black lines indicate connections between the
associated field components. The modal field
$\\{\mathbf{E}^{(m)},\mathbf{H}^{(m)}\\}$ on the focal plane $\Gamma$ serves
as the desired focal field, from which the the pupil field
$\\{\mathbf{E}_{\mathrm{inc}},\mathbf{H}_{\mathrm{inc}}\\}$ and thus the focal
field $\\{\mathbf{E}_{f},\mathbf{H}_{f}\\}$ are designed. Desired focal fields
are categorized into four groups by the time-averaged energies stored in the
transverse and longitudinal, electric and magnetic fields in the focal plane,
denoted by $\mathscr{W}_{e,\mathrm{tran}}^{(m)}$,
$\mathscr{W}_{e,\mathrm{long}}^{(m)}$, $\mathscr{W}_{h,\mathrm{tran}}^{(m)}$,
and $\mathscr{W}_{h,\mathrm{long}}^{(m)}$, respectively.

Before proceeding, we first review the forward and backward propagation in
tight focusing of vector beams. In detail, we choose a coordinate system such
that the middle plane of the nanostructure(s) coincides with the focal plane
$\Gamma$ (where $z=0$), as shown in Fig. 6. For a given eigenmode, its near
field $\\{\mathbf{E}^{(m)},\mathbf{H}^{(m)}\\}$ is evaluated from Eqs. (19)
and (20) in the focal plane $\Gamma$. For the eigenmode’s electric near field
to be mode-matched, a 3D vector field $\mathbf{E}_{f}$ is designed via the
tight focusing of a paraxial vector beam $\mathbf{E}_{\mathrm{inc}}$ at the
pupil plane (which is mapped to $\mathbf{E}_{\infty}$ on the reference
sphere).

The tight focusing process is governed by the Richards–Wolf formalism [54, 55,
56], and the focal field at a point $\mathbf{R}_{f}=(\rho,\varphi,z)$ in
cylindrical coordinates writes as

$\displaystyle\mathbf{E}_{f}(\rho,\varphi,z)$
$\displaystyle=\frac{-ik}{2\pi}\int\limits_{0}^{\theta_{\mathit{\mathrm{max}}}}\int\limits_{0}^{2\pi}\mathbf{E}_{\infty}(\theta,\phi)$
$\displaystyle\quad\times
e^{i\mathbf{k}\cdot\mathbf{R}_{f}}\sin\theta\mathrm{d}\phi\mathrm{d}\theta,$
(23)

where the wave vector $\mathbf{k}$ is along the direction from a point on the
reference sphere $(\phi,\theta)$ to the focal point and the corresponding wave
number $k=k_{0}$ (i.e., the same value as in domain $\Omega_{0}$). The maximum
converging angle $\theta_{\mathrm{max}}=\arcsin(\mathrm{NA}/n_{f})$ is limited
by the $\mathrm{NA}$ of the aplanatic system, where $n_{f}$ is the refractive
index after the reference sphere. We can rewrite Eq. (23) into Fourier
transform in Cartesian coordinates [57, 58]

$\displaystyle\mathbf{E}_{f}(x,y,z)$
$\displaystyle=\frac{-ik}{2\pi}\int\limits_{-\infty}^{+\infty}\int\limits_{-\infty}^{+\infty}\left[\frac{\mathbf{E}_{\infty}(k_{x},k_{y})e^{ikz\cos\theta}}{k^{2}\cos\theta}\right]$
$\displaystyle\quad\times
e^{-ik_{x}x}e^{-ik_{y}y}\mathrm{d}k_{x}\mathrm{d}k_{y},$ (24)

where $k_{x}=-k\sin\theta\cos\phi$ and $k_{y}=-k\sin\theta\sin\phi$ are the
$x$ and $y$ components of the wave number, respectively. Here, the surface
integral runs from $-\infty$ to $+\infty$, as the field $\mathbf{E}_{\infty}$
on the reference sphere is zero outside of the NA. In order to obtain the
necessary paraxial vector beam or the field $\mathbf{E}_{\infty}$ on the
reference sphere to generate the desired 3D mode-matching vector field
$\mathbf{E}_{f}$, we perform the inverse Fourier transform for the above Eq.
(24) except $\mathbf{E}_{f}$ is replaced by the eigenmode’s near field in the
focal plane $\mathbf{E}^{(m)}|_{\Gamma}$,

$\displaystyle\frac{\mathbf{E}_{\infty}\left(k_{x},k_{y}\right)e^{ikz\cos\theta}}{\cos\theta}$
$\displaystyle=\frac{ik}{2\pi}\int\limits_{-\infty}^{+\infty}\int\limits_{-\infty}^{+\infty}\left.\mathbf{E}^{(m)}\right|_{\Gamma}$
$\displaystyle\quad\times
e^{ik_{x}x}e^{\mathrm{i}k_{y}y}\mathrm{d}x\mathrm{d}y,$ (25)

which can be considered as the governing equation for a reversed Richards–Wolf
tight focusing process.

Ideally, the paraxial input vector beam is inversely restored from the
reference field [59, 21]

$\displaystyle\mathbf{E}_{\mathrm{inc}}(r,\phi)$
$\displaystyle=\mathbf{L}^{\mathrm{-1}}(\theta)\mathbf{E}_{\infty}(\theta,\phi),$
(26)

where the radius is $r=f\sin\theta$ with $f$ being the focal length and the
mapping matrix in cylindrical coordinates is [21]

$\mathbf{L}^{\mathrm{-1}}(\theta)=\frac{1}{f\sqrt{\cos\theta}}\sqrt{\frac{n_{f}}{n_{i}}}\begin{bmatrix}\cos\theta&0&\sin\theta\\\
0&1&0\\\ -\sin\theta&0&\cos\theta\end{bmatrix},$ (27)

where $n_{i}$ is the refractive index before the reference sphere. Above, we
followed Richards and Wolf [55] for the refraction at the aplanatic lens, the
reference field amplitude vector writes as
$\mathbf{E}_{\infty}(\theta,\phi)=f\sqrt{n_{i}/n_{f}}\sqrt{\cos\theta}[E_{\mathrm{inc},r}\hat{\theta}+E_{\mathrm{inc},\phi}\hat{\phi}+E_{\mathrm{inc},z}(\hat{\theta}\times\hat{\phi})]$.

###  III.2 Degrees of freedom

The associated magnetic fields are governed by similar formulas as in Eqs.
(23)–(26) for the electric counterparts, but for brevity we do not repeat them
here. However, it is noted that in Richards–Wolf formalism the magnetic fields
are linked to electric fields by revoking the relations
$\mathbf{H}_{\mathrm{inc}}=1/Z_{i}(\mathbf{k}_{i}/k_{i})\times\mathbf{E}_{\mathrm{inc}}$
and thus
$\mathbf{H}_{\infty}=1/Z_{f}(\mathbf{k}_{f}/k_{f})\times\mathbf{E}_{\infty}$
[54] as a result of beam-like characteristic of the vector fields on the pupil
plane and on the reference sphere. Here, $Z_{i}$ and $Z_{f}$ are the
associated wave impedances. This implies that only two components in the pupil
vector field are completely independent. In cylindrical coordinates, four
possible combinations of two field components are
$(E_{\mathrm{inc},r},E_{\mathrm{inc},\phi})$,
$(H_{\mathrm{inc},r},H_{\mathrm{inc},\phi})$,
$(E_{\mathrm{inc},r},H_{\mathrm{inc},r})$, and
$(E_{\mathrm{inc},\phi},H_{\mathrm{inc},\phi})$, and choosing any pair will
fully determine the vector field on the pupil plane, then on the reference
sphere, and finally in the focal region.

The fact that the incident vector beam is paraxial ($E_{\mathrm{inc},z}\approx
0$ and $H_{\mathrm{inc},z}\approx 0$) imposes a restriction on the radial and
longitudinal components of the reference electric [58] and magnetic fields via
the following relations:

$\displaystyle-\sin\theta E_{\infty,r}(\theta,\phi)+\cos\theta
E_{\infty,z}(\theta,\phi)$ $\displaystyle=0,$ (28) $\displaystyle-\sin\theta
H_{\infty,r}(\theta,\phi)+\cos\theta H_{\infty,z}(\theta,\phi)$
$\displaystyle=0.$ (29)

As a consequence, the focused electric field $\mathbf{E}_{f}$ bears all three
components that are nevertheless not completely independent [39]. In other
words, all three components of the electric focal field are determined
uniquely by two components of the beam-like pupil field
$\mathbf{E}_{\mathrm{inc}}$, as seen in Fig. 6(b) where the black solid lines
indicate the dependencies of field components on different planes. In order to
inversely design a focal field that is truly generated by a beam-like vector
pupil field, inherent dependencies of different field components of the
desired focal field should be carefully investigated.

###  III.3 Inverse design strategies

To inversely design a focal field that potentially matches with an eigenmode’s
field, i.e., the desired focal field, in principle all the electric and
magnetic modal field components should be taken into account. However, bearing
in mind the degrees of freedom discussed in the previous section, an
inconsistent pupil field could result from a backward propagation that ignores
the inherent dependencies of all components [39] in the desired focal field.
In contrast, one has to choose appropriate components in the desired focal
field that will unambiguously determine the beam-like pupil field and thus the
designed focal field. It is observable in Figs. 3, 4, and 5 that only some of
the electric and magnetic field components are dominant in an eigenmode’s near
field. For instance, the mode M10 of a single nanodisk in Fig. 3 shows strong
longitudinal electric field component. Moreover, in Fig. 4 modes M3–6 of a
nanodisk dimer are clearly dominated by their transverse electric and
longitudinal magnetic field components, whereas the transverse magnetic and
longitudinal electric field components are dominant in modes M7–10 of a
nanodisk dimer.

For a more quantitative analysis, we categorize the desired focal fields
according to the time-averaged energies associated with the transverse and
longitudinal electric and magnetic fields in the focal plane, i.e., the
following surface integrals of the energy densities [60],

$\displaystyle\mathscr{W}_{e,\mathrm{tran}}^{(m)}$
$\displaystyle=\frac{1}{2}\epsilon_{0}\int_{\Gamma}\left(\left|E_{\rho}^{(m)}\right|^{2}+\left|E_{\varphi}^{(m)}\right|^{2}\right)\mathrm{d}S,$
(30) $\displaystyle\mathscr{W}_{e,\mathrm{long}}^{(m)}$
$\displaystyle=\frac{1}{2}\epsilon_{0}\int_{\Gamma}\left|E_{z}^{(m)}\right|^{2}\mathrm{d}S,$
(31) $\displaystyle\mathscr{W}_{h,\mathrm{tran}}^{(m)}$
$\displaystyle=\frac{1}{2}\mu_{0}\int_{\Gamma}\left(\left|H_{\rho}^{(m)}\right|^{2}+\left|H_{\varphi}^{(m)}\right|^{2}\right)\mathrm{d}S,$
(32) $\displaystyle\mathscr{W}_{h,\mathrm{long}}^{(m)}$
$\displaystyle=\frac{1}{2}\mu_{0}\int_{\Gamma}\left|H_{z}^{(m)}\right|^{2}\mathrm{d}S.$
(33)

In principle, the surface integration is performed over an infinite focal
plane $\Gamma$, but in numerical evaluation we have to truncate the
integration surface in an area where the desired focal field remains
significant and discard the field elsewhere.

In group \tinyI⃝, the focal-plane energy associated with the transverse
electric field is larger than those of the longitudinal electric field
component as well as both focal-plane energies related to the transverse and
longitudinal magnetic field components. In group \tinyII⃝, the focal-plane
energy of the longitudinal electric field overwhelms the other three types of
energies. The magnetic counterparts of groups \tinyI⃝ and \tinyII⃝ are group
\tinyIII⃝ and group \tinyIV⃝, respectively. Concisely, we may write

$\displaystyle\mathrm{max}(\mathscr{W}_{e,\mathrm{tran}}^{(m)},\mathscr{W}_{e,\mathrm{long}}^{(m)},\mathscr{W}_{h,\mathrm{tran}}^{(m)},\mathscr{W}_{h,\mathrm{long}}^{(m)})$
$\displaystyle=$
$\displaystyle\begin{cases}\mathscr{W}_{e,\mathrm{tran}}^{(m)},&\text{in group
\tinyI⃝},\\\ \mathscr{W}_{e,\mathrm{long}}^{(m)},&\text{in group \tinyII⃝},\\\
\mathscr{W}_{h,\mathrm{tran}}^{(m)},&\text{in group \tinyIII⃝},\\\
\mathscr{W}_{h,\mathrm{long}}^{(m)},&\text{in group \tinyIV⃝}.\end{cases}$
(34)

For eigenmodes belonging to the group \tinyI⃝, the corresponding transverse
electric field is chosen as the candidate component of the desired focal
field, upon which the radial and azimuthal components of the reference field
$(E_{\infty,r},E_{\infty,\phi})$ and thus the pupil field
$(E_{\mathrm{inc},r},E_{\mathrm{inc},\phi})$ are determined. The longitudinal
electric focal field is, however, a free parameter so that the longitudinal
electric field component on the reference sphere $E_{\infty,z}$ is determined
from the radial component $E_{\infty,r}$ according to the restriction in Eq.
(28). In this regard, the pupil vector field is guaranteed to be beam-like,
and the designed focal field including the magnetic part is fully determined
by the pair of field components $(E_{\mathrm{inc},r},E_{\mathrm{inc},\phi})$.

For eigenmodes in the group \tinyII⃝, the longitudinal electric field
component of the desired focal field $E_{z}^{(m)}$ determines the longitudinal
electric reference field $E_{\infty,z}$, then the radial electric pupil field
$E_{\mathrm{inc},r}$ using the restriction in Eq. (28) where the radial
electric reference field is a free parameter. What remains to be determined is
either the azimuthal electric or the radial magnetic pupil field. But, the
azimuthal electric pupil field is only determined by the azimuthal electric
reference field, which in turn depends on the radial and azimuthal desired
electric focal field both of which are usually trivial for an eigenmode
showing strong longitudinal electric field. On the other hand, a strong
longitudinal electric field component of the desired focal field is usually
accompanied by slightly less stronger transverse magnetic field component
$H_{\rho}^{(m)}$ and/or $H_{\varphi}^{(m)}$ from which the radial magnetic
reference field $H_{\infty,r}$ and then the radial magnetic pupil field
$H_{\mathrm{inc},r}$ can be determined. The magnetic counterparts of the
aforementioned electric design strategies apply to the groups \tinyIII⃝ and
\tinyIV⃝. At this point, it is worth pointing out that the restriction on the
field components $E_{\infty,r}$ and $E_{\infty,z}$, as well as $H_{\infty,r}$
and $H_{\infty,z}$ poses a challenge to design a 3D focused vector beam, whose
electric and magnetic field components perfectly match the eigenmode’s near
field.

Besides the constraint imposed by the beam-like characteristics in the pupil
field, another limitation in the inverse design of a focal field is due to the
fact that the reference field is band-limited, i.e., spatial frequencies that
are larger than the free-space wave number are evanescent and cannot be
backward propagated onto the reference sphere. This is associated with the
fact that the tight focusing process relies on free-space beam propagation,
and the focal field is actually the far field propagated from the secondary
field on the reference sphere. In addition, back propagated waves are further
limited by the NA of the aplanatic system.

Figure 7: Scaled electric and magnetic field components [defined in Eqs. (21)
and (22)] at different stages in the inverse design for M2 in the nanodisk
dimer whose eigenmodes are shown in Fig. 4. The dashed circles represent the
nanodisk’s outlines in the desired (first column) and designed (fourth column)
focal fields, or the threshold of maximum free-space wave number $k_{0}$ in
$\mathbf{k}$ space of the desired focal field (second column) and the designed
pupil field (third column). The solid circle in the designed pupil field
indicates the wave number limited by NA with a value of 0.9 used in this
inverse design.

An example eigenmode belonging to group \tinyI⃝ is the mode M2 in the nanodisk
dimer, and its scaled electric and magnetic fields at different stages of the
inverse design are shown in Fig. 7. The first and second columns show all six
field components of the desired focal field in the real space and spatial
frequency $\mathbf{k}$ space, respectively. The inversely designed pupil
field’s components are shown in the third column, and the designed focal
field’s components can be seen in the fourth column. Both the radial and
azimuthal scaled electric field components of the designed focal field
effectively match with the desired eigenmode’s near field, even though the
desired focal field is so locally confined that a considerable amount of light
is at high spatial frequencies outside the $k_{0}$-circle and has to be
discarded in the free-space backward propagation of the inverse design.
Imposed by the beam-like characteristics of the designed pupil field, the
azimuthal scaled magnetic field is identical to the radial scaled electric
field, whereas the radial scaled magnetic field is $\pi$-phase shifted with
respect to the azimuthal scaled electric field. Four other field components
are also present as a result of tight focusing of the beam-like pupil field
and they are dependent on the radial and azimuthal electric field components
[see Fig. 6(b) for the connections of field components]. In particular, the
designed longitudinal magnetic field is $\pi$-out-of-phase with respect to the
desired longitudinal magnetic field. Nevertheless, as it will be shown later
with the expansion coefficients in the next section, this designed focal field
still effectively excites the mode M2 in the nanodisk dimer.

Figure 8: Scaled electric and magnetic field components in the inverse design
for M7 in the nanodisk tetramer (similar to Fig. 7) whose eigenmodes are shown
in Fig. 5.

Another example is shown in Fig. 8 for inversely designing a focal field to
match with the eigenmode M7 in the nanodisk tetramer. An effective match
between the designed and the desired focal fields is achieved in radial and
azimuthal electric field components, as well as the longitudinal magnetic
field component. Again, as it will be shown later with the expansion
coefficients in the next section, this designed focal field can effectively
and exclusively excite the mode M7 in the nanodisk tetramer. In this case, it
seems dropping waves with spatial frequencies higher than those limited by NA
has a trivial effect on the excitation efficiency. A full list of scaled field
components in the inverse design of focal field for all twelve eigenmodes in
the nanodisk monomer, dimer, and tetramer is compiled in the Supplementary
Material (SM).

##  IV Mode expansion coefficients

Figure 9: The mode excitation coefficients between the designed focal fields
$\mathbf{E}_{f}^{(m)}$ (along the horizontal axis) and eigenmodes
$\mathbf{E}^{(n)}$ (along the vertical axis) for the nanodisk (a) monomer, (b)
dimer, and (c) tetramer. For each designed focal field, the expansion
coefficients for the first twelve eigenmodes are normalized to the maximum
value.

Following the analysis of eigenmodes in nano-objects (Sec. II) and the inverse
design of tightly focused beam (Sec. III), it is important to see how well a
designed focal field matches with an eigenmode’s near field. A figure of merit
is the modal expansion coefficient which characterizes the overlapping between
the inversely designed vector beam and the eigenmode’s near field. Under the
SIE formulation, the modal expansion coefficient can be evaluated via the
inner product of two light fields over the nanostructures’ surfaces
$\partial\Omega_{0}$ [14],

$a^{(m,n)}=\frac{\iint_{\partial\Omega_{0}}\mathbf{x}_{L}^{(n)}\cdot\mathbf{b}^{(m)}\mathrm{d}S}{\iint_{\partial\Omega_{0}}\mathbf{x}_{L}^{(n)}\cdot\mathbf{x}^{(n)}\mathrm{d}S},$
(35)

where $\mathbf{x}_{L}^{(n)}$ and $\mathbf{x}^{(n)}$ are the left and right
eigenvectors associated with the $n$th eigenmode with an eigenvalue
$\lambda^{(n)}$, and $\mathbf{b}^{(m)}$ is the vector representation of the
focal field inversely designed from the $m$th eigenmode and evaluated on the
surface of the nanodisk(s) $\partial\Omega_{0}$. In detail, the left
eigenvector is defined by, in comparison with Eq. (17),

$\mathbf{x}_{L}^{(n)}\left(\mathcal{Z}_{0}+\mathcal{Z}_{1}\right)=\lambda^{(n)}\mathbf{x}_{L}^{(n)},$
(36)

and the vector of the designed focal field is

$\mathbf{b}^{(m)}=\left[\mu_{0}\mathbf{J}_{f}^{(m)},\epsilon_{0}\mathbf{M}_{f}^{(m)}\right],$
(37)

where $\mathbf{J}_{f}^{(m)}=\hat{n}_{0}\times\mathbf{H}_{f}^{(m)}$ and
$\mathbf{M}_{f}^{(m)}=-\hat{n}_{0}\times\mathbf{E}_{f}^{(m)}$ are the electric
and magnetic surface current densities arising from the designed focal field
over the surface $\partial\Omega_{0}$.

The eigenmode’s near field is solved directly in Eq. (17) over the surface
$\partial\Omega_{0}$. The designed focal field on $\partial\Omega_{0}$ can be
evaluated via the surface double integral in Eq. (23), where the reference
field $\mathbf{E}_{\infty}$ in the integrand is calculated through Eq. (25) by
fast Fourier transform (FFT) [61]. As a result, the reference field
$\mathbf{E}_{\infty}$ is a 2D matrix on a square grid, and we need to use Eq.
(23) rather than Eq. (24) for the focal field calculation on
$\partial\Omega_{0}$. Here, the double integral is implemented through the 2D
numerical integration based on the trapezoidal rule [62].

The expansion coefficients for the nanodisk monomer, dimer, and tetramer are
visualized in three colormap matrices in Figs. 9(a), (b), and (c),
respectively. The $m$th focal field $\mathbf{E}_{f}^{(m)}$ is inversely
designed from the corresponding $m$th eigenmode’s near field in the focal
plane. For the monomer, the designed focal fields $\mathbf{E}_{f}^{(5)}$,
$\mathbf{E}_{f}^{(8)}$, $\mathbf{E}_{f}^{(9)}$, and $\mathbf{E}_{f}^{(10)}$
largely overlap with the eigenmodes from which they are designed,
respectively. In other words, these eigenmodes can be effectively mode-matched
by tightly focused vector beams via our inverse design approach. However, the
designed focal fields $\mathbf{E}_{f}^{(m)}$ with $m=1,2,11,\text{ and }12$ do
not effectively match the corresponding desired eigenmodes and instead other
eigenmodes are strongly excited. The designed focal fields
$\mathbf{E}_{f}^{(m)}$ with $m=3,4,6,\text{ and }7$ strongly, but not
exclusively (since other eigenmodes are also excited), excite the desired
corresponding eigenmodes. In the nanodisk dimer case, the effective mode-
matching focal fields can be obtained for the eigenmodes $\mathbf{E}^{(n)}$
with $n=1\text{ to }6$, but the focal fields $\mathbf{E}_{f}^{(m)}$ with
$m=7\text{ to }12$ jump to near fields that effectively match with the
eigenmodes $\mathbf{E}^{(3)}$, $\mathbf{E}^{(4)}$, $\mathbf{E}^{(5)}$,
$\mathbf{E}^{(6)}$, $\mathbf{E}^{(5)}$, and $\mathbf{E}^{(6)}$, respectively.
For the nanodisk tetramer, the effective mode-matching focal fields are
obtained for the eigenmodes $\mathbf{E}^{(n)}$ with $n=2,3,5,6,7,8,\text{ and
}9$. However, the designed focal fields $\mathbf{E}_{f}^{(10)}$ and
$\mathbf{E}_{f}^{(11)}$ will excite the eigenmodes $\mathbf{E}^{(8)}$ and
$\mathbf{E}^{(9)}$, respectively. The designed focal fields
$\mathbf{E}_{f}^{(1)}$ and $\mathbf{E}_{f}^{(12)}$ always excite two
eigenmodes simultaneously, including the one from which the individual focal
field is inversely designed.

##  V Conclusion and discussion

We have developed an inverse design approach of generating a focused light
field that potentially matches with an eigenmode’s near field in optical
nanoantennas. We began with a rigorous analysis of eigenmodes in the the
optical nanoantennas of interest. In this work, a nanodisk monomer, dimer, and
tetramer were considered, and the first twelve eigenmodes in each
nanostructure were studied using our BEM eigenmode solver that is based on the
N-Müller formulation of the SIEs. In our inverse design strategies, the
eigenmode’s near field in the nanodisk’s middle symmetric plane (which is
chosen to coincide with the focal plane) is set as the desired focal field,
which is inversely propagated to the pupil plane where the designed pupil
field is obtained. In this backward propagation, we took into account the
beam-like characteristic of the pupil field and the fact the pupil field is
band-limited, as well as the inherent dependencies of different field
components. The eigenmodes were categorized into four groups by the time-
averaged energies associated with the transverse and longitudinal electric and
magnetic field in the focal plane, and the dominant field components were
chosen for uniquely determining the beam-like pupil field. The designed focal
field’s expansion coefficients into the first twelve eigenmodes were
evaluated. For several eigenmodes, effective and exclusive mode excitations
were achieved by designing the focal field accordingly from the eigenmodes via
our developed inverse design method. This work can have a significant impact
on optical switching, near- and far-field engineering in nano-optics.

## Acknowledgments

The authors thank Martti Kauranen for valuable advice on improving the
manuscript. This work was financially supported by the Academy of Finland
[projects 308393 and 320166 (PREIN)].

## References

  * Staude  _et al._ [2013] I. Staude, A. E. Miroshnichenko, M. Decker, N. T. Fofang, S. Liu, E. Gonzales, J. Dominguez, T. S. Luk, D. N. Neshev, I. Brener, and Y. Kivshar, Tailoring directional scattering through magnetic and electric resonances in subwavelength silicon nanodisks, ACS Nano 7, 7824 (2013).
  * Chong  _et al._ [2015] K. E. Chong, I. Staude, A. James, J. Dominguez, S. Liu, S. Campione, G. S. Subramania, T. S. Luk, M. Decker, D. N. Neshev, I. Brener, and Y. S. Kivshar, Polarization-independent silicon metadevices for efficient optical wavefront control, Nano Letters 15, 5369 (2015).
  * Baranov  _et al._ [2018] D. G. Baranov, R. Verre, P. Karpinski, and M. Käll, Anapole-enhanced intrinsic Raman scattering from silicon nanodisks, ACS Photonics 5, 2730 (2018).
  * Zhao and Reinhard [2019] X. Zhao and B. M. Reinhard, Switchable chiroptical hot-spots in silicon nanodisk dimers, ACS Photonics 6, 1981 (2019).
  * Kuznetsov  _et al._ [2016] A. I. Kuznetsov, A. E. Miroshnichenko, M. L. Brongersma, Y. S. Kivshar, and B. Luk’yanchuk, Optically resonant dielectric nanostructures, Science 354, aag2472 (2016).
  * Habteyes  _et al._ [2014] T. G. Habteyes, I. Staude, K. E. Chong, J. Dominguez, M. Decker, A. Miroshnichenko, Y. Kivshar, and I. Brener, Near-field mapping of optical modes on all-dielectric silicon nanodisks, ACS Photonics 1, 794 (2014).
  * Cambiasso  _et al._ [2017] J. Cambiasso, G. Grinblat, Y. Li, A. Rakovich, E. Cortés, and S. A. Maier, Bridging the gap between dielectric nanophotonics and the visible regime with effectively lossless GaP antennas, Nano Letters 17, 1219 (2017).
  * Camacho-Morales  _et al._ [2019] R. Camacho-Morales, G. Bautista, X. Zang, L. Xu, L. Turquet, A. Miroshnichenko, H. H. Tan, A. Lamprianidis, M. Rahmani, C. Jagadish, D. N. Neshev, and M. Kauranen, Resonant harmonic generation in AlGaAs nanoantennas probed by cylindrical vector beams, Nanoscale 11, 1745 (2019).
  * Sautter  _et al._ [2019] J. D. Sautter, L. Xu, A. E. Miroshnichenko, M. Lysevych, I. Volkovskaya, D. A. Smirnova, R. Camacho-Morales, K. Zangeneh Kamali, F. Karouta, K. Vora, H. H. Tan, M. Kauranen, I. Staude, C. Jagadish, D. N. Neshev, and M. Rahmani, Tailoring second-harmonic emission from (111)-GaAs nanoantennas, Nano Letters 19, 3905 (2019).
  * Shcherbakov  _et al._ [2014] M. R. Shcherbakov, D. N. Neshev, B. Hopkins, A. S. Shorokhov, I. Staude, E. V. Melik-Gaykazyan, M. Decker, A. A. Ezhov, A. E. Miroshnichenko, I. Brener, A. A. Fedyanin, and Y. S. Kivshar, Enhanced third-harmonic generation in silicon nanoparticles driven by magnetic response, Nano Letters 14, 6488 (2014).
  * Maier and Atwater [2005] S. A. Maier and H. A. Atwater, Plasmonics: Localization and guiding of electromagnetic energy in metal/dielectric structures, Journal of Applied Physics 98, 011101 (2005).
  * Kelly  _et al._ [2003] K. L. Kelly, E. Coronado, L. L. Zhao, and G. C. Schatz, The optical properties of metal nanoparticles: The influence of size, shape, and dielectric environment, The Journal of Physical Chemistry B 107, 668 (2003).
  * Moreno  _et al._ [2013] F. Moreno, P. Albella, and M. Nieto-Vesperinas, Analysis of the spectral behavior of localized plasmon resonances in the near-and far-field regimes, Langmuir 29, 6715 (2013).
  * Mäkitalo  _et al._ [2014] J. Mäkitalo, M. Kauranen, and S. Suuriniemi, Modes and resonances of plasmonic scatterers, Physical Review B 89, 165429 (2014).
  * Zhan [2009] Q. Zhan, Cylindrical vector beams: From mathematical concepts to applications, Advances in Optics and Photonics 1, 1 (2009).
  * Rosales-Guzmán  _et al._ [2018] C. Rosales-Guzmán, B. Ndagano, and A. Forbes, A review of complex vector light fields and their applications, Journal of Optics 20, 123001 (2018).
  * Sancho-Parramon and Bosch [2012] J. Sancho-Parramon and S. Bosch, Dark modes and Fano resonances in plasmonic clusters excited by cylindrical vector beams, ACS Nano 6, 8415 (2012).
  * Gómez  _et al._ [2013] D. E. Gómez, Z. Q. Teo, M. Altissimo, T. J. Davis, S. Earl, and A. Roberts, The dark side of plasmonics, Nano Letters 13, 3722 (2013).
  * Yanai  _et al._ [2014] A. Yanai, M. Grajower, G. M. Lerman, M. Hentschel, H. Giessen, and U. Levy, Near-and far-field properties of plasmonic oligomers under radially and azimuthally polarized light excitation, ACS Nano 8, 4969 (2014).
  * Reich  _et al._ [2020] S. Reich, N. S. Mueller, and M. Bubula, Selection rules for structured light in nanooligomers and other nanosystems, ACS Photonics 7, 1537 (2020).
  * Zang  _et al._ [2021] X. Zang, G. Bautista, L. Turquet, T. Setälä, M. Kauranen, and J. Turunen, Efficient hybrid-mode excitation in plasmonic nanoantennas by tightly focused higher-order vector beams, Journal of the Optical Society of America B 38, 521 (2021).
  * Lalanne  _et al._ [2018] P. Lalanne, W. Yan, K. Vynck, C. Sauvan, and J.-P. Hugonin, Light interaction with photonic and plasmonic resonances, Laser & Photonics Reviews 12, 1700113 (2018).
  * Powell [2014] D. A. Powell, Resonant dynamics of arbitrarily shaped meta-atoms, Physical Review B 90, 075108 (2014).
  * Fredkin and Mayergoyz [2003] D. R. Fredkin and I. D. Mayergoyz, Resonant behavior of dielectric objects (electrostatic resonances), Physical Review Letters 91, 253902 (2003).
  * Mayergoyz  _et al._ [2005] I. D. Mayergoyz, D. R. Fredkin, and Z. Zhang, Electrostatic (plasmon) resonances in nanoparticles, Physical Review B 72, 155412 (2005).
  * Gómez  _et al._ [2010] D. E. Gómez, K. C. Vernon, and T. J. Davis, Symmetry effects on the optical coupling between plasmonic nanoparticles with applications to hierarchical structures, Physical Review B 81, 075414 (2010).
  * Weiss  _et al._ [2011] T. Weiss, N. A. Gippius, S. G. Tikhodeev, G. Granet, and H. Giessen, Derivation of plasmonic resonances in the Fourier modal method with adaptive spatial resolution and matched coordinates, Journal of the Optical Society of America A 28, 238 (2011).
  * Bykov  _et al._ [2017] D. A. Bykov, E. A. Bezus, and L. L. Doskolovich, Use of aperiodic Fourier modal method for calculating complex-frequency eigenmodes of long-period photonic crystal slabs, Optics Express 25, 27298 (2017).
  * Bykov and Doskolovich [2013] D. A. Bykov and L. L. Doskolovich, Numerical methods for calculating poles of the scattering matrix with applications in grating theory, Journal of Lightwave Technology 31, 793 (2013).
  * Suryadharma  _et al._ [2017] R. N. S. Suryadharma, M. Fruhnert, I. Fernandez-Corbaton, and C. Rockstuhl, Studying plasmonic resonance modes of hierarchical self-assembled meta-atoms based on their transfer matrix, Physical Review B 96, 045406 (2017).
  * Guo  _et al._ [2012] H. Guo, B. Oswald, and P. Arbenz, 3-dimensional eigenmodal analysis of plasmonic nanostructures, Optics Express 20, 5481 (2012).
  * Paulus and Martin [2001] M. Paulus and O. J. F. Martin, Light propagation and scattering in stratified media: A Green’s tensor approach, Journal of the Optical Society of America A 18, 854 (2001).
  * Bernasconi  _et al._ [2016] G. D. Bernasconi, J. Butet, and O. J. Martin, Mode analysis of second-harmonic generation in plasmonic nanostructures, Journal of the Optical Society of America B 33, 768 (2016).
  * Powell [2017] D. A. Powell, Interference between the modes of an all-dielectric meta-atom, Physical Review Applied 7, 034006 (2017).
  * Bautista  _et al._ [2012] G. Bautista, M. J. Huttunen, J. Mäkitalo, J. M. Kontio, J. Simonen, and M. Kauranen, Second-harmonic generation imaging of metal nano-objects with cylindrical vector beams, Nano Letters 12, 3207 (2012).
  * Bautista  _et al._ [2015] G. Bautista, J. Mäkitalo, Y. Chen, V. Dhaka, M. Grasso, L. Karvonen, H. Jiang, M. J. Huttunen, T. Huhtio, H. Lipsanen, and M. Kauranen, Second-harmonic generation imaging of semiconductor nanowires with focused vector beams, Nano Letters 15, 1564 (2015).
  * Kant [2000] R. Kant, Superresolution and increased depth of focus: An inverse problem of vector diffraction, Journal of Modern Optics 47, 905 (2000).
  * Braat  _et al._ [2005] J. J. M. Braat, P. Dirksen, A. J. E. M. Janssen, S. van Haver, and A. S. van de Nes, Extended Nijboer–Zernike approach to aberration and birefringence retrieval in a high-numerical-aperture optical system, Journal of the Optical Society of America A 22, 2635 (2005).
  * Foreman  _et al._ [2008] M. R. Foreman, S. S. Sherif, P. R. T. Munro, and P. Török, Inversion of the Debye-Wolf diffraction integral using an eigenfunction representation of the electric fields in the focal region, Optics Express 16, 4901 (2008).
  * Jahn and Bokor [2013] K. Jahn and N. Bokor, Solving the inverse problem of high numerical aperture focusing using vector Slepian harmonics and vector Slepian multipole fields, Optics Communications 288, 13 (2013).
  * Chen  _et al._ [2016] Z. Chen, T. Zeng, and J. Ding, Reverse engineering approach to focus shaping, Optics Letters 41, 1929 (2016).
  * Rui  _et al._ [2016] G. Rui, J. Chen, X. Wang, B. Gu, Y. Cui, and Q. Zhan, Synthesis of focused beam with controllable arbitrary homogeneous polarization using engineered vectorial optical fields, Optics Express 24, 23667 (2016).
  * Yla-Oijala and Taskinen [2005] P. Yla-Oijala and M. Taskinen, Well-conditioned Muller formulation for electromagnetic scattering by dielectric objects, IEEE Transactions on Antennas and Propagation 53, 3316 (2005).
  * Stratton and Chu [1939] J. A. Stratton and L. J. Chu, Diffraction theory of electromagnetic waves, Physical Review 56, 99 (1939).
  * Schelkunoff [1936] S. A. Schelkunoff, Some equivalence theorems of electromagnetics and their application to radiation problems, The Bell System Technical Journal 15, 92 (1936).
  * Chew  _et al._ [2009] W. C. Chew, M. S. Tong, and B. Hu,  _Integral Equation Methods for Electromagnetic and Elastic Waves_ (Morgan & Claypool, 2009).
  * Orfanidis [2016] S. J. Orfanidis,  _Electromagnetic Waves and Antennas_ (https://www.ece.rutgers.edu/~orfanidi/ewa/, 2016).
  * Jin [2015] J.-M. Jin,  _Theory and Computation of Electromagnetic Fields_ , 2nd ed. (John Wiley & Sons, Inc., 2015).
  * Collin [1991] R. E. Collin,  _Field Theory of Guided Waves_ , 2nd ed. (John Wiley & Sons, Inc., 1991).
  * Kong [2008] J. A. Kong,  _Electromagnetic Wave Theory_ (EMW Publishing, 2008).
  * Colton and Kress [2013] D. Colton and R. Kress,  _Integral Equation Methods in Scattering Theory_ (Society for Industrial and Applied Mathematics, 2013).
  * Taskinen [2006] M. Taskinen,  _On the Implementation and Formulation of the Electromagnetic Surface Integral Equations_ , Ph.D. thesis, Helsinki University of Technology (2006).
  * Jackson [1998] J. D. Jackson,  _Classical Electrodynamics_ , 3rd ed. (Wiley, New York, 1998).
  * Wolf [1959a] E. Wolf, Electromagnetic diffraction in optical systems - I. An integral representation of the image field, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 253, 349 (1959a).
  * Richards and Wolf [1959] B. Richards and E. Wolf, Electromagnetic diffraction in optical systems, II. Structure of the image field in an aplanatic system, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 253, 358 (1959).
  * Pereira and Van de Nes [2004] S. F. Pereira and A. S. Van de Nes, Superresolution by means of polarisation, phase and amplitude pupil masks, Optics Communications 234, 119 (2004).
  * Leutenegger  _et al._ [2006] M. Leutenegger, R. Rao, R. Leitgeb, and T. Lasser, Fast focus field calculations, Optics Express 14, 11277 (2006).
  * Zhang  _et al._ [2018] G.-L. Zhang, X.-Z. Gao, Y. Pan, M.-D. Zhao, D. Wang, H.-H. Zhang, Y. Li, C. Tu, and H.-T. Wang, Inverse method to engineer uniform-intensity focal fields with arbitrary shape, Optics Express 26, 16782 (2018).
  * Wolf [1959b] E. Wolf, A scalar representation of electromagnetic fields: II, Proceedings of the Physical Society 74, 269 (1959b).
  * Snyder and Love [1983] A. W. Snyder and J. D. Love,  _Optical Waveguide Theory_ (Chapman & Hall, 1983).
  * Frigo and Johnson [2005] M. Frigo and S. G. Johnson, The design and implementation of FFTW3, Proceedings of the IEEE 93, 216 (2005).
  * Cheney and Kincaid [2007] E. W. Cheney and D. R. Kincaid,  _Numerical Mathematics and Computing_ , sixth ed. (Thomson Brooks/Cole, 2007).