Wave dynamics by a plane wave on a half-space metamaterial made of plasmonic nanospheres: a discrete Wiener-Hopf formulation

A rigorous analytical solution for the description of wave dynamics originated at the interface between a homo- geneous half-space and a half-space metamaterial made by arrayed plasmonic nanospheres is presented. The solution is cast in terms of an exact analytical representation obtained via a discretized Wiener – Hopf (WH) technique assuming that each metallic nanosphere is described by the single dipole approximation. The solution analytically provides and describes the wave species in the metamaterial half-space, their modal wavenumbers and launching coefficients at the interface. It explicitly satisfies the generalized Ewald – Oseen extinction principle for periodic structures, and it also provides a simple analytical solution for the reflection coefficient from the half-space. The paper presents a new WH formulation for this class of problems for the first time, and describes the analytical solution, which is also tested against a purely numerical technique. While only the case of orthogonal plane wave incidence and isotropic inclusions is considered here, the method can be easily generalized to the oblique incidence and anisotropic constituent (tensorial polarizability) cases. © 2011


INTRODUCTION
It is known that, when a plane wave impinges on a surface of a bulk metamaterial, it generates waves that decay away from the interface and, for what concern transmission and reflection, the metamaterial may not be considered a homogeneous material bounded by the interface itself [1][2][3][4]. In this paper, we rigorously analyze the problem of a plane wave orthogonally incident onto a metamaterial of semi-infinite extent made by a regular lattice of plasmonic nanospheres as shown in Fig. 1. With a rigorous formulation, we are able to determine the wave dynamics arising from the interface, creating reflection and transmission. However, these two simple concepts, trivial in standard homogenous materials, assume a more complicated aspect in the present problem because of the periodicity. Several waves are excited at the interfaces that travel away from it, and in most of the cases they decay exponentially. The rigorous analytical treatment allows for a precise physical description of these phenomena.
Our formulation consists in representing each sphere as a dipolar scatterer according to the single dipole approximation (SDA) [5,6], with strength proportional to the nanosphere polarizability, which exhibits a resonant behavior because of the fundamental plasmonic mode. Besides this simplified representation for the nanosphere scattered field, which is, however, well justified for small nanospheres near the plasmonic resonance, the rest of our formulation rigorously accounts for all the coupling interactions. This leads to a linear system with infinite dimensions, solved here using the discretized version of the Wiener-Hopf (WH) method [7][8][9][10][11], which specifically uses the concept of the Z transform, which is ideal for periodic problems, where the periodicity can be regarded as space sampling. The same method has been used in [11] for analyzing currents in planar semi-infinite arrays of conducting strips and their radiated fields. The reader is addressed to the Introduction of our previous paper [11] for more complete information about the discretized WH method for scattering problems in periodic structures.
In this paper, the discrete WH technique is used for the first time to study a periodic arrangement of scatterers occupying a semi-infinite space, which is regarded as a semi-infinite collection of periodically stacked planar arrays of nanospheres. The same problem was originally studied in Mahan and Obermair's seminal work [12], and in the papers thereof originated [13][14][15][16][17][18], within the so-called "nearest neighbor approximation," which, however, loses its validity when the crystal period is not small in terms of the wavelength. In [1], a comprehensive and critical review of this stream of literature is given. In [1], the problem solution validity is extended on the basis of the a generalized version of the Ewald-Oseen extinction principle, valid for crystals or periodic lattices. Namely, the Ewald-Oseen extinction principle [19][20][21], for which the polarization in a dielectric is distributed so that it cancels out the incident wave and produces the propagating wave, is extended considering Floquet harmonics in a periodic structure. Therefore, in [1], the interaction between particles is described in a complete way, and not only through the dominant Floquet wave as in the "nearest neighbor approximation." Consequently, the accuracy of [1] is extended beyond the long-wavelength limit. In [1], besides a general formulation, a detailed solution is provided analytically only for the special case of crystals formed by small scatterers, which can be treated as point dipoles with fixed orientation, and the split-ring resonators are then considered as an example. However, in [1], an analytic solution cannot be given for a crystal of metallic spheres, which can also be replaced by point dipoles but whose orientation depends on the direction of the external field. The limitation of the available analytic solutions in [1] resides in the use of the characteristic function method, according to [18]. Indeed, authors in [1] say "We hope that with some modification this method can be used for other special cases as well." The solution provided in the present paper by using the discretized WH method can be regarded as an extension of the characteristic function method to which it is strictly related. The main improvement in using the WH method is that it provides a constructive algorithm to build the solution exploiting the WH factorization and no assumption has to be made a priori. For simplicity, only the case of orthogonal incidence and isotropic polarizability constituents is here treated, but the methodology can be easily generalized to the case with oblique incidence and anisotropic scatterers. Furthermore, most of the wave dynamics, here observed with this rigorous method for the first time, would hold also for the oblique incidence case.
Assuming that the interface between the homogeneous half-space and the metamaterial is the x-y plane ( Fig. 1), it is found that a plane wave incidence would generate a number of modes propagating along the z direction, each one represented as a summation of spatial Floquet harmonics. Most of the modes are evanescent in the z direction. It is important to stress that these modes are those that are admitted in a metamaterial with infinite extent. However, the technique here proposed is able to determine exactly their launching coefficients (their weights), and these modes can be interpreted as launched by the interface. We also show that the WH method permits to determine the exact value of the plane wave reflection coefficient.
These findings can be applied within homogenization theories since they show (and quantify) exactly that besides a transmitted Bloch wave, the interface often generates other waves with significant amplitude that decay exponentially away from the interface. Indeed, these phenomena associated to the interface have posed some problems in the characterization of metamaterials in terms of effective bulk permittivity and permeability. To overcome such homogenization problem, some authors have proposed the insertion of a Drudetype transition layer (or transition sheet) [2,3] that accounts for the transition between free space and the bulk homogenized metamaterial. Since our solution rigorously describe wave dynamics also in proximity of the interface, it is intended to shed light into this class of phenomena for future better understanding. Also, the formulation and interpretation analysis could be applied to investigate the wave dynamics at an interface between free space and a photonic crystal.

