Effect of small random disorders and imperfections on the performance of arrays of plasmonic nanoparticles

In this paper, we formulate an analytical theory that quantifies the first-order effect of a small random uncontrollable disorder that is due to limitations in the realization of periodic arrays of plasmonic nanoparticles. In particular, we show how the effect of a small disorder may be quantitatively taken into account when evaluating the guidance properties of these otherwise periodic chains, and how the main effect of the small disorder consists of additional radiation losses for the guided mode. Similar quantitative analyses may be extended to the general class of periodic metamaterials, providing an idea of how disorder affects their electromagnetic response, and which types of disorder have the most effect.


Introduction
Exploring the science and applications of photonic crystals and metamaterials formed by periodic collections of inclusions has rapidly grown in interest and importance in recent years (see e.g. [1]). With a few exceptions (see e.g. [2]- [7]), however, in the majority of the theoretical contributions on this subject, the role of disorders and imperfections due to inherent fabrication limitations is neglected, while, with the present state of the art, experimental efforts in metamaterial realization and measurement inevitably involve such problems. The subwavelength scale in which the constituent inclusions of a metamaterial are typically embedded in a host medium usually allows a homogenization procedure, suggesting that the intrinsic deviation from an ideal periodic arrangement may not play a fundamental role in the metamaterial performance. This disorder, however, should be taken into account for a complete and rigorous analysis of the metamaterial response, which may describe more closely the experimental conditions. A clear example of such limitations is seen when dealing with periodic lattices of inclusions and metamaterials: it is well known that an infinite periodic arrangement of lossless inclusions cancels out the scattering losses from each individual particle, ensuring a completely lossless bulk metamaterial [8], but the measured losses in metamaterials, manifested in the imaginary parts of extracted permittivity and permeability, are often higher than expected from purely theoretical considerations. The role of fabrication imperfections and random disorders comes directly into play in this scenario, generating incoherent superposition of scattered waves from the different inclusions, which results in energy decay for a wave traveling through the metamaterial sample. This energy 'loss' has often been interpreted in the measurements as energy absorbed by the bulk sample, even though actually it is mostly related to incoherent scattering, rather than material absorption.
A different but related problem that we have recently analyzed is represented by the analytical solutions for the eigenmodes supported by one-dimensional (1D) or 3D periodic lattices of closely packed plasmonic nanoparticles [9,10], which may support highly confined guided modes and backward-wave propagation, with various potential applications for negativeindex materials and laterally confined waveguides. Also, in this case, the chain periodicity allows total cancelation of the scattering losses, supporting the propagation of low-loss eigenmodes (with loss only due to the material absorption in the particles). However, disorder is inherently associated with the realization of these chains: techniques for manufacturing large arrays of nanoparticles are currently available, but although they can ensure good control on the average dimension of the particles, it may be more challenging to ensure perfect alignment and positioning of each of them on a large scale. Since the assumption of ideally periodic arrangement of the particles has provided us with closed-form solutions for the propagation of such modes [9,10], which may also easily take into account the presence of ohmic absorption, we may ask whether an extension of these analytical results may also model the intrinsic disorder expected in the fabrication process. More broadly, this may be applied to the general analysis of metamaterials, in order to quantitatively model how the deviation from ideal periodic arrangements of the inclusions may affect their potential application in realistic setups and devices.
In the following, we explore these concepts, and in particular, we study how the presence of disorder may be modeled in simple geometries like periodically arranged nanoparticles. In particular, we show how such disorder effects, when small enough, may be equivalently described as the effect of an additional 'loss' in the materials forming the particles, directly 3 related to the statistical variance of the disorder in the system, and may characterize both quantitatively and qualitatively the effects of imperfections in the realization of these metamaterials. These results may be further generalized to generic metamaterial lattices, providing relevant physical insights into the analysis of metamaterials for practical realization.
It should be pointed out that other groups have recently analyzed the problem of disorders in metamaterials, both from theoretical and experimental points of view (see e.g. [2]- [5]), showing how these effects, although of second order, may be relevant in some applications and how they are sometimes underestimated. An extensive amount of work related to these concepts is also available for random disorder in photonic crystals and Anderson localization in random materials [11]- [16]. Here, we concentrate on the problem of periodic arrays of plasmonic nanoparticles in order to derive a quantitative expression that may model this small disorder in closed form, with a first-order perturbation approach. Unlike the other approaches presented in the literature relative to this problem, our analytical results allow us to theoretically quantify and model the small disorder in order to provide simple physical analogies and insights that could be later applied to metamaterial technology. A preliminary portion of these results has been presented in a recent conference talk [7]. We assume in the following an e −iωt time dependence.

