## Abstract

A method for determining the modes that can be guided along infinite chains of metallic nanowires when they are embedded, as in most realistic set-ups, in layered media is presented. The method is based on a rigorous full-wave frequency-domain Source-Model Technique (SMT). The method allows efficient determination of the complex propagation constants and the surface-plasmon type modal fields. Sample results are presented for silver nanowires with circular and triangle-like cross-sections lying in an air-Si-glass layered structure.

© 2011 Optical Society of America

## 1. Introduction

Metallic nano-structures for guidance of surface plasmon-polariton (SPP) type waves at optical frequencies have been the subject of intensive research in recent years. The plasmonic waveguides tightly confine propagating energy in the lateral dimension, defying the diffraction limit *λ*/(2*n*_{core}) of pure dielectric waveguides, hence offering a means of reducing the scale of integrated optical circuits. However, the large field concentration near the metallic structure results in resistive heating losses, therefore yielding a trade-off between the lateral confinement and the propagation length of plasmonic waveguides [1–3]. With prospects that reduced metal content could lower the losses, linear chains of metallic spherical nano-particles were first proposed as plasmonic waveguides with sub-wavelength lateral scale in Ref. [4]. Energy would propagate along the waveguide via inter-particle electromagnetic coupling. Complementary finite-difference time-domain (FDTD) simulations in Refs. [5, 6] suggested that chains of elongated spheroidal particles or nanorods, arranged so that the long axis is perpendicular to the direction of propagation, would be superior to chains of spherical nano-particles in terms of propagation length and confinement. The first experimental evidence of energy transfer along an array of nano-rods confirmed a propagation distance of approximately 500nm [7]. These results spur interest in 2D structures of metallic nanowire (NW) arrays not only as a rudimentary step towards plasmonic wave-guiding in 3D nano-particle chains, but in their own right.

FDTD simulations were the predominant tool for analyses of wave-guiding along chains of NWs, employed also for aperiodic configurations like the funnel feeding [8], or complex structures like the hexagonal NW lattice [9]. However, FDTD calculations yield the fields on a pre-defined grid that is constrained by artificial reflection-less boundaries, which is both computationally intensive and less suited for the analysis of an *open* structure. This invoked research and development of computational techniques more appropriate for analysing these metallic nano-structures. For example, parametric studies of the propagation distance along a single periodic chain of NWs were carried-out with the higher-order discontinuous Galerkin time-domain (DGTD) method [10] and the surface integral equation (SIE) method [11, 12]. The funnel feeding configuration was also optimized with the SIE in Ref. [12].

It is not unknown that in time-domain studies the modal content of the electromagnetic field of a wave-guiding structure is rendered by the specifics of the local external excitation [9, 11]. A more reliable excitation-unbiased characterization of wave-guiding parameters is achieved by pursuing frequency-domain modal analysis, for which the results are often presented as dispersion curves. Modal analyses of guiding structures that are invariant along the direction of propagation, such as the slab waveguide, are performed by a rigorous separation of parameters in a spatial formulation and straightforward calculations of the modal fields. Guiding structures that are periodic along the direction of propagation, such as a chain of NWs, require a more involved analysis with some approximations. In light of this complexity, many efforts have been made to develop techniques for modal analysis of periodic chains of 3D metallic nano-spheroids [13–22]. Chains of 2D NWs have been explored to a much lesser extent [23, 24], in spite of their importance.

The generalized multipole technique (GMT) for modal analysis of metallic NW chains was described in Ref. [23]. The GMT was convenient for analysis of NWs with circular cross-section, whereas Ref. [24] presented a modal analysis method based on the source-model technique (SMT) fit for NWs with arbitrary, smooth cross-section. The SMT [25] is closely related to the GMT [26] in the sense that they are both rigorous, full-wave integral equation methods, yet with an off-surface approach. Other similar methods include the method of auxiliary sources [27], the method of fictitious sources [28, 29], and the method of fundamental solutions [30].

No previous analysis of energy propagation along a chain of metallic NWs accounted for other material interfaces which may lie in proximity to the NWs, although in most cases of practical interest the NW chain would very likely reside on a substrate, or be embedded beneath a cover layer. The influence of an air-glass interface on the dispersion curve of a nano-particle array was demonstrated in Refs. [20, 21], and its effect was clearly not negligible. In this paper we build on the methodology presented in Ref. [24], and proceed to present a method for modal analysis of optical wave guiding along periodic chains of wires of general, smooth cross-section, embedded in multi-layered media, based on the SMT. This paper is organized as follows: The general structure we considered is outlined in Sec. 2 and the problem-specific method of analysis is given in detail in Sec. 3. Following the analytical work, we developed a corresponding computational tool which is robust to the choice of specific parameters of the geometry. We studied a few representative problems of NWs in layered media, and sample numerical results are presented in Sec. 4. The results are summarized in Sec. 5.

