## Abstract

Nanowire chains (NCs) are analyzed by use of a rigorous, full-wave, Source-Model Technique (SMT). The technique employs a proper periodic Green’s function which converges regardless of whether the structure is lossless or lossy. By use of this Green’s function, it is possible to determine the complex propagation constants of the NC modes directly and accurately, as solutions of a dispersion equation. To demonstrate the method, dispersion curves and mode profiles for a few NCs are calculated.

© 2009 Optical Society of America

## 1. Introduction

Linear chains of metallic nano-particles and nanowires have been the subject of intensive research recently. This effort has been motivated, primarily, by the possible use of these structures as integrated optical waveguides that can confine light transversally to sub-wavelength dimensions [1]. However, metallic structures are inherently lossy at optical frequencies, and there appears to be a correlation between losses and confinement [2]. Hence, reliable methods for determining these parameters are essential for waveguide design.

A straightforward and rigorous way to study the confinement and the effect of losses is to perform a modal analysis and determine the real and imaginary parts of the modal propagation constants, assuming a real-frequency time-harmonic variation of all fields and sources. The distance of the real part of the propagation constant from the light line is a measure of confinement [2], and the imaginary part of the propagation constant is a measure of attenuation due to losses (leaky modes, in which radiation contributes to the attenuation are encountered in the antenna literature [3,4, Ch. 19] and will not be considered here). Modal analyses of this type have been carried out for a variety of metallic optical waveguides such as the single metallic interface and metallic slab waveguides [5], and the cylindrical [6] and elliptical [7] rod waveguides. All these guiding structures are invariant along the direction of propagation. When a similar analysis is attempted for a guiding structure that is periodic along the direction of propagation, certain difficulties may arise.

One of the common ways to analyze periodic structures involves the assumption that the fields outside the periodic elements are due to periodic arrays of sources. For example, arrays of nano-particles are often modeled as dipole arrays [8]. Periodic arrays of cylinders (of arbitrary cross-section) can be modeled by use of periodic arrays of current filaments, as in [3], or periodic arrays of current strips, as in [9]. When the structure is lossy, the modal fields must decay exponentially in the direction of propagation in accordance with a complex propagation constant. However, the fields due to an infinite array of periodic sources with exponentially decreasing (increasing in the opposite direction) amplitudes must diverge everywhere, since the fields of the elementary sources decay only algebraically with distance.

A number of ways around this problem have been proposed. In [8], the infinite array is truncated to a finite length array. While this solves the problem, it has the drawback that a large number of particles must be modeled. In [10], the poly-logarithm function is used to continue analytically the values of an approximate dispersion equation from the real-line into the complex propagation-constant plane, where it would have otherwise diverged. This method is also used in [11, 12]. A different approach is to solve a scattering problem and to infer the dispersion curve from the variation of the reflection with the angle of incident light, as in [13, 14]. It is also possible to approximately determine the propagation length from the attenuation of the total fields (scattered+incident) along a finite chain as in [15]. In a scattering problem, the divergence does not arise, even if the structure is lossy and infinite, because the incident wave impinges on all the scatterers with the same amplitude.

In this paper, we focus on nanowire chains (NCs), i.e., the periodic elements are assumed to be periodic cylinders invariant in some direction transverse to the propagation direction. As this problem is 2D, it is somewhat simpler than the nano-particle chain problem. But the mentioned divergence occurs in the 2D case as well, and the techniques we propose could be applied to the 3D case.

The method of solution is the Source-Model Technique (SMT), which belongs to the Methods of Moments class of numerical techniques. The SMT, which has been used in various electromagnetic modeling scenarios [16–19], is closely related to a number of other methods: the Generalized Multipole Technique [20], the Method of Auxiliary Sources [21], the Method of Fictitious Sources [22, 23, Ch. 1.5.6], and the Method of Fundamental Solutions [24]. The earliest works are probably those by Kupradze [25] and Vekua [26], and a review can be found in [27]. The SMT was used in [9] for analysis of electromagnetic scattering from lossless periodic structures, and we adopt the formulation put forward there. In this formulation, the determination of the fields in all space is reduced to the consideration of a unit cell of the periodic structure, with periodic boundary conditions on the unit cell edges. With this view-point in mind, the modal fields can be defined as source-free solutions that obey periodic boundary conditions at the edges of the unit cell. When losses are present, the periodic boundary conditions dictate that the modal fields be the same on both edges of the unit cell, up to multiplication by a non-unit-amplitude complex scalar.

