## Abstract

We characterize the modes with complex wavenumber for both longitudinal and transverse polarization states (with respect to the mode traveling direction) in three dimensional (3D) periodic arrays of plasmonic nanospheres, including metal losses. The Ewald representation of the required dyadic periodic Green’s function to represent the field in 3D periodic arrays is derived from the scalar case, which can be analytically continued into the complex wavenumber space. We observe the presence of one longitudinal mode and two transverse modes, one forward and one backward. Despite the presence of two modes for transverse polarization, we notice that the forward one is “dominant” (i.e., it contributes most to the field in the array). Therefore, in case of transverse polarization, we describe the composite material in terms of a homogenized effective refractive index, comparing results from (i) modal analysis, (ii) Maxwell Garnett theory, (iii) Nicolson-Ross-Weir retrieval method from scattering parameters for finite thickness structures (considering different thicknesses, showing consistency of results), and (iv) the fitting of the fields obtained through HFSS simulations. The agreement among the different methods justifies the performed homogenization procedure in case of transverse polarization.

© 2011 OSA

## 1. Introduction

Plasmonic resonances can occur in metallic structures with subwavelength dimensions. This property enables the design of aggregates and arrays of nanoparticles, periodic in one (1D), two (2D) or three dimensions (3D), that possess collective resonances, with subwavelength distance between the constitutive nanoparticles. These resonances can ultimately be used to guide backward modes, create artificial dielectrics, narrow band absorption, artificial magnetism, quasi-dark modes. Therefore, a thorough understanding of the modes that can be excited in a periodic array of nanoparticles is essential to relate them to physical employment for innovative applications. We are interested in the case of periodic arrays of plasmonic nanospheres. Modal analysis in 1D and 2D periodic arrays of plasmonic nanospheres have been studied in detail in [1–8] and references therein, also providing the criteria of physical existence of the excitable modes in such periodic structures [6,7]. However, we have not found in literature a complete and exhaustive characterization description of the modes in 3D periodic arrays of plasmonic nanospheres. Therefore, the interest of this paper concerns modal analysis in such a 3D periodic array, and in particular how waves propagate inside the metamaterial. Importantly, 3D periodic arrays of nanospheres can be engineered to obtain peculiar characteristics, such as, among others, slow wave structures, and double negative materials. Mode analysis of 3D periodic arrays of metallic nanospheres have been performed in [5,9]. In [9], the authors adopted the single dipole approximation (SDA) [10,11] with the quasistatic nanosphere polarizability, according to the Clausius Mossotti approach, to design broadband effective negative index metamaterials by using the nano-transmission line network concept, and focusing mostly on the transverse polarization (with respect to the direction of propagation of the modes) which is the only one that allows for backward propagation. They considered densely packed and properly designed 3D periodic arrays because they are capable of supporting a nano-transmission line propagating mode and thus acting as an effective negative index metamaterial. However, a complete physical modal description was out of the scope of that paper and little attention has been given to the longitudinal polarization (since it does not allow for backward propagation). In [5], the authors analyzed the modes with complex wavenumber in 3D periodic arrays of silver nanospheres, accounting also for metal losses. Moreover, the authors clearly stated the difficulty of finding all possible modal solutions, which may result in missing some branches of the dispersion diagram. They reported the dominant mode for transverse polarization (see Sec. 3 where we report also a backward mode, and Sec. 6 where we show that this is not dominant), and indeed, their objective was primarily to display a representative selection of *k−β* diagrams of modal solutions (where *k* is the host medium wavenumber, and *β* is the wavenumber of the guided wave) without focusing on the physical/nonphysical properties of the found modes.

The approach described in the present paper allows for the tracking and especially for the characterization of the evolution of modes, varying frequency. We analyze optical modes with complex wavenumber in 3D periodic arrays of plasmonic nanospheres for both longitudinal and transverse (with respect to the mode traveling direction) polarization states. Each nanosphere is modeled to act as a single electric dipole, by using the SDA [10,11] and the metal permittivity is described by the Drude model. The numerical procedure developed in this paper for evaluating the complex zeros of the dispersion relation uses the Ewald representation for the *dyadic* periodic Green’s function (GF) to represent the field in 3D periodic arrays, and is based on previous scalar developments [12–14]. The Ewald representation, besides providing analytic continuation to the complex wavenumber space, results in two series with Gaussian convergence where only a handful of terms are needed [12–14].

As it will be shown in Sec. 6 dealing with a metamaterial with finite thickness, despite the presence of two modes for transverse polarization, only one is “dominant”, in the sense that it contributes mostly to the field in the array. Therefore, in case of transverse polarization, we describe the composite material in terms of homogenized effective refractive index. This is obtained by the modal analysis for transverse polarization for a cubic lattice metamaterial, is compared (i) to the one obtained from Maxwell Garnett (MG) homogenization theory, (ii) to that retrieved by scattering parameters of finite thickness structures by using the Nicolson-Ross-Weir (NRW) method, and (iii) to that retrieved directly by fitting the fields obtained through HFSS simulations. The agreement among the different methods confirms the validity of the homogeneous treatment of the composite material. We also perform parametric analyses and present a structure that is able to boost the figure of merit (defined as the ratio between the real and the imaginary part of the modal complex wavenumber).