## 2. Problem description

We consider the following geometry: a linear chain of metallic cylinders of arbitrary cross-section, characterized by a plasma-type permittivity *ɛ*_{c}, and oriented parallel to the *z* direction, as shown in Fig. 1. The linear chain is periodic in the *x* direction and the period is denoted by *L*. The material surrounding the cylinders is a dielectric slab characterized by permittivity *ɛ*_{mid}, and free-space permeability, *μ*_{0}. The slab resides on a dielectric half-space at *y*_{sub}. The dielectric lower-half-space serves as a substrate, and is characterised by *ɛ*_{sub}, *μ*_{0}. The layer atop *y*_{top} is a dielectric half-space characterised by *ɛ*_{top}, *μ*_{0}. The triple layered geometry as shown in Fig. 1 is an example of a NW chain waveguide in a multilayered structure, without loss of generality.

For all fields and sources, an exp(*jωt*) harmonic time variation is assumed and suppressed. For a given real *ω*, our aim is to determine proper modes of this structure, i.e., source-free EM fields which can propagate along the periodic structure with a constant wave-vector component *k _{x}*. The search for

*k*is in the complex plane, to account for complex phase accumulation along the direction of propagation in a structure with losses. As the modes of interest are proper, the lateral wave-vector components are required to obey the radiation condition. This means that for

_{x}*y*→ ±∞, the fields must be either zero or composed entirely of outgoing waves. A more formal mathematical statement of the radiation condition for periodic structures is somewhat involved (see Ref. [31, p. 54]).

In the following section we present an analysis only for TE to *z* modes, which have a *z*-directed magnetic field, as these are generally not attenuated as much as TM to *z* modes, which have a *z*-directed electric field. The electric field of the latter modes is tangential to the metallic interface and is virtually shorted-out. The extension to TM to *z* modes is trivial and has been integrated in the computation tool, but is not presented here.

## 3. Source-Model Technique (SMT) problem-specific solution scheme

For a periodic structure with a period *L* in the *x* direction as shown in Fig. 1, it is common practice to reduce the analysis to consider a single unit cell while enforcing periodic boundary conditions

*F*(

*x*) represents any of the EM field components, and

*k*is the propagation constant. The EM field components are usually conveniently represented in terms of an infinite sum of Floquet modes with propagation constants For a periodic structure with losses, the same periodic boundary condition is practiced, yet

_{x}*k*is allowed to take complex values to account for both phase shift and amplitude decay between the opposite walls of the unit cell. Then, an infinite sum of Floquet modes of respective complex valued

_{x}*k*

_{xn}can represent the fields of a lossy periodic structure in a concise analytical form.

Previous analyses of periodic 2D structures based on the SMT have successfully reduced the analysis to consider a single unit cell, either for solving a scattering problem [32–34], or for source-free modal analysis [24]. Once the unit cell was defined, the next step of the SMT algorithm was to separate the structure enclosed in a unit cell into subregions of homogeneous media according to the media boundaries. However, in the following analysis of a unit cell featuring a NW in a triple-layered media, the subregions are chosen slightly differently than the general guideline. We have found it possible to analyze the structure efficiently by separating the unit cell according to the NW boundary only. This yields just two subregions: inside the NW, and the rest of the unit cell.

A set of sources is arranged on a mathematical surface “outside” each subregion for the sake of approximating the fields inside that material. The fictitious sources are distributed uniformly on curves which are slightly dilated (or contracted) versions of the NW boundary curve, as depicted in Fig. 2. A simple method for generating such curves for arbitrary smooth boundary curves is described in Ref. [35]. To avoid overlap of fictitious sources with finite-volume, as their number *N* is increased the curves along which the sources are placed should approach the boundary.

Once located, the choice of radiating sources, their current distributions, or their current flow directions is greatly dependent on the fields that they are meant to approximate. For example, to approximate radiating fields with polarization TE to *z*, which is the focus of this work, it was recommended in Ref. [36] to choose magnetic-current line sources oriented parallel to the *z* axis. Following these guidelines, the fields *inside* the NW are approximated by the fields due to magnetic current filaments of yet to be determined current amplitudes, distributed on a curve which is a slightly dilated version of the cylinder’s perimeter. The filaments are regarded as radiating in a homogeneous, unbounded, obstacle-free medium that has the same material parameters as the cylinder, as shown in Fig. 2. Hence, the radiating fields are calculated from the trivial 2D Green’s function

*k*

_{0}being the free-space wave number. In Eq. (3), (

*x*

_{0},

*y*

_{0}) and (

*x,y*) denote the filament coordinates and testing point coordinates, respectively.

The fields within the unit cell and *outside* of the NW are approximated by the fields due to current sources operating in the unit cell with the cylinder removed and all other material interfaces intact. The location of the fictitious sources in the modified unit cell is shown in Fig. 3.

