1. Introduction

Surface-enhanced Raman scattering (SERS) [1,2] is a common method to increase the sensitivity of Raman spectroscopy, [3] with enhancements of Raman signals reaching $\sim 10^{10}$ for a single molecule [4–7], important for a wide variety of sensing applications [8]. In SERS, Raman-active molecules are placed in the vicinity of a textured surface (e.g. coated with metal nanoparticles [2]) that provides two multiplicative resonant enhancements: it concentrates the incoming pump field of frequency $\omega _1$ at the molecule’s location, and it also Purcell-enhances the emission at a shifted frequency $\omega _2$. In previous work, we derived general upper bounds on the Raman enhancement for arbitrarily-shaped structures given the material’s susceptibility, the size of the scatterer, and the distance to the molecule [9]. Motivated by optimistic results from these bounds, we used topology optimization (TopOpt) [10,11]—in which every “pixel” of a design is a degree of freedom (Sec. 3, below)—to inverse-design novel structures maximizing the Raman enhancement, leading to $\sim 100\times$ improvement over conventional structures [11].

This previous work only analyzed the emission of a single molecule placed at a “hot spot” of maximal electric-field intensity [10,11]. In many practical experiments, however, the molecules are distributed randomly in space, either suspended in a fluid or deposited onto a surface, [12,13] so only a small fraction of the molecules experiences the peak hot-spot enhancement. [14–16] It is an open question to determine what structures maximize average enhancement over all molecule locations. Some authors have analyzed the effect of one or two geometric parameters on averaged enhancement using a simplified metric discussed below [17,18], but neither large-scale optimization (e.g. TopOpt) nor a comprehensive theoretical approach have been developed. Also, additional nonlinear effects arise in UV Raman spectroscopy, where extremely high intensities (“too hot” hot spots) can damage the molecules and quench emission [15,16], but the implications of the effect on optimal design have never been analyzed. Spatially distributed SERS emission is challenging to model rigorously, as it naively requires running a large number of simulations for molecules at different locations, which is especially problematic for inverse design where many structures must be simulated over the course of optimization. Building on a recently developed trace formulation for optimizing incoherent emission processes [19], in this paper we propose an efficient technique for simulation and inverse design of spatially averaged SERS enhancement and analyze its results for TopOpt applied to several example problems addressing effects of hot-spot density, angular selectivity, and nonlinear damage.

In particular, we show that spatially averaged SERS emission in a single direction can be modeled with only two Maxwell solves (i.e., two numerical solutions of the discretized Maxwell equations), one for the pump process and another “reciprocal” solve for the average emission over all molecule locations (Sec. 2.3), easily generalized to support nonlinear damage (Sec. 4.3) and/or anisotropic Raman polarizability. Moreover, this formulation is straightforwardly compatible with large-scale inverse design, requiring only two additional “adjoint” simulations [20] to compute the sensitivity of the output power with respect to “every pixel” of the design (e.g. a material density at every point in TopOpt, Sec. 3). Previous authors employed a simplified $\int \vert \mathbf {E}_1\vert ^{4}$ metric for distributed Raman emission [17,18], where $\mathbf {E}_1$ is the pump electric field, and we show that this is a special case of our framework when the emission is in the same direction as the pump, the Raman shift is negligible ($\omega _1\approx \omega _2$), and the Raman molecule is isotropic. We also analyze how the $\sim \vert \mathbf {E}_1\vert ^{4}$ nonlinearity favors hot spots and field singularities (from sharp corners) in 3D, but less so in 2D (Sec. 2.4).

