## Abstract

We show that topology optimization (TO) of metallic resonators can lead to ∼10^{2} × improvement in surface-enhanced Raman scattering (SERS) efficiency compared to traditional resonant structures such as bowtie antennas. TO inverse design leads to surprising structures very different from conventional designs, which simultaneously optimize focusing of the incident wave and emission from the Raman dipole. We consider isolated metallic particles as well as more complicated configurations such as periodic surfaces or resonators coupled to dielectric waveguides, and the benefits of TO are even greater in the latter case. Our results are motivated by recent rigorous upper bounds to Raman scattering enhancement, and shed light on the extent to which these bounds are achievable.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

In this paper, we present freeform shape optimization ("topology optimization," TO [1–3]) of metallic nanoparticles for Raman scattering [4–6], and obtain non-intuitive structures $\sim 60\times$ better (in terms of emitted power) than optimized coupled-sphere [7–9] or bowtie [10–14] antennas (Sec. 3.1) for identical separation distance to the Raman molecule. Our current results are a proof-of-concept of Raman TO in 2d systems, and the resulting dramatic enhancements suggest exciting opportunities for future improvements in practical 3d Raman sensing. Stated briefly, Raman scattering consists of an incident wave at a frequency $\omega _1$ interacting with a molecule that subsequently emits electromagnetic radiation at $\omega _2$. A nanostructure can enhance both the incident-wave absorption by a focusing effect as well as the emission by a Purcell effect [15–20]. Figure 1(a) shows a schematic of this process, in which the molecule is surrounded by an unknown arrangement of metal (e.g. silver); TO is used to tailor this arrangement to maximize the Raman scattering (Sec. 2), resulting in surprising structures such as the one shown in Fig. 1(b). In addition to optimizing isolated resonators coupled with incident planewaves (Fig. 1(a)), we also optimize Raman scattering on periodic surfaces as well as resonators coupled with input/output waveguides, as depicted in Fig. 2(b) and 2(c), respectively. We show that it is important to consider the full Raman process combining *both* focusing and emission, as optimizing either emission or focusing alone sacrifices a factor of $\sim 5\times$ (Sec. 3.2). When the Raman shift $\omega _1 - \omega _2$ is more than a few percent (i.e., comparable to the bandwidth of a single plasmonic resonance), we show that this shift must be taken into account in the optimization, or one may sacrifice a factor of $\sim 5\times$ (Sec. 3.5). We find that the huge enhancements observed from TO structures compared to simple geometries are in qualitative agreement with recently discovered theoretical upper bounds to Raman scattering [21]. Quantitatively, for a material with susceptibility $\chi$, the key figure of merit for light-matter interactions is $F=\vert \chi \vert ^2/\mathrm {Im}(\chi )$ [22–24] and the Raman bounds scale as a factor of $F^3 V$ where $V$ is the volume of the scatterer: a factor of $F^2$ maximum focusing enhancement and a factor of $F$ maximum emission enhancement. The enhancement achieved by our TO designs scale roughly linearly with $V$ in agreement with the upper bounds (Sec. 3.3), but they scale roughly with $F^2$ instead of $F^3$ (Sec. 3.4) suggesting, in part, that in practice one cannot simultaneously achieve ideal focusing and emission enhancement.

In conventional Raman spectroscopy, the very small Raman cross-section of most chemicals results in Raman radiation typically on the order of 0.001% of the power of the pump signal [4]. This low efficiency is overcome in surface-enhanced Raman spectroscopy (SERS) by placing the chemicals of interest in the vicinity of a scatterer, typically a surface or collection of nanostructures, which increases the overall signal that can be collected [15–20,25]. This technique has allowed for efficiencies up to 12 orders of magnitude larger than that of traditional Raman spectroscopy, yielding detection levels down to the single molecule [26,27] and opening up applications in the fields of biochemistry, forensics, food safety, threat detection, and medical diagnostics [19,20]. While many different materials and antenna geometries have been used for SERS substrates [28–30], so far the optimization of these geometries has been limited to one or two degrees of freedom in designs such as bowtie antennas [10,31–34]. Comparison with the recently derived upper bound to Raman enhancement showed these geometries to be greatly suboptimal, with performance several orders of magnitude below the bounds [21]. Although it is possible that the bounds can be tightened (e.g. by incorporating additional physical constraints [35]), the fact that current SERS geometries are so far below these new theoretical limits made us wonder if dramatic improvements might be attainable by TO.

In this work we thus take a hitherto unexplored approach in the context of Raman scattering, in which the problem of designing metallic nanostructures to enhance the process is formulated as a mathematical optimization problem and solved using density-based topology optimization [1]. In the design process, we simultaneously optimize the focusing of the incident field and the emission from the molecule (dipole), which is shown to lead to structures with higher performance than only optimizing for one process. In brief, density-based topology optimization operates by introducing a continuous design field to control the physical material distribution, enabling the use of adjoint sensitivity analysis [36] and gradient-based optimization algorithms [37,38] to efficiently solve design problems with potentially billions of design degrees of freedom [39] . Hence, the approach provides near-unlimited design freedom, with a computational complexity dominated by the solution of the Maxwell equations, utilizing mature finite-element techniques [40]. A suite of well-understood tools, developed or matured over the last decades, are used to solve the optimization problem, interpolate material parameters, control design variations and ensure physical admissibility of the design [38,41–44]. In addition, geometric length-scale constraints are employed to ensure that all features of the final designs are above a specified size (which may be chosen to comport with fabrication constraints) [45]. Over the past 20 years, TO has been applied to an increasingly wide range of problems in electromagnetic design [2,3]. Our work on Raman optimization couples multi-frequency focusing and emission problems. Maximizing the emission alone (for a dipole source) corresponds to maximizing the local density of states (LDOS) [46,47], a formulation that has been employed for TO of optical cavities [48,49]. The focusing problem alone is related to lens and antenna design among others, for which TO has also been applied (both to small metallic particles such as those considered here and to much larger structures, such as metalenses) [44,50–54]. Nonlinearly coupled electromagnetic resonances at multiple frequencies were optimized via TO for second harmonic generation [55].

## 2. Model formulation

The Raman scattering phenomenon (sketched in Fig. 1(a)) is modelled as two sequentially-coupled, time-harmonic, classical electromagnetic problems [40,56] in a domain $\boldsymbol {\Omega }$. The first models an in-plane polarized planewave propagating through the domain, interacting with any material distributed inside. The second models the excitation of a Raman molecule in $\boldsymbol {\Omega }$ by the electric field resulting from the first problem. The Raman molecule is modelled as a dipole source, with its dipole moment given by a polarizability tensor multiplied by the electric field, excited by the incident planewave at the position of the dipole [18]. The problem may be stated formally as

The problem of enhancing the power emitted from the Raman molecule at $\lambda _2$, by tailoring the material distribution in the design domain $\boldsymbol {\Omega }_{\mathrm {D}}$, is formulated as a continuous optimization problem and solved using density-based topology optimization. In its basic form the optimization problem may be written as

Density-based topology optimization is used to solve Eqs. (4)–(5) where the physical admissibility of the final material distributions is assured using filtering, thresholding and continuation techniques developed in the mechanics community [42,43]. This allows the enforcement of a minimum length scale (6 nm) on all features in the material distribution using geometric length-scale constraints [45].

We consider the three configurations shown in Fig. 2. The first (Fig. 2(a), Sec. 3.1–Sec. 3.6) is used to benchmark the proposed approach. It consists of a Raman molecule (blue dot) placed in an air background with $\boldsymbol {\Omega }_{\mathrm {D}}$ (gray region) surrounding the molecule. An incident planewave propagates through the domain from left to right (green) and the power emitted through $\Gamma _{{\mathrm {out}}}$ by Raman scattering (red) from the molecule is sought maximized.

