Plasmonic Slot Waveguide Propagation Analysis

Plasmonic slot waveguides provide extreme light confinement with the benefits of having naturally present electrodes for switching and high thermal conductivity of the metal layers to remove excess heat. Past works relied on numerical computation for these structures, which is time-consuming and lacks physical insight. Here, we present an analytical model of plasmonic slot waveguides to determine the modal properties based on single-mode matching to continuum. The model is accurate to within 3% of rigorous numerical simulations. The theory allows for rapid design and provides physical insight into mode propagation in plasmonic slot waveguides for information processing, optical manipulation, and sensing applications.

A careful analysis of the slot waveguide geometry is required to optimize structural parameters for different applications. The trade-off between mode confinement and propagation length is a common challenge. Other design challenges include achieving high sensitivity in sensing applications and optimizing electronic response times (i.e., managing capacitance). Past works used numerical methods to study these geometries and the modes supported by them [1-3, 18, 19, 30]. Although accurate, numerical simulations do not provide physical insight and are resource intensive. For each parameter variation, a new calculation is required to determine the propagation properties. Analytical models compute physical parameters like reflection phase, which can be used in the accurate design and have been applied to other plasmonic geometries [31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48].
We present a fully analytic theory to determine the properties of modes found in plasmonic slot waveguides. An analytic expression of the reflection coefficient at the slotdielectric interface is derived by combining single-mode matching to continuum with an approximate mode shape. The reflection phase values obtained from this equation are used in the transverse resonance waveguiding condition in the geometric optics framework to calculate the mode effective index and propagation length. The results are shown to be accurate by numerical simulations performed with commercially available numerical software. The analytic approach allows for rapid calculation of the propagation properties and provides physical insight into the design space. Compared to a past work that used the effective index method as an analytic approach, here we include a more accurate formulation of the phase of reflection that incorporates the mode-shape mismatch at the boundary [49]. The approach is useful for the efficient design and optimization of plasmonic slot waveguides in applications ranging from information processing to sensing. Figure 1a shows a schematic of the slot waveguide geometry under consideration. Two gold films of thickness l and permittivity m (values taken from a past work [50]) are separated by a distance of a. They are surrounded by dielectric media of permittivity d1 at the top ( z = l∕2 ), d2 at the bottom ( z = −l∕2 ), and d inside the slot. The gap plasmon mode in the slot is quasi-TEM (transverse electromagnetic). The transverse electric field along the x−direction is represented by E x . Propagation takes place in the y − z plane, making an angle with the z−axis. The magnetic field (H), being perpendicular to the propagation direction is also directed in the y (H y ) and z (H z ) directions as shown in Fig. 1c. When exceeds the critical angle ( c ), the gap plasmon mode propagates by total internal reflection at the slot-dielectric interfaces. The two reflection coefficients r 1 and r 2 are shown in Fig. 1b. All subsequent calculations are performed at a wavelength ( ) of 1550 nm.

Effective Index Approximation
Our analysis is based on a combination of perfect electric conductor (PEC) approximation and dielectric loading, which has been previously shown to be accurate for 2D slits [51]. The effectiveness of this approach is explained by two key observations. The hyperbolic cosine dependence of the transverse electric field inside the slot on the position (x) becomes almost constant in the limit a << .
Secondly, the continuity of the normal component of the electric displacement field ( D x ) across each metal-dielectric interface (at x = −a∕2 and x = a∕2 ) gives the ratio of E x in the metal and the slot as d ∕ m . The large values of m at long wavelengths cause E x to be significantly lower in the metal relative to the slot as shown in Fig. 2a. This enables us to approximate the overall electric field by a rectangle function, shown in Fig. 2b, with unit amplitude for |x| < a∕2 and 0 otherwise.
Dielectric loading is introduced as a means of including the finite conductivity and absorption of real metal. The slot is replaced by a material of complex refractive index n MIM , which is the effective index of a gap plasmon mode given by ∕k 0 . k 0 is the free-space wave vector and is the propagation constant in a metal-insulatormetal (MIM) structure calculated by solving the following dispersion relation:

