# Time-Domain Topology Optimization of Power Dissipation in Dispersive
Dielectric and Plasmonic Nanostructures

Johannes Gedeon, Izzatjon Allayarov, Antonio Calà Lesina, and Emadeldeen
Hassan This work was supported by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) under Germany’s Excellence Strategy within the
Cluster of Excellence PhoenixD (EXC 2122, Project ID 390833453), and within
the Research Grant CA 2763/2-1 (Project ID 527470210). The computing time was
granted by the Resource Allocation Board and provided on the supercomputer
Lise and Emmy/Grete at NHR@ZIB and NHR@Göttingen as part of the NHR
infrastructure. The calculations for this research were conducted with
computing resources under the project nip00059. J.G., A.C.L. and E.H.
formulated the research question, J.G. developed the theory, implemented the
code, performed the study and wrote the manuscript. I.A. provided the Mie
analysis section. A.C.L and E.H. supervised the work and revised the
manuscript.Johannes Gedeon, Izzatjon Allayarov and Antonio Calà Lesina are
with the Hannover Centre for Optical Technologies, the Institute for Transport
and Automation Technology (Faculty of Mechanical Engineering), and the Cluster
of Excellence PhoenixD, Leibniz University Hannover, 30167 Hannover, Germany
(e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>hannover.de).Emadeldeen Hassan is with the Department of Applied Physics and
Electronics, Umeå University, SE-901 87 Umeå, Sweden.

###### Abstract

We present a density-based topology optimization scheme for locally optimizing
the electric power dissipation in nanostructures made of lossy dispersive
materials. By using the complex-conjugate pole-residue (CCPR) model, we can
accurately model any linear materials’ dispersion without limiting to specific
material classes. We incorporate the CCPR model via auxiliary differential
equations (ADE) into Maxwell’s equations in the time domain, and formulate a
gradient-based topology optimization problem to optimize the dissipation over
a broad spectrum of frequencies. To estimate the objective function gradient,
we use the adjoint field method, and explain the discretization and
integration of the adjoint system into the finite-difference time-domain
(FDTD) framework. Our method is demonstrated using the example of topology
optimized spherical nanoparticles made of Gold and Silicon with an enhanced
absorption efficiency in the visible-ultraviolet spectral range. In this
context, a detailed analysis of the challenges of topology optimization of
plasmonic materials associated with a density-based approach is given.

###### Index Terms:

absorption efficiency, adjoint method, complex-conjugate pole–residue pairs
model, FDTD method, Gold, instantaneous electric power dissipation, inverse
design, optical dispersion, plasmonics, Silicon, time domain, topology
optimization

##  I Introduction

Topology optimization (TopOpt) in nanophotonics has proven to be a powerful
method to generate novel nanostructures with desired optical functionalities.
Its inverse design approach enables the development of optical components that
far surpass those created through human intuition and physical first
principles in terms of efficiency. Topology optimization commonly refers to a
density-based approach and the use of the adjoint method to enable a gradient-
based optimization[1, 2]. The density is mapped to the material’s distribution
in space, and is iteratively updated to maximize a given objective function.
By using the adjoint method, the gradients are calculated with only two
simulations. Over the past decade, it has been applied to various engineering
challenges, from developing achromatic metalenses [3] and broadband
(subwavelength-) antennas [4, 5, 6, 7, 8], to single-photon emitters for
quantum computing [9], nonlinear multiplexer[10], multilayer waveguides
transitions [11], and small-scale particle accelerators [12].

Most of these optimization objectives are directly or indirectly linked to the
electric energy, and typically operate in a frequency regime in which the
optical material dispersion and losses are negligible. However, towards higher
frequencies and with the requirement of broadband performance, the assumption
of a dispersionless material cannot be made. Dissipation (or absorption), like
energy itself, is a fundamental quantity of the Poynting theorem and
intrinsically linked to dispersion by thermodynamics laws [13]. Absorption is
commonly optimized implicitly by tackling reflection and transmission in the
far-field [14, 15], or based on the physical properties of specific
combinations of materials, such as metal-dielectric-metal and multi-layered
architectures [16, 17, 18, 19], photonic crystals [20, 21, 22], or
incorporation of plasmonic nanoparticles into the absorbing material [23, 24].
It is timely to conduct an in-depth examination of this physical quantity in
the context of topology optimization for dispersive nanostructures.

In this work, we present a topology optimization method for locally optimizing
the electric power dissipation for arbitrary dispersive materials. Contrary to
conventional methods, our approach is explicit because the design region
overlaps/coincides with the region where the optimization objective is
specified. This opens the door to solving novel design problems where the
dissipation (or conversion of energy) in the structure itself is the objective
of interest, with possible implications in solar energy conversion in
(thermo-) photovoltaics [25], thermo-nanophotonics [26, 27], radiotherapy and
thermal ablation [28], photodetection [29, 30, 31], absorptive filters and
sensors [32, 33, 34], or plasmonic devices with a reduced power loss [35].

We use the adjoint method, relying on the formulation of Maxwell’s equations
in the time domain. Time-domain topology optimization increasingly attracts
attention as it is well suited to tackle a variety of special optimization
problems, such as the optimization of dynamic and transient effects, pulse
shaping, and time-varying materials [7, 36, 37, 38]. Since a time-domain
approach allows the entire spectral content of interest to be incorporated
into the time-dependent source of excitation, the method allows us to target
an optimized performance over a broad continuous frequency spectrum.
Furthermore, dispersive materials can be simulated with the finite-difference
time-domain (FDTD) method, e.g., via the auxiliary differential equation (ADE)
approach [39, 40, 41, 42], and can therefore be optimized across a broad
spectrum. We recently demonstrated such an approach by optimizing the field
confinement in dielectric and plasmonic nanostructures based on the complex-
conjugate pole-residue (CCPR) model [7]. The CCPR model is not restricted to a
certain physical model, as it allows the fitting of any dielectric function
[43, 42]. We build on these results and integrate this model in our equations
to enable the maximization/minimization of the dissipation for arbitrary
materials.

Starting from the Poynting theorem, we find the expression of the
instantaneous electric power dissipation of lossy dispersive media described
by the CCPR model (Sec. II). The time-average of this expression serves as our
objective function, and we present the corresponding adjoint scheme in Sec.
III-A to enable a gradient-based optimization. The discretization of the
adjoint system using the FDTD method is presented in Sec. III-C. A brief
overview of the common techniques in density-based topology optimization is
given in Sec. III-B. In this context, we provide a heuristic procedure to
determine an artificial damping term, which has a significant impact on the
optimization convergence when plasmonic materials are considered. As a study
case and to demonstrate our method, we chose the topology optimization of
spherical Silicon and Gold nanoparticles with an enhanced absorption
efficiency in the visible-ultraviolet regime. The results are presented in
Sec. IV, including an assessment of the success and challenges associated with
our optimization method.

##  II Instantaneous electric power dissipation density in lossy dispersive
media

In this section, we present the expression for the instantaneous electric
power dissipation density based on the CCPR model. The time-averaged
dissipation is defined later as our objective function. Given that, we derive
the time-domain adjoint scheme to maximize/minimize this quantity within the
topology optimization framework in Sec. III-A. In the following, we consider
an isotropic and non-magnetic medium whose dielectric function can be
described by CCPR poles with $e^{\,j\omega t}$ time-dependency as

$\varepsilon(\omega)=\varepsilon_{\infty}+\frac{\sigma}{j\omega\varepsilon_{0}}+\sum_{p=1}^{P}\left(\frac{c_{p}}{j\omega-
a_{p}}+\frac{c_{p}^{*}}{j\omega-a_{p}^{*}}\right).$ (1)

We note that the coefficients $a_{p}$ and $c_{p}$ themselves can be complex,
and $*$ represents the complex conjugation. The model can accurately describe
the complex permittivity of any lossy dispersive material as long as the
number of poles $P$ is chosen appropriately. By the proper selection of the
coefficients, this model incorporates all the standard and advanced dispersion
models commonly used, such as Debye, Drude (+ critical points), and (modified)
Lorentz. [42]. Using the auxiliary differential equation method, the full
system of the time-domain Maxwell’s equations in a source-free region reads

$\displaystyle-\nabla\times\mathbf{H}+\varepsilon_{0}\varepsilon_{\infty}\partial_{t}\mathbf{E}+\sigma\mathbf{E}+2\sum_{p=1}^{P}\Re\left\\{\partial_{t}\mathbf{Q}_{p}\right\\}$
$\displaystyle=0,$ (2a) $\displaystyle\forall
p\in\\{1,\ldots,P\\}:\partial_{t}\mathbf{Q}_{p}-a_{p}\mathbf{Q}_{p}-\varepsilon_{0}c_{p}\mathbf{E}$
$\displaystyle=0,$ (2b)
$\displaystyle\mu_{0}\partial_{t}\mathbf{H}+\nabla\times\mathbf{E}$
$\displaystyle=0,$ (2c)

where the complex auxiliary fields $\mathbf{Q}_{p}$ must be computed for all
poles $p\in\\{1,...,P\\}$. To derive the expression for the instantaneous
electric power dissipation based on the CCPR model, we follow a similar
procedure as the one presented in Ref. [44], which is therein limited to
Lorentz media only. We start from the general Poynting theorem in a source-
free domain $\Omega_{\mathrm{o}}$[45],

$-\oint_{\partial\Omega_{\mathrm{o}}}(\mathbf{E}\times\mathbf{H})\cdot\mathrm{d}\mathbf{a}=\int_{\Omega_{\mathrm{o}}}\left(\mathbf{E}\cdot\frac{\partial\mathbf{D}}{\partial
t}+\mathbf{H}\cdot\frac{\partial\mathbf{B}}{\partial
t}\right)\mathrm{d}^{3}r,$ (3)

with the ansatz

$\mathbf{E}\cdot\frac{\partial\mathbf{D}}{\partial t}=\frac{\partial
u_{e}}{\partial t}+q_{e}.$ (4)

Here, $u_{e}(t)$ represents the instantaneous electric energy density, and
$q_{e}(t)$ is the instantaneous electric dissipation density, respectively.
Using the relation $\partial_{t}\mathbf{D}=\nabla\times\mathbf{H}$, we can
calculate the product in Eq. (4) and separate the terms according to the
ansatz on the right-hand side (Supplementary Material, Sec. II). As a result,
the electric power dissipation density can be identified as

$q_{e}(t)=\sigma\mathbf{E}^{2}+2\sum_{p=1}^{P}\Re\left\\{\frac{(\partial_{t}\mathbf{Q}_{p})^{2}}{\varepsilon_{0}c_{p}}\right\\}.$
(5)

We note that this equation is valid for any time-varying electric field
$\mathbf{E}(t)$. To evaluate the dissipation over an arbitrary spectral range
of frequencies in the steady state (Appendix A), we now assume a time harmonic
field $\mathbf{E}(t):=\Re\\{\mathbf{\hat{E}}_{0}e^{j\omega t}\\}$, where the
amplitude $\mathbf{\hat{E}}_{0}$ itself is frequency independent. In that
case, the auxiliary fields from Eq. (2b) take the form

$\mathbf{Q}_{p}(t)=\,\frac{1}{2}\left(\frac{\varepsilon_{0}c_{p}}{j\omega-
a_{p}}\hat{\mathbf{E}}_{0}e^{j\omega t}+\frac{\varepsilon_{0}c_{p}}{-j\omega-
a_{p}}\hat{\mathbf{E}}_{0}^{*}e^{-j\omega t}\right).$ (6)

Substituting this expression into the Eq. (5) and averaging with respect to
time, the product terms containing factors $e^{\pm 2j\omega t}$ give zero [44,
13]. This leaves the equation of the time-averaged electric power dissipation,

$\bar{q}_{e}(\omega)=\left(\frac{1}{2}\sigma+\sum_{p=1}^{P}\Re\left\\{\frac{\omega^{2}\varepsilon_{0}c_{p}}{(j\omega-
a_{p})(-j\omega-a_{p})}\right\\}\right)|\hat{\mathbf{E}}_{0}|^{2}.$ (7)

This expression can also be referred to as the dc component of the
instantaneous electric power dissipation. The result is in agreement with the
model-independent expression derived by Landau and Lifschitz [13],

$\bar{q}_{e}=\frac{1}{2}\varepsilon_{0}\omega\varepsilon^{\prime\prime}|\mathbf{\hat{E}}_{0}|^{2},$
(8)

considering the convention
$\varepsilon=\varepsilon^{\prime}-j\varepsilon^{\prime\prime}$ of the
permittivity.

##  III Methods

###  III-A Adjoint Scheme

Here, we present the time-dependent adjoint scheme to optimize the power
dissipation in dispersive media within the density-based topology optimization
framework. The formulation builds on the method proposed in Ref. [7]. A
comprehensive derivation is provided in the Supplementary Material, Sec. III.
Density-based topology optimization (TopOpt) for inverse design refers to an
iterative design process that allows us to optimize the distribution of a
given material in a specified domain in order to optimize a certain objective
function [1]. The method requires the description of the material in the
design domain as a spatial density distribution $0\leq\rho(\bm{r})\leq 1$
which is mapped to the material’s properties, such as the permittivity given
in Eq. (1) in our case, and consequently describes the topological shape of
the photonic device as the density converges to a binary design.