We apply TopOpt to various example problems in 2D to illustrate the key tradeoffs and physical effects (Sec. 4): normal incidence and emission (Sec. 4.1), $30^{\circ }$ pump and normal emission (which performs nearly as well but with a very different design Sec. 4.2), emission with UV-like nonlinear damage (again leading to very different designs Sec. 4.3), and emission only from a material-surface coating [13] rather than in a volume/fluid coating (Sec. 4.4). By comparing TopOpt for periodic surfaces of varying periods, we observe a best density of the resulting hot spots (Sec. 4.1). We obtain optimized structures that perform about $4\times$ better than coating with optimized spheres or bowtie structures and about $20\times$ better when the nonlinear damage effects are considered. Also, we find similar optimized structures when Raman-active molecules are distributed only on the metal surface, as opposed to throughout a volumetric (fluid) coating. We believe that this framework sets the stage for future work in 3D (where field singularities are stronger), TopOpt for dielectric Raman [21] (instead of metal, trading sharper resonances for weaker localization), and related problems in scintillation detectors (where previous work optimized emission but not absorption [22]).

2. Model formulation

In this section, we provide a general mathematical framework for optimizing spatially averaged SERS enhancement. We begin with the numerical model for Raman scattering (Sec. 2.1) and then show how the trace formulation can be applied to the SERS problem in general (Sec. 2.2). Next, we consider the special case where the Raman signals are received in a single direction (Sec. 2.3). Finally, we provide some analysis of the singularities in the SERS problem (Sec. 2.4).

2.1 Numerical model for Raman scattering

Raman scattering can be modeled as a combination of two electromagnetic processes [2]: first, an incident laser (or equivalent current [23] source $\mathbf {J}_1$) produces an electric field $\mathbf {E}_1e^{-\mathrm {i}\omega _1t}$ at a frequency $\omega _1$. This solves the linear Maxwell equations $M_1\mathbf {E}_1 = \mathrm {i}\omega _1\mathbf {J}_1$, where $M_1$ is the Maxwell (vector Helmholtz) operator $M_1 = \nabla \times \mu _1^{-1}\nabla \times -\frac {\omega _1^{2}}{c^{2}}\varepsilon _1$ with $\varepsilon _1(\mathbf {x})$ and $\mu _1(\mathbf {x})$ being the relative electric permittivity and magnetic permeability at frequency $\omega _1$. Second, a molecule at position $\mathbf {x}_0$ with a Raman polarizability tensor $\alpha$ produces a dipole current density $\mathbf {J}_2 = \alpha \mathbf {E}_1(\mathbf {x}_0)\delta (\mathbf {x}-\mathbf {x}_0)e^{-\mathrm {i}\omega _2t}$ at a frequency $\omega _2$, which produces an emission field $\mathbf {E}_2$ satisfying $M_2\mathbf {E}_2 = \mathrm {i}\omega _2\mathbf {J}_2$, where $M_2$ is the Maxwell operator at the frequency $\omega _2$. The difference $|\omega _2-\omega _1|$ is the Raman shift, and usually $|\omega _2-\omega _1|\ll |\omega _1|$.

Numerically, we discretize this problem (e.g. using finite elements) into a sequence of finite-sized systems of linear equations:

Typically, we are interested in maximizing the power radiated into one or more directions/channels by $\mathbf {u}_2$ for a given incident source $\mathbf {b}_1$. This can be expressed as quadratic functions of the emission fields $\mathbf {u}_2$ via the Poynting flux. Since the power is always a real-valued quantity, it corresponds in particular to a Hermitian quadratic form

When the Raman-active molecules are distributed randomly in some region, one needs to solve for the emission field $\mathbf {u}_2$ for every single molecule (different $\alpha$ and ${A}$) and then take the the average:

2.2 Trace formulation for Raman scattering

The trace formulation [19] is motivated by previous works on thermal emission, luminescence, and related problems [25–27], where methods were developed to model incoherent emission from many molecules in a large volume as a single matrix trace operation, rather than individual matrix solves for every emitter point. Here, we briefly introduce this approach and show how it is applied in the spatially-averaged SERS problem.

The key idea is to rewrite our scalar objective Eq. (3) as a “$1\times 1$” trace, and then employ the cyclic-shift trace property [28] to group the ${A}\mathbf {u}_1$ terms together:

