Plasmon mode propagation in array of closely spaced metallic cylinders

Chain plasmonic waveguides are formed by linear arrays of metallic grains embedded in a dielectric matrix. Plasmonic structures of this kind have potential applications to subwavelength guiding, subwavelength imaging and SERS technology. We present qualitative analysis and numerical results for bound plasmonic modes propagating along the chain of closely spaced silver cylinders of subwavelength diameter. The dispersion relation and electromagnetic field structure of the modes are calculated by the cylindrical harmonic expansion method. We demonstrate that it is possible to match simultaneously both the frequency and wave number of the fundamental transverse mode and the first longitudinal mode in optical range. The application of dense chain of cylinders for optical switching between guided modes is discussed. © 2014 Optical Society of America OCIS codes: (230.7370) Waveguides; (160.3918) Metamaterials; (160.5298) Photonic crystals; (130.4815) Optical switching devices. References and links 1. D.K. Gramotnev and S.I. Bozhevolnyi, “Plasmonics beyond the diffraction limit,” Nat. Phot. 4, 83–91 (2010). 2. 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). 3. Stefan A. Maier, Pieter G. Kik, and Harry A. Atwater, “Optical pulse propagation in metal nanoparticle chain waveguides,” Phys. Rev. B 67(20), 205402 (2003). 4. 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,” Nat. Mater. 2, 229–232 (2003). 5. Y. Zhao, and Y. Hao, “Finite-difference time-domain study of guided modes in nano-plasmonic waveguides,” IEEE Trans. Antennas Propag. 55(11), 3070–3077 (2007). 6. M. Conforti, and M. Guasoni, “Dispersive properties of linear chains of lossy metal nanoparticles,” J. Opt. Soc. Am. B 27(8), 1576–1582 (2010). 7. S. M. Raeis Zadeh Bajestani, M. Shahabadi, and N. Talebi, “Analysis of plasmon propagation along a chain of metal nanospheres using the generalized multipole technique,” J. Opt. Soc. Am. B 28(4), 937–943 (2011). 8. N. A. Giannakis, J. E. Inglesfield, A. K. Jastrzebski, and P. R. Young, “Photonic modes of a chain of nanocylinders by the embedding method,” J. Opt. Soc. Am. B 30(6), 1755–1764 (2013). 9. I.L. Rasskazov, S.V. Karpov, and V.A. Markel, “Nondecaying surface plasmon polaritons in linear chains of silver nanospheroids,” Opt. Lett. 38(22), 4743–4746 (2013). 10. E. Smith and G. Dent, Modern Raman Spectroscopy: A Practical Approach (John Wiley and Sons, 2005). 11. K. Kneipp, M. Moskovits, and H. Kneipp, Surface-Enhanced Raman Spectroscopy: a Brief Perspective. In Surface-Enhanced Raman Scattering: Physics and Applications (Springer, 2006). 12. S. Kawata, A. Ono, and P. Verma, “Subwavelength colour imaging with a metallic nanolens,” Nat. Phot. 2, 438– 442 (2008). 13. V. Markel, and A. Sarychev, “Propagation of surface plasmons in ordered and disordered chains of metal nanospheres,” Physical Review B 75, 111 (2007) 14. C. R. Simovski and E. A. Yankovskaya, “Propagation of light along the waveguide of silver nano-cylinders,” Proc. SPIE 5927, Plasmonics: Metallic Nanostructures and Their Optical Properties III, 59271K (2005) 15. B. Rolly, N. Bonod and B. Stout, “Dispersion relations in metal nanoparticle chains: necessity of the multipole approach,” J. Opt. Soc. Am. B 29, 1012 (2012) 16. W. Zakowicz, “Two coupled dielectric cylindrical waveguides,” J. Opt. Soc. Am. A 14(3), 580–587 (1997). 17. R. Borghi, F. Gori, M. Santarsiero, F. Frezza and G. Schettini, “Plane-wave scattering by a perfectly conducting circular cylinder near a plane surface: cylindrical-wave approach,” J. Opt. Soc. Am. A 13(3), 483–493 (1996). 18. R. Borghi, M. Santarsiero, F. Frezza and G. Schettini, “Plane-wave scattering by a dielectric circular cylinder parallel to a general reflecting flat surface,” J. Opt. Soc. Am. A 14(7), 1500–1504 (1997). 19. K. Yasumoto, H. Toyama, and T. Kushta, “Accurate analysis of two-dimensional electromagnetic scattering from multilayered periodic arrays of circular cylinders using lattice sums technique,” IEEE Trans. Antennas Propag. 52,(10) 2603–2611 (2004). 20. S. Belan, S. Vergeles and P. Vorobev, “Adjustable subwavelength localization in a hybrid plasmonic waveguide,” Opt. Express 21(6), 7427–7438 (2013). 21. V.E. Babicheva, S.S. Vergeles, P.E. Vorobev and S. Burger, “Localized surface plasmon modes in a system of two interacting metallic cylinders,” J. Opt. Soc. Am. B 29, 1263–1269 (2012). 22. F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, Handbook of Mathematical Functions (NIST National Institute of Standards and Technology & Cambridge University Press, 2010). 23. P. Johnson, R. Christy, “Optical constants of the noble metals,” Phys. Rev. B 6, 4370 (1972). 24. E.N. Economou, “Surface Plasmons in Thin Films,” Phys.Rev. 182(2), 539 (1969). 25. D.J. Griffins, Introduction to Quantum Mechanics (Pearson Prentice Hall, 2d ed., 2005) 26. V. Twersky, “Elementary function representation of the Schlomilch series,” Arch. Rational Mech. Anal.8(1), 323–332 (1961). 27. C.M. Linton, “The Greens function for the two-dimensional Helmholtz equation in periodic domains,” Journal of Engineering Mathematics 33(4), 377–402 (1998). 28. D. V. Evans and R. Porter, “Trapping and near-trapping by arrays of cylinders in waves,” Journal of Engineering Mathematics 35(1-2), 149-179 (1999).


