Analyzing bowtie structures with sharp tips by a vertical mode expansion method

Bowtie structures of metallic nanoparticles are very effective in producing strong local fields needed in many applications. Existing numerical studies on bowtie structures are limited to those with rounded tips. Due to the field singularities at sharp edges and corners, accurate numerical solutions for bowtie structures with mathematically sharp tips are difficult to obtain. Based on an improved vertical mode expansion method (VMEM) that incorporates boundary integral equation techniques for domains with corners, we analyze bowtie structures with truly sharp tips. Numerical results are presented to reveal the effects of a few key factors including the distance between the tips, the apex angle and the substrate. © 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


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][4][5][6][7][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][10][11][12][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][15][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.

Vertical mode expansion method
A bowtie structure is shown in Fig. 1(a). It consists of two closely spaced metallic cylinders on a (a) dielectric substrate with their sharp tips facing each other. The top surface of the substrate is assumed to be the x y 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 where Ω 1 and Ω 2 are 2D domains corresponding to the cross sections of the cylinders. We let Γ 1 and Γ 2 be the boundaries of Ω 1 and Ω 2 , respectively, and denote the unbounded domain outside Ω 1 and Ω 2 by Ω 0 , and the boundary of Ω 0 by Γ 0 (thus, Γ 0 = Γ 1 ∪ Γ 2 ). Corresponding to Ω 1 , Ω 2 and Ω 0 , we have three-dimensional (3D) cylindrical regions S 1 , S 2 and S 0 , 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 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 (i) , H (i) } with a wave vector (α, β, −γ), where α and β are real, γ = (k 2 0 ε t − α 2 − β 2 ) 1/2 is positive, and k 0 is the free space wavenumber. Early versions of the VMEM are only applicable to circular cylinders in a layered background medium [9][10][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.
3. Discretize Γ 1 and Γ 2 by M 1 and M 2 points, respectively, and calculate M 1 × M 1 and M 2 × M 2 matrices that approximate the tangential differential operators T 1 and T 2 along Γ 1 and Γ 2 , respectively. The boundary Γ 0 of the exterior domain Ω 0 is discretized by 5. Set up and solve a linear system for all w (l, p) j , l ∈ {0, 1, 2}, 1 ≤ j ≤ N and p ∈ {e, h}.
6. Evaluate the electromagnetic field at desired locations and calculate quantities such as the scattered power.
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 (l, p) j 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 Ω 1 and Ω 2 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

∂G(r,r) ∂ν(r) ψ(r)ds(r),
and ν(r) is the outward unit normal vector of Γ l at r ∈ Γ l [12,14]. The NtD operator is given by N = (1 + K) −1 S. For the exterior domain Ω 0 , the BIE is slightly different [12]. 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 operatorS such thatSq = Sw. The BIE is then written as (ρ + K)v =Sq, 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 bỹ N = (ρ + K) −1S . With the graded mesh, the integral operators K andS 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Ñ (l, p) j and the tangential differential operators T 1 and T 2 , we can set up a linear system for all q (l, p) j . If Γ l is discretized by M l points, q (l, p) j 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].

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 Γ 1 of the right cylinder is given by r = [0.5G + L sin(0.5s), 0.5L tan(0.5θ) sin(s)], 0 ≤ s < 2π (1) and the parameter s is further transformed to t (for the graded mesh) as in [14,16]. The shape of this bowtie structure is chosen for its convenience in our numerical implementation. The local field around the tips and the resonant frequency are insensitive to the shape of the structure away from the tips. Our results are relevant to all bowtie structures with the same shape around the tips. In the top region, a normal incident plane wave with its electric field in the x direction is given. Using the VMEM, we calculate scattering efficiencies and local field enhancements. The scattering efficiency is defined as the scattering cross-section divided by the geometric cross-section of the bowtie structure. The field enhancement is defined as |E|/|E (i) |, where the total electric field E is evaluated at point (0, 0, D/2) in the middle of the bowtie structure. For discretizing the boundary integral operators, we use Alpert's second order hybrid Gauss-trapezoidal rule for logarithmic singularities [16,17]. Our numerical results are given for a silver bowtie structure with L = 100 nm and D = 8 nm. The refractive indices of silver at different wavelengths are interpolated from the results of Johnson and Christy [20].
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 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 Γ 1 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.
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].

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.