The Design and Optimization of Plasmonic Crystals for Surface Enhanced Raman Spectroscopy Using the Finite Difference Time Domain Method

We present computational studies of quasi three-dimensional nanowell (NW) and nanopost (NP) plasmonic crystals for applications in surface enhanced Raman spectroscopy (SERS). The NW and NP plasmonic crystals are metal coated arrays of cylindrical voids or posts, respectively, in a dielectric substrate characterized by a well/post diameter (D), relief depth (RD), periodicity (P), and metal thickness (MT). Each plasmonic crystal is modeled using the three-dimensional finite-difference time-domain (FDTD) method with periodic boundary conditions in the x- and y-directions applied to a computational unit cell to simulate the effect of a periodic array. Relative SERS responses are calculated from time-averaged electric field intensity enhancements at λexc and λscat or at λmid via GSERS4=g2(λexc)×g2(λscat) or Gmid4=g4(λmid), respectively, where g2=|E|2/|E0|2. Comparisons of GSERS4 and Gmid4 are made to previously reported experimental SERS measurements for NW and NP geometries. Optimized NW and NP configurations based on variations of D, P, RD, and MT using GSERS4 are presented, with 6× and 2× predicted increases in SERS, respectively. A novel plasmonic crystal based on square NP geometries are considered with an additional 3× increase over the optimized cylindrical NP geometry. NW geometries with imbedded spherical gold nanoparticles are considered, with 10× to 103× increases in SERS responses over the NW geometry alone. The results promote the use of FDTD as a viable in silico route to the design and optimization of SERS active devices.