Theoretical analysis
Consider an infinite linear array of identical particles, each characterized by an electric polarizability α ee and periodically displaced along the z-axis at distance d from each other (center-to-center distance). As shown in various recent papers on the problem [17]- [23], the guided eigenmodes supported by this periodic arrangement propagate along the array with an e iβz dependence, where β satisfies the following equations: whered ≡ k 0 d,β ≡ β/k 0 ,ᾱ ee ≡ k 3 0 α ee /(6 πε 0 ), k 0 ≡ ω √ ε 0 µ 0 and ε 0 and µ 0 are the background medium's permittivity and permeability, respectively. In equation (1) the nanoparticles forming the array are modeled as polarizable dipoles. The two equations refer to the two orthogonal polarizations for the eigenmodes of this system, i.e. respectively, longitudinal (figure 1(a)) and transverse (figure 1(b)) with respect to the axis of propagation [9]. Equation (1) supports the possibility of no-damping subdiffractive guided propagation in the limit of lossless plasmonic particles, since, as expected, the mathematical problem represented by (1) ensures a real-valued solution for any 1 < β < π/d under the lossless condition Im[ᾱ −1 ee ] = −1 [9]. Therefore, similar to any periodic metamaterial, this geometry in the ideally periodic limit also ensures that, despite the scattering losses from each one of the nanoparticles, the overall behavior of the chain may have zero radiation losses.
The previous equations may fully take into account the general case in which the particles forming the array are lossy. In such a case, complex solutions for β are required, with a corresponding decay in the direction of propagation. The summations in (1) are not strictly convergent in this situation [19], and an analytical continuation employing polylogarithm functions is required, obtaining the closed-form dispersion equations first introduced in [9]: being the polylogarithm function, as defined in [24].
Let us suppose now that each of the particles in the linear array is not necessarily positioned at the designated location z N = Nd, but slightly shifted as z N = Nd + δ N , where the set of real numbers δ N d represents quantitatively an imperfection or disorder in their position (see figure 1). Analogously, intrinsic fabrication limitations may generate uncontrollable small differences among the particles, and each of the polarizabilities involved in the summations (1) may be rewritten as α N = α ee (1 + ϕ N ), with ϕ N being a set of complex numbers.
In order to model the presence of disorder in this problem, we may assume that these imperfections are inherently random in nature and provide small perturbations from the ideal design values. Under these assumptions, it may be reasonably expected that such a quasiperiodic system may still be able to guide a mode with e iβz space dependence, where the disorder is anticipated to produce an analogously small perturbation in the guided wavenumber β. Here, we aim to quantify this variation, and relate it to the expected degree of disorder.
The dispersion equations (1) have been obtained for the two polarizations in the ideally periodic geometry by assuming that the induced dipole moments along the chain are in the form p N = p 0 e iβ Nd , with N being any integer (positive or negative) and p 0 being the dipole moment induced on the particle located at the origin. Each p N produces an electric field at the origin, with amplitude E N = G(Nd) · p N , where G(N d) is the dipolar dyadic Green function. Imposing that the sum of the E N would give rise to a total electric field at the origin, self-sustaining a dipole moment p 0 provides the dispersion equations (1) and (2), which may be solved for the guided wavenumber β once the array geometry is known. The presence of small disorder in the position of the particles has a perturbative effect on the original dispersion equation, which may be analyzed under a first-order approximation as the perturbation to the self-sustaining field at the location of p 0 . In this sense, the small disorder affects the evaluation of the dispersion equations in two distinct ways: the Green function is modified into G(N d + δ N ), due to the variation in the distance of each of the particles from the origin, which causes a variation in the induced field on the particle; moreover, the value of each of the dipole moments in the chain is also varied as p N = p 0 e iβ(N d+δ N ) by the variation in the particle position. The variation in the shape or electromagnetic properties of the particles may also be embedded in the estimation of the induced dipole moment of each particle, which corresponds quantitatively to the expression . Summing all these effects, each term in the summations in (1), therefore, takes the form where N goes from −∞ to ∞ and it takes into account the contribution from each of the particles in the chain. In this modified form, the sums in the dispersion relations cannot be evaluated in any closed form, since each of their terms depend separately on random occurrences of δ N and ϕ N .
Rearranging the terms in the summation and expanding in Taylor series for δ N , due to the assumption of small disorder, we may formally rewrite the dispersion relations as in equation (2) L : where the following expressions for the modified polarizabilities hold: The coefficients c N , d N and f N are the Taylor coefficients of (3). As an example, in the longitudinal polarization they take the form Equations (4) and (5) combine the effects of disorder along the chain in a single compact expression, accurate enough to evaluate the perturbation of the dispersion relations due to the random disorder in the chain in terms of a change in the 'average polarizability' of the particles forming the array, allowing a quantitative evaluation of the effects of the disorder on the guided wavenumber. In other words, the disorder along the chain is 'seen' by the propagating mode as a variation in the polarizability of the particles in an equivalent ideally periodic array. Due to the assumption of small disorder, the summation in (5) is expected to weakly perturb the solution of equation (4), and that is why the approach of summing the fields at the particle location is still applicable. In particular, in the limit of interest in which the particles are lossless or with low loss (allowing low-damping modal propagation along the chain) it is straightforward to show that the perturbation in (5) principally affects the imaginary part of equation (4), which is identically satisfied in the limit of lossless particles. This implies that a reasonably small 6 disorder weakly affects the phase velocity along the chain, but it may affect more considerably its damping coefficient (imaginary part of β).
Assuming that the random quantities δ N and ϕ N are Gaussian distributions with zero expected value and with variances σ 2 δ and σ 2 ϕ , respectively, we may evaluate this expected variation on the effective polarizability of the particles as 'seen' by the mode. The Gaussian hypothesis on the distribution of δ N and ϕ N represents a realistic assumption, since the disorder associated with technological limitations and imperfections may be reasonably assumed to be a stochastic process in which each occurrence of δ N and ϕ N is fairly independent from the others. Due to the zero mean value of δ N and ϕ N and their independence, the dominant term in the summation (5) where we have applied the properties of the polylogarithms to evaluate the summation (5) in closed form. Even though the form of these expressions is quite complex, these solutions provide some interesting insights into the effect of disorder on these plasmonic chain waveguides. It can be seen in (7) that the main perturbation on the effective polarizability is given by the variance of the disorder in the position of the particles, as shown by (7). Random (and independent) variations in the geometry of the particles appear to compensate for each other in the first-order approximation along the chain, not significantly affecting the expected value ᾱ −1 , as was also verified numerically in [4] for 3D metamaterials in a different geometry. The variance of the position disorder, however, which is a measure of the degree of disorder in the system, may significantly perturb the effective 'average' polarizability of the particles as seen by the guided mode, since the scattering losses from each of the particles cannot be completely 'canceled' by the other particles forming an ideally periodic array, producing residual scattering losses. For this reason, in the following sections, we concentrate on the analysis of random disorder in the position of the nanoparticles, assuming for simplicity that they have the same size.
Due to the nature of equation (4), when low-loss particles and low-damping guided modes are considered, the main effect in the perturbation of (7) is in fact expected to be in the imaginary part of (7), which is a measure of the 'losses' experienced by the mode. Similarly to what was done in [9], the properties of the polylogarithm functions allow a direct evaluation of this 7 imaginary part. Using the properties the following exact relations are obtained for the imaginary parts of the expected values of effective polarizabilities in the two polarizations: under the assumption that Im[β] is negligible (low-damping guided modes). This result is formally consistent with the first-order approximation for the effect of disorder in photonic crystal waveguides derived in [14]. We reiterate that equation (9) applies only in the limit of small disorder, for which these results represent small perturbations of the solution in the ideally periodic scenario. In this limit, this result is of great interest for its inherent simplicity. In this context, it is noticed that this perturbative method is consistent with iterative techniques that take into account small disorders in periodic arrays, for which a recent application to periodic metamaterials has been presented in [6]. Although our method would numerically yield consistent results for any type of random disorder (with zero mean value), the advantage of our technique when applied to Gaussian disorder is that it provides a closed-form solution for the level of additional losses caused by the disorder, directly related to its variance, as shown in equation (9). This provides relevant physical insights into the effect of small disorder in periodic structures, as we further discuss in the following.
Without resorting to any further approximation, and considering the fully dynamic interaction in the infinite disordered array of particles, equation (9) indeed states that, independent of the average distance between the particles, the disorder in the positions of particles affects the modal propagation by adding an effectively additional loss, due to the incoherent scattering radiation from the particles, proportional to the variance of the disorder. Similarly to the closed-form solution obtained in the ideally periodic chain for the imaginary part of the dispersion relations (4), which ensured power conservation along the chain [9], and for which Im[ᾱ −1 ee ] = −1, i.e. the particles are lossless and no-damping propagation is obtained, here we obtain another equally simple generalized expression that takes into account the presence of disorder along the chain and adds another form of radiation damping, quantitatively related to the degree of disorder in the chain. 8 We note that the sign of the additional terms in (9) is strictly negative, ensuring that the disorder increases the damping factor of the modes whatever be the geometry of the chain and the nature of the disorder. For passive particles, in fact, Im[ᾱ −1 ee ] = −1 −ᾱ −1 loss , with α −1 loss 0 [9]. Moreover, equation (9) is consistent with the physical expectation that an increase in the variance of the disorder would correspond to an increase in the damping factor of the guided modes, and that a tighter (and slower) propagation along the chain (increase in β) makes the propagation more sensitive to such imperfections. Another interesting feature of this result is associated with the fact that the transverse-polarization configuration, which ensures backward propagation, is slightly more sensitive to the disorder than the longitudinalpolarization propagation. This is consistent with our previous findings [9] that showed how this polarization has a slightly lower bandwidth and more sensitivity to losses than the forward-wave longitudinal propagation. Figure 2, as a first example, shows the variation of Im[β] with the standard deviation σ δ of the position disorder for a chain of silver particles, calculated using (9) and the results in [9]. Longitudinal polarization, which is the most appealing for guiding purposes, is considered here and in the following. Moreover, we concentrate on the positional disorder, which has been proven in the previous section to be the most significant mechanism of radiation loss. The chain is formed by silver spherical nanoparticles with average radius a = 10 nm and center-to-center average distance d = 22 nm. We have evaluated this chain geometry considering experimental values of the silver permittivity available in the literature [25] and considering the inherent absorption and frequency dispersion (black solid line) of silver at optical frequencies. We have also added the red dashed curve, which neglects the absorption in the particles, in order to isolate the effect of disorder on the damping of the mode. The curves in this case have been evaluated at the free-space wavelength λ 0 = 380 nm, which, following the analysis in [9], ensures the best guiding properties of the mode as the highest ratio between real and imaginary parts ofβ. At this wavelength the permittivity of silver is ε Ag = −3.8 + i 0.19 [25].

