1. Introduction

Plasmonic waveguides are particularly promising for microelectronic circuits because they have the potential to replace conventional wires, using light as an information carrier with deep subwavelength confinement instead of electrons [1]. The plasmonic waveguiding structures that have been explored to date fall into three major categories: nanoparticle chains [2], metal stripes [3], and the gap waveguides [4–6]. The most important performance parameters of these guiding structures are their propagation length and the degree of confinement they can obtain. The propagation length (L_p) is defined as the length over which the component of the time-average power in the direction of propagation of a mode decays to 1/e of its original value. The confinement is less strictly defined and often depends on the waveguide geometry, though usually it is stated to be the width or area of the region in which the power density is greater than some proportion of its maximum value [7]. The nanoparticle chains have propagation lengths and confinements on the order of hundreds of nanometers [2], and the metal stripes, supporting long-range surface plasmon polariton (LRSPP) modes [8], have propagation lengths of centimeters, but exhibit a confinement of microns [9]. Various geometries of gap waveguides have been explored, such as the v-groove [10,11], and slot or gap waveguides [4]. Gap waveguides provide a balanced tradeoff between confinement and propagation length, having confinements of hundreds of nanometers and a propagation length of several microns [4]. This paper will thus concern the gap plasmon waveguide (GPW), also known as the plasmonic slot waveguide [4].

The steady state propagation characteristics of the GPW have already been thoroughly treated [4,12]. One of the most important devices compatible with this waveguide is the directional coupler. Because this waveguide is not closed from above and below by metal, the field evanescently tails out of the gap vertically, allowing for crosstalk between two such waveguides in this direction.

The vertical coupling structure for the GPW was recently introduced and compared to other directional coupling configurations, namely horizontal couplers [13]. While coupling vertically is not as strong as horizontal coupling, coupling in the vertical dimension is an important degree of freedom for a designer, especially considering that the waveguiding metal layers could simultaneously be used to power electron-based devices embedded between them. In this work we examine in detail methods for inducing strong vertical coupling between gap plasmon waveguides. Additional previous work on vertical couplers has focused on LRSPP waveguides [14].

The rest of this paper is organized as follows. Section two introduces our geometry, restrictions of the parameter space, simulation method, and design methodology. Section three presents simulation results, and finally section four concludes the paper.

2. Design approach

2.1 Geometry specification

Figure 1(a) is a schematic of the geometry of the single GPW. It consists of a rectangular gap etched into a metal film and later backfilled by a dielectric. In this work, the core region may not necessarily be of the same material as the cladding or substrate. The gap has a width w and a height h. Propagation of the wave occurs in the +z-direction, in which the waveguide geometry is considered to be invariant for the purpose of mode solving. The refractive indices of the core, metal, and cladding are denoted via n_core, n_metal, and n_clad, respectively.

Figure 1(b) depicts the vertically-coupled geometry where one GPW is placed directly over the other with a center-to-center (CTC) spacing denoted by s. The refractive index of the middle spacer layer between the two metal films is denoted n_mid. This structure is referred to as a vertical gap plasmon coupler (VGPC). We take the geometry to be vertically-symmetric. To quantify how the total power is distributed between the two waveguides as a function of propagation distance along the coupler, we define the channel 1 region as being that above the line midway between the two gaps, and the channel 2 region as being below it, and thus this line is the boundary between the two semi-infinite waveguide cross-section domains over which the integration of time average power density is done (dashed line in Fig. 1(b)).

 

Fig. 1. Schematic cross-section of the (a) gap plasmon waveguide (GPW) and (b) vertical gap plasmon coupler (VGPC).

Download Full Size | PDF

2.2 Simulation method