Single-Mode Matching to Continuum
The single-mode matching to continuum method is used to derive an analytic expression for reflection coefficients, r 1 and r 2 . In the subwavelength regime ( a << ), it is a good approximation to assume that only a single mode can be localized in the x−direction inside the slot. The transverse electric and magnetic fields of this mode are matched with the corresponding fields of the continuum of modes outside the slot under PEC approximation.
Considering the interface at the top, the electric field in the slot is given by: Here, rect(x∕a) is the rectangle function whose value is 1 for −a∕2 < x < a∕2 and 0 everywhere else. The electric field just above the slot interface is given as: By matching these fields at the boundary (z = 0 ) and using the orthogonality of modes, we obtain the relationship between r 1 and the transmission coefficient t(u) as follows: where, u = k x ∕k 0 . Similarly, the transverse component of the magnetic field (H y ) inside the slot is given as: and just above the slot is given as: On replacing k x by uk 0 , substituting the value of t(u) obtained from Eq. 4, and applying orthogonality condition to the matched transverse magnetic fields, we obtain r 1 as: where, I 1 is: and w = a . The same steps are used to derive an expression for the reflection coefficient at the bottom of the slot ( r 2 ), the only difference in the resulting equation was that d1 is replaced by d2 . For the case of a real metal, d is replaced by n 2 MIM resulting in new values of r 1 and I 1 as: cos + I 1 and: Figure 3a shows the variation of reflection amplitude with the propagation angle for different slot widths. The material above and below the interface is ( SiO 2 ) of refractive index 1.44 ( d = d1 ). The reflection amplitudes are higher for narrower slots due to increased penetration of the electric fields into the metal. When the angle exceeds the critical angle, given by sin −1 �√ d1 ∕n MIM � and indicated by the vertical dotted lines, total internal reflection takes place but the reflection amplitude is slightly less than 1 due to metal absorption. Narrower gaps have lower critical angles due to a higher difference between refractive indices inside and outside the slot. The reflection phase shown in Fig. 3b increases with and the slot width but reaches for = ∕2 regardless of the gap widths. Figure 3c shows the reflection amplitude for a slot interface with SiO 2 inside and air outside ( d ≠ d1 ). The values of reflection amplitude for corresponding widths are slightly greater compared to those in Fig. 3a. This is due to the increased contrast between the index of the gap plasmon inside the slot and that of air outside as opposed to SiO 2 , which also results in lower values of critical angles. The reflection phases in Fig. 3d are similar in their overall trend to Fig. 3b but they appear more squeezed together due to their lower values also resulting from the higher contrast.

Geometric Optics Approach
The solutions for propagation angles of the plasmonic slot waveguide are calculated from the transverse resonance condition also known as the self-consistency condition. This approach has been previously used to determine the properties of plasmonic stripe waveguide [41]. For a waveguide to sustain a bound mode, the total phase acquired in a roundtrip must be integral multiples of 2 , written as: For all our subsequent calculations we use two configurations of plasmonic slot waveguides: symmetric and asymmetric. For the symmetric case, the two gold films are (12) n m,mode = n MIM sin m completely surrounded by SiO 2 ( d1 = d = d2 ) and therefore r 1 = r 2 . The asymmetric waveguide has air above the top interface and SiO 2 inside the slot and as the substrate Figure 4 shows the total roundtrip phase of the gap plasmon mode in the symmetric configuration for slot widths ranging from 10 nm to 50 nm and film thickness of 200 nm. This is calculated from the left-hand side of Eq. 11 by using 1 = 2 . The horizontal blue solid line corresponds to a phase of 0. Its intersections with the total phase curves, denoted by the red circles are fundamental mode solutions of the symmetric waveguide. Figure 5a shows the variation of the fundamental mode effective index and propagation length with the slot width of the symmetric waveguide. The effective index increases with a decrease in the slot width similar to a 2D slit due to the increasing amount of electric field entering the metal. As a → 0 , the effective index approaches n MIM , since the fundamental order solution → 90 • . These observations are also verified using results from numerical simulations shown by purple squares. Simulations are performed by solving the mode source in ANSYS Lumerical FDTD v2020 R2.3 within perfectly matched layer boundaries. A uniform mesh size of 0.1 nm is used inside the slot and 1 nm is used for the rest of the simulation region. Au (gold) -Johnson and Christy is chosen from the Material Database [50] to form the gold layers and a user-defined index of 1.44 is used for silica. The red solid line in Fig. 5a shows the theoretical values of propagation length calculated by the expression 1∕2[Imag(n mode k 0 ) + ]. Here, accounts for the absorption taking place during reflection computed from the following equation:

