Effective model and investigation of the near-ﬁeld enhancement and subwavelength imaging properties of multilayer arrays of plasmonic nanospheres

When a small radiating or scattering object is placed near a multilayer array of plasmonic nanospheres, on the other side the optical near ﬁeld is enhanced due to the excitation of resonant modes in the layers. For some particular frequencies, the ﬁeld behind the array is concentrated in a subwavelength region, creating a super resolution effect. Resonating layers are able to reproduce (transport) part of the evanescent spectrum to the other side of these layers which otherwise would decay rapidly. We explore the mechanism of evanescent ﬁeld transport and subwavelength ﬁeld concentration on the other side of the layered material and show the relationship between near-ﬁeldenhancement,ﬁeldconcentration,andmodaldispersioncharacteristics.Adetailedinvestigationofthese phenomena is carried out by using an effective numerical model based on the array scanning method (ASM) combined with the Ewald method to accelerate the convergence of the dyadic Green function calculation. The subwavelength-sized spheres forming the arrays are represented as single dipole radiators, and the model of their interactions takes into account all the radiative and reactive ﬁeld components.


I. INTRODUCTION
Recently a new approach to the realization of super resolution imaging devices was suggested [1][2][3][4]. It was based on the possibility to enhance evanescent fields and form subwavelength images, more properly, subwavelength field concentrations, with the use of arrays of small resonant particles such as plasmonic nanospheres. This "superlensing" approach has some advantages over other approaches which traditionally use a bulk double-negative metamaterial slab [5] or single-negative one (in devices for only one polarization) as proposed in previous studies by using one [5,6] or multiple [7] thin layers of silver. In 2004, subwavelength imaging was experimentally demonstrated in double arrays of resonant electric dipoles [1] and, in 2005, in a "superlens" formed by two parallel planar arrays of coplanar split ring resonators (SRRs) [8]. In the papers in Refs. [9][10][11] one can find further theoretical and experimental studies of such subwavelength imaging devices based on the excitation of magneto-inductive waves at radio frequencies. There, the authors theoretically and experimentally showed that they could achieve subwavelength imaging using two layers of split-ring resonators at frequencies close to the boundary of a stopband and a passband of modes traveling along the layers.
Different approaches to form subwavelength field concentration using thin layers were explored in Refs. [12][13][14][15][16]. These papers referred to the subwavelength focusing and not to the subwavelength imaging in the sense that the grid was able to generate a field concentrated in a specified subwavelength region. This subwavelength focusing is obtained by transforming an incident plane wave into a complex spatial spectrum containing evanescent components such that the phase contributions from the planar grid elements sum constructively only in a subwavelength region. In Ref. [14] the subwavelength imaging was mentioned, however, its physical mechanism is different from that considered in the present paper, as it did not use resonant response in the grids and does not allow resonant amplification of evanescent field components. Using the approaches of Refs. [12][13][14][15][16], an imaging system can be realized only using mechanical scanning.
The principle of field concentration used in the device investigated in this paper, as well as in Refs. [1][2][3][4], is based on resonant excitation of surface modes that are characterized by large tangential wave numbers. As it was shown in Ref. [17], chains of silver nanospheres may exhibit dispersion curves with rather flat regions, so that oscillations with a wide range of wave numbers can be resonantly excited by external evanescent fields. This feature of the arrays of small resonant particles can be used to design devices for subwavelength imaging or subwavelength field concentration in the optical region. In Ref. [4], such a possibility was demonstrated with the use of two parallel two-dimensional periodic arrays of metal nanospheres. The superlensing properties of the structure were demonstrated by calculating the field distribution in the image plane on the other side of the array. Note that the array had a finite extent in the cited paper and only two layers were considered [4].
In this paper we consider single-layer, double-layer, and multi-layer periodic arrays of electrically small metal spheres, with an infinite extent along two directions, excited by a single electric dipole source or two electric dipole sources on one side. This layered "artificial material" is able to provide field enhancement and subwavelength focusing effects. Both phenomena are investigated in detail and related to the modal dispersion characteristics. There are several advantages of investigating an infinite structure instead of the one with a finite extent, as in Ref. [4]. Some basic principles of such devices were not understood, and we think that it is important to simplify the problem and analyze the phenomena without the contribution of the reflected waves at the array truncations. Furthermore, as we will show in this paper, the analysis of an infinite structure (excited by a single source) can be also numerically advantageous compared to that of a finite structure. The approach suggested in this paper uses the array scanning method (ASM) [18][19][20][21][22][23][24] that expresses the field produced by a single source in a periodic environment via a spectral integral over the Brillouin zone.
In Refs. [25,26] a full-wave approach to calculate plane waves scattering by layers of two-dimensional (2D) periodic arrays of metal spheres possibly alternating by dielectric slabs has been developed. The field scattered by each sphere is represented as a sum of spherical waves. In this paper, we use the single dipole approximation (SDA), which takes into account only the field of the electric dipole spherical wave for each sphere. This approach is reasonable in our case as the spheres are much smaller than the wavelength and the electric polarizability is predominant near the plasmonic resonance. Accordingly, each sphere is modeled as a radiating dipole with a given polarizability. Since the problem is periodic in two directions, the calculation of the required periodic dyadic Green function (GF) is accelerated with the Ewald method [27][28][29], whose expressions for dyads are shown here for the first time. This approach is very efficient for the case of dipole excitation and for that of mode analysis. Besides providing an effective model suitable to analyze this kind of problem, in this paper we investigate how the field enhancement is related to the resonant modes excited in the layered array, and how these two are related to subwavelength focusing. In particular, we analyze how the field is transported by these resonant modes across the layers and how this creates field enhancement, which can be seen equivalently as an amplification of the evanescent field components. Since evanescent fields are transported away from the source, for certain frequencies they can add coherently to form a subwavelength region of strong field concentration, creating the subwavelength focusing effect.