The fields of these fictitious sources are subject to a number of boundary conditions: the periodicity condition of a lossy structure, multiple reflections from the planar interfaces of nearby media layers, and the proper radiation condition at infinity. We follow the footsteps of Refs. [24, 32–34] for the choice of appropriate sources, and employ fictitious current *strips* rather than current filaments. The fictitious strips are assembled with their expanded dimension in the direction of the periodicity, so they are “less singular” along this direction. Preferring strip sources over spatially-impulsive sources greatly improves the convergence rate of the infinite sum of Floquet modes representing the fields, without the need to resort to more complicated methods like the Ewald summation approximation [37], or various other methods described in Ref. [38]. Moreover, it was suggested in Ref. [39] that the strips could support a non-trivial magnetic current distribution, that if chosen wisely would further accelerate the evaluation of the infinite Floquet sum. The Blackman-Harris window [40] was found adequate for this purpose, because it is a smooth, real-valued function with non-zero values only in a finite interval corresponding to the width of the strip, and rapidly decaying Fourier coefficients that contribute to the speedy convergence of the field representation.

The next step after choosing properly modulated fictitious sources and source locations is to adjust their fields so that all boundary conditions within the unit cell are fulfilled. At this end, we note that in comparison with the multilayered unit cell presented in Fig. 3, the structure analyzed in Ref. [24] can be considered as a “degenerate multilayered unit cell”, consisting of a single homogeneous cladding layer and no other material boundaries. In consequence, the results of Ref. [24] can be referred to as the solution of a rudimentary private case of our study. So, we recall the formalism of the proper TM to *z* fields radiating from a magnetic-current strip source that is encapsulated in a unit cell composed of a single homogeneous material of relative permittivity *ɛ*_{mid}, and subject to lossy periodic boundary conditions. The *z* directed component of the magnetic field is given by

*k*

_{xn}is taken from Eq. (2) with a complex

*k*, to introduce losses. The Re{

_{x}*k*} can be restricted to the First Brillouin Zone (FBZ), i.e., between −

_{x}*π*/

*L*and

*π*/

*L. k*

_{yn,mid}is given by

*y*→ ±∞, which leads to Im(

*k*

_{yn, mid}) < 0. Also in Eq. (4),

*f̂*are the Fourier coefficients of the strip’s Blackman-Harris magnetic-current distribution function

_{n}*f*(

*x*), calculated via

The formulation of the radiating fields in Eq. (4) as an infinite sum of Floquet modes proved itself consequential for two reasons. First, as pointed-out in Ref. [24], this spectral analysis avoids the divergence problem that is inherent to any spatial formulation of fields radiating from an infinite periodic recurrence of sources in a structure with losses [14]. Second, it appears convenient to expand this solution of the rudimentary problem so it solves the general problem, in which the periodic strips radiate in presence of planar interfaces. The effect of the nearby planar interfaces can be straightforwardly addressed by treating each Floquet mode as a simple plane wave, and considering its multiple refractions at each interface.

We consider how the *n*^{th} plane wave in the infinite sum in Eq. (4) undergoes refraction at the interface with the substrate medium, which is characterized by permittivity *ɛ*_{mid}. Let us define the following generalizes reflection and transmission coefficients:

*k*

_{yn,mid}is taken from Eq. (6) and ${k}_{{\text{y}}_{\text{n}},\text{sub}}={\left({\varepsilon}_{\text{sub}}{k}_{0}^{2}-{k}_{{\text{x}}_{\text{n}}}^{2}\right)}^{1/2}$,

*y*

_{0}is the srip’s lateral coordinate, and

*y*

_{sub}is the location of the respective interface. The reflection and transmission coefficients between the middle and top layers, ${\Gamma}_{\text{mid}-\text{top}}^{{H}_{n}}$ and ${T}_{\text{mid}-\text{top}}^{{H}_{n}}$, are also calculated based on Eq. (8) by substituting

*ɛ*

_{sub},

*k*

_{yn,sub},

*y*

_{sub}

*ɛ*

_{top},

*k*

_{yn,sub},

*y*

_{top}, respectively. Clearly, the reflection and transmission coefficients for the opposite polarization are slightly different.

To express the fields due to the magnetic current strips in the unit cell with two cladding layers present, we have to take into account multiple internal reflections in the middle layer. Conveniently, our generalized reflection and transmission coefficients in Eq. (8a) include the phase accumulation in the lateral direction, so we may write the well-known expansion of multiple internal reflections in a finite slab simply as

*z*-directed component of the magnetic field is compactly written as