The second problem (Fig. 2(b), Sec. 3.7) models out-of-plane Raman scattering and considers a metallic surface with periodic patterning placed in a dielectric medium. The model problem is $x$-periodic with period $a$. It consists of a spatial region containing the design domain with a Raman molecule placed at its centre, truncated from below by a metallic surface. An incident planewave propagates through the domain from top to bottom and the power emitted from the molecule is sought maximized in the direction normal to the metallic surface. The Raman scattering process from the different molecules on the surface is assumed to be incoherent. Therefore a $k$-space integration technique, the array scanning method [57,58], is used to compute the power emitted from each individual Raman molecule situated in the periodic background.

The third model problem (Fig. 2(c), Sec. 3.8) concerns in-plane Raman scattering into a waveguide (output). A Raman molecule is excited by an incident field coupled into the system from a second waveguide (input). The Raman molecule is positioned at the centre of $\boldsymbol {\Omega }_{\mathrm {D}_{\mathrm {Metal}}}$ (dark gray), in which a metallic nanostructure is designed. In the remaining part of the design domain, $\boldsymbol {\Omega }_{\mathrm {D}_{\mathrm {Dielectric}}}$ (light gray), a dielectric structure, aimed at coupling the light into and out of the waveguides, is designed simultaneously. For this problem, the difference in the power emitted by the Raman molecule through $\Gamma _{{\mathrm {out}}}$ and $\Gamma _{{\mathrm {in}}}$ is sought maximized.

The model domains are truncated using different boundary conditions depending on which problem configuration is considered. For the isolated resonators (Fig. 2(a)) and the resonators coupled with the input/output waveguides (Fig. 2(c)), a perfectly matched layer is used to truncate $\boldsymbol {\Omega }$. For the Raman molecules placed on a periodic surface (Fig. 2(b)), Floquet-Bloch boundary conditions are used in plane, a perfectly matched layer truncates the top boundary and a perfect electric conduction condition is imposed on the bottom boundary.

The model problems are implemented in COMSOL Multiphysics [59] and the optimization problem is solved using the Globally Convergent Method of Moving Asymptotes [38] utilizing a maximum of 3 inner iterations per design iteration. Details concerning the numerical simulation, optimization, post-processing, and evaluation processes are given in Appendix A and B.

## 3. Results

This section details eight studies conducted to investigate various aspects of the Raman scattering enhancement (Raman enhancement) problem. Study-specific parameters are given in each of the following subsections, while a general set of parameters used across all studies may be found in Appendix B.

#### 3.1 Comparison to known solutions

To demonstrate the proposed approach, we design a silver nanostructure that enhances Raman scattering from an isolated Raman molecule (Fig. 2(a)) excited at $\lambda _{1} = 532$ nm, assuming a negligible Raman shift ($\lambda _2 = \lambda _1)$; non-negligible shifts are considered in Sec. 3.2 and Sec. 3.5. For context, the achieved enhancement is compared to enhancements generated by two simple parameter-optimized references.

Regarding the reference geometries, it is well known that placing a Raman molecule between two carefully scaled circular metallic cylinders (a 2d analogue of a coupled-sphere antenna) significantly enhances the Raman scattering process [7–9], while placing it between two identical metallic triangular structures (a bowtie antenna) further enhances the process [10–14]. Both these references are parameter-optimized to maximize their performance at the targeted wavelengths. For the coupled-cylinder antenna, the radius of the two discs is optimized over $R_{\mathrm {csa}} \in [10 \ \mathrm { nm}, 50 \ \mathrm { nm}]$, and we find the largest enhancement for $R_{\mathrm {csa}} = 48$ nm. The bowtie antenna is optimized by sweeping the tip-angle over, $\theta _{\mathrm {bta}} \in [10^{\mathrm {o}}, 180^{\mathrm {o}}]$, and the side length over $L_{\mathrm {bta}} \in [20 \ \mathrm { nm}, 100 \ \mathrm { nm}]$, and we find $\theta _{\mathrm {bta}} = 15^{\mathrm {o}}$ and $L_{\mathrm {bta}} = 70$ nm to yield the largest enhancement at $\lambda _{1} = \lambda _{2} = 532$ nm.

Figure 3 shows the enhancement of the emitted power $\big ($Eq. (4)$\big )$ versus wavelength for each of the three structures, relative to the power emitted from the Raman molecule in free space. The coupled-cylinder antenna (black) is seen to achieve a nearly wavelength-independent enhancement of $\approx 3 \cdot 10^2$, while the bowtie antenna (red) achieves a peak enhancement of $\approx 1.3 \cdot 10^4$ at $532$ nm. Interestingly, the intricate geometry of the topology-optimized nanostructure (blue), fully encapsulating the Raman molecule, achieves an enhancement of $\approx 8.0 \cdot 10^5$. This is an increase by a factor of $\sim 60\times$ relative to the bowtie antenna. The TO structure is in some sense a fusion of different features, tailored to enhance either the focusing or emission process, as will be studied more closely in Sec. 3.2. To our knowledge, no metallic structure for Raman enhancement similar to the optimized structure has previously been proposed, nor is such a structure likely to arise from applying any traditional design rules or intuition.

A finding worth noting for this study, is that starting from an initial material configuration of $\xi (\textbf {r}) = 0.0 \ \forall \textbf {r} \in \boldsymbol \Omega _{\mathrm {D}}$ results in the design process converging to a structure not encapsulating the Raman molecule. Importantly, this non-encapsulating design achieves an enhancement which is lower by more than a factor of $\sim 2\times$ compared to the fully encapsulating design in Fig. 3. In general, any structure found with gradient-based TO may only be a local optimum to the design problem. Alternative global optimization formulations are rarely practical or even feasible in large design spaces for non-convex problems [60] and do anyway not guaranty global minima. Gradient-based optimization from different starting points (a "multistart" algorithm) can explore different local optima but in this work, the main goal is to find a structure much better than what could be easily designed by hand. Comparison to theoretical upper bounds is another route to gauging global optimality of TO structures [22,61].

#### 3.2 Raman vs. focusing vs. emission

A question worth asking, as it could potentially reduce the computational effort associated with the design procedure by a factor of two, is whether it is necessary to consider the full two-step Raman scattering process in the design procedure, or if it is possible to achieve similar performance by only considering an energy-focusing $\big ($Eq. (1)$\big )$ or a dipole-emission $\big ($Eq. (3)$\big )$ problem. To answer this question using the proposed approach, we consider the following three cases. A dipole-emission case, an energy-focusing case and a case considering the full two-step process. All cases assume a 50 nm Raman shift, with $\lambda _{1} = 532$ nm and $\lambda _{2} = 582$ nm.

For the dipole-emission case, the optimization problem $\big ($Eqs. (4)–(5)$\big )$ remains unchanged, while the dipole considered in Eq. (3) has its orientation fixed along the $y$ axis and its magnitude kept constant throughout the optimization, i.e. $\textbf {f}_{2} = - \omega _{2}^2 \textbf {I} \langle 1,0\rangle \delta (\textbf {r}-\textbf {r}_{0})$, effectively removing the first model problem $\big ($Eq. (1)$\big )$ as $\textbf {f}_{2}$ no longer depends on $\textbf {E}_{1}$. This corresponds to maximizing the local density of states for the dipole [48]. For the energy-focusing case, the optimization problem is changed to one of solely maximizing $\vert \textbf {E}_{1} \vert ^2$ for an incident planewave at the position of the dipole [44,50] effectively removing the second model problem $\big ($Eq. (3)$\big )$, as the emission process is no longer considered.

