Universal scaling of local plasmons in chains of metal spheres

The position, width, extinction, and electric field of localized plasmon modes in closely-coupled linear chains of small spheres are investigated. A dipole-like model is presented that separates the universal geometric factors from the specific metal permittivity. An electrostatic surface integral method is used to deduce universal parameters that are confirmed against results for different metals (bulk experimental Ag, Au, Al, K) calculated using retarded vector spherical harmonics and finite elements. The mode permittivity change decays to an asymptote with the number of particles in the chain, and changes dramatically from 1/f 3 to 1/f 1/2 as the gap fraction (ratio of gap between spheres to their diameter), f, gets smaller. Scattering increases significantly with closer coupling. The mode sharpness, strength and electric field for weakly retarded calculations are consistent with electrostatic predictions once the effect of radiative damping is accounted for. ©2010 Optical Society of America OCIS codes: (240.6680) Surface plasmons; (250.5403) Plasmonics; (310.6628) Subwavelength structures, nanostructures. References and links 1. U. Kreibig, and M. Vollmer, Optical properties of metal clusters (Springer-Verlag, 1995). 2. P. K. Jain, W. Y. Huang, and M. A. El-Sayed, “On the universal scaling behavior of the distance decay of plasmon coupling in metal nanoparticle pairs: a plasmon ruler equation,” Nano Lett. 7(7), 2080–2088 (2007). 3. S. A. Maier, P. G. Kik, and H. A. Atwater, “Optical pulse propagation in metal nanoparticle chain waveguides,” Phys. Rev. B 67(20), 205402 (2003). 4. J. M. McMahon, A.-I. Henry, K. L. Wustholz, M. J. Natan, R. G. Freeman, R. P. Van Duyne, and G. C. Schatz, “Gold nanoparticle dimer plasmonics: finite element method calculations of the electromagnetic enhancement to surface-enhanced Raman spectroscopy,” Anal. Bioanal. Chem. 394(7), 1819–1825 (2009). 5. C. Tabor, R. Murali, M. Mahmoud, and M. A. El-Sayed, “On the use of plasmonic nanoparticle pairs as a plasmon ruler: the dependence of the near-field dipole plasmon coupling on nanoparticle size and shape,” J. Phys. Chem. A 113(10), 1946–1953 (2009). 6. A. M. Funston, C. Novo, T. J. Davis, and P. Mulvaney, “Plasmon coupling of gold nanorods at short distances and in different geometries,” Nano Lett. 9(4), 1651–1658 (2009). 7. D. Langbein, “Theory of Van der Waals attraction,” in Springer Tracts in Modern Physics (Springer, 1974), pp. 1–139. 8. 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(24), R 16356–R16359, 16359 (2000). 9. 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(9), 1714–1716 (2002). 10. S. A. Maier, M. L. Brongersma, P. G. Kik, and H. A. Atwater, “Observation of near-field coupling in metal nanoparticle chains using far-field polarization spectroscopy,” Phys. Rev. B 65(19), 193408 (2002). 11. I. Romero, J. Aizpurua, G. W. Bryant, and F. J. García De Abajo, “Plasmons in nearly touching metallic nanoparticles: singular response in the limit of touching dimers,” Opt. Express 14(21), 9988–9999 (2006). 12. N. Harris, M. D. Arnold, M. G. Blaber, and M. J. Ford, “Plasmonic resonances of closely coupled gold nanosphere chains,” J. Phys. Chem. C 113(7), 2784–2791 (2009). 13. S. Y. Park, and D. Stroud, “Surface-plasmon dispersion relations in chains of metallic nanoparticles: an exact quasistatic calculation,” Phys. Rev. B 69(12), 125418 (2004). 14. S. L. Zou, and G. C. Schatz, “Theoretical studies of plasmon resonances in one-dimensional nanoparticle chains: narrow lineshapes with tunable widths,” Nanotechnology 17(11), 2813–2820 (2006). #122462 $15.00 USD Received 8 Jan 2010; revised 19 Mar 2010; accepted 23 Mar 2010; published 26 Mar 2010 (C) 2010 OSA 29 March 2010 / Vol. 18, No. 7 / OPTICS EXPRESS 7528 15. F. J. García de Abajo, “Nonlocal effects in the plasmons of strongly interacting nanoparticles, dimers, and waveguides,” J. Phys. Chem. C 112(46), 17983–17987 (2008). 16. M. D. Arnold, and M. G. Blaber, “Optical performance and metallic absorption in nanoplasmonic systems,” Opt. Express 17(5), 3835–3847 (2009). 17. T. Hagihara, Y. Hayashiuchi, and T. Okada, “Photoplastic effects in colored KCl single crystals containing potassium metal colloids. I. Preparation of specimens enriched with potassium metal colloids,” Memoirs of Osaka Kyoiku University. Ser. 3. Natural Science and Applied Science 46, 49–56 (1997). 18. J. H. Weaver, and H. P. R. Frederikse, Optical properties of selected elements, 82 ed. (CRC Press, 2001). 19. D. W. Mackowski, “Calculation of Total Cross-Sections of Multiple-Sphere Clusters,” J. Opt. Soc. Am. A 11(11), 2851–2861 (1994). 20. D. W. Mackowski, and M. I. Mishchenko, “Calculation of the T matrix and the scattering matrix for ensembles of spheres,” J. Opt. Soc. Am. A 13(11), 2266–2278 (1996). 21. M. Ringler, A. Schwemer, M. Wunderlich, A. Nichtl, K. Kürzinger, T. A. Klar, and J. Feldmann, “Shaping emission spectra of fluorescent molecules with single plasmonic nanoresonators,” Phys. Rev. Lett. 100(20), 203002 (2008). 22. R. Fuchs, “Theory of optical properties of ionic-crystal cubes,” Phys. Rev. B 11(4), 1732–1740 (1975). 23. C. F. Bohren, and D. R. Huffman, Absorption and scattering of light by small particles (Wiley, 2004). 24. H. Raether, Surface Plasmons on Smooth and Rough Surfaces and on Gratings, Springer Tracts in Modern Physics (Springer, 1988). 25. M. A. Yurkin, and A. G. Hoekstra, “The discrete dipole approximation: an overview and recent developments,” J. Quant. Spectrosc. Radiat. Transf. 106(1-3), 558–589 (2007). 26. V. A. Markel, “Coupled-dipole approach to scattering of light from a one-dimensional periodic dipole structure,” J. Mod. Opt. 40(11), 2281–2291 (1993). 27. B. Khlebtsov, A. Melnikov, V. Zharov, and N. Khlebtsov, “Absorption and scattering of light by a dimer of metal nanospheres: comparison of dipole and multipole approaches,” Nanotechnology 17(5), 1437–1445 (2006). 28. T. J. Davis, K. C. Vernon, and D. E. Gomez, “Designing plasmonic systems using optical coupling between nanoparticles,” Phys. Rev. B 79(15), 155423 (2009). 29. R. L. Chern, X. X. Liu, and C. C. Chang, “Particle plasmons of metal nanospheres: application of multiple scattering approach,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 76(1), 016609 (2007). 30. W. Y. Chien, and T. Szkopek, “Multiple-multipole simulation of optical nearfields in discrete metal nanosphere assemblies,” Opt. Express 16(3), 1820–1835 (2008). 31. N. A. Nicorovici, R. C. McPhedran, and B. Ke-Da, “Propagation of electromagnetic waves in periodic lattices of spheres: Green’s function and lattice sums,” Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 51(1), 690–702 (1995). 32. M. Abramowitz, and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (Dover, 1964).