II. MODEL OF DIPOLE EXCITATION FOR INFINITE ARRAY OF SPHERES
The structure under investigation is shown in Fig. 1. It consists of an infinite array (of one or more layers) of metallic nanospheres, which is periodic along the x and y directions. The radius of spheres r P is considerably smaller than the wavelength λ in the host material (free space in this paper). The periods of the array are a and b with respect to x and y, respectively. The transverse positions of the spheres are ρ mn = max + nbŷ, m,n = 0, ± 1, ± 2,..., and thus the m = n = 0 sphere is located at (x,y) = (0,0). Here bold symbols denote vectors and unit vectors are denoted by a hat ( ) . In this paper we focus on multilayer configurations as in Fig. 1(b), with layers at longitudinal positions z 1 , z 2 , . . . ,z N , where N is the number of layers. The array is excited by a dipole source with the dipole moment p S placed at r S = x Sx + y Sŷ + z Sẑ .
The field excited by a single electric dipole in the proximity of an infinite periodic structure can be efficiently calculated using the array scanning method (ASM) [18][19][20][21][22][23][24]. Accordingly, the solution of the problem in Fig. 1 is retrieved from the solutions of the associated auxiliary problem in Fig. 2, where the periodic array of spheres is excited by a periodic system of linearly phased dipoles placed at positions r S,mn = r S + ρ mn , m,n = 0, ± 1, ± 2,.... Here k t = k xx + k yŷ is the transverse excitation wave vector.

016607-2
According to the ASM method, the field E(r, r S ) at an arbitrary point r produced by a single dipole source at r S is given by [23,24] E(r, r S ) = ab (2π ) 2 (2) where E ∞ (r,r S ,k t ) is the field at r, in the presence of the periodic array of spheres, produced by the periodic array of dipoles at r S,mn = r S + ρ mn , with the transverse wave vector k t , integrated over the Brillouin zone. The next section is devoted to the determination of E ∞ (r,r S ,k t ) and to its effective evaluation.