Introduction
Metallic nanostructures are known to exhibit strong subwavelength localization of the guided signals due to the excitation of the surface plasmon-polaritons (SPP) at the metal-dielectric interfaces.Various architectures of waveguides have been proposed for guiding SPP modes at optical wavelengths: metal nanorods, thin metal films, nanoholes in a metallic medium, nanogaps between metallic media, sharp metal wedges and nanogrooves, hybrid plasmonic waveguides [1].There is also a class of periodic geometries that can be used in subwavelength optics.In the so-called chain plasmonic waveguide the electromagnetic energy is transferred via the nearfield plasmon coupling between periodically arranged metal grains [2].The resulting plasmonic mode is highly confined to the vicinity of the line of the grains.The guiding of visible light using the arrays of metal nanoparticles has been demonstrated both numerically and experimentally [3,4,5,6,7,8].It has been also reported that for suitably chosen parameters the propagation of light in linear chains of metallic nanospheroids is almost free of spatial decay [9].Potential applications of the periodic arrays of metal grains are not restricted by the plasmon nanoguiding, the intensive research in electrodynamic properties of these structures is motivated also by its importance for SERS technology [10,11], subwavelength imaging and sensing [12].
The guiding properties of sufficiently sparse chains of grains can be effectively considered analytically with the aim of dipole approximation, see [13] for the case of spherical grains and [14] for the case of the ordered array of the metal wires.Clearly, the perturbation approach fails when the interwire distance is much less than the wire radius since the coupling of high-order multipoles becomes important.Using the multipole description, it was shown in [7,15] that the dipole approximation sufficiently deviates from the exact solution even for moderately dense array of spherical grains.
In this study we focus on the guiding properties of the very dense chain of infinitely long metal cylinders.We show that the typical frequency of guided modes decreases as the interwire distance diminishes and that for sufficiently small inter-wire distance the dispersion curves of two lowest modes of the chain have intersections in optical range.This feature allows to match both frequency and wave number of the bound modes with different field structures.
In our investigation, we use semi-analytical method which is based on the representation of the electromagnetic field of a guided mode as a coherent linear superposition of the localized plasmonics modes of separate wires.The approach is highly efficient for numerical solving the wave propagation and scattering problems in the composite systems [16,17,18,19,20].The accuracy of our results is verified by means of the finite element method from COMSOL.