Figure 4 reports the enhancement of the three objectives attained by introducing the nanostructures optimized for the emission (red), focusing (black) and two-step Raman scattering (blue) process. The leftmost bar group reports the emission enhancement for a dipole with unit magnitude oriented along the $y$ axis, when placed at the center of each of the three designs. The middle bar group reports the enhancement of $\vert \textbf {E}_{1} \vert ^2$ at the dipole position for each design, while the rightmost bar group reports the enhancement achieved for the full two-step Raman scattering process. From the figure, it is seen that optimizing for a given process yields the best performance for that process. In particular, by explicitly optimizing for the two-step Raman scattering process, an increase in the enhancement by a factor of more than $4.5\times$ is achieved. It is interesting to observe that the design optimized for the full Raman scattering process qualitatively consists of a fusion of geometric features found in the focusing and the emission designs, e.g the closed cavity (emission) and the protruding features on the right side of the design (focusing). This suggests that a combination—and, by extension, a trade-off— between the features is required to achieve the largest Raman enhancement. By evaluating the emission and focusing components of the enhancement bound presented in Ref. [21] for the red and black structures in Fig. 4, we discover the following. The emission enhancement achieved by the red structure ($\approx 10^3$) comes within a factor of four of the bound ($\approx 3.5 \cdot 10^3$), where difference may be explained by the length-scale imposed on the design. The focusing enhancement achieved by the black design ($\approx 10^3$) is, on the other hand, nearly three orders of magnitude from the focusing bound ($\approx 6 \cdot 10^5$). Hence, it has not been possible to reach the focusing bound, the magnitude of the discrepancy suggesting that the bound may not be tight.

#### 3.3 Size scaling

It was recently shown that an upper bound for enhancing the Raman scattering process using metallic nanostructures scales linearly with the volume (area in 2d) of the design region [21]. To investigate if this scaling is found using the proposed design approach, we perform a study considering three different outer radii of $\boldsymbol \Omega _{\mathrm {D}}$. These are $R_{\mathrm {out}} = 50$ nm (black), $R_{\mathrm {out}} = 75$ nm (red) and $R_{\mathrm {out}} = 100$ nm (blue), respectively. All other parameters are kept constant across the three cases and an assumption of a negligible Raman shift is made, i.e. $\lambda _{1} = \lambda _{2} = 532$ nm.

Figure 5 shows the Raman enhancement relative to a molecule in free space versus the area of $\Omega _{\mathrm {D}}$ with the three designs and their performance plotted in black, red and blue and a trend line of slope 1 included to guide the eye. The data shows an approximately linear scaling of the emission enhancement with volume (area) in agreement with the upper bounds derived in Ref. [21], demonstrating that the proposed approach finds the expected volume scaling.

#### 3.4 Material scaling

The upper bound for enhancing the Raman scattering process, derived in Ref. [21], scales cubically with the material-related quantity $F = \vert \chi \vert ^2 / \mathrm {Im}(\chi )$. Here a factor of $F^2$ stems from energy-focusing enhancement while the remaining factor of $F$ stems from dipole-emission enhancement. To investigate this behaviour, we conduct a materials study considering four metals, silver (Ag), gold (Au), copper (Cu) [62], and platinum (Pt) [63]. All parameters, other than the material parameter $\varepsilon (\textbf {r})$, were kept constant across the four cases and a negligible Raman shift was assumed with $\lambda _{1} = \lambda _{2} = 532$ nm.

Figure 6 shows the Raman enhancement obtained by placing the Raman molecule at the center of the optimized structures relative to free space versus $\vert \chi \vert ^2 / \mathrm {Im}(\chi )$ for the Cu- (magenta), Au- (black), Pt- (red) and Ag-design (blue). To guide the eye, two trend lines (gray) with slope 2 and slope 3 are plotted with the data. From a linear fit of the data, an approximate scaling with $F^2$ is observed. Trusting that TO indeed provides highly optimized results, as confirmed by the previous examples, this data suggests that the theoretical bounds may not be tight for two reasons. First, it may not be possible to create a design which achieves ideal focusing- and ideal emission-enhancement simultaneously, but that a trade-off must be made. The observation is supported by the data in Fig. 4, where significantly different geometries are observed when targeting the maximization of either focusing or emission, suggesting that different geometries are required to achieve either the highest attainable focusing or the highest attainable emission. Second, combining the observed $F^2$ scaling in Fig. 6 with the finding in Sec. 3.2 that the focusing bound is not reached when optimizing purely for focusing, leads to the suggestion that the bound might not be tight. In conclusion, we find that the optimized structures (while vastly superior to bowtie antennas) still fall far short of the upper bound [21], especially for large values of $F$. Future analytical work including additional constraints, e.g. coupling the two processes, may lead to a tighter bound.

#### 3.5 Accounting for the Raman shift

Depending on the molecule(s) and energy transitions of interest for a given Raman scattering problem, the Raman shift between the wavelengths $\lambda _{1}$ and $\lambda _2$ will be different. In this section, we investigate the benefits of accounting for the Raman shift in the design process. To this end, three cases are considered (Fig. 7) with Raman shifts of 0 nm (blue), 50 nm (red), and 100 nm (black), respectively. The wavelength of the incident field is kept constant at $\lambda _{1} = 532$ nm and the wavelength of the light emitted by the Raman molecule is adjusted accordingly.

Figure 7 presents the three designs along with the power-emission enhancement attained for each design when evaluated at 0 nm shift (leftmost bar group), 50 nm shift (middle bar group), and 100 nm shift (rightmost bar group). A first observation is that all three structures share the same overall geometrical features and thus in that sense look similar. However, looking more closely at each structure it is clear, that both the size and shape of each feature change across the designs. Furthermore, looking at the differences in enhancement, it is clear that these delicate feature changes are reflected in the performances of the structures. Unsurprisingly, each design exhibits the largest enhancement at the Raman shift it was optimized for. Specifically, an enhancement difference of a factor of $\sim 1.6\times$ is observed for the 50 nm Raman shift while a factor of $\sim 4.5\times$ is observed for the 100 nm Raman shift, clearly illustrating the performance benefit of accounting for the Raman shift as part of the design process.

#### 3.6 Wavelength dependence

The importance of adjusting an electromagnetic structure to the operating wavelength, such as changing the length of an antenna, is well known. In this study we demonstrate that significantly larger design adjustments than simple scaling are required to achieve the best Raman enhancement for a given operating wavelength. It is shown how the optimized nanostructure geometry, as well as the emission enhancement, changes significantly with wavelength over the interval $\lambda _{1} = \lambda _2 \in [250 \ \mathrm { nm},490 \ \mathrm { nm}]$. The study considers smaller designs with $R_{\mathrm {out}} = 30$ nm and assumes a negligible Raman shift. In a sense, it probes a similar quantity as the material-dependence study in Sec. 3.4, seeing as $\vert \chi \vert ^2 / \mathrm {Im}(\chi )$ for silver varies by orders of magnitude across the wavelength interval. However, this study considers the added complexity that the wavelength changes along with $\varepsilon (\textbf {r})$.