Introduction
Surface enhanced Raman spectroscopy (SERS) employing metallic nanostructures is a well studied phenomenon that arises predominantly from large electric field intensity enhancements when light is incident on a metal that supports surface plasmons (SPs) [1][2][3]. SPs are collective oscillations of conducting electrons with a characteristic frequency that can exist at the interface of a dielectric and a metal [4,5]. Under certain conditions, incident light resonant with the plasmon frequency can excite the SP within the metal, giving rise to sub-wavelength field confinement and intensity enhancements, defined hereafter as g 2 =|E| 2 /|E 0 | 2 , where |E| and |E 0 |are the magnitudes of the total and incident electric field, respectively. (The spatial and frequency/wavelength dependence of E is ignored for the moment for simplicity.) The SP frequency is influenced by the size, shape, and local dielectric environment of the metal and can thereby be tuned as desired [6,7]. Molecules that traverse the SP field experience not only enhanced absorption of light proportional to the enhancement factor at the excitation frequency, g 2 exc , but also enhanced Raman scattering proportional to the enhancement factor at the Raman scattered frequency, g 2 scat . The overall SERS response is then proportional to g 2 exc × g 2 scat [8]. Early SERS substrates were composed of rough metal surfaces or colloidal solutions of metallic nanoparticles whose sensitivity and reproducibility suffered from the generation of only random "hot-spots" on the metal surface or in the temporal gaps between adjacent particles. Within the last decade, nanolithography methods based on soft imprinting have enabled the inexpensive fabrication of robust SERS substrates that give rise to uniform hot-spots over large areas, with enhancement factors on the order of 10 6 to 10 8 [9][10][11][12].
The SERS substrates considered in this work are based on quasi-three-dimensional plasmonic crystals composed of dielectric supports patterned with either voids or posts upon which thin films of gold are deposited using electron beam evaporation. We do not consider sputtering, which can give rise to wall features. Nanowell (NW) geometries (Figure 1a) result from voids in the supporting dielectric [13], and nanopost (NP) geometries (Figure 1b) result from posts on the supporting dielectric [14]. Previous experimental studies of these NW and NP geometries involved fabricating an array of plasmonic crystals, each region of the array characterized by a particular diameter, D, and periodicity, P. The SERS response of a target molecule, namely benzene thiol, was measured for each region of the array. A finite-difference time-domain (FDTD) analysis of the geometries studied showed qualitative agreement (more on this below) between relative SERS responses, and a maximum in the SERS response was attributed to the excitation of a local surface plasmon resonance (LSPR) of the well or post. The present work aims to re-examine these plasmonic crystal studies as a means to compare the predictive abilities of two related FDTD approaches in estimating relative SERS responses. The FDTD method is a numerical approach in which electric and magnetic field components are propagated on a set of discrete grid points, each described by a dielectric constant consistent with the material in the structure being modeled at that point [15]. FDTD is a powerful tool for modeling the interactions of light with arbitrarily shaped, complex geometric structures and has been used in a wide variety of surface plasmon applications, including sensing, photovoltaics, imaging, optical trapping, metamaterials and optical transparency, and, of course, SERS [16][17][18][19][20][21][22][23][24]. With FDTD electric field intensities calculated at both excitation and scattered wavelengths, λ exc and λ scat , respectively, the SERS response can be calculated as where x i , y j , and z k are FDTD grid points and N is the number of grid points in the volume of air near the surface of a metallic nanostructure [23,25,26]. Equation (1) requires the computation of fields at λ exc and λ scat . For large three-dimensional grids, due either to a large spatial extent or a fine grid resolution, or for long simulation times, it is possible to reduce computational cost by approximating the SERS response as a function of a single wavelength, typically λ exc or λ mid , the midpoint between the excitation and scattered wavelenghts [13,14,27,28]. The use of λ exc is suitable if the Stokes shift is small compared with λ exc . The use of λ mid is suitable if the surface plasmon resonance is sufficiently broad over the range of excitation and scattered wavelenghts. Early studies of colloidal gold nanospheres indeed showed that particles with an LSPR parked between the excitation and scattered wavelengths were able to effectively enhance both λ exc and λ scat to produce the highest SERS response [29]. The SERS response based on this approximation using λ mid , referred to hereafter as G 4 mid , is given by The aim of this study, then, was to compare Equations (1) and (2) in estimating relative SERS responses for NW and NP plasmonic crystals, which are able to support multiple LSPRs with complicated coupling. In Section 2.1, we show that the former approach is necessary to get more quantitative relative SERS responses. With this in mind, we then set out to design and optimize SERS substrates in silico without fabrication costs. In Refs. [13,14], the dimensions studied experimentally do not allow one to distinguish between the effect of diameter and the effect of periodicity on the SERS response, both of which can affect a SP resonance. In Section 2.2 , we optimize of NW and NP geometries with respect to well depth/post height, diameter, and periodicity using Equation (1). In Section 2.3, we present two studies looking at novel plasmonic crystals based on variations of NW and NP geometries. Differences include nanoposts with square cross-sections ( Figure 1c) and nanowells with imbedded nanoparticles (Figure 1d).