In the following, we consider a design region $\Omega$ and define the
objective function on a subset $\Omega_{\mathrm{o}}\subseteq\Omega$ as

$F:=\frac{1}{T}\int_{\Omega_{\mathrm{o}}\times
I}q_{e}(t,\rho(\bm{r}))\mathrm{d}t\mathrm{d}^{3}r.$ (9)

This equation represents the time-averaged electric power loss in a subregion
$\Omega_{\mathrm{o}}$ lying in a material-interpolated design over a time
interval $I=[0,T]$. $T$ denotes the time period until the fields are decayed
after the excitation of a light pulse, which is injected outside the domain
$\Omega$. If $\Omega\equiv\Omega_{\mathrm{o}}$ and the design is considered to
be entirely surrounded by a non-absorbing medium, the objective is equivalent
to the average power transferred from the light pulse to the nanostructure.
Fig. 1 shows a schematic representation of the sets $\Omega$ and
$\Omega_{\mathrm{o}}$, which we refer to when deriving the gradients of the
objective with respect to the density on subsets.

Figure 1: Schematic illustration of the sets $\Omega$ and
$\Omega_{\mathrm{o}}$. The objective function $F$ is defined on a subset
$\Omega_{\mathrm{o}}$ (observation region) lying in the design region
$\Omega$. The design material can vary in both regions during the
optimization.

Following a linear interpolation scheme of the parameters
$\varepsilon_{\infty}$, $\sigma$ and the pole terms of the CCPR model, Eq.
(1), between background and design material, the density-dependent dissipation
term reads

$q_{e}(t,\rho)=\sigma(\rho)\mathbf{E}^{2}+2\sum_{i=1}^{2}\kappa^{(i)}(\rho)\sum_{p=1}^{P}\Re\left\\{\frac{(\partial_{t}\mathbf{Q}_{p}^{(i)})^{2}}{\varepsilon_{0}c_{p}^{(i)}}\right\\},$
(10)

where $\kappa^{(1)}(\rho):=1-\rho$ and $\kappa^{(2)}(\rho):=\rho$. The indices
$i=1,2$ represent the background and design material, respectively. We expand
the linear mixed parameter $\sigma(\rho)$ by an additional term, such that

$\sigma(\rho):=\sum_{i\in\\{1,2\\}}\kappa^{(i)}(\rho)\sigma^{(i)}+\rho(1-\rho)\gamma.$
(11)

Here, $\gamma$ is an artificial damping parameter. As we will see, this
parameter plays a crucial role for the convergence based on the objective we
have chosen, as it is directly incorporated into Eq. (10) and its gradients
with respect to the density. We describe our choices and its impact on the
example of Gold in detail in Sec. III-B2.
Defining two material-dependent fields $\mathbf{f}^{(i)}$, $i=1,2$, as

$\mathbf{f}^{(i)}:=\varepsilon_{0}\varepsilon_{\infty}^{(i)}\partial_{t}\mathbf{E}+\sigma^{(i)}\mathbf{E}+2\sum_{p=1}^{P^{(i)}}\Re\left\\{\partial_{t}\mathbf{Q}_{p}^{(i)}\right\\},$
(12)

the Maxwell’s Eqs. (2) are changed according to their density dependency in
$\Omega$ to

$\displaystyle-\nabla\times\mathbf{H}+\sum_{i=1}^{2}\kappa^{(i)}(\rho)\,\mathbf{f}^{(i)}+\rho(1-\rho)\gamma\,\mathbf{E}$
$\displaystyle=0,$ (13a) $\displaystyle\begin{subarray}{c}\text{For
}i=1,2\text{ and }\\\ \forall
p\,\in\\{1,\ldots,P^{(i)}\\}\end{subarray}\quad\partial_{t}\mathbf{Q}_{p}^{(i)}-a_{p}^{(i)}\mathbf{Q}_{p}^{(i)}-\varepsilon_{0}c_{p}^{(i)}\mathbf{E}$
$\displaystyle=0,$ (13b)
$\displaystyle\mu_{0}\partial_{t}\mathbf{H}+\nabla\times\mathbf{E}$
$\displaystyle=0.$ (13c)

The auxiliary fields $\mathbf{Q}_{p}^{(i)}$ must be computed for both
materials $i=1,2$ and each corresponding pole $p\in\\{1,...,P^{(i)}\\}$. We
refer to this system of equations in the following as forward simulation, in
which the objective function is measured.

To update the design using a gradient-based optimization method, we need the
gradient information of the objective function with respect to the density.
Therefore, we perform an additional adjoint simulation, which differs from the
forward one only in terms of excitation of the physical system. Based on the
interpolation scheme and the definition of our objective, we introduce adjoint
sources that are injected into the observation region $\Omega_{\mathrm{o}}$.
They consist of an electric source term,

$\mathbf{S}_{E}:=\phantom{-}2\,T^{-1}\sigma(\rho)\overleftarrow{\mathbf{E}},$
(14)

and auxiliary adjoint source terms,

$\mathbf{S}_{\partial_{\tau}Q_{p}}^{(i)}:=\phantom{-}2\,T^{-1}\partial_{\tau}\overleftarrow{\mathbf{Q}}_{p}^{(i)},$
(15)

for $i=1,2$ and $\forall p\in\\{1,\ldots,P^{(i)}\\}$. Here, $\tau=T-t$ denotes
the transformed variable after time reversal, and the symbol
”$\overleftarrow{}$” marks the time reversal of the corresponding fields
calculated during the forward simulation. We again define two material-
dependent fields $\tilde{\mathbf{f}}^{(i)}$, $i=1,2$ as

$\tilde{\mathbf{f}}^{(i)}:=\varepsilon_{0}\varepsilon_{\infty}^{(i)}\partial_{\tau}\tilde{\mathbf{E}}+\sigma^{(i)}\tilde{\mathbf{E}}+2\sum_{p=1}^{P^{(i)}}\Re\left\\{\partial_{\tau}\tilde{\mathbf{Q}}_{p}^{(i)}\right\\}.$
(16)

Then, the adjoint system on $\Omega_{\mathrm{o}}$ reads

$\displaystyle-\nabla\times\tilde{\mathbf{H}}+\sum_{i=1}^{2}\kappa^{(i)}(\rho)\,\tilde{\mathbf{f}}^{(i)}+\rho(1-\rho)\gamma\,\tilde{\mathbf{E}}$
$\displaystyle=\mathbf{S}_{E},$ (17a)
$\displaystyle\begin{subarray}{c}\text{For }i=1,2\text{ and }\\\ \forall
p\,\in\\{1,\ldots,P^{(i)}\\}\end{subarray}\;\;\partial_{\tau}\tilde{\mathbf{Q}}_{p}^{(i)}-a_{p}^{(i)}\tilde{\mathbf{Q}}_{p}^{(i)}-\varepsilon_{0}c_{p}^{(i)}\tilde{\mathbf{E}}$
$\displaystyle=\mathbf{S}_{\partial_{\tau}Q_{p}}^{(i)},$ (17b)
$\displaystyle\mu_{0}\partial_{\tau}\tilde{\mathbf{H}}+\nabla\times\tilde{\mathbf{E}}$
$\displaystyle=0.$ (17c)

Since the objective $F$ in Eq. (9) is defined on the subset
$\Omega_{\mathrm{o}}$ only, no adjoint sources will be injected into the
domain $\Omega\setminus\Omega_{\mathrm{o}}$. In that case, the right-hand side
of the system of equations (17) is 0. The adjoint fields
$\\{\tilde{\mathbf{E}}$, $\tilde{\mathbf{Q}}_{p}^{(i)}\\}$ together with the
forward fields $\\{\mathbf{E}$, $\mathbf{Q}_{p}^{(i)}\\}$ are used to
calculate the gradient of the objective function $F$ with respect to the
density,

$\nabla_{\rho}F:=\mathcal{A}_{|\rho\in\Omega}+\mathcal{B}_{|\rho\in\Omega_{\mathrm{o}}},$
(18)

which is a sum of two gradient terms defined on each set $\Omega$ and
$\Omega_{\mathrm{o}}$ illustrated in Fig. 1, namely

$\displaystyle\mathcal{A}_{|\rho\in\Omega}:=$
$\displaystyle\;-\int_{I}\varepsilon_{0}(\mathrm{d}_{\rho}\varepsilon_{\infty})\partial_{\tau}\tilde{\mathbf{E}}\cdot\overleftarrow{\mathbf{E}}$
(19a)
$\displaystyle-\int_{I}(\mathrm{d}_{\rho}\sigma)\tilde{\mathbf{E}}\cdot\overleftarrow{\mathbf{E}}$
$\displaystyle-\int_{I}\sum_{i=1}^{2}\sum_{p=1}^{P^{(i)}}2(\mathrm{d}_{\rho}\kappa^{(i)})\partial_{\tau}\tilde{\mathbf{E}}\cdot\Re\left\\{\overleftarrow{\mathbf{Q}}_{p}^{(i)}\right\\},$
$\displaystyle\mathcal{B}_{|\rho\in\Omega_{\mathrm{o}}}:=$
$\displaystyle\;\;\;\;\,\frac{1}{T}\int_{I}(\mathrm{d}_{\rho}\sigma)\overleftarrow{\mathbf{E}}^{2}$
(19b)
$\displaystyle+\frac{1}{T}\int_{I}\sum_{i=1}^{2}\sum_{p=1}^{P^{(i)}}2(\mathrm{d}_{\rho}\kappa^{(i)})\Re\left\\{\frac{(\partial_{\tau}\overleftarrow{\mathbf{Q}}_{p}^{(i)})^{2}}{\varepsilon_{0}c^{(i)}_{p}}\right\\}.$

From these equations, we observe simplifications for two important special
cases:

  * •

If both background and design material in $\Omega$ do not contain any CCPR
poles, i.e., they are both non-dispersive and their losses are only described
by the $\sigma$ term in Eq. (1), all auxiliary fields vanish. And since
$\mathbf{S}_{\partial_{\tau}Q_{p}}^{(i)}\equiv 0$, no auxiliary adjoint source
will be injected during the adjoint simulation. In addition, the gradients in
Eqs. (19) do not contain any terms related to the $\mathbf{Q}_{p}^{(i)}$
fields.

  * •

If the objective function in Eq. (9) defined on $\Omega_{\mathrm{o}}$ does not
depend on the density itself, i.e., the density does not vary in
$\Omega_{\mathrm{o}}$ during the optimization, the derivatives of the material
parameters with respect to the density vanish in $\Omega_{\mathrm{o}}$. And
since $\forall\rho\in\Omega_{\mathrm{o}}$:
$\mathrm{d}_{\rho}\varepsilon_{\infty}=\mathrm{d}_{\rho}\sigma=\mathrm{d}_{\rho}\kappa^{(i)}=0$,
it follows that $\mathcal{B}_{|\rho\in\Omega_{\mathrm{o}}}\equiv 0$, i.e., the
total gradient reduces to
$\nabla_{\rho}F:=\mathcal{A}_{|\rho\in\Omega\setminus\Omega_{o}}$, in
agreement with the results presented in Refs. [7, 11].

###  III-B Density-based topology optimization techniques

####  III-B1 Filtering and Projection

In density-based topology optimization, it is common to apply a filter and
projection on the density $\rho$. It can serve the purpose of eradicating the
appearance of single-pixel features, introducing a weak sense of length-scale
into the design, or curing the self-penalization issue when optimizing lossy
structures [46, 47, 48, 49]. The transformed form $\bar{\tilde{\rho}}$
describes the material interpolation and is also incorporated into the
objective function Eq. (9) in our case.
At each iteration, we average over a neighborhood around the density points
$\rho(\bm{r})$ by applying the filter operator

$\tilde{\rho}(\bm{r})=(\mathscr{F}\rho)(\bm{r}):=\frac{\int_{\mathcal{B}_{R}(\bm{r})}w\left(\bm{r},\bm{r}^{\prime}\right)\rho(\bm{r}^{\prime})\,\mathrm{d}^{3}r^{\prime}}{\int_{\mathcal{B}_{R}(\bm{r})}w\left(\bm{r},\bm{r}^{\prime}\right)\,\mathrm{d}^{3}r^{\prime}}.$
(20)

$\mathcal{B}_{R}(\bm{r})$ describes a sphere with radius $R$ around $\bm{r}$,
and $w\left(\bm{r},\bm{r}^{\prime}\right)$ is a weighting function defined as

$w\left(\bm{r},\bm{r}^{\prime}\right)=R-\left|\bm{r}-\bm{r}^{\prime}\right|.$
(21)

Subsequently, we project the filtered density by the smoothed Heaviside
function

$\bar{\tilde{\rho}}(\bm{r})=\mathscr{P}(\tilde{\rho}(\bm{r})):=\frac{\tanh(\beta\eta)+\tanh\left(\beta\left(\tilde{\rho}(\bm{r})-\eta\right)\right)}{\tanh(\beta\eta)+\tanh(\beta(1-\eta))}.$
(22)