The analysis of the triple layered geometry, as shown in Fig. 3 and examined above, is applicable for a unit cell with a general number of material layers. The “degenrate” cases of a homogeneous unit cell or double layered unit cell can be easily derived from Eq. (10), and indeed for the homogeneous unit cell one can verify that Eq. (10) reduces to Eq. (4). For a unit cell with more than three material layers, the presented analysis is conceptually equivalent to analyzing the fields in the middle layer and the two most distant cladding layers of the multilayered structure. The expansion of the field expressions in the intermediate layers, for example via transmission matrices, is trivial and is therefore not presented in our analysis.

After defining fictitious source locations and current distributions in both sub-regions of the original unit cell, the continuity of the tangential field components at the cylinder boundary is enforced at a discrete set of testing points. The testing points are uniformly distributed on the cylinder perimeter. This leads to a homogeneous matrix equation

where [*Z*] is an 2

*M*× 2

*N*impedance matrix and

*K⃗*is a column vector of unknown current amplitudes. Here, the number of testing points is

*M*, and the continuity of both the tangential electric and magnetic fields leads to 2

*M*equations. The number of current sources is 2

*N*, of which there is an equal number inside and outside of the cylinder. The structure of the impedance matrix is

*η*

_{0}) at the testing points. It is helpful to use more testing points than sources, since this spreads the error in the continuity conditions more uniformly along the boundary. It is also useful for avoiding certain spurious solutions [41].

Proper modes of the structure are recognized as pairs of *k _{x}* and

*ω*for which a non-trivial solution to Eq. (11) exists. False solutions, which are inherent to most integral-equation formulations and to an even larger extent in the SMT [42], yield a difficulty in identifying singular [

*Z*] matrices. Therefore, we avoid examining the determinant (in case the matrix is square) or the condition number (in the general case), and prefer to employ the spurious-free mode determination method of Ref. [24]. In short, the following alternative generalized eigenvalue problem:

*k*,

_{x}*ω*pairs. In Eq. (13), the dagger sign indicates conjugate-transpose. The square-root of the generalized eigenvalue, Δ

*E*(

*K⃗*), gives a physical implication as to the validity of a solution of the original eigenvalue problem in Eq. (11). In essence, Δ

*E*(

*K⃗*) is a normalized error measure for the eigen-solution

*K⃗*; it is the ratio between the absolute

*difference*of the tangential components of the fields across the cylinder boundary, vs. the absolute

*value*of the tangential components of the magnetic field at the testing points. Thus, seeking the minimal (Δ

*E*)

^{2}which solves Eq. (13) is equivalent to finding a “true” modal solution of the structure.

## 4. Numerical results

To demonstrate the capabilities of our solution scheme, a few numerical results are given. First, we consider a linear chain of NWs of circular cross-section in a triple layered structure, as illustrated in the inset of Fig. 4(a). In this and all numerical examples, the NWs are assumed to be made of silver, which is characterized by *ɛ*_{c} taken from Johnson & Christy’s experimental data [43] rather than the approximate Drude model. The radius of the NWs is taken as *R* = 20nm, and the period of the chain is *L* = 50nm. The NW chain is assumed to reside in a Silicon slab that is characterized with complex permittivity values taken from Ref. [44]. The Silicon is propped on a half-space glass substrate with *n*_{Glass} = 1.51 at *y*_{sub} = −5*R*, and covered by air. The height of the air / Si interface, *y*_{top}, is assigned varying values.

For modal analysis of the structure, the normalized frequency *ωL*/(2*πc*) is taken as the independent variable, for which complex *k _{x}* values within the FBZ are found. The values of

*ωL*/(2

*πc*) are evenly spaced. The appropriate

*k*values are sought to a tolerance of 10

_{x}^{−7}

*k*

_{0}, so that the normalized absolute error as defined in Eq. (13) is smaller than 0.02. The results are presented in the form of dispersion curves, hence the normalized

*ω*(and

*λ*) values are plot vs. their Re{

*k*} and Im{

_{x}*k*} values, in Fig. 4(a) and Fig. 4(b), respectively. This is a common form of displaying results of modal analyses, although somewhat “in reverse” to the scheme of solution used here. For reference, an identical NW chain embedded in a

_{x}*single*, homogeneous Silicon layer is also examined, and the appropriate dispersion curve is shown in Fig. 4 in blue diamonds.

The analysis is limited to plasmonic-type modes, which propagate along the metallic NW chain as evanescent surface waves. As these solutions exist outside of the light cone, only modes within the FBZ that have a propagation constant Re{*k _{x}*} > max

_{i}_{=top,mid,sub}Re{

*n*}

_{i}*k*

_{0}are inspected, as evident in Fig. 4(a). When Re{

*k*} is close to the light-line in Si, which implies that the mode propagates largely in the Silicon surrounding the NWs, the attenuation is small. As Re{

_{x}*k*} increases towards the edge of the FBZ at Re{

_{x}*k*} =

_{x}*π*/