The structure of the paper is as follows. In Sec. 2, a brief description of the simulation model to perform modal analysis in 3D periodic arrays of plasmonic nanospheres is reported. The detailed computation of the dispersion diagrams for longitudinal and transverse polarization states (with respect to the mode traveling direction) for a set of array parameters is then reported in Sec. 3, and a comparison of the performance of three different structures (parametric analysis) is in Sec. 4. We then report a comparison of reflection, transmission and absorption due to a normal incident plane wave excitation obtained using SDA with those obtained by a full-wave HFSS simulation pertaining to ten layers of arrayed plasmonic nanospheres, stacked in the direction of propagation, in Sec. 5. Finally, in Sec. 6, homogenization theory is adopted to compute the effective refractive index and the results from modal analysis are compared to the ones obtained from Maxwell Garnett theory, as well as to those retrieved from HFSS simulation either by using NRW method or directly by the fields obtained through HFSS simulations. In addition, the evaluation of the Ewald representation for the dyadic GF for 3D periodic arrays is derived for the first time, to the authors’ knowledge, in the Appendix A.

## 2. Simulation model

The structure under analysis is the 3D periodic array of metallic nanospheres reported in Fig. 1 . The monochromatic time harmonic convention, $\mathrm{exp}\left(-i\omega t\right)$, is assumed here and throughout the paper, and is therefore suppressed hereafter.

We model each nanosphere as a single electric dipole. As such, for a metallic spherical particle the induced dipole moment is

where ${\alpha}_{\text{ee}}$ is the isotropic electric polarizability of the nanosphere, ${E}^{\text{loc}}$ is the local field produced by all the nanospheres of the array except the considered nanosphere, plus the external incident field to the array, if present. Bold fonts refer to vector quantities. The SDA is a good approximation when the metallic nanospheres are used close to their fundamental plasmonic frequency, when particle dimensions are much smaller than the wavelength, and when the edge-to-edge spacing*d*between spheres is larger than the spheres’ radius

*r*(i.e., $d\ge r$). However, even for smaller distances the SDA may provide satisfactory approximated results [15]. In general, for a spacing between the spheres smaller than their radius (i.e., $0<d<r$), more accurate results would involve multipole field contributions [10,16–18]. According to Mie theory, the electric polarizability for a nanosphere is [10,11]

Consider now a 3D periodic array of nanospheres as in Fig. 1, immersed in a homogeneous background, with relative permittivity ${\epsilon}_{h}$, for which each spherical nanoparticle is placed at positions ${r}_{n}={r}_{0}+{d}_{n}$, where $n=\left({n}_{1},{n}_{2},{n}_{3}\right)$ is a triple index, and ${d}_{n}={n}_{1}a\widehat{x}+{n}_{2}b\widehat{y}+{n}_{3}c\widehat{z}$ (where a caret on top of a bold letter refers to unit vector quantities), with ${n}_{1},{n}_{2},{n}_{3}=0,\pm 1,\pm 2,\mathrm{...}$, ${r}_{0}={x}_{0}\widehat{x}+{y}_{0}\widehat{y}+{z}_{0}\widehat{z}$, and *a*, *b* and *c* are the periodicities along *x*-, *y*- and *z*-direction, respectively [11,21].

Suppose that the array is either excited or a mode (a periodic field) is present, with wavevector ${k}_{\text{B}}={k}_{x}\widehat{x}+{k}_{y}\widehat{y}+{k}_{z}\widehat{z}$. Consequently, each nanosphere will have an electric dipole moment equal to ${p}_{n}={p}_{0}{e}^{i{k}_{\text{B}}\cdot {d}_{n}}$ (the polarization direction is fixed for symmetry reasons). Then, the local electric field acting on a nanosphere at position ${r}_{0}$ is given by

Mode analysis in the 3D periodic array is performed by computing the complex zeroes of the homogeneous version of Eq. (6); i.e., when no impressed excitation is present, or ${E}^{\text{inc}}\left({r}_{0}\right)=0.$ This requires the solving of $\mathrm{det}\underset{\xaf}{A}\left({k}_{\text{B}}\right)=0$ for complex ${k}_{\text{B}}$. For a given angular frequency $\omega $, imposing the determinant of Eq. (6) to vanish, implicitly defines a dispersion relation among the three complex components of ${k}_{\text{B}}$. Once two independent components of ${k}_{\text{B}}$, say ${k}_{x}$ and ${k}_{y}$, are given, the third, say ${k}_{z}$, is imposed by the dispersion relation. For general isotropic media the dispersion relation is invariant with respect to an arbitrary rotation, and the same mode can propagate in any direction; therefore the modal analysis can be carried out in a specific direction which will be representative of all. However, in the general anisotropic case, the analysis must be carried out for any direction. Also, because of crystal symmetry and electromagnetic reciprocity, a mode which travels in a given direction can also travel in the opposite direction. The mathematical counterparts of this physical property resides in the even dependence $\underset{\xaf}{A}\left({k}_{\text{B}}\right)=\underset{\xaf}{A}\left(-{k}_{\text{B}}\right)$, for which if ${k}_{\text{B}}$ is a solution of $\mathrm{det}\underset{\xaf}{A}\left({k}_{\text{B}}\right)=0$, also $-{k}_{\text{B}}$ is. Therefore, modes are always found in couples travelling in opposite directions; however they can be easily distinguished by assuming a physical decay along the energy propagation direction. In absence of losses, instead, a four-quadrant symmetry holds, thus the correspondent complex conjugate solutions, $\pm {k}_{\text{B}}^{*}$, would be solutions as well.

