Multipolar Lattice Resonances in Plasmonic Finite-Size Metasurfaces

Abstract

Collective lattice resonances in regular arrays of plasmonic nanoparticles have attracted much attention due to a large number of applications in optics and photonics. Most of the research in this field is concentrated on the electric dipolar lattice resonances, leaving higher-order multipolar lattice resonances in plasmonic nanostructures relatively unexplored. Just a few works report exceptionally high-Q multipolar lattice resonances in plasmonic arrays, but only with infinite extent (i.e., perfectly periodic). In this work, we comprehensively study multipolar collective lattice resonances both in finite and in infinite arrays of Au and Al plasmonic nanoparticles using a rigorous theoretical treatment. It is shown that multipolar lattice resonances in the relatively large (up to 6400 nanoparticles) finite arrays exhibit broader full width at half maximum (FWHM) compared to similar resonances in the infinite arrays. We argue that our results are of particular importance for the practical implementation of multipolar lattice resonances in different photonics applications.

1. Introduction

Collective lattice resonances (CLRs) are modes inherent to plasmonic or all-dielectric nanoparticles (NPs) arranged in a regular lattice, and observed at wavelengths close to Rayleigh anomalies of the lattice. Hybridization between broad resonances localized on a single NP with discrete states associated with the lattice periodicity implies exceptionally high-quality factors of CLRs observed in a number of experimental works [1,2,3,4,5]. CLRs in arrays of NPs have a rich history since pioneering works in the early 1980s [6,7,8], with a revival in the 2000s [9,10,11,12,13] and subsequent developing of a variety of CLRs-assisted applications [14,15,16,17,18] in biosensing, lasing, color printing, fluorescence enhancement, strong coupling, and nonlinear optics.

CLRs in arrays of NPs are well-explained and understood within the coupled electric dipole (ED) [19,20], coupled ED and magnetic dipole (MD) [21], and coupled ED, MD, electric quadrupole (EQ) and magnetic quadrupole (MQ) [22] approximations. In all of these frameworks, closed-form analytical solutions are available for the extinction, scattering and absorption cross sections for infinite periodic structures. CLRs in finite arrays are almost exclusively studied within the coupled dipole approximation: for plasmonic systems, taking into account only ED interactions [23,24,25,26,27,28,29,30], and for all-dielectric arrays, taking into account ED and MD interactions [31,32,33,34]. Recent works on CLRs suggest a great promise of multipolar CLRs, obviously in arrays of high-index NPs [35] since such NPs support a richer variety of multipolar modes [36] compared to plasmonic counterparts. Nonetheless, multipolar CLRs in plasmonic nanostructures have been shown to emerge in arrays of Au [37,38,39,40], Al [41], and Au NPs [42,43] with promising applications in the Kerker effect [44] and lasing [45].

Since it is inevitable that the Q-factor of CLRs increases for higher-order multipolar resonances of the constituent scatterers involved in collective phenomena, multipolar CLRs hold a great promise. Similar to dipolar CLRs [25,29,33], the size of the illuminated area of the array is expected to play a critical role on multipolar CLRs, but there are no clear and systematic studies revealing these effects and answering a simple question: “how much is enough?” In this work, we comprehensively study multipolar CLRs in regular arrays of plasmonic NPs, both Au and Al, arranged in hexagonal and square lattices, respectively, and demonstrate how CLRs in finite arrays differ from their counterparts in infinite arrays.

2. Methods

We consider finite and infinite ($N=\infty$) 2D periodic arrays of identical NPs. In both cases, we compare respective extinction efficiencies, $Q_{\mathrm{ext}}$, which is in line with previous similar studies [25,29,30,33]. The arrays are illuminated from the top by the plane wave with normal incidence along the Z axis positive direction and polarization along the X axis; see Figure 1.