To represent the fields outside the cylinder, a proper periodic Green’s function (PPGF) is used. This function is closely related to the fields due to an elementary source operating in a unit cell, with the cylinder removed, but with the periodic boundary and radiation conditions enforced. The PPGF can be written analytically as a sum of Floquet harmonics that converges everywhere except at the source. This spectral representation converges even when the conventional spatial representation would diverge. Thus, by use of the PPGF, the modal solutions can be determined rigourously.

The remainder of this paper is organized as follows: the problem under study is specified in the next section, and the SMT is formulated in Section 3. In Section 4, the spectral representation of the PPGF, which is used in the SMT formulation, is given and its convergence is discussed. In Section 5, it is explained how the SMT formulation is used for mode determination. Numerical results for a number of NCs are given in Section 6, and the paper is summarized in Section 7.

## 2. Problem specification

We consider a linear chain of cylinders of arbitrary cross-section, all of them characterized by a relative 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 characterized by free-space permittivity,

*ε*

_{0}, and free-space permeability,

*µ*

_{0}. For all fields and sources, exp(

*jωt*) harmonic time variation is assumed and suppressed. Our aim is to determine proper modes of this structure, i.e., source-free electromagnetic fields which obey the periodicity condition

where *F*(*x*) is any of the field components and *k _{x}* is a propagation constant. As the modes are proper, they must obey the radiation condition. This means that for

*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 [28, p. 54]). For a given

*ω*, modes will exist only for certain values of

*k*. The determination of the (

_{x}*k*) pairs for which a mode exists is a central part of the mode determination process.

_{x},ωIn this paper, we will consider metallic cylinders characterized by a plasma-type permittivity. Consequently, we will only be interested in 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 therefore shorted-out.

Lastly, the mode determination problem may be reduced, by the Floquet theorem, to consideration of a single unit cell with periodic boundary conditions

at the unit cell boundaries.

## 3. Source-Model Technique (SMT)

The following formulation of the SMT for periodic structures follows [9], in which scattering from a linear chain of cylinders is analyzed. In the SMT, the fields in each homogenous region are approximated by the fields due to a set of elementary sources, of yet to be determined amplitudes, that is placed outside of each homogeneous region and is used to approximate the fields in that region. To approximate modes that are TE to *z*, all the sources are magnetic current filaments oriented parallel to the *z* axis. The fields within the unit cell and outside of the cylinder are approximated by the fields due to magnetic current filaments, of yet to be determined current amplitudes, operating in the unit cell with the cylinder removed and the periodic conditions and radiation conditions enforced. This is shown in Fig. 2. Generally, the filaments are distributed uniformly on curves which are slightly contracted and dilated versions of the media boundary curves. A simple method for generating such curves for arbitrary smooth boundary curves is described in [29]. The fields inside the cylinder are approximated by the fields due to magnetic current filaments operating in a homogeneous, unbounded, obstacle-free medium that has the same material parameters as the cylinder, as shown in Fig. 3.

Enforcing the continuity of the field components tangential to the media boundaries at a discrete set of testing points 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, *M* is the number of testing points, which are uniformly distributed on the media boundaries, and enforcing the continuity of the tangential electric and magnetic fields leads to 2*M* equations. The number of current filaments is 2*N*, of which there is an equal number inside and outside of the cylinder. To avoid linear dependence of the columns of [*Z*], as *N* is increased, the sources should approach the boundary. A simple rule of thumb used in this paper (adapted from [20, pp. 169–170]) is to make the length of the contracted and dilated curves on which the sources are placed equal to 1±2*π/N* times the length of the boundary curves.

Using more testing points than sources spreads the error in the continuity conditions more uniformly along the boundary. It is also useful for avoiding certain spurious solutions [30]. The structure of the impedance matrix is as follows