STATEMENT OF THE PROBLEM
We consider an artificial material (or metamaterial) made of small plasmonic nanospheres arranged in a Cartesian lattice (Fig. 1), embedded in the free space or in an infinite homogeneous background material. The metamaterial fills the z ≥ 0 half-space beyond the interface at z ¼ 0. The z ≤ 0 half-space is assumed to be free of particles. Nanospheres are located at r lmn ¼ lax þ mbŷ þ ncẑ, with a, b, and c denoting the lattice periodicities, while l; m ¼ 0; AE1; AE2; ::::: and n ¼ 0; 1; 2; ::::: are the associated indexes, along x, y, and z, respectively; the caret ∧ denotes unit vectors. Any plasmonic nanosphere, according to the SDA [5,6], is represented by a single electric dipole p lmn and by a scalar isotropic polarizability α, which relates it to the local electric field by We assume that the metamaterial is illuminated by a plane wave traveling in theẑ direction and impinging orthogonally on the z ¼ 0 interface. Its incident electric field is given by denoting the plane wave polarization vector and k denoting the wavenumber in the background material. The time convention expðjωtÞ is assumed throughout the paper and therefore suppressed.

FORMULATION
The scattered electric field at a generic point r produced by a single metal nanosphere, represented as single electric dipole at r 0 with the dipole moment p 0 ¼ αE loc ðr 0 Þ, where α is its polarizability [5,6] and E loc is the local electric field at r 0 , which is expressed as through the unbounded homogeneous material electric dyadic Green's function in which I is the identity dyad, ϵ is absolute dielectric permittivity of the background material, and is the scalar Green's function. Performing the differentiation in Eq. (3), we obtain where R ¼ r − r 0 andR ¼ R=R.
The local electric field E loc ðr 0 Þ represents the field produced by the incident wave E inc and all the other nanospheres, and its general expression will be clear in the following. Under orthogonal plane wave excitation E inc ðrÞ ¼ E inc 0 expð−jkzÞ, the periodicity of the problem along x, y will result in a spatial independence of m and n for the polarization p lmn ¼ p n ; i.e., all the particle in the same nth layer, n ¼ const, will have the same polarization (magnitude and phase). Therefore, it is convenient to represent the half-space metamaterial as a semiinfinite stack of layers. Accordingly, the field radiated at a generic point r by the nth metamaterial layer can be conveniently expressed as where we have defined the layer Green's function as By using the Poisson summation formula, the layer Green's function can be expressed, alternatively to its spatial representation (7), via its spectral Floquet wave representation: in which k z;rs ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi y;s q , with k x;r ¼ 2πr=a, k y;s ¼ 2πs=b, and Imfk z;rs g ≤ 0, and k AE rs ¼ k x;rx þ k y;sŷ AE k z;rsẑ . The choice of the sign in k AE rs accordingly to z≷0 ensures the series convergence and satisfaction of Sommerfeld radiation condition at jzj → ∞. For the calculation of the local field E loc ðr lmn Þ in the following it is also useful to introduce a regularized layer Green's function where the 00th term is excluded: which is therefore regular at r ¼ 0. Note that, because of odd symmetries, k x;r ¼ −k x;−r and k y;s ¼ −k y;−s induced by the normal incidence, G layer ðrÞ and Ğ layer ðrÞ are both diagonal as extradiagonal terms vanish in the summation (8).
The local electric field at the lmnth particle position, i.e., the incident field plus the field radiated by all the other particles (from all layers) at the location denoted by (l ¼ 0, m ¼ 0, n), is therefore given by G layer n−n 0 · p n 0 þ E inc n ; for n ¼ 0; 1; 2; :::::; ð10Þ for n ¼ 0; 1; 2; :::::; ð12Þ where, for the sake of simplicity, an x-polarized incident field is assumed. For convergence reasons, the expression of Ğ layer ð0Þ is evaluated by using the Ewald method [22]. By using Eq. (10) into Eq. (1), a linear system of equations (of infinite dimension) is set whose solution provides the particle dipole moments X ∞ n 0 ¼0 a n−n 0 · p n 0 ¼ αE inc n ; for n ¼ 0; 1; 2; :::: with a n ¼ δ 0n I − αG layer n . Since a n is diagonal, the vector problem can be separated into three independent scalar problems. The y and z problems admit a trivial vanishing solution as they cannot be excited by the x-polarized plane wave excitation, so that only the x component needs to be considered: X ∞ n 0 ¼0 a n−n 0 p n 0 ¼ αE inc 0 e −jkcn ; for n ¼ 0; 1; 2; :::::; ð14Þ withx · a n ·x ¼ a n and p n ¼ p nx . In particular, the coefficients are equal to a 0 ¼ 1 − αĞ layer xx ð0Þ and a n ¼ −αG layer xx ðncẑÞ, with n ≠ 0.
A. Z-Transformed Problem The linear system of equation to be solved [Eq. (14)] is in the form of a discrete convolution between the unknown unilateral sequence p n , with n ¼ 0; 1; 2; :::::, and the bilateral sequence a n , with n ¼ 0; AE1; AE2; :::::. Hence, it is convenient to introduce the (discrete Laplace) Z transforms [23,24] of the two sequences: To avoid confusion with the spatial coordinate z, differently to standard notation, ζ will denote the complex transform variable. The conformal mapping is used to establish a correspondence between the ζ and k z planes. The top (bottom) complex half k z plane is projected onto the region outside (inside) the unit circle of the complex ζ plane. Note that k z ðζ inc Þ ¼ k and k z ðζ ref Þ ¼ −k. Since a discrete sequence convolution becomes a multiplication between respective Z transforms, Eq. (14) is rewritten in the transformed domain as where RðζÞ is an unknown function that is regular inside the unit circle, while the pole singularity ζ inc ¼ expð−jkcÞ lies slightly inside the unit circle by assuming vanishing small losses in the ambient. Note that the inverse Z transform of the right-hand side (RHS; Appendix A), evaluated via residue at ζ − ζ inc , provides αE inc 0 expð−jkncÞ for n ¼ 0; 1; 2; ::::: (inside the metamaterial), which is the RHS of Eq. (14).

