Enhanced THz extinction in arrays of resonant semiconductor particles

We demonstrate experimentally the enhanced THz extinction by periodic arrays of resonant semiconductor particles. This phenomenon is explained in terms of the radiative coupling of localized resonances with diffractive orders in the plane of the array (Rayleigh anomalies). The experimental results are described by numerical calculations using a coupled dipole model and by Finite-Difference in Time-Domain simulations. An optimum particle size for enhancing the extinction efficiency of the array is found. This optimum is determined by the frequency detuning between the localized resonances in the individual particles and the Rayleigh anomaly. The extinction calculations and measurements are also compared to nearfield simulations illustrating the optimum particle size for the enhancement of the near-field. © 2015 Optical Society of America OCIS codes: (230.1950) Diffraction gratings; (240.6680) Surface plasmons; (250.5403) Plasmonics; 300.6470 Spectroscopy, semiconductors; (300.6495) Spectroscopy, terahertz. References and links 1. J. Krenn, A. Dereux, J. Weeber, E. Bourillot, Y. Lacroute, J. Goudonnet, G. Schider, W. Gotschy, A. Leitner, F. Aussenegg, and C. Girard, “Squeezing the optical near-field zone by plasmon coupling of metallic nanoparticles,” Phys. Rev. Lett. 82, 2590–2593 (1999). 2. S. Zou and G.C. Schatz, “Silver nanoparticle array structures that produce giant enhancements in electromagnetic fields,” Chem. Phys. Lett. 403, 62–67 (2005). 3. J. Aizpurua, G. W. Bryant, L. J. Richter, F. J. Garcı́a De Abajo, B. K. Kelley, and T. Mallouk, “Optical properties of coupled metallic nanorods for field-enhanced spectroscopy,” Phys. Rev. B 71, 235420 (2005). 4. H. Fischer and O.J.F. Martin, “Engineering the optical response of plasmonic nanoantennas,” Opt. Express 16, 9144–9154 (2008). 5. S. Zou, N. Janel, and G. C. Schatz, “Silver nanoparticle array structures that produce remarkably narrow plasmon lineshapes,” J. Chem. Phys. 120, 10871–10875 (2004). 6. F. J. Garcı́a de Abajo, “Colloquium: Light scattering by particle and hole arrays,” Rev. Mod. Phys. 79, 1267 (2007). 7. F. J. Garcı́a de Abajo, R. Gómez-Medina, and J. J. Sáenz, “Full transmission through perfect-conductor subwavelength hole arrays,” Phys. Rev. E 72, 016608 (2005). 8. T. W. Ebbesen, H. Lezec, H. Ghaemi, T. Thio, and P. Wolff, “Extraordinary optical transmission through subwavelength hole arrays,” Nature (London) 391, 667–669 (1998). 9. L. Martı́n-Moreno, F. J. Garcı́a-Vidal, H. J. Lezec, K. M. Pellerin, T. Thio, J. B. Pendry, and T. W. Ebbesen, “Theory of extraordinary optical transmission through subwavelength hole arrays,” Phys. Rev. Lett. 86, 1114 (2001). 10. J. Bravo-Abad, L. Martı́n-Moreno, F. J. Garca-Vidal, E. Hendry, and J. Gómez Rivas, “Transmission of light through periodic arrays of square holes: from a metallic wire mesh to an array of tiny holes,” Phys. Rev. B 76, 241102(R) (2007). #240819 Received 12 May 2015; revised 9 Aug 2015; accepted 28 Aug 2015; published 10 Sep 2015 (C) 2015 OSA 21 Sep 2015 | Vol. 23, No. 19 | DOI:10.1364/OE.23.024440 | OPTICS EXPRESS 24440 11. E. M. Hicks, S. Zou, G. C. Schatz, K. G. Spears, R. P. Van Duyne, L. Gunnarsson, T. Rindzevicius, B. Kasemo, and M. Käll, “Controlling plasmon line shapes through diffractive coupling in linear arrays of cylindrical nanoparticles fabricated by electron beam lithography,” Nano Lett. 5, 1065–1070 (2005). 12. Y. Chu, E. Schonbrun, T. Yang, and K. B. Crozier, “Experimental observation of narrow surface plasmon resonances in gold nanoparticle arrays,” Appl. Phys. Lett. 93, 181108 (2008). 13. B. Auguié and W. Barnes, “Collective resonances in gold nanoparticle arrays,” Phys. Rev. Lett. 101, 1 (2008). 14. V. Kravets, F. Schedin, and A. Grigorenko, “Extremely narrow plasmon resonances based on diffraction coupling of localized plasmons in arrays of metallic nanoparticles,” Phys. Rev. Lett. 101, 087403 (2008). 15. G. Vecchi, V. Giannini, and J. Gómez Rivas, “Surface modes in plasmonic crystals induced by diffractive coupling of nanoantennas,” Phys. Rev. B 80, 201401 (2009). 16. S. R. K. Rodriguez, A. Abass, B. Maes, O. T. A. Janssen, G. Vecchi, and J. Gómez Rivas, “Coupling bright and dark plasmonic lattice resonances,” Phys. Rev. X 1, 021019 (2011). 17. W. Zhou and T. W. Odom, “Tunable subradiant lattice plasmons by out-of-plane dipolar interactions,” Nat. Nanotech. 6, 423–427 (2011). 18. T. V. Teperik and A. Degiron, “Design strategies to tailor the narrow plasmon-photonic resonances in arrays of metallic nanoparticles,” Phys. Rev. B 86, 245425 (2012). 19. G. Pellegrini, G. Mattei, and P. Mazzoldi, “Nanoantenna arrays for large-area emission enhancement,” J. Chem. Phys. C 115, 24662–24665 (2011). 20. A. B. Evlykhin, C. Reinhart, U. Zywietz, and B. N. Chichkov, “Collective resonances in metal nanoparticle arrays with dipole-quadrupole interactions,” Phys. Rev. B 85, 245411 (2012). 21. V. G. Kravets, F. Schedin, A. V. Kabashin, and A. N. Grigorenko, “Sensitivity of collective plasmon modes of gold nanoresonators to local environment,” Opt. Lett. 35, 956–598 (2010). 22. P. Offermans, M. C. Schaafsma, S. R. K. Rodriguez, Y. Zhang, M. Crego-Calama, S. Brongersma, and. J. Gómez Rivas, “Universal scaling of the figure of merit of plasmonic sensors,” ACS Nano 5, 5151–5157 (2011). 23. W. Zhou, M. Dridi, J. Y. Suh, C. H. Kim, D. T. Co, M. R. Wasielewski, G. C. Schatz, and T. W. Odom, “Lasing action in strongly coupled plasmonic nanocavity arrays,” Nat. Nanotech. 8, 506–511 (2013). 24. A. H. Schokker and A. F. Koenderink, “Lasing at the band edges of plasmonic lattices,” Phys. Rev. B 90, 155452 (2014). 25. A. Bitzer, J. Wallauer, H. Helm, H. Merbold, T. Feurer, and M. Walther, “Lattice modes mediate radiative coupling in metamaterial arrays,” Opt. Express 17, 22108–22113 (2009). 26. B. Ng, S. M. Hanham, V. Giannini, Z. C. Chen, M. Tang, Y. F. Liew, N. Klein, M. H. Hong, and S. A. Maier, “Lattice resonances in antenna arrays for liquid sensing in the terahertz regime,” Opt. Express 19, 14653–14661 (2011). 27. J. Wallauer, A. Bitzer, S. Waselikowski, and M. Walther, “Near-field signature of electromagnetic coupling in metamaterial arrays: a terahertz microscopy study,” Opt. Express 19, 17283–17292 (2011). 28. A. Berrier, R. Ulbricht, M. Bonn, and J. G. Rivas, “Ultrafast active control of localized surface plasmon resonances in silicon bowtie antennas,” Opt. Express 18, 23226–23235 (2010). 29. A. Berrier, P. Albella, M. A. Poyli, R. Ulbricht, M. Bonn, J. Aizpurua, and J. Gómez Rivas, “Detection of deep-subwavelength dielectric layers at terahertz frequencies using semiconductor plasmonic resonators,” Opt. Express 20, 5052–5060 (2012). 30. IOFFE, http://www.ioffe.ru/SVA/NSM/Semicond/Si/basic.html. 31. Y.-S. Lee, Principles of Terahertz Science and Technology (Springer US, 2009). 32. T. Jensen, L. Kelly, A. Lazarides, and G. C. Schatz, “Electrodynamics of noble metal nanoparticles and nanoparticle clusters,” J. Clust. Sci. 10, 295–317 (1999). 33. H. C. van de Hulst, Light Scattering by Small Particles (Dover Publications, Inc., 1981). 34. V. Giannini, A. Berrier, S. A. Maier, J. A. Sanchez-Gil, and J. Gómez Rivas, “Scattering efficiency and near field enhancement of active semiconductor plasmonic antennas at terahertz frequencies,” Opt. Lett. 18, 2798–2807 (2010). 35. M. C. Troparevsky, A. S. Sabau, A. R. Lupini, and Z. Zhang, “Transfer-matrix formalism for the calculation of optical response in multilayer systems: from coherent to incoherent interference,” Opt. Express 18, 24715–24721 (2010). 36. L. Duvillaret, F. Garet, and J.-L. Coutaz, “A reliable method for extraction of material parameters in terahertz time-domain spectroscopy,” IEEE J. Sel. Top. Quant. 2, 739–746 (1996). 37. N. Ashcroft and N. Mermin, Solid State Physics (Holt-Saunders, 1976). 38. C. Jacoboni, C. Canali, G. Ottaviani and A. Alberigi Quaranta, “A review of some charge transport properties of silicon,” Solid State Electron. 20, 77–89 (1977).