## 3. Dispersion diagrams for 3D periodic arrays of plasmonic nanospheres

In this section, we analyze modes with complex wavenumber traveling along the *z*-direction, i.e., ${k}_{\text{B}}={k}_{z}\widehat{z}$, with ${k}_{z}={\beta}_{z}+i{\alpha}_{z}$, in 3D periodic arrays of silver nanospheres embedded in free space (i.e., $k={k}_{0}$), with dipole moments polarized longitudinally (L-pol) or transversally (T-pol) to the traveling direction *z*. Namely, we impose ${k}_{x}={k}_{y}=0$ for any given frequency and we solve $\mathrm{det}\underset{\xaf}{A}\left({k}_{\text{B}}\right)=0$ for complex ${k}_{z}$. Between the two possible opposite solutions, the mode with power flow toward the positive *z*-direction is selected by choosing ${\alpha}_{z}\ge 0$. The radius of the spherical nanoparticles is $r=25\text{nm}$. The adopted Drude model parameters for silver permittivity are ${\epsilon}_{\infty}=5$, ${\omega}_{p}=1.37\times {10}^{16}\text{rad/s}$ and $\gamma =27.3\times {10}^{12}{\text{s}}^{-1}$ [22,23], to fit measured data in the analyzed frequency range [24]. Based on this data, an isolated particle would exhibit a strong dipole moment around 800 THz ($ka/\pi =0.4$ assuming $a=75\text{nm}$). Modes are also characterized in terms of their polarization state by calculating the eigenvector $p$ associated to a wavenumber ${k}_{z}$. It has to be noticed that when considering propagation along a principal direction of the crystal, say *z*, the matrix $\underset{\xaf}{A}\left({k}_{\text{B}}\right)$ in Eq. (6) will present a block structure with ${A}_{xz}={A}_{zx}={A}_{yz}={A}_{zy}=0$ because of symmetry which decouples the longitudinal and the transverse polarizations. Therefore, in the L-pol case, each nanosphere has an induced dipole moment along the *z*-direction (i.e., $p={p}_{z}\widehat{z}$). In addition, when the transverse periodicities *a* and *b* along *x* and *y* are equal, yet because of symmetry, one has ${A}_{xx}={A}_{yy}$ and ${A}_{xy}={A}_{yx}$; therefore, in the T-pol case, each nanosphere has an induced dipole moment orthogonal to *z*, i.e., $p={p}_{x}\widehat{x}+{p}_{y}\widehat{y}$. Modal analysis is then performed for the three 3D periodic arrays reported in Table 1
(although we report the dispersion diagrams only for Structure I), and results are compared in Sec. 4.

#### 3.1 Transverse polarization (T-pol)

The dispersion diagrams for the Structure I outlined in Table 1 are shown in Fig. 2
for both the real and the imaginary parts of the wavenumber ${k}_{z}={\beta}_{z}+i{\alpha}_{z}$ with respect to the host wavenumber $k$ in the case of T-pol. Only modes with ${\alpha}_{z}\ge 0$, i.e., those with power flow toward the positive *z*-direction, are shown.

Mode 1 (blue line) follows a typical dispersion curve which is almost straight at low frequencies corresponding to an effective medium slightly denser than the host medium (whose ideal dispersion, the light line, ${\beta}_{z}=k$ is plotted in black dash-dotted line) with small attenuation; then, increasing the frequency, the dispersion curve bends exhibiting large phase constant (close to the edge of the Brillouin zone $-1<{\beta}_{z}c/\pi <1$), thus this mode could be employed, for example, in field concentration and imaging applications. Further increasing frequency, Mode 1 experiences a bandgap with a strong attenuation; finally, at higher frequencies Mode 1 reenters a propagation band with small attenuation. Mode 2 (green line) at low frequencies is characterized by a large phase constant at the edge of the Brillouin zone and large attenuation constant. Increasing frequency, the attenuation constant decreases (remaining yet larger than that of Mode 1 in the bandgap) in a small frequency region, with phase constant approaching very small values. At higher frequencies, Mode 2 is characterized by very small phase constant and large attenuation constant. Other modes with normalized attenuation constant larger than ${\alpha}_{z}c/\pi =2$ are present but not reported here since they dramatically decay as ${\alpha}_{z}>>k$. The dispersion diagram assuming ideal lossless nanospheres is reported in Fig. 2 as black dotted lines, and it refers to a crystal comprising spheres made of hypothetical lossless silver with Drude parameter $\gamma \to 0$. It is interesting to note how small losses keep dispersion diagram branches separated for all frequencies, while they overlap in the lossless case in some frequency range (see Fig. 3 ), thus rendering the mode frequency evolution ambiguous in some cases. The occurrence of two transversely polarized (“doubly degenerated”, that is admitting any transverse polarization) modes reveals the non-local behavior of the crystal (i.e., the presence of spatial dispersion, as also indicated in [25], which also depends on how dense the array is [26]).