Symmetric Waveguide
2 × l × tan is the path covered in a single roundtrip along the mode propagation direction ( y−direction). The red circles in Fig. 5a show the propagation length values gathered from the simulations. Figure 5b shows the variation of the effective index with metal thickness for a slot width of 20 nm. The index values increase with increasing metal thickness and approach n MIM asymptotically as l → ∞ . Conversely, the propagation length of the mode becomes shorter. Figure 6 shows the total phase for an asymmetric waveguide. The slot dimensions are the same as those presented in Fig. 4. The blue solid line represents a total phase of 0 and the red circles represent the modal solutions for the different slot widths. Figure 7a shows the influence of slot width on (13) exp (−2 × l × tan ) = |r 1 r 2 | Fig. 4 Total roundtrip phase of the EM waves in the symmetric waveguide completely surrounded by SiO 2 for slot widths of 10 nm to 50 nm and metal thickness of 200 nm. Horizontal blue line represents a total phase of 0 and the red circles are fundamental mode solutions the fundamental mode index and propagation length of the asymmetric waveguide. Similar to symmetric waveguides, the index increases with decreasing width. However, unlike a symmetric waveguide, it does not always sustain a bound mode for any combination of source and waveguide parameters. Figure 7b shows the variation of mode index with the height of metal film for a gap of 30 nm. The effective index increases with an increase in width and approaches n MIM for l → ∞ similar to the symmetric waveguide.

Asymmetric Waveguide
It is important to note that the theory is most accurate in the fundamental mode calculation when the slot's aspect ratio (ratio of metal thickness to slot width) is greater than 1. Previous work has shown that the gap plasmon mode no longer remains the fundamental mode for lower aspect ratio slots, as a result, the quasi-TEM condition does not hold [26]. The gap plasmon formulation of the plasmonic slot waveguide should therefore be used for l > a.

Higher Order Modes
Although single-mode operation is desired for most applications of plasmonic slot waveguides, we believe it is useful to validate our analysis for higher order modes as well. We use an asymmetric waveguide with the same configuration as in the "Higher Order Modes'' section for a slot width of 20 nm and an increased thickness of 750 nm. As shown in Fig. 8a we obtain zeroth order mode at an angle of 84.8 • , 1 st order mode at 65.61 • , and 2 nd order mode at 37.311 • marked by the red circles. The corresponding effective indices are calculated from Eq. 12 as 2.6308, 2.4060, and 1.6012 shown in Fig. 8b along with the numerical values for the same configuration.
We observe that the presence of higher order modes (in the z− direction) depends not just on the film thickness but also on the slot width. When the width of the slot is changed from 20 nm to 50 nm, keeping the thickness constant at 750 nm, the number of modes reduces from three to two as indicated in Fig. 8a. Their effective indices obtained from theory and numerical calculations are shown in Fig. 8b. The number of modes in the slot (localized in the z− direction) and their corresponding index values increase with an increase in the aspect ratio of the slot.