Introduction
Conducting particles can resonantly absorb and scatter electromagnetic waves by the so-called localized resonances (LRs) when their dimensions are comparable to the wavelength. The ab-sorption and scattering efficiencies depend on properties like the size, shape, material and orientation relative to the polarization of the wave. The local electric field can be increased in the vicinity of a particle due to such LRs [1][2][3][4]. When these scattering particles are separated by less than a few wavelengths they can couple radiatively. The radiative coupling can be enhanced in periodic arrays of particles through Rayleigh anomalies (RAs), which are diffraction orders in the plane of the array [5]. The resonances in periodic arrays of particles arising from the enhanced radiative coupling of LRs through diffraction are termed Surface Lattice Resonances (SLRs) [6]. SLRs have a significant influence on macroscopic optical properties of the array such as the optical extinction, and offer a new range of design parameters for tuning the spectral response of scattering structures. The extinction of light by periodic arrays is enhanced at the SLR condition as a consequence of the in-plane scattered wave by the particles in the array. This enhanced extinction has been associated through Babinet's principle to the extraordinary transmission through the complementary system, i.e. hole arrays [6,7], where surface waves on the continuous metal film (surface plasmon polaritons) are responsible of the enhanced optical transmission [8][9][10]. Aside from reshaping the extinction spectrum, the electromagnetic fields are enhanced over large volumes at specific frequencies due to SLRs. SLRs are currently under intense investigation at optical frequencies in arrays of plasmonic nanoparticles because their controllable properties that can be exploited in numerous applications such as for surface enhanced Raman scattering, bulk refractive index sensing, solid state lighting and lasing [11][12][13][14][15][16][17][18][19][20][21][22][23][24]. SLRs have been also described at THz frequencies using metallic structures [25][26][27].
In this manuscript we demonstrate experimentally, and verify theoretically the diffractive coupling at THz frequencies of LRs in semiconductor scatterers when they are ordered in a 2dimensional periodic array. In particular, we investigate arrays of silicon particles and show how the THz extinction of these arrays correlate with the frequency detuning between the LRs of the particles and the RA of the array. The extinction efficiency is maximum for a detuning close to zero and not for the largest particles as it could be incorrectly expected. We also correlate the far-field enhanced extinction with numerical simulations of the local field enhancement, illustrating an enhanced local field for a small detuning. Semiconductors are very promising materials for THz applications due to their versatile nature that can be controlled by the doping and free charge carrier concentration. Therefore, this demonstration of the diffractive coupling of LRs in Si structures defines new methods to efficiently modulate the THz extinction and near field.
The experiments and subsequential theoretical and numerical analysis presented in this manuscript consider samples where the Si structures are homogeneously surrounded by quartz, resulting in pronounced effects due to diffractive coupling and a better support of Rayleigh anomalies. For the design of optical elements based on diffractive coupling it might be desirable that the particles are accessible and placed at an interface of quartz and air. Under these conditions diffractive coupling occurs in both media, and not necessarily in phase. Although this weakens the reduced extinction at the RA condition and the enhanced extinction at the SLR frequency, a SLR will still be excited and the main mechanism of diffractive coupling will be similar to that of a sample symmetrically embedded in quartz.