For waveguides one can make the assumption that the field dependence in the propagation direction (z-direction here) goes as e ^-γz for forward propagating waves, where γ=α+iβ is the complex propagation constant, α is the attenuation constant and β is the phase constant. Thus the propagation length is L_p=1/(2α). The real part of the effective propagation index (n_eff) is obtained by dividing β by k ₀, where k ₀ is the free space phase constant. The real part of the effective propagation index and the propagation length of the fundamental gap plasmon mode are shown in Fig. 2(a) and Fig. 2(b) calculated via the finite difference frequency domain (FDFD) method [15] and compared with both Reference [12] and finite element frequency domain (FEFD) results using COMSOL Multiphysics™ (for consistency with Ref. 12, Fig. 2 uses 1.55 µm free-space wavelength, w=h=50 nm, n_metal=(-143-10i)^0.5, refractive index of 1.44 for core, cladding, and substrate) [16]. For both FDFD and FEFD formulations, γ is obtained as the eigenvalue of the respective mode. Our FDFD implementation uses a uniform orthogonal grid whose points sit on the geometry boundaries. We enforce homogeneous Dirichlet boundary conditions at the simulation edges and keep them far enough away from the waveguide so as not to have a numerically significant effect. Both the electric (E-FDFD) and magnetic (H-FDFD) field formulations were implemented as a cross-check. For the COMSOL simulations, we use the scattering boundary condition to truncate the boundary, use cubic Lagrangian elements, and keep the maximum element size along the vertical sidewalls of the GPW equal to 0.1h. There is good agreement among the four approaches. The FEFD method was found to converge on a solution to a given accuracy with fewer unknowns and computational effort compared to our FDFD implementation with a uniform grid. Further numerical computation in this paper is accomplished via FEFD.

 

Fig. 2. (a) Effective propagation index obtained by different mode solvers for the fundamental mode of the gap plasmon waveguide with h=w=50 nm, λ ₀=1.55 µm, n_clad=n_core=n_sub=n_silica=1.44, n_metal=(-143-10i)^0.5. b) Corresponding propagation length.

Download Full Size | PDF

2.3 Parameter space reductions

In the analysis that follows, the materials are chosen for their ease of processing and low loss at room temperature. Silver is the metal, with optical parameters from Hagemann, Gudat, and Kunz [17]. At 1.55 µm, n_metal≈(0.434-12.7i). The substrate, cladding and middle are composed of silica with a refractive index of 1.44. The gap cores can have refractive indices of 1.00, 1.44, 2.00, and 2.50. Materials that span this range of refractive indices include air, silicon dioxide, silicon nitride, and Group IV–VI amorphous glasses. The fabrication of air gaps has been feasibly demonstrated using multiple techniques [18,19].

We also restrict the geometry to gaps whose dimensions are on the order of 100 nm. Though it is unphysical, we take the corners of the geometry to be perfectly sharp. Previous work has shown that this simplification still leads to accuracy within 10% for the propagation constants of the supermodes [12,20]. In addition, we have found that giving the corners a radius as large as 5 nm does not significantly affect the analysis.

2.4 Analysis method

The supermode interference (SI) method is used to evaluate the VGPC since our approach is in the steady-state frequency domain [21]. In this method, the power oscillation between the waveguides is expressed in terms of the interference among the superimposed modes of the coupled structure. The interference pattern changes with propagation distance due to the different spatial phase constants of the modes. In general, many modes could be supported by these waveguides depending on the dimensions of the gap and the refractive index of the core; here we restrict ourselves to the case when the gaps are single mode. In this case, the vertically-coupled geometry supports two primary supermodes, both symmetric in the x-direction, and sometimes a third supermode, asymmetric horizontally and vertically. This third supermode is neglected in this analysis on the grounds that it would not be excitable by a horizontally-symmetric source. The two primary supermodes are either symmetric or anti-symmetric in the y-direction. The symmetry is determined based upon the polarization-defining field component, which is E_x for the coordinate system that we are considering. We denote the symmetric supermode with E_s and the asymmetric supermode with E_a. The phase constants for the symmetric and asymmetric supermodes are donated β_s and β_a, respectively.

Of central importance is determining the initial power distribution among the modes given some input condition as this will affect the extinction ratio at the output plane. Taking this into account becomes significant when one of the supermodes is poorly matched to the input excitation (i.e. when the waveguides are close together) and this leads to a decreased power transfer from one waveguide to another. We assume that the input waveguide feeding the coupler is a single GPW whose dimensions and core refractive index are identical to the gaps that comprise the coupler. As shown in Fig. 3, the input waveguide feeds directly into the channel 1 gap. We label the input mode E_i. The input plane of the coupler is our phase reference point.

 