where the entries of [*Z*
^{in}
_{H}] are the magnetic fields at the testing points due to current filaments of unit amplitude that lie inside the cylinder and those of [*Z*
^{out}
_{H}] are the magnetic fields at the testing points due to current filaments that lie outside the cylinder. The sub-matrices [*Z*
^{in}
_{E}] and [*Z*
^{out}
_{E}] are analogous to [*Z*
^{in}
_{H}] and [*Z*
^{out}
_{H}], except that their entries are the tangential electric fields (divided by *η*
_{0}) at the testing points.

The entries of the sub-matrices [*Z*
^{out}
_{H}] and [*Z*
^{out}
_{E}] are easy to evaluate as they can be obtained from the 2D Green’s function

where *H*
^{(2)}0(·) denotes the Hankel function of the second kind and zero order and ${k}_{c}=\sqrt{{\epsilon}_{c}}{k}_{0}$ is the wave number inside the cylinder, with *k*
_{0} being the free-space wave number. In (5), the filament coordinates are denoted by (*x*
_{0},*y*
_{0}) and the testing point coordinates are denoted by (*x,y*). As a function of *G*
_{out}, the fields are given by

where ${\eta}_{c}=\sqrt{{\mu}_{0}\u2044{\epsilon}_{c}}$ is the characteristic impedance inside the cylinder. The entries of the sub-matrices [*Z*
^{in}
_{H}] and [*Z*
^{in}
_{E}] are given in terms of the PPGF, as detailed in Section 4.

## 4. Proper Periodic Green’s Function (PPGF)

The PPGF for a *z* invariant geometry is defined as the solution to the following boundary-value problem

$$\mathrm{subject}\phantom{\rule[-0ex]{.2em}{0ex}}\mathrm{to}:$$

$${G}_{\mathrm{in}{\mid}_{x=-L\u20442}}={e}^{{\mathrm{jk}}_{x}L}{G}_{\mathrm{in}{\mid}_{x=L\u20442},\text{}}$$

$$\mathrm{Radiation}\phantom{\rule[-0ex]{.2em}{0ex}}\mathrm{condition}\phantom{\rule[-0ex]{.2em}{0ex}}\mathrm{for}\phantom{\rule[-0ex]{.2em}{0ex}}y\to \pm \infty .$$

This boundary-value problem, depicted in Fig. 4, is a standard Sturm-Liouville problem, which leads to a well-known spectral representation of the periodic Green’s function. The case of complex *k _{x}*, however, has not received much attention. The solution of the boundary-value problem is given by (see for example [28, p. 23])

where

The function *G*
_{in} has a denumerably infinite number of square-root-type branch points that occur at the propagation constants

corresponding to Wood’s anomalies. For each of the square-roots in (12), a choice of branch and branch-cut must be made. To comply with the radiation conditions at *y*→±∞, we must have, for *k _{yn}* purely real, Re(

*k*)≥0, and Im(

_{yn}*k*)<0 otherwise. This implies using the proper Riemann sheet with branch-cuts as shown in Fig. 5. On this Riemann sheet, Gin is periodic in

_{yn}*k*with a period of 2

_{x}*π/L*. Consequently, we can assume that

*k*is in the First Brillouin Zone (FBZ), i.e., between −

_{x}*π/L*and

*π/L*.

The fields of a unit-amplitude magnetic current filament, operating in the presence of the periodic boundary and radiation conditions, can be readily obtained from the PPGF. We have

where ${\eta}_{0}=\sqrt{{\mu}_{0}\u2044{\epsilon}_{0}}$ denotes the characteristic impedance of free-space. This current may source or sink power depending on the parameters of the boundary-value problem (9).

The series (10) converges everywhere except at *x*=*x*
_{0},*y*=*y*
_{0}, and when *y*≠*y*
_{0} the convergence is exponential. When *y*≈*y*
_{0}, however, the rate of convergence is slow. This is a well-known problem of periodic Green’s functions, for which a number of solutions have been proposed [31]. When the SMT is used, this problem has a simple solution which was proposed in [9]. The idea is to replace the current filaments with current strips of magnetic current density