Samples description
The investigated semiconductor structures consist out of 2D square arrays of 1.5 µm thick single-crystalline Si square tiles. The length d of the tile varies from d = 25 µm to d = 85 µm for the different samples, while the lattice constant Γ is kept constant at 100 µm. The fabrication has been done using standard optical lithography, as described in more detail in Refs. [28,29]. The layer containing the Si structures is bonded at both sides to an amorphous quartz substrate, which is transparent at THz frequencies. The bonding layer is benzocyclobutene (BCB) and has an estimated thickness of 3 µm. Optical microscope images of the samples are shown in Fig. 1 (a)-(e), for arrays with squares of 25 µm, 35 µm, 45 µm, 55 µm and 65 µm in length respectively. The last panel ( Fig. 1(f)) shows an unstructured Si reference layer. The images are taken after the lithography was done but before the last bonding step was carried out, in which the Si structures were bonded to the superstrate to form a symmetric sample. Intrinsic single-crystalline Si has a bandgap of 1.12 eV. At room temperature this energy exceeds the thermal energy k B T of the electrons. As a consequence most of the electrons will be confined to the valence band, and the material is transparent at THz frequencies. Under these circumstances Si behaves as a dielectric with a refractive index of n = 3.4 and is practically lossless [30]. A permanent metallic behavior has been introduced by implanting the Si with arsenic (As) atoms, with a concentration of free carriers in the order of N ≈ 10 19 cm −3 . The free charge carriers and the dimensions of the Si tiles are responsible for the LRs at THz frequencies in the structures.