Fig. 3. Side-view schematic of our excitation condition for the VGPC along the yz-plane at x=0. The input mode of the single GPW (E_i) is used to excite the coupler.

Download Full Size | PDF

Figure 4 shows E_i, E_s, and E_a for w=h=50 nm and s=100 nm, with n_core=1.44. The white lines show the boundaries of the metal and the core, the surface plot shows the time average z-propagating power density (S_avgz) profile of the modes, and the arrow plot shows the transverse electric field vector. All of the quantities are normalized to the peak value. From the power density profile, one can see that the wave’s power density peaks at the four corners of each of the gaps, demonstrating that the dominant guiding mechanism in this particular case is due to the coupled “wedge plasmon” mode [20]. The polarization of the modes is shown by the arrow plot, which shows it being dominant in the x-direction. The arrow plot also shows the symmetry of the supermodes: when E_s (Fig. 4(b)) and E_a (Fig. 4(c)) are as shown, they approximately superimpose to E_i (Fig. 4(a)). When the phase difference between them is 180°, the input mode is approximately replicated at the bottom waveguide position, at y=-50 nm.

The length over which this phase shift occurs is the coupling length (L_c=π/|β_s-β_a|) and is half the beat length between E_s and E_a. This length sets the minimum required length of any coupler, and is itself a performance parameter. Given that GPW’s are highly lossy and that their advantage rests with their ability to confine radiation in a subwavelength geometry, it is desirable that this length be as short as possible for compact, efficient devices. We take the effective end of the coupler to be at the plane z=L_c and evaluate performance parameters at this location.

The fundamental mode of the input waveguide determines the initial power distribution of the supermodes by applying transverse field continuity at the input-coupler interface and determining the transmission and reflection coefficients for guided modes (i.e. mode matching method, MMM) [22], [23]. Since we are assuming that the core dimensions and material are the same for the input waveguide as the coupler, we ignore the radiation modes. We then use these transmission coefficients as scalar coefficients weighing the superposition of E_s and E_a in the coupler, calling them t_s and t_a, respectively, normalizing each mode such that P_avgz is equal to one.

 

Fig. 4. (a) The computed transverse electric field vector (arrow plot) and time average power density, S_avgz (surface plot) for the input mode, E_i, of a square silica-filled 50 nm gap, b) the symmetric supermode, E_s, of the corresponding coupler, and c) the asymmetric supermode, E_a.

Download Full Size | PDF

The fundamental mode of the input waveguide determines the initial power distribution of the supermodes by applying transverse field continuity at the input-coupler interface and determining the transmission and reflection coefficients for guided modes (i.e. mode matching method, MMM) [22], [23]. Since we are assuming that the core dimensions and material are the same for the input waveguide as the coupler, we ignore the radiation modes. We then use these transmission coefficients as scalar coefficients weighing the superposition of E_s and E_a in the coupler, calling them t_s and t_a, respectively, normalizing each mode such that P_avgz is equal to one.

To determine the waveguide discontinuity scattering coefficients, we first consider the requirement for tangential continuity of the fields on either side of the input plane:

Here the t superscript denotes the transverse components of the field. The parameter r denotes reflection coefficient. We then apply the unconjugated (applicable for lossy modes) orthonormality condition as the overlap integral using the output modes of the waveguide to reduce (1) and (2) to scalar equations [24]:

The subscripts n and m denote modes of the same waveguide, and δ_nm is the Kronecker delta. Using Eq. (3), Eqs. (1) and (2) are transformed to a set of four equations that can be written in the following matrix form:

The i, a, and s subscripts label the input, asymmetric, and symmetric supermodes, respectively. Since this is an over-determined system, we take the pseudo-inverse of the matrix on the left of (4) to find the scattering coefficients to obtain the best (least squares) approximation [25]. Each ξ_nm is the overlap between the electric field of mode n and the magnetic field of mode m, where the modes here are from different waveguides:

2.5 Performance parameters

