Numerical investigation of mode characteristics of nanoscale surface plasmon-polaritons using a pseudospectral scheme

Chia-Chien Huang

doi:10.1364/OE.18.023711

1. Introduction

The ability for fabricating the integrated circuits is now at dimensions below 100nm. Further scaling the electronic circuits encounters some limitations such as operating bandwidth, signal delay, and energy dissipation [1,2]. For working with these problems, applying the photonic circuits offers an effective solution. However, using pure dielectrics to produce nanoscale photonic components is primarily restricted by the inherent diffraction limit. Therefore, it is difficult to make photonic integrated circuits fully compatible with the field of integrated electronic circuits. In recent years, exciting the surface plasmon polariton (SPP) is a promising way to squeeze the optical energy in nanometric cross section [3]. Theoretically, the SPP occurs at the interface between a dielectric and a metal due to the coupling of an electromagnetic wave and the free electrons of the metal [4]. Enlarging the SPP wavevector decreases the size of the resulting spatial confinement more than the conventional dielectric methods can. Many studies report various types of plasmonic structures to achieve the nanoguiding of light power [5–12]. The most essential characteristics of plasmonic waveguides include the dependences of mode confinement and propagation length on the geometric structure and the operating wavelength. Theoretically, the plasmonic waveguide structures need at least one interface composed of a dielectric and a metal to produce the SPP mode, and the SPP mode profile from the interface into the metal changes acutely (i.e., the penetration depth is only tens of nanometers). In numerical computations, an acute varying field is harder to be precisely captured than that in a pure, lossless dielectric waveguide.

Lately, many numerical schemes have been developed to analyze the SPP propagation characteristics and further to design various plasmonic integrated optical circuits. Previous published schemes include the finite-difference frequency domain (FDFD) [10,13–16], finite element method (FEM) [12,17–21], finite-difference time domain (FDTD) [8,22–25], method of lines [26–29], differential method [30,31], imaginary-distance beam propagation method [32], and density of states formulation [33]. The widely used FDFD and FEM use many discretization points to achieve accurate results, and thus lead to high memory and computation requirements. The FDTD approach is better suited to pulse propagation problems and dispersion characteristics because it scans a period of the bandwidth. Solving mode problems in a specific frequency is a time-consuming scheme that requires huge memory for three-dimensional structures. The method of lines discretizes all but one dimension to yield systems of coupled ordinary differential equations. These systems can be solved by choosing a time integrator as the time domain approach. The differential method based on the rigorous vector diffraction theories was originally dedicated to analyze periodic structures such as diffraction gratings. The unknown variables in this approach include the transverse electric and magnetic field components. Each field component is expanded by the Fourier basis. The Fourier coefficients are then substituted back to the spatial domain to find the field profiles. In the dielectric-loaded surface plasmon-polariton waveguides (DLSPPWs) [17], the differential method needs 100 terms of Fourier harmonics in each field component to yield the convergent solutions. The transverse spatial discretizations in the imaginary-distance BPM are still based on the finite difference scheme. The density of states formulation, much like the electronic density of state in the presence of a localized electronic defect, computes a Lorentzian shape near a guided wave-vector resonance. The density of states formulation also needs to calculate the Green’s dyad and using multiple-scattering description.