Experimental results
The experiments have been carried out using a THz time-domain spectrometer (THz-TDS), which uses optical rectification in ZnTe of 100 fs pulses from a Ti:sapphire oscillator for the generation of broadband THz radiation. This setup uses electro-optic sampling of the fs pulses from the same oscillator in another ZnTe crystal for the detection of THz radiation [31]. THz-TDS measures THz electric field transients by varying the optical path length between the generation and detection optical pulses. In our measurements, we detect the THz transmission through the array of Si particles. The complex transmission t(ν) is obtained by Fourier transforming the THz transients, which provides the intensity spectra I(ν) = |t(ν)| 2 and the accumulated phase φ (ν) = arg(t(ν)). A sample without any silicon, i.e. only quartz and BCB, was used as a reference. For each array the transmitted intensity is normalized to this reference, and corrected for the filling fraction f of the Si in each unit cell ( f = d 2 /Γ 2 ) to determine the extinction efficiency The extinction efficiencies are shown in Fig. 2(a) as solid curves, for the arrays with particles ranging in length from 85 µm (black curve) to 25 µm (blue curve). The vertical dotted line indicates the RA condition at 1.6 THz, which corresponds to the wave-vector that is associated with the lattice constant in the samples.
Three regimes can be identified: For large particles, 85 µm (black curve)-65 µm (magenta curve), the response of the arrays show mainly a very broad resonance. For particles in the range between 55 µm (cyan curve)-35 µm (green curve) a characteristic surface lattice resonance appears close to the RA condition. The resonance becomes broad again for the smallest particles of 25 µm (blue curve).
The extinction of the arrays is a result of the induced electromagnetic field from the scatterers interfering destructively with the incident field at the position of the detector. The phase of this field that is reradiated by the scatterers, relative to the phase of the incident field, is shown in Fig. 2(b). This differential phase delay is calculated from the complex Fourier transforms of the time-domain transients as The phase spectra show a trend consistent with the extinction spectra. For the smallest and largest particles the curves are smooth, while for the particles with lengths in the range 55 µm (cyan curve)-35 µm (green curve) there is a pronounced local maximum in the phase. This maximum in phase appears at slightly higher frequencies than its corresponding maximum in extinction, and is consistently close to the RA condition.
In the next section we compare these measurements with numerical calculations using the coupled dipole model.

