Lattice resonances in dielectric metasurfaces

Our capacity of structuring matter at the nanoscale has opened a myriad of possibilities for applications at optical frequencies. The resonant response of nanostructures at these frequencies leads to local field enhancements, increased light-matter interaction, enhanced nonlinearities and high sensitivity to small perturbations. These phenomena can be exploited for optical detection, light emission, light harvesting, and sensing.¹ A significant scientific effort has been made over the past few decades in controlling and enhancing the optical response of metallic nanoparticles.² The nanoparticles support coherent oscillations of free electrons, known as localized surface plasmon polaritons, which can locally enhance electromagnetic fields in tiny volumes.³ Intrinsic to metals are Ohmic losses due to the complex conductivity of real metals. In order to suppress these losses, more recent research has focused on low loss and high refractive index dielectric nanostructures supporting Mie resonances at optical frequencies.^4–7 These resonances, arising from displacement currents in the dielectric nanostructures, also lead to high local field enhancements and low optical losses.⁸ Moreover, compared with plasmonics, dielectric nanostructures can exhibit both electric and magnetic resonances.^9–12 

Real applications typically require field enhancements over larger volumes than those of single resonant nanoparticles. Metasurfaces and periodic lattices have emerged as systems formed by resonant nanoparticles with unusual characteristics emerging from their collective response.¹³ Arrays of metallic nanoparticles supporting collective resonances have received significant attention due to their remarkably narrow line widths (high Q-factors),^14–17 low radiation losses, and high field enhancements over large volumes.¹⁸ These collective resonances originate from the radiative coupling of localized resonances in the individual nanoparticles. This radiative coupling can be enhanced by diffraction orders in the plane of the arrays, known as Rayleigh anomalies (RAs), or by optical guided modes in thin dielectric layers. Collective resonances emerging from RAs are known as surface lattice resonances (SLRs),^19–22 while those emerging from guided modes are known as waveguide polaritons or quasiguided modes (QGMs),^23–25 although they are usually not distinguished from SLRs in the literature. The high Q-factor of collective resonances has been the reason for proposing them in applications, such as sensing,^26–28 spectroscopy,²⁹ surface-enhanced Raman spectroscopy,^30,31 solid-state lighting,^32,33 etendue reduction,³⁴ and lasing.^35–37 Similar to single nanoparticles, collective resonances in arrays of dielectric and semiconductor resonant nanoparticles have been recently proposed as an alternative to resonances in metallic arrays.^5,38–40

In this paper, we present a detailed theoretical investigation of collective resonances in arrays of resonant silicon (Si) nanoparticles using finite-difference time-domain (FDTD) simulations. Their narrow line width allows us to clearly resolve the collective resonance emerging from the diffraction enhanced radiative coupling of electric and magnetic resonances, which overlap in single nanoparticles. We illustrate the notorious differences between SLRs and quasiguided modes in Si nanoparticle arrays. Near-field simulations illustrate the strong electromagnetic field confinements that result in the collective resonances. We specially consider the field enhancement as a function of the distance to the nanoparticles. For an illumination along the normal direction to the surface, this field enhancement represents the eventual enhanced absorption that a material will experience at a particular height. It also provides the fractional local density of states or the efficiency at which the material will emit in the normal direction.

We first describe the optical resonances supported by individual Si nanoparticles. The optical constants of Si are taken from Palik.⁴¹ The nanoparticles have a cylindrical shape and are embedded in fused silica with constant refractive index $n=1.46$⁠. We select the height of the particle so that both electric and magnetic multipoles are supported, specifically the lowest order modes, i.e., the electric dipole (ED) and magnetic dipole (MD). The ED arises from the displacement current along the polarization of the incident electric field. The MD requires field retardation in the particle along the propagation direction of the incident wave, such that a displacement current loop is formed, leading to a magnetic field.^7,42 If the nanoparticle is too shallow, the field retardation is not sufficient to build a current loop. From calculations not shown here, we have determined that a height of 100 nm satisfies this requirement and, therefore, we fix the nanoparticle height to this value.