The purpose of a directional coupler is to transfer as much propagating power, defined as the time averaged-Poynting power in the z-direction (P_avgz), over from one waveguide to another. One performance parameter of the directional coupler is the extinction ratio (ER), which is the ratio of the output channel power (the time-averaged z-component of the power density, S_avgz, integrated over the channel 2 region) to the input channel (channel 1) power and is ideally infinite so that the desired output channel is isolated from the other. Since we are interested in extracting as much power in channel 2 as possible, we define another performance parameter to be the net coupled output power, P_net, which is the channel 2 power at the output plane, taking into account both the losses at the input of the coupler due to reflected and radiated power, and the absorption losses along the coupler. This intrinsically includes information about ER as well because a large ER coupler will have most of its power in channel 2. We evaluate both P_net and ER at a length along the coupler equal to L_c. We have always found that the ER is maximized at this point, though P_net is actually maximized at a slightly shorter distance. We also examine the guided reflected power (P_refl) and total estimated radiated power (P_rad) in terms of the transmission and reflection coefficients:

2.6 Design strategy

We observe that the choice of the minimum value of the CTC gap spacing, s, requires consideration that the electric field of E_a between the two metal films (to either side of the gaps) points vertically, hence the electric field resembles that of the fundamental mode (TM₀) of the two-dimensional MIM slab. Previous work has shown that the TM₀ mode has no cutoff and that its phase constant increases without bound as the slab width is reduced [26, 27]. E_a can therefore couple to the semi-infinite slab mode for small enough s [13]. Figure 5(a) is a plot of the effective index for the MIM TM₀ slab mode and the E_i, E_s, and E_a effective indices for a coupler with n_core=1.44 and w=h=50 nm obtained via FEFD simulations. Both E_a and TM₀ converge on each other for s values below about 70 nm and E_a is effectively cutoff for s values below 70 nm. At s=80 nm, with the corners of the gaps rounded with a 10 nm radius (weakening the control over the mode that the corners have), and with a simulation domain that was 32 µm wide, we find that the FWHM of S_avgz of E_a along the x-direction at y=0 is only about 2 µm (not cutoff). We do not observe any cutoff for E_s for s values greater than 55nm (corresponding to a metal-metal spacing of 5 nm), even though the transverse electric field off to the sides of the gaps resembles the TM₁ slab. This is consistent with the fact that the TM₁ slab mode is cutoff at this slab width. As expected, the phase constants of E_s and E_a approach that of E_i (the single gap), for large gap spacings. Furthermore, we find that the real part of the n_eff’s of E_s and E_a cross each other at s≈110 nm leading to an infinite coupling length, as can be seen in Fig 5(b). Also, putting the gaps either less than 110 nm apart or at about 190 nm leads to large n_eff differences and short L_c’s.

 

Fig. 5. a) Real part of the effective propagation indices (n_eff’s) for MIM-TM₀, E_s, E_a, and E_i for w=h=50 nm and n_core=1.44 in the target materials-geometry-excitation system note that the effective indices of the supermodes cross each other at s≈110 nm b) coupling length, L_c, versus vertical CTC gap spacing, s.

Download Full Size | PDF

To determine analytically how the power oscillations occur, we examine the definition of S_avgz in Eq. (8), where the “tot” suffix denotes a field quantity that is the weighted sum of the two supermodes:

The average attenuation constant is denoted α′ and β_Δ≡β_s-β_a. The power densities, S_avgz-nm, are given by

The parameters S_avgz-ss and S_avgz-aa are the mode power densities of each supermode by itself. However, S_avgz-sa and S_avgz-as are the cross power densities. The power in each channel is obtained by integrating (8) over the regions specified in Fig. 1(b).

The cross power densities show how power can oscillate as a function of distance along the coupler. The mode power densities are symmetric in the y-direction, whereas the cross power densities are asymmetric. Thus, as the phase of the cross terms is incremented, the cross power densities interfere with the self power densities in a way such that power is subtracted from the top channel and added to the lower one.