Coupled dipole model
The scattering of small particles in an ensemble can be described by a Coupled Dipole Model (CDM) calculation. In this model the scattering particles are approximated as dipoles, which can interact radiatively with each other. The properties of a scatterer are described by the polarizability α, which relates the polarization and the local field as p = αE loc . For each particle in the ensemble this local field is the sum of the incident field E inc and the field E ens which is scattered from all the other particles in the ensemble. For each particle i we can write This scattered field can be calculated from the product of the polarization vector of each particle and the dipole-dipole interaction tensor G(r) [6], via a summation over all particles j different than particle i E ens In differential form G(r) is defined as The above equations can be re-written as a linear set of self-consistent equations  where I is the identity matrix and matrix S contains all the dipolar interactions. For finite arrays these equations can be solved for P by calculating the inverse of the matrix A. The extinction efficiency is given by [32] with d the diameter of the scattering particles, k the wavenumber and n the number of particles in the lattice.
For the polarizability α we use the modified long wavelength approximation (MLWA), which holds for ellipsoidal particles of subwavelength dimension. This polarizability can be written as [32] where m = 1, 2, 3 represents the principal axes of the ellipsoid with corresponding lengths d m and α static m represents the static polarizability of the particle. The term 2 3 ik 3 α static m corresponds to the dynamic depolarization, and the term 2k 2 d m α static m to the radiative damping. The static polarizability for an ellipsoidal particle of permittivity ε p surrounded by a homogeneous medium of permittivity ε s is given by [33] α static with V the volume of the particle and L m is the shape factor taking into account the shape of the particle. This factor is given by For spherical particles the shape factor is 1/3 for the three axes. In the calculations presented in this manuscript, we consider square arrays with a lattice constant of 100 µm of 20 × 20 oblate spheroids representing disk-shaped particles. The oblate spheroids have two equal principal axes defining their diameter and a shorter one defining the height of the disk (i.e. 1.5 µm). The shape factors for the particles are for the long axes in the range 0.01 < L m < 0.03. The resulting resonance in the spectrum of these particles will be referred to as a localized resonance (LR).
As we see next, the approximation of the Si tiles as oblate spheroids particles in the CDM is sufficient to describe qualitatively the measurements.

CDM Extinction efficiency
The results of the CDM calculations are shown in Fig. 3. Figure 3(a) displays the extinction efficiencies of the arrays with solid curves, and the extinction of single particles with dashed curves. The spectra are vertically displaced for clarity as indicated by the horizontal dotted lines. The vertical dotted line indicates the Rayleigh anomaly condition at 1.5 THz. This frequency differs slightly from that of the experimental results. The difference can be attributed to the omission of the polymer bonding layer in the calculations. Results are shown for particles with lengths of 85 µm (black curve) to 25 µm (blue curve).
The resonant extinction of the single particles can be ascribed to LRs resulting from the coherent oscillation of the free charge carriers at the surface of the Si particles. The localized resonances red-shift as the length of the particle increases, while the shape of the resonance remains similar. The extinction efficiency is in first order independent of the particle length and is slightly larger than 2. The extinction spectra for the arrays are much richer. For the largest particles (black to magenta curves) the resonance is blue-shifted from the LR towards the RA condition, but remains broad. As the LR approaches the RA condition (cyan to green curves), the extinction is enhanced, and the resonance becomes narrower. The extinction efficiency is enhanced the most when the LR of the individual particles overlaps with the RA condition of the array, which corresponds to particles with dimensions between 35 µm and 45 µm (green and red curves in Fig. 3(a), respectively). These trends of enhanced extinction match the experimental results shown in Fig. 2(a). For the smallest particles (blue curve), when the LR is at much higher frequencies than the RA, the extinction efficiency at the RA condition is not so pronounced and the overall extinction of the array is not very different from that of the individual particles.