The top panel of Fig. 8 shows nine silver nanostructures (rainbow-colored) optimized and evaluated for the reported wavelengths (colored dots). The quantity $[\vert \chi \vert ^2 / \mathrm {Im}(\chi )]^2$ for silver is overlaid for easy reference (black). It is clearly observed that the enhancement factor of the optimized structures is strongly dependent on $[\vert \chi \vert ^2 / \mathrm {Im}(\chi )]^2$. Furthermore, it is seen that the optimized nanostructure geometries change significantly from left to right, starting with a reflector-like geometry for $\lambda = 250$ nm and ending with a geometry resembling the red structure in Fig. 4 optimized to maximize dipole emission. Interestingly, around $\lambda = 320$ nm (near the minimum of $[\vert \chi \vert ^2 / \mathrm {Im}(\chi )]^2$) the optimized structure is seen to experience a topological change which fundamentally changes the geometry from the reflector-like type to structures fully encapsulating the emitter. The bottom panel of Fig. 8 shows the per-wavelength max-normalized enhancement attained for each of the nine structures at each wavelength. From the panel it is seen that each of the optimized nanostructures outperformed the other structures at the wavelengths for which they were optimized. Furthermore, it is seen that as the wavelength increases the performance sensitivity increases. The data clearly shows the need for tailoring the nanostructure geometry to the particular operating conditions if the highest performance is sought.

#### 3.7 Out-of-plane Raman scattering

A common configuration for Raman-sensing applications is the distribution of Raman molecules over a surface, illuminated by some out-of-plane source with the scattered light collected out-of-plane as well. In the following study we demonstrate the strength of the proposed approach for such problem configurations (Fig. 2(b)).

We consider a periodically-patterned platinum surface (period $a = 150$ nm) submerged in a water background ($\varepsilon _{\mathrm {H}_2\mathrm {O}} \approx 1.77$ around 500 nm). For computational simplicity, the study assumes a single Raman molecule per unit cell, situated at a fixed position at the center of the 150 nm x 200 nm design domain, which is again placed immediately on top of the platinum surface. A circular region of radius $R_{\mathrm {in}} = 10$ nm, centered at the Raman molecule, is kept free of platinum throughout the optimization. The model problem is excited at $\lambda _{1} = 532$ nm with the Raman molecules emitting at $\lambda _2 = 549$ nm. The molecules in the periodic array are assumed to scatter incoherently, as this most accurately models the physical scattering process. To account for this incoherent scattering, the array scanning method is used in the design process to compute the Raman scattering from a single molecule in the periodic model using 50 $k$-points in $[-\pi /a,\pi /a]$ [57,58].

As a reference, we consider a periodic array ($a = 150$ nm) of circular platinum discs resting on top of the platinum surface, with the Raman molecules placed at the center of the gap between the discs (blue dot in Fig. 9(d)). The separation distance from the molecule to the discs is taken to be 10 nm, identical to the radius of the platinum-free circular region in the optimization case.

Figure 9 shows the optimized (left) and reference (right) surface structures and their behaviour under illumination at $\lambda _{1} = 532$ nm with a single Raman molecule placed in the center unit cell. The top row shows the periodic platinum surface structures (black) in water background (gray) with the blue dot indicating the position of the Raman molecule. The middle row shows $\vert \textbf {E}_1 \vert$ at $\lambda _{1} = 532$ nm for identical incident field strength using identical color scales. From the panels it is clear that the optimized structure provides a stronger focus of the incident field at the Raman molecule position. The bottom row presents $\vert \textbf {E}_2 \vert$ at $\lambda _{2} = 549$ nm using identical color scales, clearly illustrating the significantly enhanced emission for the optimized surface structure. Computing the power emitted by the Raman molecule for both cases reveals a Raman enhancement by a factor of $\sim 64\times$ for the optimized structure relative to the reference. While this example assumes a 2d model, it demonstrates the potential of applying the proposed approach to the design of nanostructures for full 3d out-of-plane Raman scattering problems.

#### 3.8 In-plane Raman scattering

In the final study we consider an in-plane Raman scattering problem with input and output waveguides (Fig. 2(c)). This study demonstrates the vast design freedom inherent to the proposed approach, which allows for the design of extreme enhancement structures in configurations where it would otherwise be difficult, if not impossible, to achieve good enhancement using conventional design techniques or intuition.

The objective of the design problem is to maximize the difference between the power emitted through $\Gamma _{\mathrm {out}}$ and through $\Gamma _{\mathrm {in}}$. This is achieved by designing a Pt nanostructure in $\boldsymbol \Omega _{\mathrm {D}_{\mathrm {Metal}}}$ together with a Si$_{3}$N$_{4}$ dielectric structure in $\boldsymbol \Omega _{\mathrm {D}_{\mathrm {Dielectric}}}$. Figure 10(a) shows the optimized nanostructure obtained by solving the design problem at $\lambda _{1}= 532$ nm, $\lambda _{2}= 549$ nm.

As a reference, a bowtie antenna is parameter-optimized over its tip angle $\theta _{\mathrm {bta}} \in [10^{\mathrm {o}},90^{\mathrm {o}}]$, side length $L_{\mathrm {bta}} \in [20 \ \mathrm { nm}, 100 \ \mathrm { nm}]$, and rotation angle relative to horizontal $\theta _{\mathrm {rot}} \in [-90^{\mathrm {o}},90^{\mathrm {o}}]$. The parameters leading to the largest objective value are found to be $\theta _{\mathrm {bta}} = 10^{\mathrm {o}}$, $L_{\mathrm {bta}} = 62.5$ nm, and $\theta _{\mathrm {rot}} = 25^{\mathrm {o}}$. For the dielectric material distribution in the reference a simple extension of the two waveguides towards the bowtie antenna is used. The reference nanostructure is shown in Fig. 10(b).

Figure 10(c)–10(d) shows the incident field at $\lambda _{1}= 532$ nm for the optimized and reference structure, respectively. Identical color scales are used for both plots, clearly showing the enhanced focusing of $\textbf {E}_{1}$ at the position of the Raman molecule. The rotated orientation of the reference structure, introduced to maximize the objective, means that the bowtie antenna is not capable of creating a highly localized focus between its tips. In contrast, the significantly more intricate geometry of the optimized design has no problem providing such enhancement.

Figure 10(e)–10(f) shows the field emitted by the Raman molecule at $\lambda _{2}= 549$ nm for the optimized structure and the reference structure, respectively. Again, identical color scales are used for the two plots. From these it is obvious that the optimized structure generates a significantly enhanced emission relative to the reference. In fact, the difference between the two structures in terms of emitted power through $\Gamma _{\mathrm {out}}$ is a factor of $\sim 345\times$, i.e. more than two orders of magnitude.

## 4. Conclusion

In this paper we proposed a TO-based approach for designing Raman scattering enhancing metallic nanostructures and presented a structure that increases the enhancement by a factor of $\sim 60\times$ compared to a parameter-optimized bowtie antenna (Sec. 3.1). Through a number of studies we documented: the importance of considering the full Raman scattering process in the design procedure (Sec. 3.2); that the expected linear scaling [21] of the Raman enhancement with design volume (area) is achieved (Sec. 3.3); that the Raman enhancement scales with $[\vert \chi \vert ^2/\mathrm {Im}(\chi )]^2$ rather then the predicted scaling of $[\vert \chi \vert ^2/\mathrm {Im}(\chi )]^3$ [21], possibly due to a trade off between the focusing and emission enhancement processes (Sec. 3.4); the importance of accounting for the Raman shift in the design procedure (Sec. 3.5); the importance of tailoring the nanostructure geometry to the operating wavelength as large geometric and associated performance changes occur as the wavelength is varied (Sec. 3.6). Finally, we demonstrated that the TO-based approach may be used to achieve $\sim 64\times$ enhancement for out-of-plane Raman scattering (Sec. 3.7) and $\sim 345\times$ enhancement for in-plane Raman scattering in a two-waveguide setup (Sec. 3.8), relative to simple parameter optimized reference structures.