III. PERIODIC EXCITATION OF THE INFINITE ARRAY OF NANOSPHERES
As the spheres are much smaller than the wavelength, we can treat them as point dipoles. This approximation has been already used in several other papers (e.g., in Refs. [1][2][3][4]17,[30][31][32][33][34][35][36][37]). The electric and magnetic fields at position r, produced by an elementary electric source at position r S with the dipole moment p S , are expressed in terms of the GF as Here is the symmetric electric-field dyadic GF that provides the electric field at position r produced by a dipole at r S , is the scalar GF [the time harmonic convention exp(jωt) is assumed and suppressed], I is the identity dyadic, and the symbol ∂ 2 G/∂r 2 is defined as Furthermore, ε h is the relative permittivity of the material hosting the spheres and k = ω √ ε h /c is the wave number in the host medium.

A. One layer of spheres
Suppose that a layer of spheres is located at z = 0. The dipole moment p ∞ of the sphere at the origin, (m,n) = (0,0), is given by the expression where E ∞ loc (0,r S ,k t ) is the local electric field at the origin produced by the array of sources p S,mn plus that of all the metal spherical nanoscatterers except for the one at (m,n) = (0,0), and α is the polarizability of a metal sphere. The inverse of the polarizability α is given by [1][2][3][4]17,36,37] where r p is the radius of the sphere and is the relative permittivity of metal represented here using the Drude model with the plasma radian frequency ω p and the damping radian frequency ω D , the values of which depend upon the adopted metal. The dimensionless parameter ε ∞ in Eq. (10) is chosen to better fit the real part of the material permittivity in the frequency range of interest. The local field E ∞ loc (0,r S ,k t ) is thus given by the sum is the electric-field dyadic GF for the periodically phased array of dipoles and G ∞ (r,r ,k t ) is the corresponding 2D periodic scalar GF The term G ∞ (r,r ,k t ) in Eq. (11) is the same dyadic GF as in Eq. (12) without the (m,n) = (0,0) term, and thus it is not singular at r = r . Therefore, the Mexican hat (˘) denotes here the regularized GF defined as G ∞ (r,r ,k t ) ≡ G ∞ (r,r ,k t ) − G(r,r ). After the substitution of Eq. (11) into Eq. (8) one has [4] (14) which leads to the linear system Note that the symmetric dyadic A(k t ) requires the computation of three diagonal terms and two equal elements out of its diagonal (A xy = A yx ), when represented in the Cartesian system, with a total of four independent values.
The periodic GF in Eq. (12) can be evaluated very efficiently by using the Ewald method, as shown in the Appendix. Besides its fast convergence properties, the Ewald method also has the great advantage that it converges very rapidly for complex transverse wave vectors k t , whereas the series in Eq. (12) would diverge in this case. This feature is useful when looking for complex modes that propagate along the array, which will be the subject of future studies, or when, as done in this paper, a deformation of the integration path is applied in Eq. (2).
Once the solution p ∞ of Eq. (15) for the dipole moment of each nanosphere is determined, the electric and magnetic fields at an arbitrary position r, produced by the array of exciting dipoles at r S,mn , phased with the transverse wave vector k t , can be evaluated by The field produced by a single dipole at r S is obtained by inserting E ∞ (r,r S ,k t ) or H ∞ (r,r S ,k t ) into Eq. (2). It is remarkable that for every wave vector k t in Eq.(2), only the 3 × 3 linear system in Eq. (15) needs to be solved to evaluate E ∞ (r,r S ,k t ). Therefore, even though the integral evaluation of Eq. (2) requires sampling the Brillouin zone in several points (for example, in 20 × 20 spectral points), the overall computation time is not long [23].

