1. Introduction

The last decade has witnessed an exponential increase of research in nanophotonics. Plasmonic and dielectric nanostructured materials, such as metasurfaces and metamaterials, have been developed to engineer the properties of light beyond what is allowed by bulk optical devices [1,2], thus leading to revolutionary solutions for beam structuring [3], colouring [4], biosensing [5], nanomedicine [6], tunable beam steering [7], and nonlinear generation [8,9], to name a few. Advances in nanofabrication technologies have enabled such technologies by allowing unprecedented design complexity at the nanoscale [5,10]. The opportunities offered by theoretical research and nanofabrication facilities raise the need for new methods to efficiently design and optimize such nanophotonic systems [11,12].

Advances in computing capabilities and numerical methods, such as the finite-difference time-domain (FDTD) [13] and finite-element methods (FEM) [14], empowered the design of nanophotonic devices while shortening their design cycle. Conventionally, a design cycle starts from a given layout. The layout is then parameterized and various techniques are employed to explore the parameters’ space. This offers opportunities to find new designs with improved performance or that satisfy multi-objectives. However, exploring large parameter spaces raises computational challenges. Techniques such as parameters sweep or stochastic optimization methods, e.g., genetic algorithms, are computationally intractable for large design spaces and only suitable to handle problems with few design parameters [15]. Deep learning algorithms are also an option, but require large data sets for training, and still not well established [16–18].

Topology optimization (TopOpt) is a robust inverse design approach [19,20]. It was initially introduced to optimize mechanical structures [21], then it has been successfully extended to various engineering disciplines, including acoustics [22], fluids [23], and electromagnetics [24–28]. Typically, TopOpt problems are solved using gradient-based optimization methods, where the gradient of the objective function is computed using efficient methods such as the adjoint-field method [29,30]. Optimization problems that include millions or even billions of design variables have enabled novel conceptual designs [31]. In electromagnetics, TopOpt was used to optimize non-dispersive dielectric devices in the microwave and optical regimes [32–35]. In addition, it was used to design plasmonic antennas using frequency-domain methods [36–38]. Christiansen et al. [37] proposed a non-linear interpolation scheme that was successful to enable TopOpt of plasmonic antennas near their surface plasma frequency using the FEM method.

Structural perturbations or changes in material properties, caused for example by fabrication tolerance or temperature variations, raise the demand to account for the broadband performance of optical components [39]. At optical and near-infrared wavelengths, various materials exhibit dispersion, that can be modelled via Drude, Lorentz, or critical-points functions [13,40,41]. Below their plasma frequency, the Drude model is commonly used to describe the dispersive properties of metals such as silver, gold, and aluminum [42]. The model also describes the permittivity of conductive oxides, such as indium tin oxide (ITO) [43], and conductive polymers [44], including their dielectric and epsilon-near-zero regions.

In this paper, we introduce TopOpt of dispersive optical materials in the time-domain and aim at broadband designs. We base our algorithm on the FDTD method [13] and the Drude model to describe the material dispersion. To the best of our knowledge, this is the first time the FDTD method is used for TopOpt of plasmonic devices in the time domain. To present the algorithm, we conduct the optimization of plasmonic nanoantennas, with the goal to maximize the electric energy in a specified region by finding the distribution of the electric permittivity in a design domain. To evaluate the gradient of the objective function efficiently, we employ the adjoint-field method and provide sensitivity analysis based on Maxwell’s equations in their first-order form. During the optimization, the interpolation of the design material between dielectric and metal enables the surface plasmon frequency $\omega _\text {sp}$ to fall within the frequency band of interest [45,46]. Plasmonic effects, such as field localization and enhancement, are amplified close to $\omega _\text {sp}$. This leads to hot-spots within the optimization domain that prevent the algorithm from converging to well performing designs. To overcome these convergence issues, we exploit the conductivity term in Maxwell’s equations to introduce an artificial damping. This damping counteracts the high-field localization during the optimization process and enables the algorithm to converge to good designs. The developed algorithm is demonstrated by optimizing 2D and 3D silver nanoantennas operating near plasma and infrared frequencies. In all cases, the algorithm produces novel designs with outstanding performance, thus demonstrating its use for the automatic design of dispersive optical nanostructures and metamaterials.

2. Optimization problem setup