Next, in Fig. 3, we show the evolution of the modal wavenumber ${k}_{z}={\beta}_{z}+i{\alpha}_{z}$ for varying frequency. Namely the trajectories of the complex propagation constants are tracked in the complex plane ${\beta}_{z},\text{\hspace{0.17em}}{\alpha}_{z}$ in Fig. 3(a), accounting for silver losses, and also for the ideal lossless case. In the latter case, the four-quadrant symmetry discussed in Sec. 2 is present (dotted black curves). In Fig. 3(b) we show the trajectories in the complex plane ${\beta}_{z},\text{\hspace{0.17em}}{\alpha}_{z}$, normalized to that of the host medium $k$, accounting for silver losses. When silver losses are accounted for, the symmetry with respect to the origin in Figs. 3(a) and 3(b) is a consequence of the property $\underset{\xaf}{A}\left({k}_{\text{B}}\right)=\underset{\xaf}{A}\left(-{k}_{\text{B}}\right)$ discussed in Sec. 2; in other words, the four-quadrant symmetry has been broken by the presence of losses. Notice that a mode whose power flow is toward either the positive or negative *z*-direction, for which either ${\alpha}_{z}>0$ or ${\alpha}_{z}<0$, is forward (from the phase progression point of view) when ${\alpha}_{z}{\beta}_{z}>0$, whereas it is backward when ${\alpha}_{z}{\beta}_{z}<0$. By looking at both Figs. 2 and 3, the forward Mode 1 and the backward Mode 2 could be guided in the structure in the lossy case.

Note that the presence of two transverse modes with moderately low attenuation constant ${\alpha}_{z}$ is in agreement with what previously predicted in [9] by using the nano-transmission line network concept, and analytically in [25] for an ideal lossless case.

#### 3.2 Longitudinal polarization (L-pol)

The dispersion diagrams for the Structure I outlined in Table 1 are shown in Fig. 4
for both the real and the imaginary parts of the wavenumber ${k}_{z}={\beta}_{z}+i{\alpha}_{z}$ in the case of L-pol (longitudinal with respect to the mode traveling direction *z*). Only the modes with ${\alpha}_{z}\ge 0$, i.e., whose power flow is toward the positive *z*-direction, are shown. Notice again that other modes with normalized attenuation constant larger than ${\alpha}_{z}c/\pi =2$ are present but not reported here since guided modes can travel a significant distance only when their attenuation constant is small, or ${\alpha}_{z}<<k$. The case relative to ideal lossless silver nanospheres is also reported in Fig. 4 as a black dotted line. The almost flat portion of the curve in Fig. 4(a) reveals a certain degree of non-locality, and thus spatial dispersion.

The evolution of the modal wavenumbers varying frequency is shown in the complex ${k}_{z}$ plane in Fig. 5 . In the ideal lossless case, the four-quadrant symmetry discussed in Sec. 2 holds also in this L-pol case as shown in Fig. 5(a). At low frequencies, Mode 1 (blue curve) is characterized by very small phase constant and large attenuation constant. Increasing frequency, the attenuation constant decreases in a small frequency region reaching a minimum of ${\alpha}_{z}c/\pi \approx 0.19$ at $kc/\pi \approx 0.43$. At higher frequencies, both phase and attenuation constants increase. The presence of the longitudinal mode is associated to near field dipole couplings among nanospheres.

Note again that the presence of one longitudinal mode with moderately low attenuation constant ${\alpha}_{z}$ is in agreement with what previously predicted in [9,25].

## 4. Parametric analysis: comparison of modal dispersion diagrams

#### 4.1 Transverse polarization (T-pol)

We compare in Fig. 6 the real and the imaginary parts of the wavenumber ${k}_{z}$ of Mode 1, with T-pol, computed for the three arrays in Table 1.

It can be easily noticed that increasing the periodicities in the directions orthogonal to the mode traveling direction (i.e., Structure III) highly reduces the wavenumber imaginary part, thus allowing the mode to travel longer distances for some frequency ranges, as explained in Sec. 4.3. For the case of Structure II, the phase propagation constant reaches the edge of the BZ (this does not happen for the two other structures). Furthermore, the bandgap is wider for Structure II.

#### 4.2 Longitudinal polarization (L-pol)

We compare in Fig. 7 the real and the imaginary parts of the wavenumber ${k}_{z}$ for Mode 1 (L-pol) computed for the three arrays in Table 1.

As in the previous case, it is noticed that increasing the periodicities in the directions orthogonal to the mode traveling direction (i.e., Structure III) highly reduces the imaginary part, thus allowing the mode to travel longer distances in a narrow frequency range, as explained in Sec. 4.3 (we also observe a slight frequency shift towards lower frequencies). Mode 1 in Structure II exhibits very high attenuation at all frequencies.

#### 4.3 Figure of merit

In accord to the modes shown in Secs. 4.1 and 4.2, we are here interested in quantifying how far a mode can travel in terms of *guided* wavelengths. Therefore, it is convenient to define the figure of merit