Discussion
The theoretical analysis of plasmonic slot waveguides has been done by several works using two main approaches: numerical analysis [1-3, 18, 19, 30] and effective index method (EIM) [49]. Though accurate, numerical simulations lack insight into the physics of mode propagation and are time-consuming. In contrast, the proposed analytical method determines mode properties in terms of other physical parameters. We used standard numerical integration of Eq. 10; however, approximate solutions to this equation may be attempted [40] (although we have not found them in this work). Compared to the results of a commercially available finite-difference software, we found that the numerical integration we used was two orders of magnitude faster, while the accuracy was within 3%. The accuracy of the simulations was confirmed by a convergence study.
The main advantage of the presented approach is that it allows for the efficient design of plasmonic slot waveguides while being more accurate than the analytic effective index approaches. This is because the phase of reflection is different when calculated using the EIM which does not incorporate mode-shape mismatch effects. We have shown the difference between the phase of reflection (for a 50 nm slot) obtained from EIM and the proposed analysis in Fig. 9. This is valuable for practical applications where waveguide dispersion, propagation length, and refractive index sensitivity are all more accurately calculated with the analytic approach presented here (as compared to the EIM), while not requiring time-consuming full numerical simulations that lack physical insight. For example, the propagation length is critical for photonic modulators since it provides the on-chip photon loss, and this is found to be significantly more accurate than EIM-based calculations, as shown in Fig. 10. For the dimensions studied here, and assuming that the length of a plasmon-based device is limited by the propagation length, we obtained capacitance × l p × l a values of the order of fF (where is the absolute permittivity of the MIM region and l p is the propagation length of the fundamental mode). This means the characteristic time is of the order of tens to a few hundreds of femtoseconds for 50 Ω resistance so it is unlikely to limit the bandwidth in modulation applications. Information about critical propagation angles is an important feature of the model. Both reflectivity and critical angles are dependent on the refractive index of the dielectric in the slot. The changes in the critical angle can therefore be mapped to changes in the refractive index of the slot, which has been shown to be useful in sensing applications using dielectrics [52]. The strong dependence of the reflection phase on the propagation angles, shown in Figs. 5b and 7b, can also be applied to improve the sensitivity of surface plasmon resonance sensors. We calculated the sensitivity to a 10% and 50% increase in the refractive index of the gap region (silica, 1.44), as shown in Fig. 11. Here, the sensitivity was expressed as the absolute change in the mode effective index (Δn mode ) of the slot in response to the refractive index changes.  The theory can include the gain of an amplified material in the slot providing a way to analyze the slot geometry for lasing applications [53]. Surface roughness can also play a critical role in device performance [54] and there has been recent work considering the impact in ultra-narrow gaps numerically [55]. It is possible that the present approach can be extended to include surface roughness by adding a surface roughness parameter in the future. At present, the approach works only with the rectangular slot geometry; however, past works have generalized the single-mode matching to continuum to other geometries and it is possible that the perfect electric conductor approximation used here may also be extended for rapid calculation in those geometries [46,56]. Although the proposed framework is more accurate than the effective index method and more efficient than numerical simulations, its accuracy sees a decline for extremely narrow gaps or for permittivities where there is significant field penetration into the metal. A larger proportion of the electric field entering the metals reduces the validity of the PEC approximation, limiting the model's accuracy in these scenarios.

Conclusion
In conclusion, we presented a fully analytic approach to determine the propagation properties in a 3D plasmonic slot waveguide. We derived an analytic expression for reflection at the slot's interfaces with surrounding media for different propagation angles by using PEC approximation, dielectric loading, and single-mode matching to continuum. The reflection phase values above the critical angles were substituted in the waveguide self-consistency condition to obtain values of propagation angles that supported fundamental and/or higher order modes. The angular solutions were then used to compute modal properties such as effective index and propagation length for different waveguide parameters. The theoretical results demonstrated close agreement with numerical simulation results, to within 3%.   With a growing focus on the miniaturization of circuits and devices across different areas, our approach will be a useful tool for the rapid design and optimization of plasmonic slot geometries with enhanced physical intuition.