## Abstract

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

## 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}*,

*ε*and

_{m}*ε*, respectively. The top view of the structure is shown in Fig. 1(b), where Ω

_{t}_{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*, the dielectric function

_{l}*ε*depends only on

*z*, and is denoted as

*ε*

^{(l)}(

*z*). Let

*D*be the height of the cylinders, then

*ε*

^{(l)}(

*z*) =

*ε*for

_{t}*z*>

*D*,

*ε*

^{(l)}(

*z*) =

*ε*for

_{b}*z*< 0,

*ε*

^{(0)}(

*z*) =

*ε*and

_{t}*ε*

^{(1)}(

*z*) =

*ε*

^{(2)}(

*z*) =

*ε*for 0 <

_{m}*z*<

*D*.

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, $\gamma ={\left({k}_{0}^{2}{\epsilon}_{t}-{\alpha}^{2}-{\beta}^{2}\right)}^{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–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.

- For each
*l*∈ {0, 1, 2} and the given incident wave, find the reference solution {**E**^{(l)},**H**^{(l)}} in the layered medium with*ε*=*ε*^{(l)}(*z*) for all (*x*,*y*). - Truncate
*z*by perfectly matched layers (PMLs) [18, 19], discretize*z*by*N*points, and for each*l*∈ {0, 1, 2}, calculate 1D eigenmodes $\left\{{\varphi}_{j}^{(l,p)}(z),{\eta}_{j}^{(l,p)}\right\}$ for*j*= 1, 2, ...,*N*, and*p*∈ {*e*,*h*}. Here, ${\varphi}_{j}^{(l,p)}(z)$ is the mode profile and ${\eta}_{j}^{(l,p)}$ is the propagation constant. The transverse electric (TE) and transverse magnetic (TM) modes are denoted by*p*=*e*and*p*=*h*, respectively. - 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*𝒯*_{1}and*𝒯*_{2}along Γ_{1}and Γ_{2}, respectively. The boundary Γ_{0}of the exterior domain Ω_{0}is discretized by*M*_{0}=*M*_{1}+*M*_{2}points. - For each
*l*∈ {0, 1, 2}, each*j*∈ {1, 2, ...,*N*}, each*p*∈ {*e*,*h*}, and function ${V}_{j}^{(l,p)}(x,y)$ satisfying a Helmholtz equation in Ω, use a BIE to find an_{l}*M*×_{l}*M*matrix that approximates the Neumann-to-Dirichlet (NtD) operator ${\mathcal{N}}_{j}^{(l,p)}$. The NtD operator maps the normal derivative of ${V}_{j}^{(l,p)}$ on Γ_{l}(denoted as ${w}_{j}^{(l,p)}$) to ${V}_{j}^{(l,p)}$ on Γ_{l}(denoted as ${v}_{j}^{(l,p)}$)._{l} - Set up and solve a linear system for all ${w}_{j}^{(l,p)}$,
*l*∈ {0, 1, 2}, 1 ≤*j*≤*N*and*p*∈ {*e*,*h*}. - 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}_{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 Ω_{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}_{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

*ν*(

**r**) is the outward unit normal vector of Γ

*at*

_{l}**r**∈ Γ

*[12,14]. The NtD operator is given by*

_{l}*𝒩*= (1 +

*𝒦*)

^{−1}

*𝒮*. 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 Γ

*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*

_{l}**r**=

**r**(

*t*) for Γ

*, such that the first a few derivatives of*

_{l}**r**(

*t*) vanish at the corner, and then discretize Γ

*based on a uniform sampling of*

_{l}*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 Ω

*at*

_{l}**r**∈ Γ

*divided by*

_{l}*π*. The original NtD operator is replaced by the modified NtD operator

*𝒩̃*satisfying

*𝒩̃q*=

*v*and it is given by

*𝒩̃*= (

*ρ*+

*𝒦*)

^{−1}

*𝒮̃*. 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 ${\tilde{\mathcal{N}}}_{j}^{(l,p)}$ and the tangential differential operators *𝒯*_{1} and *𝒯*_{2}, we can set up a linear system for all ${q}_{j}^{(l,p)}$. If Γ* _{l}* is discretized by

*M*points, ${q}_{j}^{(l,p)}$ is approximated by a vector of length

_{l}*M*. The linear system is obtained by enforcing the continuity of

_{l}*H*,

_{z}*E*, |

_{z}**r′**(

*t*)|

*H*and |

_{τ}**r′**(

*t*)|

*E*on

_{τ}*W*(the vertical boundaries of

_{l}*S*) for

_{l}*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 Γ_{1} of the right cylinder is given by

*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(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 Γ

_{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].

## 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.

## Funding

China Postdoctoral Science Foundation (2016M600903); Research Grants Council of Hong Kong Special Administrative Region, China (CityU 11304117); US NSF DMS-1719699.

## Acknowledgment

Part of this work was carried out by the first author at the Beijing Computational Science Research Center.

## References

**1. **S. A. Maier, *Plasmonics: Fundamentials and Applications* (Springer, 2007).

**2. **P. Nordlander, C. Oubre, E. Prodan, K. Li, and M. I. Stockman, “Plasmon hybridization in nanoparticle dimers,” Nano Lett. **4**, 899–903 (2004). [CrossRef]

**3. **D. P. Fromm, A. Sundaramurthy, P. J. Schuck, G. Kino, and W. E. Moerner, “Gap-dependent optical coupling of single ‘bowtie’ nanoantennas resonant in the visible,” Nano Lett. **4**, 957–961 (2004). [CrossRef]

**4. **A. Sundaramurthy, K. B. Crozier, G. S. Kino, D. P. Fromm, P. J. Schuck, and W. E. Moerner, “Field enhancement and gap-dependent resonance in a system of two opposing tip-to-tip Au nanotriangles,” Phys. Rev. B **72**, 165409 (2005). [CrossRef]

**5. **A. Kinkhabwala, Z. Yu, S. Fan, Y. Avlasevich, K. Müllen, and W. E. Moerner, “Large single-molecule fluorescence enhancements produced by a bowtie nanoantenna,” Nat. Photonics **3**, 654–657 (2009). [CrossRef]

**6. **X. Jiao, J. Goeckeritz, S. Blair, and M. Oldham, “Localization of near-field resonances in bowtie antennae: influence of adhesion layers,” Plasmonics **4**, 37–50 (2009). [CrossRef]

**7. **W. Ding, R. Bachelot, S. Kostcheev, P. Royer, and R. Lamaestre, “Surface plasmon resonances in silver bowtie nanoantennas with varied bow angles,” J. Appl. Phys. **108**, 124314 (2010). [CrossRef]

**8. **D A. Rosen and A. R. Tao, “Modeling the optical properties of bowtie antenna generated by self-assembled Ag triangular nanoprisms,” ACS Appl. Mater. Interfaces **6**, 4134–4142 (2014). [CrossRef] [PubMed]

**9. **S. Boscolo and M. Midrio, “Three-dimensional multiple-scattering technique for the analysis of photonic-crystal slabs,” J. Lightwave Technol. **22**, 2778–2786 (2004). [CrossRef]

**10. **L. Yuan and Y. Y. Lu, “Dirichlet-to-Neumann map method for analyzing hole arrays in a slab,” J. Opt. Soc. Am. B **27**(12), 2568–2579 (2010). [CrossRef]

**11. **X. Lu, H. Shi, and Y. Y. Lu, “Vertical mode expansion method for transmission of light through a single circular hole in a slab,” J. Opt. Soc. Am. A **31**, 293–300 (2014). [CrossRef]

**12. **H. Shi and Y. Y. Lu, “Efficient vertical mode expansion method for scattering by arbitrary layered cylindrical structures,” Opt. Express **23**, 14618–14629 (2015). [CrossRef] [PubMed]

**13. **H. Shi, X. Lu, and Y. Y. Lu, “Vertical mode expansion method for numerical modeling of biperiodic structures,” J. Opt. Soc. Am. A **33**, 836–844 (2016). [CrossRef]

**14. **R. Kress, “A Nyström method for boundary integral equations in domains with corners,” Numer. Math. **58**, 145–161 (1990). [CrossRef]

**15. **W. Lu and Y. Y. Lu, “Waveguide mode solver based on Neumann-to-Dirichlet operators and boundary integral equations,” J. Comput. Phys. **231**, 1360–1371 (2012). [CrossRef]

**16. **H. Shi, Y. Y. Lu, and Q. Du, “Calculating corner singularities by boundary integral equations,” J. Opt. Soc. Am. A **34**, 961–966 (2017). [CrossRef]

**17. **B. K. Alpert, “Hybrid Gauss-trapezoidal quadrature rules,” SIAM J. Sci. Comput. **20**, 1551–1584 (1999). [CrossRef]

**18. **J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. **114**, 185–200 (1994). [CrossRef]

**19. **H. Derudder, D. De Zutter, and F. Olyslager, “Analysis of waveguide discontinuities using perfectly matched layers,” Electron. Lett. **34**, 2138–2140 (1998). [CrossRef]

**20. **P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” Phys. Rev. B **6**, 4370–4379 (1972). [CrossRef]