The parameter $\eta$ determines the threshold, and the parameter $\beta$
controls the sharpness of the projection. For $\beta\rightarrow\infty$, the
projected density is binary, i.e., it has only the values 0 or 1. The
derivative of a functional $F[\rho]$ with respect to the original density
$\rho$ can be calculated by using the chain rule together with the Fréchet
derivatives of the filtering and projection operators,

$\frac{\delta F[\rho]}{\delta\rho(\bm{r})}=\frac{\delta
F[\bar{\tilde{\rho}}]}{\delta\bar{\tilde{\rho}}(\bm{r})}\mathscr{P}^{\prime}(\tilde{\rho}(\bm{r}))(\mathscr{F}^{\prime}\rho)(\bm{r}).$
(23)

The parameter $\beta$ value will be increased during the optimization until
the objective does not show any significant change. The measure of non-
discreteness $M_{\text{nd}}$ serves as an indicator to verify the convergence
of an optimized design to a binary solution [50],

$M_{\text{nd}}:=100\%\times
V_{\Omega}^{-1}\int_{\Omega}4\,\bar{\tilde{\rho}}(\bm{r})\left(1-\bar{\tilde{\rho}}(\bm{r})\right)\mathrm{d}^{3}r.$
(24)

This measure has a maximum if all density points have an intermediate value of
0.5, and it is minimized if the density only consists of the values 0 or 1. At
the end of the optimization, the projected density is mapped to a binary
design by a threshold with respect to the parameter $\eta$.

####  III-B2 The artificial damping term $\gamma$ for Gold

The static conductivity $\sigma$ appearing in the permittivity model Eq. (1),
has been interpolated as

$\sigma(\rho)=\sum_{i\in\\{1,2\\}}\kappa^{(i)}(\rho)\sigma^{(i)}+\rho(1-\rho)\gamma.$
(25)

The first term represents the linear interpolation between background and
design material. The second term related to $\gamma$ was originally introduced
to penalize intermediate densities of dielectric materials [51]. However, this
artificial damping term has proven to be essential for the optimization of
plasmonic materials as well [8, 7, 52]. The reference [8] in particular,
highlights the significance of a proper choice of the artificial damping
parameter $\gamma$, which has therein been heuristically observed as part of
optimizing Silver nanoantennas for field confinement in a non-dispersive
observation region.

We emphasize that the choice of the value of $\gamma$ is also essential when
maximizing the dissipation of a spherical Gold nanostructure in the visible
spectral range, as will be investigated in Sec. IV. For wavelengths above
$\approx 520\,\mathrm{nm}$, Gold nanoparticles typically show localized
surface plasmon resonances. Below this wavelength, however, the interband
transition dominates, where the Gold’s electrons absorb photons and transition
from the filled $d$-band to states above the Fermi level in the $sp$-band,
which will effectively results in a higher absorption of light. In the context
of density-based topology optimization, we must carefully choose the
artificial damping term $\gamma$ to ensure a reasonable behavior of the
creation and propagation of plasmons for intermediate density values, and to
guarantee the convergence to a binary design as a result. We will therefore
only focus on the plasmonic region in the following discussion. Although we
present a heuristic method specifically for Gold here, the underlying approach
is applicable to plasmonic materials in general. The unit for the artificial
damping term $\gamma$ is the same as that for the static conductivity
$\sigma$, which we will omit for the sake of better readability here.

We can estimate the influence of the damping term on maximizing the
dissipation of a spherical Gold structure by analyzing the propagation of
surface plasmon polaritons (SPP) along the interface between air and Gold in
1D. In that case, the prerequisites for the formation of surface plasmons are
well studied [53]. The assumption of an infinite flat layer, rather than a
complexly shaped surface, represents a drastic simplification. However, it is
reasonable as a first approximation if we interpret resonant localized surface
plasmons of a sphere as the propagation of two counterpropagating surface
plasmon waves [54], which later also serves as our initial design. We consider
a surface plasmon propagating along $x^{+}$ at the interface of air ($z>0$)
and Gold ($z<0$). Inside Gold, its magnetic and electric field components are
described as

Figure 2: Real part $\Re\\{k_{\perp}\\}$ (1st column) and negative imaginary
part $-\Im\\{\beta_{\parallel}\\}$ (2nd column) from Eqs. (27) for a material-
interpolated permittivity $\varepsilon(\bar{\rho})$ of Gold and air. The
density $\rho$ was projected according to Eq. (22) with $\eta=0.55$ and
$\beta=7$. Each row corresponds to a different artificial damping term
$\gamma$ used in the mixing of the static conductivity in Eq. (25). The
colored area marks the frequency range in which Gold mainly shows plasmonic
behavior. The damping term $\gamma=2\times 10^{5}$ has been chosen for the
optimization of the Gold nanostructure in Sec. IV.

Figure 3: Propagating surface plasmons along the interface of air and the
interpolated material $\varepsilon(\bar{\rho})$ for different artificial
damping parameters $\gamma$. All plasmons have the same frequency of
$\omega=2.8\times 10^{15}$ rad/s, and the density $\rho=0.4$ has been used for
the interpolation of the permittivities of air and Gold (with a preceding
projection), see Fig. 3. The dashed and solid lines represent the magnitude
and real part of $E_{z}$ from Eq. (26c) at a depth of $\Delta z=22$ nm within
the material, respectively.

$\displaystyle
H_{y}(x,z)=A\mathrm{e}^{-j\beta_{\parallel}x}\mathrm{e}^{k_{\perp}z},$ (26a)
$\displaystyle
E_{x}(x,z)=jA\frac{1}{\omega\varepsilon_{0}\varepsilon}k_{\perp}e^{-j\beta_{\parallel}x}\mathrm{e}^{k_{\perp}z},$
(26b) $\displaystyle
E_{z}(x,z)=-A\frac{\beta_{\parallel}}{\omega\varepsilon_{0}\varepsilon}e^{-j\beta_{\parallel}x}\mathrm{e}^{k_{\perp}z}.$
(26c)

$A$ is an arbitrary amplitude, $\varepsilon$ is the relative permittivity of
Gold, and $k_{\perp}$ and $\beta_{\parallel}$ represent the attenuation and
the propagation constants of the surface plasmon, respectively,

$\displaystyle k_{\perp}=$
$\displaystyle\;\sqrt{\beta_{\parallel}^{2}-k_{0}^{2}\varepsilon},$ (27a)
$\displaystyle\beta_{\parallel}=$
$\displaystyle\;k_{0}\sqrt{\frac{\varepsilon}{1+\varepsilon}}.$ (27b)

Now, we replace the permittivity $\varepsilon$ of Gold in Eqs. (26) & (27) by
the material-interpolated permittivity between air and Gold
$\varepsilon(\bar{\rho})$. Here, the density $\bar{\rho}=\mathscr{P}(\rho)$ is
projected according to Eq. (22) with the parameters $\eta=0.55$ and $\beta=7$.
The plot of the values $\Re\\{k_{\perp}\\}$ and $-\Im\\{\beta_{\parallel}\\}$
for different density values $\rho$ and damping parameters $\gamma$ is
presented in Fig. 3. $\Re\\{k_{\perp}\\}$ describes the wave damping inside
Gold, while $-\Im\\{\beta_{\parallel}\\}$ determines the attenuation along the
propagation direction. Fig. 3 schematically illustrates the propagation of
SPPs having the same frequency for different artificial damping parameters,
just below the interface between air and the interpolated material. For this
discussion, we included the projected density $\bar{\rho}$ in the material
interpolation, since it is also incorporated into Maxwell’s equations in due
course during the optimization. A finer tuning of the equidistance between the
$\Re\\{k_{\perp}\\}$ functions for intermediate densities can be achieved by a
suitable selection of the projection parameters and the artificial damping
parameter. The benefits of this will become clear when interpreting the
diagrams in the following.

For $\gamma=0$ (Fig. 3, 1st row), the $-\Im\\{\beta_{\parallel}\\}$ functions
are characterized by large values and sharp peaks for intermediate densities
in the plasmonic region (marked in color). At these frequencies, the SPPs are
strongly confined locally around their excitation point. By considering the
corresponding $\Re\\{k_{\perp}\\}$ values in addition, we can conclude that
they are characterized by a deep penetration into the material, which will
effectively contribute to the dissipation (Fig. 3, red). The design might
retain gray areas throughout the entire optimization process, which
contributes to high dissipation around the excitation point, but prevents the
plasmons from exploring the topology in its full spatial extent. This can
result in an unpredictable performance after the final threshold, where the
localized effects for intermediate densities that contributed to the design
evolvement are eventually erased.

Choosing a very large damping term $\gamma=1\times 10^{8}$ instead (Fig. 3,
3rd row), eradicates the large and sharp peaks of
$-\Im\\{\beta_{\parallel}\\}$, which allows an extended spatial propagation
along the surface of the SPPs, as their decay is weak. However, the
penetration into the material is drastically reduced, which can be observed
from the large values of $\Re\\{k_{\perp}\\}$ for the majority of intermediate
densities. In this case, the interpolated material tends to have the
characteristics of a perfect metal, which makes optimization less sensitive,
as a dissipation inside the interpolated material is prevented (Fig. 3, blue).

The value $\gamma=2\times 10^{5}$ (Fig. 3, 2nd row) is a compromise. The large
peaks of the $-\Im\\{\beta_{\parallel}\\}$-functions are flattened, such that
the SPPs can cover a greater distance along the propagation direction for
intermediate densities (Fig. 3, purple). Furthermore, a more uniform
equidistant pattern of the intermediate $\Re\\{k_{\perp}\\}$-functions is
given, which is bounded by the curves of pure air ($\rho=0$) and Gold
($\rho=1$). We expect this to result in more stable convergence behaviour
during optimization. It is not guaranteed that the design will not retain gray
areas, but the dominant influence of highly confined plasmons related to an
intermediate density will be significantly reduced compared to the choice
$\gamma=0$. Therefore, we can expect an improvement in the convergence towards
a binary structure. Using an unprojected density $\rho$ for the material
interpolation yield a dense clustering close to the curve with $\rho=1$ for
the majority of intermediate functions here. By its projection with the
parameters $\eta=0.55$ and $\beta=7$, as shown in the figure, a more uniform
interpolation was achieved. We chose these values as our initial projection
parameters and the material parameter $\gamma=2\times 10^{5}$ when optimizing
the spherical nanoparticle made of Gold in Sec. IV.

###  III-C FDTD Implementation

Figure 4: Schematic illustration of the discretized forward and adjoint fields
sampled on the timeline within the FDTD framework. The adjoint fields (yellow)
are shifted in time by $\Delta t/2$ compared to the forward fields (blue).

In the FDTD framework, Maxwell’s equations for both forward and adjoint
simulation must be discretized on a staggered “Yee grid” in space and time
[55]. The discretization scheme we used is presented in Appendix B, including
the update equations of electric and auxiliary fields. The extension to
anisotropic media can be found in the Supplementary Material, Sec. III 1-2.
Hereafter, we will limit our discussion to merely explaining time
discretization, and present the resulting discretized expression of the
gradients in Eqs. (17).
We consider a time interval $[0,T]$, divided into $M+1$ time steps
$t_{m}:=\frac{m}{M}T=m\Delta t$, $m\in\\{0,...,M\\}$. The fields for the
forward simulation have been discretized in the same way as presented in Ref.
[42]. Therein, Eqs. (2a) and (2b) are discretized at time steps $(m+1/2)\Delta
t$, and Eq. (2c) at time steps $m\Delta t$. Consequently, the electric and
auxiliary fields are discretized at time steps $m\Delta t$, and the magnetic
field at $(m+1/2)\Delta t$. To avoid a different update scheme for the adjoint
system than for the forward system, we discretize both Eqs. (17a) and (17b) at
time steps $m\Delta t$, and Eq. (17c) at time steps $(m+1/2)\Delta t$. That
means, the sampled adjoint equations are shifted in time by $\Delta t/2$
compared to the forward equations. As a consequence, both adjoint electric and
auxiliary fields must be sampled at time steps $(m+1/2)\Delta t$, and the
adjoint magnetic fields must be sampled at time steps $m\Delta t$. A schematic
illustration of the disposition of the discretized fields on the timeline is
demonstrated in Fig. 4.
This offset in time, compared to the forward fields, offers two further
advantages in terms of numerical accuracy and stability: (1) The injection of
the source term $\mathbf{S}_{E}$ in Eq. (17a) does not require further
numerical approximation, since it is already sampled at time steps $m\Delta t$
based on the discretization of the forward field $\mathbf{E}$. (2) The
gradients terms in Eq. (19a) that include the products
$\partial_{\tau}\tilde{\mathbf{E}}\cdot\overleftarrow{\mathbf{E}}$ and
$\partial_{\tau}\tilde{\mathbf{E}}\cdot\Re\\{\overleftarrow{\mathbf{Q}}_{p}^{(i)}\\}$
can be computed in the most accurate way when evaluating them at time steps
$m\Delta t$, while staying consistent with using the first order forward
difference approximation for the time derivative. In return, the adjoint
source term $\mathbf{S}_{\partial_{\tau}Q_{p}}^{(i)}$ must be interpolated at
time steps $m\Delta t$. Therefore, a centered difference approximation was
employed. Compared to the forward equations, the update equations for both the
adjoint electric fields and the auxiliary fields include additional terms,
coming from the injection of the adjoint sources $\mathbf{S}_{E}$ and
$\mathbf{S}_{\partial_{\tau}Q_{p}}^{(i)}$.