The simplest case is that of isotropic (scalar) Raman polarizability $\alpha$, in which case $\langle \alpha (\mathbf {x})\mathbf {E}_1(\mathbf {x})\mathbf {E}_1^{\dagger }(\mathbf {x})\alpha ^{\dagger }(\mathbf {x})\rangle _\alpha = \langle |\alpha (\mathbf {x})|^{2}\rangle _\alpha \mathbf {E}_1(\mathbf {x})\mathbf {E}_1^{\dagger }(\mathbf {x})$. Defining $|\alpha _0(\mathbf {x})|^{2} = \langle |\alpha (\mathbf {x})|^{2}\rangle _\alpha$ as the mean-square polarizability at each point (i.e. the Raman polarizability multiplied by the probability of the molecule being at that point), we obtain $C_\alpha (\mathbf {x},\mathbf {x}^{\prime })=|\alpha _0(\mathbf {x})|^{2}\mathbf {E}_1(\mathbf {x})\mathbf {E}_1^{\dagger }(\mathbf {x}^{\prime })\delta (\mathbf {x}-\mathbf {x}^{\prime })$ and consequently:

For the more general case of an anisotropic Raman polarizability tensor $\alpha$, the expression $\langle \alpha (\mathbf {x})\mathbf {E}_1(\mathbf {x})\mathbf {E}_1^{\dagger }(\mathbf {x})\alpha ^{\dagger }(\mathbf {x})\rangle _\alpha$ must be averaged over all possible orientations of the molecule, corresponding to a average of $Q \alpha Q^{T} \mathbf {E}_1 \mathbf {E}_1^{\dagger }(\mathbf {x}) Q^{T} \alpha ^{\dagger } Q$ over all possible $3\times 3$ rotation matrices $Q$. If all orientations are equally likely, this rotation average can be computed analytically with the help of formulas derived in Ref. [30].

Once the correlation matrix ${B}$ is determined, we can then apply different techniques developed in our previous work [19] to combine the trace estimation problem with the TopOpt for different scenarios depending on the number of input and output channels. In this paper, we focus on the case where the emitted Raman signals are received in a single direction/channel, which means the objective matrix ${O}$ now becomes rank 1, and the trace formulation also simplifies to two Maxwell solves, one forward and one reciprocal, as discussed below.

2.3 Single-channel simplification

For spatially incoherent Raman emission into a single direction/channel, the average power of all the emitters can be computed with a single “reciprocal” solve. This was derived in a very general setting by our previous work [19] and is closely related to the well-known Kirchhoff’s law of thermal radiation (reciprocity of emission and absorption) [31] as well as analogous results for light-emitting diodes [32] or scintillation [22]. In this section, we apply the algebraic framework of Ref. [19] to the specific case of Raman emission.

The power emitted into a single direction can be expressed as a mode overlap integral of the emitted field $\mathbf {u}_2$ and a planewave mode $\mathbf {o}$, which is algebraically of the form $\Vert \mathbf {o}^{\dagger }\mathbf {u}_2\Vert ^{2}$ [33]. In terms of the quadratic form Eq. (2), the objective matrix $O$ is now simply a rank-1 matrix ${O}=\mathbf {o}\mathbf {o}^{\dagger }$ [19], and the objective trace Eq. (4) reduces to

The formulation can be further simplified when the Raman-active molecule is isotropic. Inserting Eq. (7) into Eq. (8), we obtain

From Eq. (9), we can also see that the equation further simplifies if (i) the pump and emission directions are the same and (ii) we make an approximation of a negligible Raman shift ($\omega _1\approx \omega _2$), in which case one can take $\mathbf {E}_2^{\prime } \approx \mathbf {E}_1$. For isotropic Raman polarizability whose mean $|\alpha _0|^{2}$ is constant in some volume $V$ and zero elsewhere, this leads to the $\int _V \vert \mathbf {E}_1\vert ^{4}$ figure of merit used in several previous works [17,18] for Raman power (which is often presented heuristically, but has also been justified using reciprocity [34]).