CDM phase
The phase of the average polarization is shown in Fig. 3(b), and is calculated as φ = 1 n ∑ i arg(p i ). Solid curves represent the array including interaction, while dashed curves represent single particles. The phase of the latter undergoes a smooth transition centered around the LR resonance frequency as expected for a damped Lorentz oscillator representing the oscillation of the free charges. The phase of the arrays follows a similar trend, but has an additional modulation slightly red-shifted from the RA condition. As in the case of the experimental results (Fig. 2) the maximum in phase appears at higher frequencies than the maximum in extinction, and remains close to the RA condition.
The mechanism behind this behavior in the phase is directly related to the diffractive coupling of the particles in the lattice, and can be understood in more detail when breaking down the total scattering into the different contributing components. The results are shown in Fig. 4 where each arrow represents a field component in the complex plane The modulus of the vectors represent the field amplitudes normalized to the incident field, which has an amplitude given by the radius of the circle; The angles in Fig. 4 represent their phase φ relative to the incident field. Figure 4(a) is the calculation for the single oblate spheroid with a length of 35 µm at 1.44 THz. This frequency corresponds to maximum extinction efficiency in the case of the array (see green curve in Fig. 3(a)). The incident field (gray arrow) drives the conducting electrons in the scatterer, which scatter a fraction of this field in the forward direction (red arrow). The magnitude and relative phase of this scattered field are determined by the complex value of the particle's polarizability. The coherent sum of the incident field and the scattered field results in a reduction and change in phase of the transmitted field, which is indicated with the blue ar-  row. Figure 4(b) represents a similar calculation as in Fig. 4(a) but at 1.5 THz. This frequency corresponds to the RA condition in the case of the array. As can be appreciated in Fig. 4(a), there are no significant changes in the fields compared to those calculated at 1.44 THz. This similar behavior is due to the broad resonance and similar value of the polarizability at both frequencies.
For the periodic lattice, an additional term is present to account for the diffractive coupling. This term corresponds to the contribution of the scattered field by the ensemble of particles to the local field E ens and it is indicated by the green arrows in Figs. 4(c) and 4(d), where the fields in Fig. 4(c) are calculated at the frequency of the SLR, i.e., maximum extinction by the array, and in Fig. 4(d) correspond to the frequency of the RA (see green curve in Fig. 3). For an array that extends to infinity the total scattered field at each particle is given by E ens = Sp, where S = ∑ i G(r i ) and p = α[E inc + Sp]. The magnitude and phase of the field resulting from this diffractive coupling are related to the periodicity of the lattice and therefore the detuning of the LR from the RA condition. The superposition of this field, that is diffracted in the plane of the lattice, with the incident field constitutes the local field (black arrows in Figs. 4(c) and 4(d)) at the position of the scatterers. The scattered field by the array depends on the local field at the scatterers and their polarizability and it is given by The calculated E sca is represented by the red arrows in Figs. 4(c) and 4(d) for the SLR and RA, respectively. In both cases, the phase difference between the local field and the scattered field is determined by the phase of the polarizability. The amplitude and phase of the diffracted field at the SLR condition (green arrow in Fig. 4(c)) change the local field such, that the field scattered by the particles (red arrow) is in anti-phase compared to the incident field (gray arrow), interfering destructively with each other. As a result, the extinction is maximized. The phase that is obtained from the extinction measurements is the phase of the scattered field. At the RA condition ( Fig. 4(d)) the wavelength equals the pitch (Γ = 2π/k). At this condition the diffracted field (green arrow Fig. 4(d)) at each scatterer has been delayed by an integer number of wavelengths. As a consequence, the diffracted field has approximately the same phase as the scattered field of the single scatterer of Fig. 4(b). The amplitude is such, that when the scattered field (red arrow in Fig. 4(d)) is added to the incident field, the resulting total field (blue arrow) lies still close to the unit circle. This can be understood as an amplitude in the far field which has not changed, i.e., low extinction, but only picked up a finite additional phase.
The system as described above is characterized by a LR which is at the high-energy side of the RA, resulting in a sharp SLR. Such behavior is not observed when the LR resides at the lowenergy side of the RA. In the latter case the phase of E sca for the isolated particle will exceed π, and as a consequence the phase and magnitude of E ens will change as well. This change of phase and magnitude prevents E sca in the periodic array to be equal in magnitude and opposite in phase to E inc , therefor preventing the formation of a sharp SLR [13, 18].