*L*, the mode is indeed more surface-plasmon like. Specifically, its group velocity diminishes, the mode becomes more con-fined to the metallic NWs, and thus the attenuation increases rapidly. The results are presented in Fig. 4 only for modes with an effective propagation length L

_{−3dB}longer than 150nm. Modes with stronger attenuation were found, but are not shown here. For a NW chain in a single Silicon layer or for a triple-layered structure with

*y*

_{top}= 500

*R*more than one eigen-solution was found to exist at frequencies around

*ωL*/(2

*πc*) = 0.058. This highlights the efficiency of this computation tool at separating degenerate plasmonic-type modes, which is not guaranteed in other methods.

The multilayer structure is not symmetric, so a cut-off of the fundamental TE mode is observed in Fig. 4. This phenomenon is also well known in asymmetric dielectric slab waveguides. Furthermore, as we let either interface (in this example, the air-Si interface) approach the linear chain of particles, the dispersion curves and hence the cutoff frequency are found to be blue-shifted. This can be explained intuitively by considering the reverse scenario: as the interface is drawn away from the NWs, its’ interaction with the modal fields lessens and the mode experiences “pseudo-symmetric” boundary conditions. Thus, for larger *y*_{top} values the dispersion curves in Fig. 4 approach those of a symmetric structure, of a NW chain in a single Si layer, and specifically the cut-off frequency of the fundamental mode is red-shifted to coincide with the zero cut-off frequency of the curve of blue diamonds. The same line of reasoning can explain why for large values of Re{*k _{x}*}, which represent strongly confined modes, the dispersion curves coincide regardless of the value of

*y*

_{top}. Also, the dispersion curves of higher-order modes reach large Re{

*k*} and near-zero group velocity values rather quickly, so they are less affected by the distance of the NWs from the air-Si interface. Indeed, in Fig. 4, the dispersion curves of the 2

_{x}^{nd}and 3

^{rd}order modes in a structure with

*y*

_{top}= 500

*R*down to 5

*R*are practically indistinguishable.

The (normalized) transverse magnetic field Re{*H _{z}*} of the fundamental propagating plasmonic mode in the aforementioned structure is presented in Fig. 5 for a fixed frequency

*ωL*/(2

*πc*) = 0.05355 as the parameter

*y*

_{top}varies. The parameter

*y*

_{top}equals 2

*R*, 3

*R*and 5

*R*in Figs. 5(a), (b) and (c), respectively. The three mode profiles correspond to the three points marked by squares in Fig. 4. In the SMT, the solutions obey Maxwell’s equations exactly and the continuity conditions approximately. Therefore, the continuity of Re{

*H*}, which graphically appears perfect in Fig. 5, attests to the validity of the solutions. Also as expected, the field is not periodic with

_{z}*L*but rather undergoes a

*k*phase shift. The imaginary part of

_{x}L*k*, which represents the attenuation of the modal fields, is not apparent in these short segments of the chain but can be observed in a longer segment. Figure 5 also demonstrates the dynamics of the fundamental mode-profile evolution as the air-Si interface is drawn away from the NW chain. As could be anticipated from the asymmetry of the multilayered structure, the fundamental mode is neither odd nor even with respect to

_{x}L*y*, though the mode profile seems more symmetric for a more distant air-Si interface.

A slightly more complex geometry we consider is the structure shown in the inset in Fig. 6(a). The perimeter of the triangle-like cylinders in this structure is represented in terms of the azimuthal angle *ϕ* by:

*R*is the radius of a circle that tightly bounds the curve and

_{b}*ν*≥ 1 is a parameter which determines the curvature of the rounded corners of the triangle-shaped curve. In Fig. 6 we compare the dispersion curves of chains of NWs bounded by hypotrochoid curves with

*ν*= 2 and different

*R*values, against a reference dispersion curve of a chain of circular NWs with

_{b}*R*= 20nm. For the structures under comparison all other geometric parameters are set as follows: the chain periodicity is

*L*= 50nm, the silver NWs are found in a Si slab, lying on a glass substrate at

*y*

_{sub}= −5

*R*, and covered by air from

*y*

_{top}= 3

*R*. The purple dispersion curve corresponds to a chain of triangle-like NWs with the same perimeter as the circular NWs, and the brown dispersion curve corresponds to a chain of triangle-like NWs with the same area (per 2D unit cell) as the circular NWs. Neither the triangle-like NWs with the same perimeter nor those with the same area fit the dispersion curve of the circular NWs at all frequencies. This signifies the importance of the shape of the NW’s cross-section, even when other geometric parameters are identical.