2.4 Corner singularities and hot spots

It is common knowledge that SERS tends to favor geometries with “hot spots” where high field intensities arise from geometric singularities such as sharp tips/cusps, bowtie antennas, or touching spheres, especially for single-molecule SERS where the enhancement theoretically diverges in the limit of arbitrarily sharp tips (in the continuum macroscopic Maxwell equations with local materials). However, it is less clear whether the average Raman enhancement of many volume-distributed emitters still favors such hot spots, since the effect of a field singularity might be spatially averaged out. In this section, we analyze the effect of corner singularities on average Raman enhancement in 2D and 3D for the single-channel isotropic-Raman case of Eq. (9).

Field singularities at sharp corners are frequency-independent and can be analyzed purely using electrostatics [35], so the pump $\mathbf {E}_1$ and reciprocal $\mathbf {E}_2'$ fields have identical scaling near a sharp tip. Therefore, without loss of generality, we can analyze the simplified metric $\int \vert \mathbf {E}_1\vert ^{4}$ in the neighborhood of a sharp tip.

For a 2D sharp corner in a dielectric or metallic material $\varepsilon$ enclosing an angle $\phi < \pi$, the field singularity of the field $\mathbf {E}_1$ is a fractional power law in the distance $r$ from the tip [35]:

In 3D, the field singularity at sharp tips (e.g. cones or corners) is stronger than in 2D. For example, the fields at the tip of a 3D cone with angle $\phi < \pi$ also exhibit a singularity $\mathbf {E}_1\sim r^{t-1}$ but with a stronger power law $0 < t < 1$ (e.g. $t=\frac {1}{2\log (8/\phi )}>0$ for perfect conductors) [37]. The integral then becomes $\int r^{4(t-1)} r^{2}\, dr \sim r^{4t-1}$, which diverges for $t<1/4$ (sufficiently small $\phi$). Therefore, we expect that 3D topology optimization of incoherent Raman emission will favor arbitrarily sharp tips, limited only by the imposition of manufacturing constraints [36].

It is worth contrasting the Raman case, in which the field singularity is squared by the conjunction of pump and emission enhancement, with spontaneous emission in cases with non-optical pumping, such as light-emitting diodes [32], scintillation from high-energy particles [22], or thermal emission. In such cases, if the excitation is nearly uniform in the vicinity of a sharp tip, then the emitted power scales as $\int |\mathbf {E}_2'|^{2}$ from the reciprocal field alone. The contribution of a corner singularity is then $\sim \int r^{2(t-1)}r\,dr \sim r^{2t}$ in 2D ($t > 1/2$) and $\sim \int r^{2(t-1)}r^{2}\,dr \sim r^{2t+1}$ in 3D ($t > 0$), both of which vanish as $r\to 0$. In consequence, one does not expect arbitrarily sharp corners/tips to be favored when optimizing spatially averaged emission alone or a spatially averaged local density of states, LDOS. Indeed our previous work on topology optimization of incoherent emission [19] or scintillation [22] did not exhibit arbitrarily sharp corners, in contrast to the “bowtie antenna” singularities that typically arise when optimizing emission/LDOS from a single emitter location [38,39]. Similar considerations apply to optimization of photovoltaic cells, since maximizing absorption is equivalent to maximizing spatially averaged emission via Kirchhoff’s law.

3. Density-based topology optimization

In this section, we briefly review the technique of density-based TopOpt [40], which is used to solve the inverse design problem of tailoring the surface geometry to maximize our Raman-power objective from the previous section.