Introduction
Clusters of metal spheres exhibit localized surface plasmon (LSP) modes, some of which are red-shifted and increased in strength relative to similar isolated spheres [1].The tunability of the resonances is interesting in its own right, but in addition the associated increase in permittivity increases the strength of electric fields and hence energy transfer by forcing higher field gradients, and may also allow access to regions of the permittivity function where damping is reduced.Isolated particles may be similarly tuned by varying the shape, but it is arguably easier to dynamically tune clusters by changing the spacing between particles.These properties imply several practical applications, including nanometric rulers [2], waveguiding [3], surface enhanced Raman [4], and other plasmonics applications.
In particular, it has been suggested that the mode shifting could be employed as a nanometric ruler [2,5], by relating the position of the mode to the gap between particles.If a universal scaling law exists then interparticle distances can be measured to nanometer resolution by simply measuring the mode shift.In some papers [2,5] it has been suggested that mode shifting due to two arbitrary particles is universally exponential, however other experimental and numerical results for coupled gold rods [6] implied a different relationship, as have other studies.We shall now summarize some relevant investigations of linear chain of n spheres, including the special bisphere or dimer case where n = 2. Figure 1a shows an example of the type of geometry to be considered.
The various studies can be categorized by considering the length scales involved, which are the particle radius a, the lattice spacing d, and the wavelength λ.Using the wave-vector k = 2π/λ, we can write three important dimensionless quantities, ka, kd, and the "gap fraction" f = d / (2a) −1.Often ka<<1, meaning that the particles themselves are only weakly affected by retardation effects.Similarly in the majority of studies we consider kd<<1, so that interparticle retardation and lattice diffraction are not very significant.However the gap fraction, f, is of considerable importance here, as there is a significant difference between the weak coupling (f>>1) and strong coupling (f<<1) regimes.This can be compared to van der Waals interactions for bispheres [7], where the approximate limits for the non-retarded potential go as 1/f 6 at long range, but 1/f at short range.As these two scalings are very different and many applications operate near the cross over, any universal law should consider the two extremes.
A mode shift of order 1/f 3 characteristic of dipole-dipole coupling has been demonstrated theoretically and experimentally for relatively well-spaced chains [8][9][10].In contrast, numerical simulations imply that in addition to very large fields in the gap, closely-coupled gold bispheres exhibit a 1/f 1/2 shifting of the eigenvalue related to permittivity [11].Recently we calculated that the fractional wavelength shift of the longitudinal mode for a moderately closely-coupled chain, relative to an isolated gold sphere, is reasonably approximated by an inverse scaling [12] as a function of the gap fraction, f. Figure 1b shows an example of how the mode moves with f.Our calculations also demonstrated the dependence of mode position with number of particles in the chain, exhibiting a near-exponential mode shift with chain length.Furthermore, the decay constant for this exponential was of the order of two particles giving an indication of the localized nature of the mode.Park and Stroud considered the quasistatic dispersion of infinite chain modes at a range of spacings by first calculating universal depolarization eigenvalues and then relating those to frequency via the Drude model [13].In this example a dipole near 400nm dominates -in general we ignore higher order modes such as the quadrupole seen here emerging near 350nm.Varying the gap fraction moves the peak, as indicated by dots.In the remainder of this paper, we only plot the dipole peaks.
Our aim in this article is to determine a universal scaling for the resonances of closelycoupled chains of arbitrary length and spacing, for a number of real metals.Calculations of closely-coupled resonances are reasonably intensive due to singularity on contact [11], so it is advantageous to separate universal geometric factors from material permittivity.We will present a model to facilitate optimization of factors relating to resonance strength such as the extinction and electric field enhancement.In this article we focus on the longitudinal dipole resonance because it exhibits the largest changes compared to other polarizations and mode shapes [12].We limit our search to spheres that are small compared to the wavelength, because larger radii lead to radiative damping and diffractive resonances [14] that complicate analysis.Conversely the spheres cannot be so small that surface scattering of the electrons is an issue [1], as this would require significant modification of the permittivity.Likewise, some care is required in assessing the validity of bulk permittivity at very small separations [15].
It is of considerable practical importance to be able to predict the performance of real metals.We will explore a range of permittivity values by considering the metals Au, Ag, Al and K.These particular metals are chosen because previous studies have shown that they all have low damping for some regions, and they have quite different permittivity profiles [16].Au is used extensively in plasmonics due to its chemical inertness, and yet it has very strong interband damping around the sphere resonance, requiring modes to be red-shifted to achieve reasonable performance.Ag is also used, especially for superlensing as it has lower damping at low permittivities.Al is unique in having a very short plasma wavelength with performance that declines with red-shifting.K has a long plasma wavelength and like the other alkali metals offers excellent performance over a wide range.Obviously it is less physically amenable due to its reactivity, but this may be acceptable in some applications [17].We used dielectric tables collated in a single review [18], and we note that although there is considerable variation (notably for Ag) among different authors, our main results are universal.We have only produced results for a vacuum medium, but we do include the medium permittivity in the theoretical discussion.
Our general strategy is for each combination of metal, gap-fraction f and number of particles n, find the peak of extinction with respect to wavelength, as shown in Fig. 1b.The peak is then used to extract useful parameters from this data, including the permittivity and the mode strength.At the same wavelength we also calculate the maximum electric field E max , on the surface near the centre of the chain, as indicated in Fig. 1a.All results were calculated using fully-retarded vector spherical harmonics [19][20][21], and checked using finite element simulations of infinite chains [12].Mode parameters were compared to those predicted from the eigenmodes of a static boundary element method [22].We adapted this method to give the field enhancement parameter directly and will show how to implement it efficiently for chains of identical spheres.
Before we detail our search for a universal scaling, it is useful to compare the different metals at different wavelengths, as seen in Fig. 2 where we show the peaks that result from changing the particle separation.We only show results for the bisphere (n = 2) and infinite chains here -later results prove that these are the limiting cases.We can compare the extinction and field enhancement of the different metals in this geometry, and the results are largely consistent with that expected for simpler geometries when comparing different wavelengths.Despite the rather poor performance of gold for widely spaced particles, it can produce larger extinction other metals in the near-infrared due to lower loss than silver and stronger depolarization than potassium at those frequencies (see Fig. 2a).Field enhancement is not particularly sensitive to the material -all strongly coupled chains have high enhancements, as seen in Fig. 2b.In general the electrostatic results agree quite well with the electrodynamics for the small spheres used here, except for an apparent wavelength shift for field enhancements on Au, and some damping for infinite chains which is explained by radiative effects.It is important to understand that actual experimental results could differ quite substantially from these due to differences in permittivities and non-continuum effects as discussed earlier.(Color online).Summary of peak behavior as a function of wavelength due to changes in particle separation.Peak extinction (a), field enhancement (b), particle separation (c) and the associated mode permittivity shift (d) compared to that of an isolated sphere.Electrostatic limits (lines) are shown for the infinite chain (dashed) and bisphere (solid).The lines are limited in extent by mode behavior at short wavelengths and calculation limits at long wavelengths.Also included are VSH-Mackowski extinctions for bispheres (dots), VSH-GMM field enhancements for bispheres (crosses), and FEM results for infinite chains (circles), all damped by radiation.Details of methods can be seen in later sections.
To understand these results we also show how the chains could be used as a ruler, by plotting the gap-fractions corresponding to a range of peak wavelengths.Figure 2c shows that there is a substantial shifting of the peak at small gap-fractions.Next we will show how the dipole model allows most of the data to be collapsed onto universal curves, but to achieve this comparison of different metals it is necessary to compare the modes of the same permittivity.Specifically we compare to the permittivity shift relative to that of isolated spheres -we have plotted this function for reference in Fig. 2d.This may seem unusual but actually it is the most logical comparison since permittivity relates directly to the nature of the resonance [1,23], and we previously employed this strategy for a wide variety of plasmonic systems [16].In the following sections we first provide a review of a dipole approximation and its application to the problem at hand.This serves as a useful model since the measured resonances are dipolelike in the far-field.We briefly describe the limiting case of dipole-dipole coupling at large separations, and how to correct for radiative damping.Then we detail the numerical methods we have employed followed by results for chains of varying length and gap fraction in search of a universal law, and finally present some concluding remarks.