where *f*(*x*) is a smooth, real-valued window function that is zero outside the interval [−*s*/2, *s*/2], with s being the width of the strip. The *H _{z}* field (from which

*E*and

_{x}*E*can be derived) due to such a current strip is given by

_{y}where the Fourier coefficients, *f*̂_{n}, are given by

By choosing *f*(*x*) to have rapidly decaying Fourier coefficients the evaluation of the fields can be accelerated. In this paper, we used the Blackman-Harris [32] window, which was used in [17] and is given by

The width of the strip *s* should be as large as possible, since a larger width increases the rate of decay. Of course, the strip should be small enough so as not to intersect or even approach the boundary. We took the strip width to be a quarter of the distance from the center of the strip to the boundary.

## 5. Determination of modal solutions

A modal solution exists only for certain pairs of *k _{x}* and

*k*

_{0}for which a non-trivial solution to (3) exists. As [

*Z*] is a non-square matrix, we cannot use its determinant to test whether a solution exists for a given pair of

*k*and

_{x}*k*

_{0}. Instead, the condition number of [

*Z*] could be used and it should be infinite for a solution to exist. In practice, as the solutions are always approximate, the condition number is never infinite. Rather, it is a local maximum of the condition number as a function of

*k*and

_{x}*k*

_{0}that indicates the existence of an approximate solution. So, a naive method for determining the modes would be to search for the maxima of the condition number. However, as explained in [7], the determination of the modes can be significantly hindered by spurious solutions which are present in various integral-equation formulations and are especially detrimental in the SMT. These spurious solutions are not discrete solutions, but they form a continuum of approximate non-trivial solutions to (3). They are sets of non-zero filament amplitudes that when used to expand the fields yield fields that are very small everywhere except near the sources. A simple example of such a set is

*K*⃗=[1,-1,1,-1, …]

^{T}, where

*T*denotes transposition. These amplitudes are used to form approximations of the fields for each homogeneous region of the unit cell. When two such approximations corresponding to neighboring regions are evaluated at the testing points, which are some distance away from the sources, it is found that the fields of both approximations are very small. Hence, the absolute error in the continuity conditions associated with such a solution is small and this is manifested by a high condition number. Nevertheless, this solution is not a true solution because the absolute error in continuity conditions is not small compared with the values of the fields at the testing points.

To filter out the spurious solutions, we use the measure of singularity proposed in [7]. This measure of singularity is based on the normalized error Δ*E*(*K*⃗), which is defined as the absolute error normalized to the norm of the fields at the testing points. The idea in [7] is to exploit the fact that Δ*E*(*K*⃗) is small for a true solution and large for a spurious solution, and to devise a measure of singularity which, when applied to [*Z*], is large if and only if there exists some solution vector *K*⃗ for which Δ*E*(*K*⃗) is small. We are therefore interested in finding the smallest Δ*E*(*K*⃗) that can be obtained for a given matrix [*Z*]. Once this smallest Δ*E*(*K*⃗) is obtained, we use its reciprocal as a measure of the singularity of [*Z*]. Fortunately, finding the smallest Δ*E*(*K*⃗) is not significantly more difficult than evaluating the condition number of [*Z*]. The smallest Δ*E*(*K*⃗) can be written as

where [*Z*̃] is a matrix that maps the current amplitude vector *K*⃗ to tangential field values at the testing points. The matrix [*Z*̃] can be readily obtained from the matrix [*Z*], which maps *K*⃗ to the *difference* of the tangential fields on both sides of the media boundaries, by reversing the signs of half of its entries, and dividing by two. For TE fields, such as the ones considered in this paper, the matrix [*Z*̃] can be written as

This choice differs from that used in [7], where the fields were hybrid (neither TE nor TM). The significance of this choice is that [*Z*̃]*K*⃗ is a vector of the magnetic field values at the testing points, and it is the norm of this vector that is used to normalize the absolute error. Only the magnetic field is used because, for TE fields, the electric field due to a spurious solution is significantly larger than the magnetic field multiplied by the characteristic impedance of the medium in which the sources operate. Hence, if the electric field were included in the denominator of the expression for Δ*E*(*K*⃗), this denominator would not be small for a spurious solution, and consequently, the rejection of the spurious solutions would be less reliable.

