The Forces from Coupled Surface Plasmon Polaritons in Planar Waveguides

: We analytically investigate the forces due to Surface Plasmon Polariton (SPP) modes between ﬁnite and inﬁnitely thick metal slabs separated by an air gap. Using the Drude model and experimentally determined values of the dielectric functions of gold and silver, we study how frequency dispersion and loss in the metals affects the behavior of the SPP modes and the forces generated by them. We calculate the force using the Maxwell Stress Tensor for both the attractive and repulsive modes.


Introduction
Researchers have long held interest in converting electromagnetic energy into mechanical motion. Kepler was the first to hypothesize that solar radiation is responsible for the deflection of comet tails away from the sun. By 1903, Lebedew [1] and Nichols and Hull [2] had proved Maxwell's hypothesis that light impinging on a thin metallic disk in a vacuum would induce measurable motion. Applications for harnessing the energy of light has been seen in system ranging from the "Solar Sail" [3] to optical traps and tweezers [4,5]. Recent work has explored the nature of radiation pressure in high Q-factor microresonators [6,7], negative index systems [8,9], metamaterials [10], photonic crystals [11] and in dispersive dielectrics [9,12,13]. Additionally, studies have explored the evanescent wave bonding and antibonding between parallel dielectric optical waveguides [14,15] and microresonators [16] and the enhancement of radiation pressure in waveguides due to slow-light effects [17].
Surface Plasmon Polaritons (SPPs) offer another avenue for generating mechanical motion from light [18]. SPPs are the result of coherent coupling of photons to free electron oscillations at the boundary between a metal and a dielectric. A significant amount of work has been devoted to studying the coupling of SPPs on surfaces that are in close proximity to one another [19,20]. Long Range Surface Plasmon Polaritons (LRSPPs) [21][22][23], which result from the coupling of SPPs on opposite surfaces of a thin -on the order of the skin depth -metal slab in what is known as the Insulator-Metal-Insulator geometry (IMI), can propagate for distances up to 1 cm when excited at near-infrared frequencies [24]. SPP-induced field enhancement in gaps between metallic nanoparticles [25,26] and between large planar surfaces in the Metal-Insulator-Metal (MIM) geometry [27][28][29] have been used for Surface Enhanced Raman Spectroscopy (SERS) [30][31][32] and the creation of nanoantennas [33].
The forces on metal and dielectric nanoparticles generated by SPP excitation on planar metal surfaces have previously been studied [34][35][36][37][38]. Progress has also been made on the nature of SPP forces in metal nanoparticle clusters [39][40][41][42][43], though to our knowledge, no work thus far has addressed the forces between planar metal surfaces. In this paper we analytically investigate the forces generated by SPPs in the two-dimensional MIM and Insulator-Metal-Insulator-Metal-Insulator (IMIMI) geometries in the cases involving both "lossless" and lossy metals. This paper is structured as follows: in section 2, we derive the expressions for the dispersion of the SPP modes in idealized metal-dielectric systems. We compare the SPP dispersion using the Drude model for the dielectric function of the metal to the SPP dispersion calculated with the tabulated dielectric data of silver and gold from Palik [44]. In section 3, we calculate the forces in the IMIMI geometries, and in section 4, we discuss the characteristics of the force curves the and applications of SPP waveguide forces.