We have calculated the scattering and absorption efficiencies, i.e., the ratio of the scattering (absorption) cross section to the geometrical cross section of the particle, for cylinder diameters in the range 70–160 nm [Figs. 1(a) and 1(b), respectively]. A total-field scattered-field (TFSF) plane wave illumination (see Sec. IV) is used with a broadband spectrum, which is incident from below along the height of the particle. The scattering efficiency in Fig. 1(a) exhibits a strong peak that redshifts when the diameter is increased. This peak corresponds to the MD resonance. The peak is not symmetric and presents a shoulder at shorter wavelengths, which is attributed to the ED resonance. The spectral overlap between MD and ED resonances is a consequence of the refractive index of the surrounding medium that increases the radiative losses and broadens the resonances.⁷ For diameters larger than 120 nm, other peaks appear at shorter wavelengths, which are attributed to higher order modes.

Compared to the scattering, the absorption efficiency [Fig. 1(b)] in the silicon nanoparticles decreases as the diameter is increased. The high absorption at shorter wavelengths and small diameters is not caused by interband transitions on Si and is attributed to the interaction of the particle with the electromagnetic fields, as described by the Mie theory.⁴³ The real and imaginary components of the electric and magnetic multipole moments have different dependence on the particle size. The absorption at short wavelengths for small diameters is attributed to the magnetic and electric dipole resonances, while for larger diameters, it is due to higher order modes. Aside from these absorption peaks, we note that in the range of 600–700 nm, the absorption losses are particularly low. This makes silicon nanoparticles suitable structures for the design of lattice resonances in this spectral range.

The ED and MD resonances can be identified calculating the spatial distribution of the electric field intensity enhancement relative to the incident field intensity. We take the diameter value of $d=110$ nm, indicated in Fig. 1(a) by the white dashed line, and calculate the field distribution in the $yz$-plane for the two wavelengths corresponding to the MD and ED resonances, indicated by the black dots in Fig. 1(a). The color scale for both Figs. 1(c) and 1(d) represents the field intensity and the white arrows represent the real part of the vectorial electric field projected in the $yz$-plane. For the ED resonance, the electric field in Fig. 1(c) shows no circulation and is to a large extent oriented parallel to the incident electric field. The near-field outside the particle is distorted with respect the typical near-field created by a pure electric dipole in a disk and is a consequence of the retardation effect inside the particle. In Fig. 1(d), we can see that the electric field circulates inside the silicon particle, and thus corresponds with the magnetic dipole resonance.

SLRs arise from the enhanced radiative coupling of particle resonances by the in-plane diffraction orders of the lattice. For simplicity, we call these diffraction orders Rayleigh anomalies (RAs) in reference to the condition at which a diffraction order becomes evanescent at the interface separating two media. To investigate these resonances, we arrange the particles in an infinite periodic array with hexagonal geometry. We simulate this array embedded in a symmetric environment, i.e., with a substrate and superstrate with the same refractive index (⁠$n=1.46$⁠). This condition maximizes the radiative in-plane coupling between particles due to phase matching of the in-plane diffraction orders in the lower and upper media. A schematic representation of the array is shown in Fig. 2(a).

To tune the frequency of the SLRs to a particular spectral range, we have to consider the coupling strength between the optical resonances of the individual silicon cylinders to the RA. This coupling determines the redshift of the SLR with respect to the RA and its line width. The position of the RA can be calculated using the grating equation