In this section, we present the setup of the optimization problem that will be used to inverse design 2D and 3D nanostructures. The computational domain $\Omega$ consists of $\Omega _g \cup \Omega _{d} \cup \Omega _{s} \cup \Omega _\text {PML}$, as shown in Fig. 1(a). The domain $\Omega _g$ is an observation region where the electric energy is to be maximized. We assume that $\Omega _g$ is a dispersionless dielectric medium with relative permittivity $\varepsilon _g$ and has an area $w_g\!\times \! h_g$. Inside the design domain $\Omega _{d}\!=\!w_d\times h_d$, we aim to distribute a dispersive material to form the nanoantenna structure. We consider materials with relative permittivity described by the Drude model:

 

Fig. 1. (a) $\Omega _g$ is the observation domain, $\Omega _{d}$ is the design domain, $\Omega _{s}$ is the background medium, and $\Omega _\text {PML}$ is a perfect matched layer. The boundary $\Gamma$ is used for plane-wave injection. (b) Complex permittivity of silver [47], and fitting via the Drude model. (c) Frequency and (d) time domain plots of a truncated sinc signal with a bandwidth of 20% at half-maximum and modulating a carrier signal with a frequency of 600 THz.

Download Full Size | PDF

In this paper, we consider silver in the wavelength range $350$–$1000$ nm, where we fit its measured complex permittivity [47], via a Drude model with parameters $\varepsilon _\infty = 4.469$, $\omega _p \!=\! 1.426\times 10^{16}$ rad/s, and $\gamma _p \!=\! 4.571\times 10^{13}$ rad/s, as shown in Fig. 1(b). The background space $\Omega _{s}$ hosts the design domain $\Omega _d$, and has a constant relative permittivity $\varepsilon _s$. For simplicity, we use $\varepsilon _s=\varepsilon _g=1$, but deviation from these values is also possible. The computational domain is terminated by perfectly matched layers $\Omega _\text {PML}$. In our analysis, we use the total-field scattered-field formulation [13] to inject a plane-wave through the boundary $\Gamma$, which is located in $\Omega _{s}$. A similar setup was previously used to optimize TM nanoantennas with frequency-domain methods [36,37].

To optimize in the time-domain, the spectral content of the excitation signal determines the desired bandwidth of the structure under optimization. Ideally, we want to excite the system with a signal that has a rectangular spectrum (see Fig. 1(c)). Unfortunately, this corresponds to a sinc signal of infinite duration. Thus, to keep the simulation time reasonable, the sinc signal is truncated to a few lobes, as shown in Fig. 1(d). Such signal modulates a carrier with a frequency corresponding to the center of the spectral window of interest. In addition, we use a Hanning window to reduce the ripples in the excitation spectrum.

We formulate the conceptual optimization problem

2.1 Density-based interpolation and convergence issues

In density-based TopOpt, we use the permittivity function of the material as our design variable. For each point $x$ in $\Omega _{d}$ we want the permittivity at that point to be either that of the Drude material $\varepsilon ^{Drude} (\omega )$ or that of the background space $\varepsilon _s$. In order to interpolate between the background space and the design material, a density variable $\rho _i$ is introduced to describe the material at each edges of the computational grid in $\Omega _d$. The vector $\boldsymbol {\rho } = [\rho _1\, \rho _2\, \ldots \rho _i\, \ldots \rho _M]$ is used to hold the $M$ design variables of the optimization problem. Since we aim to use gradient-based methods to solve the TopOpt problem, the entries of the design vector are allowed to attain values between $0$ and $1$ during the optimization process. However, to obtain a manufacturable design, the final density vector must hold only the binary values $0$ or $1$. We map each design variable $\rho _i$ to the physical material parameters using:

 

Fig. 2. Interpolating between free-space permittivity ($\rho \!=\! 0$) and silver ($\rho \!=\! 1$). (a) Real and (b) Imaginary part of $\varepsilon (\omega,\rho )$ when $\sigma _\text {max} \!=\!0$ S/m. (c) Imaginary part of $\varepsilon (\omega,\rho )$ when $\sigma _\text {max} \!=\!5 \!\times \!10^{5}$ S/m (the real part is the same as (a)). Dispersion diagrams of the surface plasmon polaritons when interpolating between air and silver using (d)-(f) $\sigma _\text {max} \!=\!0$ S/m, (g)-(i) $\sigma _\text {max} \!=\!5 \!\times \!10^{5}$ S/m, and (j)-(l) $\sigma _\text {max} \!=\!5 \!\times \!10^{6}$ S/m.