B. Factorization in the Z Spectral Plane
The Z transform of a n is determined in Appendix A and given by in which ζ rs ¼ expð−jk z;rs cÞ, with ζ 00 ≡ ζ inc . Note that, because of the even symmetry a −n ¼ a n , which results from Green's function symmetries with respect to z, it holds that Að1=ζÞ ¼ AðζÞ. Furthermore, AðζÞ exhibits an infinite number of pole singularities at ζ rs , with r; s ¼ 0; AE1; AE2; :::::, inside the unit circle and an infinite set of pole singularities outside the unit circle at 1=ζ rs . The two sets of poles accumulate at ζ ¼ 0 and ζ ¼ ∞, respectively, which are therefore function singular points, too. Poles ζ rs correspond to the z-directed wavenumber of the Floquet harmonics produced by a single layer considering propagation along the positive z direction. Conversely, poles 1=ζ rs (outside the unit circle) correspond to the same wavenumbers produced by a single layer, propagating along the negative z direction. Finally, note that electrically small periods a and b imply that ζ inc ≡ ζ 00 and ζ ref ≡ 1=ζ 00 are the only propagating wavenumbers slightly inside and slightly outside the unit circle, respectively; while all other Floquet modes (produced by a single layer) represented by ζ rs and 1=ζ rs are evanescent and the associated poles lie on the positive real axis. These can still be dominant when the distance c between layers is much smaller than the wavelength and comparable to or smaller than the periods a and b. These critical points characterize the field produced by a single layer, we now determine the critical points characterizing directly fields in the three-dimensional (3D) lattice. Waves in the 3D lattice are characterized by critical points in the ζ complex plane such that AðζÞ ¼ 0. Indeed, at those "zeros," the homogeneous version (i.e., without any excitation) of Eq. (14) admits a nontrivial (nonvanishing) solution. Note that, if ζ mode q is a zero, 1=ζ mode q is also a zero because of the symmetry property discussed after Eq. (17). Such eigenvalues ζ mode q and 1=ζ mode q correspond to Bragg's modes in the metamaterial (i.e., in the 3D lattice) traveling along the AEz direction. Again, we can split them into two sets: one set of zeros ζ mode q , q ¼ 1; 2; ::::: inside the unit circle (jζ mode q j < 1) corresponding to metamaterial modes that attenuate along the positive z direction, and another set of zeros 1=ζ mode q , q ¼ 1; 2; ::::: outside the unit circle corresponding to the same respective modes but attenuating along the opposite direction, i.e., along negative z.
By using the procedure described in Appendix B, A is factorized as in which A þ accounts only for the singularities ζ rs and the zeros ζ mode q inside the unit circle and it is therefore regular and nonvanishing on and outside the unit circle, including ζ → ∞, where A þ ð∞Þ ¼ a þ 0 , with a þ 0 denoting the first (n ¼ 0) sample of the monolateral sequence a þ n counterparts of A þ ðζÞ. Conversely, A þ ð1=ζÞ ¼ A − ðζÞ accounts only for the singularities 1=ζ rs and the zeros 1=ζ mode q outside the unit circle, and it is therefore regular and nonvanishing on and inside the unit circle, including