For both finite and infinite arrays, two different geometries are considered: square array of Al NPs, and a honeycomb array of Au NPs; see Figure 1. Depending on the application and spectral range of CLRs, square arrays of plasmonic NPs are widely considered in Ag [46], Au [5,47,48], Al [41,44,49,50,51] and TiN [52] NPs, whereas honeycomb lattices are considered for Ag [53,54,55,56] and Au [57] NPs, within ED approximation in almost all cases.

For infinite arrays, we relied on the well-established and commonly used finite-difference time-domain (FDTD) method implemented in a commercial FDTD Lumerical package [58]. Absorption efficiency was calculated using a set of monitors surrounding each NP in the unit cell, while scattering efficiency was calculated as a sum of fields recorded on two monitors located below and above the total field scattered field (TFSF) region [59] inside the simulation domain. Extinction cross section, ${\sigma }_{\mathrm{ext}}$, is defined as the sum of the absorption and scattering cross sections. Periodic boundary conditions have been applied at the lateral boundaries of the simulation box, while perfectly matched layer (PML) boundary conditions were used on the remaining top and bottom sides. An adaptive mesh has been used to reproduce accurately the nanosphere shape.

For finite arrays, we used the generalized Mie theory [60,61] and considered dipolar ($\ell =1$) and quadrupolar ($\ell =2$) interactions. Extinction cross section of the entire array for the ℓ-th multipole is defined as a sum of respective contributions of each i-th constituent particle:

(1)${\sigma }_{\mathrm{ext};\ell }=\sum _{i=1}^N{\sigma }_{\mathrm{ext};\ell }^i=\frac{4\pi }{k^2}\ell (\ell +1)(2\ell +1)\sum _{i=1}^N\sum _{m=-\ell }^{\ell }\frac{(\ell -m)!}{(\ell +m)!}\Re \left(p_{m\ell }^{i,i*}{a}_{m\ell }^i+q_{m\ell }^{i,i*}{b}_{m\ell }^i\right),$

where $a_{m\ell }^i$ and $b_{m\ell }^i$ are scattering coefficients, and $p_{m\ell }^{i,i}$, $q_{m\ell }^{i,i}$ are external field decomposition coefficients [60]. The summation in the last bracers in Equation (1) is across electric, $p_{m\ell }^{i,i*}{a}_{m\ell }^i$ (ED for $\ell =1$, EQ for $\ell =2$), and magnetic, $q_{m\ell }^{i,i*}{b}_{m\ell }^i$ (MD for $\ell =1$, MQ for $\ell =2$), polarizations. Respective phases of electric and magnetic multipoles are defined as:

(2)${\phi }_{\ell }=arg\left(\sum _{m=-\ell }^{\ell }\frac{(\ell -m)!}{(\ell +m)!}{p}_{m\ell }^{i,i*}{a}_{m\ell }^i\right)\mathrm{and}{\varphi }_{\ell }=arg\left(\sum _{m=-\ell }^{\ell }\frac{(\ell -m)!}{(\ell +m)!}{q}_{m\ell }^{i,i*}{b}_{m\ell }^i\right),$

correspondingly. The total extinction cross section is the sum of cross sections for dipolar, $\ell =1$, and quadrupolar, $\ell =2$, interactions: ${\sigma }_{\mathrm{ext}}={\sigma }_{\mathrm{ext};1}+{\sigma }_{\mathrm{ext};2}$. Finally, the extinction efficiency is understood as the ratio of the extinction cross section to the area, S, of the unit cell for infinite arrays, $Q_{\mathrm{ext}}={\sigma }_{\mathrm{ext}}/S$, and to the total area, $NS$, for finite arrays, $Q_{\mathrm{ext}}={\sigma }_{\mathrm{ext}}/NS$. The simulation time and memory demand for finite arrays increase as $N^2$, so we limited our consideration up to $N=6400$, which still can be simulated in a reasonable time.

3. Results

3.1. Aluminum Square Array