In density-based TopOpt, a continuous design field (density) $\rho (\mathbf {x}) \in [0,1]$ is defined on the spatial “design” domain. This design field is first passed through a smoothing filter to regularize the optimization problem that sets a filter lengthscale $r_f$, as otherwise one may obtain arbitrarily fine features as the spatial resolution is increased. (Additional steps are required to impose strict manufacturing constraints [36]). The smoothing convolves $\rho$ with a low-pass filter to obtain a smoothed density $\tilde {\rho }$ [40]. There are many possible filtering algorithms, but in a finite-element method (FEM) setting, especially with complicated nonuniform meshes, it is convenient to perform the smoothing by solving a simple “damped diffusion” PDE, also called a “Helmholtz” filter [41]:

Next, one employs a smooth threshold projection on the intermediate variable $\tilde {\rho }$ to obtain a “binarized” density parameter $\tilde {\tilde {\rho }}$ that tends towards values of $0$ or $1$ almost everywhere [42]:

Finally, one introduces a material, described here by an electric relative permittivity (dielectric constant) $\varepsilon (\mathbf {x})$ in the Maxwell operator $M_1$ or $M_2$, given by:

Numerically, we employ a recent free/open-source FEM package, Gridap.jl [44], in the Julia language [45], which allows us to code highly customized FEM formulations efficiently in a high-level language. We discretized $\rho$ and $\{\tilde \rho,\tilde {\tilde \rho }\}$ with piecewise-constant (0th-order) and first-order elements, respectively. During optimization, one must ultimately compute the sensitivity of the objective function with respect to the degrees of freedom $\rho$. For each step outlined above (smoothing, threshold, PDE solve, etc.), we formulate a vector–Jacobian product following the adjoint method for sensitivity analysis [20] with some help from automatic differentiation [46]. Then these are automatically composed (“backpropagated”) by an automatic-differentiation (AD) system [47]. In this way, the gradient with respect to all of the degrees of freedom ($\rho$ at every mesh element) can be computed with only two additional adjoint Maxwell solves [20].

4. Results

In this section, we present various example problems building from a single-channel framework (Sec. 2.3), illustrating the key tradeoffs and physical effects. We begin with the simplest case with normal incidence and emission (Sec. 4.1), revealing a best density of the hot spots. Next, we show that the pump and emission angles can also be considered as design parameters. We provide an example where the pump field is fixed at angle $\theta _1=30^{\circ }$ and we search for the best emission angle, which turns out to be roughly normal emission and performs nearly as well as the normal-incidence pumping but with a very different design (Sec. 4.2). Then we take into consideration the effect of UV-like nonlinear damage [15] quenching the emission and optimize with this taken into account, which again leads to very different designs (Sec. 4.3). Finally, we briefly discuss the case where SERS emission is only from molecules coating the material surface rather than being distributed throughout the volume of a fluid (Sec. 4.4).

Figure 1 is a sketch of the single-channel SERS design problem in 2D. The Raman-active molecules are distributed uniformly in a fluid (water) background above a periodically patterned metal (silver) surface with period $L$. Since the Raman signal is proportional to the density of the molecules (indicated by the Raman polarizability $\alpha _0(x)$), we chose $|\alpha _0(x)|^{2}=1-\tilde {\tilde {\rho }}$ for simplicity. An incident planewave ($H_z$-polarized, $\lambda _1=532$ nm) at angle $\theta _1$ excites the molecules, and the Raman-shifted power ($\lambda _2=549$ nm) at angle $\theta _2$ is measured and optimized. The design region is an $L\times H$ (200 nm) rectangular domain, in which the material can either be fluid (i.e. with molecules) or metal (i.e. without molecules). We sweep the optimization over different periods $L$ to find a best period, which corresponds to a best density of hot spots. An infinitely thick layer of molecules would emit infinite power in the absence of water absorption and pump depletion, but since we are only interested in optimizing near-field enhancement we limit the Raman-active molecules to a half-wavelength layer of thickness $H_\mathrm {R}=(\lambda _1+\lambda _2)/(4n_\mathrm {f})$ above the design domain.

 