C. Z Transform of Polarization Solution
Using the factorization Eq. (19) in Eq. (17) leads to Since the Z transform PðζÞ of the unilateral sequence p n , vanishing for n < 0, is regular on and outside the unit circle, including ζ → ∞, where Pð∞Þ and therefore they must equal a constant, whose value αE inc 0 =A þ ð1=ζ inc Þ is calculated by evaluating the RHS at ζ ¼ ζ inc . Whence, the Z transform of the polarization sequence is obtained:

D. Solution for Polarization of Nanospheres in Metamaterial Half-Space
The polarization sequence at each nth layer, with n ¼ 0; 1; 2; :::::, is calculated via the inverse transform Eq. (A2) of Eq. (21) [23,24]: Note that, since Pðζ → ∞Þ ¼ αE inc 0 =½a þ 0 A þ ð1=ζ inc Þ, for any n < 0, the integral in Eq. (22) can be closed to infinity via Jordan's lemma (i.e., by deforming the integration path C onto a circle of infinite radius), providing a vanishing result p n ¼ 0, as expected. On the other hand, for n ≥ 0, Eq. (22) does not converge on a path at infinity but is now regular at ζ ¼ 0; hence, it is to be calculated via the Cauchy residue theorem by considering the singularities inside the unit circle C. A description of such singularities is thus in order. Note that the singularity at ζ ¼ ζ inc in Eqs. (21) and (22) is removable because A þ ðζ → ζ inc Þ ∝ ðζ − ζ inc Þ −1 → ∞, and thus by considering Eq. (18), one has Pðζ inc Þ ¼ −j2abϵE inc 0 =k. Hence, Eq. (21) exhibits only pole singularities at ζ ¼ ζ mode q , q ¼ 0; 1; 2:::::;, where A þ ðζ mode q Þ ¼ 0, which coincides with the zeros that AðζÞ has inside the unit circle. Consequently, the integral in Eq. (22) can be calculated via the Cauchy theorem as series of residues where p mode;q n ¼ C q e −jk mode z;q nc is the polarization associated to the qth mode in the semi-infinite crystal, with propagation wavenumber k mode z;q ¼ jc −1 lnðζ mode q Þ obtained via the mapping Eq. (16), and weight Note that only pole singularites ζ mode q inside the unit circle contribute to the field in the metamaterial region z > 0. No other singularities like branch points are present (cf. [11]). In light of the correspondence established, Eq. (23) is interpreted as a superposition of metamaterial modal waves transmitted in the z > 0 half-space that travel with modal wavenumbers k mode z;q and with launching coefficients C q , established at the interface. Because of the mapping Eq. (16), the only waves that can be excited at the interface and travel in the z > 0 region have wavenumber k mode z;q , with Imk mode z;q < 0. This corresponds to the physical condition that waves cannot grow when moving away from the interface. The q-modal series in Eq. (23) is usually very rapidly convergent, especially for n > 0, since one or no more than a few modes are propagating, while all the others are strongly evanescent (i.e., they attenuate along z). Note that, as stated earlier, in general, a mode with propagation expð−jkncÞ (i.e., with ζ ¼ ζ inc ) is not an admitted solution, as expected. In other words, modes inside the metamaterial propagate along z with a wavenumber different from the one of the incident wave, as explained in Section 4 in relation to the extinction principle.
In addition, it is noted that PðζÞ in Eq. (21) exhibits zeros at the poles of A þ ðζÞ; i.e., Pðζ rs Þ ¼ 0 at any rsth positive zdirected Floquet harmonic ζ rs , except for the fundamental one ζ 00 ¼ ζ inc , as stated earlier. Such a property is also a consequence of the generalized Ewald-Owseen extinction principle described in [1], as we show in Section 4.