B. A note on the computational cost
For example, suppose that S operations are needed to solve the 3 × 3 linear system (15), and the Brillouin zone is discretized in P × Q spectral points (P along k x and Q along k y ). The overall cost to determine E(r,r S ) in Eq. (2) is thus P × Q × (22 + S), where 22 is the total number of dyadic GF calculations [four elements in the matrix A(k t ), six elements in the symmetric matrix in the right-hand side of Eq. (15), and 12 elements in the symmetric matrices in Eq. (17)]. For simplicity we can assume that S = 3 3 operations are needed to solve a system with three unknowns, as, for example, needed when using the lower-upper (LU triangular matrices) inversion method. The error in E(r,r S ) in Eq. (2) resulting from the periodic sampling of the Brillouin zone is represented by the additional field contributions coming from the "virtual source images" at locations r S + mP ax + nQbŷ, with m,n = 0, ± 1, ± 2, . . . [23].
When small losses are present, in the spherical particles or in the ambient, the error can be negligible already with P and Q in the range of 20. We compare the computational cost of the ASM with P = Q = 20 spectral samples, with a case of a truncated array of 20 × 20 spheres with the total of 20 2 × 3 = 1200 unknowns. The ASM computational effort is P × Q × (22 + S) = 19,600 operations, compared to the truncated case which needs the computation of 20 2 × 6 matrix elements to fill the linear system and at least (20 2 × 3) 3 = 1.728 × 10 9 operations to solve the linear system (we recall that for simplicity we estimate S as the cube of the number of unknowns 20 2 × 3). Though in general the ASM technique is very advantageous, the case studied in this paper is one of the toughest because the superresolution achieved with the array device is produced by exciting a special mode (see Sec. V). In mathematical terms this is described by complex poles in the ASM integrand in Eq. (2), which are close to the real k x , k y axes in all quadrants. This is not a standard situation for periodic structures.
The presence of losses helps the computation in this case since the poles are slightly off the real axes. For the lossless case, poles may be on the real axes and it is recommended to perform a path deformation without enclosing poles when calculating Eq. (2) and to use adaptive integration. In all the problems treated, which follow, numerical convergence was achieved by using path deformation and adaptive integration.

C. N layers of spheres
The previous formulation is here generalized for the case of N-layer arrays [ Fig. 1(b)]. We obtain the following matrix equation for the determination of the dipole moments p ∞ n : Here p ∞ n is the dipole moment of the spheres at the (0,0) position of the layer n (n = 1,. . ., N). Once the solution p ∞ n , n = 1,. . ., N of Eq. (19) has been determined, the electric and magnetic fields at an arbitrary position r, produced by the array of exciting dipoles, phased with the transverse wave vector k t , are then evaluated by

D. Excitation by a plane wave
In case one wants to excite the layers with a plane wave to analyze transmission and reflection properties, for example, the only modifications needed are the following: one needs to replace the right-hand side of Eq. (15) with E inc (r,k t ), evaluated at r = 0 in the case of a single layer, whereas in the case of N layers one needs to replace the right-hand side of Eq. (19) with E inc (r,k t ) evaluated at r = z mẑ .

IV. MODAL PROPAGATION IN THE ARRAY OF SPHERES
Modes are found by looking for the eigensolutions of Eq. (15), which are those that satisfy Eq. (15) without the forcing term, that is, The search is made numerically. Fast evaluation of the periodic dyadic GF using the Ewald method shown in the Appendix renders the search for the eigenmodes very efficient. This Ewald GF evaluation is faster than the accelerated scheme suggested in the paper of Simovski et al. [38], due to its Gaussian convergence rate.
Dispersion curves for single-and double-layer arrays of silver nanospheres placed in a host medium with ε r = 1 are shown in Figs. 3 and 4, respectively. The radius of the nanospheres is r p = 25 nm, the periods are a = b = 73 nm, and the distance between the layers in the case of the double-layer array is h = 95 nm. The Drude model for silver  [39,40]. In this section devoted to modal analysis we assume ω D = 0 since we focus only on the determination of the dispersion diagram with real modes. Complex modes in such arrays will be the subject of future investigations.
The dispersion curves for the single-layer array are shown in Fig. 3, for 0 k x π/a and k y = 0, where we show the three supported modes. They correspond to the dipole moments p ∞ parallel to the x, y, and z axes. The mode polarized along y and traveling along x cannot be excited by the electric dipole source along z, considered in the next section. The small box in the figure represents the region where the mode polarized along z has a flat dispersion curve. This peculiarity is important to create the sought subwavelength focusing and possibly the superresolution, as was already discussed in Ref. [3].
The same kind of plot for the double-layer array is shown in Fig. 4. There are six possible modes traveling along the x direction. Two of them have dipole modes parallel to the y direction and as already stated above, they cannot be excited by a z-directed dipole. The four other modes have dipole moments p ∞ 1 and p ∞ 2 , in layers n = 1 and n = 2, with both x and z components.

