Generation of surface plasmon-polaritons by edge effects

TL;DR

This paper investigates how edges of atomically thin metamaterial sheets generate surface plasmon-polaritons, using numerical and analytical methods to understand their behavior and influence in nanophotonics.

Contribution

It introduces a combined numerical and analytical approach to model edge-induced SPPs on 2D metamaterial sheets, highlighting the impact of geometry and conductivity changes.

Findings

Edges act as sources of dominant SPPs.

Sharp geometries significantly affect plasmon behavior.

Adaptive mesh refinement improves simulation accuracy.

Abstract

By using numerical and analytical methods, we describe the generation of fine-scale lateral electromagnetic waves, called surface plasmon-polaritons (SPPs), on atomically thick, metamaterial conducting sheets in two spatial dimensions (2D). Our computations capture the two-scale character of the total field and reveal how each edge of the sheet acts as a source of an SPP that may dominate the diffracted field. We use the finite element method to numerically implement a variational formulation for a weak discontinuity of the tangential magnetic field across a hypersurface. An adaptive, local mesh refinement strategy based on a posteriori error estimators is applied to resolve the pronounced two-scale character of wave propagation and radiation over the metamaterial sheet. We demonstrate by numerical examples how a singular geometry, e.g., sheets with sharp edges, and sharp spatial changes in the associated surface conductivity may significantly influence surface plasmons in nanophotonics.

Equations

$$\displaystyle(\boldsymbol{\mathcal{E}}(\boldsymbol{x},t),\boldsymbol{\mathcal{B}}(\boldsymbol{x},t))=\text{Re}\,\big{\{}e^{-i\omega t}(\boldsymbol{E}(\boldsymbol{x}),\boldsymbol{B}(\boldsymbol{x}))\big{\}}\leavevmode\nobreak\ .$$

$$\displaystyle\begin{cases}\begin{aligned} -i\omega\boldsymbol{B}^{\text{sc}}+\nabla\times\boldsymbol{E}^{\text{sc}}\;&=\;0,\\[1.00006pt] i\omega\tilde{\varepsilon}\boldsymbol{E}^{\text{sc}}+\nabla\times\big{(}\mu^{-1}\boldsymbol{B}^{\text{sc}}\big{)}\;&=\;0,\\[1.00006pt] \nabla\cdot\boldsymbol{B}^{\text{sc}}\;&=\;0,\\[1.00006pt] \nabla\cdot\big{(}\tilde{\varepsilon}\boldsymbol{E}^{\text{sc}}\big{)}\;&=\;0.\end{aligned}\end{cases}$$

$$\displaystyle\begin{cases}\begin{aligned} \{\boldsymbol{\nu}\times\big{[}\big{(}\mu^{-1}\boldsymbol{B}\big{)}^{+}-\big{(}\mu^{-1}\boldsymbol{B}\big{)}^{-}\big{]}\}\Big{|}_{\Sigma}\;&=\;\{\sigma^{\Sigma}(\boldsymbol{x})\,[(\boldsymbol{\nu}\times\boldsymbol{E})\times\boldsymbol{\nu}]\}\Big{|}_{\Sigma}\leavevmode\nobreak\ ,\\[1.99997pt] \{\boldsymbol{\nu}\times\big{(}\boldsymbol{E}^{+}-\boldsymbol{E}^{-}\big{)}\}\Big{|}_{\Sigma}\;&=\;0\leavevmode\nobreak\ ,\end{aligned}\end{cases}$$

$$\underset{|\mathbf{x}|\to \infty }{\lim }\left\{{\mathbf{B}}^{\text{sc}}\times \mathbf{x}-c^{-1}|\mathbf{x}|{\mathbf{E}}^{\text{sc}}\right\}=0\text{leavevmode},\underset{|\mathbf{x}|\to \infty }{\lim }\left\{{\mathbf{E}}^{\text{sc}}\times \mathbf{x}+c|\mathbf{x}|{\mathbf{B}}^{\text{sc}}\right\}=0(\mathbf{x}\notin \Sigma )\text{leavevmode}.$$

$$\left|\frac{\omega \mu {\sigma }^{\Sigma }}{k}\right|\ll 1\text{leavevmode},$$

$$u(x)=u^{\mathrm{in}}(x)+\frac{i\omega \mu {\sigma }^{\Sigma }}{k}\left(\frac{d^2}{d{x}^2}+1\right){\int }_0^{\infty }d{x}^{\prime }\mathfrak{K}(x-x^{\prime })u(x^{\prime })x>0\text{leavevmode},$$

$$u(x)-\frac{i\omega \mu {\sigma }^{\Sigma }}{k}\left(\frac{d^2}{d{x}^2}+1\right){\int }_{-\infty }^{\infty }d{x}^{\prime }\mathfrak{K}(x-x^{\prime })u(x^{\prime })=u^{\mathrm{in}}(x)+g(x)x\in \mathbb{R}\text{leavevmode}.$$

$$u(x)\approx$$

$$-e^{-{\hat{Q}}_{+}(\sin \alpha )+{\hat{Q}}_{+}(k_m/k)}{e}^{i(k_m/k)x}\text{leavevmode}.$$

$$\widehat{Q}_{+}(\xi)=\frac{\xi}{\pi i}\int_{0}^{\infty}\frac{d\zeta}{\zeta^{2}-\xi^{2}}\,\ln\Biggl{(}1+\frac{\omega\mu\sigma^{\Sigma}}{2k}\sqrt{1-\zeta^{2}}\Biggr{)}\leavevmode\nobreak\ ,\qquad{\rm Im}\xi>0\leavevmode\nobreak\ ,$$

$${\hat{Q}}_{+}(0)\approx 0.0005+0.05i,{\hat{Q}}_{+}(k_m/k)\approx 0.346574+0.392699i\text{leavevmode}.$$

$$u(x)\approx \frac{2}{\sqrt{\pi }}\sqrt{\frac{2k}{\omega \mu {\sigma }^{\Sigma }}}{e}^{-i\pi /4}\sqrt{x}\text{mbox}asx\searrow 0\text{leavevmode}.$$

$$\mathbf{x}\to k\mathbf{x},\nabla \to \frac{1}{k}\nabla ,k_{m,r}=\frac{k_m}{k},$$

$$\mu \to {\mu }_r=\frac{1}{{\mu }_0}\mu ,\tilde{\epsilon }\to {\tilde{\epsilon }}_r=\frac{1}{{\epsilon }_0}\tilde{\epsilon },{\sigma }^{\Sigma }\to \mathrm{Space}\text{leavevmode}Sci.\text{leavevmode}Rev.=\sqrt{\frac{{\mu }_0}{{\epsilon }_0}}{\sigma }^{\Sigma }.$$