Although a layered structure with triangle-like NWs cannot be modeled with any of the approximate coupled-dipole / quadrupole solution schemes for nanowire chains, with our method it is solved straightforwardly and with no more complexity than the circular nanowire case. Also, as the SMT is a full-wave technique, it is straightforward to plot the normalized mode profiles of the modes which propagate along the structure with NWs of any cross-section. Figure 7 shows the magnitude of the magnetic field |**H**| in color and the power-density flow Re{**S**} as black arrows for the three pairs of points marked by *circles* and *triangles* in Fig. 6, corresponding to a chain of circular or triangle-like NWs of equal perimeter, respectively. For improved visibility the arrows are scaled by a non-linear transform, hence their sizes are relative, not absolute. Therefore, it is hard to notice the attenuation in this figure. It is evident in Fig. 7 that the time-average power density inside the cylinders flows opposite to the direction of the modal front, hence in direction −*x*. This is expected, because turbulence of the power flow density in the vicinity of plasmonic metallic interfaces is a well-known phenomenon [45].

The vertically adjacent sub-figures in Fig. 7 present the “same” mode in circular and triangle-like NWs for identical frequencies, material and geometric properties. Figures 7(a),(d) display the fundamental mode, which is recognized at the quasi-static limit as a coupled dipolar mode. Figures 7(b),(e) display the 2^{nd} order mode, which is analogous at the quasi-static limit to the quadrupolar mode. Figures 7(c),(f) display the 3^{rd} order mode, in which the field magnitude looks to have a number of local extrema along the circumference of each element in the chain, and exhibits strong destructive interference patterns in each unit-cell with a strong dependency on the shape of the cross-section of the NWs. The interference patterns seem relate to a geometric resonance along the circumference of each NW, thus implying an analogy to the quasi-static limit of particle multipolar modes is less appropriate for 3^{rd} and higher-order modes. We wish to stress that the higher-order modes, which have been explored to a lesser extent by other methods, can be easily identified and plotted with the SMT analysis tool.

To give a sense of the convergence of the numerical procedure, we calculate the propagation constant along chains of circular and triangle-like NWs at a frequency matching that of Figs. 7(b),(e) with an increasing number of sources, *N*. The results for a Matlab^{®} 2010b implementation of the proposed method, running on a 32-bit instruction set on an Intel^{®} Core 2 Duo CPUs @3GHz each, are shown in Table. 1. Here, the propagation constant is detailed to an accuracy of 6 significant digits, unlike Figs. 7(b),(e), which present a rounded value with only 2 significant digits. For this calculation, the summation of Floquet modes in Eq. (10) is truncated when the partial sum is at least 10^{5} times greater (in absolute value) than the difference between consecutive terms. *T* is the computation time in seconds for identifying the complex *k _{x}* which fits the chosen frequency at the desired level of precision, and Δ

*E*is the normalized error in Eq. (13).

As expected, for an increase in *N* the computation time increases and Δ*E* is seen to decrease. It is also evident that the run-time is independent of the cross-section of the NW, while the accuracy of the modal solution for triangle-like NWs is (at worst) an order of magnitude less than that of the circular NWs. Taking into consideration the rapid field variations at the sharp bends of the triangle-like shape, the average error at evenly-spaced testing points is indeed expected to assume a higher value than for the circular NW case.

## 5. Summary

This paper expands a full-wave SMT-based modal analysis of metallic NW chains in unbounded free-space [24] to encompass the case of metallic NW chains embedded in layered media. The SMT is rigorous and versatile, and can be used as frequency-domain method, thereby rendering it more suitable than time-domain methods for modeling frequency dependent phenomena like modal solutions. During the course of validating the solution, a few representative cases of NW chains in triple-layered media, such as nanowires of circular and triangle-like cross sections were studied. The solution scheme was found to be robust to the choice of specific parameters of the geometry, and non-circular NWs were treated with very much the same efficiency as circular ones. Standard plasmon-type modes, such as the dipolar and quadrupolar modes, were identified as easily as higher-order modes. The numerical results demonstrated that the shape of the NW cross-section clearly affects the dispersion curves of the NW chain. The proximity of the NWs to a planar material interface was shown to alter the symmetry of the mode-profile, and to blue-shift the cut-off frequency of the plasmonic-type modes.

## Acknowledgments

The authors wish to thank A. Hochman for valuable discussions and for providing supplementary code sections which have not yet been included in the online version of the SMT-Package [42]. This work was supported in part by the Israel Science Foundation and by the Goldberg Fund For Electronics Research.

## References and links

**1. **R. Zia, M. D. Selker, P. B. Catrysse, and M. L. Brongersma, “Geometries and materials for subwavelength surface plasmon modes,” J. Opt. Soc. Am. A **21**, 2442–2446 (2004). [CrossRef]

**2. **P. Berini, “Figures of merit for surface plasmon waveguides,” Opt. Express **14**, 13030–13042 (2006). [CrossRef] [PubMed]

**3. **R. Buckley and P. Berini, “Figures of merit for 2d surface plasmon waveguides and application to metal stripes,” Opt. Express **15**, 12174–12182 (2007). [CrossRef] [PubMed]