We start with square arrays of Al spherical NPs, see Figure 1a. The extinction spectra of arrays for a different amount N of constituent NPs are shown in Figure 2a. Two distinct resonances were observed at $\lambda =421$ nm and $\lambda =468$ nm. To understand the origin of these resonances, we plotted the decomposition of the spectra into multipoles (ED, EQ, MD, MQ) in Figure 2d. It can be seen that the main contribution at the $\lambda =468$ nm resonance was due to ED, while the $\lambda =421$ nm peak appearde due to MD and EQ with a minor contribution of ED interaction. We noticed that $n_h{h}_a=420$ nm, which means that the peak at $\lambda =421$ nm corresponded to a multipolar (hybrid MD+EQ) CLR coupled to the $[1,0]$ Rayleigh anomaly. The ED resonance in the finite array ($N\le 400$) at $\lambda =468$ nm had a lower amplitude compared with the infinite array; it slightly shifted to longer wavelengths and quickly converged to the infinite array already for $N=1600$. On the other hand, the multipolar CLR at $\lambda =421$ nm did not converge to the respective CLR in the infinite lattice, even for $N=6400$ NPs.

To obtain deeper insights on the resonances of the arrays, we plotted the spatial distribution of the phase of the scattered field and the extinction efficiency for each constituent NP of $40\times 40$ array at $\lambda =421$ nm (ED, MD, EQ modes) and at $\lambda =468$ nm (ED mode), see Figure 2b,c. Indeed, it can be observed that extinction of EQ and MD modes was the same for almost the entire array, except for a thin boundary region, where the ED contribution was quite large. Moreover, we emphasize that it is insufficient to track only “local extinction” [29] to infer the collective lattice origin of the resonance, but it is also instructive to consider the phase of the scattered field at each NP, as shown in Figure 2b,c. We can conclude that the field at MD and EQ resonances was in the same phase for almost every NP in the array, which highlights even more the collective nature of lattice resonances.

The evolution of the spatial distribution of the extinction efficiency and respective phases with N is shown in Figure 3. The pattern of spatial distribution of extinction per particle for the MD and EQ modes did not change with the increase of N. On the contrary, for the ED mode, the region with high extinction values shifted from the center of the array to the edges as the number of NPs increased, while for $N=6400$ the extinction and phases distributions became homogeneous. At the same time, the spatial distribution was characterized by the different pattern features. For MD and EQ it was due to the difference between the edges and the internal part of the array, while for ED there were the inhomogeneities on the scales of several lattice periods. It should be mentioned that features of the extinction distribution in the form of stripes were oriented along the polarization of the external electric field for ED and perpendicular to the polarization for the EQ and MD modes. For different multipoles, the spatial distribution became homogeneous at different sizes of the array: while for the ED mode the pattern became homogeneous only for $N=6400$, for the MD and EQ modes it was already homogeneous for $N=900$. The behavior of the phase distribution pattern with the increase of NPs number was similar to the extinction pattern. The latter can be explained by the fact that both CLRs in the considered structures emerged due to phase matching so that the phase distribution played a critical role in the formation of high amplitude patterns. This differs from the behavior of disordered systems of interacting particles, such as fractal aggregates, where the local field enhancement is due to the anisotropy of the spatial distribution of the particles [62].

3.2. Au Honeycomb Array

The extinction spectra of honeycomb arrays with different N are shown in Figure 4a, where two distinct resonances at $\lambda =760$ nm and at $\lambda =830$ nm can be observed. While the long-wavelength resonance had a pure ED nature, as shown in Figure 4b,d, the short-wavelength resonance was a superposition of MD and EQ resonances; see Figure 4c,d. As opposed to a square array of Al NPs, the ED interaction is completely suppressed at the short-wavelength resonance [42] in the honeycomb structure, cf. Figure 2d and Figure 4d. Notice that the MD contribution at $\lambda =760$ nm was surprisingly significant and cannot be neglected, as was suggested in Reference [42]. Similar to the case of a square array, a multipolar CLR in array with relatively large $N=5904$ did not converge to the respective CLR in the infinite array, while the ED resonance resembled an infinite lattice response for sufficiently small N.