FDTD simulations
To get a better insight in the interpretation of the experiments, simulations have been performed using the finite-difference in the time-domain method (FDTD). The layout of the unit cell and dimensions of the particles are taken from the microscope images shown in Fig. 1. The permittivity of the Si is the same as for the CDM calculations. The surrounding medium is considered as quartz with minor losses ε s = 4 + 0.04i. For single particles, the simulation volume is surrounded by perfectly matched layers (PMLs) and the extinction cross sections simulated with FDTD have been normalized to the geometric cross sections of the Si particles, resulting in an extinction efficiency. The results for the periodic arrays are simulated using periodic boundary conditions (PBC) in the directions perpendicular to the propagation direction of the incident field and the extinction has been normalized to the filling fraction of the particles. The resulting extinction efficiencies for the arrays are shown in Fig. 5(a), while those for the single particles are shown in Fig. 5(b).
The extinction spectra for the single particles show a pronounced red-shift as the particle length increases. The magnitude of the extinction efficiency is nearly independent of the particle length. The finite extinction at high frequencies is attributed to higher order resonances. This contribution to the extinction from multipolar resonances is obviously not present in the dipole calculations of the previous section. For the periodic structures the results agree very well with the experimental data. For the large particles (85 µm-65 µm) the spectral response is broad, and there is no pronounced resonance around the RA condition. Also for the smaller particles the resonance is weak. For the arrays with particle lengths in the range 55 µm-35 µm a pronounced enhancement of the extinction is visible, slightly red-shifted from the Rayleigh anomaly. The FDTD simulations overestimate the magnitude and underestimate the width of the resonances compared to the measurements. This occurs mostly for the smaller particles and can be attributed to the finite bandwidth and cross-section of the experimental THz beam, while the FDTD simulations have been carried out with plane-wave illumination.

Discussion
A resonance is characterized by its central frequency, magnitude and width. These parameters have been extracted from the experimental data (Fig. 2) and the FDTD simulations (Fig. 5) and are shown in Fig. 6. The red triangles in this figure correspond to the experimental data of the SLRs in the arrays, while the blue circles are the simulated SLRs and the crosses are the simulated LRs. The SLR frequency, defined as the frequency of maximum extinction and shown in  This small change of the resonance frequency compared to that of the single particles (represented in the figure inset) when the size of the particles is varied, indicates that the diffractive coupling strongly modifies the collective response of the localized resonances. The extinction efficiency of the arrays (Fig. 6(b)) has a maximum value of 4 at a particle length of 35 µm, decreasing rapidly for larger and smaller particle lengths. The extinction efficiency of the single particles does not change significantly with the particle length. It is not practical to assign a full width at half maximum to the resonances due to the asymmetric line shape of the extinction spectra, with an extinction that remains significant at high frequencies in both simulation and experiment. The extinction does however drop to zero at low frequencies, and therefore the lower half width at half maximum (HWHM) is used to quantify the widths of the resonances. As can be seen in Fig. 6(c), for both the FDTD and the experimental results the width of the SLRs increases as the particle length is increased. The width of the LRs of the single particles has the opposite behavior.
In order to describe these dependencies of the SLRs with the particle length and their physical origin, we introduce a dimensionless parameter, δ , that defines the frequency detuning between the RA and the LRs For the FDTD simulations this detuning has been calculated for each particle size using the data presented in Fig. 6(a). The calculated frequency detunings are displayed for clarity in the upper horizontal axes of Fig. 6. As can be appreciated in Fig. 6(b) the maximum extinction efficiency for the arrays occurs for a detuning close to zero, i.e. δ 0. Also when this condition is satisfied, the HWHM of the SLRs is minimum. These results illustrate that the SLRs leading to the enhanced extinction can be regarded as the result of the coupling of LRs in the individual particles to the RAs or diffracted orders in the plane of the array. This coupling is optimum for small detuning leading to the maximum in the SLRs extinction. SLRs are also characterized by narrow line widths, which are mainly the result of the reduced radiative damping by destructive interference of the scattered fields [5].
The same analysis has been done for the SLRs calculated with the coupled dipole model. This is shown in Fig. 7, where the results of a 20 × 20 lattice (solid curves) are compared with those of single scatterers (dashed curves). Even though the CDM approximated the rectangular tiles as the dipolar response of oblate spheroids, and ignored the presence of the bonding layer, there is a qualitative good agreement between the experimental and numerical results. The CDM predicts in Fig. 7(a) a frequency shift of the SLR which scales with the particle length (solid curve). This shift is much more pronounced for the resonance frequency of the single particles (dashed curve). For small particles the LR is at high frequencies and the RA condition asymptotically keeps the SLR at 1.5 THz. In the experiments and FDTD the larger particles are no longer in the limit of truly dipolar behavior. The CDM neglects any contribution from multipolar modes in these larger particles that couple through diffraction at the RA condition at 1.5 THz.
The dependence of the extinction of the SLR on the particle length shown in Fig. 7(b) with the solid curve agrees very well with the experiments, reproducing the maximum around the particle length of 35 µm at which δ 0. The extinction efficiency of the single particles is practically independent of the particle length (dashed curve). The small increase of the extinction efficiency for smaller particles can be attributed to the material dispersion.
The HWHM of the LR (dashed curve in Fig. 7(c)) scales inversely with the particle length. The solid curve in Fig. 7(c), describing the array, displays a completely different behavior. A reduction in the length of the particle reduces the width of the resonance, which is a direct result of the SLR approaching the RA condition, reaching a minimum HWHM when δ 0. This is in agreement with the experimental results.