Comparison of FDTD Simulated SERS Responses for Nanowell and Nanopost Plasmonic Crystals
Time-averaged electric field intensity enhancements were calculated for λ = 785, 821, and 857 nm using the finite-difference time-domain (FDTD) method for a series of nanowell (NW, Figure 1a) and cylindrical nanopost (NP, Figure 1b) geometries. For the NW geometries, an SU-8 support with refractive index n = 1.59 was used with a relief depth RD = 360 nm and gold metal thickness MT = 40 nm [13]. For the NP geometries, an NOA support with refractive index n = 1.56 was used with a relief depth RD = 200 nm and gold metal thickness MT = 24 nm [14]. All parameters were consistent with fabricated arrays. Periodic boundary conditions were implemented in the xand y-directions to simulate the array. Periodicities and well/post diameters considered are listed in Table 1.
SERS responses G 4 SERS and G 4 mid were calculated using Equations (1) and (2), respectively.  Experimental Raman intensity measurements reveal a maximum at a diameter D = 456 nm (P = 730 nm) at 55,000 counts. Experimental SERS spectra were collected for a 15 mM solution of benzene thiol using a dispersive Raman microscope. The units are not directly comparable to G 4 SERS and G 4 mid , which are based on unitless electric field intensity enhancements. The FDTD SERS responses, therefore, are scaled such that the maximum experimental and theoretical values coincide. It is clear from Figure 2a that G 4 SERS , using Equation (1), reproduces the experimental SERS responses better than G 4 mid , using Equation (2). Differences between the G 4 SERS and G 4 mid results can be understood in terms of the exictation of different LSPRs for D = 456 nm and D = 514 nm at 785 nm, 821 nm, and 857 nm, and predicted SERS response based on λ mid alone can be misleading.
In a similar comparison, Figure 2b contains a plot for the NP geometries comparing the FDTD simulated SERS responses (dashed lines with open symbols) with experimental measurements (solid line with filled circles) reported in Ref. [14]. Experimental Raman intensity measurements reveal a maximum at a diameter D = 224 nm (P = 584 nm) at 9000 counts, and the FDTD SERS responses are again scaled such that the maximum experimental and theoretical values coincide. (It should also be noted that the units for the experimental SERS response for the NWs and the NPs are not exactly the same, and it turns out that NP arrays produce higher SERS responses than similar NW geometries [14].) While the agreement between G 4 SERS and G 4 mid are more similar in terms of predicted SERS responses, the G 4 SERS values are slightly better in a least squares sense. The better agreement in the NP case can be understood as less coupling possible between the disc and the metal film, in contrast to the NW case where the void, film, and lower disc can couple, as evidenced by the differences in the transmission spectra seen in Figure 3.

Optimization of FDTD Simulated SERS Responses for Nanowell and Nanopost Plasmonic Crystals
The better quantitative agreement between G 4 SERS and G 4 exp allows us to optimize relative SERS responses for the NW and NP geometries by scanning geometry parameter space, allowing the relief depth (RD), diameter (D), periodicity (P), and metal thickness (MT) to vary and calculating the resulting relative SERS response. We define the NW geometry with D = 456 nm, P = 730 nm, RD = 360 nm, and MT = 40 nm as the NW "control". We also define an optimization factor such that any O.F. > 1 results in an enhancement of the expected SERS response. Figure 4 contains a plot of optimization factors as the nanowell RD, D, P, and MT were varied sequentially. The final optimization factor is 6.0, indicating a nanowell plasmonic crystal with P = 730 nm, D = 500 nm, RD = 160 nm, and MT = 70 nm would give rise to a 6× increased SERS signal over the experimentally optimum structure. Figure 5 shows a comparison of electric field enhancements at λ = 785 nm and λ = 857 nm for the control (Figure 5a) and optimized (Figure 5b) arrays, respectively. Increased field intensities in the well region lead to a higher SERS response.

Optimization of FDTD Simulated SERS Responses for Novel Plasmonic Crystals
In this section, we describe calculations aimed at increasing relative SERS responses even further by considering new geometries that should support either increased field intensities and/or additional coupling at 785 nm and 857 nm. To this end, we considered a square nanopost (Figure 1c) array and a nanowell array with imbedded nanoparticles (Figure 1d). The motivation for considering the square post was to excite local surface plasmons at the sharper corners, and such an array could easily be fabricated using the same soft-nanolithography techniques used to construct the cylindrical nanoposts considered above. The square NPs were positioned such that edges of length D were oriented along the xand y-directions. In order to couple into surface plasmons at the corners, we used incident light polarized at 45 o . We considered the following range of values: D = 130 − 500 nm, RD = 150 − 260 nm, P = 500 − 800 nm, and MT = 6 − 20 nm. An optimum square NP was found to have an optimization factor of O.F. = 6.3 with D = 150 nm, RD = 190 nm, P = 730 nm, and MT = 24 nm. Figure 7 depicts the structure and fields at 785 nm. (Field distributions at 857 nm were similar to 785 nm and are not presented here.) For the plasmonic crystals considered thus far, optimum SERS responses are associated with the coupling of local surface plasmons of the film and either of the void (NW) or post (NP), as indicated by periodicity and shape affect on the electric field intensities. We also considered imbedding spherical gold nanoparticles in NW geometries. The aim of this study was to increase electric field intensities by allowing coupling between the particle and the perforated film as well as between the particle and the well disc. Beginning with the optimum NW geometry from Section 2.2, we calculated G 4 SERS for a range of particle diameters, from 80 nm to 300 nm. Optimization factors are plotted in Figure 8. The control NW geometry is the same as used in Section 2.2, namely D = 456 nm, P = 730 nm, RD = 360 nm, and MT = 40 nm. It is not surprising that G 4 SERS (and hence O.F.) increases with increasing particle diameter. While isolated solid gold nanospheres do not absorb strongly in the near infrared, scattering increases for increasing particle size. Figure 9 shows the extinction efficiency (extinction cross section per geometrical cross section) for a 300 nm gold sphere. A broad peak is seen above 700 nm. The increased scattering as well as the coupling to the nanowell gives rise to increasingly large enhancements.