While all studies in this work consider 2d model problems, the demonstration of extreme enhancements compared to well known reference geometries is promising. Hence, an important next step is the extension of the proposed approach to three dimensions. Applying the approach for Raman scattering problems modelling realistic operating conditions may help reveal hitherto undiscovered nanostructures for extreme Raman enhancement. If such structures are found, they may serve to improve existing Raman scattering based tools significantly as well as enable the development of novel tools.

An important topic for future work is to incorporate variability in the location of the Raman molecule. In some experimental situations, the Raman molecules are suspended in a fluid and the objective is to maximize the *average* emission over all possible molecule locations [32,64,65]. Naively, this would require solving the model problem for a vast number of different dipole locations and averaging the emission, making such a task computationally infeasible for all but the smallest of problems. However, efficient "trace" formulations have been developed for related problems involving a distribution of random emitters—thermal emission [66,67] and spontaneous emission [67]—and we believe that similar techniques are applicable to the Raman problem.

Another future step is the investigation of the sensitivity of the optimized structures towards perturbations of their geometry, such as those encountered in fabrication. Results in previous studies on designing plasmonic nanostructures using TO, e.g. [52], as well as several results in this work, suggest that the optimized nanostructures are sensitive to geometric variation, as small variations lead to large performance changes. By employing well known robust optimization techniques [43,68] it may be possible to reduce the geometric sensitivity of the nanostructures, hereby bridging a gap between simulations and experiment.

## Appendix A. Optimization, post processing and evaluation procedure

The physics is modelled in COMSOL Multiphysics [59] and the optimization problem is solved using the Globally Convergent Method of Moving Asymptotes (GCMMA) [38].

The following stopping criterion is used to terminate the optimization:

Here $i$ denotes the current optimization iteration, $i_{\mathrm {min}} = 70$ denotes the minimum number of iterations taken. $\beta$ denotes the thresholding strength used in the filtering procedure and $\beta _{\mathrm {max}} = 32$. $\Phi _{i}$ denotes the objective function value at the $i$’th iteration and $n \in \mathbb {N}^{+}$.

After the design process is completed a post-processing procedure is performed to extract the final nanostructure from the optimized material distribution. In this procedure the final $\epsilon (\textbf {r})$-field is sampled in a Cartesian (x,y)-grid with 0.1 nm resolution. The sampled distribution is smoothed using a simple cone-shaped kernel [41] with 1.5 nm filter radius to remove all sub-nanometer kinks left by the 1 nm topology optimization mesh. The final geometry is then extracted as the $0.5$-contour of the smoothed field. Finally, the geometry is rescaled to have an inner radius of $r_{\mathrm {inner}} = 10$ nm for consistency across designs.

The reported performances and fields are evaluated using the final post processed design, discretized using an unstructured triangular body fitted mesh with 1 nm side length for the structure. Quadratic continuous Lagrange basis functions are used to resolve the physics during the evaluation step.

Note that the exact position and magnitude of the Raman enhancement peak depend strongly on the final geometry. This means that a size scaling of a structure by a few percent may result in a shift of several nm in the emission peak position as well as a change in the overall enhancement.

## Appendix B. Study parameters

This appendix lists the parameters used in setting up the model and associated optimization problems. Unless specified otherwise below or in the individual sections detailing each study, the parameters listed in the following are used.

## B. 1. Model domain

For the model domain sketched in Fig. 2(a) the following values are used: The model domain $\boldsymbol {\Omega }$ is a square with side length, $L_{\boldsymbol {\Omega }} = 600$ nm. $\boldsymbol {\Omega }$ is surrounded by a perfectly matched layer with a depth of $D_{\mathrm {PML}}= 300$ nm. The design domain $\boldsymbol {\Omega }_{\mathrm {D}}$ is centered in $\boldsymbol {\Omega }$ and taken to be a perfect 2d torus with inner radius, $R_{\mathrm {in}} = 10$ nm and outer radius, $R_{\mathrm {out}} = 100$ nm. The Raman molecule is modelled as a point dipole source placed at the center of $\boldsymbol {\Omega }_{\mathrm {D}}$. The power emitted by the dipole at $\lambda _{2}$ is computed by evaluating Eq. (4) over $\boldsymbol {\Gamma }_{\mathrm {RE}},$ a circular curve centered on the dipole with radius: $R_{\boldsymbol {\Gamma }_{\mathrm {RE}]}} = 250$ nm.

For the model domain sketched in Fig. 2(b) the following values are used: The model domain $\boldsymbol {\Omega }$ is a rectangle of width, $\mathrm {a} = W_{\boldsymbol {\Omega }} = 150$ nm and height $H_{\boldsymbol {\Omega }} = 600$ nm. $\boldsymbol {\Omega }$ is truncated from above by perfectly matched layer with a depth of $D_{\mathrm {PML}}= 300$ nm. $\boldsymbol {\Omega }$ is truncated from below using a perfect electric conductor boundary condition: $\textbf {n} \times \textbf {E} = 0$. $\boldsymbol {\Omega }$ is truncated by a Floquet Bloch boundary condition on the left and right boundaries to model the $x$-periodicity. The top 500 nm of $\boldsymbol {\Omega }$ is taken to contain water ($\varepsilon _r \approx 1.77$). The bottom 100 nm of $\boldsymbol {\Omega }$ is taken to contain platinum [63]. The design domain $\boldsymbol {\Omega }_{\mathrm {D}}$ is placed immediately on top of the platinum region. The design domain $\boldsymbol {\Omega }_{\mathrm {D}}$ has a width of $W_{\boldsymbol {\Omega }_{\mathrm {D}}} = 150$ nm and a height $H_{\boldsymbol {\Omega }_{\mathrm {D}}} = 200$ nm. The design domain $\boldsymbol {\Omega }_{\mathrm {D}}$ has a region fixed to contain water at its center with radius $R_{\mathrm {in}} = 10$ nm. The Raman molecule is modelled as a point dipole source placed at the center of $\boldsymbol {\Omega }_{\mathrm {D}}$. The power emitted by the dipole at $\lambda _{2}$ is computed by evaluating Eq. (4) over $\boldsymbol {\Gamma }_{\mathrm {RE}},$ a straight horizontal line 50 nm below the top boundary of $\boldsymbol {\Omega }$.

For the model domain sketched in Fig. 2(c) the following values are used: The model domain $\boldsymbol {\Omega }$ is a square with side length, $L_{\boldsymbol {\Omega }} = 1200$ nm centered at (0 nm,0 nm). $\boldsymbol {\Omega }$ is surrounded by a perfectly matched layer with a depth of $D_{\mathrm {PML}}= 300$ nm. The design domain $\boldsymbol {\Omega }_{\mathrm {D}_{\mathrm {Metal}}}$ is centered at (-300 nm,-300 nm). $\boldsymbol {\Omega }_{\mathrm {D}_{\mathrm {Metal}}}$ is a square of side length $L_{\boldsymbol {\Omega }_{\mathrm {D}_{\mathrm {Metal}}}} = 200$ nm with a circular hole of radius $R_{\mathrm {in}} = 10$ nm at its center. The Raman molecule is modelled as a point dipole source placed at the center of $\boldsymbol {\Omega }_{\mathrm {D}}$. The design domain $\boldsymbol {\Omega }_{\mathrm {D}_{\mathrm {Dielectric}}}$ is centered at (-200 nm,-200 nm). $\boldsymbol {\Omega }_{\mathrm {D}_{\mathrm {Dielectric}}}$ is the difference between a square of side length $L_{\boldsymbol {\Omega }_{\mathrm {D}_{\mathrm {Dielectric}}}} = 400$ nm and $\boldsymbol {\Omega }_{\mathrm {D}_{\mathrm {Metal}}}$. The input waveguide runs horizontally and is centered at $y = -300$ nm. It starts at the left edge of $\boldsymbol {\Omega }$ and runs until $\boldsymbol {\Omega }_{\mathrm {D}_{\mathrm {Dielectric}}}$. The output waveguide runs vertically and is centered at $x = 300$ nm. It starts at the top edge of $\boldsymbol {\Omega }$ and runs until $\boldsymbol {\Omega }_{\mathrm {D}_{\mathrm {Dielectric}}}$. The width of both waveguides is 150 nm. The power emitted by the dipole at $\lambda _{2}$ into the input and output waveguides is computed by evaluating Eq. (4) over $\boldsymbol {\Gamma }_{\mathrm {in}}$ and $\boldsymbol {\Gamma }_{\mathrm {out}}$ respectively. $\boldsymbol {\Gamma }_{\mathrm {in}}$ is a vertical line through the input waveguide at $x = -300$ nm. $\boldsymbol {\Gamma }_{\mathrm {out}}$ is a horizontal line through the output waveguide at $t = 300$ nm.