The higher the figure of merit *F* is, the longer the traveled distance in terms of guided wavelength is, relatively to the attenuation. We show in Fig. 8
the figure of merit versus frequency of the transverse and longitudinal cases analyzed in Secs. 4.1 and 4.2, comparing *F* of the Mode 1 for both T-pol and L-pol, for the three structures in Table 1. In both polarizations, the modes that show the best performance are the ones in Structure III in Table 1, for which $a=b=2c$: for the T-pol case, the peak of *F* is about $7\times {10}^{4}$, whereas it is about 19 for the L-pol case. Also notice how the figure of merit for Structure III for T-pol is large for a wide frequency range analyzed (i.e., about 100-650 THz), whereas for L-pol is large only for a narrow frequency range (i.e., 790-850 THz). In Structure III, modes with T-pol and L-pol can travel longer distances in terms of guided wavelengths than in the other two arrays, as can be inferred by the graphs in Fig. 8.

## 5. 3D lattice with finite thickness along z: transmission, reflection and absorption

We analyze propagation in Structure I in Table 1 (Fig. 1) made of 10 layers of silver nanospheres stacked along the *z*-direction (Fig. 9
), by using (i) HFSS and (ii) the SDA with the Mie polarization expression in Eq. (2). The stack is illuminated by a normally incident plane wave traveling toward + *z*, and the magnitude of transmission *T* and reflection *R* of the stack are shown in Fig. 9, together with the absorption $A=1-{\left|T\right|}^{2}-{\left|R\right|}^{2},$ evaluated in dB as 10Log *A*. Results in Fig. 9 show good agreement between the HFSS full-wave simulation and the SDA theoretical results (also the phase of *R* and *T* obtained from the two methods are in good agreement, not shown). The strong stop-band between 730 THz ($kc/\pi =0.365$) and 890 THz ($kc/\pi =0.445$) is due to the large imaginary part of the attenuation constant ${\alpha}_{z}$ of Mode 1 shown in Fig. 2. Notice also that absorption *A* is particularly strong at the edges of the stop-band. In the evaluation of the SDA results we have used the Ewald representation of the dyadic GF for 2D periodic arrays reported in [27].

## 6. Homogenization: effective refractive index

In general, considering only modes with a moderately low attenuation constant ${\alpha}_{z}$, a plane wave impinging on the composite material in Fig. 9 could excite any of the permitted modes, i.e., two modes with T-pol and one mode with L-pol, with different amplitudes (as also briefly discussed in [9]). However, for a normal incident plane wave, only modes with T-pol could be excited. We will show in Sec. 6.2 that even though the analyzed structure has a certain degree of spatial dispersion (as mentioned in Sec. 3), Mode 1 (T-pol) in Sec. 3.1 is “dominant”. Accordingly, a wave in the composite material, in the case of transverse polarization, could be described with good approximation as a TEM wave in a homogeneous material, which can in turn be represented by an effective refractive index. Though, strictly speaking, the effective permittivity characterizing the crystal as a homogenized material has a small degree of non-locality (spatial dispersion) [9,25].

#### 6.1 NRW retrieval method and comparison with theoretical results

Transmission and reflection coefficients for a stack of layers are here used to retrieve the effective refractive index of the metamaterial by using the Nicholson-Ross-Weir (NRW) method [28–33]. Treating the composite slab as a uniform continuous medium with same thickness *t*, according to NRW, the complex effective refractive index can be retrieved by

*q*and +/− in Eq. (8).

The comparison of the retrieved effective refractive index of Structure I is then performed using three different methods: (i) by normalizing the modal wavenumber obtained via SDA shown in Sec. 3, as ${n}_{\text{eff}}={k}_{z}/{k}_{0}$, (ii) by using the Maxwell Garnett formulas [34,35] with the polarization expression from Mie theory (the comparison with the quasistatic polarizability expression (Clausius Mossotti) has been reported in [36]), and (iii) by using the NRW algorithm in Eq. (8), for *N* layers of arrayed metallic nanospheres, stacked along the *z*-direction, with *T* and *R* computed via HFSS simulation. The three methods provide comparable results as shown in Fig. 10
(where *N* = 10). Note that large values of refractive index are obtained near *f* = 750 THz ($kc/\pi =0.375$), as well as very small ones (close to zero) near *f* = 840 THz ($kc/\pi =0.42$), showing that this material, besides providing small guided wavelengths, can also provide an epsilon near zero (ENZ) artificial medium in a narrow frequency band. Losses can be rather low, in the high refractive index region, as also shown by the figure of merit in Fig. 8, and in the ENZ region.

In Fig. 11
we show the results obtained by using the NRW method (with HFSS simulations) for *N* = 4 and *N* = 6 layers of metallic nanospheres stacked along the *z*-direction in comparison to the result in Fig. 10. Observe that both real and imaginary parts agree well, although the imaginary part of the NRW-HFSS method is smaller around the bandgap for increasing number of layers. This is indeed in the region where $\left|T\right|$ is close to 0. In this region, the HFSS results show some discrepancies probably due to the sensitivity of Eq. (8) to small perturbations of $\left|T\right|$, as $\left|T\right|\to 0$ [37] (which are mainly due to the HFSS numerical precision). To understand if this is true, in Fig. 12
we show the calculations for the same stack of layers by using the SDA, which is more stable numerically than HFSS, especially for small $\left|T\right|$. Results from the NRW-SDA method are in good agreement to those obtained with mode analysis for both 4 (where the lowest $\left|T\right|\approx -70\text{dB}$) and 6 layers (where the lowest $\left|T\right|\approx -110\text{dB}$) (in Fig. 11 the simulation with 6 layers was already in disagreement). However, in the case of 10 layers, in the frequency region of extremely low value of $\left|T\right|$, also the NRW-SDA provides inaccurate results in the bandgap, due to the limited SDA numerical precision in the numerical algorithm for the estimation of the transmission coefficient and to the high sensitivity of Eq. (8) explained above.