V. SUPERRESOLUTION AND FIELD ENHANCEMENT
We show here that two main effects are provided by a single-layer or a multilayer array of nanospheres excited by a nearby dipole source: subwavelength focusing and near-field enhancement. Now, the source is situated at the coordinate origin, its dipole moment is directed along the z axis, and its magnitude is chosen as p S,z = 10 −30 C/m. The periods of the array are a = b = 73 nm. The first layer is placed at z 1 = 95 nm. In the case of a multilayer array of nanospheres, the interlayer distance is h = 95 nm and the relative permittivity of the media hosting the nanospheres is ε h = 1. The field is calculated in the "image plane" z = (N + 1)95 nm, at the distance 95 nm above the last layer (the number N of layers may vary, see Fig. 2). The radius of nanospheres is r p = 25 nm. In Fig.  5 the near-field enhancement |E z (r obs )| 2 / |E inc z (r obs )| 2 , defined as the ratio between the field intensity with and without the array at the position r obs = (0,0,(N + 1)95 nm), is shown versus frequency for different numbers of layers N = 1,2,3,4. The field without the array corresponds to the field E inc z of an isolated dipole. The near-field enhancement is observed in a wide frequency range for any number of layers, and increases when increasing the number of layers.
The frequencies relative to the strongest near-field enhancement for one and two layers in  spatial spectrum, which is reconstructed at the observation point at a specified frequency due to the resonant array. The presence of the array results in the excitation of a strong field in the vicinity of the layers since several spectral components (waves with different wave numbers) are excited by the source and are able to transmit their evanescent field to the observation point. The field enhancement is particularly strong for a large number of layers. Therefore the device made by layers of nanospheres can be used to amplify certain evanescent field components, with respect to what happens in a homogeneous medium. Figure 6 shows the field distribution moving along the x direction, for y = 0, in the "image" plane z obs = (N + 1)95 nm for a double-layer array (N = 2), at the two  frequencies (f = 787 THz and f = 808 THz) yielding the maximum of the near-field enhancement shown in Fig. 5. The half-power width of the "image" is marked by arrows and the result in Fig. 6 clearly shows that subwavelength field localization (denoted as "focusing" by some researchers) is achieved at the frequencies which correspond to the maximum near-field enhancement. For the frequency f = 787 THz, the half-power width of the "image" is significantly smaller than a quarter wavelength. This frequency corresponds to the flattest part of a dispersion curve in Fig. 4, which means that the largest part of evanescent spectral components excited at that frequency are enough to reconstruct a subwavelength field localization. For a reference we have also reported in Fig. 6 the field that would result without the array, which corresponds to the incident field (dashed line). As it is much weaker than the field with the array, this is another confirmation of the near-field enhancement already seen in Fig. 5. Figure 7 shows similar results for three-and four-layer arrays, showing for the first time that subwavelength focusing can be achieved for more than two layers, and that subwavelength focusing is obtained at more frequencies. Field localization with a half-power width, approximately equal to 0.15λ, is achieved, for example, at f = 790 THz for four layers.
To better illustrate the physical mechanism, in Fig. 8 we plot the field across the array layers for the frequency f = 808 THz, which corresponds to the maximum field enhancement for two layers, and is very close to that for three and four layers of nanospheres. Figures 8(a) and 8(b) represent the x and z components, respectively, of the E field along the line parallel to the z axis with x = a/2 and y = 0. In all considered cases the field component E y vanishes for symmetry reasons. It is interesting to note that the field has a similar behavior for any number of layers. It decreases more slowly than the incident field, for increasing z, because of the presence of the layers of nanospheres. The plots clearly show that the presence of the layers permits to maintain the field strength at larger distances than without these layers. The incident field decreases as 1/r (r is the distance from the dipole) and rapidly becomes very weak. Instead, when the nanosphere layers are present, a small incident field (to any layer) is able to excite the resonant mode, which is thus responsible for producing a strong scattered field near the layer. This field decays away from each layer, but is able to strongly excite the subsequent resonant layer.
Only above the last layer of nanospheres does the field decrease more rapidly than the incident field. It decays exponentially since it is related to a mode in the array propagating in the transverse direction, and the "spatial" component excited by the source is already too weak (we stress here that a source near a periodic structure excites the modes and a "spatial" field that collects the continuous spectrum of wave numbers, as shown in Sec. V of Ref. [21], and in Ref. [41]). The "maintenance" of the evanescent field across the layers can be interpreted as "enhancement" or "amplification" when compared to the field without the layers (the incident field), which decreases rapidly. Since the field maintained (or transferred) across the layers is made of the evanescent part of the excited spectrum, it has large transverse wave-number components (small spatial wavelengths) and is thus able to reconstruct the subwavelength field localization far from the source.
The results of the investigation of the superresolution possibility are shown in Fig. 9. The double-layer array is excited by two z-directed dipole sources with dipole moments The frequency is f = 789 THz (λ = 380 nm). Thus, the distance to the observation plane from the dipole sources is 0.6λ, and the distance between the sources is 0.29λ = 110 nm. In Fig. 9 the horizontal dash-dotted line corresponds to the half-power level of the field in the presence of the array. Note that considering the half-power level in Fig. 9(a) with the two layers of nanospheres, the "images" of the two sources at y = ±55 nm are well separated in the observer plane at z = 229 nm. Instead, without the two-layer array of nanospheres the two field contributions, each produced by a source, are not distinguishable one from the other, in the observer plane at z = 229 nm [ Fig. 9(a), dashed line].
In Fig. 9(b) we also show the field distribution along the x direction, with y = 0 and z = 229 nm (note that the two sources are displaced along y). There are strong interference maxima which are due to the superposition of the fields of two coherent contributions. Indeed, as already seen in Figs. 6 and 7, the field produced by a single source in the observation plane has several maxima. From a field analysis near each source, not shown here (an example could be seen in Ref. [42]), the field is strong in a ring region of the image plane (in the x,y plane, for z fixed) around each of the sources. Therefore, interference maxima produced by the two sources may arise at the points of intersection of the corresponding two rings. In Fig. 9(b) these maxima almost reach the half-power level of the field maximum in Fig. 9(a). Their level can be reduced by optimization. The three-dimensional (3D) Fig. 9(c) illustrates the distribution of the electric field maxima in the image plane. The results clearly show strong enhancement of evanescent fields as well as subwavelength focusing, where the two main peaks corresponding to two nearly placed sources are clearly resolved. On the other hand, the presence of parasitic interference maxima at the level of some 70% of the main peak shows limitations for the imaging applications in our current design.