## B. 2. Numerical solution

For the numerical solution of all model problems the geometry is discretized using an unstructured first order finite element mesh [40]. The design domain is discretized using triangular elements with a uniform side length of $h = 1$ nm, dictating the resolution of the design. The remaining model domain is discretized using a non-uniform mesh of triangular elements with a smallest side lengths of $h=1$ nm near the design domain and a maximum side length of $h = 1/16$ effective wavelength ($\lambda / n$), to ensure the accuracy of the model.

The design problems considered in this work were solved on a laptop with 16 GB RAM and an Intel(R) Core(TM) i5-7200 CPU @ 2.50 GHz processor. The approximate wall time spent per design iteration for the design problem in Sec. 3.1 was 90 seconds, resulting in approximately 15 hours spent executing the full design process. It is noted that these numbers are naturally strongly dependent on the number of design iterations, number of elements in the mesh, and by extension on the size of the model and design domain, and finally that computational speed was not the focus of this work.

For the model problems sketched in Fig. 2(a)–2(b) the incident-field problem $\big ($Eq. (1)$\big )$, is solved using the scattered-field approach [40], where the background field is taken to be a perfect planewave.

For the model problems sketched in Fig. 2(c) the incident-field problem is solved using the total-field approach [40]. First the lowest-order mode confined to the input waveguide is computed. This mode is then used to excite the system and a total-field problem is solved.

For all model problems the emitter problem $\big ($Eq. (3)$\big )$ is solved for the total field.

## B. 3. Physics related

The physics related parameters are chosen as follows: The wavelength of the incident field is taken to be $\lambda _{1} = 532$ nm. The wavelength of the emitted field is taken to be $\lambda _{2} = 532$ nm. The relative permittivity and relative permeability of air are taken to be $\varepsilon _{\mathrm {air}} = \mu _{\mathrm {air}} = 1.0$. The relative permittivity and relative permeability of water are taken to be $\varepsilon _{\mathrm {water}} = 1.77$, $\mu _{\mathrm {water}} = 1.0$ at the operating wavelengths. The relative permittivity and relative permeability of silver (Ag), gold (Au) and copper (Cu) are taken from [62]. The relative permittivity and relative permeability of platinum (Pt) are taken from [63]. The relative permittivity and relative permeability of Si$_3$N$_4$ are taken to be $\varepsilon _{\mathrm {Si}_3\mathrm {N}_4} = 2.05$, $\mu _{\mathrm {water}} = 1.0$ at the operating wavelengths. The speed of light is taken to be $c = 3 \cdot 10^8$ m/s.

## B. 4. Design related

For the designs considering the model problem sketched in Fig. 2(a) mirror symmetry is imposed on the designs normal to the horizontal axis due to the nature of the model problem, where the planewave in Eq. (1) propagates through the domain along the horizontal axis.

For the designs considering the model problem sketched in Fig. 2(b) mirror symmetry is imposed on the design normal to the vertical axis due to the nature of the model problem, where the planewave in Eq. (1) propagates through the domain along the vertical axis.

## B. 5. Optimization related

For the model problem sketched in Fig. 2(a) the following values are used: As an initial guess for all optimizations a full metal disc is used, i.e. $\xi _{\mathrm {initial}}(\textbf {r}) = 1 \ \forall \ \textbf {r} \in \boldsymbol {\Omega }_{\mathrm {D}}$. The filter radius in the smoothing operation is, $r_f = 10$ nm. ($r_f = 5$ nm was used for the study detailed in Sec. 3.6.) The thresholding level in the thresholding operation is, $\eta = 0.5$. The initial thresholding strength is, $\beta _{\mathrm {initial}} = 8$. The thresholding strength is increased gradually during the optimization through the values, $\beta = 8,16,32$. A minimum length scale of all features in the designs is ensured using geometric length-scale constraints with $c_{\mathrm {LS}} = 625$, $\eta _{e} = 0.75$, $\eta _{d} = 0.25$. The length-scale constraints are only enforced for $\beta = 32$ to allow the design to form freely without limiting topology changes in the earlier stages of the optimization.

For the model problem sketched in Fig. 2(b) the following values are used: As an initial guess a full metal region is used, i.e. $\xi _{\mathrm {initial}}(\textbf {r}) = 1 \ \forall \ \textbf {r} \in \boldsymbol {\Omega }_{\mathrm {D}}$. The filtering radius used in the smoothing operation is, $r_f = 10$ nm. The thresholding level used is in the thresholding operation is, $\eta = 0.5$. The initial thresholding strength used in the thresholding operation is, $\beta _{\mathrm {initial}} = 8$. The thresholding strength is increased gradually during the optimization through the values, $\beta = 8,16,32$. A minimum length scale of all features in the designs is ensured using geometric length-scale constraints using $c_{\mathrm {LS}} = 625$, $\eta _{e} = 0.75$, $\eta _{d} = 0.25$. The length-scale constraints are only enforced for $\beta = 32$ to allow the design to form freely without limiting topology changes in the earlier stages of the optimization.

For the model problem sketched in Fig. 2(c) the following values are used: As an initial guess for $\boldsymbol {\Omega }_{\mathrm {D}_{\mathrm {Metal}}}$ a full metal region is used, i.e. $\xi _{\mathrm {initial}}(\textbf {r}) = 1 \ \forall \ \textbf {r} \in \boldsymbol {\Omega }_{\mathrm {D}_{\mathrm {Metal}}}$. As an initial guess for $\boldsymbol {\Omega }_{\mathrm {D}_{\mathrm {Dielectric}}}$ a full dielectric region is used, i.e. $\xi _{\mathrm {initial}}(\textbf {r}) = 1 \ \forall \ \textbf {r} \in \boldsymbol {\Omega }_{\mathrm {D}_{\mathrm {Dielectric}}}$. The filtering radius in the smoothing operation is: $r_f = 10$ nm. The thresholding level in the thresholding operation is: $\eta = 0.5$. The initial thresholding strength used is in the thresholding operation is: $\beta _{\mathrm {initial}} = 8$. The thresholding strength is increased gradually during the optimization through the values: $\beta = 8,16,32$. A minimum length scale of all features in the designs is ensured using geometric length-scale constraints using $c_{\mathrm {LS}} = 625$, $\eta _{e} = 0.75$, $\eta _{d} = 0.25$. The length-scale constraints is only imposed for $\beta = 32$ to allow the design to form freely without limiting topology changes in the earlier stages of the optimization.

## Funding

Villum Fonden (8692); U.S. Army Research Office (W911NF-13-D-0001); National Science Foundation (1453218, 1709212).

## Acknowledgments

