Inverse design of nanoparticles for enhanced Raman scattering

We show that topology optimization (TO) of metallic resonators can lead to ∼102 × 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


Introduction
In this paper, we present freeform shape optimization ("topology optimization," TO [1][2][3]) of metallic nanoparticles for Raman scattering [4][5][6], and obtain non-intuitive structures ∼ 60× better (in terms of emitted power) than optimized coupled-sphere [7][8][9] or bowtie [10][11][12][13][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 ω 1 interacting with a molecule that subsequently emits electromagnetic radiation at ω 2 . A nanostructure can enhance both the incident-wave absorption by a focusing effect as well as the emission by a Purcell effect [15][16][17][18][19][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 ∼ 5× (Sec. 3.2). When the Raman shift ω 1 − ω 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 ∼ 5× (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 χ, the key figure of merit for light-matter interactions is F = | χ| 2 /Im( χ) [22][23][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. Raman scattering from a molecule (point dipole) placed at the centre of a (gray) topologyoptimized silver structure at λ 1 = λ 2 = 532 nm. The color scheme is saturated due to the diverging field at the source. 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][16][17][18][19][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][29][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][32][33][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][42][43][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][51][52][53][54]. Nonlinearly coupled electromagnetic resonances at multiple frequencies were optimized via TO for second harmonic generation [55].

Model formulation
The Raman scattering phenomenon (sketched in Fig. 1(a)) is modelled as two sequentiallycoupled, time-harmonic, classical electromagnetic problems [40,56] in a domain Ω. 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 Ω 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 Here Ω ∈ R 2 represents the modelling domain and r and r 0 denote points in Ω. The material parameters µ and ε(r) are the magnetic permeability and the electric permittivity, respectively (non-magnetic materials are assumed in the following, i.e. µ r ≡ 1). Further, ω i = c/λ i denotes the angular frequency with c being the speed of light and λ i the wavelength, E i the electric field and f i the excitation sources for i ∈ 1, 2. Finally α denotes the polarizability tensor, which in this work is taken to be the identity αE 1 (r 0 ) = IE 1 (r 0 ) . The problem of enhancing the power emitted from the Raman molecule at λ 2 , by tailoring the material distribution in the design domain Ω 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 Here S 2 = Re 1 2 E 2 × H * 2 denotes the time averaged Poynting vector computed from E 2 , which in turn is obtained by solving Eqs. (1)-(3) for a given ξ(r) and (•) * denotes the complex conjugate. The vector n is the outward pointing normal vector for the integration curve Γ out ( Fig. 1(a)). The permittivity distribution ε(r) is determined by the value of ξ(r) through the interpolation [44], Here M i denotes the materials considered (metal or air) and n and κ denote the refractive index and extinction coefficient of M i , respectively. Finally, i denotes the imaginary unit. 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 Ω D (gray region) surrounding the molecule. An incident planewave propagates through the domain from left to right (green) and the power emitted through Γ 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 Ω D Metal (dark gray), in which a metallic nanostructure is designed. In the remaining part of the design domain, Ω D 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 Γ out and Γ 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 Ω. 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.

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.

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 λ 1 = 532 nm, assuming a negligible Raman shift (λ 2 = λ 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][8][9], while placing it between two identical metallic triangular structures (a bowtie antenna) further enhances the process [10][11][12][13][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 csa ∈ [10 nm, 50 nm], and we find the largest enhancement for R csa = 48 nm. The bowtie antenna is optimized by sweeping the tip-angle over, θ bta ∈ [10 o , 180 o ], and the side length over L bta ∈ [20 nm, 100 nm], and we find θ bta = 15 o and L bta = 70 nm to yield the largest enhancement at λ 1 = λ 2 = 532 nm. Figure 3 shows the enhancement of the emitted power Eq. (4) 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 ≈ 3 · 10 2 , while the bowtie antenna (red) achieves a peak enhancement of ≈ 1.3 · 10 4 at 532 nm. Interestingly, the intricate geometry of the topology-optimized nanostructure (blue), fully encapsulating the Raman molecule, achieves an enhancement of ≈ 8.0 · 10 5 . This is an increase by a factor of ∼ 60× 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 ξ(r) = 0.0 ∀r ∈ Ω 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 ∼ 2× compared to the fully encapsulating design in Fig. 3. In Fig. 3. Raman enhancement as a function of wavelength, for a molecule placed at the center of different silver nanostructures (dark blue dot) relative to a molecule placed in free space. A topology-optimized structure (blue), a bowtie antenna (red) and coupled-cylinder antenna (black) are considered, all designed to maximize enhancement at λ = 532 nm. 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].

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 Eq. (1) or a dipole-emission Eq. (3) 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 λ 1 = 532 nm and λ 2 = 582 nm.
For the dipole-emission case, the optimization problem Eqs. (4)-(5) 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. f 2 = −ω 2 2 I 1, 0 δ(r − r 0 ), effectively removing the first model problem Eq. (1) as f 2 no longer depends on 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 |E 1 | 2 for an incident planewave at the position of the dipole [44,50] effectively removing the second model problem Eq. (3) , 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 |E 1 | 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× 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 (≈ 10 3 ) comes within a factor of four of the bound (≈ 3.5 · 10 3 ), where difference may be explained by the length-scale imposed on the design. The focusing enhancement achieved by the black design (≈ 10 3 ) is, on the other hand, nearly three orders of magnitude from the focusing bound (≈ 6 · 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.

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 Ω D . These are R out = 50 nm (black), R out = 75 nm (red) and R 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. λ 1 = λ 2 = 532 nm. Figure 5 shows the Raman enhancement relative to a molecule in free space versus the area of Ω 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.

Material scaling
The upper bound for enhancing the Raman scattering process, derived in Ref. [21], scales cubically with the material-related quantity F = | χ| 2 /Im( χ). 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 ε(r), were kept constant across the four cases and a negligible Raman shift was assumed with λ 1 = λ 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 | χ| 2 /Im( χ) 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.

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 λ 1 and λ 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 λ 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 ∼ 1.6× is observed for the 50 nm Raman shift while a factor of ∼ 4.5× 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.

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 λ 1 = λ 2 ∈ [250 nm, 490 nm]. The study considers smaller designs with R out = 30 nm and assumes a negligible Raman shift. In a sense, it probes a similar quantity as the materialdependence study in Sec. 3.4, seeing as | χ| 2 /Im( χ) for silver varies by orders of magnitude across the wavelength interval. However, this study considers the added complexity that the wavelength changes along with ε(r).
The top panel of Fig. 8 shows nine silver nanostructures (rainbow-colored) optimized and evaluated for the reported wavelengths (colored dots). The quantity [| χ| 2 /Im( χ)] 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 [| χ| 2 /Im( χ)] 2 . Furthermore, it is seen that the optimized nanostructure geometries change significantly from left to right, starting with a reflector-like geometry for λ = 250 nm and ending with a geometry resembling the red structure in Fig. 4 optimized to maximize dipole emission. Interestingly, around λ = 320 nm (near the minimum of [| χ| 2 /Im( χ)] 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.

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-ofplane 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 (ε H 2 O ≈ 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 in = 10 nm, centered at the Raman molecule, is kept free of platinum throughout the optimization. The model problem is excited at λ 1 = 532 nm with the Raman molecules emitting at λ 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 [−π/a, π/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 λ 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 |E 1 | at λ 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 |E 2 | at λ 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 ∼ 64× 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.

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 Γ out and through Γ in . This is achieved by designing a Pt nanostructure in Ω D Metal together with a Si 3 N 4 dielectric structure in Ω D Dielectric . Figure 10(a) shows the optimized nanostructure obtained by solving the design problem at λ 1 = 532 nm, λ 2 = 549 nm. The parameters leading to the largest objective value are found to be θ bta = 10 o , L bta = 62.5 nm, and θ rot = 25 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 λ 1 = 532 nm for the optimized and reference structure, respectively. Identical color scales are used for both plots, clearly showing the enhanced focusing of 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 λ 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 Γ out is a factor of ∼ 345×, i.e. more than two orders of magnitude.  (Fig. 2(c)). (a) Optimized platinum (black) and Si 3 N 4 structure (white) with Si 3 N 4 input/output waveguides (white) and water background (gray). (b)) Reference platinum bowtie antenna (black) input/output waveguides (white) and water background (gray). (c)-(d)) |E 1 |-field at λ 1 = 532 nm resulting from the interaction between the incident field through the lower left waveguide and the (c) optimized and (d) reference structure. (e)-(f) |E 2 |-field emitted at λ 2 = 549 nm from a Raman molecule positioned at the center of (e) the optimized and (f) reference platinum nanostructure.

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 ∼ 60× 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 [| χ| 2 /Im( χ)] 2 rather then the predicted scaling of [| χ| 2 /Im( χ)] 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 ∼ 64× enhancement for out-of-plane Raman scattering (Sec. 3.7) and ∼ 345× 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 min = 70 denotes the minimum number of iterations taken. β denotes the thresholding strength used in the filtering procedure and β max = 32. Φ i denotes the objective function value at the i'th iteration and n ∈ 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 (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 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 Ω is a square with side length, L Ω = 600 nm. Ω is surrounded by a perfectly matched layer with a depth of D PML = 300 nm. The design domain Ω D is centered in Ω and taken to be a perfect 2d torus with inner radius, R in = 10 nm and outer radius, R out = 100 nm. The Raman molecule is modelled as a point dipole source placed at the center of Ω D . The power emitted by the dipole at λ 2 is computed by evaluating Eq. (4) over Γ RE , a circular curve centered on the dipole with radius: For the model domain sketched in Fig. 2(b) the following values are used: The model domain Ω is a rectangle of width, a = W Ω = 150 nm and height H Ω = 600 nm. Ω is truncated from above by perfectly matched layer with a depth of D PML = 300 nm. Ω is truncated from below using a perfect electric conductor boundary condition: n × E = 0. Ω is truncated by a Floquet Bloch boundary condition on the left and right boundaries to model the x-periodicity. The top 500 nm of Ω is taken to contain water (ε r ≈ 1.77). The bottom 100 nm of Ω is taken to contain platinum [63]. The design domain Ω D is placed immediately on top of the platinum region. The design domain Ω D has a width of W Ω D = 150 nm and a height H Ω D = 200 nm. The design domain Ω D has a region fixed to contain water at its center with radius R in = 10 nm. The Raman molecule is modelled as a point dipole source placed at the center of Ω D . The power emitted by the dipole at λ 2 is computed by evaluating Eq. (4) over Γ RE , a straight horizontal line 50 nm below the top boundary of Ω.
For the model domain sketched in Fig. 2(c) the following values are used: The model domain Ω is a square with side length, L Ω = 1200 nm centered at (0 nm,0 nm). Ω is surrounded by a perfectly matched layer with a depth of D PML = 300 nm. The design domain Ω D Metal is centered at (-300 nm,-300 nm). Ω D Metal is a square of side length L Ω D Metal = 200 nm with a circular hole of radius R in = 10 nm at its center. The Raman molecule is modelled as a point dipole source placed at the center of Ω D . The design domain Ω D Dielectric is centered at (-200 nm,-200 nm). Ω D Dielectric is the difference between a square of side length L Ω D Dielectric = 400 nm and Ω D Metal . The input waveguide runs horizontally and is centered at y = −300 nm. It starts at the left edge of Ω and runs until Ω D Dielectric . The output waveguide runs vertically and is centered at x = 300 nm. It starts at the top edge of Ω and runs until Ω D Dielectric . The width of both waveguides is 150 nm. The power emitted by the dipole at λ 2 into the input and output waveguides is computed by evaluating Eq. (4) over Γ in and Γ out respectively. Γ in is a vertical line through the input waveguide at x = −300 nm. Γ 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 (λ/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 Eq. (1) , 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 Eq. (3) 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 λ 1 = 532 nm. The wavelength of the emitted field is taken to be λ 2 = 532 nm. The relative permittivity and relative permeability of air are taken to be ε air = µ air = 1.0. The relative permittivity and relative permeability of water are taken to be ε water = 1.77, µ 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 ε Si 3 N 4 = 2.05, µ water = 1.0 at the operating wavelengths. The speed of light is taken to be c = 3 · 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. ξ initial (r) = 1 ∀ r ∈ Ω 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, η = 0.5. The initial thresholding strength is, β initial = 8. The thresholding strength is increased gradually during the optimization through the values, β = 8, 16, 32. A minimum length scale of all features in the designs is ensured using geometric length-scale constraints with c LS = 625, η e = 0.75, η d = 0.25. The length-scale constraints are only enforced for β = 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. ξ initial (r) = 1 ∀ r ∈ Ω D . The filtering radius used in the smoothing operation is, r f = 10 nm. The thresholding level used is in the thresholding operation is, η = 0.5. The initial thresholding strength used in the thresholding operation is, β initial = 8. The thresholding strength is increased gradually during the optimization through the values, β = 8, 16, 32. A minimum length scale of all features in the designs is ensured using geometric length-scale constraints using c LS = 625, η e = 0.75, η d = 0.25. The length-scale constraints are only enforced for β = 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 Ω D Metal a full metal region is used, i.e. ξ initial (r) = 1 ∀ r ∈ Ω D Metal . As an initial guess for Ω D Dielectric a full dielectric region is used, i.e. ξ initial (r) = 1 ∀ r ∈ Ω D Dielectric . The filtering radius in the smoothing operation is: r f = 10 nm. The thresholding level in the thresholding operation is: η = 0.5. The initial thresholding strength used is in the thresholding operation is: β initial = 8. The thresholding strength is increased gradually during the optimization through the values: β = 8, 16, 32. A minimum length scale of all features in the designs is ensured using geometric length-scale constraints using c LS = 625, η e = 0.75, η d = 0.25. The length-scale constraints is only imposed for β = 32 to allow the design to form freely without limiting topology changes in the earlier stages of the optimization.