As mentioned above, we are interested in finding the solution to (21), namely, the smallest Δ*E*(*K*⃗) that can be obtained for a given [*Z*]. It can be shown [33, Ch. 7] that the smallest Δ*E*(*K*⃗) is the square-root of the smallest generalized eigenvalue of the following generalized eigenvalue problem

where the dagger sign indicates conjugate-transpose. Denoting the smallest generalized eigenvalue by ξ_{min}, the smallest Δ*E*(*K*⃗) is simply $\sqrt{\xi min},$, and the measure of singularity is $1\u2044\sqrt{\xi min}.$.

Lastly, to determine the pairs of *k*
_{0} and *k _{x}* for which a mode exists, we search for the local maxima of $1\u2044\sqrt{{\xi}_{min}({k}_{x},{k}_{0})}$, and this search can be carried out effectively because $1\u2044\sqrt{{\xi}_{min}({k}_{x},{k}_{0})}$ is not influenced by the spurious solutions. The search for the maxima is carried out by use of the adaptive search algorithm described in [7].

## 6. Numerical results

To demonstrate the capabilities of the proposed modeling scheme, a few numerical results are given. Dispersion curves for a NC of circular cylinders, of radius *R*=50 nm and period *L*=120 nm, are shown in Fig. 6. The parameter varied in these curves is *λ*, and the values shown are the real and imaginary parts of the effective refractive index *n*
_{eff}=*k _{x}/k*

_{0}. In this and all numerical examples, the cylinders are assumed to be made of silver, which is characterized by a Drude-model permittivity function given by [34]

with *τ*=1.45×10^{−14}s and *ω _{p}*=1.32×10

^{16}rad/s. As can be readily observed in Fig. 6, at long wavelengths, the dispersion curve of the fundamental mode approaches the dispersion curve of a mode propagating along a single flat interface between silver and air. However, as the wavelength decreases, the dispersion curve of this mode, and those of higher order modes, differ significantly from the single interface dispersion curve. As expected, there is a trade-off between attenuation and confinement to the NC. Specifically, when Re(

*n*

_{eff}) ≈1, which implies that the mode propagates largely in the air surrounding the NC, the attenuation is small. As Re(

*n*

_{eff}) increases, which implies that the mode becomes more confined to the NC, the attenuation increases as well. Also, the rapid increase in attenuation near the edge of the FBZ is related to the vanishing group velocity there [35].

A few mode profiles, corresponding to four points on the dispersion curves of Fig. 6, are shown in Fig. 7. These points are marked A-D in Fig. 6; to avoid clutter, the points are marked only in Fig. 6(a). In the SMT, the solutions obey Maxwell’s equations exactly and the continuity conditions approximately. Therefore, the continuity of Re(*H _{z}*), which graphically appears perfect in Fig. 7, attests to the validity of the solutions. As could be anticipated from the symmetry of the NC, the modes are either odd or even with respect to

*y*. The attenuation is not apparent in these short segments of the NC, but it can be observed in a longer section shown in Fig. 8. Also evident in Fig. 8, is that the

*x*-directed time-average power flow density assumes negative values inside the cylinders. This is a well-known phenomenon in plasmonic waveguides [36].

We turn now to the analysis of a pair of NCs, which are separated by a vertical distance equal to the period *L*=120 nm. The dispersion curves for the NC pair are shown in Fig. 9. The *H _{z}* of the fundamental mode, shown in Fig. 10(a), is an even function of y, while the next mode, shown in Fig. 10(b), is an odd function of

*y*. In the case of the even mode, the power flows mostly in the air region between the NCs, whereas in the odd mode, the power flows mostly in the exterior air region.

The last structure considered is a NC shown in the inset in Fig. 11. The cylinders in this structure are of a more general cross-section defined by hypotrochoid curves, which have the following parametric representation

where *ϕ* is a parametric variable ranging from 0 to 2*π*, and *R _{b}* is the radius of a circle that tightly bounds the curve, which is a hypotrochoid. For