Based on this time discretization of forward and adjoint equations, the
gradients in Eqs. (19) can be computed as follows:

$\displaystyle\mathcal{A}_{\rho\in\Omega}:=$ (28a) $\displaystyle-$
$\displaystyle\sum^{M}_{m=0}\varepsilon_{0}\mathrm{d}_{\rho}\varepsilon_{\infty}\mathbf{E}^{M-m}\cdot\left(\tilde{\mathbf{E}}^{m+\frac{1}{2}}-\tilde{\mathbf{E}}^{m-\frac{1}{2}}\right)$
$\displaystyle-$
$\displaystyle\sum^{M}_{m=0}\mathrm{d}_{\rho}\sigma\mathbf{E}^{M-m}\cdot\frac{1}{2}\left(\tilde{\mathbf{E}}^{m+\frac{1}{2}}+\tilde{\mathbf{E}}^{m-\frac{1}{2}}\right)\,\Delta
t$ $\displaystyle-$
$\displaystyle\sum^{M}_{m=0}\sum_{i=1}^{2}\sum_{p=1}^{P^{(i)}}2\mathrm{d}_{\rho}\kappa^{(i)}\Re\left\\{\mathbf{Q}^{M-m,(i)}_{p}\right\\}\cdot\left(\tilde{\mathbf{E}}^{m+\frac{1}{2}}-\tilde{\mathbf{E}}^{m-\frac{1}{2}}\right),$
$\displaystyle\mathcal{B}_{\rho\in\Omega_{\mathrm{o}}}:=$ (28b)
$\displaystyle\sum^{M}_{m=0}\frac{\mathrm{d}_{\rho}\sigma}{T}\left(\mathbf{E}^{M-m}\right)^{2}\Delta
t$ $\displaystyle+$
$\displaystyle\sum^{M}_{m=0}\sum_{i=1}^{2}\sum_{p=1}^{P^{(i)}}\mathrm{d}_{\rho}\kappa^{(i)}\Re\left\\{\frac{\left(\mathbf{Q}^{M-1-m,(i)}_{p}-\mathbf{Q}^{M-m+1,(i)}_{p}\right)^{2}}{T\,2\,\Delta
t\,\varepsilon_{0}\,c^{(i)}_{p}}\right\\}.$

From these equations and the definition of the adjoint source terms, we
observe that the electric field $\mathbf{E}$, as well as the auxiliary fields
$\mathbf{Q}_{p}^{(i)}$, must be stored in space and time during the forward
simulation. This implies a significant increase in memory consumption, the
more poles are used for an accurate material description of both background
and design materials.
For the spatial discretization across $N$ computational cells, we assign
density values $\rho_{n,k}$ on each position of three electric field
components $k\in\\{1,2,3\\}$ at the edges of the $n$th Yee cell. The gradients
in Eq. (28) are calculated per component such that a filtering of the original
density as well as of the gradients is essential to obtain a physical
meaningful structure [7].

Algorithm 1 outlines the time-domain gradient-based TopOpt routine using the
FDTD method.

Algorithm 1 Time-domain gradient-based TopOpt algorithm to maximize/minimize
the time-averaged power dissipation using the FDTD method

1:Initial density $\rho$, filter radius $R$, initial and maximum penalization
parameter $\beta_{0}$ and $\beta_{\mathrm{max}}$, threshold value $\eta$, the
maximum number of iterations $\text{itr}_{\text{max}}$ and time steps $M$.

2:

  * •

$M$ and $\text{itr}_{\text{max}}$ can be determined according to whether the
objective function $F$ and the design $\rho$ no longer show any significant
change, respectively.

  * •

Set a condition $\texttt{cond}_{\,\beta}$ for the increment of $\beta$.

3:

4:for $\text{itr}=1,\ldots,\text{itr}_{\text{max}}$ do

5: Compute the filtered and projected density $\bar{\tilde{\rho}}$

6: procedure Forward Simulation

7: for $m=0,\ldots,M$ do

8: Compute $\mathbf{H}^{m+\frac{1}{2}}$

9: Store $\mathbf{E}^{m}$ and all $\mathbf{Q}_{p}^{m,(i)}$

10: Compute $\mathbf{E}^{m+1}$

11: Compute all $\mathbf{Q}_{p}^{m+1,(i)}$

12: Update objective $F$

13: end for

14: end procedure

15: procedure Adjoint Simulation

16: for $m=0,\ldots,M$ do

17: Compute $\tilde{\mathbf{H}}^{m}$

18: Store current $\tilde{\mathbf{E}}^{m-\frac{1}{2}}$

19: Compute $\tilde{\mathbf{E}}^{m+\frac{1}{2}}$

20: Update $\nabla_{\bar{\tilde{\rho}}}F$

21: Compute all $\tilde{\mathbf{Q}}_{p}^{m+\frac{1}{2},(i)}$

22: end for

23: end procedure

24: Compute $\nabla_{\rho}F$

25: Update $\rho$

26: if $\texttt{cond}_{\,\beta}$ is true then

27: Increase $\beta$

28: end if

29:end for

30:Threshold $\bar{\tilde{\rho}}$ with respect to $\eta$
($\beta\rightarrow\infty$)

##  IV Spherical nanoparticles with enhanced absorption efficiency

Figure 5: Complex relative permittivity
($\varepsilon=\varepsilon^{\prime}-j\,\varepsilon^{\prime\prime}$) of (a)
Silicon and (b) Gold using the CCPR model in Eq. (1). Two poles were used for
the fit of Silicon, and three poles for Gold. The experimental data is
reported in the references [56, 57], respectively. The corresponding CCPR
parameters can be found in Ref. [7].

We applied our method to the inverse design of spherical nanoparticles for an
enhanced absorption efficiency in the visible and near ultraviolet regime. To
put the generic formulation of the optimization for arbitrary classes of
materials under test, we have chosen Gold and Silicon as representatives for
metallic and dielectric/semiconductive materials. The particles are embedded
in a non-absorbing background medium, which has been chosen as air.

This optimization problem was selected as our study case due to the following
reasons. The absorption efficiency $Q_{\text{abs}}(\omega)$ is proportional to
the absorbed power $\mathscr{W}_{\text{abs}}$ of the structure transferred
from an incident wave with intensity $I_{\text{inc }}$,

$Q_{\text{abs}}=\frac{\mathscr{W}_{\text{abs}}}{I_{\text{inc
}}\sigma_{\text{geom }}}.$ (29)

The geometrical cross section $\sigma_{\text{geom }}=4^{-1}d^{2}\,\pi$ is
determined by the radial expansion of the spherical particle of a diameter
$d$, and the absorbed power $\mathscr{W}_{\text{abs}}$ is the integral over
the time-averaged electric power dissipation density in Eq. (7) over the
design’s volume (considering $\Omega\equiv\Omega_{\mathrm{o}}$). Thus, the
value of $Q_{\text{abs}}$ over a broad range of frequencies reflects the
success of achieving a broadband performance of dissipation in a direct way.
All frequencies of interest can be covered in the source of excitation in our
time-domain optimization routine (Sec. IV-A).
We restrict the design to have a spherical boundary to enable a comparison to
analytical solutions using Mie theory and ensure sufficient accuracy of all
FDTD simulations. In addition, we use a sphere filled with the design material
both as an initial design and as a reference to compare the absorption
efficiencies of the final optimized designs. The given setup will highlight
the features evolving in 3D, which contribute to the absorption based on
different physical effects appearing in dielectric and plasmonic materials.
Therefore, we did not enforce explicit manufacturability constraints to allow
the design to develop freely in terms of its topological complexity.
Moreover, both Silicon and Gold are ideal candidates for the chosen wavelength
range of 300-700 nm with an imaginary part $\varepsilon^{\prime\prime}\gg 0$
in the dielectric function (Fig. 5), such that an absorbing behavior is
expected for both materials in general. The increasing imaginary part towards
the violet and near ultraviolet regime for both materials is mainly caused by
the direct band transitions in Silicon and the interband transitions in Gold,
respectively. The optimization for Gold is considered challenging and
represents an ultimate test for broadband optimization in time domain, since
both plasmonic effects and penetrating dissipative effects (due to interband
transitions) occur in different spectral ranges, respectively. We will see in
the following how the structures form topologically to best utilize their
complex material properties over the broad spectral range to enable efficient
absorption of energy.

Figure 6: Normalized objective functions (1st row) and measure of non-
discreteness (2nd row) per iteration while topology optimizing the (a) Silicon
and (b) Gold nanostructures. The normalization constant $F_{\text{Sphere}}$ is
the time-averaged power loss calculated for a sphere to whose boundaries the
designs were restricted in their expansion. The threshold was applied at the
final iteration, marked by a red circle, resulting in a binary design. Before
applying the threshold, the measured non-discreteness $M_{\text{nd}}$ of the
densities was $3.25\%$ for Silicon and $1.24\%$ for Gold.

Figure 7: (a.1) - (a.2): Smoothed illustrations of the optimized
nanostructures made of Silicon (left) and Gold (right) for enhanced
absorption. Both designs were restricted to evolve within a spherical volume
with a fixed diameter during the optimization. The systems are excited by a
$z$-polarized plane wave propagating along the $y$-axis. (b.1) - (b.2):
Absorption efficiencies of the design (in red) compared to a sphere with the
same diameter (black dashed) over the spectral range of interest. They were
computed by integrating the time-averaged electric power dissipation density
$\bar{q}_{e}$ in Eq. (7) over the design’s volume, and by using the Fourier-
transformed electric fields $\mathbf{E}(\omega)$ from the FDTD simulation. In
addition, the absorption efficiencies of a sphere predicted by Mie theory is
shown (black dotted). It served as verification of sufficient numerical
accuracy before commencing the optimization with identical simulation
parameters. (c.1) - (c.2): Spatial distributions of the amplitude of the
electric field $\|\mathbf{E}\|$ (1st row) and the time-averaged electric power
dissipation density $\bar{q}_{e}$ (2nd row) for the wavelength at which the
design has its absorption maximum. The amplitude of the electric field is
normalized to that of the incident field. The unit of the dissipation density
is given in $\mathrm{W}/\mathrm{m}^{3}$. A comparison with the spatial
distributions for a sphere is shown too.

Figure 8: Cross-sections through the optimized designs made of Silicon and
Gold, showing the amplitude of the electric field $\|\mathbf{E}\|$ and the
time-averaged electric power dissipation density $\bar{q}_{e}$ for the
wavelength at which the corresponding design has its absorption maximum. The
amplitude of the electric field is normalized to that of the incident field.
The unit of the dissipation density is given in $\mathrm{W}/\mathrm{m}^{3}$.
The planes were selected based on where the dissipation density shows
particularly strong local hotspots; see Figs. 7 (c.1)-(c.2).

We optimized a spherical Silicon and Gold nanoparticle surrounded by air
($\varepsilon^{(1)}=1$) for having an enhanced absorption efficiency
$Q_{\text{abs}}$ in the spectral range 300-700 nm. A summary of the parameters
used for the FDTD simulations and topology optimization can be found at the
end of this section (Sec. IV-A). The projection parameters and artificial
damping term $\gamma$ for Gold were selected based on the preliminary
considerations from Sec. III-B2. Both materials have been modelled by the CCPR
parameters listed in Ref. [7] and the corresponding fits are shown in Fig. 5.
Two poles were used for the fit of Silicon, and three poles for Gold. The
systems were excited by a $z$-polarized plane wave, covering the frequencies
of interest. We restricted the designs to a spherical volume with diameter
$d$, which was chosen to be $d=300$ nm for Silicon and $d=100$ nm for Gold.
Here, we compromised between highlighting material-dependent features evolving
due to different physical effects appearing in dielectric and plasmonic
structures, while ensuring numerical accuracy and a reasonable simulation time
and memory usage when storing the fields during the optimization. A sufficient
numerical accuracy was validated by simulating the response of a (binary)
sphere with the same setup employed in the design optimizations, followed by a
comparison of the computed dissipation via Eq. (7) against the analytical
results from Mie theory [58, 59].

The convergence of both structures is illustrated in Figs. 6 (a)-(b). The
figures show the normalized objective function with respect to that of a
sphere sharing the same diameter, and the measure of non-discreteness per
iteration step. In both cases, we started with a filtered and projected sphere
as the initial geometry. We let the Silicon and Gold design evolve for 50 and
100 iterations, respectively. A stable convergence behavior for both materials
can be observed, and they tend to evolve to a binary structure. In the case of
Gold, we recognize a jump in performance after the final threshold - and from
iteration 80 onwards, the design shows no significant improvement in terms of
performance and its tendency towards a binary structure. In fact, we observed
that the Gold structure retains localized gray areas on the inner surface of
the centrally formed cavity. Increasing the projection parameter $\beta$ could
not remedy this. As previously described in Sec. III-B2, we can attribute this
to the optimization driven by propagating, spatially confined but deeply
penetrating surface plasmons, and after the threshold, the localized surface
plasmon resonances come into full effect.