SCATTERED FIELD
Once the dipole moment at each n layer is calculated, it is straightforward from Eq. (6) to derive the exact expression of the field scattered by the semi-infinite crystal at an arbitrary observation point r ¼ xx þ yŷ þ zẑ inside or outside the crystal half-space by using the superposition principle as where the series can be formally extended to n ¼ −∞, since p n ¼ 0 for n < 0. Again, because of the periodicity of the problem, we use the z ransforms whose standard theory was developed in [23] for sampled functions like p n . Nonetheless, by resorting to the advanced Z transform [24], both the scattered field E scat x ðrÞ and the layer Green's function G layer xx ðr − ncẑÞ, which are defined continuously as a function of any observation point r within the periodic cell, can be treated within the same framework. In this way, analogously to Eq. (14), the scattered field Eq. (25) can also be interpreted as a convolution product of sequences and calculated via the inverse Z transform x;r Þe −jðk x;r xþk y;s yþk z;t δ z Þ ζ Δ z =c k 2 z;rs − ½k z;t − ðj ln ζÞ=c 2 ; ð27Þ with k z;t ¼ 2πt=c, is the extended Z transform of the layer Green's function Eqs. (7) and (8) derived in Appendix C. In Eq. (26), the observation point z coordinate is decomposed as z ¼ nc þ Δ z , in which n ¼ ⌊z=c⌋ and Δ z ∈ ½0; cÞ denotes the so called "delay parameter" corresponding to a displacement of the observation point within the periodic cell with respect to the lattice points z ¼ nc. Note that, similar to AðζÞ, in the ζ complex plane, Eq. (27) exhibits pole singularities at ζ rs inside the unit circle and at 1=ζ rs outside the unit circle. The branch cut on the ζ-negative real semiaxis, introduced in each term of the series by the logarithm and the fractional power, is then canceled in the summation over t; therefore, it is fictitious and not present in the entire function.

A. Reflected Field
When observing the field scattered back into the homogeneous half-space, z < 0 implies n < 0 so that the integrand of Eq. (26) provides a vanishing contribution when deforming the integration path on a circle of infinite radius. In such a deformation, Eq. (26) is calculated via Cauchy theorem as the series of the residues of the integrand outside the unit circle C. Indeed, these are only the pole singularities at ζ ¼ 1=ζ rs of G layer xx ðζ; x; y; z 0 Þ, as PðζÞ is regular outside C. Hence, Eq. (26) is reduced to where W scat x;rs ¼ ðk 2 − k 2 x;r ÞPðζ −1 rs Þ=ðj2abϵk z;rs Þ is the weight of the rsth Floquet wave traveling in the negative z direction. In the case that only the dominant Floquet mode is propagating (a; b; c < λ=2), far from the interface (z → −∞) only the fundamental 00th wave harmonic reflected by the interface is retained in Eq. (28) and where Γ denotes the field reflection coefficient given by

