1. Introduction

Many applications in the field of plasmonics are related to localized surface plasmon resonances of metallic nanoparticles [1]. The resonances induce strong local fields around the nanoparticles, which are useful in sensing, Raman scattering, near-field imaging, photovoltaics, nonlinear optics, etc. Nanoparticle pairs are of significant importance, because large field enhancements can occur in the narrow gaps between the nanoparticles [2]. Plasmonic resonances depend strongly on the size and shape of the nanoparticles, and for nanoparticle pairs the distance between them. Bowtie structures with the tips of two nanoparticles facing each other are particularly interesting [3–8], since sharp tips can further increase the local field enhancement. For practical applications, it is important to understand how the scattering spectra, resonant frequencies and field enhancements depend on geometric and material parameters of nanoparticle pairs, such as the apex angle and the tip-to-tip distance of a bowtie structure. All existing numerical studies on bowtie structures assume the tips are rounded with a small radius of curvature, probably because bowtie structures with mathematically sharp tips are not practical. However, when the gap between the tips is small, the field between the tips is very sensitive to the shape of the tips. If the tips are assumed to be locally circular with a small radius of curvature, then the field changes significantly as the radius is varied [8]. It is thus important to analyze idealized bowtie structures with truly sharp tips. The theoretical results for idealized structures can serve as useful guidelines when bowtie structures are designed for practical applications.

General numerical methods, such as the finite-difference time-domain (FDTD) method and the frequency-domain finite element method (FEM), can certainly be used to analyze metallic nanoparticles, but they require extremely small mesh sizes when the nanoparticles have sharp edges and corners. On the other hand, due to existing nano-fabrication techniques, the metallic nanoparticles are often tiny cylindrical objects, with their top and bottom surfaces parallel to the substrate and their side boundaries perpendicular to the substrate. The vertical mode expansion method (VMEM) is a special numerical method for cylindrical objects embedded in a layered background [9–13]. It is based on the observation that if the physical space can be divided into a number of cylindrical regions where the material properties depend only on a “vertical” variable z, then the electromagnetic field can be expanded in each region in one-dimensional (1D) eigenmodes with unknown “coefficients” that are two-dimensional (2D) functions. The 1D eigenmodes are functions of z and the 2D unknown functions satisfy scalar Helmholtz equations. The method was originally developed for a slab structure with circular holes [9]. In [12,13], we extended the VMEM to cylindrical structures with arbitrary cross sections, using a boundary integral equation (BIE) method to handle the 2D unknown functions. The BIEs are formulated on curves corresponding to the boundaries of the 2D cross sections of the cylinders.

The VMEM presented in [12,13] is suitable for cylindrical structures with smooth boundaries. Without a proper modification, the method cannot be used to analyze cylindrical bowtie structures with mathematically sharp tips, because the tips correspond to corners in the cross sections. Fortunately, there are a number of well-established methods for BIEs on 2D domains with corners [14–16]. In this paper, we present an improved VMEM suitable for cylindrical structures with sharp edges, based on a particular BIE method for 2D domains with corners [16]. The new VMEM involves a scaling for the unknown functions, a graded mesh technique for discretizing the boundaries, and a quadrature scheme for discretizing the boundary integral operators. As a demonstration of the application of the new method, we present computational results that provide accurate predictions on the effect of some key parameters on field enhancements and scattering efficiencies for a silver bowtie structure. Our findings are compared with known observations. They further demonstrate that sharp tips help to increase the field intensity in between tips.

2. Vertical mode expansion method

A bowtie structure is shown in Fig. 1(a). It consists of two closely spaced metallic cylinders on a dielectric substrate with their sharp tips facing each other. The top surface of the substrate is assumed to be the xy plane. The dielectric constants of the substrate, the metallic cylinders and the medium above the structure are ε_b, ε_m and ε_t, respectively. The top view of the structure is shown in Fig. 1(b), where Ω₁ and Ω₂ are 2D domains corresponding to the cross sections of the cylinders. We let Γ₁ and Γ₂ be the boundaries of Ω₁ and Ω₂, respectively, and denote the unbounded domain outside Ω₁ and Ω₂ by Ω₀, and the boundary of Ω₀ by Γ₀ (thus, Γ₀ = Γ₁ ∪ Γ₂). Corresponding to Ω₁, Ω₂ and Ω₀, we have three-dimensional (3D) cylindrical regions S₁, S₂ and S₀, given by S_l = {(x, y, z) | (x, y) ∈ Ω_l, −∞ < z < ∞}. In S_l, the dielectric function ε depends only on z, and is denoted as ε^(l)(z). Let D be the height of the cylinders, then ε^(l)(z) = ε_t for z > D, ε^(l)(z) = ε_b for z < 0, ε⁽⁰⁾(z) = ε_t and ε⁽¹⁾(z) = ε⁽²⁾(z) = ε_m for 0 < z < D.

 