The pseudospectral scheme [34,35] with high-order accuracy and fast convergence was originally developed to solve fluid dynamics problems. In later years, the author combined the scheme with the multidomain approach to study the mode problems of lossless dielectric optical waveguides [36–39]. The calculated results of the multidomain pseudospectral scheme demonstrate that it is indeed reliable and efficient to find accurate effective indices and mode field profiles for various dielectric optical waveguides. Compared with the results calculated by the FDFD and FEM, the multidomain pseudospectral scheme achieves better accuracy and performance in reducing memory capacity because it requires fewer matrix equations [36–39]. This study is the first to use the multidomain pseudospectral scheme to investigate the SPP waveguide structures, including the DLSPPW and inverted metal slot waveguide (IMSW). In previous reports, the multidomain pseudospectral scheme require only 10 terms of basis functions in each subdomain to achieve high accuracy while calculating the lossless dielectric waveguides. In contrast, the plasmonic waveguide structures require extremely fine discretization points to capture the fast decay field within a narrow scale. For the FDFD and FEM approaches, the computational efforts increase significantly while solving the SPP problems. Consequently, this study aims to circumvent these massive computation requirements using the multidomain pseudospectral scheme, and offers an efficient and accurate approach to analyze the mode characteristics of nanoscale plasmonic structures. The remainder of this paper is organized as follows. Section II formulates the wave equations for the three-dimensional waveguide structures. Section III presents the related issues of the numerical scheme. Section IV calculates the numerical examples to demonstrate the accuracy of the proposed scheme, and discusses the mode characteristics of DLSPPW and IMSW in detail. Finally, Section V draws conclusions.

2. Mathematical formulations

Assuming a monochromatic optical wave with a time dependence of exp(iωt) and propagating along the z direction in a medium with refractive index n(x,y,z), the vector wave equation derived from Maxwell’s equations based on the magnetic field vector H is

\nabla^{2} H + k_{0}^{2} n^{2} H + \frac{\nabla n^{2}}{n^{2}} \times (\nabla \times H) = 0

where k ₀ = 2π/λ ₀ and λ ₀ the wavelength in vacuum. For solving a mode problem, the waveguide structure in the longitudinal propagation direction is invariant. The refractive index profile only depends on the transverse directions, that is, n = n(x,y). Here, the wave with z dependence of the form exp(-ißz) is considered, where β = k ₀ n_eff is the propagation constant and n_eff is the mode effective index. The general waveguide structures involve distinct materials with a high-contrast refractive index in the computational domain, especially for the plasmonic waveguides consisting of a dielectric medium and a metal medium with negative permittivity. This approach partitions the computational domain into a few subdomains with a uniform refractive index. The last term of Eq. (1) equals zero due to the homogeneous refractive index profile in each subdomain. This reduces the vector wave equation to two decoupled wave equations for the transverse magnetic field components H_x and H_y as follows:

\frac{\partial^{2} H_{s}}{\partial x^{2}} + \frac{\partial^{2} H_{s}}{\partial y^{2}} + k_{0}^{2} (n^{2} - n_{e f f}^{2}) H_{s} = 0, (s = x, y) .

However, the polarization dependence and the coupling effect can be recovered by incorporating the interfacial boundary conditions. These boundary conditions include the continuous normal and tangential components of magnetic fields at each intra-element boundary. Moreover, the continuities of the longitudinal components H_z and E_z provide the coupling effect of H_x and H_y through the relations

\nabla \times H =

jωε ₀ n ²(x,y)E and

\nabla \cdot H = 0

. For a horizontal interface, since the derivatives of H_x and H_y with respect to x on both sides of the interface are equal, the continuity of E_z yields

n_{y +}^{2} \frac{\partial H_{x}}{\partial y} |_{y -} - n_{y -}^{2} \frac{\partial H_{x}}{\partial y} |_{y +} = (n_{y +}^{2} - n_{y -}^{2}) \frac{\partial H_{y}}{\partial x}

and the continuity of H_z yields

\frac{\partial H_{y}}{\partial y} |_{y +} = \frac{\partial H_{y}}{\partial y} |_{y -}

where y ⁺ and y ^- refer to the locations infinitesimally above and below the horizontal interface, respectively. Likewise, for a vertical interface,

n_{x +}^{2} \frac{\partial H_{y}}{\partial x} |_{x -} - n_{x -}^{2} \frac{\partial H_{y}}{\partial x} |_{x +} = (n_{x +}^{2} - n_{x -}^{2}) \frac{\partial H_{x}}{\partial y}

and