B. Field Inside the Metamaterial Half-Space and Relation to the Ewald-Oseen Extinction Principle
On the other hand, when observing the scattered field inside the metamaterial half-space, z > 0 one has n ≥ 0 so that the integrand of Eq. (26) is regular at ζ ¼ 0. Therefore, Eq. (26) is now calculated via the Cauchy theorem as the sum of the residues at the poles inside C. Now, the integrand exhibits pole singularities at ζ mode q , which are the poles of PðζÞ, which lead to the q-indexed modes in the crystal where k mode rst;q ¼ k x;rx þ k y;sŷ þ ðk mode z;q þ 2πt=cÞẑ denotes the wave vector of the rstth 3D crystal Floquet harmonic and denotes its weight. Harmonics are represented by three indexes according to a 3D periodicity. The singularities at ζ rs , except for rs ¼ 00, are removable because the pole singularities ζ rs of G layer xx ðζ; x; y; z 0 Þ inside the unit circle C are canceled out by the zeros PðζÞ has at the same points. Indeed, Pðζ rs Þ ¼ 0, except for rs ¼ 00 for which Pðζ 00 ¼ ζ inc Þ ¼ −j2abϵE inc 0 =k, so that only the pole of G layer xx ðζ; x; y; z 0 Þ at ζ inc gives a contribution to Eq. (26). Hence, the scattered field (26) is evaluated as It is worth noting that the first term of the above expression of the field E scat x scattered in the metamaterial half-space exactly contains the negative of the incident field, in accordance with the Ewald-Oseen extinction principle [19,20]. Consequently, the total field E x ðrÞ ¼ E inc x ðrÞ þ E scat x ðrÞ in the metamaterial (i.e., the sum of the incident and the scattered field) only comprises the superposition of crystal modes, i.e., E x ðrÞ ¼ P q C q E mode x;q ðrÞ, where C q expresses the complex amplitude of each field crystal mode launched at the interface.
Consistent with [1], we have here rigorously proved that the field inside the half-space crystal (or metamaterial) is provided in terms of crystal modes, each of which is represented in terms of spatial harmonics [Eq. (31)]. We would like to stress that the exact incident field cancellation by the scattered field is even a more general property than what is shown here and it happens in a more complicated fashion in other structures. For example, in [25,26], it has been shown that, when a planar array is excited by a nearby localized source (not by plane wave like in this paper), the scattered field along the array plane, provided by the array itself, cancels the incident (direct) field provided by the localized source. Indeed, in the geometry in [25,26], where a localized source excites a planar array, the electromagnetic modes in the array cannot cancel the incident field because these are two different wave species with different phase velocities and exponential/algebraic decays. Therefore, besides the array modes, a different wave species, called "spatial wave" or "space wave" is generated in terms of Floquet harmonics, whose fundamental one cancels the incident field along the array plane.
In Section 5, some examples illustrate the presented discretized WH procedure and obtained electric field solutions for various crystal geometries.

ILLUSTRATIVE EXAMPLES
In all the illustrative examples that follow, we consider silver nanospheres with various radii and lattice constants. The silver permittivity is described by the Drude relative dielectric function ϵ m ¼ ϵ ∞ − ω 2 p ½ωðω þ iγÞ −1 , where ϵ ∞ ¼ 5 is the background permittivity of the metal, ω p ¼ 1:37 × 10 16 rad=s is the plasmon radian frequency, and γ ¼ 27:3 × 10 12 s −1 is the damping frequency accounting for metal losses [6,27]. This parameterization provides a reasonably accurate description of the dielectric properties of silver across the optical range. In Eq. (1) and subsequent formulas, we use the Mie expression of the nanosphere electric polarizability α, shown in [6,27].
A. Example 1 First, we consider a half-space filled with a cubic lattice a ¼ b ¼ c ¼ 100 nm of silver nanospheres with radius of 30 nm in free space, illuminated by a unit plane wave E inc 0 ¼ 1 V=m at frequency f ¼ 600 THz.