Fig. 1. Sketch of the 2D Raman scattering design problem. The Raman-active molecules are distributed uniformly in fluid (water, $n_\mathrm {f}=1.33$) background near a periodically patterned metal (silver) surface with period $L$. The incident planewave ($H_z$-polarized, $\lambda _1=532$ nm, green) at angle $\theta _1$ excites the molecules, and the Raman-shifted signal ($\lambda _2=549$ nm, blue) at angle $\theta _2$ is measured and optimized. A half-wavelength molecular layer of thickness $H_\mathrm {R}=(\lambda _1+\lambda _2)/(4n_\mathrm {f})$ is placed above the design domain of height $H=200$ nm.

Download Full Size | PDF

4.1 Normal incidence and emission

In this example we consider the case when the incident pump and measured emission are both normal to the surface ($\theta _1=\theta _2=0^{\circ }$). We maximize this power $P$ of the emitted Raman signal at different periods $L$, normalized to a baseline power for a flat metal surface (metal half-space). Since this optimization problem is non-convex, TopOpt may easily converge to different local optima from different initial geometries [20]. In Fig. 2, we plot only the local optima with largest power we found for each period (from 100 nm to 500 nm with a sampling of 50 nm spacing) for 20 different random starting structures. We find that the largest local optima have very similar performance, within $\sim 10$%, giving us some confidence that there are unlikely to be dramatically better local optima yet to be found.

As shown in the inset of Fig. 2, the optimized patterns share a small “notch” feature that creates a hot spot (localized resonance), and some of them exhibit spontaneous symmetry breaking: the resulting pattern is asymmetric although the problem is mirror-symmetric [48]. This “notch” hot spot is different from the field singularity arising from sharp corners where the field theoretically diverges. Here, the minimum length-scale of those notches is about 10 nm (note that this is 2 times of the filter size of 5 nm, thus the notch size is not constrained by filter size). The best density of those hot spots is found to be at period $L=400$ nm, which is nearly the wavelength $\lambda _\mathrm {vacuum} \sqrt {\frac {1}{\varepsilon _\mathrm {f}} + \frac {1}{\varepsilon _\mathrm {metal}}} = 375$ nm of surface plasmons at a flat silver–air interface [49]. We also did a similar optimization in the UV regime with $\lambda _1=400$ nm and $\lambda _2=437$ nm (not shown here), and the best period $L=300$ nm was also found to be close to the surface-plamson wavelength $241$ nm. Intuitively, periodic perturbations with this wavelength implement a grating coupler between normal-incident radiation and a surface-plasmon resonance [49,50], but of course the period changes as the surface is deformed substantially. We also performed an optimization for a doubled period of $L=800$ nm (not shown), and unsurprisingly found that it converges to two hot spots in each unit cell with similar performance to the single-hotspot $L=400$ nm design.

 

Fig. 2. Emitted power of Raman-active molecules situated on the optimized surfaces, compared to that of a flat surface, over different periods for silver at normal incidence pump ($\lambda _1=532$ nm, $\theta _1=0^{\circ }$) and emission ($\lambda _2=549$ nm, $\theta _2=0^{\circ }$). The insets show the optimized patterns and the pump fields $\vert E_1\vert$ (normalized to incident planewave) in a unit cell for typical periods.

Download Full Size | PDF

To gain a better understanding of how the optimized surface performs, we also compare to a surface coated with optimized spheres – optimized over both the sphere diameter and period. As shown in Fig. 3(a), the optimized spheres are of diameter 310 nm and have period equal to their diameters (i.e. touching). The performance of our TopOpt surface is about $4\times$ better than that of the optimized spheres. A similar optimization was carried out for a bowtie structure, and the optimized bowtie performs slightly worse than the optimized spheres. From the pump field displayed below the pattern in Fig. 3, we can also see that the notch hot spot in the optimized surface is more spread out than the singular hot spot produced by tangent spheres—as predicted in Sec. 2.4, 2D distributed SERS emission does not favor singularities at points or cusps.

 