Qualitative description
Generally the guided modes of chain waveguide can be considered as a coherent excitation of plasmonic modes of individual cylinders.Say, in the main approximation the fundamental mode of sufficiently sparse chain is described by the interacting dipole harmonics associated with transversal oscillations of electrical charge.However, when the interparticle distance gets smaller, the dipolar model becomes inappropriate because of the significant contribution of the high-order harmonics.The optical properties of closely packed particles are very different from those of distant particles.In the limit of extremely small inter-wire spacing the fundamental mode of the chain has nothing to do with the superposition of oscillating dipole harmonics of individual wires.The guiding mechanism in the very dense chain is the coherent excitation of strongly localized plasmonic modes in the narrow gaps between the cylinders, see Fig. 1.It is possible to estimate the typical frequency of a plasmonic mode, making a start from the solution for the localized plasmonic modes in the system of two metal cylinders separated by a small gap [21].The analysis is quite simple if the cylinders' radius R is much smaller than the wavelength λ , R λ .Then, the electromagnetic field near the cylinders can be treated in the quasi-static approximation.The field is concentrated inside the plane part of the gap region of length ∼ √ Rh where standing plasmon wave is excited, see Fig. 1.In this case, the dispersion of the mode stems mainly from the frequency dependence of the dielectric permittivities ε m and ε d for the metal and the surrounding dielectric media correspondingly.The frequency scaling with the large parameter R/h 1 of both the fundamental and the first longitudinal modes is estimated from the relation up to some factors of the order of unity.The factors can be found from the exact solution of the quasi-static equation for each mode, thus they are determined by the field structure on the scale L of the whole period.In particular this means that the result generally depends on quasimomentum q.Below we demonstrate that estimate (1) is valid for wave numbers which are far enough from the light line.

Semi-analytical approach
Let us consider an infinite periodic array of identical metal cylinders embedded in a homogeneous dielectric matrix.The adjacent cylinders are separated by the gap of width h.The radius of the cylinders is R, so the period of the array is L = 2R + h.We choose the Cartesian reference system as it is shown in Fig. 2: z-axis is directed along the cylinders, whereas x-axis is directed normally to the plane of the cylinders.We are interested in a plasmon-polariton TM-mode of frequency ω and the wave number q propagating along y-axis.In this case magnetic field H z is parallel to the cylinder axis and the only non-zero electric field components E x and E y are determined by H z .The guided mode of the chain can be represented as superposition of localized surface plasmonic modes of separate cylinders.Let m-th cylinders has centre mL in the y-direction.We represent the magnetic field in the outer dielectric as a double series with Hankel functions of the first kind and field inside a cylinder m as series of the modified Bessel functions of the first kind where (r m , ϕ m ) are the polar coordinates centered in the m-th cylinder, k = √ ε d ω/c and κ = √ −ε m ω/c.These expansions satisfy the two-dimensional Helmholtz equation (∆ 2 + εω 2 /c 2 )H z = 0 inside the homogeneous areas where permittivity is constant.The boundary conditions are the continuity of the field components tangential to each cylinders boundary.Note that each term in series (3) and (2) corresponds to the electromagnetic field produced by electric multipole inside the particular cylinder.
The Blochs periodic boundary condition (PBC) H z (x, y+ jL) = H z (x, y)e iqL j must be satisfied for the guided mode of the wave number q, where j is integer.Thus, the electromagnetic field of guided mode is completely determined by set of coefficients A ν ≡ A 0 ν and B ν ≡ B 0 ν .To eliminate the unknowns B ν associated with the internal field we impose the continuity conditions on the cylindrical surface using use Grafs addition theorem [22].As a result, we arrive to an infinite set of linear homogenious equations for the outer coefficients A ν .If the Ohmic losses inside the cylinders are negligible, the dispersion law is determined by the eigenvalue problem where Cn is determined in Eq. ( 15) and Fnν = Im F nν , see Eq. (8).During a numerical solution, one should keep first terms in sum (4) which resolve angles order of h/R with desired accuracy.The estimate follows from the analytical consideration of plasmonic modes for two closed Frequaency, eV Fig. 3.The metal permittivity at resonance frequency for the first two guided modes of the chain waveguide in comparison with the resonant permittivity in system of two cylinders, see Eq. ( 1).We choose q = π/L to suppress as much as possible the retardation effects.Note, that the surface charge distribution has two nodes for the transverse mode and four nodes for the longitudinal one at q = π/L, see Fig 1 .It follows from quasi-static solution [21] that the difference in the charge distribution leads to the longitudinal mode should have indeed higher frequency than the transverse one.
metallic wires [21].Due to symmetry x → −x, the linear system (4) is uncoupled for transverse and longitudinal modes, which correspond to amplitudes t n = A n + A −n and l n = A n − A −n respectively.More details of calculations and comparison with dipole approximation can be found in Appendix.