Factorization
In Fig. 2, the Z-transform function AðζÞ is plotted in the ζ complex plane; the amplitude is expressed in decibels (left), while its phase is expressed in radians (right). Since the lattice period is about λ=5, only one Floquet wave couple is propagating away from each layer, represented by poles at ζ 00 ¼ ζ inc , and 1=ζ 00 ¼ ζ ref close to the unit circle C; all the others are hardly visible close to the origin. Analogously, only the first couple of zeros, ζ mode 0 and 1=ζ mode 0 (representing the modes in the metamaterial half-space), are close to C, while all the other are close to the origin. In Fig. 3, the factorizing function A þ ðζÞ is plotted in the ζ complex plane. Now only singularities and zeros inside C are present while the function is regular outside C.
Next, in Fig. 4, the Z transform PðζÞ of the dipole moment sequence is plotted in the ζ complex plane. Namely, we plot PðζÞζ −1 , which has to be integrated to calculate the interface layer dipole moment p 0 via Eq. (22). For the calculation of inner layer dipole moments p n , the ζ n factor will add an nth-order zero at the origin. Note that the plotted amplitude of the dimensional quantity PðζÞζ −1 has been normalized to the polarization of an isolated nanosphere under the same illumination αE inc 0 . It is apparent that PðζÞ exhibits only pole singularities inside the unit circle. The dominant one is at ζ mode 0 , corresponding to the medium dominant mode, while all the others are strongly evanescent, very close to the origin, and hardly visible even for PðζÞζ −1 because they are nearly cancelled by zeros at ζ rs , also very close to the origin. The dominant mode wavenumber, via k mode z;0 ¼ n eff k, corresponds to an effective refraction index n eff ¼ 1:283 − j0:003, corresponding to an artificial dielectric with small losses.

Polarization of Particles
Finally, in Fig. 5, the profile of the induced dipole moment in the layers is plotted versus the layer number. The WH solution (red thick line) is compared against a purely numerical method of moment (MoM) solution of Eq. (14) for a large (N ¼ 2000) but finite number of layers, and the agreement is almost perfect; a hardly noticeable deviation reveals an almost negligible stationary wave in the numerical solution due to a reflection at the second interface. Such discrepancy might be further reduced by increasing the number of layers considered, thanks to the presence of losses. The linear phase profile (left) reveals the presence of a propagating mode inside the metamaterial (for better readability the phase has been unwrapped), which slightly decays exponentially because of losses, as appears from the amplitude profile (right).

B. Example 2
As a second example, we consider a smaller period cubic lattice a ¼ b ¼ c ¼ 40 nm of silver nanospheres with radius of 10 nm, embedded in a silica background ϵ h ¼ 2:2, illuminated by a unit plane wave Now the topology of the critical points of AðζÞ is similar but the metamaterial mode wavenumber k mode z;0 is smaller than the ambient wavenumber k, as it appears by the reciprocal position of ζ mode 0 and ζ 00 in Fig. 6 (zeros of AðζÞ are on the right of the poles). Indeed, the metamaterial-dominant mode exhibits an effective relative refraction index n eff ¼ n 0 eff − jn 00 eff ¼ 0:476 − j0:036 with n 0 eff < 1 and significant attenuation, as confirmed by the polarization profile in Fig. 7. Again, the solution is successfully checked against a numerical reference for which now just the N ¼ 200 layer is sufficient because of higher mode attenuation.

C. Example 3
The same configuration in Example 2 (Subsection 5.B) is now analyzed at f ¼ 735 THz, which is at the edge between a passband (at higher frequencies like in Example 2) and a stopband (at lower frequency as in Example 4). Now, the spectral point ξ 0 associated with the metamaterial-dominant mode would approach the point ζ ¼ 1, if the occurrence of strong losses had not displaced it deeper inside the unit circle (Fig. 8). Assuming that k mode z;0 ¼ n eff k, the resulting effective relative refraction index n eff ¼ n 0 eff − jn 00 eff ¼ 0:110 − j0:270 reveals an extreme n 0 eff ≈ 0 behavior though associated with important attenuation, as revealed by the considerable decaying rate of the layer polarizability (Fig. 9).

D. Example 4
The same structure analyzed in Examples 2 and 3 is also considered here. When decreasing the frequency, the spectral point ξ 0 bends onto the real axis and move toward ζ ¼ 0. At f ¼ 720 THz, the ζ complex plane topology is shown in         10. Now the metamaterial mode is dramatically decaying with an effective relative refraction index n eff ¼ n 0 eff − jn 00 eff ¼ 0:0989 − j0:8250 that is almost purely imaginary since the metamaterial mode is in its cut-off regime. Figure 11 shows the nanosphere dipole moment versus layer number: dipoles exhibit a strong exponential decay and a very small phase variation moving away from the interface.