**4. **M. Quinten, A. Leitner, J. R. Krenn, and F. R. Aussenegg, “Electromagnetic energy transport via linear chains of silver nanoparticles,” Opt. Lett. **23**, 1331–1333 (1998). [CrossRef]

**5. **S. A. Maier, P. G. Kik, and H. A. Atwater, “Observation of coupled plasmon-polariton modes in au nanoparticle chain waveguides of different lengths: Estimation of waveguide loss,” Appl. Phys. Lett. **81**, 1714–1716 (2002). [CrossRef]

**6. **S. A. Maier, P. G. Kik, and H. A. Atwater, “Optical pulse propagation in metal nanoparticle chain waveguides,” Phys. Rev. B **67**, 205402 (2003). [CrossRef]

**7. **S. A. Maier, P. G. Kik, H. A. Atwater, S. Meltzer, E. Harel, B. E. Koel, and A. A. Requicha, “Local detection of electromagnetic energy transport below the diffraction limit in metal nanoparticle plasmon waveguides,” Nature Materials **2**, 229–232 (2003). [CrossRef] [PubMed]

**8. **S. K. Gray and T. Kupka, “Propagation of light in metallic nanowire arrays: Finite-Difference Time-Domain studies of silver cylinders,” Phys. Rev. B **68**, 045415 (2003). [CrossRef]

**9. **W. Saj, “FDTD simulations of 2D plasmon waveguide on silver nanorods in hexagonal lattice,” Opt. Express **13**, 4818–4827 (2005). [CrossRef] [PubMed]

**10. **X. Ji, W. Cai, and P. Zhang, “High-order DGTD methods for dispersive maxwell’s equations and modelling of silver nanowire coupling,” Int. J. Numer. Meth. Eng. **69**, 308–325 (2007). [CrossRef]

**11. **H. S. Chu, W. B. Ewe, and E. P. Li, “Optical properties of a single-chain of elliptical silver nanowires,” in *Proceedings of the 7th IEEE International Conference on Nanotechnology* (IEEE2007), pp. 850–853.

**12. **H. S. Chu, W. B. Ewe, E. P. Li, and R. Vahldieck, “Analysis of sub-wavelength light propagation through long double-chain nanowires with funnel feeding,” Opt. Express **15**, 4216–4223 (2007). [CrossRef] [PubMed]

**13. **M. L. Brongersma, J. W. Hartman, and H. A. Atwater, “Electromagnetic energy transfer and switching in nanoparticle chain arrays below the diffraction limit,” Phys. Rev. B **62**, R16356–R16359 (2000). [CrossRef]

**14. **W. H. Weber and G. W. Ford, “Propagation of optical excitations by dipolar interactions in metal nanoparticle chains,” Phys. Rev. B **70**, 125429 (2004). [CrossRef]

**15. **D. S. Citrin, “Coherent excitation transport in metalnanoparticle chains,” Nano Lett. **4**, 1561–1565 (2004). [CrossRef]

**16. **A. F. Koenderink and A. Polman, “Complex response and polariton-like dispersion splitting in periodic metal nanoparticle chains,” Phys. Rev. B **74**, 033402 (2006). [CrossRef]

**17. **A. Alù and N. Engheta, “Theory of linear chains of metamaterial/plasmonic particles as subdiffraction optical nanotransmission lines,” Phys. Rev. B **74**, 205436 (2006). [CrossRef]

**18. **K. H. Fung and C. T. Chan, “Plasmonic modes in periodic metal nanoparticle chains: a direct dynamic eigenmode analysis,” Opt. Lett. **32**, 973–975 (2007). [CrossRef] [PubMed]

**19. **A. A. Govyadinov and V. A. Markel, “From slow to superluminal propagation: Dispersive properties of surface plasmon polaritons in linear chains of metallic nanospheroids,” Phys. Rev. B **78**, 035403 (2008). [CrossRef]

**20. **T. Yang and K. B. Crozier, “Dispersion and extinction of surface plasmons in an array of gold nanoparticle chains: influence of the air/glass interface,” Opt. Express **16**, 8570–8580 (2008). [CrossRef] [PubMed]

**21. **E. Simsek, “Full analytical model for obtaining surface plasmon resonance modes of metal nanoparticle structures embedded in layered media,” Opt. Express **18**, 1722–1733 (2010). [CrossRef] [PubMed]

**22. **M. Conforti and M. Guasoni, “Dispersive properties of linear chains of lossy metal nanoparticles,” J. Opt. Soc. Am. B **27**, 1576–1582 (2010). [CrossRef]

**23. **N. Talebi and M. Shahabdi, “Analysis of the propagation of light along an array of nanorods using the Generalized Multipole Techniques,” J. Comput. Theor. Nanosci. **5**, 711–716(6) (2008). [CrossRef]