Download Full Size | PDF

The presence of metal and dielectric in the design domain leads to localization and enhancement of the electric field due to surface plasmon modes arising during the optimization process. This hinders the convergence of the objective function. To understand the reason for such convergence problem, we use, as an illustrative model, the dispersion relation of surface plasmon polaritons (SSPs) at a flat interface between a dielectric $\varepsilon _s$ and our design material $\varepsilon (\omega,\rho )$ [45,46]:

With the aim to reduce the peaks in the dispersion diagrams at intermediate values of $\rho$, we estimate the value of $\sigma _\text {max}$ based on a parameter sweep. Plots for our chosen $\sigma _\text {max}$ value (i.e., $\sigma _\text {max} = 5\times 10^{5}$ S/m) are shown in Figs. 2(g)-(i). A too-large value of $\sigma _\text {max}$ (e.g., $\sigma _\text {max} \!=\! 5\times 10^{6}$ S/m) leads to an intermediate material with conductivity higher than silver, as shown in Fig. 2(j)-(l). To conclude, a too-small value of $\sigma _\text {max}$ is not enough to introduce sufficient damping to counteract the plasmonic effects, while a too-large value prevents the algorithm from converging to black and white designs since the structure becomes less lossy.

3. Methods

3.1 Governing equations and sensitivity analysis

Inside the analysis domain, that we assume source-free and non-magnetic with permeability $\mu _0$, the time-dependent electric field $\boldsymbol {\mathcal {E}}$ and magnetic field $\boldsymbol {\mathcal {H}}$ are governed by Maxwell’s equations

In our case, we use the design permittivity model

We solve the design problem by using a gradient-based optimization method. To do so, we need derivatives of the objective function with respect to the design variables $\varepsilon _\infty$, $\psi$, and $\sigma$. We use the adjoint-field method to derive expressions for such derivatives. Assuming that the design variables are perturbed by $\delta \varepsilon _\infty$, $\delta \psi$, and $\delta \sigma$, the corresponding first variations of $W$ are

To find explicit expressions for Eqs. (12), we use the system governing Eqs. (10) and employ the adjoint-field method [29,30]. In the following, we derive an explicit relation for $\delta _\psi \! W$. To simplify the derivation, we drop the domain $\Omega _\text {PML}$ and consider $\Gamma$ as the external boundary of the analysis domain $\Omega$, see Fig. 1(a). In addition, we drop the differential $dt$ and $d\Omega$ since they can be inferred from the limits of the integrals. We consider the initial-boundary-value-problem,

We define the adjoint fields $\boldsymbol {\mathcal {E}}^{*}$, $\boldsymbol {\mathcal {J}}^{*}$, and $\boldsymbol {\mathcal {H}}^{*}$. We perform the scalar product of Eqs. (14a),  (14b), and  (14c) with $\boldsymbol {\mathcal {E}}^{*}$, $\frac {\boldsymbol {\mathcal {J}}^{*}}{\varepsilon _0 \psi }$, and $\boldsymbol {\mathcal {H}}^{*}$, respectively. We add the result of multiplication, integrating over the whole analysis domain $\Omega$ and the observation interval $(0,T)$, and applying integration by parts, we obtain

We assume that the adjoint fields satisfy the terminal conditions $\boldsymbol {\mathcal {E}}^{*}\!=\!\boldsymbol {\mathcal {J}}^{*}\!=\! \boldsymbol {\mathcal {H}}^{*}\!=\! \boldsymbol {0}$ at $t\!=\!T$. By arranging the terms in Eq. (15), utilizing Eq. (14d), and adding and subtracting $\delta _\psi W$, we obtain

Now, if we require

Expression (18) is the directional derivative of $W$ when $\psi$ is perturbed by $\delta \psi$. The gradient of $W$ with respect to $\psi$ can be identified as the integral kernel

The adjoint system  (17) is a terminal-value-problem which, by changing the time variable (i.e., $t\!=\!T-\tau$) and the sign of the magnetic field $\boldsymbol {\mathcal {H}}^{*}$ (i.e., to preserve the direction of the Poynting vector), can be written as