\frac{\partial H_{x}}{\partial x} |_{x +} = \frac{\partial H_{x}}{\partial x} |_{x -}

where x ⁺ and x ^- refer to the locations infinitesimally to the right and to the left of the vertical interface, respectively. Note that the spurious modes do not satisfy the divergence-free magnetic fields [40] and then the spurious modes are completely eliminated in [41] by imposing the constraint

\nabla \cdot Η =0

. Importantly, the constraint

\nabla \cdot Η =0

is explicitly imposed by the proposed scheme to patch the subdomains to effectively prohibit the spurious modes.

3. Numerical approaches

The pseudospectral approach first expands the dependent variables in each subdomain using a set of orthogonal basis functions [34,35]. The product of basis functions and Fourier coefficients are then transferred to the individual Lagrange-type interpolation functions and the grid points of fields. This process transforms the unknown Fourier coefficients in the frequency domain to the field points in the physical domain to obtain the mode profiles directly. The pseudospectral scheme requires that the wave equations must be satisfied exactly at specific collocation points depending on the basis functions. The differential equations are then converted to system of linear equations. For the mode problems, the effective refractive index is found by solving a matrix eigenvalue equation. To preserve the exponential convergence behavior of the pseudospectral scheme, the computational window is partitioned into several subdomains with homogeneous refractive indices. Finally, the subdomains are assembled by patching the physical interfacial boundary conditions to the entire computational domain.

In subdomain r, the transverse magnetic fields (H_x and H_y) are represented as a product of suitable (described in the next subsection) basis function θ(x) in the x-direction, the basis function ψ(y) in the y-direction, and the H_x and H_y field values denoted as $H_{x, p q}^{r}$ and $H_{y, p q}^{r}$ at (n_x + 1) × (n_y + 1) collocation points associated by the corresponding basis functions as follows:

H_{s}^{r} (x, y) = \sum_{p = 0}^{n_{x}} \sum_{q = 0}^{n_{y}} θ_{p}^{r} (x) ψ_{q}^{r} (y) H_{s, p q}^{r}, (s = x, y)

where

θ_{p}^{r} (x_{m}) = δ_{m p}

,

ψ_{p}^{r} (y_{m}) = δ_{m p}

, and δ_mp denotes the Kronecker delta. Note that the number of term of basis function expanded and collocation point is identical. Substituting Eq. (4) into Eq. (2), and then Eq. (2) perfectly satisfies these (n_x-1) × (n_y-1) collocation points in subdomain r. The differential equations are then converted to a matrix eigenvalue equation

[\begin{array}{l} P^{r} \begin{matrix}  \end{matrix} \begin{matrix} \begin{matrix} 0 \end{matrix} \end{matrix} \\ 0 \begin{matrix}  \end{matrix} P^{r} \end{array}] \begin{matrix}  \end{matrix} [\begin{array}{l} H_{x}^{r} \\ H_{y}^{r} \end{array}] = {(k_{0}^{} n_{e f f}^{})}^{2} [\begin{array}{l} H_{x}^{r} \\ H_{y}^{r} \end{array}],

where the operator P ^r is

P^{r} H_{s}^{r} = [\frac{\partial^{2} H_{s}^{r}}{\partial x^{2}} + \frac{\partial^{2} H_{s}^{r}}{\partial y^{2}} + k_{0}^{2} n_{r}^{2} H_{s, p q}^{r}] |_{x = x_{i}, y = y_{j}}, (s = x, y)

= \sum_{i = 1}^{n_{x} - 1} \sum_{j = 1}^{n_{y} - 1} [\sum_{p = 0}^{n_{x}} \sum_{q = 0}^{n_{y}} {θ_{p}^{r (2)} (x) ψ_{q}^{r} (y) + θ_{p}^{r} (x) ψ_{q}^{r (2)} (y) + k_{0}^{2} n_{r}^{2} (x, y) θ_{p}^{r} (x) ψ_{q}^{r} (y)}] |_{x = x_{i}, y = y_{j}} \cdot [H_{s, p q}^{r}]