Numerical results and discussion
Next we apply the approach to investigate light guiding by a dense grid of silver cylinders, where the width of the gap is much smaller than the wires radius, h R. The frequency de-  pendence of silver permittivity in optical range is extracted from experimental data [23].First, we show that the estimate ( 1) is indeed correct.For the purpose, we plot the dependence of the resonance value of the metal permittivity as a function of cylinder radius R at fixed h = 1nm for the fundamental transverse and first longitudinal modes, see.Fig. 3.We take quasi-momentum q = π/L at the boundary of Brillouin zone to suppress retardation effects in the simulation.Second, we investigate the dispersion curves for the fundamental and first longitudinal modes.We track the dependence of dispersion curves on the thickness of the inter-wire gap at fixed wire radius R = 25nm and show that there is some critical value h c ≈ 5nm of the width, for which two dispersion curves touch each other, see Fig. 4(a).This means, that it is possible to match both the frequency and the wave vector of the transverse and longitudinal modes of the chain waveguide in optical range.For the smaller values of the width, there is the band of quasi-momentums, for which the dispersion curve of the longitudinal mode lies below that of the fundamental one, see Fig. 4(b).The phenomena has a qualitative explanation in terms of quasi-static approach.Indeed, at quasi-momentums qL 1 the surface charge distribution has only two nodes for longitudinal mode and four for transverse one in contrary with the opposite case qL π, see Fig. 3. Hence, for qL 1 the transversal mode corresponds to smaller absolute value of ε m , as it is for higher bound modes in a system of two closed metallic wires [21].
The found dispersion curves correspond to frequency range where the loss tangent of silver is indeed small, tg < 0.05, that justifies our approximation of lossless metal.Note that waveguide is subwavelength since the period L ∼ 60 nm is much smaller than the free-space wavelength λ which varies from ∼ 400 nm to ∼ 600 nm along the curves.In Fig. 5 we show two-dimensional plots of the electromagnetic field on the xy plane at wave number q = π/L.Thus, in the case of dense grid the field is strongly localized within the nanoscale gap between cylinders.We also plot schematic distribution of surface charges at the metallic boundary to reveal the vector structure of the electric field distribution, see Fig. 1.
It is informative to compare guiding along the chain of cylinders with signal guiding by the metal film of width 2R (see Fig. 2(b)) to show that waveguide modes in a dense array have substantial difference in comparison with those in a metal film [24].For both these systems the fundamental mode which has zero cut-off frequency is transverse one, while the second guided mode is longitudinal, Fig. 4(b).The difference is that the modes in the chain waveguide have lower frequency which generally is determined not solely by the thickness of the structure but also by the thickness of the inter-wire gap.This means that chain waveguide provides stronger mode confinement than metal film of the same thickness.
The intersection of dispersion curves is an unique feature of dense chain of metal cylinders.The plasmonic waveguides with such dispersion property seem to be unknown in literature.We believe that this feature can find practical applications.For example, if there is a mean to disturb the symmetry x → −x of the waveguide by an external force, then the considered waveguide can be exploited as an optical switcher between two lowest modes when both the frequencies and the wavenumbers of the modes coincide.