By differentiating system  (13) with respect to $\varepsilon _\infty$ and $\sigma$, and following similar procedures as before, we obtain the same adjoint system  (20) and the following gradient expressions

Thus, to evaluate the gradient of the objective function, we solve the forward system  (13) and the adjoint system  (20), then we use Eqs. (21), (22), and (23) to form the full gradient components.

3.2 Numerical treatments and optimization algorithm

We solve numerically the system of governing equations, discussed in the previous section, using the FDTD method [13]. We adopt the auxiliary differential equation approach to implement the Drude model in the FDTD method, and the uniaxial perfectly matched layer (UPML) is used to simulate the open-space radiation boundary condition [13,40]. The computational domain is discretized into uniform square (2D) or cubical (3D) Yee cells with spatial steps $\Delta x \!=\! \Delta y \!=\! \Delta z$ in all Cartesian directions.

In the forward system, we use the FDTD method and discretize the electric field at full-time indices and the magnetic field at half-time indices. In the adjoint system, however, the FDTD discretization of the electric- and magnetic fields are performed at half-time indices and full-time indices, respectively [25,30]. The discretized objective function is

To avoid mesh-dependency or self-penalization issues, in density-based TopOpt it is common to filter the design variables [48–52]. That is, instead of using $\rho _i$ in Eq. (5), we replace it with $\tilde {\rho }_i$, where the filtered design vector $\boldsymbol {\tilde {\rho }}$ is obtained through the mapping

In this work, we use an open-close, nonlinear filter operator $\mathcal {F}(\cdot )$ that consists of a cascade of four $fW$-mean filters [52]. The filter has two tuning parameters that determine its size and the level of nonlinearity. Here, we fix the filter size to a constant value of $5\Delta x$ and only employ the nonlinearity parameter to smoothly decrease the level of greyness in the design during the optimization process [52,53]. Using the chain rule, the derivative of the discrete objective function with respect to the design variable $\rho _i$ is evaluated by

We compared the derivatives computed using Eq. (27) against those evaluated by finite differences. The comparison showed more than $4$ digits match in precision between the two methods. We write the discrete version of the optimization problem as

 

Fig. 3. Flowchart of the optimization algorithm.

Download Full Size | PDF

4. Results

In this section, we demonstrate our method through some design examples of 2D and 3D plasmonic nanostructures. In 2D, only the case of TM excitation is presented. In fact, TE excitation does not lead to plasmonic effects, and the optimization algorithm shows regular convergence to traditional focusing mirror structures. To enable fast simulation, we implement the FDTD method to execute on graphics processing units (GPUs). The computations are executed on nodes equipped with NVidia V100 GPUs and 64 GB of memory. Based on the problem size, one call to the Maxwell solver uses a simulation time between less than a minute and few minutes.

4.1 2D, TM nanoantennas

We chose a design domain $\Omega _d$ with dimensions $w_d \!=\! h_d\!=\!100\,\text {nm}$, and the observation domain $\Omega _g$ is centered within $\Omega _d$ and has dimensions $w_g \!=\! h_g\!=\!10 \,\text {nm}$, see Fig. 1(a). We use a space-step $\Delta x=0.5\,\text {nm}$, and a time-step $\Delta t$ satisfying the Courant stability criterion. The simulation domain is truncated by 15 UPML cells placed 30 cells away from $\Omega _d$. The excitation is an $\mathcal {H}_{x}$ polarized plane-wave propagating towards the positive $y$ axis. The excitation spectrum has a bandwidth of 20% at half-maximum and is centered at $413$ nm, as shown in Fig. 4(c) along with the results of optimization. The design variables are mapped to the in-plane permittivity components associated with the Yee edges where the electric field components $\mathcal {E}_y$ and $\mathcal {E}_z$ are located. Excluding the observation region, the design domain includes 79 560 design edges.

 

Fig. 4. (a) Progress of the objective function and some samples showing the development of AntTM1 (see Visualization 1). (b) AntTM1 topology. (c) Average field enhancement in $\Omega _g$ together with the excitation spectrum. (d)-(j) Field distribution at $\lambda =375$, $410$, $430$, and $460$ nm.

Download Full Size | PDF

