Characterization of complex plasmonic modes in two-dimensional periodic arrays of metal nanospheres

Two-dimensional periodic arrays of noble metal nanospheres support a variety of optical phenomena, including bound and leaky modes of several types. The scope of this paper is the characterization of the modal dispersion diagrams of planar arrays of silver nanospheres, with the ability to follow individual modal evolutions. The metal spherical nanoparticles are described using the single dipole approximation technique by including all the retarded dynamic field terms. Polarizability of the nanospheres is providedby the Mie theory. Dispersion diagrams for both physical and nonphysical modes are shown for a square lattice of Ag nanospheres for the case of lossless and lossy metal particles, with dipole moments polarized along the x , y , and z directions. Though an array with one set of parameters has been studied, the analysis method and classification are general. The evolution of modes through different Riemann sheets and analysis of guidance and radiation are studied in detail. © 2011 Optical


INTRODUCTION
The plasmonic properties of noble metal nanoparticles (for example, nanospheres composed of gold or silver) have recently been studied and adopted extensively at infrared and optical frequencies because they offer a new range of possible applications, such as the enhancement of Raman spectra [1,2], their use for biosensor devices [3,4], and their potential aptitude to produce artificial magnetism [5][6][7][8][9][10][11]. Furthermore, periodic arrangements of metal nanoparticles present peculiar properties that could be applied to produce enhancement of evanescent fields in subwavelength regions [12][13][14][15][16][17][18], among other innovative applications. All the above possible applications require, as a first step, a good understanding of mode propagation through periodic arrangements of plasmonic nanoparticles in both two dimensions (2D periodic planar arrays), which is the main goal of the present investigation, and one dimension (linear chains, waveguides, 1D periodic). This knowledge is fundamental to fully control and associate determinate characteristics with a particular set of excitation wavelengths.
In recent years, mode propagation and other related phenomena in 2D periodic arrays of nanospheres have been the subject of study by different groups from the analytical, computational, and experimental points of view [19][20][21][22][23][24]. Broader analyses of 2D periodic arrays of both nanospheres and isotropic nanoscatterers have been carried out in [11,[25][26][27][28]. Likewise, mode propagation in 1D periodic arrays of nanospheres has been analyzed in [21,[29][30][31][32][33][34][35], and broader analyses of 1D periodic arrays of isotropic nanoscatterers and radiators in [36][37][38][39][40]. A review of the state of the art of linear and planar arrays of plasmonic nanospheres has recently been presented in [41]. Despite these extensive studies, several difficulties are still encountered in the analysis and, especially, the interpretation and the characterization of mode propagation in periodic plasmonic structures composed of metal nanospheres. In particular, though it is feasible to determine modes in 1D and 2D periodic arrays, information from previously published works does not provide the knowledge of which mode is physical (i.e., that can be excited) and which one is not, and often also when a mode is proper or improper (proper/improper modes decay/grow away from the array, respectively; this concept is clarified in Section 2). Moreover, it also has been problematic to distinguish the evolution of different modes in the so-called dispersion diagram k − β, where β is the wavenumber of the traveling wave and k is the host medium wavenumber (proportional to the frequency). Thus, to the authors' knowledge, a complete physical modal characterization of 2D periodic arrays of metal nanospheres, understanding mode evolution in the dispersion diagram for complex traveling waves, has never been done. The ultimate goal of this work is, then, the characterization of the complex modes in a 2D periodic array of plasmonic nanospheres, with particular focus in the physical waves that can be excited, i.e., launched into the array by a source in proximity of the array, a defect or a truncation. Different methodologies have been followed in the analysis of this kind of problem. Thus, for the study of the optical properties and resonances of a 2D periodic array of noble metal nanospheres on a dielectric slab, the authors of [22] implemented the spherical-wave expansion and solved boundary conditions problems that are applicable to nanoparticles with either real or complex dielectric values. A semianalytical description of traveling waves on 2D and 3D periodic arrays of lossless magnetodielectric spheres has been presented in [19,20]. There, a spherical-wave and scattering-matrix formulation has been adopted, although the analysis was limited to the dispersion relation of traveling waves propagating along the array axes. The authors of [19,20] considered a lossless array of scatterers made of a perfect electric conductor or silver metal, and they provided dispersion relations for the longitudinal and transverse modes. In [20], Mie theory was applied to the computation of the dispersion relation diagram, kd − βd, of 2D arrays of lossless scatterers, providing plots of the dispersion diagram for the longitudinal and transverse waves with respect to the axes of the array for perfect electric conductor, silver, and diamond nanospheres embedded in either vacuum or glass. Modal analysis of 2D periodic systems of nanospheres has also been investigated in [23], where a point-dipole description was used for each nanoparticle after introducing the concept of effective polarizability to obtain the dispersion relation that describes the frequency response of the array with respect to the excitation wavenumber. A limited discussion of leaky waves was also given in that paper.
The most comprehensive review of the computational and analytical methods for modeling 1D, 2D, and 3D periodic arrays of magnetodielectric nanospheres using polylogarithmic functions was done in [21]. In that work, the analytical study provided equations that relate β to the corresponding freespace wavenumber k of the traveling waves with real and complex wavenumbers, including all proper and improper waves (though without providing a distinction between them) in the direction parallel to a specified array axis. Moreover, the authors of [21] clearly stated the difficulty of finding all possible modal solutions, which may result in missing some branches of the dispersion diagram. Indeed, their objective was primarily to display a representative selection of kd − βd diagrams without providing information on the proper/ improper and physical/nonphysical properties of the found modes. As a matter of fact, the physical set of waves is a subset of all the possible waves of the system. The authors of [32] focused instead on some general properties of 1D periodic arrays and on the determination of the required polarizability of each nanoparticle (not necessarily nanosphere) to have a certain dispersion diagram.
In [21,32,38], polylogarithmic functions were been used to provide an analytic continuation of the propagation wavenumber β into the complex domain. However, as previously discussed, that analytical approach, as used in those papers, provides modal solutions without giving any condition to distinguish between physical from nonphysical and proper from improper wave solutions. Moreover, the dispersion diagrams presented in those works are made of a single continuous curve.
In contrast to what has been done previously, in the present work, we are able to follow single mode evolution in such a way that the dispersion diagram is given as a comprehensive superposition of the evolution of all the modes in the system. As shown here, in some cases, varying frequency, an improper mode could become proper or, vice versa, a proper mode could become improper, transitioning from a physical to a nonphysical one, for example. For this purpose, an analytical formulation in the spectral-wavenumber domain is presented, useful to obtain all the possible modes (real and complex) existing in a 2D periodic array of nanospheres. We use a standard Floquet-wave (space harmonic) expansion [42] that makes use of the single dipole approximation [43] (summarized in [29]) and Mie theory [43] to describe the nanosphere polarizability. Comprehensive numerical results of the dispersion diagrams of real/complex modes are provided for transverse and longitudinal polarizations with respect to the array axes. Following the general classification of the modes in the periodic structures summarized in [42,44], we provide a full characterization of the modes in terms of their direction of propagation (forward/backward), in terms of their guidance and radiation properties (bound/leaky), in terms of the position of the wave vector on the Riemann sheets (proper/ improper), and also in terms of their actual physical excitation in the structure (physical or nonphysical modes). In this way, we are capable of identifying the subset of physical modes that are allowed in the array as well as all the subset of nonphysical ones. Given the knowledge of which modes can be excited, we also discuss which one can be used for guiding and radiating purposes. It is noteworthy that this method can as well be applied successfully to 1D periodic arrays of nanoscatterers using a periodic Green's function (GF) that can be extended analytically for a complex β, as that in [45].
This paper is organized as follows. Section 2 introduces the problem studied and the structure under analysis. Section 3 provides the mathematical formulation and assumptions needed to perform the modal study. A discussion on the physical validity of the different modes is given in Section 4. Then, results of the mode analysis in a lossless 2D periodic array with square lattice are shown in Section 5 (without loss of generality, the numerical results are restricted to complex modes traveling along the principal axes). The lossy case is analyzed in Section 6, and the results are compared with those of the lossless case. A discussion regarding the applicability of the lossless model for β − k diagrams in place of the lossy one is briefly provided as well. A summary of descriptions, reasoning, and interpretations of the obtained results is reported in Section 7, regarding the guidance and radiation properties of such a structure. Some brief conclusions are then given in Section 8.

STATEMENT OF THE PROBLEM
In this analysis, a time dependence of the type expð−iωtÞ is assumed and suppressed throughout this paper. Moreover, in the following equations, bold letters refer to vector quantities, a caret on top of a bold letter refers to unit vector quantities, and a bar under a bold letter refers to dyadic quantities.
Our goal in this paper is the study of the plasmonic wave dispersion in a 2D periodic array of metal nanospheres. For this purpose, let us consider a 2D array of metal nanospheres that are periodically located in a homogeneous medium along the x and y directions with periods a and b, respectively, as is shown in Fig. 1. For periodic arrays excited by a plane wave or for a mode in the array, the dipole moment, p mn , of the ðm; nÞth nanosphere is expressed as where p 00 is the dipole moment of the nanosphere located at the origin of the coordinate system; is the Bloch wave vector of the incident wave or of the modal solution, with β and α its real and imaginary vector components and r mn ¼ max þ nbŷ (m, n ¼ 0; AE1; AE2; …) denotes the location of the ðm; nÞth nanosphere. In general, a mode does not necessarily have its phase vector, β, parallel to its attenuation vector, α. Each mode in the structure is represented in terms of an infinite sum of Floquet waves (also called spatial harmonics) as The ðp; qÞth Floquet wave vector of the guided mode is written as with β pq ¼ β þ ð2πp=aÞx þ ð2πq=bÞŷ. Away from the array plane (positioned at z ¼ 0), each Floquet wave behaves as expðik z;pq jzjÞ, where is the z-directed ðp; qÞth Floquet wavenumber and k is the wavenumber of the host medium.
A proper wave has α z;pq > 0; i.e., it decays away from the array. An improper wave has α z;pq < 0; i.e., it grows away from the array. A proper/improper mode has all its wave harmonics proper/improper. Note that the apparent contradiction that an improper mode grows indefinitely does not exist, because each physical improper wave has a bounded domain of existence, as is well explained in [46], Ch. 5.
In order to obtain a careful description of real/complex modes in an array of metal nanospheres, we summarize here a classification to physically characterize the behavior of the various modes existing in the periodic structure. The structure in Fig. 1 supports different guided modes depending on whether the dipole moment of the spheres is parallel or orthogonal to the plane containing the array, and they can be of two types [42,44]. If all the Floquet waves are slow waves (phase velocity is slower than the speed of light), the mode is a surface-wave mode or a bound mode. Bound modes are nonradiating modes: all the Floquet waves decay exponentially away from the structure (α z;pq > 0), being jβ pq j > k (outside the visible region) for all ðp; qÞ indices, so that the only decay is associated with the metal or host losses [42,47,48]. If at least one of the Floquet waves is faster than the speed of light [49] [namely, at least one spatial harmonic is in the visible region, i.e., jβ pq j < k for some ðp; qÞ], the mode is a radiating mode or leaky mode; i.e., at least one of the spatial harmonics is radiating or leaking energy into space (see, for instance, Table 12.1 in [42] for a summary of this classification of complex modes). In uniform structures, this radiation takes place in the forward direction (β pq · α > 0), and the wave grows unbounded for growing jzj. This unbound nature is why this mode is called improper (namely, its wavenumber is found on the improper sheet of the corresponding Riemann surface) [49]. In contrast, in periodic structures, the radiation can take the form of either a forward or backward leaky wave. The forward leaky waves in these structures are always found to be improper, whereas radiation in the form of backward (β pq · α < 0) leaky waves is always proper (i.e., their field is bounded in the vertical direction) [49].
The above discussion has pointed out that radiation in periodic structures can take the form of bound and unbounded waves in the z direction, which opportunely raises the question of the physical existence of these radiating modes and, in general, of all the modes [42]. This very relevant question will be treated in the next sections after the dispersion equation of the structure is derived.

FORMULATION
The electromagnetic modal analysis of the structure shown in Fig. 1 follows the single dipole approximation summarized in [29,43]. This approximation is an effective tool to model periodic collections of nanoscatterers, and gives a very good description of small metal particles [23,50] when r p ≤ d=3, where r p is the radius of the particle and d ¼ minða; bÞ. If the host medium is isotropic, the induced electric dipole of the ðm; nÞth nanosphere is characterized by where E loc mn is the local field at the position of this particle [E loc mn is produced by the incident field plus that scattered by all the nanospheres of the array except the (m; n)th nanosphere itself], and α ee is the polarizability of the nanospheres (isotropic and identical for every spherical nanoparticle). For α ee we use the following expression obtained from the Mie theory [43]: where ε 0 is the free-space permittivity, ε r;h is the relative permittivity of the host medium, n r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ε r;m =ε r;h p is the relative refractive index, with ε r;m being the relative permittivity of the metal nanosphere, s 1 ¼ n r kr p , s 2 ¼ kr p , and ψ 1 ð·Þ, ξ 1 ð·Þ are the Ricatti-Bessel functions [51] [ψ 0 1 ð·Þ and ξ 0 1 ð·Þ are their derivatives with respect to the argument]. The relative permittivity of the metal as a function of frequency is provided by the Drude model: where ε r;∞ is a convenient fitting parameter to match the experimental data, ω p is the plasma frequency, and γ is the Drude damping term. The local field E loc ðr 00 ; k B Þ at the position r 00 is expressed as [29] E loc ðr 00 ; k B Þ ¼ Ğ ∞ ðr 00 ; r 00 ; k B Þ · p 00 þ E inc ðr 00 Þ; ð9Þ where Ğ ∞ is the regularized dyadic GF that accounts for the field contributions produced by all the (m; n) nanospheres except the ðm; nÞ ¼ ð0; 0Þ nanosphere [29]. The computation of this regularized dyadic GF is discussed in Appendix A.
Because we are interested in the modal analysis of the 2D periodic arrays, we assume that E inc ¼ 0. Combining Eqs. (6) and (9), it is then found that the polarization of the ðm; nÞ ¼ ð0; 0Þ nanosphere is given by This last expression leads to the following linear system of equations: where A is a 3 × 3 matrix given by with I being the identity dyadic. The modes in the 2D periodic array of nanospheres considered in Fig. 1 can now be determined from the solutions of the homogeneous matrix of Eq. (11); namely, the modal wavenumbers correspond to the complex zeros of the determinant of Aðk B Þ. Thus, the dispersion relation of the structure can formally be written as Now taking into account that the nanospheres are arranged in the x-y plane, it is possible to write Aðk B Þ as [23] from which the dispersion relation of the structure can finally be expressed as

PHYSICAL EXISTENCE OF THE MODAL SOLUTIONS
First it should be noted that the periodic nature of the array in the x and y directions causes the function det½Aðk B Þ in Eq. (13) to define a Riemann surface with an infinite number of (p; q) branch points given by the (k x ; k y ) values that satisfy k 2 x;p þ k 2 y;q ¼ k 2 , which correspond to If we assume for simplicity that the mode travels along x (i.e., k y ¼ 0), the branch points are as shown in Fig. 2, where k 00 ¼ k and k pq ¼ k − p2π=a − q2π=b, as described by Eq. (14) in [48] or in [42]. Figure 2 shows how these branch points are periodically located in the k x plane. Each (p; q) branch point defines two Riemann sheets: the (p; q) top Riemann sheet, where Imðk z;pq Þ > 0 (the so-called proper sheet), and the (p; q) bottom sheet, where Imðk z;pq Þ < 0 (improper sheet). A detailed study on the nature of the resulting Riemann surface reveals that all the (p; q) proper sheets are actually periodic and that they overlap to give place to just one top (proper) sheet with the expected appearance of different Brillouin zones (BZs) accounting for the periodic nature of the structure under study. However, the improper sheets are not found to be periodic, and no sheet overlapping is observed. Thus, the Riemann surface has one common top (proper) periodic sheet and an infinite number of improper sheets disconnected among themselves. Nevertheless, the multivalued Riemann surface as a whole is expected to present some sort of periodicity congruent with the existence of BZs (and ultimately with the actual periodicity of the array). This feature is preserved by the fact that the part of the improper sheet with respect to the (p; q) branch point, which is exactly within the limits of the (p; q) BZ, is exactly the same for all the branch points. In this way, every BZ has exactly the same structure in the top sheet as well as in the bottom sheet that is accessible through the branch points located in this particular BZ.
The previous knowledge on the nature of the Riemann surface is key to discerning the physical/nonphysical nature of the different modes of the periodic structure. Thus, in order to identify the modes that are physical in the present case, we start with the appropriate definition of the inverse Fourier transform that has to be used to represent any field excited by a realistic point source [42,48,52]. Taking into account the periodic nature of the Riemann surface and the different poles and branch-point singularities, the adequate inversion contour of the Fourier transform should be the thicker gray solid line shown in Fig. 2 (there are more details in [42]). This inversion contour runs entirely on the real axis of the k x complex plane for the case of a lossy medium, but it detours around the realaxis singularities in the case of a lossless medium. Assuming a source at the origin, and an observer along the positive x direction, it is possible to apply a further path deformation toward the positive imaginary half plane. The deformation consists of an infinite periodic set of the steepest descent paths and residues of those poles encountered in the deformation [42,46,52]. Each encircled pole of the periodic set represents a physical mode. The poles associated with physical modes can be located on both the proper and improper sheets of the Riemann surface. Hence, it is the possibility that the physical modes are located on the improper sheets of the Riemann surface that actually justifies and requires a precise knowledge on the nature of the Riemann surface as well as the study of the evolution of the modes in the improper sheets.
Thus for the particular case shown in Fig. 2, inside the first BZ [defined as −π=a < Reðk x Þ < π=a], has been plotted the wavenumbers of two physical modes: a bound mode and a leaky mode. The leaky mode with 0 < Reðk x Þ < k 00 ¼ k is physical if it belongs to the bottom Riemann sheet with respect to the ð0; 0Þ branch point, which would cause the ð0; 0Þth Floquet wave to be a forward improper wave. The bound mode characterized by Reðk x Þ > k 00 , however, is physical if it is located on the upper Riemann sheet (forward proper wave) [42,47]. If we now consider the ð1; 0Þ branch point (which belongs to the p ¼ 1 BZ in Fig. 2), for the improper leaky mode to be physical, it has to correspond to the ð1; 0Þth Floquet wave (namely, it has to belong to the bottom Riemann sheet with respect to the ð1; 0Þ branch point). And for the p ¼ −1 BZ, the physical mode is that corresponding to the ð−1; 0Þth spatial harmonic. Further plots in next sections show the evolution of the wavenumbers of these modes (real and imaginary parts) in the complex k x plane. In summary, it can be said that the physical modes are those whose corresponding poles in the complex k x plane (assuming k y ¼ 0) are captured in the path deformation shown in Fig. 2 (see also [42,47,48,52]), and, therefore, for an observer along the positive x axis, they have Imðk x Þ > 0 [solutions with Imðk x Þ ¼ 0 can be valid in the lossless case]. In the p ¼ 0 BZ, it implies that the physical modes for an observer along the positive x direction are in one of these four categories: (i) improper forward leaky modes in the ð0; 0Þth bottom Riemann sheet with 0 < Reðk x Þ < k, (ii) proper forward bound modes (top Riemann sheet) with Reðk x Þ > k, (iii) proper backward bound modes with Reðk x Þ < −k, and (iv) proper backward leaky modes with −k < Reðk x Þ < 0. When losses are present, the physical validity of the modes can be tracked back to the corresponding lossless case.

MODES IN LOSSLESS STRUCTURES
In this section, we present dispersion curves for a 2D periodic array of silver nanospheres of radius r p ¼ 25 nm in a lossless host medium with ε r;h ¼ 1 arranged in a square lattice with a ¼ b ¼ 73 nm. The parameters in the Drude model for silver permittivity are ε r;∞ ¼ 5, ω p ¼ 1:37 × 10 16 rad=s, and γ ¼ 0 for lossless nanospheres [5,53]. The modal wavenumbers in the array of nanospheres are computed by searching for the complex zeros of Eq. (13), under the assumption that the complex wave vector is parallel to the x axis (namely, k B ¼ k xx , with k x ¼ β þ iα) and that the supported modes have the dipole moment p 00 parallel to one of the Cartesian axes. Because of the periodicity of the structure, we restrict our study to the first BZ, with proper modes (upper Riemann sheet) plotted in solid lines and improper modes (bottom Riemann sheet) in dotted lines. Furthermore, in the absence of losses, there is a four-quadrant complex modal wavenumber symmetry. In other words, if k x is a complex solution, then −k x , k Ã x and −k Ã x are also solutions (where Ã indicates the complex conjugate). If k x is a real solution, we also have a mode that can propagate with the −k x wavenumber. In what follows, our discussion on the physical validity of the modes reaching an observer along the positive x axis will be restricted to those modes whose wavenumber k x has positive or null imaginary part, Imðk x Þ ≥ 0, because otherwise their corresponding poles could not be captured by the integration path deformation shown in Fig. 2 or in [42]. Because of the reciprocity of the structure under study, for an observer along the negative x direction, the corresponding mode with wavenumber −k x will have the same nature.
The evolution with frequency of the real and imaginary parts of the wavenumber k x (namely, the dispersion diagram) for modes with dipole moments polarized along x, y, and z is shown in Figs. 3 and 4. The physical modes for an observer placed along the positive x direction (that is, those that are excited and effectively traveling with either positive or negative phase velocity and reaching the observer) are tagged in Figs. 3 and 4 by circular markers (○) and the nonphysical modes by crossed circular markers (⊗). These markers will also be used in some of the following figures.
For the case of dipoles polarized along x, a nonphysical real improper leaky mode [red dotted curve in Figs. 3(a) and 4(a), β > k) transitions to a physical real proper bound mode [red solid curves in Figs. 3(a) and 4(a)] at the normalized frequency ka=π ≈ 0:38. This forward-bound mode has a growing propagation constant as the frequency increases. Because of the periodicity, at the boundary of the BZ [Reðk x Þ ¼ π=a], this mode meets another related backward-bound mode (which arrives at the boundary of the BZ from the right-side adjacent BZ), and they join together to give rise to a pair of complex conjugate proper modes [blue solid curve in Figs. 3(a) and 4(a)] with phase constant β ¼ Reðk x Þ ¼ π=a and attenuation constant α ¼ AEImðk x Þ (certainly, only the complex mode with α > 0 is physical and bound, as explained in Section 4). Figure 3(a) also shows an improper complex mode (dotted green curve) that is physical when its phase constant is located in the leaky region above the light line in the first quadrant of the complex k x plane (that is, when 0 < β < k, α > 0), until it reaches the light line to transition to a nonphysical improper bound mode (β > k).
The evolution of the modal wavenumbers is more apparent in Fig. 5, which shows the excursion of the wavenumbers in the complex k x plane as the frequency varies (namely, the loci of the k x modal solutions). In this figure it can easily be appreciated that the solutions have a symmetry with respect to the Imðk x Þ ¼ 0 axis (i.e., conjugate values are also solutions in the lossless case). The arrows in this figure, and in other following similar figures, indicate the direction of increasing frequency. The details of how the improper real solution (not visible in Fig. 5) transitions into the real proper one and then gives rise to a pair of complex proper solutions are depicted in Fig. 6: there, the two forward and backward proper real solutions (red) meet at kx ¼ π=a and become two proper complex solutions (blue) when the frequency increases.
The dispersion diagrams of the modes having dipoles oriented along the y direction are shown in Figs. 3(b) and 4(b), and the loci of their wavenumbers are shown in Fig. 7. In the detailed picture shown in Fig. 8, we can observe in the neighborhood of point A (separating two BZs) that two real proper modes (red solid curves) moving in opposite directions with increasing frequency meet at the limit of the BZ to turn into a pair of complex conjugate proper modes with Reðk x Þ ¼ π=a. Each one of this pair of complex proper modes (purple curves in Fig. 8) can be seen to evolve along the limit line of the BZ up to meet another complex proper mode [green solid curve, β ¼ Reðk x Þ ¼ π=a] to give rise at points B and B 0 to two new pairs of complex proper modes (blue solid curves), which start to move horizontally toward the right-/left-hand sides of the aforementioned meeting points.
All the modes shown in Fig. 8 with Imðk x Þ ≥ 0 are physical according to the previous discussion in Section 4, and this explains the tags that appear in Figs. 3(b) and 4(b). In particular, the presence of two pairs of physical complex proper bound modes is noteworthy, plotted with blue solid curves in these figures for ka=π > 0:37. Only two of the four solutions  corresponding to these complex proper bound (more specifically, those with jβj > k and α > 0) are tagged as physical in Fig. 3(b). This pair of bound waves always appear together; i.e., when they are physical, they both exist. These modes transition into proper leaky modes for jβj < k. However, the proper forward leaky solution with β > 0, α > 0 is nonphysical above the light line (when 0 < β < k), whereas the proper leaky solution in the range −k < β < 0, α > 0 is physical and backward. Figures 3(b) and 4(b) also show a real improper mode (pink dotted curve), never physical, approaching the light line. This mode lies behind the proper real (red curve) in Fig. 7. Figures 3(c) and 4(c) show the dispersion diagrams of the modes polarized along z. Two real forward proper modes with Reðk x Þ > 0 (red solid curve and cyan solid curve with negative slope) meet at the normalized frequency ka=π ≈ 0:39 to give rise to a pair of complex conjugate proper modes (blue solid curves). The analogous counterparts, modal curves with negative Reðk x Þ < 0, meet at ka=π ≈ 0:39. The physical modes for Reðk x Þ > 0 and Reðk x Þ < 0 are tagged by (○). In the complex k x plane shown in Fig. 9, it can be observed that this pair of complex conjugate proper modes (blue solid curve) move toward the origin of the complex k x plane as the frequency increases. They are both bound and physical (those with α > 0) when jβj > k, and they always exist in pairs. At a certain frequency [see Fig. 3(c)], they cross the light line and the mode with 0 < β < k becomes a nonphysical forward proper leaky mode, while the one with −k < β < 0 becomes a physi-cal backward proper leaky mode. Increasing the frequency, at the origin in Fig. 9, these modes (blue solid curves) go through the branch cuts (positioned as in Fig. 2) to reach the improper Riemann sheet, where they are identified as improper complex modes (dotted curves in Fig. 9). Up to the frequency where these improper complex modes cross the light line, the improper mode with 0 < β < k is physical. As the frequency increases, the pair of improper complex modes becomes nonphysical and again approaches the real axis. There, they meet together to transition into a pair of nonphysical improper real modes. The details of the above transitions can be observed in Fig. 10.

MODES IN LOSSY STRUCTURES
In this section, we account for metal losses in the array of silver nanospheres previously studied in Section 5 assuming the following Drude damping term [5,53]: γ ¼ 27:3 × 10 12 s −1 . Further comparison between the lossless and lossy cases will make apparent the situations for which the simpler lossless case is a good model for the more realistic lossy case. The presence of metal losses always adds certain imaginary part to the wavenumbers of the previous real modes studied in the lossless case in Section 5. Thus, the wavenumber of all the     Figs. 3(b) and 4(b). At the normalized frequency ka=π ¼ 0:37, two real forward and backward proper modes [solid red curve in Fig. 3(b)] meet at the boundary of the BZ, generating at A a pair of complex conjugate proper modes (purple curve). These modes evolve up and down the BZ limit line until they meet other complex proper modes [solid green curve in Fig. 3(b)] at points B and B 0 . Starting from these points, two new pairs of complex modes [solid blue curve in Fig. 3(b)] move horizontally to the right-and left-hand sides. modes studied in this section are complex, independently of whether they come from real or complex solutions in the lossless case previously considered in Figs. 3 and 4. Furthermore, the presence of losses breaks the four quadrant symmetry. The only symmetry retained in the lossy case comes from reciprocity: if k x is a solution, then −k x is also a solution. Because of losses, the conjugate k Ã x is no longer a solution; however, if losses are not very significant, there can be a solution close to AEk Ã x . Indeed, note that each complex curve in Figs. 3 and 4 becomes a pair of almost superimposed curves in Figs. 11 and 12. The above fact makes the transition or meeting points discussed in the lossless case not appear now.
For x polarization, the dispersion diagrams in Figs. 11(a)  and 12(a) show how a complex improper mode (red dotted curve) becomes a complex proper mode (red solid curve) at ka=π ¼ 0:38. In contrast with the lossless case plotted in Fig. 5, in Fig. 13 it can be seen that this complex proper mode does not meet any other mode in its evolution to the edges of the BZ. This last figure clearly shows that this proper complex mode continues evolving in the first BZ as frequency increases, with its imaginary part continuously growing. Figure 13 also shows the evolution of the improper modes (green dotted curve). As a consequence of the losses, the phase and attenuation constants are slightly different for the two modes in the fourth and first complex quadrants (the modes in the third and second quadrants are the opposite ones with respect to the first and fourth quadrants, respectively). As mentioned above, there is no longer a pair of complex conjugate modes (as happened in the lossless case shown in the green dotted curve in Fig. 5), although this fact is almost unnoticeable in Fig. 13 due to the small amount of losses. For the same reason, the physical validity of the modes in the lossy case follows practically the same pattern as the one shown in Figs. 3(a) and 4(a) for the lossless case. For the sake of clarity, the tags on the physical/nonphysical modes are also shown in Fig. 13.
In the case of y polarization, Figs. 11(b) and 12(b) show the dispersion diagrams of two complex proper modes (solid blue and red curves). This evolution is depicted in the complex wavenumber plane in Fig. 14. Unlike the behaviors plotted in Figs. 7 and 8 for the analogous modes in the lossless case, there is not any splitting point in this case: both proper modes evolve separately. Moreover, in the lossy case, it is clear that there is not symmetry with respect to the real axis; namely, the blue and red curves are not conjugated solutions (however, the blue solution is close to the complex conjugate of the red solution). The tags in Figs. 11(b), 12(b), and 14 clearly show which one of the set of complex solutions is physical. In order to clarify the behavior of these modes, a detail of the evolution in the complex k x plane is shown in Fig. 15, which can be compared with that for the lossless case in Fig. 8. In the neighborhood of point A in Fig. 15, a complex proper mode (red solid curve) travels toward the limit line of the BZ in the forward direction with an attenuation constant α ¼ Imðk x Þ⪆0. Above the normalized frequency ka=π ≈ 0:37 (in the neighborhood of point A), the magnitude of its attenuation constant starts to increase, whereas its phase constant starts to decrease up to reach zero at ka=π ≈ 0:48 (point B). At this normalized frequency, this complex proper mode goes through the branch cut to become an improper complex mode [red dotted curve in Figs. 11(b), 12(b), and 15]. Another complex proper mode (blue solid curve) moves down close to the limit line of the second BZ as the frequency increases. However, as the normalized frequency approaches the value ka=π ≈ 0:37 (neighborhood of point A), the magnitudes of both its phase and attenuation constants start to increase. The evolution of this mode in the second BZ is reflected in the evolution of the complex proper mode plotted in the second quadrant of Fig. 14, which shows that this mode is physical in all the considered frequency ranges. Also, as happened in the lossless case shown in Figs   nonphysical, propagates toward the light line as the frequency increases. This mode appears in Fig. 14 close to the real axis.
In Figs. 11(c) and 12(c) are plotted the dispersion diagrams of the modes for z polarization along with their corresponding physical/nonphysical tags. The apparent multiplicity of modes shown in these figures disappears in Fig. 16, which shows the frequency evolution of these modes in the complex wavenumber plane. This figure also makes more apparent the physical/ nonphysical character of the modes. The complex proper modes (red and blue solid curves) show a similar evolution as that shown in Fig. 14 for the case of y polarization. A clear difference observed in Fig. 16 is that, at ka=π ≈ 0:42, these two complex proper modes go through the branch cut (through either the real or the imaginary axes) to become improper complex modes (dotted red and blue curves). At this normalized frequency, it can be seen that the physical complex backward proper mode coming from the second quadrant (blue solid curve) transitions into a physical complex forward improper leaky mode (blue dotted curve), which eventually becomes nonphysical forward improper bound when Reðk x Þ > k [see Figs. 11(c) and 16]. The other complex proper mode coming from the first quadrant (red solid curve) was already nonphysical before transitioning to a nonphysical complex improper mode in the fourth quadrant.
To conclude this analysis, having both lossless and lossy dispersion diagrams available, it is possible to observe whether the lossless approximation for the computation of β − k diagrams for the lossy case provides a satisfactory result  (a fact that is especially important for complex modes whose physics can hardly be analyzed within the lossy model). Analyzing the x polarization, the lossless approximation could be adopted without making a large error (as can be noticed by comparing the graphs in Figs. 5 and 13). For y polarization, by comparing the graphs in Figs. 7 and 14, this approximation would produce a larger error than the previous case for the proper modes, especially for the region close to the boundary of the BZ. For z polarization, looking at the graphs in Figs. 9 and 16, this approximation would produce the largest error with respect to the other two polarizations, and as such should not be used to approximate the modes in a lossy structure.

ANALYSIS OF GUIDANCE AND RADIATION
A. Physical Bound Modes Physical bound modes travel a long distance in terms of the free-space wavelength when their attenuation constant is small, i.e., when α ≪ k (low decay). As was previously stated, our analysis is restricted to modes along the x axis, such that β and α are both parallel to the x axis, k x ¼ β þ iα. Taking into account the condition of physical existence, and considering the harmonic in the fundamental BZ, a hypothetical observer along the positive x direction would be reached by a bound mode as outlined at the end of Section 4. A more general condition for bound and physical modes is given in [42,44].
In the planar array of nanospheres under study, the bound modes that travel without large decay are discussed next (the description is performed looking at the more realistic lossy case dispersion diagrams shown in Figs. 11 and 12; a similar discussion can be provided for the lossless case in Figs. 3 and 4). For x polarization, the proper forward mode [solid red curve in Figs. 11(a), 12(a), and 13] is bound and physical. As shown in Fig. 12(a), its attenuation constant α is small for a narrow frequency region. For y polarization, the best bound mode is the proper forward mode [solid red curve in Figs. 11(b), 12(b), and 14]; it has a very small imaginary part near the light line for a wide frequency range. For increasing frequency, the imaginary part grows. The proper backward mode (solid blue curve) is also bound, but its attenuation constant is large. For z polarization, the best bound mode is the proper forward mode [solid red curve in Figs. 11(c), 12(c), and 16], because its attenuation constant is small for a wide frequency range. Also the backward proper mode (solid blue curve) is bound, but its attenuation constant is not small. None of the improper modes can be bound.
As a last remark concerning z polarization, we point out that the real mode dispersion analysis carried out in [14,16] shows a mode whose slope transitions from positive to negative. That mode has been associated with the capability of producing superresolution with arrayed nanospheres. The nanospheres in [14,16] were assumed lossless, and only real    modes were shown. The same dispersion analysis has been found in Fig. 3(c), and we have shown that the mode considered in [14,16] is actually made of two distinct modes: a physical forward proper real mode (the red solid curve with β > 0) and a nonphysical backward proper real mode (cyan solid curve with β > 0). Furthermore, a third complex proper mode (solid blue curve with β > 0 and α > 0) is intersecting the previous two in Fig. 3(c) where the slope vanishes. Indeed, we also show that this complex mode in the second quadrant of the complex k x plane is physical and must be included in the explanation of the complex behavior of fields generated by the arrayed nanospheres. The intersection of these three different mode types is well represented in the proper sheet of Figs. 9 and 10. The more realistic case of mode dispersion analysis for an array that includes losses is shown in Figs. 11(c) and 16.

B. Physical Radiating (Leaky) Modes
Physical leaky modes provide very directive radiation when their attenuation constant is small, i.e., α ≪ k. Again, our analysis is restricted to modes along the x axis, such that β and α are both parallel to the x axis, k x ¼ β þ iα. Taking into account the physical condition, and considering the harmonic in the fundamental BZ, a hypothetical observer along the positive x direction would be reached by a radiating mode excited at the origin as outlined at the end of Section 4. Considering the cases represented in Figs. 11-16, the only modes that can provide a directive radiation (for which they should have a low attenuation constant in order to travel a long distance before leaking all of their power) are z polarized and, more specifically, the improper forward mode (dotted blue curve) and the proper backward mode (solid blue curve) in Figs. 11(c), 12(c), and 16. In Fig. 12(c) it can be seen that there is a frequency range where the attenuation constants of these modes are rather small, which means that they can travel a large distance from a source, thus providing a very large "equivalent radiating aperture." Note that, in principle, these modes can also radiate at broadside (namely, orthogonal to the x-y plane), because their phase constants approach the origin [47]. However, as these modes are z polarized, each element has a radiation null along the z direction (the so-called element factor in array theory). This is observed by noticing that the imaginary part, α, tends to vanish as the real part, β, approaches zero. In other words radiation losses vanish as β → 0. Also the x-polarized improper forward mode [dotted green curve in Figs. 11(a), 12(a), and 13] is radiating, but its attenuation constant is very large and therefore cannot provide directive radiation. For y polarization, the only radiating mode is the proper backward (solid blue curve), but it cannot provide directive radiation, because its attenuation constant is rather large.

CONCLUSIONS
In this work we have presented a thorough study of the modes of a 2D periodic array of metal nanospheres at optical frequencies. We have paid special attention to the mathematical continuation of the different modal solutions in their excursions across the different sheets of the Riemann surface defined by the corresponding dispersion equation. A comprehensive study of this mathematical continuation is key to avoid the eventual loss of some solutions and also because it does provide essential information about the physical nature of the modes. Although the reported method allows for the determination of all the possible modes, we have shown only those with a small imaginary part (i.e., low decay when propagating along the array), and we have discussed the nature of the modes from different perspectives: real/complex, proper/improper, physical/nonphysical, forward/backward, and bound/radiating. This is the first time (to our knowledge) such a complete physical characterization in a 2D periodic array of metal nanospheres has been presented showing which modes can actually be excited by a localized source and how a mode can evolve when the frequency is increasing; e.g., it can stop/start being physical, or it can start radiating. The presented analysis is also useful because it shows when lossy nanospheres could be approximated by lossless ones and when that approximation fails to be valid. The above discussion and classification of the modes allowed us to have a very complete knowledge of the characteristics of the different modes. This knowledge is certainly crucial for a further exploration of practical applications of 2D periodic arrays of metal nanospheres.

ACKNOWLEDGMENTS
The work of Francisco Mesa has been supported by the Spanish Ministerio de Ciencia e Innovación, the European Union FEDER funds (projects TEC2010-16948 and Consolider Ingenio 2010 CSD2008-00066), and by the Spanish Junta de Andalucía (project TIC-4595). Ana L. Fructos would like to acknowledge the financial support from the Spanish Junta de Andalucía (mobility grant) during her stay at the University of California at Irvine, under the supervision of Professor F. Capolino.