Fig. 3. Comparison of a) optimized spheres and b) optimized surface for the normal pumping and emission case. The optimized spheres are found to have a diameter of 310 nm and are adjacent to each other. The optimized surface is found to have a period of $L=400$ nm. The pumping fields $\vert E_1\vert$ are displayed below the pattern. The reciprocal fields $\vert E_2^{\prime }\vert$ are similar because of the small Raman shift.

Download Full Size | PDF

4.2 Oblique incidence

In this example (Fig. 4), we consider the case where the incident pump field is fixed at an angle $\theta _1=30^{\circ }$ and we search for a best surface and emission angle $\theta _2$ to maximize the power (per unit cell) of the emitted signal for a fixed period of $L=400$ nm.

 

Fig. 4. Optimized surface with a fixed incident angle at $\theta _1=30^{\circ }$. a) Pattern of the optimized surface for the best emission angle: normal emission ($\theta _2=0^{\circ }$). b) Emitted power (per unit cell) of the optimized surface with $\theta _1=30^{\circ }$ as a function of the emission angle $\theta _2$, normalized by the normal emission from a flat surface. c) Reciprocal fields $\vert \mathbf {E}^{\prime }_2\vert$ at selected emission angles.

Download Full Size | PDF

In Fig. 4(a), we show the optimized surface, which we find to be best by scanning the power over all emission angles $\theta _2$ (Fig. 4(b)). As it happens, the best emission angle for this structure is nearly 0$^{\circ }$. The structure again has small notches that create hot spots, but with a very different design compared to the normal-incident pump in the previous section. The power near $\theta _2=-30^{\circ }$, inverse to the incident direction, is also large, which is expected since the reciprocal field is similar (for small Raman shift) to the pump field when $\theta _1 = -\theta _2$. On the other hand, the emitted power at $\theta _2=30^{\circ }$ is very low. From Fig. 4(c), we can see that this occurs because the reciprocal fields $\mathbf {E}_2^{\prime }$ excite a localized resonance for both $\theta _2=-30^{\circ }$ and $0^{\circ }$ but not at $\theta _2=30^{\circ }$. After all, the resonant frequencies of a periodic surface depend on angle, corresponding to a Bloch wavevector.

4.3 Nonlinear damage

In this example, we consider how to design surfaces that enhance spatially-averaged SERS spectra without creating very strong localized fields that will damage molecules and quench emission, a familiar nonlinear phenomenon observed experimentally, especially for UV SERS spectroscopy [15,16]. Here we assume a threshold magnitude $\vert \mathbf {E}_\mathrm {th}\vert$, above which the pump fields will damage the molecules and quench emission. We model this phenomenon by treating the mean Raman susceptibility $|\alpha _0|^{2}$ in Eq. (9) as nonlinear, exponentially decreasing for pump fields larger than the threshold. We replace $|\alpha _0|^{2}$ with

Figure 5(a) shows the optimized surface for a nonlinear damage threshold $\vert E_\mathrm {th}\vert =10$ and a pump planewave with $\vert E\vert =1$, assuming normal incidence at $\lambda _1=532$ nm and normal emission at $\lambda _2=549$ nm. We can see that the pump-field pattern is much more spatially spread out than the highly concentrated hot spots of the previous sections. For comparison, if we use the previous structure optimized without a threshold constraint from Sec. 4.1, we find that its performance rapidly degrades in the presence of nonlinear damage. This is shown in Fig. 5(b): the emitted power falls rapidly with decreasing damage threshold for the surface and spheres optimized without this threshold because it is easy for the pump field to damage the molecules in their hot spots. However, the emitted power from the design with nonlinear damage taken into account remains constant for powers above the design threshold $\vert E_\mathrm {th}\vert =10$ and only drops for thresholds lower than this. Moreover, it performs about $20\times$ better than the optimized spheres at $\vert E_\mathrm {th}\vert =10$. Because for the optimized spheres, a large region around the hot spots above the threshold will go dark and the fields change more rapidly. However, for the optimized structure with threshold, the fields are more evenly distributed below the target threshold.

 