The authors thank Juejun Hu for insightful discussions.

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **M. P. Bendsøe and O. Sigmund, * Topology Optimization* (Springer, 2003).

**2. **J. S. Jensen and O. Sigmund, “Topology optimization for nano-photonics,” Laser Photonics Rev. **5**(2), 308–321 (2011). [CrossRef]

**3. **S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vuckovic, and A. W. Rodriguez, “Inverse design in nanophotonics,” Nat. Photonics **12**(11), 659–670 (2018). [CrossRef]

**4. **D. A. Long, * Raman spectroscopy* (McGraw-Hill, 1977).

**5. **G. Turrell and J. Corset, * Raman Microscopy: Developments and Applications* (Academic, 1996).

**6. **N. B. Colthup, L. H. Daly, and S. E. Wiberley, * Introduction to Infrared and Raman Spectroscopy* (Academic, 1990).

**7. **W. Huang, W. Qian, P. K. Jain, and M. A. El-Sayed, “The effect of plasmon field on the coherent lattice phonon oscillation in electron-beam fabricated gold nanoparticle pairs,” Nano Lett. **7**(10), 3227–3234 (2007). [CrossRef]

**8. **W. Zhu, M. G. Banaee, D. Wang, Y. Chu, and K. B. Crozier, “Lithographically fabricated optical antennas with gaps well below 10 nm,” Small **7**(13), 1761–1766 (2011). [CrossRef]

**9. **W. Rechberger, A. Hohenau, A. Leitner, J. R. Krenn, B. Lamprecht, and F. R. Aussenegg, “Optical properties of two interacting gold nanoparticles,” Opt. Commun. **220**(1-3), 137–141 (2003). [CrossRef]

**10. **E. Hao and G. C. Schatz, “Electromagnetic fields around silver nanoparticles and dimers,” J. Chem. Phys. **120**(1), 357–366 (2004). [CrossRef]

**11. **S. Dodson, M. Haggui, R. Bachelot, J. Plain, S. Li, and Q. Xiong, “Optimizing electromagnetic hotspots in plasmonic bowtie nanoantennae,” J. Phys. Chem. Lett. **4**(3), 496–501 (2013). [CrossRef]

**12. **M. Kaniber, K. Schraml, A. Regler, J. Bartl, G. Glashagen, F. Flassig, J. Wierzbowski, and J. J. Finley, “Surface plasmon resonance spectroscopy of single bowtie nano-antennas using a differential reflectivity method,” Sci. Rep. **6**(1), 23203 (2016). [CrossRef]

**13. **W. Yue, Z. Wang, J. Whittaker, F. Lopez-Royo, Y. Yang, and A. V. Zayats, “Amplification of surface-enhanced Raman scattering due to substrate-mediated localized surface plasmons in gold nanodimers,” J. Mater. Chem. C **5**(16), 4075–4084 (2017). [CrossRef]

**14. **J. Zhang, M. Irannejad, and B. Cui, “Bowtie Nanoantenna with Single-Digit Nanometer Gap for Surface-Enhanced Raman Scattering (SERS),” Plasmonics **10**(4), 831–837 (2015). [CrossRef]

**15. **M. Moskovits, “Surface-enhanced spectroscopy,” Rev. Mod. Phys. **57**(3), 783–826 (1985). [CrossRef]

**16. **A. Campion and P. Kambhampati, “Surface-enhanced Raman scattering,” Chem. Soc. Rev. **27**(4), 241–250 (1998). [CrossRef]

**17. **K. Kneipp, M. Moskovits, and H. Kneipp, * Surface-Enhanced Raman Scattering* (Springer, 2006).

**18. **E. Le Ru and P. Etchegoin, * Principles of Surface-Enhanced Raman Spectroscopy: and Related Plasmonic Effects* (Elsevier, 2008).

**19. **P. L. Stiles, J. A. Dieringer, N. C. Shah, and R. P. Van Duyne, “Surface-Enhanced Raman Spectroscopy,” Annu. Rev. Anal. Chem. **1**(1), 601–626 (2008). [CrossRef]

**20. **C. L. Haynes, A. D. McFarland, and R. P. Van Duyne, “Surface-Enhanced Raman Spectroscopy,” Anal. Chem. **77**(17), 338A–346 A (2005). [CrossRef]

**21. **J. Michon, M. Benzaouia, W. Yao, O. D. Miller, and S. G. Johnson, “Limits to surface-enhanced raman scattering near arbitrary-shape scatterers,” Opt. Express **27**(24), 35189–35202 (2019). [CrossRef]

**22. **O. D. Miller, A. G. Polimeridis, M. T. H. Reid, C. W. Hsu, B. G. DeLacy, J. D. Joannopoulos, M. Soljačić, and S. G. Johnson, “Fundamental limits to optical response in absorptive systems,” Opt. Express **24**(4), 3329–3364 (2016). [CrossRef]

**23. **R. D. Averitt, S. L. Westcott, and N. J. Halas, “Linear optical properties of gold nanoshells,” J. Opt. Soc. Am. B **16**(10), 1824 (1999). [CrossRef]

**24. **C. F. Bohren and D. R. Huffman, * Absorption and Scattering of Light by Small Particles* (Wiley, 1998).

**25. **A. Ahmed and R. Gordon, “Single molecule directivity raman scattering using nanoantennas,” Nano Lett. **12**(5), 2625–2630 (2012). [CrossRef]

**26. **S. Nie and S. R. Emory, “Probing Single Molecules and Single Nanoparticles by Surface-Enhanced Raman Scattering,” Science **275**(5303), 1102–1106 (1997). [CrossRef]

**27. **K. Kneipp, Y. Wang, H. Kneipp, L. T. Perelman, I. Itzkan, R. R. Dasari, and M. S. Feld, “Single molecule detection using surface-enhanced raman scattering (SERS),” Phys. Rev. Lett. **78**(9), 1667–1670 (1997). [CrossRef]

**28. **B. Sharma, R. R. Frontiera, A. I. Henry, E. Ringe, and R. P. Van Duyne, “SERS: Materials, applications, and the future,” Mater. Today **15**(1-2), 16–25 (2012). [CrossRef]

**29. **S. Fateixa, H. I. Nogueira, and T. Trindade, “Hybrid nanostructures for SERS: Materials development and chemical detection,” Phys. Chem. Chem. Phys. **17**(33), 21046–21071 (2015). [CrossRef]

**30. **P. A. Mosier-Boss, “Review of SERS substrates for chemical sensing,” Nanomaterials **7**(6), 142 (2017). [CrossRef]

**31. **J. P. Camden, J. A. Dieringer, Y. Wang, D. J. Masiello, L. D. Marks, G. C. Schatz, and R. P. Van Duyne, “Probing the structure of single-molecule surface-enhanced Raman scattering hot spots,” J. Am. Chem. Soc. **130**(38), 12616–12617 (2008). [CrossRef]

**32. **E. C. Le Ru, E. Blackie, M. Meyer, and P. G. Etchegoint, “Surface enhanced raman scattering enhancement factors: A comprehensive study,” J. Phys. Chem. C **111**(37), 13794–13803 (2007). [CrossRef]

**33. **D. A. Genov, A. K. Sarychev, V. M. Shalaev, and A. Wei, “Resonant Field Enhancements from Metal Nanoparticle Arrays,” Nano Lett. **4**(1), 153–158 (2004). [CrossRef]

**34. **A. Sundaramurthy, K. B. Crozier, G. S. Kino, D. P. Fromm, P. J. Schuck, and W. E. Moerner, “Field enhancement and gap-dependent resonance in a system of two opposing tip-to-tip Au nanotriangles,” Phys. Rev. B: Condens. Matter Mater. Phys. **72**(16), 165409 (2005). [CrossRef]