$$\displaystyle\nabla\times\big{(}\mu_{r}^{-1}\,\nabla\times\boldsymbol{E}\big{)}-\tilde{\varepsilon}_{r}\,\boldsymbol{E}\;=\;0.$$

$$\begin{gathered} {\int }_{\Omega }({\mu }_r^{-1}\nabla \times {\mathbf{E}}^{\text{sc}})\cdot (\nabla \times \bar{\mathbf{\phi }})\mathrm{d}x-{\int }_{\Omega }({\tilde{\epsilon }}_r{\mathbf{E}}^{\text{sc}})\cdot \bar{\mathbf{\phi }}\mathrm{d}x \\ -i{\int }_{\Sigma }(\mathrm{Space}\text{leavevmode}Sci.\text{leavevmode}Rev.{\mathbf{E}}_T^{\text{sc}})\cdot {\bar{\mathbf{\phi }}}_T\mathrm{d}{o}_x-i{\int }_{\partial \Omega }\sqrt{{\mu }_r^{-1}{\tilde{\epsilon }}_r}{\mathbf{E}}_T^{\text{sc}}\cdot {\bar{\mathbf{\phi }}}_T\mathrm{d}{o}_x= \\ i{\int }_{\Sigma }(\mathrm{Space}\text{leavevmode}Sci.\text{leavevmode}Rev.{\mathbf{E}}_T^{\text{in}})\cdot {\bar{\mathbf{\phi }}}_T\mathrm{d}{o}_x\text{leavevmode}. \end{gathered}$$

$$\displaystyle\boldsymbol{X}(\Omega)=\Big{\{}\boldsymbol{\varphi}\in\boldsymbol{L}^{2}(\Omega)\;:\;\nabla\times\boldsymbol{\varphi}\in L^{2}(\Sigma)^{3},\;\boldsymbol{\varphi}_{T}\big{|}_{\Sigma}\in L^{2}(\Sigma)^{3},\;\boldsymbol{\varphi}_{T}\big{|}_{\partial\Omega}\in L^{2}(\partial\Omega)^{3}\Big{\}}\leavevmode\nobreak\ ,$$

$$\displaystyle\{\boldsymbol{\nu}\times\big{(}\boldsymbol{E}^{+}-\boldsymbol{E}^{-}\big{)}\}\Big{|}_{\Sigma}=0$$

$$\displaystyle\{\boldsymbol{\nu}\times\big{[}\big{(}\mu^{-1}\boldsymbol{B}\big{)}^{+}-\big{(}\mu^{-1}\boldsymbol{B}\big{)}^{-}\big{]}\}\Big{|}_{\Sigma}\;$$