Decompositions of the extinction on each NP and corresponding phases are shown in Figure 5. The spatial extinction distribution pattern was also different compared to the square array of Al NPs for both resonances. The solid region with a high extinction of ED resonance acquired a periodic pattern with a constant period that doid not depend on N, with decreasing contrast of the periodic pattern with the increase of NPs number. The behavior of hybrid multipolar EQ–MD resonance also differed from the square array of Al NPs. The single region with NPs having high extinction split into four separate regions for a large amount of NPs in the array, $N>2488$. In all cases, the phase distribution pattern determined the extinction pattern, since the CLR in our case were connected to Rayleigh anomalies and phase matching played crucial role in emergence of the resonances.

4. Discussion

For practical purposes, it is instructive to estimate the full width at half maximum (FWHM) and Q-factor of the respective resonances, which are shown in Figure 6 for multipolar hybrid CLRs. Points correspond to $\Delta \lambda$ (FWHM) and $\lambda /\Delta \lambda$ (Q-factor) obtained directly from spectra presented in Figure 2a and Figure 4a. The approximation as $\Delta \lambda =\Delta {\lambda }_{\inf }+k/\sqrt{N}$ [33] for the Al square array ($\Delta {\lambda }_{\inf }$ = 1.9 nm, $k=269$ nm) and Au honeycomb array ($\Delta {\lambda }_{\inf }=6.5$ nm, $k=184$ nm) were used to fit the data. The asymptotes of these approximations, which correspond to the respective quantities for infinite arrays, are shown with dashed lines. The asymptotic $\Delta {\lambda }_{\inf }=1.9$ nm for a square array is consistent with $\Delta {\lambda }_{\inf }=2.12$ nm obtained directly from the FDTD calculations for the infinite array. For the honeycomb array, the corresponding $\Delta {\lambda }_{\inf }=6.5$ nm was used for approximation and $\Delta {\lambda }_{\inf }=8.3$ nm was used for FDTD calculations. These results are in agreement with the lower limit for FWHM, which follows from the uncertainty relation for the electromagnetic wave ([33], Equation 5)). Figure 6 shows that the honeycomb lattice is more suitable for harnessing multipolar CLRs, since the values of FWHM and Q-factor in finite arrays are closer to those of infinite arrays compared to the case of a square lattice. Figure 6 makes it clear that the best answer to the question “how much is enough” is nothing but “∞”. However, we can use $\Delta \lambda =\Delta {\lambda }_{\inf }+k/\sqrt{N}$ approximation to estimate the N value required for the CLR in the finite array to resemble CLR in the infinite array. For instance, $N=8\times {10}^6$ for an Al square array and $N=3.2\times {10}^5$ for an Au honeycomb array provide FWHM within $5\%$ discrepancy from $\Delta {\lambda }_{\inf }$. This corresponds to ≈$0.8\times 0.8$ mm${}^2$ and ≈$0.4\times 0.36$ mm${}^2$ overall sizes of finite arrays.

Finally, the electric field distribution at the distance of $z=65$ nm of the NPs centers for MD–EQ CLRs at $\lambda =421$ nm for a square array and at $\lambda =760$ nm for a honeycomb array are shown in Figure 7. It can be seen that the inhomogeneities of the field over the finite structure are consistent with the extinction shown in Figure 3 and Figure 5.

5. Conclusions