and

θ_{p}^{(h)} (x)

and

ψ_{q}^{(h)} (y)

indicate the h-th order derivatives of θ_p(x) to x and ψ_q(y) to y, respectively. After obtaining the matrix eigenvalue equation in Eq. (5) for each subdomain, the next step is to assemble all the matrix elements into a global computational domain. If the computational domain is divided into m subdomains, the pattern distribution of the matrix elements forms the following block diagonal matrix equation:

[\begin{matrix} Q_{1} & 0 & 0 & 0 \\ 0 & Q_{2} & 0 & 0 \\ 0 & 0 & ⋱ & 0 \\ 0 & 0 & 0 & Q_{s} \end{matrix}] [\begin{matrix} H_{1} \\ H_{2} \\ ⋮ \\ H_{s} \end{matrix}] = {(k_{0} n_{e f f})}^{2} [\begin{matrix} H_{1} \\ H_{2} \\ ⋮ \\ H_{s} \end{matrix}]

where

Q^{r} = [\begin{matrix} P^{r} & 0 \\ 0 & P^{r} \end{matrix}], Η^{r} = [\begin{matrix} H_{x}^{r} \\ H_{y}^{r} \end{matrix}], (r = 1, 2, 3... m) .

In addition to the elements inside each subdomain, the final matrix must include other contributions at interfacial boundary points (remainder collocation points) satisfying the physical conditions in Eq. (3)a)-(3d) between different materials. Enforcing the boundary conditions from Eq. (3)a)-(3d) recovers the exponential convergence rate of pseudospectral scheme [35]. After adding the interfacial boundary conditions Eq. (3)a)-(3d), the pattern of the final matrix is no longer a block diagonal form, but remains a moderate sparse matrix. For SPP mode problems, n_eff is complex. The real part characterizes the wave vector (which is approximately inversely proportional to the mode size) and the imaginary part represents the propagation length (the power drops to 1/e with respect to the input power) of the SPP mode.

The other issue of the proposed scheme is to determine the basis functions to represent the mode fields. Distinct basis functions depending on the field behaviors are chosen for different subdomains in this work. Chebyshev polynomials are suitable to expand the optical fields in interior subdomains with finite boundaries because of their robustness to non-periodic structures. In contrast, Laguerre-Gaussian functions (LGFs) expand the exponential field profiles of guided modes in exterior subdomains with a semi-infinite boundary. For Chebyshev polynomials, the explicit form of the Lagrange-type interpolation function is as follows [35],