Generalized quasistatic dipole model
It is useful to consider the general properties of LSP of a small system that may in fact be a cluster of individual particles.The system response at an arbitrary frequency can be described by the polarizability α, and a convenient model in the electrostatic limit, assuming the system size is much smaller than the wavelength, is [23]: where V is the volume of metal in the system, ε is the frequency-dependent permittivity of the metal, ε b is the permittivity of the background, L is the depolarization that affects the position of the resonance, and A is a strength factor.L & A are independent of the metal but are geometry dependent.The strength factor was introduced by Fuchs [22], where it was shown that in a multimode system the strength factors add to unity.We omit the sum because we are mainly interested in the dominant dipole mode which is often unaffected by near-by modes, at least for low-loss metals.
The polarizability α can be used to predict useful quantities in the quasistatic (small-size) limit, including the extinction cross-section, and scattering cross-section ( ( Generally, when considering small metal particles, absorption C a = C e -C s dominates extraction of energy from the resonance, and scattering is small.Note that Eq. ( 3) is correct for dipole fields of a small system, but some care is required in applying it to a cluster such as a chain, since the cluster is not necessarily "small".We will discuss this point further in the next section.
The polarizability is also related to the near field enhancement of particles [24], that is the ratio of the surface electric field E to the incident field E i : ( ) Here we have introduced a new coefficient B, which is a property of the mode fields at a given observation point, and A is essentially a dipole moment of B. The coefficients mentioned have well-known trivial values for simple geometries, and these serve as useful limiting cases.For example, L = 1/3 for an isolated sphere, whereas for any small isolated ellipsoid L varies depending on shape but A = 1 and B = 1.However, in the case of a chain of spheres all three parameters vary and must be determined numerically, which we present in our results.We now determine the position and strength of the peaks associated with the dipole mode, since this is usually of most interest in a practical sense.First we must determine the peak position, which occurs at a particular permittivity ε = ε' +iε′′ .Note that, in contrast to diffractive resonances, the LSP resonance is only indirectly related to wavelength via the permittivity function of the metal.Assuming ε ε ′ ′′ − >> , the ellipsoid peak occurs when the real part of the permittivity approaches the pole at ( ) and at the exact maximum of C e we find .
Other quantities such as scattering and field enhancement have slightly different peak positions but in the quasistatic limit they align.
Of particular interest in the present work is the pole shift relative to that of an isolated sphere ( s ε ), and we express this quantity as 1.
In general we find Eq.( 7) useful because it partially compensates for finite size effects, and is more practical than the related depolarization used in Eq. ( 5).It is worth connecting the various descriptions of mode shift for sphere chains in the weak (L = 1/3) and close ( 0 L → ) coupling limits, including the Drude model for low-loss metals, which implies Hence depolarization and permittivity shifts are inversely related across a wide range, but the relationship to wavelength shift changes.However it is important to note that real metals have more complicated wavelength relationships which can deviate significantly from the Drude model, especially in the weak coupling region.Hence we must emphasize that wavelength is not particularly useful here and permittivity is better for determining universal scalings.
Substituting the peak described by Eq. ( 6) back into Eq.( 1), and assuming ε ε ′ ′′ − >> , the strength of the strongest possible dipole at a given metal permittivity is approximately Equation ( 8) can be used to estimate the optimum peak extinction and field enhancements if A & B are known, or vice-versa.Finally, it has been shown that the peak sharpness of a single resonance in the limit of zero scattering is absorption limited.This can be quantified by considering the frequency dependence of the resonance quality factor Q (centre frequency to bandwidth ratio), which at low damping is given by [1]: Note that the peak sharpness is only indirectly related to geometric factors via the resonance permittivity.We now consider some corrections to the model and what sort of predictions can be made by coupled dipole models.

Radiative corrections and weak coupling models
It is important to consider the effect of radiative damping, because all systems of finite size are affected.Radiative damping in dipoles may be approximated through energy considerations which give the so-called radiative reaction (RR) correction to the polarizability [25]: This expression gives a good approximation to the peak height for single spheres, even when scattering dominates.Unfortunately it grossly underestimates radiative pole shift, which can be a significant problem when the permittivity changes rapidly, as it does in real materials.
Equation ( 11) should be interpreted to mean that the non-retarded extinction e C can be estimated from the equivalent retarded quantity by a correction involving the ratio of the retarded extinction RR e C and absorption RR a C , provided that the comparison is made at the peak.In the results section we use this equation to make a fair comparison of retarded calculations to the "universal" non-retarded curves.It is also reasonable to expect that a similar damping applies to the Q of the resonance and to the electric enhancement, whereas the scattering cross-section goes like Equation ( 11) and ( 12) are useful because they are expected to apply to clusters in general, provided retardation across the cluster is not too great.
On the other hand, it is instructive to consider the specific case of clusters with weak coupling, to find a limiting case for the mode permittivity.Chains with small particles that are weakly coupled can be adequately modeled [14] using the coupled dipole approximation.Note especially that because the chains that we consider in this section are well separated and composed of very small particles they are only weakly coupled by radiation, so interference effects are negligible.Most of the dipoles (except those at the ends) are observed to be similar in strength [12].Thus it is reasonable to consider average per-particle polarizability so that chains of different length can be easily compared.The weak coupling limit of polarizability per-particle is given by: where α 1 is the isolated-particle polarizability, and S n is the lattice sum of the free-space Green's function.This approach has been covered quite extensively in the literature, as summarized in Zou [14], so we only discuss the main points here.In systems with many large particles very narrow diffractive (global) resonances may be observed [14], but they will not feature in this article because they involve strongly retarded interactions that we wish to avoid.Using quasistatic limits it can be shown that the lattice sum required for Eq. ( 13) is given by ( ) where p is a polarization dependent factor which is 2 for longitudinal and −1 for transverse, and s n is particle-number dependent and rises from 1 for bispheres up to 2ζ(3) ≈2.4 for infinite chains [26].Equation ( 14) implies the permittivity shift is approximated by 6 3 n πa S .In general, the coupled dipole model is not fully applicable to close clusters because it fails to account for higher multipoles that dominate the interaction [27], so we must consider more general methods.

Methods applicable to strong coupling
Numerical methods are required to account for the near-field interactions, and these all have advantages and disadvantages.Some studies have used arbitrary discretization techniques such as boundary elements [11], finite elements [12] and a discrete dipole approximation (DDA) [2], but these methods usually require very significant computational resources.DDA is known to have particularly poor convergence for the system under study [12], due to sensitivity to surface discretization.The geometry in question is relatively amenable to spherical harmonic expansions, either electrostatically [13] or fully retarded vector spherical harmonics (VSH) [12], however the mathematics of this approach is quite complicated.All methods struggle near contact due to large field gradients [11].In this article we have drawn on the strengths of a surface integral method and VSH, along with a commercial implementation of the finite element method, and we now detail those methods.
One approach we have used is to extend an electrostatic surface integral [22] to directly calculate mode parameters L, A and B. While some authors have approached coupled systems from the point of view of coupling eigenmodes of arbitrary individual particles [28], our spherical particles have rather trivial eigenmodes (see VSH below) and we have instead calculated eigenmodes of the system as a whole.Our implementation is described in detail in the next section.The main accuracy concerns are sufficient sampling of the integral, correctly integrating the singularities, and calculating the lattice sum for the infinite chain.We were able to calculate some results sampling up to 3000 points in the integral, sufficient to converge results to much better than 10 −3 relative error, although calculations with less favorable particle counts limited those calculations to sparser sampling and hence reduced accuracy.Surface integral results are presented later, in parallel with other results for real metals.
We also used Mackowski's [19,20] well-tested public-domain implementation of coupled vector spherical harmonic (VSH) expansion -also known as the T-matrix method -to give a more realistic result incorporating the effects of retardation.The implementation does not give mode parameters directly -instead we estimated them from the cross-sections using equations developed in the model section.Despite the numerical advantages of VSH, the presence of nearly singular fields means that a large number of multipoles N is required to correctly predict closely-coupled longitudinal modes [27].In this work we used an adaptive convergence test on N that required the relative absorption cross-section to converge to within at least 10 −3 for at least two consecutive values of N. It is also possible to calculate electric fields using VSH [29,30], although this capability is not provided in Mackowski's code, so we used the publicly available GMM-FIELD [21], with a slight correction to the field normalization factors.We used a similar convergence criterion, and found that in general the fields converge more slowly than the cross-sections.We focused on the field on the surface along the chain axis, which is expected to be maximum.Near close-contact the field across the gap is relatively invariant.
As an additional check we also used a commercial implementation of the finite element method, for which the details are published elsewhere [12].Again, symmetry was exploited wherever possible, which was especially important given the expensive numerical scaling of this method.The implementation incorporates retardation and gives electric fields directly with cross-sections derived via a surface integration.We only calculated the infinite chain results this way, because other geometries are less symmetric and have "open" boundaries parallel to the scattering direction, requiring much greater resources.All implementations (including dipole-dipole) were validated against each other using test cases where at least one of the methods was known to be correct.
We compared the four chosen metals (Au, Ag, Al, K) using all the retarded methods.In our retarded calculations we fixed the sphere diameter (2a) at 15nm for Au, Ag and K, but we reduced the Al diameter to 3nm to reduce retardation at the short wavelengths associated with this material.We used bulk dielectrics for all calculations due to the difficulty of including size-dependent effects for all methods and to avoid complicating the analysis.In isolation the spheres are large enough that electron scattering would be not be too great [1], and are sufficiently small that radiative damping is low (<10% contribution).However, in chains with very small separations the field singularity may induce significant electron scattering which would reduce the peak shift [15], and we also observed a significant increase in radiative scattering for Au, Ag & K.Where possible we tested separations in the range 1/30<d/(2a)-1<1 and chain lengths from 2 to 100 particles, but convergence often became impractically slow at small separations and large particle numbers, especially in regions of low metal loss.In every case we have focused on the most red-shifted (dipole longitudinal) mode, as this is the most fundamental.

Surface Integral Method
Surface integral methods are particularly effective for simple shapes such as the sphere chains presented here, and casting as an eigenvalue problem is a powerful framework for characterizing resonances.Fuchs used this approach to describe resonances of single cubes in the quasistatic limit [22] -here we apply it to sphere chains and elucidate how to effectively calculate the field enhancement coefficient B. We have ignored the effect of the host medium to simplify the notation.
The self-consistent equation for the internal electric field is written where P = χE , 1 4 ε = + πχ , dS is a surface element and E 0 is the applied field.The primed coordinate is the source, and non-primed is observation.Using surface normal fields ( ) ( ) ( ) , the eigenproblem corresponding to Eq. ( 15) can be written where 2π are the eigenvalues of the interaction matrix and , where (i) and (j) represent the indices of the unprimed and primed discretized spaces respectively.The eigenvalues implied by Eq. ( 16) are directly related to permittivity poles 2π .2π Equation ( 18) relates directly to ( 5) and ( 6) and can be used to determine depolarization L.
The strength of the m th mode in the longitudinal (z direction) can be found via which requires calculating all significant eigenvectors, or alternatively finding the complex limit 0 lim , ( ) where double prime indicates the imaginary part.Note that the limit given by Eq. ( 20) can be adversely affected by degenerate modes, so we calculated the eigenvalues directly via Eq.( 19).Typically we calculated all of the eigenvalues because the number required near contact is substantial and it was not much more expensive to calculate all of them.The normal component of the external electric field can be determined from the internal field, which in turn can be written as a superposition of eigenmodes, and by inspection we have derived that Equation ( 21) seems obvious but was not mentioned in the original reference [22].Now we discuss some efficient strategies for evaluating Eq. ( 17) for periodic linear chains of spheres.Although there is a vast amount of detailed literature on lattice sums and Green's functions, including some that applies to sphere arrays (e.g [31].),we could not find any specific reference that could be rapidly adapted to the current framework.Hence we performed all the derivations using obvious strategies and detail the results here.The longitudinal mode is axisymmetric, therefore the problem can be reduced to the dimension along the chain, by calculating the elliptic integral over the azimuthal coordinate.For a single sphere this can be expressed as where θ is the angular coordinate in the longitudinal plane, and K is the elliptic integral as per the Abramowitz convention [32].More generally, the interaction between spheres n sites apart is where the sphere to sphere distance / q r = qd a , q is the integer sphere-sphere separation and K and E are the elliptic integrals with the same argument As with all Green's functions, there is a singularity at r = r', which must be treated carefully in the numerical integration.One way to do this is to do a series expansion of Eq. ( 22) about the singularity, and then analytically integrate over the source coordinate: Here the line on the surface has been discretized by angle with N intervals of ∆θ, which is the main limit on convergence near contact.We find N = 360 is sufficient to converge f = 0.033 to about 10 −3 , but typically we used between 2 to 8 times as many points depending on machine constraints.
In the general n body problem the matrix becomes nN on a side, greatly increasing computational demands.Due to the periodicity of the chain, the entire matrix can be deduced from the first particle row block (N by nN).It is also useful to note the rotational symmetry of each N by N block.These facts can be used to avoid recalculating the elliptic integrals.
However, generally matrix operations such as diagonalization or inversion require the entire matrix.Nevertheless, the symmetry of the bisphere and infinite chain cases can be exploited, by considering a single observation particle and collapsing source particles on top while accounting for symmetries.In the bisphere case the lattice summation of the Green's function is readily accomplished observing symmetries: and likewise, in the infinite case: ( ) We have inserted the subscript Σ in Eq. ( 26) and ( 27) to distinguish the total interaction matrix from the interactions of individual spheres as given in Eq. ( 22) through (25).Note that the coordinate reversal is achieved through a trivial rearrangement of the non-reversed matrix.Further, calculation of the full lattice sum for the infinite chain can be truncated by splicing in a series expansion of the lattice sum at long-range.We conservatively truncated the full sum at q = 9, and used a third order series approximation to account for more distant neighbors, which reduces the relative norm error of G Σ∞ to about 10 −6 in the worst case.Explicitly, we calculated the summation in (27) by estimating the distant part of the sum as: ( ) ( ) where the series coefficients over the index p are: ( ) Equation ( 28) and ( 29) should be combined with exact near-neighbor terms as given by ( 23) and (24) to complete the estimate of the infinite lattice sum in Eq. (27).
We used this method to calculate electrostatic limits for the universal parameters, and we now present the results in conjunction with those from other methods.

Universal scaling results
Firstly we consider the mode position, as shown in Fig. 3, where we have calculated the shift of the mode permittivity.As mentioned earlier, the mode position is significantly affected by just two parameters -the number of particles n, and the gap fraction f = d/(2a)-1.The effect of these two parameters is substantially independent.Our earlier work demonstrated an approximately exponential decay with particle number [12].The dependence on gap fraction is more complicated, but nevertheless, all materials fall on a similar curve that goes from 1/f 3 at large separations to 1/f 1/2 at very small separations.Any differences between the different metals are probably due to retardation effects associated with the non-zero size of the chain.While it is not obvious from Fig. 3, it is interesting to note that the effect of chain length eventually reduces as the particles come into contact, which is expected as the near-fields dominate the interaction.Secondly we consider the mode sharpness Q, as shown in Fig. 4. The data is reasonably well predicted by Eq. ( 9) when radiative damping is accounted for.Since the quality factor is inversely related to energy loss, it is reasonable to expect that Q e /Q a = C a /C e , as mentioned in Section 3.There are a few significant deviations due to the unusual permittivity of Au near the isolated sphere resonance and the presence of a neighboring mode in the case of Ag.With the outliers explained, we can have confidence that the radiative correction is appropriate.In Fig. 5 we consider the mode strength A, which describes how the polarizability compares to that of an equivalently depolarized ellipsoid.As discussed in Section 3, we correct for radiative damping by using C e 2 /C a per particle for the left hand side of Eq. ( 2), and then we find that the electrodynamic results are reasonably well bound by the electrostatic bisphere and infinite chains results, if the data is plotted as a function of peak shift (Fig. 5b).There are some deviations evident for large gap fractions (Fig. 5a) that may be attributable to interband transitions and neighboring modes.Deviations for infinite chains are probably due to remaining retardation effects.It is interesting that the infinite chain strength reduces least rapidly as the peak shifts, in other words the chain suppresses the generation of higher order (non-radiative) modes compared to the bisphere., which shows the universal scaling of peak extinction.Electrostatic limits are plotted as lines, for both the bisphere (solid) and infinite chain (dashed).Retarded calculations, corrected for scattering by using Ce 2 /Ca on the left hand side of Eq. ( 2), are overlaid using markers.This includes VSH results for all chain lengths (dots), and finite elements for infinite chains only (circles).
We also test how well Eq.( 3) predicts scattering, as shown in Fig. 6, by considering the polarizability of the chain as a collective entity.This approach is reasonable, provided that the chain length is short compared to the wavelength.Widely spaced particles result in a long chain, but this case does not greatly interest us.On the other hand, a large number of particles could also lead to a long chain, however we did not investigate this limit fully, because resources did not allow and to avoid chain-end resonances.Note that for long chains the actual scattering is less than that predicted by (3), because the coupling between particles reduces, and the spheres begin to act individually.For shorter chains, the most significant deviations (~20%) are observed for Au, n = 2, and the closest coupling -this is probably influenced by higher-order modes.The maximum field enhancement coefficient B, on axis just outside the surface of the sphere nearest the middle of the chain, was calculated using the surface integral and compared to retarded results for real metals using Eq.(4). Figure 7 shows reasonable agreement, with very substantial enhancements compared to equivalently depolarized ellipsoids.As with the dipole strength there are some minor variations for infinite chains which could be explained by the effect of neighboring modes or radiation.Interestingly the field is weaker for longer chains at equivalent depolarization, which may be the result of increased symmetry.Electrostatic surface integral limits are shown as lines with infinite chains (dashed) and bisphere (solid).Retarded results corrected for radiative damping are shown as markers.This includes VSH-GMM for all chain lengths (crosses), and finite elements for infinite chains only (circles).

Conclusion
Our calculations applicable across strong and weak coupling regimes demonstrate that while the behavior of the mode position as a function of chain geometry is indeed nearly universal across different materials, there is not a single scaling rule that applies across all coupling regimes.Our calculations certainly indicate that for the metals studied and within the validity of our model the scaling is not exponential with the particle gap, and does lie between 1/f 3 at long range and 1/f 1/2 at close range.The model presented here is separated into geometric factors that require intensive calculation yet are universal (depolarization L, dipole strength A, field coefficient B), and simple functions of the particular metal permittivity used.Provided that radiative damping is considered, the model accurately predicts important parameters like absorption cross-section and electric field enhancement.Extinction is larger for longer chains, but field enhancement is slightly weaker, although this is a small effect compared to the very large enhancements possible in depolarized systems.

Fig. 1 .
Fig. 1. (Color online) An example geometry showing a chain of n = 3 spherical metal particles (a), and the optimization process (b).Only longitudinally polarized excitation is considered.Extinction as a function of wavelength (blue line) is shown for a chain with 3 Ag spheres of radius a = 7.5nm and d = 16.5nmcenter-spacing, giving a gap fraction f = d/(2a)-1 = 0.1 corresponding to the diagram on the left.In this example a dipole near 400nm dominates -in general we ignore higher order modes such as the quadrupole seen here emerging near 350nm.Varying the gap fraction moves the peak, as indicated by dots.In the remainder of this paper, we only plot the dipole peaks.

#
122462 -$15.00USD Received 8 Jan 2010; revised 19 Mar 2010; accepted 23 Mar 2010; published 26 Mar 2010 (C) 2010 OSA Fig. 2. (Color online).Summary of peak behavior as a function of wavelength due to changes in particle separation.Peak extinction (a), field enhancement (b), particle separation (c) and the associated mode permittivity shift (d) compared to that of an isolated sphere.Electrostatic limits (lines) are shown for the infinite chain (dashed) and bisphere (solid).The lines are limited in extent by mode behavior at short wavelengths and calculation limits at long wavelengths.Also included are VSH-Mackowski extinctions for bispheres (dots), VSH-GMM field enhancements for bispheres (crosses), and FEM results for infinite chains (circles), all damped by radiation.Details of methods can be seen in later sections.

#(Fig. 3 .
Fig. 3. (Color online) Shift of the longitudinal dipole mode permittivity, (a) as a function of the gap fraction and (b) as a function of the chain length.The VSH results for finite spheres are indicated with markers, where different colors indicate different metals.Electrostatic limits are shown as lines.The dashed line in (a) is the infinite chain result, and the solid line is the bisphere.In (b) the lower line has a gap fraction d/2a-1 = 1, and the upper is 1/30.

Fig. 4 .
Fig. 4. (Color online) Peak quality factor compared to that predicted by Eq. (9), as a function of the measured contribution of metal absorption to energy extraction.The line shows the expected relationship, and all VSH (Mackowski) results have been included as dots, covering all separations and chain lengths.This partially confirms the veracity of the damping correction used for later results.

Fig. 5 .
Fig. 5. (Color online) Dipole peak strength for different chain lengths n, as a function of gap fraction (a), and as a function of the corresponding peak shift (b), which shows the universal scaling of peak extinction.Electrostatic limits are plotted as lines, for both the bisphere (solid) and infinite chain (dashed).Retarded calculations, corrected for scattering by using Ce 2 /Ca on

Fig. 6 .
Fig. 6. (Color online) Comparison of peak scattering predicted by Eq. (3) as a function of phase length of the chain, where Eq. (3) used cluster α based on (1) and (2) with uncorrected Ce.The curved dotted line is inversely proportional to the x axis.All VSH (Mackowski) results,

Fig. 7 .
Fig. 7. (Color online) Maximum electric field enhancement coefficient, for different lengths of chains, as a function of gap fraction (a), and as a function of the corresponding peak shift (b).Electrostatic surface integral limits are shown as lines with infinite chains (dashed) and bisphere (solid).Retarded results corrected for radiative damping are shown as markers.This includes VSH-GMM for all chain lengths (crosses), and finite elements for infinite chains only (circles).
Nevertheless, we can cast RR in a different form that gives some physical insight as to the expected magnitude of damping.At resonance ( ) ( ) #122462 -$15.00USD Received 8 Jan 2010; revised 19 Mar 2010; accepted 23 Mar 2010; published 26 Mar 2010 (C) 2010 OSA