#### 6.2 Field fitting retrieval method

From the HFSS full-wave simulation of 10 layers of plasmonic nanospheres shown in Sec. 5, illuminated by a normally incident plane wave traveling toward +*z*, with electric field polarized along *y*, we extract the *y* component of the electric field (1 point per layer, at the center of each sphere) at two representative frequencies of the curves reported in Fig. 10: first at 745 THz ($kc/\pi =0.3725$), corresponding to large $\mathrm{Re}\left[{n}_{\text{eff}}\right]$ and low $\mathrm{Im}\left[{n}_{\text{eff}}\right]$, and then at 875 THz ($kc/\pi =0.4375$), corresponding to low $\mathrm{Re}\left[{n}_{\text{eff}}\right]$ and low $\mathrm{Im}\left[{n}_{\text{eff}}\right]$.

At each of the two frequency points the total field in each sphere is represented as the superposition of a direct (${E}_{+}$) and a reflected (${E}_{-}$) wave, pertaining to a *single mode* (T-pol) with complex wavenumber ${k}_{z}$, traveling along the ± *z* directions (see Fig. 9) as follows

*single mode*. Moreover, from the fitting, we extract ${k}_{z}/{k}_{0}\approx 2.15+0.07i$ at 745 THz (which is compared with the mode solution in Fig. 10 where ${k}_{z}/{k}_{0}\approx 2.13+0.07i$) and${k}_{z}/{k}_{0}\approx 0.40+0.07i$ at 875 THz (which is compared with the mode solution in Fig. 10 where ${k}_{z}/{k}_{0}\approx 0.40+0.07i$). A good agreement is then observed. Also, looking at the phase information in Figs. 13(b) and 14(b), it can be easily noticed that the phase accumulated by the wave at 745 THz after the 10 layers is much larger than the one at 875 THz, showing that the latter travels very slowly inside the structure, confirming that at this frequency we have an ENZ composite material. Moreover, the ratio between the reflected and the direct wave $\Psi \left(n\right)=\Psi \left(0\right){e}^{-2i{k}_{z}\left(n-{\scriptscriptstyle \frac{1}{2}}\right)c}$ (with $\Psi \left(0\right)={E}_{-}/{E}_{+}$) has been verified to be consistent for every

*n*, where $\Psi \left(0\right)=0.13{e}^{i0.078}$ at 745 THz and $\Psi \left(0\right)=0.06{e}^{i1.55}$ at 875 THz. As a consequence, we can conclude that Mode 1 for Structure I, reported in the dispersion diagrams in Fig. 2, is dominant.

## 7. Conclusion

This work provides a description of modal analysis in 3D periodic arrays of metallic nanospheres. We observe the presence of one longitudinal mode and two transverse modes, one forward and one backward. The performance of the guided Mode 1 for both T-pol and L-pol in the three arrays in Table 1 have been compared, and it has been noticed that the modes can travel longer distances in terms of guided wavelengths when the periodicities in the directions transversal to the mode traveling direction are larger than the periodicity along the mode traveling direction (i.e., Structure III in Table 1). In case of T-pol, Mode 1 has been found to be dominant, allowing for the description of the composite material in terms of homogenized effective refractive index. This has been indeed retrieved by four different methods, in good agreement with each other. The method shown in this paper can be applied to a variety of composite materials made of 3D periodic arrays of nanoparticles.

## Appendix A: Ewald representation for the dyadic GF for 3D periodic arrays

The 3D array regularized dyadic GF reported above Eq. (4) can be computed as

where ${\stackrel{\u2323}{G}}^{\infty}\left(r,{r}_{0},{k}_{\text{B}}\right)$ is the regularized scalar GF, not singular at $r={r}_{0}$, defined as

with $n=\left({n}_{1},{n}_{2},{n}_{3}\right)$ a triple index, ${d}_{n}={n}_{1}a\widehat{x}+{n}_{2}b\widehat{y}+{n}_{3}c\widehat{z}$, ${R}_{n}=r-{r}_{0}-{d}_{n}$, ${R}_{n}=\left|{R}_{n}\right|$, and $\nabla \nabla $ denotes the Hessian dyadic differential operator. According to the Ewald representation [38,39] the dyadic GF is split into the sum

of a spectral and a spatial dyadic part, which are obtained by using in Eq. (10) an analogous split of the scalar GF ${\stackrel{\u2323}{G}}^{\infty}\left(r,{r}_{0},{k}_{\text{B}}\right)={G}_{\text{spectral}}\left(r,{r}_{0},{k}_{\text{B}}\right)+{\stackrel{\u2323}{G}}_{\text{spatial}}\left(r,{r}_{0},{k}_{\text{B}}\right)$. The expressions of these latter two scalar terms, the spectral ${G}_{\text{spectral}}\left(r,{r}_{0},{k}_{\text{B}}\right)$ and the regularized spatial ${\stackrel{\u2323}{G}}_{\text{spatial}}\left(r,{r}_{0},{k}_{\text{B}}\right)$ can be found in [13,14], and in the particular case where $r={r}_{0}=0$, as required in Eq. (6), it follows that

