Complex bound and leaky modes in chains of plasmonic nanospheres

: Bound and leaky modes with complex wavenumber in chains (linear arrays) of plasmonic nanospheres are characterized for both longitudinal and transverse polarization states (with respect to the array axis). The proposed method allows for the description of each mode evolution when varying frequency. As a consequence, full characterization of the guided modes with complex wavenumber is provided in terms of propagation direction, guidance or radiance, proper or improper, and physical or nonphysical conditions. Each nanosphere is modeled according to the single dipole approximation, and the metal permittivity is described by the Drude model. Modal wavenumbers are obtained by computing the complex zeroes of the homogeneous equation characterizing the field in the one dimensional periodic array. The required periodic Green’s function is analytically continued into the complex wavenumber space by using the Ewald method. Furthermore, a parametric analysis of the mode wavenumbers is performed with respect to the geometrical parameters of the array.


Introduction
Periodic arrays of plasmonic nanospheres have been studied because of their interesting performance in several practical applications. For example, to achieve near-field enhancement and super-resolution at optical wavelengths [1][2][3][4][5]; to achieve directive radiation [6,7]; to perform microscopy and spectroscopy using SERS [8,9]; to design biosensors [10,11], and DNA detectors [12]. However, despite these innovative applications, it is still not clear how to associate determinate characteristics with a particular set of excitation wavelengths. Several difficulties are still encountered in the analysis and interpretation of mode propagation in one dimensional (1D-periodic, linear chains) and two dimensional (2D-periodic) arrays of metallic nanospheres when varying frequency, for example. Based on what it has been published so far it is somehow hard to discern if a mode is physical (excitable) or not. The scope of this paper is the characterization of the bound (non radiating) and leaky (radiating) modes with complex wavenumber in linear chains of plasmonic nanospheres (Fig. 1), accounting for metal losses, as we recently did in [13] for 2D-periodic arrays.
Modes of linear chains of metallic nanospheres, object of this paper, have been studied in [7,[14][15][16][17][18][19][20][21][22][23][24][25][26], assuming either lossless or lossy metal and using the single dipole approximation (SDA), described in [27,28]. In particular, in [14] it has been shown that the dispersion relation k-β (where k is the host medium wavenumber, and β is the wavenumber of the guided wave) for a linear chain of gold nanospheres, obtained from finite-difference time-domain (FDTD) calculations, closely agreed with that predicted by the SDA in case of real wavenumbers. In [15], the authors have calculated the dispersion relation for a linear chain of gold spherical particles in a uniform host in the case of particle dimensions much smaller than the wavelength, including all the multipolar terms within the quasistatic approximation, finding that the dipole approximation was reasonably accurate for a periodicity larger than three times the nanosphere radius, becoming increasingly inaccurate for smaller periodicity. In [16], the SDA with the expression for the nanosphere polarizability taken from Mie theory has been applied to compute the dispersion diagram (the real part of the wavenumber) of lossless chains of gold nanospheres for longitudinal and transverse polarizations, comparing the results with the ones reported in [14,15]. In [7,[18][19][20][21]26], the authors used polylogarithmic functions to analytically continue the validity of the dispersion relation into the complex wavenumber plane, and they also accounted for metal losses. However, conditions to classify proper and improper modes (proper/improper modes are defined as those that vanish/grow getting away from the 1D array axis, and this concept will be better clarified in Sec. 2.2), as well as to distinguish the subset of physical waves allowed in the linear chain from the complete mathematical set of modes, have not been developed yet when using polylogarithmic functions. In particular in [26], the authors derived the 1D dispersion relation by considering all multipole field contributions, showing that in the case of a linear chain of plasmonic nanospheres the dispersion relation is well approximated by the one obtained through the SDA with the expression for the nanosphere polarizability taken from Mie theory, as we do in this paper, assuming small nanospheres with respect to the wavelength. Experimentally, modes have been analyzed in linear chains of gold nanoparticles comparing the dispersion diagram results from measurements (real part of the wavenumber) with corresponding ones obtained from theoretical fully retarded solutions and numerical FDTD simulations [22,23]. In [24], the effect of small random disorder (e.g., due to fabrication limitations) in linear chains of metallic nanoparticles has been analyzed theoretically. Specifically, guidance properties of these aperiodic arrays (in fact, the periodicity has been broken by a small disorder) have been analyzed, resulting into additional radiation losses for the guided mode. A review of the state of the art of studies performed for both 1D-and 2Dperiodic arrays of plasmonic nanospheres, including mode analysis, fabrication possibilities, and possible applications, has been outlined in [25]. To conclude, we want to emphasize that all the aforementioned excellent referenced work did not provide yet a complete characterization and classification of the modes in a linear chain of plasmonic nanospheres, in terms of leaky/bound, proper/improper, forward/backward, including the conditions for their physical existence, which is the goal of the present paper. We want to stress that for a physical mode we mean a mode that can be excited (launched) into the array by a localized source in proximity of the array, or by a defect or by an array truncation. Other modes that are not classified as physical here may, however, be excited by much more complicated source distributions, though this study is not within the scope of this paper. The analysis in this paper is a necessary step towards a thorough understanding of modal description in periodic arrays. Several results and observations provided in this paper have not been provided elsewhere.
The novelty of the approach presented here is that it allows for the tracking and especially for the characterization of the evolution of all the modes varying frequency. Each nanosphere is modeled as a single dipole, by using the SDA [25,27,28] and the metal permittivity is described by the Drude model. The computation of the complex zeroes of the dispersion relation of the 1D-periodic array is performed by using a periodic Green's function (GF), analytically continued into the complex wavenumber space by using the Ewald representation. The computation of the complex roots of complex valued functions is a difficult problem in its own rights. We have (i) implemented the Muller's method [29,30] (which can be used to find complex zeroes of analytic functions) and (ii) used the built-in IMSL function ZANLY for Fortran [31] (which uses Muller's method) for verification purposes. The two obtained solutions are in excellent agreement with each other. The use of the periodic GF in the computation of the complex zeroes will be explained in Sec. 2.1. As a consequence, full characterization of the modes with complex wavenumber in the linear chain is provided, in terms of guidance/radiance, proper/improper, physical/nonphysical, and propagation direction conditions. In this paper, a new numerical procedure that uses the Ewald representation for the periodic dyadic GF for linear arrays has been introduced, based on previous scalar developments [32][33][34]. Our method is alternative to the one used in [18,20,21,26] based on polylogarithms, and future studies will be devoted to relate these two techniques.
The structure of the paper is as follows. The mathematical description for the simulation model adopted to perform the analyses described in this paper is briefly summarized in Sec. 2, including a discussion regarding conditions of guidance, radiance and physical existence of modes and their constituent Floquet waves. Then, the characterization of the complex bound and leaky modes in linear chains of plasmonic nanospheres is reported in Sec. 3, considering nanospheres made of silver, and accounting for metal losses. Though our examples relate to an array with a precise set of representative parameters, the analysis, wave classification and results are general. The analysis of guidance and radiance is reported in Sec. 4.
Parameterization of modes with respect to the geometrical parameters of the array is shown in Sec. 5. Criteria for physical existence (excitability) of modes are presented in Appendix A. The dyadic form of the Ewald representation of the periodic GF evaluation used in this paper is provided in Appendix B. ε . The radius of each nanosphere is r, and the array spatial period is d.

Theoretical background
The monochromatic time harmonic convention ( ) is assumed. 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.

Simulation model
As stated in Sec. 1, according to SDA, each plasmonic nanosphere of the linear chain, at optical frequencies, acts as a single electric dipole and can be described by the induced dipole moment loc ee , n n α = p E (1) where ee α is the electric polarizability of the nanosphere, and loc n E is the local field produced by all the nanospheres of the array except the considered n-th nanosphere plus the external incident field to the array [25,28]. SDA is a good approximation when the metallic nanospheres are used close to their fundamental plasmonic frequency and particle dimensions are much smaller than the wavelength, when 3 d r ≥ ; however, for small period comparable with the spheres' radius (i.e., 2 3 r d r < < ), more accurate results would involve multipole field contributions [27,[35][36][37]. The polarizability of spherical nanoparticles can be expressed, according to Mie theory [27], by the following equation [25,27,28,38,39] ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) where r is the nanosphere radius, is the relative refractive index, and m ε is the relative permittivity of the metal, which is described by the Drude model fitting parameter, p ω is the plasma frequency of the metal (expressed in rad/s) and γ is the damping factor (expressed in 1/s).
Consider now a linear chain of nanospheres as in Fig. 1, for which each nanosphere is placed at positions n z nd = , with 0, 1, 2, n = ± ± … , and d is the spatial period of the chain. Suppose that the array is immersed in a homogeneous background, with relative permittivity h ε , and that the array is then excited by a plane wave or by a quasi-periodic excitation with wavevector Bẑ k = k z [28]. Consequently, each spherical nanoparticle will have a dipole moment equal to (1) by using n = 0. Then, the local field acting on a nanosphere at position 0 0 z = is given by represents the regularized electric field dyadic GF for the phased periodic array of nanospheres, defined as and n G is the dyadic GF of the background medium. Since the array is placed in a homogeneous medium, the GF (field produced at 0 z = by a unit dipole at z nd = ) simply is where I is the identity dyad. Substituting the expression for the local field given in Eq.
Mode analysis in the linear chain is then performed by computing the complex zeroes of the homogeneous linear system in Eq. (6), after having imposed no excitation source (i.e., ). In other words, the complex z k -zeroes of the determinant of ( ) ) using the methods described in Sec. 1. The series inside the GF definition in Eq. (4) is slowly convergent, and it is even divergent for complex values of the wavenumber z k , and some technique must be adopted to allow analytic continuation of the GF into the complex z k plane. The GF representation adopted in this work makes use of Ewald's method for linear arrays, which, besides providing analytic continuation, it also provides series with Gaussian convergence and only a handful of terms are needed [32][33][34]. The dyadic form of the Ewald representation of the periodic GF evaluation used in this paper is provided in Appendix B.

Floquet waves representation
The electric field in the linear chain can be described as ( where d is the spatial period of the periodic structure, 0, 1, 2, n = ± ± … , z all its spatial harmonics proper, whereas an improper mode has at least one spatial harmonic improper. Note that certain leaky waves are improper, and physical, but they never grow indefinitely because they exist only in certain regions of space [41,42] (in other words, their domain of existence is limited). This means that the mode itself is not present in certain regions, though the radiation caused by a leaky mode is present everywhere as explained in [42]).

Physical excitation conditions by a source, a defect or a truncation
Among all the mathematical wavenumber solutions of ( ) only a subset represents physical waves, i.e., waves that can be excited by a localized source, or defect, or array truncation, as described in Appendix A. The physical modes that can be excited in linear chains of plasmonic nanospheres, according to Appendix A and what follows, are summarized in Table 1, similar to what was previously reported for the 2D-periodic case in [41,43], based on the complex wavenumber of the fundamental harmonic (with p = 0). Table 1 applies to the case of chain period shorter than half wavelength, since it implies that only wavenumbers in the fundamental BZ may represent fast waves.
Assume that a point source is near the origin (Fig. 2) and an observer is located along the positive z axis of the array. The field at the observer is represented by a spectral path integral on the real axis as shown in the figure in Appendix A. Deformation onto the steepest descent paths encounters one or more periodic sets of poles that represent the wavenumbers of the physical FWs [41,44,45]. A mode is represented by an infinite set of wavenumbers, according to Eq. (7). Periodic sets of poles that are not captured in the path deformation do not contribute to the final field value and therefore represent modes that are not excited. In order to be captured in the path deformation in Appendix A, a pole representing a wave with complex wavenumber z z z k i β α = + should have 0 z α > . Therefore, considering the modal harmonic in the fundamental BZ (often the one with p = 0), a physical forward wave belongs to the I quadrant of the complex z k plane (i.e., 0 z z β α > ), whereas a physical backward one to the II quadrant (i.e., 0 z z β α < ). Moreover, to be physical, not only the wave needs to have a wavenumber with 0 z α > , but also its wavenumber has to respect the conditions stated in Table 1, which are also valid in the case of an observer along the negative z axis for which physical waves have 0 z α < . Further discussion about integration path and field representation is provided in Appendix A and references therein. Figure 2 shows how each physical wave outlined in Table 1 evolves in the analyzed linear chain, for the case in which the observer is along the positive z-direction. If the wave is proper, as defined in the previous section, it vanishes for ρ → ∞ (as can be observed in Figs. 2(A), 2(B) and 2(D), where the wave decays along the radial direction). If the wave is improper, instead, it grows for ρ → ∞ , as shown in Fig. 2(C). Note that in this case, the leaky wave existence domain is limited in space as described at the end of Sec. 2.2, and therefore it does not violate any physical principle [41,42]. If the wave is forward, the phase travels towards the positive z-direction, as highlighted by the red arrows in Figs. 2(A) and 2(C), and attenuates while it travels. Whereas, if the wave is backward, the phase travels towards negative z-direction, as highlighted by the red arrows in Figs. 2(B) and 2(D), but still attenuates along positive z (for conservation of energy, a wave cannot grow along positive z). If the wave is bound, then no energy is radiated away (as in Figs. 2(A) and 2(B)); radiation phenomena are present when the wave is leaky, as shown in Figs. 2(C) and 2(D) by the black arrows.

Conditions for guidance and radiance
An important part of this work lies in the fact that we show which modes are physical (i.e., that are excitable by a localized source), and consequently employable in practical applications. For this reason, the conditions to discriminate physical bound and radiating (leaky) modes are here stated based on Table 1 and the discussion in Appendix A.
A bound mode is a mode that has all its Floquet wavenumbers outside the "radiating spectral region", i.e., the spectral region for which , β < − and z α is positive (backward wave), i.e., the wavenumber is either in the first or in the second quadrant of the complex z k plane. Improper bound modes cannot be launched by a source in Fig. 2.
A leaky mode is a mode that has at least one Floquet wavenumber within the radiating spectral region, or A guided mode (leaky or bound) travels a long distance when the imaginary part of its wavenumber is small with respect to the host wavenumber, or z k α << (i.e., low decay).
For completeness, we want to emphasize that the higher order Floquet harmonics outside the radiating spectral region (with , z p k β > , 0 p ≠ ) can be either proper or improper non radiating waves, and only the ones represented by poles captured from the path deformation described in Appendix A can be excited by a source in Fig. 2, and thus be physical. Also, as mentioned in Sec. 2.3, in the following sections we restrict our analysis to the case of chain period shorter than half wavelength, where only wavenumbers in the fundamental BZ may represent fast (leaky) waves.

Mode analysis
In this section, complex modes in a linear chain of silver nanospheres embedded in free space (i.e., which are the only ones that can be captured in the path deformation described in Appendix A when observing along the positive z axis. Note that this subset has the same information as the complete set of modes (i.e., a complete basis for the field) in the linear chain because of the symmetry property stating that z k and z k − are both solutions. Figure 3 shows the dispersion diagrams for the real and imaginary parts of the wavenumber of the FW harmonic in the fundamental BZ of the modes with T-pol and relatively small imaginary part. Figure 4 shows the wavenumber evolutions in the complex z k plane varying frequency, using two different normalizations. We want to clarify that the reported curves in Figs. 3 and 4 are relative to both physical and nonphysical solutions (because together they form a complete basis for the field); physical ones are then distinguished by using Table 1.

Transverse polarization (T-pol)
Modes 'Proper 1', 'Proper 2', 'Improper 1' and 'Improper 2' in Figs. 3 and 4 are computed using the method described in [33], where only m = 0 and −1 are considered, and adapted to the present case as in Appendix B. When the wavenumber dispersion curves change complex quadrant, they can be continued with dispersion curves computed with different m-values (see appendix B). We also show the solution 'Improper 3′, computed with m = +1, as a continuation of the m = 0 red curve into the II and IV quadrants (Fig. 4), which also represents the wavenumber close to the complex conjugate of the 'Improper 1', according to the symmetry properties of the modal solutions described in Sec. 3. Though these waves computed with m values different from 0 and −1 are correct mathematical solutions, their physical validity requires further investigation, as explained in Appendix B. For an observer along the positive z-direction, according to Table 1  whereas the mode 'Improper 2' is always nonphysical, being in the II quadrant (remember that improper modes in the II quadrant cannot be captured, according to Appendix A). From the point of view of modal wavenumber evolution when increasing frequency, the wavenumber of the blue mode 'Proper 1' goes from the II to the I quadrant of the complex z k plane (or from the IV to the III due to the symmetry of the modal wavenumber solutions with respect to the origin stated in Sec. 3, i.e., z k and z k − are both solutions, even in lossy cases), thus it crosses the imaginary axis (exhibiting its lowest z α before this transition) and so a branch cut positioned there (as shown in Appendix A). Therefore, this modal wavenumber passes from the top to the bottom Riemann sheet, i.e., it becomes the one of the blue mode 'Improper 1'. This transition is well observed in Fig. 4. It is noteworthy that due to this change in the Riemann sheet of the modal wavenumber, the physical backward leaky mode 'Proper 1' (II quadrant), when frequency increases, transitions into the physical forward leaky mode 'Improper 1' (I quadrant). The wavenumber of the red mode 'Proper 2' is almost over the light line z k β = at low frequencies. Increasing frequency, it departs from the light line in the bound region until its trajectory is inverted (Fig. 4), and then it goes from the I to the II quadrant (or from the III to the IV due to symmetry), and thus, similarly to what described for the mode 'Proper 1', crosses a branch cut, transitioning into the wavenumber of the mode 'Improper 3′. With smaller losses, the wavenumber of the mode 'Proper 2' would be able to depart farther away from the light line in the I quadrant, reaching large values of z β , and thus causing short plasmonic wavelengths.  Figure 5 shows the dispersion diagrams for the real and imaginary parts of the wavenumber of the FW harmonic in the fundamental BZ of the modes with L-pol and relatively small imaginary part. Figure 6 shows the wavenumber evolutions in the complex z k plane varying frequency, using two different normalizations. We want to clarify again that the reported curves in Figs. 5 and 6 are relative to both physical and nonphysical solutions; physical ones are then distinguished by using Table 1. Modes 'Proper 1' and 'Improper 1' are computed using the method described in [33], and adapted to the present case as in Appendix B, where only m = 0 and −1 are considered. When the wavenumber dispersion curves change complex quadrant, they can be continued with different m-values (see appendix B). We also show the curve related to the 'Improper 2' mode obtained with m = +1, which represents the continuation of the 'Proper 1' solution into the complex z k plane quadrants II and IV (Fig. 6). This wavenumber also represents the wavenumber close to the complex conjugate solution symmetric to the 'Improper 1' solution, according to the symmetry properties described in Sec. 3. The physical validity of the waves found with m-values different than 0 and −1, requires further investigation, as explained in Appendix B. For an observer along the positive z-direction, according to Table 1 and Sec. 2.4, and looking at quadrants I and II in the upper half part of the complex z k plane in Fig. 6, the mode 'Proper 1' is a physical forward bound mode in the I quadrant when z k β > ; similarly, the mode 'Improper 1' is a physical forward leaky mode in the I quadrant when 0 z k β < < . The mode 'Improper 1' in the II quadrant is always nonphysical. From the point of view of modal wavenumber evolution, at high frequencies, the wavenumber of the blue mode 'Proper 1' is in the bound region in the I quadrant, as in Fig.  6(b). Decreasing frequency, the real part of its wavenumber increases while its imaginary part decreases up to about / 0.45 kd π ≈ , then both the real and imaginary parts of this modal wavenumber decrease. Around / 0.39 kd π ≈ , this wavenumber enters in the radiating spectral region (i.e., z k β < ). There, its imaginary part transitions from positive to negative values (see the inset in Fig. 6). This transition translates into the wavenumber passing from the I to the IV quadrant of the complex z k plane (or from the III to the II due to symmetry), thus crossing the

( )
Re z k axis and so a branch cut located parallel to the

( )
Re z k axis (in Appendix A, this branch cut is slightly above the axis when small ambient losses are present; it would be exactly on the axis in absence of ambient losses), and it changes Riemann sheet, i.e., it becomes the wavenumber of the blue mode 'Improper 2'. This transition is well observed in Fig. 6. The wavenumber of the forward green mode 'Improper 1' starts in the I quadrant of the complex plane in Fig. 6. Increasing frequency, its imaginary part goes to zero while its real part remains under the light line (i.e., z k β > ) as in the inset in Fig. 6, thus it does not encounter any branch cut, as can be seen in Appendix A, and remains improper with negative growing imaginary part (this modal wavenumber transitions from the I to the IV quadrant, or from the III to the II due to symmetry).
Notice for example that in [21] modes 'Improper 1' and 'Proper 1' have been described as the same mode, discontinuous for a small frequency range; whereas in [26] they have been described as one continuous longitudinal mode.

Bound modes for transverse and longitudinal polarizations
Accounting for the presence of losses in the silver nanospheres, the only bound mode for Tpol that travels without large decay is the mode 'Proper 2' in Figs. 3(a) and 3(b). This mode is bound, forward and physical (as described in Sec. 3 ). In the entire domain of its physical existence, this mode can be guided very well, in the sense that it can travel long distances from the excitation source. The other physical bound mode for T-pol is the 'Proper 1' that has the imaginary part of its wavenumber which is large in comparison to the host wavenumber k . This mode is bound, backward and physical; it presents first decreasing imaginary part while frequency increases, a minimum at / 0.4 , and then increasing imaginary part. The only bound mode for L-pol that travels without large decay is the mode 'Proper 1' in Figs. 5(a) and 5(b). This mode is bound, forward and physical (as described in Sec. 3.2); it presents growing imaginary part for increasing frequency, and can be excited in all its domain of physical existence, although it would be slowly decaying while traveling only when it has very small attenuation constant (e.g., near the light line at about / 0.39

Leaky modes for transverse and longitudinal polarizations
The radiating leaky modes for T-pol are described next. The physical modes found in the array of silver nanospheres here analyzed radiate, but they cannot provide very directive radiation since the imaginary part of their wavenumber is large with respect to the host wavenumber k [48]. In Figs. 3 and 4, the mode 'Proper 1' is a leaky backward mode, whereas the mode 'Improper 1' is a leaky forward mode. We will show in Sec. 5 the possibility of modifying the imaginary part of the wavenumber of these two modes by varying geometrical parameters of the array. The only radiating leaky mode for L-pol is the mode 'Improper 1' in Figs. 5 and 6. This mode is leaky and forward. However, notice that this mode lies very close to the light line; since the dipoles approximating the nanospheres are longitudinally polarized (i.e., along the zdirection), their far-field radiation pattern would have the classical sinθ shape, with a farfield null corresponding to the direction along which the dipoles are aligned. Indeed, in the case of the analyzed linear chain, when z k β = the radiation should take place along the zdirection and this is not possible due to the aforementioned null, whereas when z k β < for z k β ≈ , a very small radiation (still due to the dipole sinθ radiation pattern) could take place.
A detailed analysis of how to adopt these leaky modes to obtain directive radiation at optical frequencies has been recently shown in [7].

Parametric analysis of bound and leaky modes
In this section, we analyze the evolution of bound and leaky modes by varying the radius r of each nanosphere and the spatial period d of the array.  Fig. 7. The linear chain under analysis here is the same as the one in Sec. 3 where we vary one parameter at the time. In Fig. 7(a), d = 75 nm and r varies from 15 nm to 35 nm; whereas in Fig. 7(b), r = 25 nm and d varies from 55 nm to 100 nm. It is observed in Fig. 7(a) that increasing the radius of each nanosphere to 35 nm (blue lines), both the 'Proper 1' and 'Improper 1' modes present lower losses for a wider frequency range than the array in Sec. 3 (green lines). Also, notice that by increasing the nanospheres' radius, the mode 'Proper 1' is "shifted" towards lower frequencies, whereas the mode 'Improper 1' is "shifted" towards higher frequencies in Fig. 7(a). It is noteworthy that for r = 35 nm the backward proper mode has a very small imaginary part in correspondence of / 0.38 kd π ≈ (where this mode is bound); if excited at this particular frequency, this mode would travel long distances from the excitation source. Moreover, z β is quite large, thus causing short plasmonic wavelengths. In Fig. 7(b), we show the mode dependence on the array period: it can be noticed that the imaginary part z α increases with increasing periodicity. Notice that this is in agreement with what it has been previously shown in [20].
In the L-pol case, the dependence with respect to the array geometrical parameters of the forward 'Improper 1' (physical leaky only for 0 z k β < < , nonphysical bound for z k β > ) is illustrated in Fig. 8. The performed parameterization is the same as in Fig. 7. As for the T-pol case previously analyzed, it can be observed in Fig. 8(a) that increasing the radius of each nanosphere to 35 nm (blue lines), the 'Improper 1' mode presents lower losses for a wider frequency range than the array in Sec. 3 (green lines). Furthermore, notice that for large nanospheres' radius the frequency range for which the mode 'Improper 1' is physical is wider (this range is about 0.01 for the red curve and about 0.21 for the blue curve, using the scale / kd π in Fig. 8(a)) and "shifted" towards lower frequencies. In Fig. 8(b), we modify the periodicity of the array with respect to the case in Fig. 5, and it is noticed that the imaginary part z α decreases for decreasing d. Again, this result is in agreement with what it has been previously shown in [20]. Also, notice that, similar to the case in Fig. 8(a), for small periodicity d, the frequency range for which the mode 'Improper 2' is physical is wider than the other cases with larger periodicity.

Conclusion
This work provides an extensive and comprehensive description of modal analysis in a representative case of a 1D-periodic array of metallic nanospheres in terms of bound/leaky, proper/improper, physical/nonphysical wave conditions. The subsets of physical and nonphysical waves have been identified and the general principle of their identification has been provided. A physical wave is one that can be excited by a realistic localized source, and it is classified in terms of path deformation in the complex wavenumber domain. Though the theory is involved, a simple Table summarizing the classification in terms of longitudinal and transverse wavenumbers, based on spectral theory, has been given. This is the first time such a complete physical characterization in a chain of plasmonic nanospheres has been presented, and to the authors' knowledge, this analysis for complex modes is the most comprehensive carried out so far. The formulation and classification are general. The high subwavelength regime and resonance of the array nanospheres are responsible for the interesting modal dispersion described in this paper; for example, it is very interesting to note that the T-pol physical backward leaky wave "Proper 1" in Fig. 4, increasing frequency, transitions into a physical forward leaky wave 'Improper 1'. Knowledge of complex mode propagation is useful for the development of innovative sensors, field confinement and field enhancement structures, and highly directive antennas. Mode diagrams can be engineered by modifying the array parameters, thus allowing the design of waveguiding structures by using physical bound modes (such as 'Proper 2' in T-pol, and 'Proper 1' in L-pol) and radiating structure by using physical leaky modes (such as 'Improper 1' in L-pol). The physical validity of these last improper solutions related to non-principal Riemann sheets of the Hankel function involved in cylindrical problems need further analysis, therefore they are not shown in this paper and will be the focus of future work. Nevertheless, here we present a detailed physical description of the modes previously reported (for example in [21,26]), distinguishing between proper and improper ones, forward and backward, bound and leaky, in a well defined framework. We would like to point out that the study of longitudinal complex modal propagation in open or non-uniform cylindrical domains, like a dielectric rod or a waveguide with an inner dielectric rod (simpler cases than the present one because they are not periodic), is an old subject [50-53]. However, not much has been done in the study of modal solutions in cylindrical problems involving other Riemann sheets of the Hankel function. Several papers have been discussing the physical validity (individual excitability) of complex modes, even reaching contradicting conclusions [50-53], as discussed in [54]. Recently, the authors of [54] claimed they sorted out old misconceptions regarding the physical validity of certain complex modes. Here we have provided a complete description regarding the principal Riemann sheet of the Hankel function but as already said further studies are required to analyze the solutions obtained from non-principal Riemann sheets. These are necessary to establish the wavenumber symmetry in the z k complex plane, and indeed, according to [55], the complex conjugate symmetry of improper solutions can be obtained only considering solutions in adjacent Riemann sheets of the Hankel function (i.e., with different m values). What described in [55] has been confirmed here for the mode 'Improper 1' in T-pol in Figs. 3 and 4 (where we also found the wavenumber close to the complex conjugate solution correspondent to the mode 'Improper 2' by using m = +1, not shown), and for the mode 'Improper 1' in L-pol in the I and III quadrants of the complex z k plane in Figs. 5 and 6 (where we also found the wavenumber close to the complex conjugate solution correspondent to the mode 'Improper 1' in the II and IV quadrants of the complex z k plane by using m = +1, not shown). In this paper, in accordance with [55], symmetry of the proper solutions has been found in the proper plane (i.e., with m = 0).