**24. **A. Hochman and Y. Leviatan, “Rigorous modal analysis of metallic nanowire chains,” Opt. Express **17**, 13561–13575 (2009). [CrossRef] [PubMed]

**25. **Y. Leviatan and A. Boag, “Generalized formulations for electromagnetic scattering from perfectly conducting and homogeneous material bodies-theory and numerical solution,” IEEE Trans. Antennas Propag. **36**, 1722 –1734 (1988). [CrossRef]

**26. **C. Hafner, *The Generalized Multipole Technique for Computational Electromagnetics* (Artech House, Norwood, MA, 1990).

**27. **O. M. Bucci, G. D’Elia, and M. Santojanni, “Non-redundant implementation of Method of Auxiliary Sources for smooth 2D geometries,” Electronics Letters **41**, 1203–1205 (2005). [CrossRef]

**28. **G. Tayeb and S. Enoch, “Combined fictitious-sources-scattering-matrix method,” J. Opt. Soc. Am. A **21**, 1417–1423 (2004). [CrossRef]

**29. **D. Maystre, M. Saillard, and G. Tayeb, “Special methods of wave diffraction,” in *Scattering*,, P. Sabatier and E. Pike, eds. (Academic Press, 2001), chap. 1.5.6.

**30. **G. Fairweather and A. Karageorghis, “The Method of Fundamental Solutions for elliptic boundary value problems,” Adv. Comput. Math. **9**, 69–95 (1998). [CrossRef]

**31. **R. Petit, *Electromagnetic Theory of Gratings* (Springer-Verlag, 1980).

**32. **A. Boag, Y. Leviatan, and A. Boag, “Analysis of two-dimensional electromagnetic scattering from a periodic grating of cylinders using a hybrid current model,” Radio Sci. **23**, 612–624 (1988). [CrossRef]

**33. **A. Boag and Y. Leviatan, “Analysis of two-dimensional electromagnetic scattering from nonplanar periodic surfaces using a strip current model,” IEEE Trans. Antennas Propag. **37**, 1437 –1446 (1989). [CrossRef]

**34. **A. Ludwig and Y. Leviatan, “Analysis of bandgap characteristics of two-dimensional periodic structures by using the Source-Model Technique,” J. Opt. Soc. Am. A **20**, 1553–1562 (2003). [CrossRef]

**35. **M. Szpulak, W. Urbanczyk, E. Serebryannikov, A. Zheltikov, A. Hochman, Y. Leviatan, R. Kotynski, and K. Panajotov, “Comparison of different methods for rigorous modeling of photonic crystal fibers,” Opt. Express **14**, 5699–5714 (2006). [CrossRef] [PubMed]

**36. **Y. Leviatan and A. Boag, “Analysis of TE scattering from dielectric cylinders using a multifilament magnetic current model,” IEEE Trans. Antennas Propag. **36**, 1026 –1031 (1988). [CrossRef]

**37. **F. Capolino, D. Wilton, and W. Johnson, “Efficient computation of the 3D Green’s function for the Helmholtz operator for a linear array of point sources using the Ewald method,” J. Comput. Phys. **223**, 250 – 261 (2007). [CrossRef]

**38. **G. Valerio, P. Baccarelli, P. Burghignoli, and A. Galli, “Comparative analysis of acceleration techniques for 2-D and 3-D Green’s functions in periodic structures along one and two directions,” IEEE Trans. Antennas Propag. **55**, 1630–1643 (2007). [CrossRef]

**39. **A. Boag, Y. Leviatan, and A. Boag, “Analysis of diffraction from doubly periodic arrays of perfectly conducting bodies by using a patch-current model,” J. Opt. Soc. Am. A **7**, 1712–1718 (1990). [CrossRef]

**40. **F. Harris, “On the use of windows for harmonic analysis with the discrete Fourier transform,” Proceedings of the IEEE **66**, 51–83 (1978). [CrossRef]

**41. **W. Schroeder and I. Wolff, “The origin of spurious modes in numerical solutions of electromagnetic field eigenvalue problems,” IEEE Trans. Microw. Theory **42**, 644 –653 (1994). [CrossRef]

**42. **A. Hochman and Y. Leviatan, “Efficient and spurious-free integral-equation-based optical waveguide mode solver,” Opt. Express **15**, 14431–14453 (2007). [CrossRef] [PubMed]

**43. **P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” Phys. Rev. B **6**, 4370–4379 (1972). [CrossRef]

**44. ** Virginia Semiconductor Inc., “Optical properties of Silicon,” http://www.virginiasemi.com/pdf/OpticalPropertiesofSilicon71502.doc.

**45. **B. Prade and J. Y. Vinet, “Guided optical waves in fibers with negative dielectric constant,” J. Lightwave Technol. **12**, 6–18 (1994). [CrossRef]