To conclude, we have systematically studied multipolar collective lattice resonances in infinite and finite arrays of plasmonic nanoparticles using rigorous theoretical treatment. We showed that narrow hybrid EQ–MD multipolar resonances in finite arrays are broader than they are predicted to be for perfectly periodic (effectively infinite) arrays, while the broad ED resonance appears to be the same both in infinite lattices and in relative large finite arrays. After closer inspection of the role of each particle of the array on its total electromagnetic response, it was shown that phase distribution for the EQ mode has significant inhomogeneities, while for the MD and ED modes it is relatively uniform for large arrays, which explains observed spectral behavior for finite and infinite arrays. We believe that our results are important for experimental setups where the uniform electromagnetic response across nanoparticles in the array is essential, for example in surface-enhanced Raman scattering [63], sensing [64] and other applications. Moreover, we anticipate that similar research would be of great interest for arrays with more complex unit cells [65,66,67]. The obtained results set the ground for optimization problems, with a primary goal to harness narrowband infinite-like behavior in compact finite-size structures with engineered local geometrical parameters.

Author Contributions

Conceptualization, A.E.E., V.S.G., S.V.K., I.L.R.; methodology, A.E.E. and V.S.G., S.P.P.; software, A.E.E. and V.S.G.; validation, A.S.K.; investigation, A.S.K.; writing—original draft preparation, A.S.K.; writing—review and editing, V.S.G., I.L.R., S.P.P.; visualization, A.S.K.; funding acquisition, A.E.E. All authors have read and agreed to the published version of the manuscript.

Funding

The reported study was funded by the Russian Science Foundation project number 19-72-00066.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available upon reasonable request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Schematic representation of (a) square array of Al nanoparticles (NPs) with ha=280 nm, RAl=60 nm, medium refractive index nh=1.5 and (b) honeycomb array of Au NPs with hb=600 nm, RAu=120 nm, nh=1.45. The red rectangle in (a) and red hexagon in (b) denote primitive cells for corresponding structures. The blue rectangle in (b) highlights the simulation domain with periodic boundary conditions for finite-difference time-domain (FDTD) calculations. Notice that the simulation domain and primitive cell coincide for the square array in (a). The unit cell for the honeycomb array consist of two NPs but in case of an infinite array simulated using the FDTD method, four NP comprising unit cells were used to reproduce infinity with orthogonal translation vectors (b).

Figure 2. (a) Extinction spectra for arrays of N spherical Al NPs with RAl=60 nm arranged in a square lattice with ha=280 nm. (b,c) Extinction efficiency per particle, corresponding phases of scattered wave and its decomposition at λ=421 nm (ED, EQ, MD) and λ=468 nm (ED) for N=1600. (d) Total extinction spectra and the respective contributions of multipolar modes (ED, MD, EQ, MQ) for arrays with different N, as labeled in the legend. ED: electric dipole; MD: magnetic dipole; EQ: electric quadrupole; MQ: magnetic quadrupole.

Figure 3. Evolution of the extinction efficiency per NP and scattered wave phases for different N for ED, MD and EQ resonances. Highlighted N=1600 case corresponds to data shown in Figure 2b,c.

Figure 4. (a) Extinction spectra for honeycomb arrays of N spherical Au NPs with RAu=120 nm and period (distance between NPs centers) hb=600 nm. (b,c) Extinction efficiency per particle, corresponding phases of scattered wave and its decomposition for resonances at λ=760 nm (ED, MD, EQ) and λ=830 nm (ED) for N=2488. (d) Total extinction spectra and the respective contributions of multipolar modes (ED, MD, EQ, MQ) for the arrays with different N, as labeled in the legend.

Figure 5. Evolution of the extinction efficiency per NP and scattered wave phases for different N for ED, MD, and EQ resonances. Highlighted N=2488 case corresponds to data shown in Figure 4b,c.

Figure 6. (a) Full width at half maximum (FWHM), Δλ, and (b) quality factor, Q=λ/Δλ, of a multipolar collective lattice resonance (CLR) at λ=421 nm for an Al square array and at λ=760 nm for an Au honeycomb array. Dashed lines show the estimate of the respective quantities for infinite arrays ([33], Equation (5)).

Figure 7. Electric field distribution at x–y plane with z=65 nm for (a) Al square array of N=1600 NPs and (b) Au honeycomb array with N=2488 NPs.