Calculation of the dispersion of SPP waveguides
SPPs are transverse magnetic (TM) polarized modes that exist at the interface of two materials when the real part of the electric permittivity, ε(ω), changes sign across the interface. The most common example of such a system is the boundary between a metal and a dielectric at optical frequencies. The field profile of an SPP at an interface is a solution of the wave equation, permeability of free space, μ 0 , for nonmagnetic materials at optical frequencies, and E(r,t) is the electric field. When ε(ω) changes sign across an interface, the continuity of the normal component of the displacement vector, D(r,t), implies a solution with evanescent fields on both sides of the interface. Using the coordinate system of Fig. 1, we express the electric field as E(r,t) = E 0 exp(ik · r − iωt), where k = k 0 nr = k zẑ + κŷ is the wavevector of SPPs, −ẑ is the direction of propagation and In Eq. (1), k 0 = ω/c, n is the refractive index of the medium, k z is wavevector in the direction of SPP propagation, which is conserved across the interface. We can write k z as β + iα, where β is the propagation constant and α is the loss factor. Im{κ} > Re{κ} for SPPs. For convenience, we define k y ≡ iκ so that for SPPs we can rewrite Eq. (1) as Using these conventions, we can calculate the field profiles for SPPs supported in the two geometries shown in Fig. 1. The subscripts 1 and 2 will always be used to denote the metallic and dielectric regions, respectively in the equations throughout this paper, as labeled in Fig. 1. Fig. 1. The Metal-Insulator-Metal (MIM, (a)) and Insulator-Metal-Insulator-Metal-Insulator (IMIMI, (b)) geometries. ε 1 is the electrical permittivity of the metal and ε 2 is the permittivity of the dielectric. The roman numerals in the IMIMI geometry correspond to the regions defined in Eq. (3). In both geometries, the origin is placed at the center of the dielectric gap of width 2w, and SPP propagation is in the -z-direction in the calculations.
The SPP fields in the IMIMI geometry can be expressed by the following set of equations: where A ...J are the field amplitudes which satisfy the boundary conditions for the fields, k y1 , k y2 are the y-components of the k-vectors in the two materials, and the factor exp [−i(ωt + k z z)] has been dropped from the expressions for clarity. The equations for the MIM geometry are obtained by taking d → ∞.
At interfaces between nonmagnetic materials, E z and H x are continuous. Because of the symmetry of the IMIMI structure, there are two independent solutions which satisfy its boundary conditions: one corresponding to D = F and one to D = −F . We chose to define the overall mode symmetry in terms of the parallel electric field component, E z , which matches the symmetry of the charge distribution in the structure. Thus, D = −F corresponds to symmetric modes and D = F corresponds to antisymmetric modes. Solving the system defined by Eq. (3)-(5) and the boundary conditions for the antisymmetric modes, we find the following relation: When combined with Eq. (2) for each of the two media, Eq. (6) gives a transcendental equation for the dispersion, ω(k z ). The dispersion relation for the symmetric modes is given by replacing tanh(k y2 w) with coth(k y2 w) in the left hand side of Eq. (6). We find that Eq. (6) and its symmetric counterpart each give rise to two fundamental solutions corresponding to a SPP mode. These four IMIMI geometry modes -two symmetric and two antisymmetric -represent the couplings between the four metal-dielectric interfaces in the geometry. To understand these modes, it is helpful to treat the geometry as coupled IMI SPP waveguides, as shown in Fig. 2. The E z field profiles of the two LRSPP modes ( Fig. 2(a)) are antisymmetric with respect to the center of the metal slabs, and can couple together symmetrically ( Fig. 2(b)) and antisymmetrically ( Fig. 2(c)). We refer to these modes as S a and A a , respectively, where the capital character denotes the overall symmetry of the mode and the subscript corresponds to the symmetry of the constituent IMI waveguide modes.
IMI waveguides also support Short Range Surface Plasmon Polaritons (SRSPP), which have shorter propagation lengths due to a larger mode overlap with the metal slabs and have symmetric E z field profiles with respect to the center of the metal slab ( Fig. 2(d)). Two SRSPP waveguide modes will couple symmetrically ( Fig. 2(e)) and antisymmetrically ( Fig. 2(f)). These modes are referred to as S s and A s , respectively.
In the MIM limit (d → ∞), S s and S a are degenerate, so only one symmetric mode exists ( Fig. 2(g)). We refer to it here as S 0 , where the subscript 0 implies this degeneracy. Similarly, the MIM geometry supports only one antisymmetric mode, A 0 ( Fig. 2(h)). At this point, we can use the field symmetry to find the sign of the force generated by each of our modes. Since a symmetric mode corresponds to symmetric charge oscillations, we expect modes with symmetric profiles to generate repulsive forces between the slabs. Likewise, we expect the antisymmetric modes to be attractive.
In order to solve the dispersion relation, we need to model the dielectric function of the meal When these waveguides are brought in proximity to one another, LRSPP 1 and LRSPP 2 will couple symmetrically (b) and antisymmetrically (c). The symmetric Short Range Surface Plasmon Polariton (SRSPP) modes supported by the IMI waveguide (d) will also couple symmetrically (e) and antisymmetrically (f). The MIM geometry supports only two modes, known here as S 0 (g) and A 0 (h). and the insulator. By letting the insulator be air, we can set ε 2 = ε 0 . The simplest model for the metal is the Drude model, which allows us to write the dielectric function as: The Drude model treats a metal as a damped free electron gas where ω p = 2πν p = Ne 2 /ε 0 m 0 is the plasma frequency and γ = Ne 2 /σ 0 m 0 is the damping coefficient. In these expressions, N is the density of free electrons in the metal, e is the electron charge, m 0 is the electron mass and σ 0 is the DC conductivity of the metal. The damping coefficient, γ, is very small compared to ω p for lightly damped systems like noble metals. We find that we can simplify things further by ignoring the loss and taking only the real part of Eq. (7), maintaining the key characteristics of the model and noting that below ω p , ε ε . By substituting the real part of Eq. (7) into Eq. (6), we solve for the dispersion relations of the modes described in Fig. 2, and plot β (ω) in Fig. 3 for both the MIM and the IMIMI geometries, for gap widths of 25 nm and 100 nm. The S 0 mode (cyan lines, Fig. 3(a) and 3(c)) exhibits a cutoff and does not exist at optical frequencies for values of w of interest to us, i.e., w < π/β , where π/β is equal to half of the SPP wavelength. For this reason it will not be discussed in this paper. The A 0 mode wavevector (red lines) increases asymptotically toward a cutoff frequency, ν p / 15 Hz for both gap widths plotted, though it approaches the asymptote more quickly for larger gap widths. Frequency, ν (Hz) The values for silver do not differ from these values enough to produce plots that are distinguishable from those shown here. The thicknesses of the metal slabs in the IMIMI geometry are held constant at 20 nm.
lines) -for gap widths of 2w = 30 nm and 2w = 100 nm, respectively, and a slab thickness of 20 nm in the IMIMI geometry. Comparing Fig. 3(b) to Fig. 3(a) reveals that the IMIMI S a and A s modes have dispersive properties similar to those of the MIM S 0 and A 0 modes, respectively, particularly at small gap widths. For this reason we will also not discuss S a in this paper. A s and the remaining IMIMI modes all approach the ν p / √ 2 asymptote. A a exhibits the least dispersion at low frequencies, as evidenced by the fact that below 10 15 Hz, the wavevector remains close to the light line.
As the gap width increases, we see that the A a and S a modes and the A s and S s modes each approach degeneracy ( Fig. 3(c) & 3(d)). At large gap widths, the interaction of the SPPs between the two slabs weakens, so the remaining two degenerate modes are those of the IMI LRSPP and SRSPP.
The bulk plasmon appears in red in all four panels of Fig. 3 above the light line (β > ω/c) and above the plasma frequency (ν p = 2.18 × 10 15 Hz), where metals experience ultraviolet transparency. Figure 4 shows the dispersion relations for the MIM (Fig. 4(a), 4(c)) and IMIMI ( Fig. 4(b), 4(d)) geometries modeled with the tabulated date for gold from Ref. dielectric function -particularly the imaginary part -of amorphous, polycrystalline and single crystal metals will be notably different from one another, with a variance of up to 20% [45], so it is important to realize that this tabulated data will not precisely match the actual dielectric function of a fabricated metal film. We note that in both of these figures, we have only plotted the A 0 mode in the MIM geometry (red lines, panels (a) and (c)) and the A s (red lines, panels (b) and (d)) and S s (blue lines, panels (b) and (d)) modes in the IMIMI geometry. Once again, the MIM S 0 and the IMIMI S a modes do not exist at optical frequencies for these gap widths and the IMIMI A a mode has such weak dispersion that the force generated by it will be at least an order of magnitude smaller than the A s and S S modes. We have also included the Drude model dispersion (gray dots) for the two modes in both figures for comparison.  Figure 4 shows that the Drude model is an excellent approximation for gold below 4×10 14 Hz (λ 0 ≈ 750 nm), but becomes increasingly worse above that frequency. The reason is that the free electron model for a metal does not account for interband absorption, which begins for gold around the aforementioned frequency, and for silver around 6e14 Hz (λ 0 ≈ 500 nm). When absorption increases to the point that ε (ν) = ε (ν), the SPP mode switches from having normal dispersion to having anomalous dispersion. We refer to this frequency, where dβ /dν = ∞, as the turnaround frequency, ν t , which for gold is approximately ≈ 6×10 14 Hz. For silver, ν t ≈ 9×10 14 Hz.
In both the MIM and IMIMI geometries, the analysis of the modes in Fig. 3 for Drude metals applies to gold and silver. The A 0 wavevector between the slabs is larger below ν t when the gap width is small (Fig. 4(a), 5(a)) than when it is large (Fig. 4(c), 5(c)). In the IMIMI geometry, as the frequency increases toward ν t , the wavevectors of both modes become significantly larger than predicted by the Drude model. However, they still behave in the same way. The S s wavevector at a given frequency decreases as the gap width decreases, while the A s wavevec- tor increases. While gold is more dispersive than silver below ν ≈ 6×10 14 Hz, silver exhibits significant dispersion between ν t,Au and ν t,Ag . Extremely large wavevectors are achievable in small-gap width ( Fig. 4(a), 4(c)) silver-insulator plasmonic structures. In the IMIMI geometry, both A s and S s are extremely dispersive beneath ν t at small and large gap widths, while once again, these modes approach degeneracy at large gap widths (Fig. 4(d), 5(d)).  Figure 6 further illustrates this point. The wavevectors of the three modes are plotted using the Drude model and the Palik data for gold and silver as a function of gap width, β (2w), at a freespace wavelength of λ 0 = 600 nm. The effective mode index, n e f f = β c/ω, is plotted along the right y-axis.
The wavevectors for the MIM A 0 mode ( Fig. 6(a)) calculated using the different models differ only slightly, with the wavevector calculated with Palik's data for gold being predictably larger due to the proximity of the operating frequency to ν t . By contrast, the IMIMI A s wavevector ( Fig. 6(b)) calculated with Palik's gold data is significantly larger than the wavevectors calculated with the Drude mode and Palik's data for silver. Comparing (a) to (b), however, reveals that A s behaves like A 0 , especially when dispersion and loss are low, as is true for silver and Drude metals at λ 0 = 600 nm. The similarity of these modes implies that the A s mode represents a strong coupling between the inner surfaces of the thin metal slabs of the IMIMI geometry, and that the mode's behavior is only weakly dependent on the thickness of the slabs.
Independent of the metal model, the wavevector of these two modes increases exponentially, meaning the group velocity, v g = c(n e f f + ωdn e f f /dω) −1 , decreases exponentially as the gap width between the slabs decreases. Thus these two modes, for extremely small gap widths, can generate slow light, as well as the enhanced field "hot spots" at optical frequencies that has been described in previous MIM waveguide studies [28][29][30][31].
The S s wavevector (Fig. 6(c)) decreases as the gap width decreases for all dielectric models of the metal. This behavior, in contrast to that of the antisymmetric modes, is asymptotic, not exponential. The wavevector, β , is largest between two gold slabs due to the proximity of the operating frequency to ν t . As the gap width approaches zero, the S s wavevector approaches the value of an IMI SRSPP waveguide of thickness 2d, implying that this mode corresponds to a depletion of optical energy from in between the two slabs, in contrast to the enhancement from the antisymmetric modes. We will discuss this further in section 4 when we analyze the energy profiles of the modes. We can also see clearly that at large gap widths -approaching the region of weak coupling between the two gold slabs -the A s and S s wavevectors approach the same value -the value of the SRSPP wavevector in the IMI waveguide geometry.