Fig. 1 (a): A bowtie structure of a pair of metallic cylinders placed on a substrate. (b): The cross section of the bowtie structure in the xy plane.

Download Full Size | PDF

To analyze the scattering and local field enhancement of electromagnetic waves by the bowtie structure, we need to solve the frequency-domain Maxwell’s equations for electric field E and magnetic field H. In the top region, we specify an incident plane wave {E⁽ⁱ⁾, H⁽ⁱ⁾} with a wave vector (α, β, −γ), where α and β are real, $\gamma ={\left(k_0^2{\epsilon }_t-{\alpha }^2-{\beta }^2\right)}^{1/2}$ is positive, and k₀ is the free space wavenumber. Early versions of the VMEM are only applicable to circular cylinders in a layered background medium [9–11]. In [12], we developed a more general VMEM for cylinders with arbitrary cross sections. The method was presented for a single cylinder with a smooth boundary, but the extension to multiple cylinders is straightforward. For two cylinders with smooth boundaries, the VMEM takes the following steps.

More details about the method are given in our previous publications [11,12]. Here, we give a few remarks to highlight the main features of the method. Step 1 is needed, since the total field containing the incident wave, is not outgoing as z → +∞, but the mode expansions used in the VMEM are only applicable to outgoing fields. This is related to Step 2 where z is truncated by PMLs. The PML technique is widely used to model outgoing wave fields [18]. Since the 1D eigenmodes are calculated with z truncated by PMLs [11,19], the expansions based on these 1D eigenmodes can only approximate outgoing wave fields. In the VMEM, the mode expansions are applied to the difference between the total field and the reference solutions. In Steps 3 and 4, the tangential differential operators and the NtD operators are approximated by matrices. They are needed in Step 5 to set up the linear system, because the mode expansions involve not only $V_j^{(l,p)}$ but also their tangential and normal derivatives. The linear system in Step 5 is obtained from the continuity of H_z, E_z, H_τ and E_τ on the vertical boundaries of the cylindrical regions, i.e., W_l = {(x, y, z) | (x, y) ∈ Γ_l, −∞ < z < ∞}. Here, H_τ and E_τ are the horizontal tangential components.

In the following, we discuss the changes needed when Ω₁ and Ω₂ have corners. These changes are mostly related to Steps 4 and 5. To simplify the presentation, we drop the subscript j and superscript (l, p) for all related quantities such as $V_j^{(l,p)}$, ${\eta }_j^{(l,p)}$, $v_j^{(l,p)}$ and $w_j^{(l,p)}$. The function V satisfies the Helmholtz equation ${\partial }_x^{2}V+{\partial }_y^{2}V+{\eta }^{2}V=0$ in Ω_l. The related free-space Green’s function is $G(\mathbf{r},\tilde{\mathbf{r}})=\left(\text{i}/4\right)H_0^{(1)}\left(\eta \left|\mathbf{r}-\tilde{\mathbf{r}}\right|\right)$, where r = (x, y), r̃ = (x̃, ỹ), and $H_0^{(1)}$ is the Hankel function of the first kind and zeroth order. For l = 1 or 2, we have a BIE (1 + 𝒦)v = 𝒮w, where 𝒮 and 𝒦 are boundary integral operators given by

If the domain Ω_l has a corner, w may blow up due to possible field singularities. To overcome this difficulty, we replace w by a new unknown function, use a graded mesh technique to discretize Γ_l and a proper quadrature rule to approximate the boundary integral operators. The graded mesh technique was originally developed by Kress [14]. It uses a special parametric representation r = r(t) for Γ_l, such that the first a few derivatives of r(t) vanish at the corner, and then discretize Γ_l based on a uniform sampling of t. As in [15], we introduce a new unknown function q = |r′(t)|w and an integral operator 𝒮̃ such that 𝒮̃q = 𝒮w. The BIE is then written as (ρ + 𝒦)v = 𝒮̃q, where ρ = ρ(r) is the inner angle of Ω_l at r ∈ Γ_l divided by π. The original NtD operator is replaced by the modified NtD operator 𝒩̃ satisfying 𝒩̃q = v and it is given by 𝒩̃ = (ρ + 𝒦)⁻¹𝒮̃. With the graded mesh, the integral operators 𝒦 and can be discretized using the same method for boundary integral operators on smooth curves. The method developed by Kress is highly accurate for Helmholtz equations with a real η [14]. For us, Alpert’s hybrid Gauss-trapezoidal rule is preferred, since η may have a large imaginary part. More details on the discretization of the boundary integral operators are given in [16].