Numerical examples and validation
As predicted by the previous analysis, increasing the disorder along the chain induces a perturbation ofβ, mainly reflected in its imaginary part. Small deviations from the ideal position of the particles do not considerably perturb the guiding properties of the chain, and the absorption is usually dominated by the inherent ohmic absorption of metallic particles. Still, the sensitivity to disorder is well modeled by the theoretical considerations in this geometry. Figure 3 shows a similar case, but for larger silver particles and correspondingly larger separation between the centers of the particles, highlighting the difference in the effect of disorder with a change in the geometry of the chain. In this case a = 15 nm, d = 33 nm and λ 0 = 390 nm, for which ε Ag = −3.9 + i 0.19.
It is evident that in this second scenario the effect of disorder is somewhat lower. This is due to the fact that in this second scenario the particles are larger and therefore less affected by the inherent ohmic losses and effect of disorder, since the electric field is less concentrated in their volume. Despite the mode being more confined around the chain (higher Re[β]) this second example ensures more robustness to ohmic losses and to disorder, consistent with its inherent robustness reported in [9]. Variation of imaginary (a) and real (b) parts ofβ as a function of the standard deviation of disorder along a chain of silver nanoparticles, considering (solid) and neglecting (dashed) ohmic absorption. The radius of the particles is a = 10 nm and the average distance between any two neighboring particles (center to center) is d = 22 nm at the free-space wavelength λ 0 = 380 nm. Here and in the following, the longitudinal polarization is considered. Figure 4, as a confirmation of the validity of the previous theoretical results, reports the numerical simulations of a linear chain of particles with some disorder in which the location and size of each particle has been perturbed from the original periodic geometry of figure 2. In particular, we have assumed a Gaussian random noise in the position of the particles with standard deviation σ δ = 0.5 nm and a deviation from the ideal radius of each particle a = 10 nm equal to σ a = 0.5 nm. The figure reports the calculated Im[β] for the ideal chain without disorder and for the chain with added disorder, both evaluated using equation (4), corresponding to the black solid line and the red dashed line, respectively. Moreover, we have reported the calculations for a random chain, obtained by applying (1) and (3) and truncating the summation to the first 20 elements (that is, evaluating the eigenwave number considering the contributions from the nearest 40 nanoparticles around the origin). As we already mentioned, the summation in (1) is not strictly convergent when complex propagation constants are considered, due to the fact that the mode amplitude grows exponentially in one direction of the chain. This is consistent with any complex eigenmode solution (leaky waves or lossy geometries), and this problem is usually solved with proper analytical continuation of the series in the complex domain, as in equation (4).
Since now the summation is constituted by random variables, and it cannot be solved in closed form, we have considered its solution obtained after proper truncation, which is a good approximation to the exact analytical continuation derived in (2). In this scenario, where each particle has a random deviation in position and size, the closed-form analytical solution cannot be derived, but the summation truncated to the first 20 elements provide good agreement with our theoretical results that embed the effect of disorder on the closed-form expression (4) and (5), as is evident from the figure. The residual oscillation with frequency, evident in figure 4, is related to the truncation effect, as shown by evaluating equations (1)-(3) in the case of no disorder and of increased losses using equation (4) (black and red dotted lines, respectively). It is evident that the line corresponding to the random chain (dotted blue line) agrees very well with the truncated dispersion relation obtained using the previous theoretical results. We have verified similar agreement when we change the number of elements considered in the truncated series. Figure 5 shows the full-wave numerical simulations, performed with commercial electromagnetic software [26], simulating the chain geometry of figure 2 with (figure 5(a)) and without ( figure 5(b)) considering the disorder in the position of the particles, in the amount of σ δ = 2 nm. These simulations refer to the wavelength λ 0 = 400 nm, for which the silver permittivity is ε Ag = −4.56 + i 0.22 [25]. In particular, the figure reports a snapshot of the instantaneous magnetic field distribution orthogonal to the plane of the chain. It is evident from the figure how the guided mode along the chain, although supported in both cases, decays faster due to the presence of disorder, consistent with the results of the previous theory and with the presence of an effective additional loss factor in the particles forming the chains. In the figure, we have also reported the plot of the longitudinal electric field amplitude sampled on a line parallel to the chain axis at a distance of 80 nm. Despite the sharp variations due to the granularity of the spheres, a clear exponential decay is visible in both cases in the figure, which respectively, corresponding to the black (darker) and red (lighter) lines. As expected, these extracted values of Im[β] are consistent with the values predicted by equation (9), despite the simplified dipolar approximation of our analytical model.