Calculations of SPP Forces in the MIM and IMIMI Geometries
With the values of the SPP wavevectors obtained with Eqs.
The symmetric solutions can be obtained by replacing cosh(k y2 w) with sinh(k y2 w) in Eqs. (8)- (10). We can relate the field amplitudes for each of the modes to the power flowing in them along the z-axis: where S = (1/2)E × H * is the Poynting vector for complex fields and * denotes the complex conjugate. We can rewrite Eq. (11) as the power per unit waveguide width (see Fig. 1) using the relationship between E y and H x expressed in Eq. (4) as then solve Eq. (12) for |D 2 | in terms of P: whereX ≡ X /D for X = A , B, C . Additionally, k y j ≡ k y j + ik y j , and ε j ≡ ε j + iε j where j = 1, 2, and '±' corresponds to the antisymmetric and symmetric mode solutions, respectively. We will hold P constant at 1 mW/μm throughout this paper.
With the field amplitudes in terms of power, we can solve for the force using the Maxwell Stress Tensor (MST). Starting with microscopic Maxwell's Equations, we can calculate the macroscopic dielectric properties of our system by representing the materials as an ensemble of dipole resonators and taking the average of response. From this, a statement of conservation of momentum can be obtained [46, 47]: where is the MST, EE represents the outer product of the two vectors, P is the polarization vector, with D = ε 0 E + P, ↔ denotes a second rank tensor and ← → I is the identity tensor. The first term on the right hand side of Eq. (14) can be expressed in terms of the momentum of the electromagnetic field, G field , as which is equal to zero when averaged over one period of oscillation. The second term on the right hand side of Eq. (14) represents the mechanical force, and in a sourceless geometry (ρ = 0, J = 0) is written as: where ... denotes the time average. We can see from this equation that the force can be expressed in terms of the local, oscillating charges and currents that result from the polarizability of the material. However, since we do not care about the distribution of the force density throughout our volume, we can use the left hand side of Eq. (14) to find the force in the y-direction, between the metal slabs. We can write the force for the symmetric mode as and the antisymmetric mode forces as where Eq. (19) becomes the negative of Eq. (18) in the lossless limit. Previous work [14,17] showed that one could equivalently calculate the force between dielectric waveguides by taking the spatial gradient of the electromagnetic energy: where U = Nhω and N is the photon density in the mode, and the derivative is taken at constant wavevector, k z , due to translational invariance of the modes. This method is not accurate in plasmonic systems for two reasons. First, translational invariance along the direction of propagation as well as conservation of the adiabatic invariant U/ω, which is proportional to N, cannot be assumed any longer due to optical losses. Secondly, a change in ω and the corresponding change in ε(ω) will lead to a shift in k z , making Eq. (20) nonphysical. Therefore, we must rely on the Stress Tensor to calculate forces. The forces generated by the A s and S s modes in the IMIMI geometry are plotted in Fig. 7(a) and 7(b) for the freespace wavelength λ 0 = 600nm. We plot the force between 20nm thick gold (green lines), silver (blue lines) and Drude metal (red lines) slabs. In Fig. 7(c) and 7(d), we plot the A s and S s mode forces between silver slabs at three freespace wavelengths: λ 0 = 450 nm (cyan lines), λ 0 = 600 nm (blue lines), and λ 0 = 1000 nm (magenta lines). We plot the force in units of pN/μm 2 and note that 1 pN/μm 2 = 1 Pa. We note that the modes in (a) and (b) of this figure correspond directly to the modes plotted in Fig. 6(b) and 6(c). Additionally, we have only plotted the magnitudes of the forces, noting that the S s mode is repulsive and the A s mode is attractive. There are two distinct coupling regimes for the two modes. The first, at large gap widths, is characterized by weak mode coupling and weak forces. The magnitudes of the forces generated by both the A s and S s modes in this regime are identical, as seen by comparing the force curves in Fig. 7(a) to those in Fig. 7(b) and the curves in Fig. 7(c) to those in Fig. 7(d). The gap width at which the forces generated by the A s and S s modes are no longer identical is on the order of 1μm.
The gap width at which the attractive and repulsive modes begin to behave differently depends on the penetration depth of the mode in the dielectric, δ 2 = 1/k y2 , for large w. This value is directly related to the point where the SRSPP modes on the two slabs begin to overlap with each other. At λ 0 = 600nm (Fig. 7(a) and 7(b)), δ 2 is largest for Drude metal slabs and smallest for gold slabs. For silver slabs (Fig. 7(c) and 7(d)), δ 2 is largest at λ 0 = 1000nm and smallest at λ 0 = 450nm. In both of these cases, δ 2 is largest when the operating frequency is closest to the turnaround frequency of the metal, ν t . This agrees with what we would expect by looking at Eq. (2), where we can see that δ 2 should vary inversely with β . We can also see that the force in this regime at a given gap width is stronger when δ 2 is larger.
The second coupling regime is at small gap widths, where the coupling between the two slabs is strong and the attractive and repulsive modes behave quite differently. The force generated by the attractive, A s , mode increases exponentially, but at a slower rate than when the coupling was weak. The strength of the A s mode force is only weakly dependent on the dielectric function ( Fig. 7(a)) and the freespace wavelength ( Fig. 7(a)), with the stronger force occurring when the wavevector, β , is largest.
The force generated by the repulsive, S s , mode peaks at the boundary between weak and strong coupling and decreases as the slabs are brought closer together. At λ 0 = 600nm ( Fig.  7(b)), the force peak is highest for gold and smallest for Drude metals. For silver slabs (Fig.  7(d)), the force peak is highest at λ 0 = 450nm and smallest at λ 0 = 1000nm. The peak is highest when β is largest, which occurs at frequencies closest to ν t for the metal being used.
The strength of the repulsive force in the strong coupling regime corresponds directly to the change that the wavevector undergoes as the gap width changes, as shown in Fig. 6(c). As the gap width 2w → 0, the wavevector asymptotically approaches the value corresponding to a geometry where the two metal slabs are in contact, effectively creating an IMI structure with a metal thickness 2d = 40nm. For silver and Drude metals, the change in wavevector β is small, but it is significantly larger for gold slabs at λ 0 = 600nm.