The results of the optimization are presented in Figs. 7 and 8, illustrating
the topologically designed structures, their performance compared to a sphere
with the same diameter $d$, and the electric field and local dissipation
distribution for the wavelengths of maximum absorption. The absorption
efficiencies $Q_{\text{abs}}$ from Eq. (29) have been computed by integration
of Eq. (7) over the design’s volume together with the electric fields from the
FDTD simulation (via discrete Fourier transformation), considering an
intensity of the incident wave of 1 $\;\mathrm{W/m^{2}}$. The geometrical
cross section was determined by the given spherical restriction of the
designs.

The spherical structure made of Silicon is characterized by a series of
interlocked material-free segments that enveloped at its surface, Fig. 7
(a.1). As we can see from the performance plot in Fig. 7 (b.1), the design
demonstrates a significant increase in absorption compared to that of a sphere
in the range 300-550 nm. We note that both absorption peaks of a sphere,
corresponding to the electric quadrupole and magnetic octupole modes in the
range 630-670 nm [60], are suppressed in the case of the topologically
optimized structure. They are correlated to resonances characterized by a high
quality factor, which are associated with a long-lasting duration of the
oscillations (Appendix C). As a consequence, it will effectively contribute to
a higher absorption of energy over time. We can explain the evolution of the
design primarily favoring the increase in absorption for $\lambda<$ 550 nm by
recalling the following: The maximization of the dissipation in the time
domain is equivalent to maximizing the integral value over the entire spectrum
of frequencies contained in the time-dependent source of excitation. The
evolving design is therefore not restricted to keep or produce sharp
resonances at specific wavelengths, as long as it yields a high absorption
over a broadband range in return. Furthermore, the imaginary part of the
permittivity of Silicon gets larger towards the ultraviolet regime. Thus, the
optimizer seems to exhibit a propensity to leverage this characteristic for
the maximization of the objective function. Reproducing single sharp
resonances associated with higher absorption is limited due to the nature of
the computational approach itself. In this case, the limitation is imposed by
the simulation time. That can be understood by noticing the discrepancy
between the numerically obtained peak at $\lambda\approx$ 645 nm of the sphere
and that predicted by Mie theory, Fig. 7 (b.1). We recall that an increase in
simulation time yields a more accurate result but also increases the time per
iteration and memory consumption due to the storage of fields (in space and
time) during the optimization. In order to better understand the optical
characteristics of the optimized Silicon design, we analyzed the resonance
spectrum of the design compared to a sphere in Appendix C. We observed that
our design’s topology supports the formation of a cascade of high-quality-
factor resonances in the wavelength range $<$ 550 nm, which get effectively
absorbed.
Fig. 7 (c.1) gives an insight into the spatial profile of the electric field
and dc component of the dissipation, Eq. (7), for the wavelength
$\lambda_{0}=$ 404 nm. Here, a comparison with a sphere is provided in
addition. The electric field’s amplitude is confined in the material-free
notches on the surface of the designed structure. By examining the spatial
dissipation distribution within the structure itself, we observe numerous
diffusely distributed hotspots that collectively enhance the absorption of the
electric field’s energy. A clearer picture of this phenomenon is provided by
Fig. 8, which shows cross-sections depicting the localization of particularly
strong hotspots.

The Gold design shown in Fig. 7 (a.2) appears as a more compact structure.
Three distinct holes on the surface on either side with respect to the
$z$-axis form the entrance to an inner cavity. From the performance plot in
Fig. 7 (b.2) we observe a $\approx$ 40% improvement at the absorption peak
corresponding to the electric dipole mode resonance, which is slightly shifted
towards larger wavelengths. The design yields an amplification of
$Q_{\text{abs}}$ across almost the entire spectrum compared to that of the
sphere. The spatial electric field and the dissipation profile for
$\lambda_{0}=$ 538 nm in Fig. 7 (c.2) and their cross-sections in Fig. 7 (c.2)
provide insight into the topological formation of the cavity with regard to
the plasmonic effects. The electric field is mainly confined at the surface,
where thin walls separate the outside from the inner cavity. In these walls,
the local dissipation has its maximum. This is understandable if we recall the
physical nature of plasmons. The topology enables a balance between their
formation and longevity along the surface, and their penetration into the
material itself, which is limited by the skin depth inside the metal.

We conclude that the optimization worked successfully for both material types,
and a broadband performance could be achieved in each case. We note that the
choice of the artificial damping parameter $\gamma$ for Gold was crucial, as a
too small or large value impairs the convergence or leads to the design
remaining gray. Despite the remarkable broadband performance of the optimized
Gold structure, we observed that the corresponding density retained thin gray
areas along the surface of the inner cavity. Neither an increase in iteration
steps nor modulation of the projection parameter $\beta$ could circumvent it.
In the case of Silicon, we performed additional optimizations with different
imposed diameters of the spherical restriction (Appendix C). We found that as
the diameter was increased, the optimized designs yield an improvement in
performance towards larger wavelengths. This is plausible as more space for
the design to develop is provided. Thus, the evolving structure can utilise
the dissipative effects related to resonances at longer wavelengths.

###  IV-A Optimization and simulation parameters

Figure 9: (a) Time and (b) frequency domain plots of the excitation source’s
amplitude. The area in yellow marks the frequency range of interest, bounded
by the values $\omega_{1}=6.28\times 10^{15}\,\mathrm{rad/s}$ and
$\omega_{2}=2.7\times 10^{15}\,\mathrm{rad/s}$ of the half bandwidth, which
correspond to the wavelengths of 300 nm and 700 nm, respectively.

TABLE I: TopOpt and FDTD configuration settings used for designing the
spherical nanostructures presented in Sec. IV

Parameter |  Silicon |  Gold
---|---|---
Simulation domain |  $125\times 125\times 125$ cells |  $100\times 100\times 100$ cells
Design domain $\Omega$ |  $86\times 85\times 86$ cells |  $61\times 60\times 61$ cells
Observation domain $\Omega_{\mathrm{o}}$ |  $86\times 85\times 86$ cells |  $61\times 60\times 61$ cells
CPML |  25 cells |  20 cells
d$x$, d$y$, d$z$ |  $4$ nm |  $2$ nm
d$t$ |  6.67 as |  3.34 as
Filter radius $R$ |  12 nm |  10 nm
Time steps per sim. $M$ |  22000 |  12000
$\eta$ |  0.5 |  0.55
$\beta_{0}$ |  10 |  7
$\beta_{\mathrm{max}}$ |  20 |  42
$\gamma$ |  0 |  $2\times 10^{5}$
TopOpt iterations |  50 |  100

Table I shows the configuration of both TopOpt and FDTD parameters used for
the optimizations presented in the previous section. The simulation domain was
surrounded by convolutional perfectly matched layers (CPML) to mimic open
boundaries by minimizing the reflection from the edges of the simulation grid.
Design domain $\Omega$ and observation region $\Omega_{\mathrm{o}}$ were
placed in the center. They were both cubic and shared the same dimensions,
i.e., $\Omega\equiv\Omega_{\mathrm{o}}$. We chose a homogeneous sphere with
$\rho\equiv\eta$ as our initial geometry, which was subsequently filtered and
projected. We started with an initial projection parameter $\beta_{0}$ and let
the value increase gradually during the optimization until it reached a
maximum value of $\beta_{\mathrm{max}}$. During the forward simulations, the
system was excited by a $z$-polarized plane wave carrying a truncated sinc
signal covering the frequencies of interest, Fig. 9. The number of time steps
$M$ per forward and adjoint simulation (Algorithm 1) was determined by whether
the absorption spectrum of a sphere matches that of the Mie theory with
sufficient accuracy (see Fig. 7 (b.1)-(b.2)). In addition, it was ensured that
there was no significant change in the objective function during the
optimization after this time period (due to the decay of the fields). We let
the Silicon and Gold design evolve for 50 and 100 iterations, respectively.
After each iteration, the density was updated using the Method of Moving
Asymptotes (MMA) [61] and subsequently stenciled out via an indicator
function, $\rho\cdot\mathbf{1}_{\mathcal{B}_{d/2}(\bm{r}_{0})}$, to restrict
the evolving design to a spherical volume with a fixed diameter $d$. Here,
$\bm{r}_{0}$ denotes the center point of the domain $\Omega$, and
$\mathcal{B}_{d/2}(\bm{r}_{0})$ describes a sphere with radius $d/2$ around
$\bm{r}_{0}$. The diameter was chosen to be $d=300$ nm for Silicon and $d=100$
nm for Gold.
The optimizations were performed with our parallelized in-house FDTD software
with TopOpt implementation [7] on the supercomputer HLRN-IV-System Emmy in
Göttingen, Germany, provided by the North German Supercomputing Alliance as
part of the National High Performance Computing (NHR) infrastructure. The
simulation domain has been decomposed into multiple chunks [62] in order to
accelerate the optimization and to ensure an efficient memory distribution
when storing the electric and auxiliary fields in space and time. The
optimization of the Silicon structure presented in Sec. IV took $\approx$ 11
hours on 2304 CPUs, while the optimization for the Gold structure took
$\approx$ 10 hours using 1152 CPUs (Cascade 9242).

##  V Conclusion

We introduced a novel time-domain gradient-based topology optimization method
to locally optimize the electric power dissipation for dispersive materials
over a broadband range of frequencies. The method is formulated based on the
CCPR model to accurately describe the material’s dispersion within any desired
spectral range. By using the auxiliary differential equation (ADE) method, we
derived the expression for the instantaneous electric power dissipation based
on the CCPR model, and formulated an adjoint scheme to obtain the gradient
information of the time-averaged dissipation defined as our objective
function. The discretization and integration of the adjoint system into the
finite-difference time-domain (FDTD) framework was explained in detail. Our
method was successfully demonstrated on the optimization of spherical
nanoparticles made of Silicon and Gold for an enhanced absorption efficiency
in the visible and ultraviolet regimes.

The adjoint scheme presented in Sec. III-A is kept in a general form and is
thus not restricted to the specific setup used in our study case. It further
allows for other configurations, such as imposing periodic boundary conditions
to enable the design of broadband absorbing metasurfaces, considering both
background and design materials as dispersive to optimize hybrid
nanostructures, or maximizing dissipation in a specific subvolume
$\Omega_{\mathrm{o}}\subset\Omega$ to generate localized heat sources (e.g.,
for thermal ablation). We note that despite the success, the optimization
comes with challenges, as it was shown by the theoretical consideration of the
plasmonic effects of Gold during optimization and its convergence behavior in
practice. We can attribute this to the density-based approach in the context
of maximizing dissipation per se, which requires appropriate tuning with
respect to a material parameter for intermediate densities. We assume that a
refinement of the material interpolation can lead to even better results. It
is also possible that the adaptation and implementation of the latest
achievements in shape optimization could be promising for tackling this
problem [63], which completely dispenses with the concept of material
interpolation.

Our contribution holds great potential for tackling inverse design problems,
where dispersive materials are included, and a reduced or enhanced absorption
over a broadband range is desired. Areas of application include solar energy
harvesting, by maximizing light absorption in photovoltaic solar cells or
improving thermal radiation absorption in thermophotovoltaics; broadband
metasurface absorbers, selective plasmonic sensors and absorptive filters;
anti-reflective or high-absorption coatings; devices with thermal tuning
exploiting thermal effects to vary the refractive index; or designing low-loss
devices like waveguides, couplers, and antennas with improved efficiency.

##  Appendix A Time average of the electric power dissipation density for the
time-harmonic case

In the time-harmonic case, the instantaneous electric power dissipation
density in Eq. (5) can be split into a dc and a sinusoidal ac term [44]. To
measure the broadband performance of an optimized device in the steady state
regarding the time-average of that dissipation density, we extract its dc
component here.
Consider a time harmonic field
$\mathbf{E}(t):=\Re\\{\mathbf{\hat{E}}_{0}e^{j\omega_{0}t}\\}$ oscillating
with a frequency $\omega_{0}$, where $\hat{\mathbf{E}}_{0}$ itself is
frequency independent. The Fourier transform yields

$\hat{\mathbf{E}}(\omega)=\frac{1}{2}\left(\hat{\mathbf{E}}_{0}\,\delta\left(\omega_{0}-\omega\right)+\hat{\mathbf{E}}_{0}^{*}\,\delta\left(\omega_{0}+\omega\right)\right).$
(30)

The auxiliary field components in the frequency domain based on the CCPR model
in Eq. (1) are related to the electric field components via [42]

$\hat{\mathbf{Q}}_{p}(\omega)=\frac{\varepsilon_{0}c_{p}}{j\omega-
a_{p}}\hat{\mathbf{E}}(\omega).$ (31)

Thus, based on the time-harmonic field, the auxiliary field in time domain and
its derivative read