where $k→∥d$ is the wave vector for the in-plane diffracted orders, $k→∥i$ is the in-plane projection of the incident wave vector $k→i$⁠, and $G→$ is the reciprocal lattice vector of the hexagonal array. If we choose a periodicity $a$ = 485 nm, we find that the first order RA happens at 613 nm for normal incidence (⁠$k∥i=0$⁠). The optical resonances of single Si cylinders with diameters between 80 nm and 110 nm are sufficiently detuned from the RA [see Fig. 1(a)]. For these diameters, we use FDTD simulations (see Sec. IV) to calculate the optical extinction of the array, defined as $1−T$⁠, where $T$ is the zero order transmittance for a plane wave incident normal to the array plane from below. The results are displayed in Fig. 2(b) for the wavelength range 400–700 nm. A series of peaks around 600–650 nm can be seen, in addition to the optical resonances of the individual particles at shorter wavelengths. These spectral features, not present in the spectrum of individual silicon particles, correspond to the SLRs and are shown in more detail in Fig. 2(c). We can see several peaks that redshift and broaden as the diameter of the nanoparticles increases. This behavior corresponds to a stronger coupling between the single particle resonances and the RA as their detuning is decreased—note that changing the particle diameter shifts the localized resonances of the individual nanoparticles, but does not change the RA, which depends on the lattice constant. The most relevant feature in the extinction spectra is the emergence of a second SLR for $d>90$ nm, which increases in intensity and redshifts with the diameter in a similar fashion as the first SLR at longer wavelengths. The formation of two collective resonances in the array has its origin in the electric and magnetic dipole resonances of the single nanoparticles. This result was already demonstrated by Evlyukhin et al. using a coupled-dipole model for both the electric and magnetic dipoles.⁵ The total induced electric (magnetic) dipole moment of an individual nanoparticle in the array is the sum of the induced electric (magnetic) dipole moment by the incident field in this particle and by the rest of the particles in the array. At a particular wavelength, this leads to the coherent and collective magnetic response of the array. Along with the formation of two SLRs, we note their highly symmetric line profile. Typically, SLRs exhibit an asymmetric line shape—Fano resonance⁴⁴—which results from the interference between the broad resonance in the individual nanoparticles and the narrow resonance given by the RA. In our system, this interference is small due to the large detuning between these two resonances.

To identify the nature of the SLRs in Fig. 2(c), we investigate the spatial distribution of the electric and magnetic field intensity enhancements relative to the incident field intensity. We focus on the array of particles with diameter $d$ = 110 nm, as this array exhibits both SLRs. The electric and magnetic field profiles on the $yz$-plane are shown in Figs. 3(a) and 3(c) for the SLR at 617 nm, and in Figs. 3(b) and 3(d) for the SLR at 627 nm. The color scale represents the magnetic and electric field intensity enhancements, in logarithmic scale, while the black arrows are the real components of the vectorial electric and magnetic fields projected in the $yz$-plane. At $λ$ = 617 nm, the circulation of the electric fields results in strong magnetic fields inside the silicon particle [Fig. 3(c)]. These fields indicate that the SLR at 617 nm arises from the collective coupling of magnetic dipoles in the nanoparticles, and we label this resonance as MD-SLR. In Fig. 3(b), the electric field is aligned with the incident field and the intensity shows a dipolar field pattern. Therefore, the SLR formed at 627 nm corresponds to the collective coupling of electric dipoles in the nanoparticles, and we label it as ED-SLR. Note that induced electric and magnetic dipole moments inside the Si particle for the ED-SLR and MD-SLR, respectively, are oriented orthogonal to each other along the direction of the incident fields. This result is also observed in the couple-dipole model for infinite arrays and implies an orthogonality relation between the ED-SLR and the MD-SLR at normal incidence.⁵ Compared to the local fields of the ED and MD resonances in the single particle illustrated in Figs. 1(c) and 1(d), the electric and magnetic fields for both ED-SLR and MD-SLR show a strong intensity enhancement, which is confined near the array plane. This in-plane confinement is attributed to the diffracted orders grazing to the sample which are responsible for the enhanced radiative coupling of the nanoparticles.