$$\{\begin{array}{l}\begin{array}{rl}{\mu }_r^{-1}& ⟶{\breve{\mu }}_r^{-1}=B{\mu }_r^{-1}A\text{leavevmode},\\ {\tilde{\epsilon }}_r& ⟶{\breve{\epsilon }}_r=A^{-1}{\tilde{\epsilon }}_r{B}^{-1}\text{leavevmode},\\ \mathrm{Space}\text{leavevmode}Sci.\text{leavevmode}Rev.& ⟶{\breve{\sigma }}_r^{\Sigma }=C^{-1}\mathrm{Space}\text{leavevmode}Sci.\text{leavevmode}Rev.B^{-1}\text{leavevmode}.\end{array}\end{array}$$

$$\displaystyle A=T^{-1}_{\boldsymbol{e}_{x}\boldsymbol{e}_{r}}\text{diag}\,\Big{(}\frac{1}{\bar{d}^{2}},\frac{1}{d\bar{d}},\frac{1}{d\bar{d}}\Big{)}T_{\boldsymbol{e}_{x}\boldsymbol{e}_{r}}\leavevmode\nobreak\ ,\quad B=T^{-1}_{\boldsymbol{e}_{x}\boldsymbol{e}_{r}}\text{diag}\,\big{(}d,\bar{d},\bar{d}\big{)}T_{\boldsymbol{e}_{x}\boldsymbol{e}_{r}}\leavevmode\nobreak\ ,$$

$$\displaystyle C=T^{-1}_{\boldsymbol{e}_{x}\boldsymbol{e}_{r}}\text{diag}\,\Big{(}\frac{1}{\bar{d}},\frac{1}{\bar{d}},\frac{1}{d}\Big{)}T_{\boldsymbol{e}_{x}\boldsymbol{e}_{r}}\leavevmode\nobreak\ ,$$

$$d=1+is(r)\text{leavevmode},\bar{d}=1+i/r{\int }_{\rho }^{r}s(\tau )\text{d}\tau \text{leavevmode},$$

$${\eta }_K^2={\rho }_K^2{\omega }_K^2+{\rho }_{{\partial }_K}^2{\omega }_{\partial K}^2\text{leavevmode}.$$

$${\rho }_K^2$$

$${\rho }_{\partial K}^2$$

$$\displaystyle\omega_{K}^{2}=\int_{K}\big{\|}z-I_{h}z\big{\|}^{2}\,\text{d}x,\quad\omega_{\partial K}^{2}=\int_{K}\big{\|}z-I_{h}z\big{\|}^{2}\,\text{d}o_{x}\leavevmode\nobreak\ .$$

$$\displaystyle\mathcal{J}(\boldsymbol{E}):=\frac{1}{2}\int_{\Sigma}\big{\|}\nabla\times\boldsymbol{E}\,\big{\|}^{2}\rm{d}o_{x}$$

$$\begin{gathered} {\int }_{\Omega }({\mu }_r^{-1}\nabla \times \mathbf{\phi })\cdot (\nabla \times \bar{\mathbf{z}})\mathrm{d}x-{\int }_{\Omega }({\tilde{\epsilon }}_r\mathbf{\phi })\cdot \bar{\mathbf{z}}\mathrm{d}x-i{\int }_{\Sigma }(\mathrm{Space}\text{leavevmode}Sci.\text{leavevmode}Rev.{\mathbf{\phi }}_T)\cdot {\bar{\mathbf{z}}}_T\mathrm{d}{o}_x \\ -i{\int }_{\partial \Omega }\sqrt{{\mu }_r^{-1}{\tilde{\epsilon }}_r}{\mathbf{\phi }}_T\cdot {\bar{\mathbf{z}}}_T\mathrm{d}{o}_x={\int }_{\Sigma }(\nabla \times {\mathbf{E}}_h)\cdot (\nabla \times \bar{\mathbf{\phi }})\mathrm{d}x\text{leavevmode}. \end{gathered}$$

$$E_{x,r}^{\text{sc}}(x)\simeq E_0\left[\frac{1}{1+\frac{\mathrm{Space}\text{leavevmode}Sci.\text{leavevmode}Rev.}{2}}-1\right]-E_0{e}^{-{\hat{Q}}_{+}(\sin \alpha )+{\hat{Q}}_{+}(k_{m,r})}{e}^{i{k}_{m,r}x}\text{leavevmode},$$

$$E_{x,r}^{\mathrm{sc}}(x)\approx E_0\left[\frac{1}{1+\frac{\mathrm{Space}\text{leavevmode}Sci.\text{leavevmode}Rev.}{2}}-1\right]-A{E}_0{e}^{-{\hat{Q}}_{+}(\sin \alpha )+{\hat{Q}}_{+}(k_{m,r})}{e}^{i(k_{m,r}x+\phi )}\text{leavevmode},x>0\text{leavevmode},$$

$$\displaystyle\varphi(d)=\frac{A\,\pi}{\gamma}\,\bigl{[}1+(d-d_{0})^{2}/\gamma^{2}\bigr{]}^{-1}\leavevmode\nobreak\ .$$

Full text

{AMS}

65N30, 78M10, 78M30, 78A45

1 Introduction

Recently, advances in nanophotonics have been made possible through the design of atomically thick materials, e.g., graphene, with tunable, novel electronic structure and macroscopic conductivity [34, 18]. A far-reaching goal is to precisely control coherent light at the nanoscale, in the terahertz frequency range [40]. This objective can be pursued by manipulation of the microscopic parameters of conducting sheets. Under suitable conditions, the sheet may behave as a metamaterial, exhibiting a dielectric permittivity with a negative real part as a function of frequency. By appropriate current-carrying sources and geometry, one may generate electromagnetic waves that propagate with relatively short wavelength along the sheet [4]. A celebrated type of wave, linked to various technological applications, is the surface plasmon-polariton (SPP) [4, 16, 27].

The present paper focuses on the computation of time-harmonic SPPs intimately related to geometric effects, specifically the presence of edges on conducting films, in the frequency regime in which the metamaterial character of the sheet is evident. Our goal is to demonstrate numerically how spatial changes in the morphology and surface conductivity of the sheet may amplify, or spoil, the observation of SPPs. To this end, we combine computational tools that include: (i) an implementation of the finite element method for the discontinuity of the tangential magnetic field across the sheet which is inherent to a class of thin conducting metamaterials; and (ii) an approximate solution by the method of Wiener and Hopf for an integral equation describing the electric field tangent to a semi-infinite sheet (“reference case”) [29].

Our main results can be summarized as follows.

• –

We formulate a two-dimensional (2D) model for fine-scale SPPs induced by edge discontinuities. A few prototypical geometries are investigated numerically. In particular, we describe SPPs generated by: the edge discontinuity of a semi-infinite sheet; the gap between two co-planar, semi-infinite sheets; and a resonant, finite conducting strip.

• –

For the numerics, we adapt a variational framework for the finite element treatment of wave propagation [26] to the 2D setting of SPP generation by edge effects. The underlying weak formulation embodies a discontinuity of the tangential magnetic-field component across the sheet and point singularities at material edges or discontinuities of the surface conductivity.

• –

We validate our numerical treatment by comparison of numerics to an approximate solution for the reference case of the semi-infinite sheet. We verify the analytically predicted singular behavior of the electric field near the edge.

• –

Based on our numerical simulations, we predict means of enhancing the observation of SPPs. Specifically, we demonstrate that the presence of a highly conducting material in the gap between two co-planar, semi-infinite metamaterial sheets may lead to an increase of the SPP amplitude. We also show that in this setting the dependence of the SPP amplitude on the gap width is distinctly different from the case with an “empty gap”, i.e., when the material of the gap is the ambient medium. In addition, we characterize resonances for SPPs on a finite strip numerically, and demonstrate that the SPP maximal magnitude as a function of the gap width is described by a Lorentzian function.

It has been well known that, because of phase matching in electromagnetic wave propagation, SPPs cannot be excited by plane waves incident upon an (idealized) infinite conducting sheet [4]. The relatively large wavenumber of the desired SPP cannot be directly matched at the interface. This limitation is usually remedied by introduction of suitably tuned gratings or localized current-carrying sources such as dipole antennas near resonance [4, 1].

In this paper, we use a finite element approach to explore an alternate scenario: the generation of SPPs via diffraction of waves by material defects such as sharp edges of conducting sheets. We point out that such geometric singularities, edges, can induce SPPs under radiation by appropriately polarized plane waves. Our numerical simulations involve hypersurfaces embedded in a 2D space for the sake of computational ease. However, our key conclusions can be generalized to higher dimensions. One of these conclusions is that the induced SPP may dominate over the diffracted field in a range of distances away from the edge.

Our work here forms an extension of the approach in [26] to more realistic, singular geometries. In [26], we introduced a numerical framework for tackling propagation of SPPs on planar hypersurfaces via the finite element method. Our formulation relies on a variational statement that incorporates the weak discontinuity of the tangential component of the magnetic field across the sheet. This approach offers the advantage of local refinement based on a-posteriori error control. The computational work reported in [26] is restricted to idealized geometries, aiming to introduce a general platform rather than address realistic applications.

In the present paper, we take a major step closer to applications by extending our finite element approach [26] to settings with sharp edges and gaps on metamaterial, conducting films. Our ultimate purpose is to model and simulate the effect of defects on the generation and propagation of SPPs. We validate our numerics via comparison to an approximate analytical solution for the semi-infinite conducting sheet with recourse to the Wiener-Hopf method of factorization [29]. This solution analytically reveals that the algebraically singular field near the edge transcends to a slowly decaying, fine-scale SPP away from the edge.

1.1 Scope: Excitation of SPPs

The excitation of SPPs on interfaces between metals and dielectrics has been conceived as a means of confining and manipulating coherent light in the infrared spectrum [16]. Surface plasmons form the macroscopic manifestation of the resonant interaction of electromagnetic radiation (photons) with electrons in the plasma of the metal surface. SPPs, in particular, are loosely defined macroscopically as waves that have a relatively short wavelength and decay slowly along the interface. SPPs may offer significant technological advantages; for example, improvement of the emission of light from corrugated metallic surfaces or nanoparticles [19, 31].

The macroscopic theory of surface plasmonics usually relies on Maxwell’s equations, specifically notions of classical electromagnetic wave propagation near boundaries. This subject has been the focus of systematic investigations in the case with radiowaves propagating over the earth or sea; see, e.g., [23, 8]. A central concept is the surface or lateral wave, which is confined closely to the boundary [23]. Typically, at radio frequencies the surface wave between air and earth has a phase velocity approximately equal to the phase velocity in air.

In surface plasmonics in the terahertz frequency regime, however, the constitutive relations and geometry of the associated materials become more intricate. Hence, the character of the ensuing surface wave is substantially different from that at radio frequencies [27]. Here, we view the SPP as a type of surface wave. Notably, the wavelength of the SPP with transverse-magnetic polarization can plausibly be made much smaller than the wavelength of a plane wave in free space at the same frequency [4, 7].

An emerging question is the following. How can the geometry or surface conductivity of a low-dimensional material be controlled to generate a desirable SPP? A satisfactory answer to this question requires understanding through reliable computations how geometry affects solutions to a class of boundary value problems for Maxwell’s equations. In these problems, an atomically thick material, e.g., graphene, may introduce a jump proportional to a local property such as conductivity in components of the electromagnetic field tangent to the sheet. The coupling of this type of boundary condition to singular geometries in 2D is the subject of the present paper.

1.2 Past computational approaches

Computational methods for plasmonics appear to be tailored to specific applications. Next, we provide a brief, non-exhaustive summary of the main computational tools found in the existing literature.

First, we comment on analytical treatments. In relatively simple settings with planar boundaries, dispersion relations of SPPs have been sought analytically through the use of plane-wave excitations; see, e.g., the recent review in [4]. Interestingly, the corresponding geometries lack singular regions, e.g., sharp edges. A notable exception concerns a metallic contact modeled as a strip of fixed width on a substrate in 2D [4, 36]. For this geometry, an analytical solution has been found via solving an integral equation in the Fourier domain for the electric field tangent to the strip in the simplified case with an electrically small strip width [4]. This solution is not applicable to conducting sheets with edges, since in the latter case the edge may not be treated as a perturbation. A numerical study of electrically large metal contacts by spectral methods is provided in [36]; however, the transition of the field from a singular behavior near the edge to an SPP is not discussed in this work [36].

Other settings in plasmonics consist of dipole sources over infinite planar boundaries; then, analytical and semi-analytical solutions are developed via the Fourier transform of the field [20]. In some exceptional cases, when the dipole and observation point lie on the sheet, exact evaluation of field components is possible [28]. These approaches have the merit of yielding features of SPPs inherent to the nature of the point source; but they convey little or no information about how realistic geometric effects related to the finite size of the sheet may influence the SPP.

In more complicated geometries, the excitation of plasmons has been studied under the “quasistatic approximation” in which the typical size of the scatterer is much smaller than the wavelength of the incident radiation field [27]. However, this approximation is expected to be questionable as the frequency becomes higher or the size of the scatterer is comparable to the free-space wavelength [27].

A variety of numerical methods for nanophotonics have been reported in the literature; for a recent review, see, e.g., [17]. These approaches include: the volume integral-equation approach [21, 24, 22], which exploits the integral form of Maxwell’s equations along with the respective Green’s function in closed form in the frequency domain; and the akin boundary element method [14, 30, 37], which makes use of Dirac masses as test functions and often employs the scalar and vector potentials. Numerical methods of different nature are the finite-difference time-domain method [32, 25] and the finite-difference method [15], which invoke space discretization and the respective approximation of the electromagnetic field by piecewise-constant functions. This type of approximation may be challenged in plasmonics, where the electromagnetic field may vary appreciably over short distances. An improvement has been offered by the finite element method [33, 39, 8] with certain choices of basis functions, for example, divergence-free functions. An issue of importance is to satisfy the radiation condition far away from the source. This can be accomplished in the numerical scheme via the notion of the perfectly matched layer (PML) [5, 9]. In passing, we should also mention the discrete dipole approximation [13, 11] by which the (continuum) scatterer is replaced by a finite array of polarizable particles or dipoles; this technique can be viewed as the outcome of discretization of the volume integral equation.

1.3 Our computational treatment

Recently, we formulated a variational framework for the finite element treatment of wave propagation along metamaterial conducting sheets embedded in spaces of arbitrary dimensions [26]. Our formulation incorporates a weak discontinuity of the tangential component of the magnetic field; this jump is responsible for the fine scale of the SPP. The corresponding discretization scheme was implemented with a modern finite element toolkit [2], and accounts both for the small wavelength of the SPP as well as for the radiation condition at infinity via a PML. This approach allows for error control and, thus, an efficient numerical approximation of the underlying boundary value problem. In our past work [26], we validated this approach by numerical simulations restricted to settings with dipoles over infinite planar geometries in 2D.

With the present paper, our goal is to adapt the aforementioned approach to more realistic settings and real-world applications. Hence, we extend our finite element treatment to more complicated geometries, especially conducting sheets with edges and gaps in 2D. We demonstrate via simulations that our variational framework and numerical implementation by appropriate curl-conforming Nédélec-elements correctly produces the edge singularity of the electric field. Our numerics for finite sheets show how the diffracted electric field transcends from a singular behavior near each edge to the SPP away from the edge. We validate our numerical approach by use of a semi-infinite sheet: in this case, our numerical results are in excellent agreement with the analytical prediction [29]. This study places on a firm foundation our finite element approach for SPPs generated and sustained by defects and finite-size effects.

1.4 Pending issues

Our work here leaves a few open problems for near-future study. For example, we do not solve the related boundary value problem for Maxwell’s equations in three spatial dimensions (3D). This important topic is the subject of work in progress. Furthermore, in our model we use homogeneous and isotropic material parameters. The case with spatially varying conductivity of the sheet, where the SPP may be the result of homogenization [7], lies beyond our present scope. In a similar vein, we have not made any attempt to study more complicated forms of discontinuities for the electromagnetic-field components across the sheet, say, in the presence of a magneto-electric effect [38]. We should also mention the experimentally appealing case with a receiving antenna lying on graphene [1]. This problem can be described by our variational approach, and is the subject of future research as well.

Note that surface plasmonics comprise a class of multiscale problems: the electronic structure of low-dimensional materials should be linked to the phenomenology of SPPs. Here, we invoke a macroscopic model of SPPs via Maxwell’s equations. The emergence of the related phenomenology from microscopic principles needs to be further understood.

1.5 Paper organization

The remainder of this paper is organized as follows. In section 2, we introduce the relevant boundary value problem and summarize analytical results for the reference case [29]. Section 3 focuses on our numerical approach; in particular, the use of a PML and adaptive local refinement for good resolution of the SPP is discussed. In Section 4, we present and discuss computational results for the following geometries: a semi-infinite sheet (reference case); two co-planar, semi-infinite sheets with a gap; and a resonant, finite conducting strip. Section 5 concludes our paper with a summary of the main findings.

2 Model and analytical results

In this section, we introduce the main ingredients of the model and analytical approach: a boundary value problem for Maxwell’s equations, which incorporates a discontinuity for the magnetic-field component tangent to an arbitrary conducting sheet in 2D (section 2.1); and an analytical formula for the SPP on a semi-infinite sheet (section 2.2). Our formula manifests the fine scale of the SPP if the sheet conductivity satisfies a certain condition consistent with the metamaterial nature of the sheet. In this vein, we describe solutions for the electric field far from and near the edge on a semi-infinite sheet [29] (section 2.2).

2.1 Boundary value problem and geometry

Following [26], we formulate a boundary value problem for Maxwell’s equations by including a conducting sheet in an unbounded space; see Figure 1. The starting point is the strong form of Maxwell’s equations for the time-harmonic electromagnetic field $(\boldsymbol{\mathcal{E}}(\boldsymbol{x},t),\boldsymbol{\mathcal{B}}(\boldsymbol{x},t))$, viz.,

[TABLE]

By splitting the electric and magnetic fields into the primary (incident) field and the scattered field, $\boldsymbol{E}=\boldsymbol{E}^{\text{in}}+\boldsymbol{E}^{\text{sc}}$, $\boldsymbol{B}=\boldsymbol{B}^{\text{in}}+\boldsymbol{B}^{\text{sc}}$, Maxwell’s equations (in a source-free form) for the scattered field outside the sheet are written as [8]

[TABLE]

Here, Gauss’ law, expressed by the last two equations, is redundant. In (1), the material parameters are time independent. The second-rank tensors $\mu(\boldsymbol{x})$ and $\tilde{\varepsilon}(\boldsymbol{x})$ denote the effective magnetic permeability and complex permittivity of the unbounded medium; in particular, $\tilde{\varepsilon}(\boldsymbol{x})=\varepsilon(\boldsymbol{x})+i\sigma(\boldsymbol{x})/\omega$, where $\varepsilon(\boldsymbol{x})$ and $\sigma(\boldsymbol{x})$ are the usual permittivity and conductivity of the medium.

Equations (1) are enforced in suitable unbounded regions of the Euclidean space $\mathbb{R}^{n}$ ($n=2,\,3$), by exclusion of the conducting sheet. Hence, system (1) must be complemented with the appropriate boundary condition across the sheet which is represented by an oriented hypersurface $\Sigma$, $\Sigma\subset\mathbb{R}^{n}$, with unit normal $\boldsymbol{\nu}$ and surface conductivity $\sigma^{\Sigma}(\boldsymbol{x})$ [4]. For a conducting sheet, the electromagnetic field satisfies the conditions [4]

[TABLE]

where $\boldsymbol{E}^{\pm},\,\boldsymbol{B}^{\pm}$ is the restriction of the field to each side ($\pm$) of the sheet and $\sigma^{\Sigma}(\boldsymbol{x})$ is a second-rank tensor with domain on $\Sigma$.

In the presence of compactly supported external sources, e.g., electric or magnetic dipoles, (1) and (2) must be complemented with a radiation condition at infinity. For an isotropic medium, this condition takes the form of a Silver-Müller condition [8]. Specifically, if $c$ is the speed of light at infinity, we require that

[TABLE]

Motivated by the setting with an infinite planar sheet [4], we will pay particular attention to material properties and external sources that allow for the generation of a short-wavelength SPP on the conducting sheet. Specifically, for an isotropic and homogeneous sheet we expect that a necessary condition for such an SPP is $\text{Im}\,\sigma^{\Sigma}>0$. In addition, the sheet needs to be radiated by an appropriately (transverse-magnetic) polarized wave. We conclude this subsection with a practically appealing definition [7]:

Definition 2.1**.**

For scalar $\sigma^{\Sigma}$, the nonretarded frequency regime is characterized by

[TABLE]

*i.e., a surface resistivity ($1/\sigma^{\Sigma}$) that is much larger in magnitude than the intrinsic impedance of the ambient medium. *

By recourse to an analytical solution for the reference case (section 2.2), we will see that (4) along with the condition $\text{Im}\,\sigma^{\Sigma}>0$ imply that an SPP is present on a semi-infinite sheet and has a wavenumber, $k_{m}$, with $|k_{m}|\gg k$. This property in turn yields two distinct length scales of wave propagation in this problem.

2.2 Reference case: Explicit formulas

In the case with a semi-infinite sheet, $\Sigma=\{(x,0)\in\mathbb{R}^{2}\,:\,x\geq 0\}$ [see Figure 2(a)], Maxwell’s equations (1) with jump condition (2) admit a closed-form solution for the electric field, $E_{x}(x,0)$, tangential to the sheet. This solution is obtained by application of the Wiener-Hopf method to an integral equation [29]. In the following we briefly motivate and describe our formalism and results; for a detailed discussion we refer the reader to [29].

The key ingredient of the analytical approach is the formulation of (1) and (2) in terms of an integral equation for the (non-dimensional) $x$-component, $u(x)=E_{x}(k^{-1}x,0)/E_{0}$, of the electric field on the sheet, where $E_{0}$ is a typical amplitude of the incident electric field. The governing integral equation reads

[TABLE]

where $u^{\rm in}$ corresponds to the incident field, $x$ is scaled by $1/k$ where $k$ is real, and $\mathfrak{K}$ is an appropriate kernel, which corresponds to a Green’s function for the scalar Helmholtz equation; see [29, Sec. 1]. By formally extending $u$ and $u^{\rm in}$ to the whole real axis through setting $u(x)=u^{\rm in}(x)\equiv 0$ for $x<0$, and introducing an unknown function, $g$, for consistency ($g(x)\equiv 0$ if $x>0$), we write

[TABLE]

Equation (6) is solved explicitly in terms of a Fourier integral by use of the Wiener-Hopf method [29]. To this end, we consider a plane wave as the incident field, $E_{x}^{\rm in}=E_{0}e^{ix\,\sin\alpha}$, for some incidence angle, $\alpha$. In the following, we summarize the main findings of this approach.

Field far from the edge

In applications of plasmonics, the region consisting of points of the sheet sufficiently far from the edge deserves some attention. This region can be characterized by the condition $x\gg 1$. In this region, the SPP may dominate over the diffracted field; the latter is defined as the scattered field after removal of the direct reflection of the incident wave. The solution to (6) can be written as

[TABLE]

A few remarks on (7) are in order. The first term describes the sum of the incident field and its direct reflection from an infinite sheet, while the second term is the SPP. The effect of the edge is expressed by values of the split function $\widehat{Q}_{+}(\xi)$, which is defined by the contour integral [29]

[TABLE]

and can be evaluated approximately, in closed form, in the nonretarded regime (Definition 2.1) [29]. In principle, this $\widehat{Q}_{+}(\xi)$ enters the Fourier integral for the exact solution, $u$, for $x>0$ [29]. In the region under consideration (if $x\gg 1$), $\widehat{Q}(\xi)$ needs to be computed for particular values of $\xi$ in order to compare formula (7) against numerical results of the finite element method (section 4). For example, for $\alpha=0$, we compute $\widehat{Q}_{+}(\xi)$ for $\xi=0$ and $\xi=k_{m}/k$ approximately (see [29, Sec. 3]), viz.,

[TABLE]

We should add that the SPP wavenumber, $k_{m}$, obeys the dispersion relation $k_{\perp}:=\sqrt{k^{2}-k_{m}^{2}}=-2k^{2}/(\omega\mu\sigma^{\Sigma})$, which furnishes the wavenumber, $k_{\perp}$, of propagation transverse to the sheet; thus, $k_{m}/k=\sqrt{1-4k^{2}/(\omega\mu_{0}\sigma^{\Sigma})^{2}}$, and $|k_{m}/k|\gg 1$ in the nonretarded regime (Definition 2.1) [29]. By imposing $\text{Im}\,k_{\perp}>0$ according to the radiation condition at infinity vertically to the sheet, we obtain $\text{Im}\,\sigma^{\Sigma}>0$ which expresses the metamaterial character of the sheet.

An additional contribution to $u(x)$, not included in (7), is the radiation field due to the edge. This contribution is negligible in the region considered here and, thus, is omitted.

Near-edge field

Another region of significance is the vicinity of the edge, where $|(k_{m}/k)x|\ll 1$. In this region, the SPP interferes with the incident, reflected and radiation fields to yield a vanishing $x$-component of the (total) electric field on the sheet. Hence, as $x\searrow 0$, the solution, $u$, of integral equation (6) cannot be separated into the distinct (physical) contributions that are evident in the field far from the edge. Instead, we analytically compute the asymptotic behavior [29]

[TABLE]

This formula manifests the singular behavior of the field at the edge.

3 Numerical approach

In this section, we describe an implementation of the finite element method in order to solve the boundary value problem of section 2.1 based on curl-conforming Nédélec-elements [35]. First, we state boundary-value problem (1) in terms of a weak formulation that embodies discontinuity (2) of the tangential magnetic-field implicitly by means of an interior interface integral. Second, we discuss a number of numerical aspects to solve the weak formulation; in particular, the use of a specifically tuned PML, as well as a local refinement strategy to resolve the SPP.

3.1 Variational formulation and discretization

By a standard manipulation in electromagnetic theory, the substitution of the first equation of (1) into the second one yields a second-order partial differential equation for the vector-valued electric field, $\boldsymbol{E}$. To rescale the resulting equation to a desired, dimensionless form, we choose to set $k=1$, or, equivalently, scale the vector position $\boldsymbol{x}$ by $1/k$, in the spirit of section 2.2. Thus, we introduce the following related scalings [26]:

[TABLE]

Accordingly, the equation for $\boldsymbol{E}$ now reads

[TABLE]

The multiplication of the last equation with a test function $\boldsymbol{\varphi}$ and subsequent integration by parts yield [26]

[TABLE]

where $E^{\text{sc}}$ is the scattered field and $E^{\text{in}}$ denotes the incident field (section 2.1). An appropriate trial and test space for the weak formulation is [26],[33, Theorem 4.1]

[TABLE]

where $L^{2}$ denotes the space of measurable and square integrable functions. By choosing this space, the formerly strong interface condition (2) is now naturally embedded in the variational formulation. In more detail, the statement

[TABLE]

is a consequence of $\boldsymbol{E}\in\boldsymbol{H}(\text{curl};\Omega)$, and the jump condition

[TABLE]

is enforced by the term $-\,i\,\int_{\Sigma}{\rm Space\leavevmode\nobreak\ Sci.\leavevmode\nobreak\ Rev.}\boldsymbol{E}_{T}\cdot\bar{\boldsymbol{\varphi}}_{T}\,{\mathrm{d}}o_{x}$ in the variational formulation [26].

We implement (12) by using curl-conforming Nédélec-elements [35] with the help of the finite element toolkit deal.II [2]. The computational domain $\Omega$ was discretized with a quadrilateral mesh. In order to allow for local refinement, we use the well-known concept of hanging nodes (see, e. g., [6] for an overview) to relax the usual mesh regularity assumptions [10]. The resulting system of linear equations is solved with the direct solver Umfpack [12].

3.2 Perfectly matched layer for SPPs

Next, we discuss a construction of a PML [5, 9] for the rescaled Maxwell equations with a jump condition in connection to the boundary problem of section 2.1 [26]. The concept of a PML was pioneered by Bérenger [5] and can be viewed as a layer with modified material parameters ($\tilde{\varepsilon}_{r}\,$, $\mu_{r}\,$) placed near the boundary of the computational domain; cf. Figure 2. The core idea is to tune the material parameters inside the PML in such a way that all outgoing electromagnetic waves decay exponentially with no artificial reflection due to truncation of the domain. The PML is an indispensable tool for truncating unbounded domains for time-harmonic Maxwell’s equations, and other, akin partial differential equations, and is often used in the numerical approximation of scattering problems [33, 9, 5].

We follow the approach to a PML for time-harmonic Maxwell’s equations discussed in [9]. The idea is to use a formal change of coordinates from the computational domain $\Omega\subset\mathbb{R}^{3}$ with real-valued coordinates to a domain $\acute{\Omega}\subset\{z\in\mathbb{C}\>:\>\text{Im}\,z\geq 0\}^{3}$ with complex-valued coordinates (and non-negative imaginary part) [33]; and then transform back to the real-valued domain. For simplicity, we assume that the interface $\Sigma$ is parallel to the unit vector $\boldsymbol{e}_{r}$ within the PML (i.e., the normal $\boldsymbol{\nu}$ is orthogonal to $\boldsymbol{e}_{r}$). For details, we refer the reader to [26]. This procedure results in the following modified material parameters $\tilde{\varepsilon}_{r}\,$, $\mu_{r}\,$ and ${\rm Space\leavevmode\nobreak\ Sci.\leavevmode\nobreak\ Rev.}$ within the PML:

[TABLE]

In the above, we introduced the $3\times 3$ matrices

[TABLE]

where

[TABLE]

for an appropriately chosen scaling factor $s(\tau)$ that will be defined later. Note that $T_{\boldsymbol{e}_{x}\boldsymbol{e}_{r}}$ is the matrix that rotates $\boldsymbol{e}_{r}$ onto $\boldsymbol{e}_{x}$, and $\tau$ is the distance from the origin. The PML is assumed to be a spherical shell starting at distance $\rho$ from the origin, as shown in Figure 2.

3.3 Adaptive local refinement

By the assumption that $\text{Im}\,{\rm Space\leavevmode\nobreak\ Sci.\leavevmode\nobreak\ Rev.}>0$ in the nonretarded frequency regime (Definition 2.1), the SPP has a wavelength much smaller than the one in the ambient medium (at the same frequency, $\omega$). Thus, wave propagation along the metamaterial sheet, $\Sigma$, has a pronounced two-scale character, being characterized by length scales of the order of $1/(\text{Re}\,k_{m})$ and $1/k$; here, $\text{Re}\,k_{m}\gg k$.

In our numerical simulation, we use typical values of ${\rm Space\leavevmode\nobreak\ Sci.\leavevmode\nobreak\ Rev.}$ in the nonretarded regime for which the SPP wavelength is one to two orders of magnitude smaller than the wavelength in the surrounding medium. This poses a challenge because, on the one hand, the minimal computational domain (that still has a well-controlled error in slow oscillating modes) is limited by the free-space wavelength, $2\pi/k$; on the other hand, the minimal resolution necessary to resolve SPPs scales with $1/|k_{m}|$. Accordingly, in order to minimize computational cost while ensuring that the SPP is sufficiently resolved, we use an adaptive, local refinement strategy based on the dual weighted residual method [3].

Next, we outline the basics of our strategy. Starting from a relatively coarse mesh, the resolution is successively improved by a number of iterative refinement steps where a subset of cells is chosen for refinement. The selection of cells for refinement is made with the help of local (per cell $K$) error indicators, $\eta_{k}$, that are given by

[TABLE]

The cell-wise residuals, viz., the integrals [3, 26]

[TABLE]

are multiplied by the weights

[TABLE]

In the above, $\boldsymbol{E}_{h}$ denotes the finite element approximation on $\boldsymbol{E}$. The weights, $\omega_{K}$ and $\omega_{\partial K}$, are in turn computed with the help of: the solutions, $z_{K}$, of a “dual” problem and their respective interpolants, $I_{h}z$, in the finite element space. The rationale of using a dual problem for computing the requisite weights is that these can be tuned to a quantity of interest in the form of a (possibly non-linear) functional [3]. In our particular application, we choose to use

[TABLE]

as the quantity of interest. Hence, the dual problem reads

[TABLE]

4 Validation of numerical method and further numerics

In this section, we focus on numerical computations by our finite element method. These computations have a two-fold purpose: validation of our numerical method and extraction of further predictions. First, we address the prototypical geometry of the semi-infinite metamaterial sheet in order to validate and verify our numerical approach by comparison of simulations against the analytical description of section 2.2. Second, we numerically simulate wave propagation under an incident plane wave in a number of realistic geometries which are relevant to nanophotonics applications; to our knowledge, no analytical results are available for these geometries in the existing literature. For the numerical experiments, we use $E_{x}^{\text{i}n}=i$ (i.e., $E_{0}=i$ and $\alpha=0$) throughout.

Figure 2 depicts the geometries used in our numerical tests. We examine the following configurations.

• –

The reference case, i.e., the semi-infinite metamaterial sheet [Figure 2(a)]. Our numerical simulations for this setting are compared against an analytical solution. In particular, we verify the quality of approximation (7) for the $x$-component of the electric field.

• –

Two co-planar, symmetrically placed, semi-infinite metamaterial sheets with a gap [Figure 2(b, c)]. In this configuration, the incident plane wave excites an SPP on each sheet. If the edges of the sheets are sufficiently close to each other, the induced SPPs may interfere destructively. We numerically examine the relative amplitude of the resulting SPP (compared to the reference case) on one sheet as a function of the gap width, $d$. We consider two different cases for the material of the gap. First, the gap is empty (Figure 2b); and, second, the gap is filled with a highly conducting material (Figure 2c). Accordingly, we show numerically that the dependence of the SPP amplitude on $d$ is dramatically different in these cases.

• –

A finite strip of metamaterial [Figure 2(d)]. For our numerics, the width, $d$, of the strip is chosen to have small to intermediate values compared to $1/k$. Because of a standing wave formed on the strip, we expect that a strong resonance effect can occur for suitable values of $d$. We numerically verify this resonance and quantify the (maximal) SPP magnitude as a function of $d$.

4.1 Semi-infinite strip

First, we consider the geometry of the semi-infinite strip, $\Sigma=\{(x,0)\in\mathbb{R}^{2}\,:\,x\geq 0\}$, depicted in Figure 2(a). The primary purpose of our numerical simulation in this setting is to validate the analytical results discussed in section 2.2. In particular, by subtraction of the incident field, approximate formula (7) yields a corresponding expression for the $x$-component of the scattered electric field, $E^{\text{sc}}_{x}$. In terms of the rescaled quantities (section 3.1), the scattered field reads

[TABLE]

where the first term (with brackets) represents the directly reflected field and the second term is the SPP. Here, for $\alpha=0$, the requisite values of $\widehat{Q}_{+}(\xi)$ in the exponent are $\widehat{Q}_{+}(0)\approx 0.0005+0.05i$, and $\widehat{Q}_{+}(k_{m}/k)\approx 0.346574+0.392699i$. Note that in the case with the surface conductivity ${\rm Space\leavevmode\nobreak\ Sci.\leavevmode\nobreak\ Rev.}=2.0\cdot 10^{-3}+0.2i$ (${\rm Space\leavevmode\nobreak\ Sci.\leavevmode\nobreak\ Rev.}=8.89\cdot 10^{-4}+0.133i$), the SPP wavenumber is $k_{m,r}\;=\;10.0489+0.0994937i$ ($k_{m,r}\;=\;15.0701+0.100288i$).

In Figure 3, we compare graphically the outcome of formula (18) to results of our numerical method for the corresponding scattered field. Specifically, the real and imaginary parts of the scattered electric field in the $x$-direction, $E_{x}^{\rm sc}$, are plotted as a function of the spatial coordinate, $x$, along the sheet, $\Sigma$. Evidently, the numerical and analytical results are in excellent agreement outside the PML and for (roughly) 2-3 SPP wavelengths away from the origin ($2\leq x\leq 20$ in our numerics).

Furthermore, we test the finite element numerical simulations against analytical prediction (9) for the asymptotic behavior of $E_{x}$ along the sheet near the edge, as $x$ becomes sufficiently small. The predicted behavior is confirmed numerically in Figure 4, in which the real part of this field component is plotted versus $x$ closely enough to the edge of the sheet.

4.2 Co-planar, semi-infinite sheets with a gap

Next, we examine the (symmetric) configuration that consists of two co-planar, semi-infinite sheets with a gap of width $d$, $\Sigma=\{(x,0)\in\mathbb{R}^{2}\,:\,x<-d/2\ \mbox{or}\ x>d/2\}$; see Figure 2(b,c). This setting is of fundamental interest in applications because experimental setups often involve arrays of strips of metamaterials such as graphene [4]. By using the present geometry, we are able to isolate the influence of the gap width on the surface plasmon of the metamaterial, since each sheet is infinitely long. In particular, resonances associated with standing waves on the sheet are absent from this setting (see section 4.3 for the resonant, finite strip).

We consider two cases for the material of the gap. In one case, the gap is empty, or, in other words, the material of the gap is identified with the medium of the ambient space [Figure 2(b)]; and in another case the gap is filled with a highly conducting material [Figure 2(c)]. Because of the proximity of the two edges, the SPPs induced on each sheet are in principle expected to interfere. Our numerics show that the closer the sheets are to each other, i.e., as the gap width, $d$, decreases, the stronger the effect of destructive interference of the SPPs is. For example, this trend is demonstrated in Figure 5, which depicts the SPP versus $x$ near the empty gap for three values of $d$.

We perform a parameter study in order to examine the strength of the SPP as a function of the gap width. We examine the cases with an empty gap and a gap filled with a material of (scaled) surface conductivity $\sigma_{r,2}^{\Sigma}$ such that $\text{Re}\,\sigma_{r,2}^{\Sigma}\gg 1$, $\text{Im}\,\sigma_{r,2}^{\Sigma}=0$. In this study, we carry out the fitting of the finite element-based numerical solution for the $x$-component of the scattered electric field to an analytical formula constructed from (18). Specifically, this formula reads

[TABLE]

which is supplemented with two parameters, namely, the extra factor, or relative amplitude, $A$ for the SPP and a corresponding phase shift, $\varphi$. These parameters, $A$ and $\varphi$, are determined via fitting of the last expression for $E_{x,r}^{\rm sc}$ to our numerics.

Interestingly, in the two distinct cases of the gap material mentioned above, we observe two fundamentally different scaling laws for $A$, which hold for a range of values of $d$. These behaviors are shown in a semi-logarithmic plot of the relative amplitude, $A$, versus the gap width, $d$; see Figure 6. In the case with an empty gap, we observe that the relative amplitude, $A$, exhibits a logarithmic behavior with $d$, viz., $A\sim-\log\,d$. In contrast, in the case with a highly conducting material, the corresponding scaling is approximately linear, viz., $A\sim d$. Both scaling laws hold at least for $2^{-8}\leq d\leq 2^{-4}$ (Figure 6). Outside this region for $d$, our numerical fitting yields results consistent with the limits $\lim_{d\to\infty}A=A_{\infty}$ and $\lim_{d\to 0}A=0$, where $A_{\infty}$ is the (constant) relative amplitude of the SPP for an infinitely large gap, when the two sheets tend to be isolated from each other.

4.3 Finite strip

Next, we consider the geometry of a finite metamaterial strip, $\Sigma=\{(x,0)\in\mathbb{R}^{2}\,:\,-d/2\leq x\leq d/2\}$ [Figure 2(d)]. In this case, if the length $d$ of the strip is sufficiently large, the two edges are expected to cause resonances to the observed SPP. In other words, the maximal magnitude of the SPP generated in this structure should exhibit peaks as a function of the strip width, $d$ (for fixed $k$).

In our numerical simulations, we compute the resulting, maximal magnitude of the SPP as a function of $d$; see Figure 7. For these numerics, we choose the surface conductivity to be equal to ${\rm Space\leavevmode\nobreak\ Sci.\leavevmode\nobreak\ Rev.}=2.0\cdot 10^{-3}+0.2i$. We select two peaks among those observed in the chosen range of values of $d$. These resonances are depicted in separate plots here [Figure 7(a,b)]. In order to quantify the two selected resonances, we carry out the fitting of the numerical results in each case to the Lorentzian function

[TABLE]

The parameters $d_{0}$, $\gamma$ and $A$ of this formula are determined via the fitting procedure. Specifically, we obtain the following values: $d_{0}=1.372\,\lambda_{m}$, $\gamma=0.015\,\lambda_{m}$, $A=1.34$ for the first resonance [Figure 7(a)]; and $d_{0}=10.375\,\lambda_{m}$, $\gamma=0.11\,\lambda_{m}$, $A=1.49$ for another (tenth) resonance [Figure 7(b)] in the chosen range of values of $d$. Here, $\lambda_{m}=2\pi/(\text{Re}\,k_{m})$ is the SPP wavelength on the semi-infinite sheet of the same conductivity (see section 2.2).

5 Conclusion

In this paper, we numerically studied the generation of SPPs on atomically thick metamaterial sheets by edge effects in different geometries. Our chosen configurations included: a semi-infinite sheet; two co-planar, semi-infinite sheets with a gap of variable width; and a finite strip of variable width. In our computations, we used an adaptive finite element method with curl-conforming Nédélec-elements in order to resolve the fine scale of the SPP propagating along the sheet in the presence of edges. Our numerical approach here forms an extension of the method introduced in [26] to more realistic geometries.

We validated our numerical treatment by comparison of the finite element-based numerics to an analytical solution for the semi-infinite sheet [29]. By further numerical simulations, we demonstrated that the presence of a highly conducting material in the gap between two co-planar, semi-infinite metamaterial sheets can increase the SPP magnitude, and leads to a distinctly different dependence of the SPP on the gap width in comparison to the case with an empty gap. In addition, we numerically characterized SPP resonances on a finite strip, and demonstrated that the SPP maximal magnitude (as a function of gap width) is well described by a Lorentzian function.

Our results point to a few open problems. For instance, our computations have focused on 2D, although our approach is applicable to arbitrary spatial dimensions. An emerging question concerns the use of truly three-dimensional (3D) geometries, which may contain metamaterial sheets with corners, conical singularities or arbitrarily curved hypersurfaces. The generation and propagation of SPPs in 3D settings is the subject of ongoing work.