{'fraction_non_alphanumeric': 0.03913307159500782, 'fraction_numerical': 0.01659804015448945, 'mean_word_length': 4.672641680246243, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Adaptive finite element simulations of waveguide configurations involving parallel 2D material sheets

We discuss analytically and numerically the propagation and energy transmission of electromagnetic waves caused by the coupling of surface plasmon polaritons (SPPs) between two spatially separated layers of 2D materials, such as graphene, at subwavelength distances. We construct an adaptive finite-element method to compute the ratio of energy transmitted within these waveguide structures reliably and efficiently. At its heart, the method is built upon a goal-oriented a posteriori error estimation with the dual-weighted residual method (DWR). Further, we derive analytic solutions of the two-layer system, compare those to (known) single-layer configurations, and compare and validate our numerical findings by comparing numerical and analytical values for optimal spacing of the two-layer configuration. Additional aspects of our numerical treatment, such as local grid refinement, and the utilization of perfectly matched layers (PMLs) are examined in detail.


Introduction
Graphene is a two-dimensional carbon allotrope with one-atomic thickness that is arranged in a honeycomb lattice structure [1]. It has a wide potential for applications in nanophotonics due to a number of desirable electronic and optical features, such as extreme confinement, low losses, and tunability [1,2]. In the infrared regime, the electric surface conductivity of such a 2D material is characterized by being complex-valued with a dominant positive imaginary part. This allows for the propagation of surface plasmon polaritons (SPPs), which are highly confined to the 2D material and slowly decaying electromagnetic waves. Within the frequency domain of interest, the wavelength of a SPP is up to two orders of magnitude smaller than that of the exciting ambient light.
Waveguide structures that enable subwavelength confinement of the optical modes are of great importance in nanophotonics [3,4,5,6]. Traditionally, waveguides have been implemented through local modulation of the shape and/or refractive index profile of the optical dielectric medium [7]. These dielectric waveguides, however, are restricted by the diffraction limit of light, λ 0 /n, where λ 0 is the wavelength in free space and n is the refractive index [4]. SPPs, on the other hand, can be confined within a very small area beyond the diffraction limit of light, and can be used as an information carrier for highly-integrated photonic circuits [8]. Thus, a number of plasmon-based waveguides have been proposed in the past decade, such as metallic nanowires [9,10], metallic nanoparticle arrays [5,11], hybrid plasmonic waveguides [6], and gain assisted plasmonic waveguides [12].
Although some studies that investigated the plasma modes and optical SPP modes of a double-layer graphene were published [13,14,15], the majority of research on graphene has been on a single-layer systems consisting of a single, planar sheet of graphene [16,17,18]. By introducing a second, parallel sheet, placed at a small but finite distance to the other sheet, it is possible to drastically change the confinement and propagation characteristics of SPPs. The purpose of this paper is two-fold.
First, we numerically investigate an infinite 2D waveguide by computing a finite element approximation for the solution of the corresponding scattering problem governed by time-harmonic Maxwell's equations. This prototypical geometry is motivated by proposed waveguide configurations that include graphene layers, or carbon nanotubes as an integral part of their design [19]. By adjusting the confinement of the two-layer system, we find an optimal spacing for which the coupling of the SPPs is maximal. To this end, a goal-oriented mesh refinement strategy and a perfectly matched layer are utilized.
Secondly, we derive and discuss an integral equation describing the timeharmonic electromagnetic field of a double-layer system. We examine the contributions from a pole of the scattered field solution, which are responsible for the generation of SPPs. We demonstrate that our findings are in accordance with those from a single-layer system by observing the evolution of the scattered field under different interlayer spacings. In conclusion, we note that the maximal spacing found in the numerics is in agreement with the value found in the analytical expressions.

Related works
The SPP dispersion of a single graphene layer and a single graphene layer deposited on dielectric substrates has been extensively investigated by many authors [20,16,18]. Additional confinement and a change of propagation characteristics can be achieved by stacking another layer of these 2D materials on top of a single-layer [21]. By extension, theoretical aspects of bilayer graphene [15,22,13], multilayer graphene [23], and intercalated graphite [24] have been studied recently. For example, [15] finds a frequency for which reflection in a double-layer graphene system, with both equal and different surface conductivities, is zero, leading to exponentially amplified transmitted modes. [25,26] shows that for plasmonic crystals, which consist of stacked metallic layers arranged periodically with subwavelength distance, embedded in a dielectric medium, the TM polarized waves experience an effective dielectric function that combines a bulk energy of the microstructure of the ambient dielectric medium and surface average of the surface conductivity of each sheet. Homogenization of layered structures and extension to a general hypersurface are also discussed.
However, a rigorous numerical and analytical treatment of waveguide configurations involving time-harmonic Maxwell's equations is not of primary interests in these publications. This paper aims to address three points. First, we introduce a reliable and efficient numerical method to readily compute propagation characteristics of the SPPs in a double-layer structure. Second, we validate the numerical findings against an analytical solution. Third, we demonstrate that our numerical approach is easily extensible to different computational domains.

Paper organization
The paper is organized as follows. In Section 2, we derive the variational formulation that serves as the basis for our numerical and analytical investigation. In Section 3 we develop the numerical framework, including a goal-oriented mesh adaptation based on the dual weighted residual (DWR) method, and a perfectly matched layer (PML). We perform a direct numerical simulation of a prototypical two-layer system in Section 4, and identify optimal spacings of two-layer systems for maximal transmission. In Section 5, we derive analytic solutions for the two-layer system and validate our numerical findings against them. Finally, Section 6 concludes the paper with a summary of our results and an outlook.

Variational formulation
In this section, we lay out the variational formulation for time-harmonic Maxwell's equations with an interface condition. We introduce a rescaling for time-harmonic Maxwell's equations that will ease the numerical computation of SPPs [27].

Preliminaries: Maxwell's equations
Time-harmonic Maxwell's equations with an electric current density read [28,29]: Here, E(x) and B(x) denote the electric and magnetic field, respectively. µ(x) and ε(x) are complex-valued rank 2 tensor quantities, where µ(x) denotes the magnetic permeability and ε(x) denotes the electric permittivity. The vectorvalued quantity J a (x) is an externally applied, electric current density. We assume a general x dependence of all quantities with some (weak) regularity conditions to ensure unique solvability that will be stated later. The (constant) temporal frequency ω > 0 arises from the time-harmonic nature of the solution fields, i. e., a solution of (1) is a special solution of the general time-dependent Maxwell's equations by rewriting all vector-valued components F by We are interested in simulating waveguide configurations that feature 2D material sheets. The sheets are modeled as an idealized hypersurface Σ with an effective surface conductivity σ(x) defined on Σ [30,31,2]. In general, Σ shall consist of two parallel, possibly curved, conducting sheets separated by a fixed distance d; see Figure 1. The discontinuity along the surface due to the conductivity leads to a jump condition in the tangential component of the magnetic field [30,28,2]: where ν is a fixed normal vector field associated with Σ; the symbol [ . ] Σ denotes the jump over Σ with respect to ν, and the subscript T denotes the tangential part of the respective vector, F T = (ν×F )×ν. Under appropriate conditions on σ(x), jump condition (3) generates SPPs on the interface [31,2]. We make the assumption that ε(x) and µ(x) become homogeneous and isotropic for large |x| and impose the Silver-Müller radiation condition at infinity [27], viz., Here, c = 1/ √ εµ denotes the speed of light. The explicit inclusion of this condition is omitted in our numerical simulation by incorporating an appropriate boundary condition and a PML.

Rescaling and variational formulation
Numerical values in SI units for solutions of (1) are many orders of magnitude apart. Further, the typical length scale of SPP is one to two orders of magnitude smaller than the free-space wavelength k 0 [2]. These discrepancies Ω Σ ν ν J a n ∂Ω Figure 1: Schematic of the computational domain, Ω, with boundary ∂Ω, outer normal n, and normal field ν defined on the waveguide Σ. An electric Hertzian dipole, Ja, is placed inside a prescribed waveguide structure, Σ.
in the magnitude of the length scales makes the direct numerical simulation of the SPPs difficult. As a remedy we use a rescaling to dimensionless units that normalizes the length scale by the free-space wavenumber, k 0 , [27]: Here, ε 0 and µ 0 denote the vacuum permittivity and permeability, respectively. This leads to two distinctly separated scales: one related to the free-space wavenumber ∼ 1, and another for the SPP wavenumber, k SPP ∼ 1/σ r , on the conducting sheets [27]. The rescaled, dimensionless form of time-harmonic Maxwell's equation (1) reads with the jump condition Let the computational domain Ω ⊂ R 3 be bounded, simply-connected, and Lipschitz-continuous with piecewise smooth boundary ∂Ω. Assume that Σ is a Lipschitz-continuous, piecewise smooth boundary. The Silver-Müller radiation condition (5) reads Multiplying (6) with the complex conjugate of a test function ϕ and subsequent integration by parts using (7) and (8) recovers the corresponding weak formulation: In the above, L 2 (·) 3 denotes the space of vector-valued square integrable functions, and H(curl) is the subspace of L 2 (·) 3 consisting of square integrable functions whose (distributive) curl admits a representation by a square integrable function. Equation (9) will serve as a starting point for a finite element discretization. For elementary results on existence and uniqueness, we refer to [27,28,32,33].

Numerics: Computational domain and discretization scheme
In this section, we briefly present the geometry and numerical tools used in the computations. A PML is introduced and its role in negating the undesired effects of the absorbing boundary condition of the plasmon modes is described. Additionally, a local adaptive mesh refinement strategy is presented that captures the highly oscillatory behaviors of the plasmons near the interfaces. The variational formulation (9) is discretized on a non-uniform quadrilateral mesh with higher-order, curl-conforming Nédélec elements [34]. Such a choice is ideal where the weak jump condition is naturally treated by aligning with the mesh. Let X h (Ω) ⊂ X(Ω) be a finite element subspace spanned by Nédélec elements. Then under a sufficiently refined initial mesh, the variational formulation is uniquely solvable for E h ∈ X h (Ω) and for all ϕ ∈ X h (Ω).

Geometry
In this paper we will study a prototypical geometry consisting of two flat, conducting layers in a square domain. The two layers are arranged parallel to each other with distance d apart; see Figure 3. This prototypical geometry is motivated by proposed waveguide configurations that include for example concentrically arranged carbon nanotubes as integral part of their design [19] (see Figure 2. Even though our waveguide configuration is quite simple in comparison, we make the claim that due to the dominance of the SPP interaction of the two layers we actually capture the quantitative behavior of the two-layer interaction quite well. Our computational framework has thus the potential of guiding the design of more complicated waveguide structures in the future. Two different current sources are considered: a single vertical Hertzian dipole placed at the midpoint of the two sheets (see  configuration with one dipole placed directly on each sheet (see Figure 3c). The use of such dipole sources in order to excite the desired SPP modes in the waveguide configuration is again an idealization. However, due to the large magnitude of the excited SPP modes, the influence of the dipole is neglibible already after a short distance from the source. Only the dominant travelling mode is observed, which enables us to study essentially "free" coupled SPP structures in the waveguide.
In the numerical computation we regularize dipoles to a small but finite thickness. This necessitates a change of sign in the double dipole configuration when crossing over the interface Σ. See Section 4.1 for details. In the case of the double dipole excitation, the proximity of the sources to the interface renders the SPP mode that is excited to dominate [30]. For the single dipole excitation, however, the strength of the SPP vanishes exponentially with the distance d. For this reason, we will primarily focus on numerical results obtained with double dipole excitation, 2(b).
The values x tr and x re are the transmission and reception locations, respectively, at which the tangential component of the energy is to be measured. In principle, a number of choices for the quantity of interest, J (E), are possible. In the following, we use a non-linear quantity of interest given by an energy transmission ratio: The numerator computes the transmitted energy, measured at a vertical strip located sufficiently far from the source, and the denominator calculates the received energy. The integrands in J (E) are modified by weight functions that localize the integral to a vertical strip where we measure the field intensity for transmission and reception. The choice (11) for the quantity of interest leads to a localized right-hand side J of the dual problem that is sensitive to the highly oscillatory SPPs associated with the electric field, E.
Remark. The functional J (E) is continuously differentiable as long as the denominator is different from zero. This is indeed the case for the solution E and all approximations E h for our choice of geometry and dipole excitation.

Perfectly Matched Layer
A perfectly matched layer (PML) is a truncation procedure motivated from electromagnetic scattering problem in the time domain. The underlying idea of a PML is to surround the computational domain with an artificial sponge layer such that all outgoing electromagnetic waves decay exponentially with minimal artificial reflection [35,36,28].
As outlined in [27,28], we carry out a change of coordinates from the computational domain with real-valued coordinates to a domain with complex-valued coordinates. Projecting back to the real coordinates yields again system (1) with (3), but with modified material parameters (ε r , µ −1 r , σ r ) inside the PML. We refer the reader to [27] for details.
The PML can be implemented by suitably replacing (ε r , µ −1 r , σ r ) within the PML. For a spherical absorption layer we define the matrices Here, r denotes the distance to the origin, s(τ ) is an appropriate nonnegative scaling function that will be defined later, T exer is the rotation matrix that rotates e r onto e x . The material parameters are hence transformed inside the PML as follows:

A posteriori error estimation and local refinement
One of the computational challenges of our problem is the need for a much finer mesh refinement near the interfaces Σ 12 and Σ 23 in order to resolve all small scale SPP structures. We discuss now an efficient adaptive refinement scheme utilizing an a posteriori error estimator based on the dual weighted residual (DWR) method [37].
Consider the following dual problem: Find a solution Z ∈ H(curl; Ω) such that for all ϕ ∈ X(Ω), where J (E) is a quantity of interest mapping The dual solution Z encodes how the target error quantity depends on local properties of the data [37]. Next, we define local error indicators with the help of the solutions E and Z of the primal problem (9) and dual problem (15), respectively [27], [37, Prop. 2.1]: Here, ρ Q and ρ * Q denote the primal and dual cell-wise residual, respectively, associated with variational equations (9) and (15): where A( . ) is given in (10). Here, χ Q denotes the indicator function associated to Q, that is, χ Q (x) is 1 inside the cell Q, and 0 otherwise. The local error indicator η Q given by (17) can now be approximated and used in a local refinement strategy [37].
Remark. The remainder term R is cubic in the error E − E h and can therefore generally be neglected [37]. More precisely, for our particular choice of quantity of interest (11) a lengthy calculation reveals provided that the numerator in (11) is smaller than the denominator (which is true for our choice of geometry). Given the fact that we place the measurement position x tr close to the source we conclude that |R| is well controlled and small in our case. Our goal is an optimal local refinement for the numerical simulation of energy propagation of the SPPs at the location of our choosing. Consequently, the weight Z − Z h in residual (18) is generally large near the interface and at points where the influence of the solution on quantity (11) is high.
In practice, the numerical evaluation of (18) and (19) is typically done with the use of a higher-order approximation for the dual solution Z and the primal solution E. However, such a calculation of a higher-order approximation is computationally costly. We therefore use a patch-wise projection π (2) 2H Z h to a higher-order space on a coarser mesh level [38]:

Direct numerical simulation
In this section we present computational results for the two-layer system that was introduced above. We demonstrate numerically how the (effective) wavenumber of SPP structures depends on the interlayer distance d and investigate the functional relationship of the energy transmission ratio to the interlayer distance d. We determine the optimal spacing, which will later be used to compare against the analytical findings from Section 5. We validate our local refinement strategy by comparing the convergence rates with uniform refinement, and demonstrate the effectiveness of the numerical tools discussed in Section 3. All numerical computations are carried out with the finite element library deal.II [39].

Setup and discretization parameters
We consider a vertical electric dipole positioned at a 1 = (−0.7, 0) (for single dipole excitation), and at a 2/3 = (−0.7, ±d/2) (for the double dipole excitation). The current density J a is thus given by for single, or double dipole excitation, respectively. We use two values for the surface conductivities, σ Σ r,12 = σ Σ r,23 = σ r , σ r = 0.002 + 0.2i, and 0.002 + 0.15i, that are within realistic parameter ranges [27]. The computational domain, Ω, is chosen to be a square with edge length 4. A spherical PML is enforced for ρ > 1.6. Following [28,27], we regularize the Dirac deltas in the current density as follows, for x − a i 2 < r d , and 0 otherwise. The signum function ensures that the regularized dipole changes sign whenever the regularization crosses the conducting layer Σ 12 , or Σ 23 . We choose a fixed value of r d = 10·2 −12 throughout the paper. We set the position at which we evaluate (11) to x tr = −0.65 and x re = 0.75. This choice maximizes the distance |x re − x tr | while ensuring that evaluation points are sufficiently far away from the regularized dipole sources and the PML. We use the following scaling function s(ρ) for the PML [27] and set the free parameter to s 0 = 0.05 in our computations.

Validation of local refinement strategy
We validate our numerical framework by comparing values for the quantity of interest (11) obtained by extrapolating numerical values under uniform, and under local refinement. The computations were performed for d = 0.20. The results are shown in Figure 4. The data is fitted to the curve f (x) = a + b x c and extrapolated. For local refinement, the fit parameters we obtain are a local = 0.205333, b local = 1.25229×10 8 , c local = 1.57811, whereas for uniform refinement we get a unif = 0.206216, b unif = 85.8656, c unif = 0.612326. As expected [27], we obtain a much faster convergence rate (c ≈ 1.6) in the quantity of interest for local refinement as opposed to uniform refinement (c ≈ 0.6). We conclude that our computation of the energy transmission ratio with J (E) ≈ a local is reliable within 1%.

Optimal spacing
Next, we preform a parameter study of the energy transmition ratio for varying interlayer spacing d ranging from d min = 4 · 10 −12 to d max = 0.2, where d max corresponds to about 1/3 of the single-layer SPP wavelength, 2π/Re(k m,r ), or 1/30 of the free space wavelength [27]. The (interpolated) results are shown in Figure 5. We used a releatively coarse initial mesh for all computations with around 20 thousand degrees of freedom. After 12 local refinement cycles

1/DoF
Convergence rate uniform refinement local refinement uniform refinement fit local refinement fit Figure 4: Convergence of the energy tranmission ratio for uniform and local refinement for d = 0.20. The relative error of the transmission ratios obtained with both refinement strategies, respectively, are plotted against the inverse of the number of degrees of freedom. We observe a convergence order of about c ≈ 0.6 for uniform refinement and c ≈ 1.6 for local refinement. The reference value J E * was obtained by taking a weighted average of the asymptotic transmission ratios, c local and c uniform .
using the adaptive refinement procedure outlined in Section 3.3, we reached roughly 2 million degrees of freedom on the finest mesh. We make a qualitative comparison of three representative cases: Figure 6 shows the real part of the computed scattered electric field in x-direction, Re(E sc x ), for SPPs on a singlelayer system, and a two-layer system with d = d opt = 0.05245 and d = 0.1805, respectively. In the case of optimal spacing, d = d opt , the wavenumber of the excited SPPs is roughly twice as large as the one obtained for the singlelayer case, cf. Figures 6a and 6b. For large enough d, for example d ≈ 0.18, we observe that the wavenumber of the excited SPP approaches the singlelayer case, cf. Figures 6a and 6c. This indicates that the two-layer system is converging to the single-layer setting, analogous to the behavior observed for the energy transmission ratio.

Analytic solution and validation
In this section, we derive an analytic solution for (1) and (3) for an (idealized) infinite two-layer system with single dipole excitation; see Figure 3b and Figure 7). We identify the limiting behavior for the case of large interlayer spacing d, compute the effective wavenumber of the dominant SPP mode and  validate our numerical findings with these results. For better readability and ease of comparison [30], we revert back from the rescaled version (6) to the original form of Maxwell's equations (1). Let us consider two planar sheets in R 2 situated at y = 0 and y = d, respectively; see Figure 7. The conducting sheets separate R 2 into three regions: Region 1 ({y > d}) has wavenumber k 1 and shares the boundary with region 2 ({0 < y < d}), whose wavenumber is given by k 2 . Region 3 ({y < 0}) shall have wavenumber k 3 , where k 2 j = ω 2 ε j µ. Here, ε j denotes a complex-valued permittivity (j = 1, 2, 3). Let a vertical electric Hertzian dipole be positioned at (0, a) in between the interfaces, viz., J a = δ(x)δ(y − a)e y . Define the Fourier transform,F (ξ, y) of the vector-valued fields (F ≡ B, E) through the integral formula Applying the Fourier transform to Maxwell's equations (1) gives A number of elementary algebraic manipulations of (24) yield the differential equation  where we have set β 2 j = k 2 j − ξ 2 . We now make the following solution ansatz for the magnetic field obeying the Sommerfeld radiation condition [30]: y > d, a 2 e iβ2y + b 2 e −iβ2y − ξµ 2β2 e iβ2|y−a| , 0 < y < d, b 3 e −iβ3y , y < 0.
The remaining electric field components can be derived from the relationŝ Next, we determine closed expressions for the coefficients in (26) by matching with boundary conditions (3) on each interface: Ω: x y Region 1 (k 1 ) conducting sheet (conductivity σ 12 ) Region 2 (k 2 ) Region 3 (k 3 ) conducting sheet (conductivity σ 23 ) Here, the constants R ij,m are given by By substituting back into (26) and undoing the Fourier transform, all field components of E and B can be expressed as analytic integrals. In particular, we are interested in E 2x (x, y), the electric field component in x-direction between the two-layers:

Approximation of the pole contribution
Next, we obtain the scattered electric field in x-direction, E sc 2x , observed at y = 0 by subtracting the incident field ωµ from (28). Some additional minor rearrangement yields 1 − R 12 R 23 e 2iβ2d e iξx =: (I) + (II) + (III). (29) Each term of the integrand contains SPP contributions stemming from different conducting sheets: Term (I) and (II) arise from the SPP situated at Σ 23 and Σ 12 , respectively. Term (III) is a mixed term due to the interlayer coupling of SPPs.
We now discuss the role of simple poles in the evaluation of integral (29). For the sake of simplicity, let us now assume that k ≡ k 1 = k 2 = k 3 and σ ≡ σ 23 = σ 12 . Thus, R 12 = R 23 ≡ R and β j ≡ β, where Waveguide modes correspond to single poles of the integrand in integral (29) [30]. Inspecting (29) we see that these are exactly given by the condition Re iβd = ±1. The solution for the branch with the plus sign recovers even waveguide modes, and, correspondingly, the minus sign recovers odd waveguide modes [40].
For d sufficiently small, there is only a single dominant mode. We analyze this case further. The common prefactor of term (I) and (II) is given by .
The TM surface plasmon corresponds to the residue contribution to the electromagnetic field from the pole ξ = k B m , where k B m is a solution of the transcendental relationship for the mode, 2k 2 + ωµσβ(1 − e iβd ) = 0 [30]. Now, Here, the subscript p denotes evaluation at the pole. Each of I, II, III consists of the branch-cut contribution and the pole contribution. We omit the discussion of the branch-cut in this paper and focus only on the simple pole, ξ = k B m . This is because for an infinite conducting sheet, the SPP is identified with the part of the electromagnetic field equal to the contribution to the Fourier integrals of the simple pole that solves the above transcendental relationship. For a more thorough discussion on the branch-cut and its computation, we refer readers to [27]. By the residue theorem, .  And similarly for the second integrand term, .
The interlayer pole contribution is calculated in the same fashion. .

Limiting behavior and effective SPP wavenumber
For d 1 and fixed dipole position a, we expect the solution of the doublelayer system to approach the solution of the single-layer system. We justify this claim by observing that as d → ∞, which corresponds to the single-layer solution [27]. Further, the remaining terms (II) p and (III) p converge to zero due to the presence of e iβpd . Under the rescaling introduced in Section 2, equation (31) becomes We numerically solve the above via a root finding algorithm and plot the real and the imaginary part of the wavenumber k B m,r as a function of distance d r ; see Figure 8. For d r small, Im(k B m,r ) is greater than the single-layer counterpart. However, upon entering a regime where the interlayer coupling dominates, i.e., where k B m,r becomes less lossy than k M m,r , we observe that there is an optimal distance, d opt r ≈ 0.05538, at which the scattered field solution in the x-direction attains its maxium. For d r large enough, only the contribution from the bottom interface remains (term (II) p and (III) p vanish), and the wavenumber converges to that of the single-layer case.

Comparison and validation of numerical results
Finally, we compare the numerical results obtained in Section 4 to the analytical solution derived above.
In particular, we expect to observe that the contribution from the interlayer coupling term (III) p of the SPP dominates in the energy transmission ratio. In this vein, our direct numerical simulation is performed for the double-dipole excitation on the interfaces so that the (III) p dominates [27]; see Section 4. We postulate that the complex-valued wavenumber k B m,r given by (37) and associated with the SPP mode (III) p describes the effective transmission behavior of the two-layer system.
The optimal distance obtained from the analytical solution is computed to be d opt r ≈ 0.05538; see Figure 8. In order to test the validity of our numerical method, we compare this value against the optimal distance, d opt = 0.05245, obtained by numerically computed the energy transmission ratio; see Figure 5. Both values are in very good agreement. We attribute the small discrepancy of both values to the different current sources that were used.

Conclusion
In this paper, we extended a variational framework for the numerical simulation of the SPPs excited by a current-carrying source on an infinite conducting sheet to the SPPs generated by single/double excitations in a waveguide configuration. The conducting sheets, e.g. graphene, are modeled as idealized hypersurfaces that naturally takes into account the jump condition of the magnetic field.
We demonstrate that the interlayer coupling of the SPPs present in the twolayer system is responsible for higher confinement and less losses than the singlelayer system. We computed optimal interlayer spacings using two approaches. First, we compute via finite element simulations with double-dipole excitation. Second, we derive the pole contribution of the x-directed scattered field solution and solve for the wavenumber of the SPPs (37) numerically. The computed values are in agreement with one another. The numerical results on the energy transmission ratio are in very good agreement with analytic results obtained for the SPP mode of the double-layer structure.
Our numerics admit several generalizations and extensions. In particular, our variational framework can be readily used without modification to model any geometric configuration that is meshable by quadrilaterals. This flexibility enables numerical simulations of curved waveguides, systems consisting of multiple layers, or complicated optical devices in the near future.