The simulated arrays present a MD-SLR blue-shifted with respect to the ED-SLR. This blue-shift indicates that the coupling strength between the RA and the magnetic dipole resonance is smaller than in the case of the electric dipole resonance. We can also observe this different coupling in the spatial distribution of the magnetic field intensity in Fig. 3(c): Although the intensity is much higher than in Fig. 3(d), it is mainly located inside the silicon nanoparticle. However, the magnetic and electric fields of the ED-SLR are extended outside the nanoparticles with larger intensities. It is important to stress that the coupling strength can be tuned with both the periodicity (i.e., the RA) and the geometry of the nanoparticles (i.e., the Mie resonances). In rectangular arrays, the periodicities in the x- and y-directions allow tuning separately the coupling for ED-SLR and MD-SLR, as their orthogonality relation means that they couple along perpendicular directions. Therefore, it is possible to achieve a crossover between the ED-SLR and MD-SLR and invert the spectral position of the ED-SLR respect to the MD-SLR.^12,40,45

To demonstrate the potential of SLRs in dielectric arrays to interact with thin layers of optically active materials, we have calculated the total intensity enhancement (IE), also known as excitation enhancement.⁴⁶ The IE is defined as the integral of the field intensity enhancement over the volume or area occupied by the material, normalized to a reference without the nanoparticle array. Specifically, we aim to illustrate the dependence of the IE with the distance to the particle array plane, which corresponds to a situation in which the SLRs interact with planar and thin-layer materials.⁴⁷ Therefore, we discretize the IE as a surface integral in the $xy$-plane that is evaluated as a function of the position $z$ and wavelength,

where $E(x,y,z,λ)$ is the electric field at wavelength $λ$ and position $z$ when the nanoparticle array is present, and $Eref(x,y,z,λ)$ corresponds to the electric field in the absence of the array. The integrals are evaluated in a unit cell of the array with area $A$⁠. The same method can be used to calculate the IE of the magnetic field. However, from the simulations shown in Fig. 3, we concluded that most of the MD-SLR enhancement is inside the Si nanoparticles, which makes it inaccessible to magnetic field active materials. Consequently, we focus on the electric field. For optically active materials, such as molecular dyes and quantum dots, with transition dipole moments randomly oriented, the modulus of $|E(x,y,z,λ)|$ includes the three field components, $Ex$⁠, $Ey$⁠, and $Ez$⁠. We have also considered the coupling of the SLRs to two-dimensional materials (e.g., 2D transition metal dichalcogenides or quantum wells), for which the optical transitions have dipole moments confined in the plane. In this case, only the projected components of $E(x,y,z,λ)$ in the $xy$-plane are evaluated in the integral. Further, the integral excludes the area of the unit cell occupied by the silicon particles.

Simulations of the IE for the array with $a=485$ nm, $h=100$ nm, and $d=110$ nm are shown for the MD-SLR and the ED-SLR in Fig. 4. The red curves are calculations including the three components of the electric field in the IE, while the blue curves include only the in-plane components. The gray area indicates the points where the integral excludes the electric fields inside the nanoparticles. The nanoparticle upper boundary leads to the small peaks observed in Figs. 4(a) and 4(b) at $z=100$ nm. We can see for the MD-SLR in Fig. 4(a) that the IE is mainly for the $z$-component of the electric field. This response is attributed to the diffraction enhanced coupling between magnetic dipole moments oriented along the $x$-axis. Since the electric field is given by the curl of the magnetic field, this creates a field strongly oriented in the direction perpendicular to the $xy$-plane. In contrast, the ED-SLR in Fig. 4(b) has the IE mainly in the $xy$-plane, which corresponds with the diffraction enhanced coupling between electric dipole moments in the $xy$-plane.

By virtue of this particular field orientation of the MD-SLR and the ED-SLR, we note that materials with different orientations in transition dipole moments will couple preferentially to one of the two SLRs. In addition, we can see from these calculations, that such an IE would enhance the absorption of the material if it overlaps spectrally and spatially with the SLR. Further, the IE also provides information about the directivity of the emission of emitters coupled to SLRs: Using Lorentz reciprocity theorem, the emission probability of a dipole in certain direction, i.e., the fractional radiative local density of states (FLDOS), is proportional to the local field enhancement produced by an illuminating beam incident in that direction.^39,47,48 Therefore, layers of emitting materials positioned at the height of the maximum field enhancement will radiate in the direction normal to the surface.