*ν*≥1, the hypotrochoid resembles a triangle with rounded corners, and the radius of curvature of the rounded corners increases with increasing

*ν*. Modeling a sharp triangular shape would require careful placement of the sources near the corners, such as was done in [37] or [34]. Hence, we used

*ν*=2 to obtain a shape which is triangular but not too sharp at the corners.

The dispersion curves of this NC are quite similar to those of the circular cylinder NC (shown in Fig. 6). The mode profiles at points A and B, shown in Fig. 12, are also quite similar to their counterparts in Fig. 7. It may be interesting to look at the influence of the shape of the NC on some parameters of interest, such as the real and imaginary parts of the effective indices. In Fig. 13, we show the difference, Δ*n*
_{eff}≜*n*∘_{eff}−*n*
^{▷}
_{eff}, between the effective indices of the fundamental modes of the hypotrochoidal-cylinder NC of Fig. 11 and a circular-cylinder NC. To try and isolate the effect of the shape from the effect of the size of the cylinders, we choose them to have the same circumference. As can be inferred from Fig. 13, the attenuation rates can differ quite significantly, while the phase velocities differ very little. At longer wavelengths, the fundamental mode of the circular-cylinder NC is attenuated slightly more rapidly than the fundamental mode of the hypotrochoidal-cylinder NC, while at shorter wavelengths it is the other way around. The difference becomes quite significant as the edge of the FBZ is approached.

#### 6.1. Accuracy and computational resources

To give an idea of the trade-off between accuracy and computational resources, we calculated the effective index at point A in Fig. 6, with an increasing number of sources, *N*. The calculation involves searching for the complex effective index that yields the smallest error in the continuity conditions, where each evaluation of this error function requires the formation of the matrix [*Z*] and the solution of the generalized eigenvalue problem (23). The search begins with an adaptive sampling of the error function at real values of *n*
_{eff} ranging from 1 to *λ*/(2*L*). Once a minimum is found, a small region of the complex plane near this minimum is searched by a standard optimization algorithm (see [7] for more details regarding this search algorithm).

The results of the calculation with increasing *N* are shown in Table. 6.1, where *T* is the computation time for a MATLAB implementation running on a 3GHz PC. The number of sources

is doubled from one line in the table to the next, and Δ*E* is seen to decrease by about an order of magnitude. The computation time increases at first linearly, and for *N*=40, a bit more sharply. Lastly, the error in continuity conditions appears to be a reliable indicator of the accuracy of the solution. Throughout this work, it was ensured that the error Δ*E* did not exceed 2%. This required *N*=20 for the analysis of the circular-cylinder NC, *N*=30 for the NC pair, and *N*=40 for the hypotrochoidal-cylinder NC. The Green’s function summations in (10) were truncated when the partial sum was at least 10^{5} times greater (in absolute value) than the difference between consecutive terms.

## 7. Summary

A method for rigorous modal analysis of plasmonic NCs was described. The method is based on the SMT and it employs a spectral representation of the periodic Green’s function, which converges for real and complex propagation constants. Three representative NC structures were analyzed: a circular-cylinder NC, a pair of circular-cylinder NCs placed close together, and a hypotrochoidal-cylinder NC. It was shown that in the fundamental mode of the NC pair most of the power is confined to the region between the NCs. Also, it was shown that the attenuation rate of the fundamental mode of the hypotrochoidal-cylinder NC can differ significantly from that of the fundamental mode of a circular-cylinder NC with the same circumference. Lastly, the trade-off between accuracy and computational resources was outlined.

## Acknowledgements

We would like to thank the anonymous reviewers for their valuable comments and suggestions. This work was supported in part by the Israel Science Foundation and the Russell Berrie Nanotechnology Institute (NTU-Technion grant). AH gratefully acknowledges the financial support of the Andrew and Erna Finci Viterbi Foundation.

## References and links

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

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

**3. **F. Capolino, D. R. Jackson, and D. R. Wilton, “Fundamental properties of the field at the interface between air and a periodic artificial material excited by a line source,” IEEE Trans. Antennas Propag. **53**, 91–99 (2005). [CrossRef]