θ_{p} (x) = \frac{{(- 1)}^{p + 1} (1 - x^{2}) T_{v}^{'} (x)}{c_{p} n^{2} (x - x_{p})}, c_{p} = {\begin{cases} 2, if p =0, N \\ 1, if 1 \leq p \leq N -1 \end{cases}

where T_v(x) is the general Chebyshev polynomial of order v and x_p’s are the collocation points for Chebyshev polynomials. For LGFs, the explicit form is as follows [35],

θ_{p} (α x) = \frac{x L_{v} (α x)}{α (x - x_{p}) [x L_{v}^{'} (α x)] |_{x = x_{p}}} e^{- α (x - x_{p}) / 2},

where L_v(αx) is the Laguerre polynomial of order v and x_p’s are the collocation points for LGFs. Parameter α, called the scaling factor, and affects the accuracy for a given term of LGFs. The definite procedure to determine α has been derived in earlier research [37].

4. Simulation results and discussion

This section investigates two plasmonic waveguides. First, the SPP modes in a DLSPPW structure are solved to demonstrate the computational merits of the proposed scheme compared with those of two reliable numerical schemes, the FEM [17] and differential method (DM) [30]. The mode confinement and figure of merit for DLSPPWs have also already been extensively studied in [42] by combining numerical simulations and leakage radiation microscopy. Moreover, the dependences of the mode confinement and propagation length by varying the geometric parameters of DLSPPW are extensively investigated. Secondly, this section studies the mode characteristics of a low-loss IMSW [43] applied to nanoscale circuits in detail, and discusses the figure of merit (FOM) [44].

4.1 The dielectric-loaded surface plasmon-polariton waveguide (DLSPPW)

The mode characteristics of DLSPPW have been verified experimentally by near-field optical microscopy and leakage radiation microscopy, which are able to guide efficiently the SPP modes [45]. A variety of photonic components based on the DLSPPW structure involving ring resonators [46], multimode interference splitter [47], active Bragg reflector [48], and thermo-optics control Mach-Zehnder interferometers [49] have been successively reported because the advantages including high localization, moderate propagation length (tens of micrometers), thermo-, electro-, all-optical functionalities, and easily integrated to optoelectronic devices. For the conventional integrated optics components; however, the efficiency of coupling has not been investigated yet. It would be useful to provide examples of more complex plasmonic circuits. The DLSPPW is constructed by depositing a dielectric polymer ridge with thickness t and width w on a gold film with thickness h covering a semi-infinite dielectric substrate. Figure 1(a) depicts the cross section of the DLSPPW. The refractive indices of materials are air n_a = 1, core n_p = 1.535, gold film n_g = 0.55-11.5i (the form of complex refractive index depends on the plane-wave solution form chosen [50], and this study uses the exp[i(ωt-βz)] solution), and substrate n_d = 1.6 at an excitation wavelength of optical communication λ = 1.55μm. This study investigates the mode characteristics by varying the width and thickness of the polymer ridge at a fixed thickness of gold film h = 0.1μm, and then analyzes the wavelength dependence.

Fig. 1 (a) The cross section of DLSPPW with refractive indices of core n_p, gold film n_g, substrate n_d, and air n_a.(b) the division of computational domain for the DLSPPW structure.

w(μm)
n_x ( = n_y )	0.8	0.7	0.6	0.5	0.4	0.3	0.2
5	1.34830	1.32360	1.29197	1.25117	1.19885	1.13447	1.06509
10	1.34787	1.32287	1.29111	1.25038	1.19851	1.13485	1.06606
15	1.34741	1.32226	1.29036	1.24954	1.19754	1.13395	1.06547
20	1.34742	1.32237	1.29050	1.24974	1.19782	1.13434	1.06582

w(μm)
n_x ( = n_y )	0.8	0.7	0.6	0.5	0.4	0.3	0.2
5	39.3	39.8	40.7	42.5	46.0	53.9	74.4
10	42.2	42.9	43.9	46.0	49.6	58.1	79.4
15	42.9	43.3	44.7	46.2	50.1	58.6	80.3
20	42.2	43.1	44.4	46.1	50.2	58.2	80.2

w(μm)	0.8	0.7	0.6	0.5	0.4	0.3	0.2
FEM	1.348	1.323	1.291	1.250	1.198	1.133	1.064
DM	1.348	1.323	1.291	1.250	1.198	1.133	1.064
This work	1.3474	1.3224	1.2905	1.2497	1.1978	1.1343	1.0658

w(μm)	0.8	0.7	0.6	0.5	0.4	0.3	0.2
FEM	42.2	42.8	44.4	46.4	50.1	59.1	81.3
DM	42.2	42.8	44.0	46.0	50.1	59.1	80.4
This work	42.2	43.1	44.4	46.1	50.2	58.2	80.2

Abstract

1. Introduction

2. Mathematical formulations

3. Numerical approaches

4. Simulation results and discussion

4.1 The dielectric-loaded surface plasmon-polariton waveguide (DLSPPW)

4.2 The inverted metal slot waveguide

5. Conclusion

Acknowledgements

References and links

Cited By

Figures (17)

Tables (4)

Equations (14)

Optics Express