VI. INFLUENCE OF THE DISTANCE h BETWEEN TWO LAYERS
Figures 10 through 12 illustrate the influence of the distance h between layers, when only two layers are considered, on the nanospheres' excitations and on the resulting near-field enhancement. All the other parameters were fixed: The dipole source is situated at r S = 0, the layers have periods a = b = 73 nm, and the bottom one is 73 nm above the source. The observation point is placed at a distance of 73 nm above the upper layer, at r obs = (2a+h)ẑ. The distance h between the layers varies. In Fig. 10 it is interesting to note that the nanosphere dipole moment p 1 , representing the dipole strength of the nanosphere at (x,y) = (0,0) in the layer closest to the source, is larger than the exciting dipole moment p S . In other words |p 1 |/|p S | > 1 in a wide range of frequencies and distances h. Note that when h is large, the excitation of p 1 in the first layer of nanopsheres does not depend on h anymore, as shown by the parallel vertical lines in Fig. 10, since the interference between the two layers of nanospheres becomes weak. The strongest near-field enhancement corresponds to the strongest excitations of both layers, and especially of the second layer p 2 . However, the maxima in Figs. 10 and 11, showing |p 1 |/|p S | and |p 2 |/|p S |, exhibit similar trends (darkest regions), which correspond to the strongest near-field FIG. 11. As in Fig. 10, with the difference that the strength of the dipole moment p 2 of the second layer is shown. enhancement, as shown in Fig. 12, in the presence of the two layers of nanospheres. These results are important to show that at any specific frequency there is a corresponding optimal distance h that maximizes the field enhancement. For each small value of h there are two frequencies giving maximum near-field enhancement. Figure 11 shows that to have a strong excitation of the second layer the distance h should not exceed twice the period of the array because otherwise the second layer could not be excited strongly enough by the evanescent field from the first layer. In general, strong excitation of the second layer implies strong near-field enhancement, as shown in Fig. 12. However, as also shown in Fig. 12, larger interlayer distances h may also yield to strong near-field enhancement because, even though in this case the second layer would not be strongly excited, the incident field (the direct field, in absence of the layers) would be very weak because of to the increased observation-source distance |r obs − r S | = 2a + h. As an example, one should note the 30 dB near-field enhancement at λ = 375 nm for h = 190 nm (i.e., for |r obs − r S | = 336 nm).