**4. **R. E. Collin and F. J. Zucker, *Antenna Theory, Part II*, McGraw-Hill (1969).

**5. **J. Burke, G. Stegeman, and T. Tamir, “Surface-polariton-like waves guided by thin, lossy metal films,” Phys. Rev. B **33**, 5186–5201 (1986). [CrossRef]

**6. **L. Novotny and C. Hafner, “Light propagation in a cylindrical waveguide with a complex, metallic, dielectric function,” Phys. Rev. E **50**, 4094–4106 (1994). [CrossRef]

**7. **A. Hochman and Y. Leviatan, “Efficient and spurious-free integral-equation-based optical waveguide mode solver,” Opt. Express **15**, 14,431–14,453 (2007). [CrossRef]

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

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

**10. **D. Citrin, “Plasmon-polariton transport in metal-nanoparticle chains embedded in a gain medium,” Opt. Lett. **31**, 98–100 (2006). [CrossRef]

**11. **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]

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

**13. **Q. H. Wei, K. H. Su, S. Durant, and X. Zhang, “Plasmon Resonance of Finite One-Dimensional Au Nanoparticle Chains,” Nano Lett. **4**, 1067–1072 (2004). [CrossRef]

**14. **T. Yang and K. 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]

**15. **H. Chu, W. Ewe, W. Koh, and E. Li, “Remarkable influence of the number of nanowires on plasmonic behaviors of the coupled metallic nanowire chain,” Appl. Phys. Lett. **92**, 103,103 (2008). [CrossRef]

**16. **Y. Leviatan, A. Boag, 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]

**17. **A. Boag, Y. Leviatan, and A. Boag, “Analysis of electromagnetic scattering from linear periodic arrays of perfectly conducting bodies using a cylindrical-current model,” IEEE Trans. Antennas Propag. **39**, 1332–1337 (1991). [CrossRef]

**18. **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]

**19. **A. Hochman and Y. Leviatan, “Analysis of strictly bound modes in photonic crystal fibers by use of a source-model technique,” J. Opt. Soc. Am. A **21**, 1073–1081 (2004). [CrossRef]

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

**21. **O. M. Bucci, G. D’Elia, and M. Santojanni, “Non-redundant implementation of method of auxiliary sources for smooth 2D geometries,” Electronics Letters **41**(22), 1203–1205 (2005). [CrossRef]

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

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

**24. **G. Fairweather and A. Karageorghis, “The method of fundamental solutions for elliptic boundary value problems,” Advances in Computational Mathematics **9**(1), 69–95 (1998). [CrossRef]

**25. **V. D. Kupradze, “About approximate solutions of a mathematical physics problem,” Success of Mathematical Sciences **22**(2), 59–107 (1967).

**26. **I. N. Vekua, Reports of the Academy of Science of the USSR 44(6), 901–909 (1953).

**27. **D. I. Kaklamani and H. T. Anastassiu, “Aspects of the Method of Auxiliary Sources (MAS) in computational electromagnetics,” IEEE Antennas and Propag. Mag. **44**, 48–64 (2002). [CrossRef]

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

**29. **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]

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

**31. **A. Peterson, L. Scott, and R. Mittra, *Computational Methods for Electromagnetics*, IEEE Press, (1998).

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

**33. **R. Bellman, *Introduction to matrix analysis*, McGraw-Hill, (1970).

**34. **E. Moreno, D. Erni, C. Hafner, and R. Vahldieck, “Multiple multipole method with automatic multipole setting applied to the simulation of surface plasmons in metallic nanostructures,” J. Opt. Soc. Am. A **19**, 101–111 (2002). [CrossRef]

**35. **K. C. Huang, E. Lidorikis, X. Jiang, J. D. Joannopoulos, K. A. Nelson, P. Bienstman, and S. Fan, “Nature of lossy Bloch states in polaritonic photonic crystals,” Phys. Rev. B **69**, 195,111 (2004). [CrossRef]

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

**37. **S. Eisler and Y. Leviatan, “Analysis of electromagnetic scattering from metallic and penetrable cylinders with edges using a multifilament current model,” IEE-Proc. H **136**, 431–438 (1989).