Figure 4(a) shows the progress of the normalized objective function $\tilde {W}$ versus the iteration numbers. We start the algorithm with a uniform initial distribution $\rho _i \!= \!0.5$ for all design variables. Included in the same figure are some snapshots to show the development of the design. The black color indicates silver ($\rho \! =\!1$) and the white color indicates air ($\rho \!=\!0$). Notably, the main topology of the antenna evolves after only a few tens of iterations. However, most of the late iterations are used to remove the intermediate material and form crisp boundaries, while the objective function keeps increasing monotonically. The design algorithm converged after $308$ iterations to a design with a grayness level $\zeta < 0.2$%. We threshold this design around $\rho _\text {th}=0.5$ and show the final design, AntTM1, in Fig. 4(b). On the side facing the incident wave, the topology of AntTM1 developed as a flared horn, backed by a small cavity region. On the other side, we see two slightly tilted vertical arms. This optimized topology shares similarities with results reported in the literature using the FEM method [37].

Fig. 4(c) shows the average field enhancement of AntTM1 at $\Omega _g$, which correlates well with the excitation spectrum shown in the same figure. We cross-validate our computations with the commercial software package Ansys Lumerical FDTD [55]. Slight differences between the two computations are attributed to differences in geometry descriptions. Inside $\Omega _g$ and within the excitation window, AntTM1 exhibits more than $20$-fold field enhancement. The peak of the performance, $(\overline {|E|/|{E}_\text {in}|})_{\Omega _g} \!\!=\! 27.8$, occurs at the wavelength $\lambda \!=\!430$ nm, which resides at long wavelengths in the excitation window. Figs. 4(d)-(j) show the electric field distribution of AntTM1 at four wavelengths, marked in Fig. 4(c). The electric field is maximum at $\Omega _g$, marked by the box at the center. However, we observe a field localization and enhancement at the device’s boundaries for short wavelengths, which justifies the decrease of the energy enhancement in $\Omega _g$ at short wavelengths. We attempt to improve the optimization results further by exploring the following investigations.

4.1.1 Effect of the design domain size

The design obtained in Fig. 4 hits the boundary of the design domain, which suggests the need for a larger design space. We double the size of the design domain to $200\!\times \!200$ nm$^{2}$, and we solve the optimization problem, which now includes $319\,960$ design variables. The algorithm converged in $297$ iterations to the topology shown in Fig. 5(a). The new design, AntTM2, has more topological features compared to AntTM1. We observe an additional vertical arm that evolved in the rear-side of the device, and the arm around $\Omega _g$ appears straight. The performance of AntTM2 has improved at short wavelengths, and it attains a nearly flat response within the excitation window, as shown in Fig. 5(b). Fig. 5(c) shows the field distribution at $\lambda \!=\!440$ nm.

 

Fig. 5. Topology, average field enhancement in $\Omega _g$, and field distribution at wavelength of maximum performance for (a)-(c) AntTM2 (see Visualization 2), (d)-(f) AntTM3 (see Visualization 3), and (g)-(i) AntTM4 (see Visualization 4), respectively. The colored insets in (d) shows the design at iteration #1.

Download Full Size | PDF

4.1.2 Effect of a fixed gap geometry

Maximizing the energy in $\Omega _g$ suggests that the objective function is more sensitive to design variables close to $\Omega _g$ compared to those away from it. We attempt to relax such a non-uniform sensitivity distribution and investigate its effect on the results. We fix the geometry region below and above $\Omega _g$ to silver with the same area as $\Omega _g$. Figs. 5(d)-(f) shows AntTM3 optimized over a design domain with size $100\times 100$ nm$^{2}$. For AntTM3, fixing the area around the gap boosts the average field enhancement to a maximum value of $40.6$ at the wavelength $435$ nm; AntTM1 has a maximum value of $27.8$. The new nanoantenna exhibits better performance at long wavelengths, however, the performance at short wavelengths stays essentially the same as in the previous cases. These results indicate the challenges that plasmonic effects pose on optimizing nanoantennas near the surface plasmon frequency.

4.1.3 Wideband optimization

We use the previous setup and attempt to optimize over a wider spectrum covering the wavelength window $375$–$900$ nm. That is, we use an excitation signal with a half-maximum bandwidth of $82$% centered around $637.5$ nm. Figs. 5(g)-(i) show the topology, the average field enhancement, and the field distribution at the wavelength of maximum performance for the new design. The new nanoantenna, AntTM4, shows a wideband performance. An average field enhancement above $10$-fold is possible to achieve. Moreover, the average field enhancement hits a maximum of $41.5$ at the wavelength $810$ nm. Here, we still observe the performance bias of the optimized nanoantennas towards long wavelengths.