Conclusion
We have investigated a plasmonic waveguide based on ordered array of closely spaced metallic cylinders whose radius is deep below the diffraction limit.Interaction between adjacent cylinders in the structure leads to coherent energy transport along the array.The typical frequency of the lowest guided modes is determined by the surface plasmon resonance inside gap between two adjacent wires and thus decreases when the interwire distance decreases.It is shown that the dispersion curves for two lowest bound modes which are transverse and longitudinal have intersections when the inter-wire distance is small enough.In particular, we found the parameters for silver cylinders for which the dispersion curves touch each other.We argued that the intersections are the result of the confinement of the bound modes inside the narrow gaps between the adjacent cylinders.
The intersection of dispersion branches distinguishes the dense chain waveguide from the sparse one and can potentially be used for optical switching: if one can break symmetry x → −x of the waveguide applying an external force, then longitudinal mode is converted into transverse one during propagation along the waveguide and vise versa.The system of two propagating modes can be treated in the terms of two-level system, where the distance along the waveguide plays the role of the time and the energies of two states are equal due to equality of the quasi-momentums.The external force produces non-zero off-diagonal element associated with the transition between these states.To achieve the most effective conversion the length of the waveguide should correspond to half of Rabi cycle, see e.g.[25].The touch of the dispersive curves as in Fig. 4 produces widening of the frequency band where an approximate spatial phase synchronization can be achieved and thus the conversion may be implemented.Note that Ohmic losses lead to finite propagation length of the modes and thus the losses possesses restriction on the minimum value of the external force to meet a requirement of the conversion efficiency.
Infinite cylinders which we have considered here can be viewed as an idealized model of long rods or prolate spheroids directed along z-axis.In the case of long rods of the length 2b, we estimate that our 2D model gives well approximation for the 3D system when b L. Indeed, we showed, that the exact dispersion law of considered modes is determined by the field structure on the scale of whole period L. In the case of prolate ellipsoids which was considered in particular in [9] the criteria changes to b L R/h, where b is the longest semi-axis of the spheroids.The criteria follows from the condition the size of the flat area of the gap in z-direction should be large in comparison with the period of the lattice.We checked these estimates via direct numerical simulation of wave propagation in 3D chains of spheroids and rods, see Fig. 6(b).

Dispersion
Here we derive the characteristic equation ( 4) for the eigenproblem of the chain waveguide.In accordance with Grafs addition formula for the cylindrical functions where we have written H n for the Hankel function of the first kind H n .This expansion is valid for r < |m|L.Then the outer magnetic field close to the surface of cylinder "0" (for r < L) can be expressed solely in terms of in its local polar coordinates (r, ϕ) The unknown amplitudes A µ n can be expressed in the terms of the set of amplitudes A 0 n in accordance with BPC.Then, the last expression can be rewritten as where for convenience we have put A 0 n = A n and The other non-zero components or electromagnetic field are The continuity conditions for H z and E ϕ on surface of cylinder "0" lead to an infinite system of linear homogeneous equations After excluding coefficients B n we obtain where We can write F nν = e −iπ(ν−n)/2 S ν−n where the so called lattice sum S n is defined as Obviously, the relations S −n = (−1) n S n and F −n,−ν = F n,ν are valid.A problem with Eq. ( 13) is that the series over Hankel functions converge very slowly.However, it was shown by Twersky [26] that lattice sum can be rewritten in an alternative form which is rapidly convergent.In the lossless case (i.e., when Im ε m = 0) the bound plasmon mode has real wave number q under the light line, q > k = √ ε d ω/c.Let us consider the first Brillouin zone, 0 ≤ k < q < π/L.It can be seen from Eqs. ( 24)-( 26) that F nν = −δ nν + i Fnν , where Fnν is real [28].Thus, the characteristic equation (11) reduces to the real form where Cn = κY n (kR)I n (κR) + kY n (kR)I n (κR) κJ n (kR)I n (κR) + kJ n (kR)I n (κR) , J n and Y n are the Bessel functions of the first and second kind correspondingly and Fnν = Im F nν , see Eq. ( 8).
Next we should consider separately the symmetric and antisymmetric modes.Let us pass to coefficients The linear system ( 14) is uncoupled for these new variables Cn Solving the problem numerically it is useful to pass to the equation where Ã = | C| 1/2 A, and truncate it to some finite size N that controls numerical precision of the scheme.The characteristic determinant of the resulting set of equation vanishes along the dispersion curve ω(q) When h R the electromagnetic field of guided mode is strongly localized in the region where the gap can be considered as approximately plain.The transversal size of the plain part of the gap is evaluated as ∼ 2 √ 2Rh that corresponds to the central angel ∼ 2 2h/R.The number of cylindrical harmonics needed to resolve the field structure of such angular size is N ∼ π R/2h.Thus, the rate of convergence of proposed scheme is determined by the parameter R/h.In particular this means that our method does not work in the case of touching cylinders, h = 0.In numerical calculations the convergence analysis was carried out to ensure that the resonant permittivity ε m (ω) varied by less than 1% with the changing of truncation order.
Figure 6(a) demonstrates the fail of the electric dipole approximation in the limit of closely spaced cylinders, R = 25nm and h = 1nm.The dispersion relation of the fundamental propagating mode is calculated by numerically solving Eq. ( 21) with truncated orders N = 1 and N = 15.Dipole approximation gives wrong prediction of dispersion relation over the main part of the Brillouin zone, where the dispersion curve sufficiently deviates from the light line.
In order to test the semi-analytical approach describing here we firstly reproduced some of known results for sufficiently sparse grid of nm-size cylinders embedded in air, R = 25 nm and h = 25 nm.Following [5] we approximated the frequency dispersion of metal permittivity by the lossless Drude formula ε m (ω) = 1 − ω 2 p /ω 2 where the plasma frequency for silver is ω p = 6.18 eV .Our results showed good agreement with those obtained in [5] by means of a conformal dispersive finite-difference time-domain (FDTD) method.We also verified the results by means of finite-element package FEMLab from COMSOL.The numerical model is based on the resonance excitation of the guided plasmonic modes by evanescent electromagnetic field.
In this study the imaginary part ε m of metal permittivity ε m = ε m + iε m was neglected, but the approach could be generalized to the case of non-zero material losses.For lossy cylinders the guided mode has complex wave vector q = q + iq and decays during propagation.Dispersion relation can be found by fixing frequency ω and finding numerically complex wave vector corresponding to a zero of the complex characteristic determinant.The procedure is straightforward but leads to some additional computational efforts.The height of cylinders is 2b.The transversal geometrical parameters was chosen to be R = 25nm and h = 5nm.Thus, whereas our 2D model is good approximation for chain of cylinders at moderate aspect ratios, it is weaker approximation for chain of prolate spheroids.Solution for 3D model was obtained using COMSOL.