Discussion
In this work, we set out to use the FDTD method to calculate electric field intensity enhancements at λ exc , λ mid , and λ scat for a series of plasmonic crystals to compare G 4 SERS and G 4 mid , as defined in Equations (1) and (2), in their ability to reproduced experimental measurements. By using the word relative, we mean that maximum G 4 SERS and G 4 mid are scaled to match maximum experimental values. The scaling translates the unitless electric field intensity enhancement, g 2 , into units consistent with experimental SERS measurements, and as structure parameters are varied, increases or decreases in the SERS response can be predicted. For both geometries, G 4 SERS understandably gives rise to better quantitative agreement with experimental SERS measurements. For the NP geometry, G 4 mid is a good approximation to G 4 SERS due to the broad spectral feature in the λ exc -λ scat range ( Figure 3). However, for the NW geometry which supports more complex coupling between multiple SP resonances in this range, G 4 mid failed to produce even qualitative agreement for some arrays. These results are not surprising; after all, G 4 mid is an approximation, and when underlying assumptions break down, the accuracy of an approximation also suffers. However, they do suggest caution in using G 4 mid for predictive applications involving complex nanostructures for SERS. It is important to note that the better quantitative agreement of G 4 SERS comes at a computational cost given simulations at both excitation and scattering wavelengths must be considered.
It is also important to note that alternate strategies could be used to estimate SERS responses. While Equations (1) and (2) represent an average over all grid points corresponding to air in the vicinity of the metal surface, one could envision using only maximum values of G 4 mid or G 4 SERS or averaging over only values above some threshold. While we do not present the results here, we also performed estimations of SERS responses based on these two alternate approaches and found that neither case reproduced relative experimental as well as an average over all grid points in a volume of air near each NW or NP.
We also set out to design and optimize plasmonic crystals in silico by varying structure parameters, such as relief depth (RD), well/post diameter (D), periodicity (P), metal thickness (MT), and post shape to maximize excitation and/or coupling of local surface plasmons. While only modest improvements were seen for the cylindrical nanopost, we saw a 6× increase in the NW geometry and the square NP geometry over their respective control structures. For the cylindrical NW and NP geometries and for the square NP geometry, we used a "sequential" optimization scheme, which involved varying only one parameter at a time, updating that parameter with the new optimum value before proceeding to the next parameter. This sequential approach enabled us to look at a finer grid of values for each parameter, but it is important to note that better geometries might have been excluded, as changing one parameter can affect optimum values of another previously determined. One improved approach would be to perform a more self consistent approach where we cycle through sequentially until convergence of parameters. Another approach would be to use the sequential method to identify potential optimum values for each parameter and then consider all combinations of that courser set of parameters (as was done in the case of imbedded nanoparticles in this Section 2.3 for the NW with imbedded nanoparticles.) It is likely an order of magnitude improvement over the control would be seen for the the NW and square NP geometries. In regards to fabrication, it is important to note that the novel square NP geometry presented herein should both be easily be fabricated using soft-nanolithograpy [9].
By far the most dramatic increases in predicted SERS responses arise from the addition of 300 nm gold nanoparticles to the NW geometry, with optimization factors ranging from 10× to 2000×, with an average of 600× increase in the predicted SERS response. In regards to fabrication, particles up to 300 nm can be synthesized with a uniform spherical shape and narrow size distribution using a combination of chemical reduction and annealing [30]. The challenge, though, is an even distribution of particles within the wells and not on the top metal film. We believe it is possible for the following reasons. As seen in Table 2, the optimum NW geometry with imbedded particles consisted of a 700 nm periodicity and a 650 nm well. The wells, then, correspond to almost 70% of the total surface area. By controlling the concentration of 300 nm particles, configurations with a single particle in a well should be statistically favored.
Finally, while the structures herein were optimized for the benzene thiol band at 1073 cm −1 , it is important to note that plasmonic crystals can be designed and optimized using FDTD for any particular fingerprint Raman band for any analyte of interest. The significance of this work is in showing the potential of using FDTD to design and optimize novel plasmonic crystals in lieu of a costly trial-and-error fabrication approach.