The output power of the coupler is related to its efficiency and the extinction ratio between the two channels. The shorter the coupling length, the less attenuation occurs, so the higher the efficiency. The efficiency can also be improved by decreasing the attenuation constant, which can be accomplished by making the gap larger or decreasing the refractive index inside and locally around the gaps. On the other hand, the cross-products of the cross power densities need to be as close as possible to the self-power terms so that power can be transferred effectively (a more complete interference), leading to a large ER. This implies that the transverse field component profiles of the supermodes need to be as similar as possible in shape and relative magnitude (determined by the transmission coefficients). Since the propagation constants tend to follow the shape of the modes, these requirements imply that designing an efficient coupler is an optimization problem where perturbing the mode shapes leads to a smaller ER but shorter L_c and somewhere P_net is maximized.

3. Optimization

The purpose of this section is to optimize the coupler geometry with respect to net coupled output power. Here we explore the effects of changing n_core, the CTC gap spacing, s, and the gap width, w. We vary s from 60 to 500 nm, take w=25, 50, or 100 nm. All other parameters are as specified in section 2.3 above.

3.1 Square gaps

Figure 6 shows the results for the case where w=h=50 nm. The coupling length in general decreases as the waveguides are brought closer together, but there is a singularity (denoted s_inf) where the supermode phase constants cross each other. There is a local minimum in L_c for s values larger than s_inf. An infinite coupling length, however, does not necessarily imply complete isolation. Even when the phase constants are the same, power is transferred by the more rapid attenuation of one supermode versus another (E_a has always been observed to have a higher loss constant for these analyses), in which case, at the limit of large coupling length, almost half of the total power would be transferred. The ER is shown in Fig. 6(b). ER also decreases with s, exhibiting a local minimum near s_inf where a close to unity extinction ratio is obtained. Strictly speaking the ER is not defined at s_inf, so only ER values away from s_inf are considered. Also shown is P_rad in Fig. 6(c), which decreases monotonically as the waveguides are separated further apart, corresponding to a better match of the supermodes to the single input mode. In all cases tested it was found that P_refl was less than 0.5%.

The conditions for maximum net coupled output power is shown in Fig. 6(d). This maximum occurs for s values greater than s_inf. For lower s values, the ER and the modal loss overcome further decreases in the L_c. It is observed that there is a close correspondence between the center-to-center (CTC) spacing that maximizes the net coupled output power, denoted s_opt, and the relative minimum in the L_c, for each n_core tested. This is due to the high loss of these waveguides dominating the output characteristics for the range of s values above s_inf, where the ER is still appreciable, and causing shorter L_c’s to be favored over large ER’s. The P_net is maximized at s_opt=112, 136, 188, and 296 nm for n_core=1.00, 1.44, 2.00, and 2.50, respectively. Figure 6(e) is a plot of the average supermode power propagation length, L_p-avg.

Increasing n_core produces a smaller minimum L_c, larger ER, smaller s_opt and greater P_net. This can be explained by the fact that the increased ER and decreased L_c are able to overcome the decrease in the propagation length. Figure 6(d) shows that the increase in the maximum P_net increases at a decreasing rate with n_core. Because increasing n_core also leads to more highly confined modes, P_rad and P_refl also decrease.

In Fig. 6(f) we validate the supermode interference/mode matching method (SI-MMM) by observing the power oscillations along the length of a coupler using the finite-difference time-domain method (FDTD), as implemented by RSoft [28]. We choose the case where n_core=2.5 and s=s_opt=112 nm. The power in channel 1 and 2 is denoted P_ch1 and P_ch2, respectively. The FDTD simulation consists of an input waveguide feeding into the coupled structure, with the input mode field profile generated by COMSOL as the excitation. Our simulation setup was exactly as depicted in Fig. 3, using the input mode as a CW source. We extracted the power values along the coupler length by projecting the total FDTD field onto the weighted sum of the coupler supermodes obtained from COMSOL, using overlap integrals (similar to those used to obtain Eq. (4) from Eqs. (1) and (2)) to obtain the weights and normalizing the weights to the input mode power at coupler input plane. For this case, |t_s|=0.8024 and |t_a|=0.5687 (both with negligible imaginary parts). Our results show that not only does the SI method given here give the correct propagation length, but also that the MMM, using the pseudo-inverse of the over-determined matrix to find the scattering parameters at the discontinuity, correctly calculates initial power distribution of the supermodes at the coupler input. A benefit of the SI-MMM approach is that it transforms a 3D problem into a 2D one, allowing one to analyze the structure in a matter of minutes instead of days (on a desktop computer).

 