$\displaystyle\mathbf{Q}_{p}(t)=$
$\displaystyle\int_{-\infty}^{+\infty}\hat{\mathbf{Q}}_{p}(\omega)e^{j\omega
t}d\omega$ (32) $\displaystyle=$
$\displaystyle\,\frac{1}{2}\left(\frac{\varepsilon_{0}c_{p}}{j\omega_{0}-a_{p}}\hat{\mathbf{E}}_{0}e^{j\omega_{0}t}+\frac{\varepsilon_{0}c_{p}}{-j\omega_{0}-a_{p}}\hat{\mathbf{E}}_{0}^{*}e^{-j\omega_{0}t}\right),$
$\displaystyle\partial_{t}\mathbf{Q}_{p}(t)=$
$\displaystyle\,\frac{1}{2}\left(\frac{j\omega_{0}\varepsilon_{0}c_{p}}{j\omega_{0}-a_{p}}\hat{\mathbf{E}}_{0}e^{j\omega_{0}t}+\frac{-j\omega_{0}\varepsilon_{0}c_{p}}{-j\omega_{0}-a_{p}}\hat{\mathbf{E}}_{0}^{*}e^{-j\omega_{0}t}\right).$
(33)

Substituting $\partial_{t}\mathbf{Q}_{p}$ into our formula for the
instantaneous electric power dissipation density,

$q_{e}(t)=\sigma\mathbf{E}^{2}+2\sum_{p=1}^{P}\Re\left\\{\frac{(\partial_{t}\mathbf{Q}_{p})^{2}}{\varepsilon_{0}c_{p}}\right\\},$
(34)

yields an expression consisting of both the ac and dc terms. Since we consider
the time-average of the dissipation as our objective, we are only interested
in its dc part. Using the identity

$\Re\\{\mathbf{W}\\}^{2}=\frac{1}{4}\left(\mathbf{W}+\mathbf{W}^{*}\right)^{2}=\frac{1}{2}|\mathbf{W}|^{2}+\frac{1}{2}\Re\left\\{\mathbf{W}^{2}\right\\},$
(35)

that holds for any complex vector field $\mathbf{W}$, and neglecting all the
terms containing the factors $e^{\pm 2j\omega t}$, leaves the dc component,

$\displaystyle\bar{q}_{e}(\omega_{0})=\left(\frac{1}{2}\sigma+\sum_{p=1}^{P}\Re\left\\{\frac{\omega_{0}^{2}\varepsilon_{0}c_{p}}{(j\omega_{0}-a_{p})(-j\omega_{0}-a_{p})}\right\\}\right)|\hat{\mathbf{E}}_{0}|^{2}.$
(36)

We can further express this result in terms of the imaginary part of the
permittivity $\varepsilon$. By noticing that

$\displaystyle\,\frac{1}{2j}\left(\frac{c_{p}}{j\omega_{0}-a_{p}}-\frac{c_{p}}{-j\omega_{0}-a_{p}}\right)+\mathrm{c.c.}$
(37) $\displaystyle=$
$\displaystyle\,\frac{1}{2j}\left\\{\left(\frac{c_{p}}{j\omega_{0}-a_{p}}+\frac{c_{p}^{*}}{\phantom{-}j\omega_{0}-a_{p}^{*}}\right)-\mathrm{c.c.}\right\\},$
(38)

and using the identities $\Re\\{z\\}=\frac{1}{2}(z+z^{*})$ and
$\Im\\{z\\}=\frac{1}{2j}(z-z^{*})$ for complex numbers $z$, the prefactor in
brackets in Eq. (36) can be written as

$\displaystyle\;\frac{\sigma}{2}+\sum_{p=1}^{P}\Re\left\\{\frac{\omega_{0}^{2}\varepsilon_{0}c_{p}}{(j\omega_{0}-a_{p})(-j\omega_{0}-a_{p})}\right\\}$
(39) $\displaystyle=$
$\displaystyle\;\frac{\sigma}{2}-\varepsilon_{0}\omega_{0}\sum_{p=1}^{P}\Re\left\\{\frac{1}{2j}\left(\frac{c_{p}}{j\omega_{0}-a_{p}}-\frac{c_{p}}{-j\omega_{0}-a_{p}}\right)\right\\}$
(40) $\displaystyle=$
$\displaystyle-\frac{1}{2}\varepsilon_{0}\,\omega_{0}\,\Im\left\\{\varepsilon(\omega_{0})\right\\}.$
(41)

From that, we observe that Eq. (36) is equivalent to the model-independent
expression derived by Landau and Lifschitz in Ref. [13],

$\bar{q}_{e}(\omega_{0})=\frac{1}{2}\varepsilon_{0}\,\omega_{0}\,\Im\left\\{\varepsilon(\omega_{0})\right\\}|\hat{\mathbf{E}}_{0}|^{2},$
(42)

considering the convention
$\varepsilon=\varepsilon^{\prime}-j\varepsilon^{\prime\prime}$.

##  Appendix B Discretization of the adjoint equations

Here, we present the update equations for the adjoint fields in time,
considering the time-averaged dissipation in Eq. (9) to be the objective
function. The update scheme for the forward equations is identical, except for
the time shift of $\Delta t/2$ and the absence of the adjoint source terms. In
that case, the equations are reduced to those presented in Ref. [42]. We again
consider isotropic media for both background and design material, and denote
these materials with the indices $i=1,2$, respectively. The material
parameters $\varepsilon_{\infty}$, $\sigma$ appearing in the following
equations are assumed to be material interpolated by the density $\rho$, and
the prefactor $\kappa^{(i)}$ is defined as $\kappa^{(1)}(\rho):=1-\rho$ and
$\kappa^{(2)}(\rho):=\rho$. The extension to anisotropic media can be found in
the Supplementary Material.

We discretize Eq. (17a) and all Eqs. (17b) at time steps $m\Delta t$. Since
the electric and auxiliary fields of the forward system are considered to be
discretized at time steps $m\Delta t$ too (Fig. 4), we can discretize the
adjoint source term $\mathbf{S}_{E}$ depending on the forward electric field
as

$\displaystyle\mathbf{S}^{m}_{E}=2T^{-1}\sigma\overleftarrow{\mathbf{E}}\left[m\Delta
t\right]=2T^{-1}\sigma\mathbf{E}^{M-m},$ (43)

and each adjoint source term $\mathbf{S}_{\partial{\tau}Q_{p}}^{m,(i)}$
depending on the time-derivative of the corresponding forward auxiliary field
as

$\displaystyle\mathbf{S}_{\partial{\tau}Q_{p}}^{m,(i)}=$
$\displaystyle\;2T^{-1}\partial_{\tau}\overleftarrow{\mathbf{Q}}_{p}^{(i)}\left[m\Delta
t\right]$ (44) $\displaystyle=$
$\displaystyle\;2T^{-1}\;\frac{\left(\mathbf{Q}_{p}^{M-m-1,(i)}-\mathbf{Q}_{p}^{M-m+1,(i)}\right)}{2\Delta
t},$ (45)

where the centered difference approximation was employed. The update equation
for each adjoint auxiliary field is then

$\displaystyle\tilde{\mathbf{Q}}_{p}^{m+\frac{1}{2},(i)}=$
$\displaystyle\;\frac{2+a_{p}^{(i)}\Delta t}{2-a_{p}^{(i)}\Delta
t}\tilde{\mathbf{Q}}_{p}^{m-\frac{1}{2},(i)}$ (46) $\displaystyle+$
$\displaystyle\;\frac{\varepsilon_{0}c_{p}^{(i)}\Delta t}{2-a_{p}^{(i)}\Delta
t}\left(\tilde{\mathbf{E}}^{m+\frac{1}{2}}+\tilde{\mathbf{E}}^{m-\frac{1}{2}}\right)$
$\displaystyle+$ $\displaystyle\;\frac{2\Delta t}{2-a_{p}^{(i)}\Delta
t}\mathbf{S}_{\partial{\tau}\mathbf{Q}_{p}}^{m,(i)}.$

Using these expressions, the update equation of the electric field reads

$\tilde{\mathbf{E}}^{m+\frac{1}{2}}=\alpha^{-1}\mathbf{g}^{m-\frac{1}{2}},$
(47)

where $\alpha$ is a density-dependent constant defined as

$\displaystyle\alpha:=\frac{\varepsilon_{0}\varepsilon_{\infty}}{\Delta
t}+\frac{\sigma}{2}+\frac{2}{\Delta
t}\sum_{i=1}^{2}\sum_{p=1}^{P^{(i)}}\kappa^{(i)}\Re\left\\{\frac{\varepsilon_{0}c_{p}^{(i)}\Delta
t}{2-a_{p}^{(i)}\Delta t}\right\\},$ (48)

and $\mathbf{g}^{m-\frac{1}{2}}$ changes at every time step according to

$\displaystyle\mathbf{g}^{m-\frac{1}{2}}=$
$\displaystyle\;(\nabla\times\tilde{\mathbf{H}})^{m}+\left(\frac{2\varepsilon_{0}\varepsilon_{\infty}}{\Delta
t}-\alpha\right)\tilde{\mathbf{E}}^{m-\frac{1}{2}}$ (49) $\displaystyle-$
$\displaystyle\;\frac{2}{\Delta
t}\sum_{i=1}^{2}\sum_{p=1}^{P^{(i)}}\kappa^{(i)}\Re\left\\{\frac{2a_{p}^{(i)}\Delta
t}{2-a_{p}^{(i)}\Delta t}\tilde{\mathbf{Q}}_{p}^{m-\frac{1}{2},(i)}\right\\}$
$\displaystyle-$ $\displaystyle\;\frac{2}{\Delta
t}\sum_{i=1}^{2}\sum_{p=1}^{P^{(i)}}\kappa^{(i)}\Re\left\\{\frac{2\Delta
t}{2-a_{p}^{(i)}\Delta
t}\mathbf{S}_{\partial{\tau}\mathbf{Q}_{p}}^{m,(i)}\right\\}$ $\displaystyle+$
$\displaystyle\;\mathbf{S}^{m}_{E}.$

##  Appendix C Mie analysis of the Silicon TopOpt design

For a better understanding of the optical properties of the optimized Silicon
design from Sec. IV, we examined the resonances of the optimized structure
compared to the sphere (Fig. 10). To identify the resonances of both
structures, we illustrate not only their scattering efficiency
$Q_{\mathrm{sca},\,\mathrm{lossy}}$ and absorption efficiency
$Q_{\mathrm{abs}}$, but also the scattering efficiency
$Q_{\mathrm{sca},\,\mathrm{lossless}}$ for the artificial undamped case, i.e.,
setting $\Im\\{\varepsilon\\}=0$. Here, the spectral range is limited to
400-700 nm, in which it was observed that the optimized design tends to
outperform the sphere in terms of absorption efficiency towards smaller
wavelengths (see Fig. 7 (b.1)). From the
$Q_{\mathrm{sca},\,\mathrm{lossless}}$ plots, it can be observed that our
design exhibits a denser distribution of resonances in the range of 400-550 nm
compared to the sphere. The sharp peaks (high-$Q_{f}$ resonances) are
associated with higher quality factors $Q_{f}$, which correspond to a longer
lifetime $\tau=2Q_{f}(\omega)/\omega$ of the resonances [64]. As a
consequence, more energy can be absorbed in this period of time if we include
the damping, which effectively contributes to increasing our objective
function in Eq. (9). Considering the damping of Silicon, we see a stronger
reduction of scattering and an increase of absorption compared to the sphere,
which indicates an efficient absorption of the resonances. Figure 10 also
shows why the optimized design exhibits relatively weak absorption in the
630-700 nm range. High-$Q_{f}$ resonances, as in the sphere case, are
completely absent here.

From these observations, we can conclude that the topological shaping of the
design leads to the formation of a cascade of high-$Q_{f}$ resonances, which
in turn enhance the absorption. It is therefore reasonable to analyze the
performance of the optimized designs in relation to the design volume.
Restricting the spherical design to a smaller diameter reduces the degrees of
freedom available for the development of the topology, which is needed for the
emergence of multiple high-Q resonances across the entire spectral range of
300-700 nm. To verify that, we performed topology optimizations of spherical
Silicon particles for different diameters. To ensure an optimal trade-off
between numerical precision, computational time, and memory consumption, we
restricted our optimizations to diameters up to 300 nm. The results are
illustrated in Fig. 11. We can clearly see that an increase in the diameter
leads to an improvement in broadband performance. The absorption curve
corresponding to a diameter of 300 nm is the performance of the design
presented in Sec. IV. We can expect that increasing the radius further will
result in even better broadband absorption, as the increased spatial extent
and topological complexity allow the appearance of resonances at longer
wavelengths too.

Figure 10: Scattering efficiency $Q_{\mathrm{sca},\,\mathrm{lossy}}$ (orange)
and absorption efficiency $Q_{\mathrm{abs}}$ (red) of both Silicon (a) sphere
and (b) TopOpt design presented in Sec. (IV) for the wavelength range 400-700
nm. To identify the resonances, the scattering efficiency
$Q_{\mathrm{sca},\,\mathrm{lossless}}$ (blue) is also shown here for the
undamped case, i.e., setting $\Im\\{\varepsilon\\}=0$.