Unlike SLRs, which are favored when the array of nanoparticles is embedded in a homogeneous dielectric, quasiguided modes (QGMs) need a waveguide structure with a higher refractive index than the surrounding medium. Without the periodic array of scatterers, guided modes lie outside the lightcone, decoupled from the radiation modes as a consequence of the momentum mismatch. An array of particles acts as a grating coupler, providing to radiative modes the required momentum to couple to guided modes.^23,25 This makes possible the excitation of guided modes, but by the same argument, such modes are now leaky and couple out to free space propagating modes. Hence, the term “quasiguided”.

Although quasiguided modes have been observed in plasmonic particle arrays, their excitation is possible using dielectric scatterers as the silicon cylinders studied here.^25,38 The proposed structure is shown schematically in Fig. 5(a) and consists of a hexagonal array with $a$ = 485 nm, where the Si nanoparticles, with $h$ = 100 nm, are placed on top of a fused silica substrate (⁠$n$ = 1.46). On top of the substrate and embedding the particles, there is a slab waveguide with $n$ = 1.59 and a thickness of 700 nm, and on top of this waveguide, there is an infinite layer of air (⁠$n$ = 1). The value of the waveguide thickness is considered so that there are quasiguided modes supported in the wavelength range of 600–700 nm. As with the RAs, the coupling strength between guided modes and the single particle resonances is determined by the detuning between them. Consequently, and following the same argument than in the Sec. II B, we vary the diameter of the Si particles in the range 80–110 nm. The extinction results from the FDTD simulations are illustrated in Figs. 5(b) and 5(c) for a plane wave like illumination incident from the substrate along the normal direction (see Sec. IV). The different spectra exhibit a rich diversity of resonances. If we focus on Fig. 5(b) in the range of 400–600 nm, where the Mie resonances of individual silicon particles are excited, we now see several peaks that shift and increase as the diameter increases. These resonances are attributed to high order quasiguided modes that couple to the Mie resonances. Considering that radiative and absorption losses increase for the silicon particles at shorter wavelengths [Figs. 1(a) and 1(b)], we focus on the range of 600–700 nm, where we find several collective resonances. We can see in Fig. 5(c) two main resonances with large extinction next to weaker resonances with very narrow line widths. These four peaks redshift and broaden as the diameter of the particles is increased, with a pronounced increase in extinction for the weaker resonances. Further, the two largest peaks have an asymmetric profile not observed when the environment has a homogeneous refractive index [see Fig. 2(c)]. The reason for this asymmetry lies in the interference between the collective mode and light that is reflected back and forth inside the waveguide.

The resonances in a grating-waveguide system can be predicted analytically by combining Eq. (1) with the phase condition for constructive interference in the waveguide. However, this model assumes the grating elements to be small and with low scattering efficiency, which is not the case for our system. The presence of the nanoparticles changes drastically how light propagates inside the waveguide. In addition, the different refractive index of the particle also changes the optical path of a light ray compared to the waveguide without the array. Therefore, to understand the nature of the four resonances observed in the spectral range of 600–700 nm in Fig. 5, we simulate the spatial distribution of the electric and magnetic fields. We focus again on the array of Si nanoparticles with $d$ = 110 nm. In Fig. 6, we show the electric and magnetic field intensity enhancements (relative to the incident field) for the four resonances in the $yz$-plane of the simulation unit cell, and for the same illumination conditions as in Fig. 5. The color scale represents the intensity enhancement of the electric and magnetic fields in a logarithmic scale. The black arrows represent the real part of the vectorial fields projected in the $yz$-plane. The shape of the Si cylinder is indicated with the white rectangle and the white horizontal lines indicate the boundaries of the waveguide. We can see that both the electric [Figs. 6(a)–6(d)] and magnetic [Figs. 6(e)–6(h)] fields are strongly confined inside the waveguide. The narrow resonances at $λ$ = 617 nm [Figs. 6(a) and 6(e)] and $λ$ = 654 nm [Figs. 6(c) and 6(g)] show circulation of the electric field, a signature of a magnetic mode, which is confirmed by the strong magnetic fields observed in Figs. 6(e) and 6(g). Compared with the resonance at $λ$ = 654 nm, we can observe for $λ$ = 617 nm in Fig. 6(a) that there are two regions around which the electric field has circulation, one in the bottom of the waveguide, along the array plane and the other in the upper part of the waveguide. Moreover, the circulation in the upper part is opposite to the circulation in the bottom. As a result, the magnetic field created is antisymmetric in the waveguide along the $z$-axis. Given that the magnetic field is strongly oriented in the $x$-axis, we can identify this resonance at $λ$ = 617 nm with a first order quasiguided TM mode, while the resonance at $λ$ = 654 nm corresponds to the fundamental quasiguided TM mode.⁴⁹ We label them as TM0 and TM1.