4.2 3D antennas

We use the developed method to optimize antennas in a 3D setup. We extend the problem model (Fig. 1) to include a $30$ nm thickness in the $x$-direction. The design space has a volume $\Omega _d=200\!\times \!200\!\times \!30$ nm${}^{3}$ with $\Omega _g \!=\!12\!\times \!12\!\times \!30$ nm${}^{3}$, and we use $\Delta x\!=\! 2$ nm. Similar to the 2D case, here we fix the geometry region below and above $\Omega _g$ to silver with the same area as $\Omega _g$. The discretized design domain includes $469\,341$ design variables associated with its interior edges. We impose symmetry along the $x$-axis to enable antennas producible by current technologies. That is, we optimize 3D antennas and aim for planar structures. The excitation is an $\mathcal {E}_z$ polarized plane-wave propagating in the positive $x$ axis. We use the same setup and solve the optimization problem for three different wavelength excitation windows. The first and second excitation spectra, shown in Fig. 6(a), have a half-maximum bandwidth of $20$% centered around $413$ nm and $513$ nm, respectively. The third excitation spectrum, shown in Fig. 6(b), has a half-maximum bandwidth covering the spectral window $375$–$900$ nm. Each excitation spectrum results in a different topology which we name Ant3D1, Ant3D2, and Ant3D3, see Figs. 6(c)-(e). Interestingly, we observe the increase of the figure-of-eight void area of the three nanoantennas, around $\Omega _g$, as the excitation spectrum includes long wavelengths.

 

Fig. 6. Average field enhancements inside $\Omega _g$ together with the excitation spectrum for (a) Ant3D1 and Ant3D2, (b) Ant3D3. Topologies of (c) Ant3D1 (see Visualization 5), (d) Ant3D2 (see Visualization 6), and (e) Ant3D3 (see Visualization 7). Field enhancements at the middle layer of (f)-(i) Ant3D1 at wavelengths marked in Fig. 6(b).

Download Full Size | PDF

Fig. 6(a) shows the average field enhancement of Ant3D1 and Ant3D2, and Fig. 6(b) shows the performance of Ant3D3. Inside $\Omega _g$, the nanoantennas Ant3D1, Ant3D2, and Ant3D3 exhibit a maximum average field enhancement of $50.5$, $52.8$, and $46.7$ at the wavelength $445$ nm, $545$ nm, and $615$ nm, respectively. Ant3D3 exhibits another peak of $74.0$ at $970$ nm, which resides slightly outside the intended excitation spectrum. The optimized structures tend to exhibit a better field enhancement at long wavelengths. Figs. 6(f)-(i) show the field distribution of Ant3D3 at wavelengths marked in Figs. 6(b). The optimized structures are capable to maximize the electric energy at the observation domain $\Omega _g$. At short wavelengths, however, we observe high-field localization near nanoantenna’s boundaries, which indicates strong plasmonic effects that are responsible for the decrease in the achieved performance. Further investigations are needed to obtain a balanced performance over the wavelength window of interest.

5. Conclusion

We introduced a density-based TopOpt approach to design plasmonic dispersive nanoantennas. Our approach is based on Maxwell’s equations in the time-domain, and we use the Drude model, which can fit the material dispersion of metals and conductive polymers, as well as epsilon-near-zero materials, such as conductive oxides. For the TM and the 3D setups, the interpolation between metallic and dielectric phases results in high field-localization associated with plasmonic effects, which prevent the algorithm from converging to well-performing designs. Guided by dispersion diagrams of metal-dielectric interfaces, we proposed an artificial damping approach to suppress the field-localization during the optimization process, which enables the algorithm to converge to good designs. For the TE setup, artificial damping is not needed and the algorithm encounters no convergence issues. Various setups for narrowband and wideband optimization are presented, resulting in novel 2D and 3D nanoantenna designs with outstanding performances. Our method opens new opportunities for the automatic design and optimization of dispersive nanophotonic structures with broadband optical response for applications in nonlinear plasmonics, integrated optics, plasmonic colouring and absorbers, epsilon-near-zero photonics, thermoplasmonics, biosensors, and ultrafast optics.