Discussion and Conclusions
To understand the difference between in behavior of the A s and S s modes more concretely, it is helpful to look at how the distribution of energy changes in the mode as the gap width changes. The electromagnetic energy density in a linear, dispersive material has been thoroughly discussed theoretically [49][50][51] and can be expressed as We can solve this equation in the metal region using the Drude model and in the dielectric region where dispersion is negligible (dε /dω = 0) for the two modes at an operating wavelength (λ 0 = 450 nm) and plot the crossections for a range of gap widths in Fig. 8. We choose Drude metal slabs in Fig. 8 because they most clearly illustrate the key features of the energy distribution within the modes. This analysis applies independent of the metal or frequency, however, as long as it is below ν t . The energy density inside the metal slabs is plotted as having negative value for clarity. Figure 8(a) shows the energy crossections for the A s mode and Fig. 8(b) shows the crossection for the S s mode. Note that the colormaps in the two panels are not of the same scale. At large gap widths, the SPPs on the two metal slabs are essentially uncoupled. The value of the energy density at the inner and outer surface of each metal slabs is approximately equal for both modes, displaying little mode overlap and little interaction between the two modes on the IMI waveguides.
At small separations, the strong coupling across all four metal-dielectric interfaces is evident. In Fig. 8(a), as the gap width decreases below 100 nm, the energy density in the A s becomes concentrated in the space between the slabs, and becomes more than an order of magnitude larger than the energy density outside the slabs. This redistribution of energy, from outside to inside the slabs, is due to the antisymmetric surface charge distribution across the gap, and explains the attractive nature of the A s mode. Additionally, the amount of energy carried in the metal increases at small gap widths, corresponding to the higher n e f f and β seen in Fig. 6(b). Fig. 8. IMIMI energy density crossections at λ 0 = 450nm for geometries using Drude metals. The plots show the energy density of the modes for gap widths between 10 and 400 nm. In (a), the crossections for the A s mode. In (b), the crossections for the S s mode. Note that the colormaps in the two panels are not of the same scale.
Conversely, the energy density of the S s mode Fig. 8(b) decreases to zero as the gap width decreases, while the amount of energy outside of the slabs increases. This redistribution is due to the symmetric surface plasmon charge oscillations across the gap, and corresponds to the repulsive nature of this mode. Furthermore, the energy carried in the metal slabs simultaneously decreases, resulting in the smaller n e f f and β seen in Fig. 6(c).
Nanomechanical forces will play important roles in future devices, both as an avenue for discovery and as a hindrance. SPPs offer an on-chip, optical, solution, to actuation and detection of motion in a nanomechanical resonator, for charge and mass sensing and switching applications. Furthermore, at the length scales relevant to optical forces, one will also have to contend with the Casimir force. For comparison, the Casimir force -−hcπ 2 /3840w 4 -between parallel ideal metal plates separated by 100 nm is 13 pN/μm 2 and decreases to ≈ 5.5pN/μm 2 between two 20 nm thick gold slabs. This value is only slightly smaller than the optical forces presented here at the power level assumed in this paper. We can imagine a system, however, where we can control the power level of our excitation source, and selectively excite the repulsive S s mode to cancel out the attractive Casimir interaction. The ability to generate a net-neutral interaction between supported metallic nanostructures offers new directions for preventing stiction in nanomechanical devices.
We have shown detailed calculations of the dispersion of SPP modes in two geometries: Metal-Insulator-Metal and Insulator-Metal-Insulator-Metal-Insulator. We have treated the metals using the Drude model and with tabulated data for silver and gold from Ref. [44]. From the wavevector dispersion, we have calculated the field profiles, energy, and forces for the modes of these two geometries. Because of the significant dispersion of gold at green-to-red visible frequencies, SPP mode forces are significantly larger than seen with the Drude model and the tabulated data for silver. While the MIM geometry supports attractive forces, the IMIMI geometry support modes with both attractive and repulsive characteristics, making it potentially desirable for many nanomechanical applications.