**35. **S. Molesky, W. Jin, P. S. Venkataram, and A. W. Rodriguez, “Bounds on absorption and thermal radiation for arbitrary objects,” arXiv:1907.04418.

**36. **D. A. Tortorelli and P. Michaleris, “Design sensitivity analysis: Overview and review,” Inverse Probl. Eng. **1**(1), 71–105 (1994). [CrossRef]

**37. **J. Nocedal and S. J. Wright, * Numerical Optimization - Second Edition* (Springer Science+Business Media LLC, 2006).

**38. **K. Svanberg, “A class of globally convergent optimization methods based on conservative convex separable approximations,” SIAM J. Control **12**(2), 555–573 (2002). [CrossRef]

**39. **N. Aage, E. Andreassen, B. S. Lazarov, and O. Sigmund, “Giga-voxel computational morphogenesis for structural design,” Nature **550**(7674), 84–86 (2017). [CrossRef]

**40. **J.-M. Jin, * The Finite Element Method in Electromagnetics - Third Edition* (Wiley-IEEE, 2014).

**41. **B. Bourdin, “Filters in topology optimization,” Int. J. Numer. Meth. Eng. **50**(9), 2143–2158 (2001). [CrossRef]

**42. **J. K. Guest, J. H. Prévost, and T. Belytschko, “Achieving minimum length scale in topology optimization using nodal design variables and projection functions,” Int. J. Numer. Meth. Eng. **61**(2), 238–254 (2004). [CrossRef]

**43. **F. Wang, B. S. Lazarov, and O. Sigmund, “On projection methods, convergence and robust formulations in topology optimization,” Struct. Multidicip. O. **43**(6), 767–784 (2011). [CrossRef]

**44. **R. E. Christiansen, J. Vester-Petersen, S. P. Madsen, and O. Sigmund, “A non-linear material interpolation for design of metallic nano-particles using topology optimization,” Comput. Method. Appl. M. **343**, 23–39 (2019). [CrossRef]

**45. **M. Zhou, B. S. Lazarov, F. Wang, and O. Sigmund, “Minimum length scale in topology optimization by geometric constraints,” Comput. Method. Appl. M. **293**, 266–282 (2015). [CrossRef]

**46. **W. R. Frei, H. T. Johnson, and K. D. Choquette, “Optimization of a single defect photonic crystal laser cavity,” J. Appl. Phys. **103**(3), 033102 (2008). [CrossRef]

**47. **J. Lu, S. Boyd, and J. Vuckovic, “Inverse design of a three-dimensional nanophotonic resonator,” Opt. Express **19**(11), 10563–10570 (2011). [CrossRef]

**48. **X. Liang and S. G. Johnson, “Formulation for scalable optimization of microcavities via the frequency-averaged local density of states,” Opt. Express **21**(25), 30812–30841 (2013). [CrossRef]

**49. **F. Wang, R. E. Christiansen, Y. Yu, J. Mørk, and O. Sigmund, “Maximizing the quality factor to mode volume ratio for ultra-small photonic crystal cavities,” Appl. Phys. Lett. **113**(24), 241101 (2018). [CrossRef]

**50. **E. Wadbro and C. Engström, “Topology and shape optimization of plasmonic nano-antennas,” Comput. Method. Appl. M. **293**, 155–169 (2015). [CrossRef]

**51. **Y. Deng, Z. Liu, C. Song, P. Hao, Y. Wu, Y. Liu, and J. Korvink, “Topology optimization of metal nanostructures for localized surface plasmon resonances,” Struct. Multidicip. O. **53**(5), 967–972 (2016). [CrossRef]

**52. **J. Vester-Petersen, R. Christiansen, B. Julsgaard, P. Balling, O. Sigmund, and S. Madsen, “Topology optimized gold nanostrips for enhanced near-infrared photon upconversion,” Appl. Phys. Lett. **111**(13), 133102 (2017). [CrossRef]

**53. **J. Vester-Petersen, S. P. Madsen, O. Sigmund, P. Balling, B. Julsgaard, and R. E. Christiansen, “Field-enhancing photonic devices utilizing waveguide coupling and plasmonics - a selection rule for optimization-based design,” Opt. Express **26**(18), A788–2001 (2018). [CrossRef]

**54. **H. Chung and O. D. Miller, “High-na achromatic metalenses by inverse design,” arXiv:1905-09213v1.

**55. **Z. Lin, X. Liang, M. Loncar, S. G. Johnson, and A. W. Rodriguez, “Cavity-enhanced second-harmonic generation via nonlinear-overlap optimization,” Optica **3**(3), 233–238 (2016). [CrossRef]

**56. **D. J. Griffiths, * Introduction to Electrodymanics - Fourth Edition* (Pearson Education Limited, 2014).

**57. **R. A. Sigelmann and A. Ishimaru, “Radiation from periodic structures excited by an aperiodic source,” IRE Trans. Antennas Propag. **13**(3), 354–364 (1965). [CrossRef]

**58. **F. Capolino, D. R. Jackson, D. R. Wilton, and L. B. Felsen, “Comparison of methods for calculating the field excited by a dipole near a 2-d periodic material,” IRE Trans. Antennas Propag. **55**(6), 1644–1655 (2007). [CrossRef]

**59. **“Comsol multiphysics^{®} v. 5.4. www.comsol.com,”.

**60. **O. Sigmund, “On the usefulness of non-gradient approaches in topology optimization,” Struct. Multidicip. O. **43**(5), 589–596 (2011). [CrossRef]

**61. **G. Angeris, J. Vuckovic, and S. Boyd, “Computational bounds for photonic design,” arXiv:1811.12936v3 43, 589–596 (2011).

**62. **P. B. Johnson and R. W. Christy, “Optical constants of noble metals,” Phys. Rev. B **6**(12), 4370–4379 (1972). [CrossRef]

**63. **W. S. M. Werner, “Optical constants and inelastic electron-scattering data for 17 elemental metals,” J. Phys. Chem. Ref. Data **38**(4), 1013–1092 (2009). [CrossRef]

**64. **S.-Y. Ding, E.-M. You, Z.-Q. Tian, and M. Moskovits, “Electromagnetic theories of surface-enhanced raman spectroscopy,” Chem. Soc. Rev. **46**(13), 4042–4076 (2017). [CrossRef]

**65. **H. K. Lee, Y. H. Lee, C. S. L. Koh, G. C. Phan-Quang, X. Han, C. L. Lay, H. Y. F. Sim, Y.-C. Kao, Q. An, and X. Y. Ling, “Designing surface-enhanced raman scattering (sers) platforms beyond hotspot engineering: emerging opportunities in analyte manipulations and hybrid materials,” Chem. Soc. Rev. **48**(3), 731–756 (2019). [CrossRef]

**66. **A. W. Rodriguez, M. T. H. Reid, and S. G. Johnson, “Fluctuating surface current formulation of radiative heat transfer for arbitrary geometries,” Phys. Rev. B **86**(22), 220302 (2012). [CrossRef]

**67. **A. G. Polimeridis, M. T. H. Reid, W. Jin, S. G. Johnson, J. K. White, and A. W. Rodriguez, “Fluctuating volume-current formulation of electromagnetic fluctuations in inhomogeneous media: Incandescence and luminescence in arbitrary geometries,” Phys. Rev. B **92**(13), 134202 (2015). [CrossRef]

**68. **R. E. Christiansen, F. Wang, B. S. Lazarov, and O. Sigmund, “Creating geometrically robust designs for highly sensitive problems using topology optimization - acoustic cavity design,” Struct. Multidicip. O. **52**(4), 737–754 (2015). [CrossRef]