where ${\gamma}_{n}^{2}={\left|{k}_{\text{B}}+{k}_{n}\right|}^{2}-{k}^{2}$, with ${k}_{n}=\left(2\pi {n}_{1}/a\right)\widehat{x}+\left(2\pi {n}_{2}/b\right)\widehat{y}+\left(2\pi {n}_{3}/c\right)\widehat{z}$, and

Here $\text{erfc}\left({\beta}^{\pm}\right)$ denotes the complementary error function of argument ${\beta}^{\pm}={R}_{n}E\pm ik/\left(2E\right)$, where $E$ is the Ewald parameter as in [40], which is in general chosen as

Also, in Eq. (14) a prime (${f}^{\prime}$) denotes a derivative of $f$ with respect to its argument.

Differentiating Eqs. (13) and (14) in rectangular and spherical coordinates, respectively, leads to the explicit representations for the Hessian dyads

where

with ${\widehat{R}}_{n}={R}_{n}/{R}_{n}$.

In Eqs. (19) and (18), a double and a triple prime (${f}^{\u2033}$ and ${f}^{\u2034}$) denote the second and third derivatives of *f* in Eq. (15) with respect to its argument, respectively.

As a last remark, we want to mention that there is a difference between the scalar and the dyadic GF in terms of number of elements needed for convergence in Eqs. (13)-(14) and (17)-(18): namely, we have observed that it was sufficient to use 11 terms for each index (i.e., ${n}_{1},{n}_{2},{n}_{3}=-5,\mathrm{...},5$) to obtain a converging dyadic GF, whereas 3 terms for each index (i.e., ${n}_{1},{n}_{2},{n}_{3}=-1,0,1$) were usually enough for the scalar GF to converge.

## Acknowledgment

The authors acknowledge partial support from the National Science Foundation, award NSF-CMMI #1101074, and from the European Commission FP7/2008, research area “NMP-2008-2.2-2 Nanostructured meta-materials”, grant “METACHEM”, no. 228762. The authors also thank Ansys for providing HFSS that was instrumental in this analysis.

## References and links

**1. **A. Alù and N. Engheta, “Theory of linear chains of metamaterial/plasmonic particles as subdiffraction optical nanotransmission lines,” Phys. Rev. B **74**(20), 205436 (2006).

**2. **R. A. Shore and A. D. Yaghjian, “Traveling waves on two- and three-dimensional periodic arrays of lossless scatterers,” Radio Sci. **42**(6), RS6S21 (2007).

**3. **R. Sainidou and G. F. de Abajo, “Plasmon guided modes in nanoparticle metamaterials,” Opt. Express **16**(7), 4499–4506 (2008).

**4. **D. Van Orden, Y. Fainman, and V. Lomakin, “Optical waves on nanoparticle chains coupled with surfaces,” Opt. Lett. **34**(4), 422–424 (2009).

**5. **R. A. Shore and A. D. Yaghjian, *Complex Waves on 1D, 2D, and 3D Periodic Arrays of Lossy and Lossless Magnetodielectric Spheres* (Air Force Research Laboratory, Hanscom, AFB, MA, 2010; revised, 2011).

**6. **S. Campione, S. Steshenko, and F. Capolino, “Complex bound and leaky modes in chains of plasmonic nanospheres,” Opt. Express **19**(19), 18345–18363 (2011).

**7. **A. L. Fructos, S. Campione, F. Capolino, and F. Mesa, “Characterization of complex plasmonic modes in two-dimensional periodic arrays of metal nanospheres,” J. Opt. Soc. Am. B **28**(6), 1446–1458 (2011).

**8. **S. Campione and F. Capolino, “Linear and planar periodic arrays of metallic nanospheres: Fabrication, optical properties and applications,” in *Selected Topics in Photonic Crystals and Metamaterials*, A. Andreone, A. Cusano, A. Cutolo, and V. Galdi, eds. (World Scientific, 2011), pp. 143–196.

**9. **A. Alù and N. Engheta, “Three-dimensional nanotransmission lines at optical frequencies: A recipe for broadband negative-refraction optical metamaterials,” Phys. Rev. B **75**(2), 024304 (2007).

**10. **C. F. Bohren and D. R. Huffman, *Absorption and Scattering of Light by Small Particles* (Wiley, 1983).

**11. **S. Steshenko and F. Capolino, “Single dipole approximation for modeling collections of nanoscatterers,” in *Theory and Phenomena of Metamaterials*, F. Capolino, ed. (CRC Press, 2009), 8.1.

**12. **P. Myun-Joo, P. Jongkuk, and N. Sangwook, “Efficient calculation of the Green's function for the rectangular cavity,” IEEE Microw. Guid. Wave Lett. **8**(3), 124–126 (1998).

**13. **I. Stevanoviæ and J. R. Mosig, “Periodic Green's function for skewed 3-D lattices using the Ewald transformation,” Microw. Opt. Technol. Lett. **49**(6), 1353–1357 (2007).