Fig. 5. a) Pattern of the optimized surface with target threshold $\vert E_\mathrm {th}\vert =10$ for a pump planewave with $\vert E\vert =1$ and the corresponding pump field $\vert E_1\vert$. b) Power from the prior normal incidence surface (blue) and spheres (black) optimized without threshold, and the new surface optimized with $\vert E_\mathrm {th}\vert =10$ (orange), as a function of the damage threshold $\vert E_\mathrm {th}\vert$.

Download Full Size | PDF

4.4 Surface emitters

In some experiments, there is a monolayer of Raman-active molecules deposited on the metal surface [13,15], which increases their exposure to the resonant enhancement. In this case, the average emission should be computed as a integral over the surface instead of the volume, as in the examples above. One simple technique to model this is to make the mean Raman polarizability $|\alpha _0|^{2}$ proportional to $\tilde {\tilde {\rho }}(1-\tilde {\tilde {\rho }})$, which is $\approx 0$ except near the metal surface. (Here, $\tilde {\tilde {\rho }}$ is the thresholded and filtered density field from Sec. 3. More rigorously, one can employ a double filtering technique to achieve exact identification of the surface [51–53]). This allows us to take the surface geometry into account during sensitivity analysis and optimization.

Figure 6 shows the optimized surface for spatially-averaged SERS emission and the corresponding pump fields. We can see that the surface also has “notch” structures, but with a smaller minimum lengthscale of about 5 nm compared to the volume-averaged TopOpt surface in Sec. 4.1, and performs about $3\times$ better than the optimized spheres when averaged over the surface. One thing to be noted is that the surface integral could diverge in principle at sharp corners, as can be seen from the $\int \vert \mathbf {E}_1\vert ^{4}$ metric (Sec. 2.4). Here, the design does not seem to exhibit sharp corners, probably because we imposed a soft minimum lengthscale by setting a smoothing filter radius of $r_f=2$ nm. However, one might ultimately need to impose manufacturing constraints and/or nonlinear damage thresholds to prevent inverse design from favoring singular structures for surface emission, even in 2D.

 

Fig. 6. Optimized (silver) surface for surface-average instead of volume-average at period $L=400$ nm, for an pump field at $\lambda _1=532$ nm and emission field at $\lambda _2=549$ nm.

Download Full Size | PDF

5. Conclusion

We presented a general framework for optimizing spatially-averaged SERS enhancement, requiring only two Maxwell solves per optimization step for emission in a single direction. We explored this technique with a number of 2D examples to illustrate the computational technique as well as the basic phenomena of best hot-spot densities, angular dependence, and the effect of nonlinear damage. The next step is to carry these techniques into 3D, where the same computational principles apply but radically different structures may arise due to the stronger singularities at sharp 3D tips (Sec. 2.4). In 3D, these singularities mean that the imposition of manufacturing constraints [36] and/or nonlinear thresholding will play a key role. Because only a small range of geometries have previously been explored for this problem, it is possible that substantial practical improvements may be uncovered by TopOpt for 3D SERS, especially in less-explored circumstances such as distinct input/output directions, nonlinear damage, or even integrated SERS with waveguide channels [10].

Another important complementary problem is the development of theoretical upper bounds to distributed SERS emission, generalizing earlier work bounding emission at a single location [9], as well as related efforts to bound the “density” of resonant modes (e.g. for solar cells [54]). Computationally, there are a wide variety of nonlinear-optics problems that may potentially be optimized using techniques involving coupled linear Maxwell solves, from scintillation processes [22] to harmonic generation [55]. Balancing the tradeoffs between multiple physical processes is precisely where large-scale optimization has the greatest advantages over intuitive human design.