Conclusions
The results presented in this paper confirm that the disorder in plasmonic waveguides and metamaterials may give rise to additional scattering losses that may be quantified as a function of the degree of disorder present in the periodic lattice. The geometry analyzed here, a linear array of plasmonic nanoparticles, has been selected for two reasons: the possibility of a closedform analytical solution, generalizing the results of [9], and the inherent property that this geometry shares with the more general class of periodic metamaterials, i.e. supporting lowloss propagation by the means of cancelation of the radiation losses that an ideally periodic lattice ensures. It has been proven in various experimental setups that subdiffraction optical waveguides, like the one considered here, indeed suffer from the presence of losses, and signals cannot travel over a few wavelengths in realistic scenarios with subwavelength cross-sectional size. The results presented here allow a proper quantification of the effects of small random disorder on the expected propagation distance and on the further increase in damping due to radiation losses. We believe this may be useful for the design and fabrication of these waveguides and of optical metamaterials based on similar concepts. With similar arguments, these results may be extended to the 3D scenario of closely packed nanoparticles forming backward-wave metamaterials [10] and to the more general class of periodic metamaterials and photonic crystals. We may quantify the amount of expected scattering losses in a realistic plasmonic array or metamaterial affected by disorder associated with technological limitations and imperfections, justifying how the measured losses in such setups are often larger than those expected from purely theoretical predictions obtained under the assumptions of ideally periodic configurations. It is interesting to underline that our theoretical results show quantitatively how, to a first-order approximation, small variations in the shape or electromagnetic properties of the inclusions are less important than an analogous disorder in their position, as regards the overall electromagnetic behavior of the metamaterial setup and in particular its damping and absorption properties. To conclude, these results represent an important step in quantifying the radiation and scattering losses induced by disorder in plasmonic waveguides and periodic metamaterials.