Near fields
The enhanced extinction efficiency in the particle arrays is closely related to a redistribution of energy in the near field in the proximity of the particles. This is illustrated in this section with FDTD simulations. We compare single particles to periodic lattices of particles at their respective resonance frequencies. These frequencies are indicated in Fig. 5 with the colored circles for each of the particle lengths under consideration. Cross-sections of the near fields are shown in Fig. 8 for particles of 25 µm, 35 µm, 45 µm and 55 µm in length. The cuts are taken through the center of each particle in the plane defined by the wavevector and the polarization vector of the incident field. Panels 8(a)-8(d) in Fig. 8 show the near-field intensity enhancements in one unit cell of the periodic lattices for increasing particle length. The horizontal dark bars at z = 0 µm represent the particles. We see that the fields are enhanced up to two orders of magnitude at the edges of the particles, as a result of the accumulation of the charges that are driven by the incident THz electric field. The fields however are enhanced not only in the proximity of the particle, but also extend in the region below the plane of the array. This effect is the largest in panel 8(b) where the local field intensity of the particle of 35 µm in length is shown. This is also in agreement with the extinction measurements and simulations discussed in the preceding sections, where for the 35 µm particles the largest extinction efficiency was observed. The regions of reduced intensity enhancement in the upper halves of the figures are a result of destructive interference between the incident field and the fraction of the THz radiation that is back scattered from the particles. For the single particles at their respective LR frequencies (panels 8(e)-8(h) of Fig. 8) the field is enhanced mostly at the edges of the particle, but this enhancement is not as pronounced as for the particles in the array. Figure 9 shows the averaged intensity enhancement (AIE), defined as the simulated intensity enhancement at the resonance frequencies integrated and averaged over a full unit-cell (dimensions 100 µm × 100 µm × 100 µm), as a function of the length of the particles. Figure 9(a) shows the AIE for the arrays of particles at each SLR frequency. The AIE of the arrays reproduces the behavior of the far-field extinction shown in Fig. 6(b), with a maximum AIE exceeding 2.5 for the array of particles with a length of 35 µm. The simulation for the single particles, for which the results are shown in Fig. 9(b), does not show an optimum particle size. There is only a slight increase of the AIE as the particle length increases, which can be explained by the larger volume over which the field in enhanced. For comparison, the FDTD simulations dis-cussed above are repeated for particles made out of a perfect electric conductor (PEC). Results for the AIE at the SLR and LR frequencies are shown in Figs. 9(c) and 9(d), respectively. Due to the suppression of Ohmic losses in PECs and the larger scattering efficiency compared to the Si particles, the AIE is larger for the PEC particles. Using semiconductors with a larger carrier mobility than Si, like GaAs or InSb, should provide a larger AIE performing even better than PECs due to the plamonic behavior of these semicondutors in the THz frequency range [34].

Conclusions
We have experimentally investigated surface lattice resonances in periodic arrays of intrinsically doped Si particles at THz frequencies. The mechanism of the coupling of localized resonances in the particles to diffractive orders in the plane of the arrays (Rayleigh anomalies) has been studied as a function of the size of the particles. A maximum extinction efficiency is found for a small frequency detuning between the localized reasonances and the Rayleigh anomaly. The experimental results are reproduced by numerical calculations using the coupled dipole model and electrodynamic simulations using the Finite Difference in Time Domain method. Simulations also show that the optimum extinction efficiency is correlated with a maximum near field intensity enhancement.