Figure 11: Absorption efficiencies $Q_{\text{abs}}$ of topology optimized
nanostructures made of Silicon for different diameters, by which the designs
were restricted in their spatial extension. The absorption efficiency
corresponding to a diameter of 300 nm corresponds to the design presented in
Sec. IV.

## References

  * [1] M. P. Bendsøe and O. Sigmund, _Topology Optimization_. Springer Berlin Heidelberg, 2004.
  * [2] J. Jensen and O. Sigmund, “Topology optimization for nano-photonics,” _Laser & Photonics Reviews_, vol. 5, no. 2, pp. 308–321, 2011.
  * [3] L. Zhang, C. Wang, Y. Wei, Y. Lin, Y. Han, and Y. Deng, “High-efficiency achromatic metalens topologically optimized in the visible,” _Nanomaterials_ , vol. 13, no. 5, 2023.
  * [4] T. Nomura, K. Sato, K. Taguchi, T. Kashiwa, and S. Nishiwaki, “Structural topology optimization for the design of broadband dielectric resonator antennas using the finite difference time domain technique,” _International Journal for Numerical Methods in Engineering_ , vol. 71, no. 11, pp. 1261–1296, 2007.
  * [5] A. Erentok and O. Sigmund, “Topology optimization of sub-wavelength antennas,” _IEEE Trans. Antennas Propag._ , vol. 59, no. 1, pp. 58 –69, Jan. 2011.
  * [6] J. Wang, X.-S. Yang, X. Ding, and B.-Z. Wang, “Antenna radiation characteristics optimization by a hybrid topological method,” _IEEE Trans. Antennas Propag._ , vol. 65, no. 6, pp. 2843–2854, Jun. 2017.
  * [7] J. Gedeon, E. Hassan, and A. Calà Lesina, “Time-domain topology optimization of arbitrary dispersive materials for broadband 3D nanophotonics inverse design,” _ACS Photonics_ , vol. 10, no. 11, pp. 3875–3887, 2023.
  * [8] E. Hassan and A. Calà Lesina, “Topology optimization of dispersive plasmonic nanostructures in the time-domain,” _Opt. Express_ , vol. 30, no. 11, pp. 19 557–19 572, May 2022.
  * [9] O. Yesilyurt, Z. A. Kudyshev, A. Boltasseva, V. M. Shalaev, and A. V. Kildishev, “Efficient topology-optimized couplers for on-chip single-photon sources,” _ACS Photonics_ , vol. 8, no. 10, p. 3061–3068, Oct. 2021.
  * [10] T. W. Hughes, M. Minkov, I. A. D. Williamson, and S. Fan, “Adjoint method and inverse design for nonlinear nanophotonic devices,” _ACS Photonics_ , vol. 5, no. 12, pp. 4781–4787, Dec. 2018.
  * [11] E. Hassan, B. Scheiner, F. Michler, M. Berggren, E. Wadbro, F. Röhrl, S. Zorn, R. Weigel, and F. Lurz, “Multilayer topology optimization of wideband SIW-to-waveguide transitions,” _IEEE Transactions on Microwave Theory and Techniques_ , vol. 68, no. 4, pp. 1326–1339, 2020.
  * [12] N. V. Sapra, K. Y. Yang, D. Vercruysse, K. J. Leedle, D. S. Black, R. J. England, L. Su, R. Trivedi, Y. Miao, O. Solgaard, R. L. Byer, and J. Vučković, “On-chip integrated laser-driven particle accelerator,” _Science_ , vol. 367, no. 6473, pp. 79–83, 2020.
  * [13] L. Landau, E. Lifshitz, and L. Pitaevskii, _Electrodynamics of Continuous Media: Volume 8_ , ser. Course of theoretical physics. Elsevier Science, 1995.
  * [14] Y. Wang, G. Wu, Y. Wang, and J. Liu, “Inverse design of a broadband tunable terahertz metasurface absorber,” _Optics Communications_ , vol. 540, p. 129526, 2023.
  * [15] V. Ganapati, O. D. Miller, and E. Yablonovitch, “Light trapping textures designed by electromagnetic optimization for subwavelength thick solar cells,” _IEEE Journal of Photovoltaics_ , vol. 4, no. 1, pp. 175–182, 2014\.
  * [16] A. Azad, W. Kort-Kamp, M. Sykora, N. Weisse-Bernstein, T. Luk, A. Taylor, D. Dalvit, and H.-T. Chen, “Metasurface broadband solar absorber,” _Scientific Reports_ , vol. 6, 09 2015.
  * [17] X. Jiang, H. Yuan, D. Chen, Z. Zhang, T. Du, H. Ma, and J. Yang, “Metasurface based on inverse design for maximizing solar spectral absorption,” _Advanced Optical Materials_ , vol. 9, no. 19, Jun. 2021.
  * [18] S. Shrestha, Y. Wang, A. C. Overvig, M. Lu, A. Stein, L. D. Negro, and N. Yu, “Indium tin oxide broadband metasurface absorber,” _ACS Photonics_ , vol. 5, no. 9, p. 3526–3533, Aug. 2018.
  * [19] H. Liu, M. Xie, Q. Ai, and Z. Yu, “Ultra-broadband selective absorber for near-perfect harvesting of solar energy,” _Journal of Quantitative Spectroscopy and Radiative Transfer_ , vol. 266, p. 107575, 2021.
  * [20] A. Starczewska and M. Kepinska, “Photonic crystal structures for photovoltaic applications,” _Materials_ , vol. 17, no. 5, p. 1196, Mar. 2024.
  * [21] V. Rinnerbauer, Y. Shen, J. D. Joannopoulos, M. Soljačić, F. Schäffler, and I. Celanovic, “Superlattice photonic crystal as broadband solar absorber for high temperature operation,” _Opt. Express_ , vol. 22, no. S7, pp. A1895–A1906, Dec 2014.
  * [22] J. Dong, X.-H. Deng, P. Zhang, and J. Yuan, “Multi-channel THz perfect absorber using graphene-based fibonacci photonic crystals,” _J. Opt. Soc. Am. B_ , vol. 40, no. 3, pp. B28–B34, Mar 2023.
  * [23] M. Omrani, R. Keshavarzi, M. Abdi-Jalebi, and P. Gao, “Impacts of plasmonic nanoparticles incorporation and interface energy alignment for highly efficient carbon-based perovskite solar cells,” _Scientific Reports_ , vol. 12, no. 1, Mar. 2022.
  * [24] A. Paris, A. Vaccari, A. Calà Lesina, E. Serra, and L. Calliari, “Plasmonic scattering by metal nanoparticles for solar cells,” _Plasmonics_ , vol. 7, 09 2012.
  * [25] L. Mascaretti, Y. Chen, O. Henrotte, O. Yesilyurt, V. M. Shalaev, A. Naldoni, and A. Boltasseva, “Designing metasurfaces for efficient solar energy conversion,” _ACS Photonics_ , vol. 10, no. 12, p. 4079–4103, Dec. 2023\.
  * [26] G. Baffou, F. Cichos, and R. Quidant, “Applications and challenges of thermoplasmonics,” pp. 946–958, 9 2020.
  * [27] G. P. Zograf, M. I. Petrov, S. V. Makarov, and Y. S. Kivshar, “All-dielectric thermonanophotonics,” _Advances in Optics and Photonics_ , vol. 13, p. 643, 9 2021.
  * [28] S. Cilla, A. Ianiro, C. Romano, F. Deodato, G. Macchia, M. Buwenge, N. Dinapoli, L. Boldrini, A. G. Morganti, and V. Valentini, “Template-based automation of treatment planning in advanced radiotherapy: a comprehensive dosimetric and clinical evaluation,” _Scientific Reports_ , vol. 10, no. 1, Jan. 2020.
  * [29] K.-T. Lin, H.-L. Chen, Y.-S. Lai, and C.-C. Yu, “Silicon-based broadband antenna for high responsivity and polarization-insensitive photodetection at telecommunication wavelengths,” _Nature Communications_ , vol. 5, no. 1, Feb. 2014.
  * [30] K. Chen, R. Adato, and H. Altug, “Dual-band perfect absorber for multispectral plasmon-enhanced infrared spectroscopy,” _ACS Nano_ , vol. 6, no. 9, p. 7998–8006, Aug. 2012.
  * [31] D. Shrekenhamer, W. Xu, S. Venkatesh, D. Schurig, S. Sonkusale, and W. J. Padilla, “Experimental realization of a metamaterial detector focal plane array,” _Phys. Rev. Lett._ , vol. 109, p. 177401, Oct 2012.
  * [32] J.-B. Brückner, J. L. Rouzo, L. Escoubas, G. Berginc, O. Calvo-Perez, N. Vukadinovic, and F. Flory, “Metamaterial filters at optical-infrared frequencies,” _Opt. Express_ , vol. 21, no. 14, pp. 16 992–17 006, Jul 2013.
  * [33] Y. Yao, Z. Liao, Z. Liu, X. Liu, J. Zhou, G. Liu, Z. Yi, and J. Wang, “Recent progresses on metamaterials for optical absorption and sensing: a review,” _Journal of Physics D: Applied Physics_ , vol. 54, no. 11, p. 113002, Jan. 2021.
  * [34] J. Hong, H. Son, C. Kim, S.-E. Mun, J. Sung, and B. Lee, “Absorptive metasurface color filters based on hyperbolic metamaterials for a CMOS image sensor,” _Opt. Express_ , vol. 29, no. 3, pp. 3643–3658, Feb 2021\.
  * [35] E. Hassan, A. B. Evlyukhin, and A. Calà Lesina. (2023) Anapole plasmonic meta-atom enabled by inverse design for metamaterials transparency. arXiv: 2303.00478v1. Accessed on March 1, 2023. [Online]. Available: https://doi.org/10.48550/arXiv.2303.00478
  * [36] H. Yasuda and S. Nishiwaki, “A design method of broadband metalens using time-domain topology optimization,” _AIP Advances_ , vol. 11, no. 5, May 2021\.
  * [37] J. Baxter and L. Ramunno, “Inverse design of optical pulse shapes for time-varying photonics,” _Opt. Express_ , vol. 31, no. 14, pp. 22 671–22 684, Jul 2023.
  * [38] R. J. Tang, S. W. D. Lim, M. Ossiander, X. Yin, and F. Capasso, “Time reversal differentiation of FDTD for photonic inverse design,” _ACS Photonics_ , vol. 10, no. 12, p. 4140–4150, Nov. 2023.
  * [39] J. H. Greene and A. Taflove, “General vector auxiliary differential equation finite-difference time-domain method for nonlinear optics,” _Opt. Express_ , vol. 14, no. 18, pp. 8305–8310, Sep 2006.
  * [40] K. P. Prokopidis and D. C. Zografopoulos, “Investigation of the stability of ADE-FDTD methods for modified lorentz media,” _IEEE Microwave and Wireless Components Letters_ , vol. 24, no. 10, pp. 659–661, 2014.
  * [41] O. Ramadan, “Stability-improved ADE-FDTD implementation of Drude dispersive models,” _IEEE Antennas and Wireless Propagation Letters_ , vol. 17, no. 5, pp. 877–880, 2018.
  * [42] K. P. Prokopidis and D. C. Zografopoulos, “Time-domain studies of general dispersive anisotropic media by the complex-conjugate pole–residue pairs model,” _Applied Sciences_ , vol. 11, no. 9, p. 4, 2021.
  * [43] M. Han, R. Dutton, and S. Fan, “Model dispersive media in finite-difference time-domain method with complex-conjugate pole-residue pairs,” _IEEE Microwave and Wireless Components Letters_ , vol. 16, no. 3, pp. 119–121, 2006\.
  * [44] W. Shin, A. Raman, and S. Fan, “Instantaneous electric energy and electric power dissipation in dispersive media,” _J. Opt. Soc. Am. B_ , vol. 29, no. 5, pp. 1048–1054, May 2012.
  * [45] F. D. Nunes, T. C. Vasconcelos, M. Bezerra, and J. Weiner, “Electromagnetic energy density in dispersive and dissipative media,” _J. Opt. Soc. Am. B_ , vol. 28, no. 6, pp. 1544–1552, Jun 2011.
  * [46] B. Bourdin, “Filters in topology optimization,” _International Journal for Numerical Methods in Engineering_ , vol. 50, no. 9, pp. 2143–2158, 2001.
  * [47] E. Hassan, E. Wadbro, and M. Berggren, “Patch and ground plane design of microstrip antennas by material distribution topology optimization,” _Progress In Electromagnetics Research B_ , vol. 59, p. 89–102, 2014.
  * [48] F. Wang, B. S. Lazarov, and O. Sigmund, “On projection methods, convergence and robust formulations in topology optimization,” _Structural and Multidisciplinary Optimization_ , vol. 43, pp. 767–784, 2011.
  * [49] E. Hassan, E. Wadbro, and M. Berggren, “Topology optimization of metallic antennas,” _IEEE Trans. Antennas Propag._ , vol. 62, no. 5, pp. 2488–2500, May 2014.
  * [50] O. Sigmund, “Morphology-based black and white filters for topology optimization,” _Struct. Multidiscip. Optim._ , vol. 33, no. 4–5, pp. 401–424, 2007.
  * [51] J. S. Jensen and O. Sigmund, “Topology optimization of photonic crystal structures: a high-bandwidth low-loss T-junction waveguide,” _Journal of the Optical Society of America B_ , vol. 22, no. 6, p. 1191, Jun. 2005.
  * [52] J. Andkjær, S. Nishiwaki, T. Nomura, and O. Sigmund, “Topology optimization of grating couplers for the efficient excitation of surface plasmons,” _J. Opt. Soc. Am. B_ , vol. 27, no. 9, pp. 1828–1832, Sep 2010\.
  * [53] S. A. Maier, _Surface Plasmon Polaritons at Metal / Insulator Interfaces_. Springer US, 2007, p. 21–37.
  * [54] Y.-R. Fang and X.-R. Tian, “Resonant surface plasmons of a metal nanosphere treated as propagating surface plasmons,” _Chinese Physics B_ , vol. 27, no. 6, p. 067302, jun 2018.
  * [55] A. Taflove and S. C. Hagness, _Computational electrodynamics: the finite-difference time-domain method_ , 3rd ed. Norwood: Artech House, 2005.
  * [56] C. Schinke, C. Peest, J. Schmidt, R. Brendel, M. R. Vogt, I. Kröger, S. Winter, A. Schirmacher, S. Lim, H. Nguyen, and D. Macdonald, “Uncertainty analysis for the coefficient of band-to-band absorption of crystalline silicon,” _AIP Advances_ , vol. 5, p. 067168, 06 2015.
  * [57] P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” _Phys. Rev. B_ , vol. 6, pp. 4370–4379, Dec 1972.
  * [58] G. Mie, “Beiträge zur optik trüber medien, speziell kolloidaler metallösungen,” _Annalen der Physik_ , vol. 330, no. 3, p. 377–445, Jan. 1908.
  * [59] S. Prahl, “miepython: Pure python implementation of Mie scattering,” Aug. 2023\. [Online]. Available: https://doi.org/10.5281/zenodo.8218010
  * [60] ITMO University. (2024) Mie calculator. Accessed on May 22, 2024. [Online]. Available: https://physics.itmo.ru/en/mie#/spectrum
  * [61] K. Svanberg, “The method of moving asymptotes—a new method for structural optimization,” _International Journal for Numerical Methods in Engineering_ , vol. 24, no. 2, pp. 359–373, Feb. 1987.
  * [62] A. Vaccari, A. Calà Lesina, L. Cristoforetti, and R. Pontalti, “Parallel implementation of a 3D subgridding FDTD algorithm for large simulations,” _Progress In Electromagnetics Research_ , vol. 120, pp. 263–292, 01 2011.
  * [63] C. Yu, C. Brakensiek, H. Schumacher, and K. Crane, “Repulsive surfaces,” _ACM Trans. Graph._ , vol. 40, no. 6, dec 2021.
  * [64] B. E. Saleh and M. C. Teich, _Fundamentals of photonics_. John Wiley & Sons, 2019.