**14. **G. Lovat, P. Burghignoli, and R. Araneo, “Efficient evaluation of the 3-D periodic Green's function through the Ewald method,” IEEE Trans. Microw. Theory Tech. **56**(9), 2069–2075 (2008).

**15. **A. Vallecchi, M. Albani, and F. Capolino, “Collective electric and magnetic plasmonic resonances in spherical nanoclusters,” Opt. Express **19**(3), 2754–2772 (2011).

**16. **J. D. Jackson, *Classical Electrodynamics* (Wiley, 1998).

**17. **V. A. Markel, V. N. Pustovit, S. V. Karpov, A. V. Obuschenko, V. S. Gerasimov, and I. L. Isaev, “Electromagnetic density of states and absorption of radiation by aggregates of nanospheres with multipole interactions,” Phys. Rev. B **70**(5), 054202 (2004).

**18. **D. W. Mackowski, “Calculation of total cross sections of multiple-sphere clusters,” J. Opt. Soc. Am. A **11**(11), 2851–2861 (1994).

**19. **M. Abramowitz and I. A. Stegun, *Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables* (Dover Publications, 1965).

**20. **A. D. Rakic, A. B. Djurisic, J. M. Elazar, and M. L. Majewski, “Optical properties of metallic films for vertical-cavity optoelectronic devices,” Appl. Opt. **37**(22), 5271–5283 (1998).

**21. **S. Campione, M. Albani, and F. Capolino, “Complex modes and near-zero permittivity in 3D arrays of plasmonic nanoshells: loss compensation using gain [Invited],” Opt. Mater. Express **1**(6), 1077–1089 (2011).

**22. **A. Alù, A. Salandrino, and N. Engheta, “Negative effective permeability and left-handed materials at optical frequencies,” Opt. Express **14**(4), 1557–1567 (2006).

**23. **I. El-Kady, M. M. Sigalas, R. Biswas, K. M. Ho, and C. M. Soukoulis, “Metallic photonic crystals at optical wavelengths,” Phys. Rev. B **62**(23), 15299–15302 (2000).

**24. **P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” Phys. Rev. B **6**(12), 4370–4379 (1972).

**25. **M. G. Silveirinha, “Generalized Lorentz-Lorenz formulas for microstructured materials,” Phys. Rev. B **76**(24), 245117 (2007).

**26. **A. Alù, “First-principles homogenization theory for periodic metamaterials,” Phys. Rev. B **84**(7), 075153 (2011).

**27. **S. Steshenko, F. Capolino, P. Alitalo, and S. Tretyakov, “Effective model and investigation of the near-field enhancement and subwavelength imaging properties of multilayer arrays of plasmonic nanospheres,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **84**(1 Pt 2), 016607 (2011).

**28. **A. M. Nicolson and G. F. Ross, “Measurement of the intrinsic properties of materials by time-domain techniques,” IEEE Trans. Instrum. Meas. **19**(4), 377–382 (1970).

**29. **W. B. Weir, “Automatic measurement of complex dielectric constant and permeability at microwave frequencies,” Proc. IEEE **62**(1), 33–36 (1974).

**30. **A. H. Boughriet, C. Legrand, and A. Chapoton, “Noniterative stable transmission/reflection method for low-loss material complex permittivity determination,” IEEE Trans. Microw. Theory Tech. **45**(1), 52–57 (1997).

**31. **D. R. Smith, S. Schultz, P. Markoscaron, and C. M. Soukoulis, “Determination of effective permittivity and permeability of metamaterials from reflection and transmission coefficients,” Phys. Rev. B **65**(19), 195104 (2002).

**32. **C. R. Simovski, “On the extraction of local material parameters of metamaterials from experimental or simulated data,” in *Theory and Phenomena of Metamaterials*, F. Capolino, ed. (CRC Press, 2009), 11.1.

**33. **S. A. Ramakrishna and T. M. Grzegorczyk, *Physics and Applications of Negative Refractive Index Materials* (CRC Press and SPIE Press, 2009).

**34. **A. Sihvola, *Electromagnetic Mixing Formulas and Applications* (IEEE Publishing, 1999).

**35. **A. Sihvola, “Mixing rules,” in *Theory and Phenomena of Metamaterials*, F. Capolino, ed. (CRC Press, 2009), 9.1.

**36. **S. Campione, S. Steshenko, M. Albani, and F. Capolino, “Characterization of the optical modes in 3D-periodic arrays of metallic nanospheres,” in *URSI General Assembly and Scientific Symposium* (Istanbul, Turkey, 2011).

**37. **X. Chen, T. M. Grzegorczyk, B.-I. Wu, J. Pacheco Jr, and J. A. Kong, “Robust method to retrieve the constitutive effective parameters of metamaterials,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **70**(1 Pt 2), 016608 (2004).

**38. **P. P. Ewald, “Die Berechnung optischer und elektrostatischer Gitterpotentiale,” Ann. Phys. (Berlin) **369**(3), 253–287 (1921).

**39. **K. E. Jordan, G. R. Richter, and P. Sheng, “An efficient numerical evaluation of the Green's function for the Helmholtz operator on periodic structures,” J. Comput. Phys. **63**(1), 222–235 (1986).

**40. **A. Kustepeli and A. Q. Martin, “On the splitting parameter in the Ewald method,” IEEE Microw. Guid. Wave Lett. **10**(5), 168–170 (2000).