We turn now our focus to the resonances at $λ$ = 628 nm [Figs. 6(b) and 6(f)] and $λ$ = 667 nm [Figs. 6(d) and 6(h)]. The electric field shows a strong enhancement for $λ$ = 628 nm [Fig. 6(b)] in the waveguide and around the nanoparticle. From the $yz$-plane projection of the electric field, represented by the black arrows, we can observe that the electric field is mainly oriented along the $y$-axis inside the waveguide and that has an antisymmetric behavior along the $z$-axis. This resonance resembles a first order TE quasiguided mode that we label as TE1. The resonance at $λ=667$ nm should correspond then with the fundamental quasiguided TE0 mode. However, the electric field is mainly confined around the nanoparticles, resembling instead an ED-SLR. This result is a consequence of the perturbation introduced in the spatial distribution of the fields by the scattering of the nanoparticles. Silicon particles scatter light in the waveguide in many directions, affecting the phase relationship of the fields and their interference, which has to include now the scattered fields as well. For particles with a large scattering cross section, as the Si particles with diameter $d=110$ nm, the perturbation is high, and the condition for in-phase accumulation of fields in the waveguide is not satisfied anymore. Radiative coupling between the particles, which depends on the phase relationship between in-plane scattered fields, is not affected and the result is that an ED-SLR-like mode is excited, that we label as TE0-SLR. The transition between TE0 and TE0-SLR depends on the scattering cross section of the particle. Decreasing the scattering cross section will reduce the perturbation and bring back the constructive interference in the waveguide. Such a transition is possible since both TE0 and ED-SLR have the same symmetry and are transverse-electric modes. The effect of the particle scattering on the spatial distribution of fields in QGMs has not been addressed in detail in the literature.^25,47

Next, we investigate the field enhancement (IE) associated to the quasiguided modes shown in Fig. 6 for the particle with diameter $d=110$ nm. We have calculated the IE as a function of the vertical distance to the array at the resonance wavelengths using Eq. (2) and excluding the field inside the nanoparticles. Both the array and waveguide are considered as the resonant structure, so the reference consists of a homogeneous medium with an identical refractive index than the substrate. To compare the results with the IE calculated for the homogeneous array, we focus on the electric field. We note, however, that the magnetic field enhancement is not negligible inside the waveguide for the TM0 and TM1 modes [see Figs. 6(e) and 6(g)]. The results of the IE are shown in Fig. 7. The red and blue curves represent the IE considering all the electric field components and only the in-plane components, respectively. The gray (yellow) area indicates the $z$-positions, where Si particles (waveguide) is present. The TM0 and TM1 modes in Figs. 7(a) and 7(c) have a significant enhancement of the $Ez$ component, as expected for TM modes where the magnetic component is oriented along the $x$-axis. The IE for the TM0 is remarkably high and can be attributed to the high Q-factor of this resonance. In contrast, in the TE0-SLR and TE1 modes of Figs. 7(b) and 7(d), the IE shows a larger contribution of in-plane electric fields, as expected from the electric dipole excited in the nanoparticles. Compared with the IE of the MD-SLR in Fig. 4(a), we note that the IE in Fig. 7(a) reaches a maximum inside the waveguide, in agreement with our description of the resonance as a TM0 mode. First order modes present two maxima, which can be seen for the TM1 and the TE1 in Figs. 7(c) and 7(d). Only the TE0-SLR mode in Fig. 7(b) shows a different IE profile as a function of $z$⁠. Due to the perturbation introduced by the particles, the fields are only enhanced in the array plane, which leads to an ED-SLR-like IE [see Fig. 4(b)].