Field structure
The magnetic field in the vicinity of the cylinder "0" is where r < L. The mode is symmetric about the line joining the cylinder centres (which is x = 0) when A −n = A n .In the case of antisymmetric mode A 0 = 0 and A −n = −A n .

Lattice sums
According to [26] and [27] S 0 = −1− 2i where B m is the Bernoulli polynomial, γ = 0.577 is Eulers constant and for the sake of brevity we define

Fig. 1 .
Fig. 1.Charge distribution at wave number q = π/L for a) the fundamental mode and b) the first longitudinal mode.

Fig. 2 .
Fig. 2. a) The periodic grid metal cylinders.b) The continuous metal film of the same thickness.

Fig. 4 .
Fig. 4. The dispersion diagram for a dense periodic chain of silver cylinders embedded in air.The frequency dependence of silver permittivity is taken from experimental data [23].a) The set of parameters R = 25 nm, h = 5 nm corresponds to touching of two dispersion curves.b) The case of very dense array, R = 25 nm, h = 1 nm.The point of intersection of the dispersion branches corresponds to frequency ω = 2.96 eV (λ = 418 nm) and wave number q = 0.426π/L.

Fig. 5 .
Fig. 5. Spatial distribution of electromagnetic field for first two guided mode of dense chain of silver cylinders at wave number q = π/L.The radius of the cylinders is R = 25 nm and the gap width is h = 1 nm.(a) and (c) normalised distributions of total electric and magnetic field of transverse mode; (b) ans (d) normalised distributions of total electric and magnetic field of longitudinal mode.

Fig. 6 .
Fig. 6. a) The dispersion relation of fundamental mode calculated by numerical solving the dispersion equation with different values of truncated index: N = 1 (dipole approximation) and N = 15 (full multipole calculations) b)The frequency of the longitudinal mode at qL/π = 0.4: comparison between the 2D model of infinite cylinders (dashed line), 3D chain of rods (asterisks) and 3D chain of spheroids (dots).The longest semi-axis of spheroids is b.The height of cylinders is 2b.The transversal geometrical parameters was chosen to be R = 25nm and h = 5nm.Thus, whereas our 2D model is good approximation for chain of cylinders at moderate aspect ratios, it is weaker approximation for chain of prolate spheroids.Solution for 3D model was obtained using COMSOL.