Practical Explicitly Invertible Approximation to 4 Decimals of Normal
Cumulative Distribution Function Modifying Winitzki’s Approximation of erf

Alessandro Soranzo∗, Emanuela Epure⋆

∗ Dipartimento di Matematica e Geoscienze – University of Trieste – Trieste –
Italy

email<EMAIL_ADDRESS>
⋆ Esteco S.p.A. (Optimization Software) – Area Science Park – Trieste – Italy

email<EMAIL_ADDRESS>
Abstract. We give a new explicitly invertible approximation of the normal
cumulative distribution function:
$\Phi(x)\simeq\frac{1}{2}+\frac{1}{2}\sqrt{1-{e}^{-x^{2}\frac{17+{x}^{2}}{26.694+2x^{2}}}}$,
$\forall x\geq 0$, with absolute error $<4.00\cdot 10^{-5}$, absolute value of
the relative error $<4.53\cdot 10^{-5}$, which, beeing designed essentially
for practical use, is much simpler than a previously published formula and,
though less precise, still reaches 4 decimals of precision, and has a
complexity essentially comparable with that of the approximation of the normal
cumulative distribution function $\Phi(x)$ immediatly derived from Winitzki’s
approximation of erf$(x)$, reducing about $36\%$ the absolute error and about
$28\%$ the relative error with respect to that, overcoming the threshold of 4
decimals of precision.

2010 Mathematics Subject Classification: 33B20 , 33F05 , 65D20 , 97N50.

Keywords: normal cdf, $\Phi$, error function, erf, Winitzki, normal quantile,
probit, erf-1, approximation, non-linear fitting.

This paper is devoted to approximate some special functions, in particular the
normal cumulative distribution function.

Though computers now allow to compute them with arbitrary precision, such
approximations are still worth for several reasons, including to catch the
soul of the considered functions, allowing to understand at a glance their
behaviour. Let’s add that, despite technologic progress, those functions – of
wide practical use – not always are available on pocket calculators. In
practice, the ancient numerical tables are still widely used; but they give
approximations only for some values and, if linearly interpolating to
approximate intermediate values, both precision and simplicity of use reduce.
And surely by computers you may obtain graphs of those functions, but for a
mathematician the meaning content of formuls is greater.

Furthermore, here we produce an explicitly invertible (and, in fact, simply)
approximation, which allow to keep coherence working contemporarily with the
considered function and its inverse.

For the [6] [8] special functions

$\Phi(x):=\int_{-\infty}^{x}\frac{1}{\sqrt{2\pi}}e^{\frac{-t^{2}}{2}}dt,$ (1)

and the related error function $erf(x)=2\Phi(x\sqrt{2})-1$ $\big{(}\forall
x\in I\\!\\!R\big{)}$, there are several approximations; in particular see
classical [1] [3] [4] and recent [2] [5] [7] [9]; approximations for $x\geq 0$
are sufficient because of the symmetry formula $\Phi(-x)=1-\Phi(x)$
$\big{(}\forall x\in I\\!\\!R\big{)}$.

Restricting now our attention only to those approximations which are simply
explicitly invertible – in the sense, explicitly invertible without requiring
to solve cubic or quartic equations – the most precise appears to be [5]

$\Phi(x)\cong\frac{1}{2}+\frac{1}{2}\sqrt{1-{e}^{\frac{-1.2735457x^{2}-0.0743968x^{4}}{2+0.1480931x^{2}+0.0002580x^{4}}}}\quad\left\\{\begin{array}[]{ll}|\varepsilon(x)|<1.14\cdot
10^{-5}\\\ |\varepsilon_{r}(x)|<1.78\cdot 10^{-5}\end{array}\forall x\geq
0\right.$ (2)

which is an improvement preserving (despite the adding of the quartic
monomial) the simple explicit invertibility (essential solving a biquadratic
equation after obvious substitutions) of this approximation of $\Phi$

$\Phi(x)\simeq\frac{1}{2}+\frac{1}{2}\sqrt{1-{e}^{-\frac{{x}^{2}\,\left(\frac{4}{\pi}+0.0735\,{x}^{2}\right)}{2\,\left(1+0.0735\,{x}^{2}\right)}}}\quad\left\\{\begin{array}[]{ll}|\varepsilon(x)|<6.21\cdot
10^{-5}\\\ |\varepsilon_{r}(x)|<6.30\cdot 10^{-5}\end{array}\forall x\geq
0\right.$ (3)

immediately derived by
$\Phi(x)=\frac{1}{2}+\frac{1}{2}erf\Big{(}\frac{x}{\sqrt{2}}\Big{)}$
$\big{(}\forall x\geq 0\big{)}$ from this [7] Winitzki’s Approximation of erf

$erf(x)\cong\sqrt{1-e^{-x^{2}\frac{\frac{4}{\pi}+0.147x^{2}}{1+0.147x^{2}}}}\quad\left\\{\begin{array}[]{ll}|\varepsilon(x)|<1.25\cdot
10^{-4}\\\ |\varepsilon_{r}(x)|<1.28\cdot 10^{-4}\end{array}\forall x\geq
0\right..$ (4)

In this note we give this new (simply) explicitly invertible approximation of
the normal cumulative distribution function

$\Phi(x)\simeq\frac{1}{2}+\frac{1}{2}\sqrt{1-{e}^{-x^{2}\frac{17+{x}^{2}}{26.694+2x^{2}}}}\quad\left\\{\begin{array}[]{ll}|\varepsilon(x)|<4.00\cdot
10^{-5}\\\ |\varepsilon_{r}(x)|<4.53\cdot 10^{-5}\end{array}\forall x\geq
0\right.$
---
(5)

which, beeing designed essentially for practical use,

$\bullet$ is much simpler than (2) and, though less precise, still reaches 4
decimals of precision;

$\bullet$ has a complexity essentially comparable with that of (3) reducing
about $36\%$ the absolute error and about $28\%$ the relative error with
respect to that, overcoming the threshold of 4 decimals of precision.

Instead, the corresponding approximation of erf is not so worth, because,
though reduces about $36\%$ the absolute error of (4), it remains with the
precision of 3 decimals, and furthermore the absolute value of the relative
error $|\varepsilon_{r}(x)|<1.79\cdot 10^{-4}$ is quite greater than in (4).

See below the graphs (made by Mathematica ®) of the approximation, of the
absolute error and of the absolute value of the relative error for $0\leq
x\leq 7$. For $x\geq 7$ the trivial approximation $\Phi(x)\simeq 1$ has
absolute error and absolute value of the relative error highly less than
$4\cdot 10^{-5}$ and $4.53\cdot 10^{-5}$ respectively. (Nevertheless, if
interested in a formal proof of the majorization $|\varepsilon(x)|<4\cdot
10^{-5}$ for $x\geq 7$, you may follow [5]).

$\begin{array}[]{cc}\includegraphics[width=158.99377pt,height=108.405pt]{Wini26x694xGraph.pdf}&\includegraphics[width=158.99377pt,height=108.405pt]{Wini26x694x0x00004-AbsErrZoom2.pdf}\\\
\rm{Fig.\>1.\>The\>new\>approximation\>(\ref{eq:ourPhiSimple})}&\rm{Fig.\>4.\>Second\>zoom\>of\>Fig.\>2.}\\\
\includegraphics[width=158.99377pt,height=108.405pt]{WiniArxivSecondAbsErr.pdf}&\includegraphics[width=158.99377pt,height=108.405pt]{Wini26x694xRelErr.pdf}\\\
\rm{Fig.\>2.\>Absolute\>error}&\rm{Fig.\>5.\>Absolute\>value\>of\>relative\>error}\\\
\includegraphics[width=158.99377pt,height=108.405pt]{Wini26x694x0x00004-AbsErrZoom1.pdf}&\includegraphics[width=158.99377pt,height=108.405pt]{Wini26x694x0x0000453-RelErrZoom.pdf}\\\
\rm{Fig.\>3.\>First\>zoom\>of\>Fig.\>2}&\rm{Fig.\>6.\>Zoom\>of\>Fig.\>5.\par}\end{array}$

## References

  * [1] Abramowitz, M. Stegun, I.A. (Eds.), (1972). Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th printing New York: Dover, pagg. 932 and 299.
  * [2] Dyer, S.A. Dyer, J.S. (2007). Approximations to error function, Instrumentation $\&$ Measurement Magazine, IEEE 10, no.6:45 – 48.
  * [3] Hart, J.F. et al. (1968). Computer Approximations, SIAM series in applied mathematics, John Wiley & Sons, Inc., New York - London - Sydney, 140, pagg. 288 – 289.
  * [4] Johnson, N. Kotz, S. Balakrishnan, N. (1994). Continuous Univariate Distributions, Vol. 1, 2nd ed. Boston, MA: Houghton Mifflin.
  * [5] Soranzo, A. Epure, E. (2012). Simply Explicitly Invertible Approximations to 4 Decimals of Error Function and Normal Cumulative Distribution Function, www.intellectualarchive.com (selecting Mathematics) and http://arxiv.org/abs/1201.1320v1.
  * [6] Wikipedia. http://en.wikipedia.org/wiki/Error$\\_$function (read 2011, September)
  * [7] Winitzki, S. A handy approximation for the error function and its inverse, http://docs.google.com/viewer?a=v$\&$pid=sites$\&$srcid=ZGVmYXVsdGRv bWFpbnx3aW5pdHpraXxneDoxYTUzZTEzNWQwZjZlOWY2 http://sites.google.com/site/winitzki http://sites.google.com/site/winitzki/sergei-winitzkis-files (2008) (read 2011, December)
  * [8] Wolfram Research, Inc. http://functions.wolfram.com/GammaBetaErf/Erf/10/01/ (read 2011, September)
  * [9] Zogheib, B. Hlynka, M. (2009). Approximations of the Standard Normal Distribution, University of Windsor, Dept. of Mathematics and Statistics. (Legally) available at http://web2.uwindsor.ca/math/hlynka/zogheibhlynka.pdf (read 2012, March).

Remark. This text has been explicitly placed by the Authors in the Creative
Commons Public Domain.