Materials and Methods
To model the interactions of light with the nanowell/nanopost plasmonic crystals, we used Meep [31], a full three-dimensional finite-difference time-domain (FDTD) solver for Maxwell's curl equations. In the FDTD method, electric and magnetic fields are represented on a discrete, staggered grid and are propagated in time in using a leap-frog algorithm [15]. The material composition at each grid point is specified by a dielectric constant consistent with the spatial layout of each plasmonic crystal. Dielectric supports were modeled with a constant (non-dispersive) refractive index in line with photopolymers used in fabricated arrays: SU-8 (n = 1.59) [13] for NW geometries and NOA (n = 1.56) [14] for NP geometries. Gold was modeled using a dispersive Drude plus two Lorentzian fit of Johnson and Christy's dielectric data [32] over 300-1000 nm: with ∞ = 4.05665, ω D = 8.75866 eV, γ D = 0.04062 eV, ω 1 = 2.92893 eV, γ 1 = 0.83423 eV, σ 1 = 1.09960, ω 2 = 4.05286 eV, γ 2 = 1.74845 eV, and σ 3 = 2.2195. Periodic boundary conditions in the xand y-directions were used to simulate square arrays. Meep's perfectly matched layers (PML) were used as absorbing boundary conditions at the top and bottom of the computational cell to prevent spurious reflections [33]. The computational box size in the z-direction was 1500 nm with a grid spacing of 2 nm and 60 nm of PML in both directions. Light was propagated in the z-direction from air from a line source at the air/PML interface. With θ src defined such that θ src = 0 o corresponds to an x-polarized light source and θ src = 90 o corresponds to a y-polarized light source, θ src = 0 o was used for all NW and cylindrical NP geometries, and θ src = 45 o was used for the square NP geometries. Electric fields for λ = 785, 821, and 857 nm were propagated for 100 fs to ensure complete coupling into surface plasmons. Electric field intensities, with and without plasmonic crystals present, were averaged over the last two periods (c/λ) to get the time-averaged electric field intensity enhancement, g 2 =|E| 2 /|E 0 | 2 , at each wavelength. We calculated G 4 SERS and G 4 mid (Equations (1) and (2)) by averaging over g 2 only for points in air within 100 nm above the NW or NP. Points in the gold film or dielectric support were excluded from the average since they would not be accessible to a target molecule, and the 100 nm range from the NW or NP allowed for electric fields corresponding to surface plasmons to sufficiently decay. Zero-order transmission (T) spectra were calculated by taking the ratio of the transmitted power to the incident power integrated over a plane within the dielectric support. The extinction cross section for the 300 nm gold sphere was calculated using Mie theory [5,34]. Plots of electric field intensities were generated using ParaView [35,36].