After computing the modified NtD operators ${\widetilde{𝒩}}_j^{(l,p)}$ and the tangential differential operators 𝒯₁ and 𝒯₂, we can set up a linear system for all $q_j^{(l,p)}$. If Γ_l is discretized by M_l points, $q_j^{(l,p)}$ is approximated by a vector of length M_l. The linear system is obtained by enforcing the continuity of H_z, E_z, |r′(t)|H_τ and |r′(t)|E_τ on W_l (the vertical boundaries of S_l) for l = 1, 2. The process is similar to that given in [13].

3. Numerical results

A bowtie structure consisting of two identical drop-shaped cylinders with their tips facing each other is shown in Fig. 1(a). A Cartesian coordinate system is chosen so that the bottom of the bowtie structure is in the horizontal plane z = 0, the bisectors of the structure are the x and y axes, and the x axis passes through the two tips. The drop-shaped cylinders are assumed to have a tip-to-base length L, a thickness D, and an apex angle θ. The distance between the two tips is denoted as G. The cross-section boundary Γ₁ of the right cylinder is given by

First, we compare a bowtie structure embedded in air with one placed on a dielectric substrate with ε_b = 2.1609. The apex angle and the gap between the tips are θ = π/3 and G = 5 nm, respectively. The calculated field enhancements and scattering efficiencies are shown in Fig. 2(a) and 2(b), respectively. The blue and red curves correspond to bowtie structures in air and on a substrate, respectively. These two cases show significant differences in both peak wavelengths and peak magnitudes. The peak wavelengths of the blue and red curves in both figures are about 756 nm and 882 nm, respectively. The red curves also have smaller peak magnitudes. Therefore, the dielectric substrate causes the resonant wavelength to red-shift and the near-field intensity to decrease. The same conclusion was obtained earlier by Jiao et al. [6]. Our results are also consistent with previous FEM results for a bowtie structure with rounded corners in air [8]. For our bowtie structure in air, there are three peaks at 756 nm, 513 nm and 441 nm, and the field enhancement at 756 nm is 865. Rosen and Tao [8] studied a silver bowtie structure with the same thickness D = 8 nm, the same apex angle θ = π/3 and the same gap G = 5 nm, but with rounded corners of radius 1 nm, and they obtained a peak at 765 nm for a field enhancement of 479, and two smaller peaks at 526 nm and 462 nm. Therefore, the field enhancement is much lower if the tips are rounded. In our calculations, the z variable is truncated to (−100, 108) nm, and discretized by N = 78 or 108 points for wavelength below or above 500 nm, respectively. The boundary Γ₁ is discretized by M 1 = 50 points. Reflection symmetries of the bowtie structure are utilized to reduce the size of the final linear system by a factor of 4.

 

Fig. 2 Field enhancements and scattering efficiencies for bowtie structures: (a,b) with or without a substarte; (c,d) with different tip-to-tip gaps; (e,f) with different apex angles.

Download Full Size | PDF

Next, we consider the effect of varying gap size to the silver bowtie structure (with apex angle θ = π/3) on the substrate (ε_b = 2.1609). Our numerical results are shown in Fig. 2(c) and 2(d) for G = 2.5, 5, 10, 20, 40 nm. It is clear that as the gap is decreased, the local field enhancement increases rapidly, the resonant wavelength increases slightly, and the peak scattering efficiency also increases slightly. These results are consistent with previous experimental and numerical results reported in [3] and [4]. Finally, we study the effect of varying apex angle to the bowtie structure on the substrate with the fixed gap G = 5 nm. Our results for θ = 60°, 70°, 80° and 90° are shown in Fig. 2(e) and 2(f). We observe that as the apex angle is decreased, the resonant wavelength increases, the field enhancement and scattering efficiency also increase. Our results agree with the experimental and numerical results of [7].

4. Conclusion

Bowtie structures of metallic nanoparticles are extremely useful in applications that require strong local fields. The field enhancement of a bowtie structure depend not only on the distance between the two tips but also on the shape of the tips. Existing numerical studies have all assumed the tips are rounded with a small radius of curvature. We study idealized bowtie structures with mathematically sharp tips, based on a new version of the VMEM incorporating techniques for BIEs on 2D domains with corners. Our numerical results indicate that sharp tips help to increase the field intensity in between the two tips. The VMEM is a special method for cylindrical structures in a layered background, and it can resolve the singular electromagnetic fields at the two tips without using excessive computing resources.