Fig. 6. Performance parameters for the VGPC versus n_core and s: in Figs. (a)–(e), the thick solid line is for n_core=1.00, the thin solid line is for n_core=1.44, the thick dotted line is for n_core=2.00, and the thin dotted line is for n_core=2.50. a) Coupling Length, L_c, which tends to decrease with increasing n_core and decreasing s, but reaches a singularity where the supermode phase constants cross each other b) Extinction Ratio, ER, which increases when either n_core or s is increased c) Radiated power, P_rad, which monotonically decreases with increasing s, d) Net coupled output power, P_net, which is large where L_c is small, and reaches a maximum where s=112 nm for n_core=2.5 e) Average supermode propagation length, L_p-avg, which is a weak function of s, f) Comparison of supermode interference with the mode matching method (SI-MMM) and FDTD in predicting the power oscillations along the coupler for the case where n_core=2.50. Excellent agreement is observed.

Download Full Size | PDF

It is important to note that P_ch2 is actually maximized a little bit before the coupling length. This is due to the high loss of the waveguide causing the total propagating power of the waveguide to decay before the amount transferred to channel 2 has maximized. However, one can see that the difference in powers and coupler lengths is small, so we simply evaluate P_ch2 at L_c. In all cases that we studied, we found that the ER along the length of the coupler was maximized at z=L_c as well, and this could be advantageous for noise reduction.

3.2 Comparison of square gaps to rectangular gaps

In this section, we compare the optimal performance parameters of the square gap with gaps with w=25 and 100 nm (h=50 nm). The relationships of the performance parameters with respect to n_core and s, shown in Fig. 6, mainly carry over to the cases where the gap is no longer square. Because we have found that P_net is maximized at larger values for larger n_core, we summarize the results obtained in Table 1, with n_core=2.50.

Increasing w causes the modes to be less-highly confined in the gaps, and so E_a is more easily perturbed by the proximity of another metal layer. Thus, the s_opt values tend to increase with w. The coupling length also increases, doubling from about 1.5 µm to 3 µm with w increasing from 25 to 100 nm, this being due to the magnitude of the propagation constants being smaller. The ER does decrease by a smaller amount, and with the L_p-avg increasing significantly, this causes P_net to increase by 9 %, even with the longer coupling length. This means that a designer may decide to use a wider gap to increase the output power at the small expense of a decreased extinction ratio. The radiated power stays around 3%, and the reflected power is negligible.

Table 1. Performance parameters for n_core=2.50 for different core widths, w.

View Table

4. Conclusions

In this work, the effects of changing the core refractive index, gap width, and center-to-center spacing on the performance characteristics of vertically-coupled gap plasmon waveguides are examined using the supermode interference technique along with the mode matching method to determine the power distribution among the supermodes at the coupler input. Coupling length, extinction ratio, net coupled output power, radiated power, and reflected power have been identified as coupler performance monitors. Supermode interference in conjunction of mode matching shows excellent agreement with three dimensional finite difference time domain simulation results. It was found that for smaller center-to-center spacings, the phase constants of the supermodes cross each other leading to an infinite coupling length. As the waveguides get closer together, the modal shape mismatch becomes greater as more of the asymmetric mode’s power is conveyed laterally, leading to a decreased ER and resulting in a maximum in the net coupled output power. This maximum corresponds well to the local minimum seen in the coupling length versus center-to-center spacing curves and also occurs at smaller center-to-center spacing values for higher core refractive indices. Making the width of the core larger from 25 to 100 nm makes the modes less confined, increasing the coupling length and the net coupled output power from 1.58 to 3.19 and from 43.7% to 52.0%, respectively, but decreasing the extinction ratio from 7.39 to 6.57. Supermode propagation length can be increased by making the core refractive index low which can be understood in terms of making the magnitude of the complex propagation constant smaller. This structure could be used as a power splitter at half the coupling length. Even though the losses are high, long propagation lengths are not necessary as the gap plasmon waveguide has a high enough confinement that it could be used to connect micron-sized low-loss dielectric waveguides to deep submicron transistor circuits. The tradeoff for confinement is reduced efficiency.