The application of QGMs for light-matter interaction calls for more precaution compared with the SLRs. In addition to the different orientation of the field components, we have to consider the spatial dependence of the IE. The fields for the TE and TM have maxima and minima inside the waveguide. Only materials that are spatially overlapping the regions with high IE will experience an increase in excitation rate and, by Lorentz reciprocity, strong emission in the normal direction. Materials positioned at $z=200$ nm will be effectively quenched in the TE1 mode. This effect was measured for metallic array-waveguides, where it was shown that the photoluminescence enhancement (PLE) changed for emitters situated in different positions within the waveguide.⁴⁷ 

We have investigated numerically the diffraction enhanced radiative coupling of Mie resonances in arrays of dielectric nanoparticles. In particular, we have investigated hexagonal arrays of Si nanopillars with large scattering cross sections. We differentiate two systems depending on the media surrounding the nanoparticles: On the one hand, the array supports surface lattice resonances arising from the enhanced radiative coupling of Si nanoparticles if the array is homogeneously surrounded by a dielectric. On the other hand, when the array is embedded in a waveguide, it supports quasiguided modes. Common to these resonances is the significant line narrowing due to the reduction of radiation losses. We also determine the spatial dependence of the electric field enhancement in these systems, finding remarkable differences between SLRs and QGMs. A correct and detailed description of these field enhancements is of utmost importance for the design and optimization of applications relying on high local fields over large areas, such as light emission, lasing, sensing, or for photovoltaics.

The optical properties of the individual nanoparticles and the periodic structures are simulated using a commercial FDTD software (Lumerical). The simulations of the scattering and absorption efficiencies were done using a total-field scattered-field (TFSF) source with a broadband (400–700 nm) beam incident along the longitudinal axis of the silicon cylinder, which is defined as the $z$-axis. The TFSF source divides the simulation region into two concentric volumes, one central around the particle with the total fields and another external where only the scattered fields propagate. Power transmission monitors are positioned around and inside the TFSF source to calculate the absorption and scattering cross sections, respectively. The efficiency is then calculated dividing the former quantities by the geometrical cross section, i.e., the area of the cylinder perpendicular to the propagation vector $k$ of the incident field. Perfectly matched layer (PML) boundaries are implemented for every direction. Nonuniform meshes are used for the simulations and a 1 nm refinement mesh around the particle when monitoring the near-field intensities.

The simulations of the arrays were performed using periodic boundary conditions in the $x$- and $y$-directions, and PML boundaries in the $z$-axis. FDTD only deals with Cartesian grids and for a hexagonal array, the unit cell is nonprimitive. The illumination consists of a broadband (400–700 nm) beam, approximated by a plane wave, which is incident normal to the array plane (the $xy$-plane) from the substrate. Nonuniform meshes are used, and a 3 nm refinement mesh is used in the array plane, extending over a length equal to the height of the nanoparticle. Due to the very low absorption of silicon between 600 and 700 nm, the lattice resonances are long-lived. Therefore, we used long simulation times (5000 fs), especially for the QGMs (7000 fs) with auto shutoff levels of $10−6$ and even $10−7$⁠. This ensures the convergence of the simulations. To extract the transmission and the fields, several monitors are placed in the $xy$-plane at different $z$ positions.

To model the optical constants of Si, Lumerical uses a multicoefficient model with a polynomial fit of the real and imaginary components of the complex dielectric function in a specific wavelength range that is consistent with the Kramers–Kronig relations.

This work was financially supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) through the Gravitation program “Research Center for Integrated Photonics” and the Innovational Research Incentives Scheme (Vici Grant No. 680-47-628). We would like to thank M. Ramezani, R. Thottungal-Valapu, and A. G. Curto for useful discussions.