VII. CONCLUSION
The near-field enhancement and subwavelength focusing properties of layers of infinite arrays of metallic nanospheres have been investigated at optical frequencies for a dipole source oriented orthogonally to the array plane. Earlier, similar investigations were performed for arrays of a finite extent and only for one and two layers. The main advantages of the present formulation are that certain physical phenomena are not hindered by truncation effects, and the possibility to use more efficient algorithms to numerically facilitate the investigation of multilayer arrays, based on the array scanning method (ASM) combined with the Ewald method, whose dyadic expressions have been provided here for the first time.
The near-field enhancement and subwavelength focusing properties are investigated for a single layer and multilayer arrays of metal nanospheres. It is clearly shown that the frequencies providing these phenomena correspond to the flat parts of the dispersion curves for the modes polarized orthogonally to the array. This property provides a way to predict and optimize the near-field enhancement and subwavelength focusing of an array based on its dispersion characteristics.
The influence of the distance between the layers of the array on the resonant excitation of each layer and on the frequency of the maximal near-field enhancement has been investigated and certain optimal regions have been determined. Tunability of the operational frequency, imposed by the plasmon resonance of the chosen metal, can be achieved by using dielectric coremetal shells.
In summary, the explored device made by layers of metal nanospheres can be used for two purposes: near-field enhancement and subwavelength focusing. This paper analyzes and stresses certain properties which show how to control such phenomena. However, we should stress that this is not an imaging device in the usual sense of this term because by measuring the field distribution in the "image plane" it is not possible to uniquely reconstruct the distribution of possibly volumetrically distributed sources. The presence of interference maxima at some lower level also make the superresolution application still challenging.

APPENDIX: ACCELERATION OF THE PERIODIC GREEN'S FUNCTION BY USING THE EWALD METHOD
We show here how G ∞ (r,r ,k t ),G ∞ (r,r ,k t ), ∇G ∞ (r,r ,k t ), ∂ 2 G ∞ ∂r 2 (r,r ,k t ), and ∂ 2G∞ ∂r 2 (r,r ,k t ) discussed in the text are calculated with the Ewald method [27][28][29] at point r = ρ + zẑ, due to a dipole source at r = ρ + z ẑ. This approach is significantly faster than that based on the analytical extraction of the slowly convergent part of the representing series [36].