E. Example 5
In the previous examples, the wave dynamics inside the metamaterial was dominated by just one mode. In such cases the metamaterial behaves very similarly to a homogeneous medium that supports a transmitted plane wave and no peculiar wave dynamics is observed. To illustrate a different wave dynamics, we now consider a parallelepiped lattice with a ¼ b ¼ 90 nm and c ¼ 30 nm of silver nanospheres with radius of 10 nm, embedded in a silica background ϵ h ¼ 2:2, illuminated by a unit plane wave E inc 0 ¼ 1 V=m, at frequencies f ¼ 720 and 735 THz. Compared to the previous cases, one can readily notice in Fig. 12 that AðζÞ has more critical points. Two modes are present, corresponding to ζ mode 0 and ζ mode 1 , associated with modal wavenumbers k mode z;0 and k mode z;1 , respectively, provided by the mapping Eq. (16). Both the zeros ζ mode 0 and ζ mode 1 have magnitudes smaller than the pole ζ 00 , associated with the background wavenumber k, implying that both modes attenuate exponentially away from the interface. However, while at f ¼ 735 THz, jζ mode 1 j ≪ jζ mode 0 j ≈ 1; i.e., the higher-order mode ζ mode 1 rapidly decays and only the dominant mode ζ mode 0 is significantly propagating inside the metamaterial, at f ¼ 720 THz, jζ mode 1 j ≈ jζ mode 0 j < 1 and the two modes contribute to the field inside the metamaterial. Figure 13 shows the dipole moments versus layer number at two different frequencies. Note that at f ¼ 720 THz (Fig. 13, left) the magnitude profile is irregular because of the interference of the two modes, which have comparable attenuation rates. Instead, at f ¼ 735 THz (Fig. 13, right) only one mode is dominant. Another important observation is that, at f ¼ 735 THz, even when the profile seems regular, the simple exponential decay behavior excludes the first layer (n ¼ 0) and therefore a simple homogenization theory would not work, and more sophisticated "transition layers" may be necessary, as, for example, hypothesized in [2,3]. Again, there is perfect agreement with the WH solution and the MoM, as expected.

CONCLUSION
Wave dynamics originated by the interface between a homogeneous half-space and a half-space metamaterial made of arrayed plasmonic nanospheres is studied in this paper for the case of plane wave orthogonal incidence. The solution, obtained with a discrete WH method, is cast in terms of an exact analytic representation of the dipole moment induced in the nanospheres and of the field at any point. The field inside the metamaterial is expressed as a superposition of the modes calculated as though it was unbounded. It is shown that more than one mode can be excited at the interface and propagates inbound the metamaterial. The WH method determines the modal wavenumbers and the modal launching coefficients at the interface. The method can be generalized to a variety of other cases, including oblique incidence. While this paper is devoted only to the fundamental aspects of this new approach for studying waves originated by the interface with a halfspace metamaterial, the exact analytic formulation may also lead to deeper understanding of the wave dynamics in these half-space problems. A future paper will be devoted to the skew incidence case, investigating refraction, anisotropy, and birefringence there occurring.

APPENDIX A
The (discrete Laplace) Z transform of a sequence f n is defined as [11,23,24] FðζÞ ¼ and can be inverted by where C denotes a counterclockwise closed path encircling the origin and entirely in the region of convergence. To calculate the Z transform AðζÞ in Eq. (18), the generic term of the sequence a n , for n ≠ 0, is expressed through Eqs. (11) and (8) as a n ¼ jα 2abϵ while a 0 ¼ 1 − αĞ layer xx ð0Þ. Therefore, by splitting the definition Eq. (A1) for A into the negative and the positive semi-infinite series and exploiting the symmetry a n ¼ a −n , Fig. 13. (Color online) Particle dipole moments p n in the metamaterial. Amplitude (normalized to αE inc 0 , decibels, left) and unwrapped phase (radians, right). The WH solution (red circles) is compared against a numerical MoM solution for a large but finite number of layers N ¼ 20 (black dots).

APPENDIX B
Traditionally, A AE ðζÞ is defined as where C þ (C − ) denotes a counterclockwise (clockwise) path on the unit circle that is suitably deformed to leave outside (inside) the pole at ξ ¼ ζ if jζj ≤ 1 (jζj ≥ 1). Note that A AE factorize both poles and zeros of A, which are all singularities of ln AðζÞ. In [11], where the reader is referred for more details, an alternative definition is proposed, which is more suitable for numerical integration since the integration path is on the unit circle independent of the location of ζ, as the pole singularity is removed. As a matter of fact, the integrand in Eq. (B2) is intended to be analytically continued to its limit value d dζ ln AðζÞ þ 1 2ζ AðζÞ at ξ ¼ ζ. Furthermore, Eq. (B2) automatically imposes A − ðζÞ ¼ A þ ð1=ζÞ when Að1=ζÞ ¼ AðζÞ by properly setting an arbitrary multiplication constant in the factorization Eq. (19).