Plasmons in two-dimensional lattices of near-field coupled nanoparticles

TL;DR

This paper develops an analytical framework to understand collective plasmonic modes in two-dimensional nanoparticle lattices, accounting for topological properties, environmental coupling, and retarded interactions, with implications for experimental observation.

Contribution

It provides exact solutions for quasistatic dispersion relations in arbitrary 2D nanoparticle lattices, including non-Bravais and topologically nontrivial structures, incorporating environmental effects.

Findings

Analytical expressions for plasmonic dispersion in various lattice geometries.

Closed-form formulas for radiative and nonradiative damping rates.

Insights into how retarded interactions modify band structures.

Abstract

We consider plasmonic metasurfaces constituted by an arbitrary periodic arrangement of spherical metallic nanoparticles. Each nanoparticle supports three degenerate dipolar localized surface plasmon (LSP) resonances. In the regime where the interparticle distance is much smaller than the optical or near-infrared wavelength associated with the LSPs, the latter couple through the dipole-dipole interaction and form collective plasmonic modes which extend over the whole metasurface. Within a Hamiltonian model which we solve exactly, we derive general expressions which enable us to extract analytically the quasistatic plasmonic dispersion for collective modes polarized within the plane and perpendicular to the plane of the metasurface. Importantly, our approach allows us not only to consider arbitrary Bravais lattices, but also non-Bravais two-dimensional metacrystals featuring nontrivial topological properties, such as, e.g., the honeycomb or Lieb lattices. Additionally, using an open quantum system approach, we consider perturbatively the coupling of the collective plasmons to both photonic and particle-hole environments, which lead, respectively, to radiative and nonradiative frequency shifts and damping rates, for which we provide closed-form expressions. The radiative frequency shift, when added to the quasistatic dispersion relation, provides an approximate analytical description of the fully retarded band structure of the collective plasmons. While it is tempting to make a direct analogy between the various systems which we consider and their electronic tight-binding equivalents, we critically examine how the long-range retarded and anisotropic nature of the dipole--dipole interaction may quantitatively and qualitatively modify the underlying band structures and discuss their experimental observability.

Figures

Equations

$$H=H_{\mathrm{pl}}+H_{\mathrm{ph}}+H_{\mathrm{eh}}+H_{\mathrm{pl}\text{-}\mathrm{ph}}+H_{\mathrm{pl}\text{-}\mathrm{eh}},$$

$$H_{\mathrm{pl}}=H_{\mathrm{pl}}^0+H_{\mathrm{pl}}^{\mathrm{int}}$$

$$H_{\mathrm{pl}}^0=\sum _s\sum _{{\mathbf{R}}_s}\sum _{\sigma }\left\{\frac{{\left[{\Pi }_s^{\sigma }({\mathbf{R}}_s)\right]}^2}{2N_{\mathrm{e}}{m}_{\mathrm{e}}}+\frac{N_{\mathrm{e}}{m}_{\mathrm{e}}}{2}{\omega }_0^2{\left[h_s^{\sigma }({\mathbf{R}}_s)\right]}^2\right\},$$

$$H_{\mathrm{pl}}^{\mathrm{int}}=$$

$$\times \frac{{\delta }_{\sigma {\sigma }^{\prime }}-3(\hat{\sigma }\cdot {\hat{\rho }}_{s{s}^{\prime }})({\hat{\sigma }}^{\prime }\cdot {\hat{\rho }}_{s{s}^{\prime }})}{{\left|{\mathbf{R}}_s-{\mathbf{R}}_{s^{\prime }}^{\prime }\right|}^3},$$

$$b_s^{\sigma }({\mathbf{R}}_s)=\sqrt{\frac{N_{\mathrm{e}}{m}_{\mathrm{e}}{\omega }_0}{2\hslash }}{h}_s^{\sigma }({\mathbf{R}}_s)+\mathrm{i}\frac{{\Pi }_s^{\sigma }({\mathbf{R}}_s)}{\sqrt{2N_{\mathrm{e}}{m}_{\mathrm{e}}\hslash {\omega }_0}}$$

$$H_{\mathrm{pl}}=$$

$$+\frac{\hslash \Omega }{2}\sum _{s{s}^{\prime }}\sum _{\begin{array}{c}{\mathbf{R}}_s{\mathbf{R}}_{s^{\prime }}^{\prime }\\ ({\mathbf{R}}_s\ne {\mathbf{R}}_{s^{\prime }}^{\prime })\end{array}}\sum _{\sigma {\sigma }^{\prime }}\left[b_s^{\sigma }({\mathbf{R}}_s)+{b_s^{\sigma }}^{†}({\mathbf{R}}_s)\right]$$

$$\times \left[b_{s^{\prime }}^{{\sigma }^{\prime }}({\mathbf{R}}_{s^{\prime }}^{\prime })+{b_{s^{\prime }}^{{\sigma }^{\prime }}}^{†}({\mathbf{R}}_{s^{\prime }}^{\prime })\right]$$

$$\times \frac{{\delta }_{\sigma {\sigma }^{\prime }}-3(\hat{\sigma }\cdot {\hat{\rho }}_{s{s}^{\prime }})({\hat{\sigma }}^{\prime }\cdot {\hat{\rho }}_{s{s}^{\prime }})}{{\left(|{\mathbf{R}}_s-{\mathbf{R}}_{s^{\prime }}^{\prime }|/d\right)}^3},$$

$$\Omega =\frac{{\omega }_0}{2}{\left(\frac{a}{d}\right)}^3.$$

$$H_{\mathrm{ph}}=\sum _{\mathbf{k},{\hat{\lambda }}_{\mathbf{k}}}\hslash {\omega }_{\mathbf{k}}{a_{\mathbf{k}}^{{\hat{\lambda }}_{\mathbf{k}}}}^{†}{a}_{\mathbf{k}}^{{\hat{\lambda }}_{\mathbf{k}}},$$

$$H_{\mathrm{pl}\text{-}\mathrm{ph}}=$$

$$+\frac{N_{\mathrm{e}}{e}^2}{2m_{\mathrm{e}}{c}^2}{\mathbf{A}}^2({\mathbf{R}}_s)],$$

$${\mathbf{\Pi }}_s({\mathbf{R}}_s)=$$

$$\mathbf{A}({\mathbf{R}}_s)=\sum _{\mathbf{k},{\hat{\lambda }}_{\mathbf{k}}}{\hat{\lambda }}_{\mathbf{k}}\sqrt{\frac{2\pi \hslash c^2}{\mathcal{V}{\omega }_{\mathbf{k}}}}\left(a_{\mathbf{k}}^{{\hat{\lambda }}_{\mathbf{k}}}{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot {\mathbf{R}}_s}+{a_{\mathbf{k}}^{{\hat{\lambda }}_{\mathbf{k}}}}^{†}{\mathrm{e}}^{-\mathrm{i}\mathbf{k}\cdot {\mathbf{R}}_s}\right)$$

$$H_{\mathrm{eh}}=\sum _s\sum _{{\mathbf{R}}_s}\sum _i{ϵ}_{{\mathbf{R}}_{s}i}{c}_{{\mathbf{R}}_{s}i}^{†}{c}_{{\mathbf{R}}_{s}i}^{},$$

$$H_{\mathrm{pl}\text{-}\mathrm{eh}}=$$

$$\times ⟨{\mathbf{R}}_{s}i|\sigma |{\mathbf{R}}_{s}j⟩c_{{\mathbf{R}}_{s}i}^{†}{c}_{{\mathbf{R}}_{s}j}^{},$$

$$\Lambda =\sqrt{\frac{\hslash m_{\mathrm{e}}{\omega }_0^3}{2N_{\mathrm{e}}}}.$$

$${\beta }_{\tau }^{\epsilon }(\mathbf{q})=\sum _{s\sigma }\left[u_{\tau s}^{\epsilon \sigma }(\mathbf{q})b_s^{\sigma }(\mathbf{q})+v_{\tau s}^{\epsilon \sigma }(\mathbf{q}){b_s^{\sigma }}^{†}(-\mathbf{q})\right]$$

$$H_{\mathrm{pl}}=\sum _{\mathbf{q}}\sum _{\tau \epsilon }{H}_{\tau }^{\epsilon }(\mathbf{q}),H_{\tau }^{\epsilon }(\mathbf{q})=\hslash {\omega }_{\tau }^{\epsilon }(\mathbf{q}){{\beta }_{\tau }^{\epsilon }}^{†}(\mathbf{q}){\beta }_{\tau }^{\epsilon }(\mathbf{q}),$$

$$\displaystyle\big{[}\omega_{0}$$

$$\displaystyle+\Omega\sum_{s^{\prime}\sigma^{\prime}}\big{[}u_{\tau s^{\prime}}^{\varepsilon\sigma^{\prime}}(\mathbf{q})-v_{\tau s^{\prime}}^{\varepsilon\sigma^{\prime}}(\mathbf{q})\big{]}f_{s^{\prime}s}^{\sigma^{\prime}\sigma}(\mathbf{q})=0$$

$$\displaystyle-\big{[}\omega_{0}$$

$$\displaystyle+\Omega\sum_{s^{\prime}\sigma^{\prime}}\big{[}u_{\tau s^{\prime}}^{\varepsilon\sigma^{\prime}}(\mathbf{q})-v_{\tau s^{\prime}}^{\varepsilon\sigma^{\prime}}(\mathbf{q})\big{]}f_{s^{\prime}s}^{\sigma^{\prime}\sigma}(\mathbf{q})=0.$$

$$f_{s{s}^{\prime }}^{\sigma {\sigma }^{\prime }}(\mathbf{q})=$$

$$\times \left[{\delta }_{\sigma {\sigma }^{\prime }}-3(\hat{\sigma }\cdot {\hat{\rho }}_{s{s}^{\prime }})({\hat{\sigma }}^{\prime }\cdot {\hat{\rho }}_{s{s}^{\prime }})\right]$$

$${\varphi }_{\tau }^{\epsilon }(\mathbf{q})=\arccos \left(\left|{\hat{u}}_{\tau }^{\epsilon }(\mathbf{q})\cdot \hat{q}\right|\right).$$

$${\omega }^z(\mathbf{q})={\omega }_0\sqrt{1+2\frac{\Omega }{{\omega }_0}{f}^{zz}(\mathbf{q})},$$

$${\omega }^{{\epsilon }_{\parallel ,\pm }}(\mathbf{q})=$$

$$\displaystyle\pm\sqrt{\left[f^{xx}(\mathbf{q})-f^{yy}(\mathbf{q})\right]^{2}+4\left[f^{xy}(\mathbf{q})\right]^{2}}\bigg{]}\bigg{\}}^{1/2}$$

$$\omega_{\tau}^{z}(\mathbf{q})=\omega_{0}\sqrt{1+2\frac{\Omega}{\omega_{0}}\big{[}f_{11}^{zz}(\mathbf{q})+\tau|f_{12}^{zz}(\mathbf{q})|\big{]}},\quad\tau=\pm 1.$$

$${\omega }_{\tau }^z(\mathbf{k})\simeq {\omega }_0-\Omega |f_{11}^{zz}(\mathbf{K})|+\tau v^z|\mathbf{k}|,\tau =\pm 1,$$

$${\mathcal{H}}_{\mathbf{k}}^{\mathrm{eff}}=\left[\hslash {\omega }_0-\hslash \Omega |f_{11}^{zz}(\mathbf{K})|\right]1_2-\hslash v^z\mathbf{\sigma }\cdot \mathbf{k},$$

$${\sigma }_x=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right){\sigma }_y=\left(\begin{array}{cc}0& -\mathrm{i}\\ \mathrm{i}& 0\end{array}\right){\sigma }_z=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)$$

Full text

I Introduction

The interaction of light with a small metallic particle results in collective electronic modes termed localized surface plasmons (LSPs) Bertsch and Broglia (1994); Kreibig and Vollmer (1995). In the case where the wavelength of the incoming light is much larger than the nanostructure itself, the LSP corresponds to a dipolar collective oscillation of the electronic cloud against the inert ionic background. While such a phenomena was empirically discovered centuries ago by late Romans Freestone et al. (2007), the underlying physics was only theorized by Mie at the beginning of the 20th century, who solved Maxwell’s equations for a metallic sphere embedded in a dielectric medium Mie (1908); Born and Wolf (1980); Bohren and Huffman (2004). Since then, plasmonic nanostructures have attracted a surge of interest due to their ability to perform subwavelength optics by confining the electromagnetic field to nanometric regions Barnes et al. (2003); Maier (2007); Gramotnev and Bozhevolnyi (2010). Due to the extreme sensitivity of the LSP resonance frequency to the nanoparticle size, shape, material, and dielectric environment Kreibig and Vollmer (1995); Kelly et al. (2003), a wealth of appealing technological applications have risen from the field of nanoplasmonics, such as, e.g., chemical Mayer and Hafner (2011) and biological Anker et al. (2008) sensors.

When two metallic nanoparticles are positioned in close vicinity of each other (i.e., separated by a distance typically smaller than the LSP wavelength) so that they form a dimer, another factor crucially influencing the resonance frequencies of the latter is the Coulomb interaction between the LSPs. The quasistatic dipole–dipole interaction, which decays with the interparticle distance $d$ as $1/d^{3}$, gives rise to coupled modes, often coined “hybridized” modes Jain and El-Sayed (2010), which correspond to symmetric (in-phase) or antisymmetric (out-of-phase) configurations of the oscillating electric dipolar moments on each nanoparticle. For transverse-polarized modes (with respect to the axis formed by the dimer), the high- (low-) energy plasmonic state corresponds to an in-phase (out-of-phase) configuration. Conversely, for longitudinal modes, the low- (high-) energy state corresponds to aligned (antialigned) dipole moments. The splitting in frequencies between these hybridized modes scales with the interparticle distance as $1/d^{3}$, and can be spectroscopically resolved as long as the linewidth (which is of both radiative and nonradiative nature) of the two resonance peaks is somewhat smaller than the above-mentioned splitting. The picture above is valid as long as the two nanoparticles are not too close to each other, so that higher multipolar modes do not mix with the dipolar ones Prodan et al. (2003); Prodan and Nordlander (2004); Nordlander et al. (2004); Park and Stroud (2004), and the quantum tunneling of electronic charges between the two particles can be disregarded, such that so-called charge transfer plasmons are irrelevant Esteban et al. (2012); Savage et al. (2012); Scholl et al. (2013); Esteban et al. (2015); Rossi et al. (2015); Zhu et al. (2016). Since the pioneering work by Ruppin Ruppin (1982), who extended Mie’s theory Mie (1908); Born and Wolf (1980); Bohren and Huffman (2004) to two nearby metallic spheres embedded in a dielectric medium, hybridized plasmonic modes in nanoparticle dimers have been investigated in numerous experimental Jain and El-Sayed (2010); Tamaru et al. (2002); Rechberger et al. (2003); Olk et al. (2008); Chu et al. (2009); Koh et al. (2009); Barrow et al. (2014); Savage et al. (2012); Scholl et al. (2013) and theoretical works Gérardy and Ausloos (1983); Nordlander et al. (2004); Dahmen et al. (2007); Bachelier et al. (2008); Zuloaga et al. (2009); Esteban et al. (2012); Zhang et al. (2014); Brandstetter-Kunc et al. (2015, 2016); Downing et al. (2017); Savage et al. (2012); Scholl et al. (2013); Esteban et al. (2015); Rossi et al. (2015).

In periodic arrays of near-field coupled nanoparticles, the dipolar interaction between the LSPs leads to collective modes that are extended over the whole lattice. In one-dimensional (1D) chains of regularly spaced nanoparticles, such collective plasmons were extensively studied both at the theoretical Brongersma et al. (2000); Maier et al. (2003); Park and Stroud (2004); Weber and Ford (2004); Citrin (2004); Simovski et al. (2005); Citrin (2006); Koenderink and Polman (2006); Markel and Sarychev (2007); Fung and Chan (2007); Lee et al. (2012); Compaijen et al. (2013); del Pino et al. (2014); Compaijen et al. (2015); Petrov (2015); Brandstetter-Kunc et al. (2016); Downing et al. (2018); Compaijen et al. (2018) and experimental Krenn et al. (1999); Maier et al. (2002, 2003); Koenderink et al. (2007); Crozier et al. (2007); Apuzzo et al. (2013); Barrow et al. (2014); Gür et al. (2018) levels, since these systems may serve as plasmonic waveguides where plasmon–photon hybrid modes (so-called plasmon polaritons) are laterally confined to subwavelength scales and can possibly propagate over macroscopic distances. The importance of retardation effects in the dipolar interaction, which become relevant for nanoparticles in the chain spaced by a distance of the order of the LSP wavelength, was put forward in Refs. Weber and Ford (2004); Citrin (2004); Simovski et al. (2005); Citrin (2006); Koenderink and Polman (2006); Markel and Sarychev (2007); Fung and Chan (2007); Compaijen et al. (2015); Petrov (2015); Downing et al. (2018); Compaijen et al. (2018). In particular, it was shown that retardation leads to a pronounced discontinuity in the dispersion relation of the collective plasmons polarized transversely to the chain for wave vectors corresponding to the intersection of the light cone with the quasistatic band structure. The longitudinal dispersion relation remains however continuous, but presents discontinuities in its derivative in the vicinity of the above-mentioned crossing. Similar effects also occur in related condensed matter systems, such as exciton-polaritons in quantum wires Citrin (1992). The crucial role played by radiative and absorption losses on the propagation of plasmonic waves along the nanoparticle chain was also studied in detail in the previous works of Refs. Brongersma et al. (2000); Maier et al. (2003); Park and Stroud (2004); Weber and Ford (2004); Citrin (2004); Simovski et al. (2005); Citrin (2006); Koenderink and Polman (2006); Markel and Sarychev (2007); Fung and Chan (2007); Lee et al. (2012); Compaijen et al. (2013); del Pino et al. (2014); Compaijen et al. (2015); Petrov (2015); Compaijen et al. (2018); Downing et al. (2018). Notably, Ref. Brandstetter-Kunc et al. (2016) showed that the nonradiative Landau damping, that is, the disintegration of the collective plasmons into particle-hole pairs, is of primarily importance as it dominates the plasmon linewidth for nanoparticles of only a few nanometers in size.

Recently, dimerized Ling et al. (2015); Downing and Weick (2017); Gómez et al. (2017); Pocock et al. (2018); Downing and Weick (2018); Pocock et al. (2019); Downing et al. (2019) as well as zig zag chains of nanoparticles Poddubny et al. (2014); Sinev et al. (2015) were proposed as a plasmonic analog of the celebrated Su-Schrieffer-Heeger (SSH) model Su et al. (1979, 1980); Heeger et al. (1988) presenting nontrivial topologically protected edge states. In particular, the robustness of such topological states against the long range retarded dipolar interactions was discussed in Refs. Pocock et al. (2018); Downing and Weick (2018); Pocock et al. (2019).

The extension of the concepts introduced above to two spatial dimensions offers new exciting possibilities. Metasurfaces, that is, two-dimensional (2D) periodic arrangements of subwavelength metallic nanostructures, are indeed at present a very active field of research, as they enable one to tailor light in a way that goes far beyond what can be achieved with conventional optics. Thus far, the vast majority of the litterature on plasmonic metasurfaces (see, e.g., the review articles of Refs. Meinzer et al. (2014); Glybovski et al. (2016); Kravets et al. (2018); Cherqui et al. (2019) and references therein) focused on the regime where the separation distance between each resonant element is of the order of the LSP wavelength, as this can be experimentally achieved with nowadays nanofabrication techniques. In this regime, the diffractive electromagnetic far fields generated by the essentially noninteracting nanoparticles of the array interfere and give rise to so-called surface lattice resonances (SLRs). The latter are of particular interest since they lead to much narrower absorption lines as that of the individual constituents of the metasurface, as well as angle-dependent dispersions, as was theoretically predicted in Refs. Carron et al. (1986); Zou et al. (2004) and later experimentally verified in Refs. Kravets et al. (2008); Auguié and Barnes (2008); Chu et al. (2008). Further works have demonstrated the use of SLRs in tailoring frequency stop gaps Rodriguez et al. (2011) and are of relevance to applications in light emission Rodriguez et al. (2012); Giannini et al. (2010). Genuinely quantum-mechanical effects Tame et al. (2013), such as the exciting perspective of lasing Zhou et al. (2013); Schokker and Koenderink (2014); Hakala et al. (2017); Guo et al. (2019), as well as Bose–Einstein condensation Martikainen et al. (2014); Hakala et al. (2018), have also been demonstrated in 2D plasmonic lattices. Notably, the works of Refs. Humphrey and Barnes (2014); Guo et al. (2017) combining modeling and experiments have brought attention to the role of the geometrical arrangement of the nanostructures composing the metasurface on the SLR properties.

In the present work, we focus on the less explored case of near-field coupled nanoparticles supporting dipolar LSPs in metasurfaces. In this case, the stronger dipolar coupling between LSPs, as compared to the weak diffractive couplings encountered in SLRs, can exhibit potentially interesting analogies with atomically thin, 2D materials, such as graphene Castro Neto et al. (2009) or transition metal dichalcogenides Mak and Shan (2016), where the electronic band structures are usually well described by tight-binding calculations. Metasurfaces composed of near-field nanoparticles may indeed present appealing nontrivial features in their band structure, paving the way to topological photonics performed with subwavelength elements Lu et al. (2014); Ozawa et al. (2019); Rider et al. (2019). For instance, it was theoretically demonstrated that a honeycomb lattice of plasmonic nanoparticles that are near-field coupled present chiral massless Dirac-like bosonic collective excitations Han et al. (2009); Weick et al. (2013); Sturges et al. (2015) which behave as electrons in graphene Castro Neto et al. (2009). Such a honeycomb lattice further hosts topologically protected edge states Wang et al. (2016). Exotic, so-called type-II Dirac plasmon-polaritons presenting a fully tunable tilted conical dispersion, were also recently unveiled in Ref. Mann et al. (2018). Straining the honeycomb metasurface and modifying the electromagnetic environment was recently shown to induce tunable pseudo-magnetic fields for polaritons Mann et al. .

Due to the vast number of possible 2D lattices of near-field coupled plasmonic nanoparticles with potentially interesting properties in their band structure, here we develop a general theoretical framework which enables us to consider the plasmonic properties of arbitrary metasurfaces. Our open quantum system approach, which builds on previous works on plasmonic dimers Brandstetter-Kunc et al. (2015); Downing et al. (2017) and chains Brandstetter-Kunc et al. (2016); Downing et al. (2018), allows us to unveil analytical expressions for the quasistatic plasmonic dispersion relations for collective modes polarized parallel and perpendicular to the plane of the metasurface. By considering the coupling of the purely plasmonic modes to photons of the electromagnetic vacuum, we also consider perturbatively the effects of retardation in the light–matter interaction, and we show that such retardation effects play a crucial role on the plasmonic band structure. Our approach, that straightforwardly include Ohmic losses, further gives access to the radiative lifetime of the plasmonic modes, which we evaluate analytically. Importantly, we also consider the decay of the collective plasmons into electron-hole pairs and show that the resulting Landau damping can be as significant as it is in single nanoparticles which are only a few nanometers in size.

The analytical treatment of the light–matter interaction has been shown Downing et al. (2018) in the case of chains to provide a rather good account of the results stemming from laborious numerical calculations based on fully retarded solutions to Maxwell’s equations for both the plasmonic dispersion relation and the radiative linewidth. In particular, the discontinuity presented by the dispersion relation (by its derivative) for the transverse (longitudinal) modes is qualitatively captured by the perturbative open quantum system approach. When we compare the results of the present study for 2D arrays of nanoparticles to existing numerical classical calculations for the square Zhen et al. (2008) and the honeycomb lattices Han et al. (2009), we further obtain good agreement.

Our model is not only applicable to plasmonic nanoparticles, but is also of relevance for any 2D system of regularly spaced subwavelength resonators coupled through the dipole–dipole interaction within the near-field regime, such as, e.g., dielectric nanoparticles where Ohmic losses are notably reduced Du et al. (2011); Slobozhanyuk et al. (2015); Bakker et al. (2017), microwave resonators Mann et al. (2018), magnonic microspheres Pirmoradian et al. (2018), cold atoms Perczel et al. (2017), or supercrystals made of semiconducting quantum dots Baimuratov et al. (2013).

This paper is organized as follows: Section II presents our Hamiltonian model and the open quantum system approach which we use to study collective plasmons in a generic 2D array of interacting spherical metallic nanoparticles. In Sec. III we present the diagonalization procedure of the purely plasmonic Hamiltonian, which gives access to the quasistatic dispersion relation of the collective modes. The latter is then extensively discussed, including the cases of Bravais and non-Bravais lattices. In Secs. IV and V, we consider the effects of the photonic environment alone which we treat perturbatively and present our results for the dispersion relation of the plasmonic modes including retardation effects (Sec. IV), as well as the corresponding radiative lifetimes (Sec. V). Section VI then analyzes the (possibly crucial) role played by the electronic environment onto the collective modes and presents our results for their Landau damping decay rates, as well as their associated electronic frequency shifts. In Sec. VII, we discuss the experimental observability of the plasmonic modes, before we conclude in Sec. VIII.

II Model

We consider an ensemble of interacting spherical metallic, nonmagnetic nanoparticles of radius $a$ forming an arbitrary 2D Bravais lattice with a basis. The array is characterized by the vectors $\mathbf{R}=n\mathbf{t}_{1}+m\mathbf{t}_{2}$ forming the Bravais lattice. Here, $\mathbf{t}_{1}$ and $\mathbf{t}_{2}$ are the primitive lattice vectors, while $n\in[0,\mathcal{N}_{1}]$ and $m\in[0,\mathcal{N}_{2}]$ are integer numbers with $\mathcal{N}_{1}$ ($\mathcal{N}_{2}$) the number of unit cells in the $\mathbf{t}_{1}$ ($\mathbf{t}_{2}$) direction. The array is composed of $\mathcal{S}$ sublattices, and the nanoparticles belonging to the sublattice $s=1,\ldots,\mathcal{S}$ are located at $\mathbf{R}_{s}=\mathbf{R}+\mathbf{d}_{s},$ where $\mathbf{d}_{s}$ is the vector belonging to the $xy$ plane and connecting the sublattice $s$ to $\mathbf{R}$ (see Fig. 1). By convention, we set $\mathbf{d}_{1}=0$ in the remainder of the paper.

Each nanoparticle forming the array supports three degenerate orthogonal dipolar LSP resonances polarized along the $\sigma=x$, $y$, and $z$ directions and characterized by the frequency $\omega_{0}$. Neglecting the effect of the embedding medium, as well as the screening of the valence electrons by the core electrons, one has $\omega_{0}=\omega_{\mathrm{p}}/\sqrt{3}$ Kreibig and Vollmer (1995), where $\omega_{\mathrm{p}}=\sqrt{4\pi n_{\mathrm{e}}e^{2}/m_{\mathrm{e}}}$ is the plasma frequency (throughout this paper, we use cgs units). Here, $-e$ ($<0$) is the electron charge, $m_{\mathrm{e}}$ its mass, and $n_{\mathrm{e}}$ the electronic density of the metal.

The dipolar LSPs supported by the nanoparticles in the array interact with their neighbors through the Coulomb interaction. Such a coupling gives rise to collective plasmonic modes which extend over the whole metasurface Meinzer et al. (2014). In this work we focus on the subwavelength regime where $d\ll k_{0}^{-1}$, with $d$ the center-to-center nearest neighbor interparticle distance and $k_{0}=\omega_{0}/c$, with $c$ the speed of light in vacuum. Furthermore, we assume that the interparticle distance $d\gtrsim 3a$. In this parameter regime, Park and Stroud Park and Stroud (2004) have shown by means of an exact diagonalization procedure based on a generalized tight-binding approach including multipoles up to $\ell=80$ (with $\ell$ the angular momentum) that the dominant interaction between two nanoparticles is the near-field quasistatic electric dipole–dipole ($\ell=1$) interaction which decays as $1/d^{3}$. For smaller interparticle distances ($2a\leqslant d\lesssim 3a$), higher multipolar interactions with $\ell>1$ Park and Stroud (2004) as well as tunneling phenomena Zhu et al. (2016) play an important role for the description of the collective modes, and we do not enter this regime of parameters in the remainder of the paper. Notice also that, since we consider nonmagnetic particles, and since our model does not incorporate an external magnetic field, magnetoplasmonic modes (see, e.g., Ref. Weick and Weinmann (2011)) which would result in oscillating magnetic dipolar moments are irrelevant to the present study.

Similarly to individual LSPs, the extended plasmonic modes supported by the metasurface are coupled to both a photonic and an electronic environment. The collective plasmons are indeed coupled to vacuum electromagnetic modes through the light–matter interaction, giving rise to finite radiative lifetimes as well as radiative frequency shifts, stemming from the retarded part of the dipole–dipole interaction Downing et al. (2017, 2018). Moreover, such collective modes are coupled to electron-hole excitations inside the nanoparticles, leading to the nonradiative Landau damping and an additional frequency renormalization.

We write the full Hamiltonian of the system as

[TABLE]

where $H_{\mathrm{pl}}$ describes the purely plasmonic degrees of freedom, while $H_{\mathrm{ph}}$ and $H_{\mathrm{eh}}$ correspond to the photonic and electronic environments, respectively. In Eq. (1), $H_{\mathrm{pl\textrm{-}ph}}$ and $H_{\mathrm{pl\textrm{-}eh}}$ are the coupling Hamiltonians of the plasmonic subsystem to photons and electron-hole pairs, respectively.

Within the Coulomb gauge Cohen-Tannoundji et al. (1992); Craig and Thirunamachandran (1984), the plasmonic Hamiltonian

[TABLE]

describing the near-field coupled LSPs is characterized by a noninteracting and an interacting term Weick et al. (2013); Brandstetter-Kunc et al. (2016). The noninteracting part related to individual nanoparticles reads

[TABLE]

with $h_{s}^{\sigma}(\mathbf{R}_{s})$ the $\sigma$ component of the displacement field associated with the dipole moment $\mathbf{p}(\mathbf{R}_{s})=-eN_{\mathrm{e}}\sum_{\sigma}h_{s}^{\sigma}(\mathbf{R}_{s})\hat{\sigma}$ of a single LSP located at position $\mathbf{R}_{s}$, while $\Pi_{s}^{\sigma}(\mathbf{R}_{s})$ is the momentum conjugated to $h_{s}^{\sigma}(\mathbf{R}_{s})$. Here and in what follows, hats designate unit vectors. In Eq. (3), $N_{\mathrm{e}}$ is the total number of valence electrons in each nanoparticle. The interacting term in Eq. (2) arises from the quasistatic dipole–dipole interaction and reads

[TABLE]

where $\boldsymbol{\rho}_{ss^{\prime}}=\mathbf{R}^{\prime}_{s^{\prime}}-\mathbf{R}_{s}$.

As we deal with nanoparticles of small sizes, quantum-size effects (such as Landau damping) can be important. Moreover, a quantum treatment of the plasmonic degrees of freedom provides a self-contained description of the light–matter interaction Downing et al. (2017, 2018). In view of the analysis of these effects, we hence write the plasmonic Hamiltonian (2) in terms of the bosonic ladder operators

[TABLE]

that annihilate an LSP at position $\mathbf{R}_{s}$ on sublattice $s$ with polarization $\sigma=x,y,z$ and its adjoint ${b_{s}^{\sigma}}^{\dagger}(\mathbf{R}_{s})$. Together with Eqs. (3) and (4), Eq. (2) thus takes the form

[TABLE]

with the coupling constant

[TABLE]

Note that $\Omega\ll\omega_{0}$ since we consider interparticle distances $d\geqslant 3a$.

The Hamiltonian (II) displays some similarities with a tight-binding Hamiltonian of an electronic 2D system Grosso and Parravicini (2014). The first term on the right-hand side of Eq. (II) [$\propto{b_{s}^{\sigma}}^{\dagger}(\mathbf{R}_{s})b_{s}^{\sigma}(\mathbf{R}_{s})$] corresponds to a fixed on-site energy, while the resonant terms [$\propto{b_{s}^{\sigma}}^{\dagger}(\mathbf{R}_{s})b_{s^{\prime}}^{\sigma^{\prime}}(\mathbf{R}^{\prime}_{s^{\prime}})$] in the second term describe the creation of an LSP at the lattice site $\mathbf{R}_{s}$ together with the destruction of an LSP at the lattice site $\mathbf{R}^{\prime}_{s^{\prime}}$, similary to a hopping term.

There are however important differences between the plasmonic Hamiltonian (II) and an electronic tight-binding Hamiltonian. Firstly, plasmons correspond to bosonic excitations, which do not have a finite chemical potential. Secondly, the dipole–dipole interaction responsible for the existence of the collective plasmons is quite different compared to the hopping amplitude in a tight-binding model Grosso and Parravicini (2014). On the one hand, the long-range dipolar interaction scales with $1/|\mathbf{R}_{s}-\mathbf{R}^{\prime}_{s^{\prime}}|^{3}$, whereas the hopping amplitude decreases exponentially with the distance. Thus, the dipolar interaction beyond the first neighbors can have important effects. On the other hand, the dipole–dipole interaction depends on the polarization of the excitations, contrarily to those in tight-binding models. Thirdly, there are additional nonresonant terms [$\propto{b_{s}^{\sigma}}^{\dagger}(\mathbf{R}_{s}){b_{s^{\prime}}^{\sigma^{\prime}}}^{\dagger}(\mathbf{R^{\prime}}_{s^{\prime}})+\mathrm{h.c.}$] in Eq. (II) which do not conserve the number of quasiparticles and play a crucial role for physical quantities depending on the plasmonic eigenstates, e.g., the collective mode damping rates (Brandstetter-Kunc et al., 2016). How the above-mentioned differences may crucially affect the plasmonic band structure is extensively discussed in Sec. III.

The Hamiltonian (1) further describes the coupling of the collective plasmons to vacuum electromagnetic modes in a volume $\mathcal{V}$ described by the Hamiltonian

[TABLE]

where $a_{\mathbf{k}}^{\hat{\lambda}_{\mathbf{k}}}$ (${a_{\mathbf{k}}^{\hat{\lambda}_{\mathbf{k}}}}^{\dagger}$) annihilates (creates) a photon with wave vector $\mathbf{k}$, transverse polarization $\hat{\lambda}_{\mathbf{k}}$ (i.e., $\mathbf{k}\cdot\hat{\lambda}_{\mathbf{k}}=0$), and dispersion $\omega_{\mathbf{k}}=c|\mathbf{k}|$. In the long-wavelength limit $|\mathbf{k}|a\ll 1$, the minimal-coupling Hamiltonian between plasmons and photons in Eq. (1) reads Cohen-Tannoundji et al. (1992)

[TABLE]

with

[TABLE]

and where

[TABLE]

is the vector potential evaluated at the nanoparticle centers. Note that since we consider interparticle separation distances much smaller than the wavelength associated with the LSP resonances, we neglect Umklapp processes in Eqs. (8) and (9). Importantly, within the Coulomb gauge, the first term on the right-hand side of Eq. (9) contains the retardation effects stemming from the finite velocity of light Craig and Thirunamachandran (1984).

In addition to the photonic environment, the collective plasmons are coupled to electron-hole excitations described by the Hamiltonian Weick et al. (2005, 2006)

[TABLE]

where $c^{\phantom{\dagger}}_{\mathbf{R}_{s}i}$ ($c^{\dagger}_{\mathbf{R}_{s}i}$) annihilates (creates) an electron in the nanoparticle located at $\mathbf{R}_{s}$ associated with the one-body state $|\mathbf{R}_{s}i\rangle$ with energy $\epsilon_{\mathbf{R}_{s}i}$ in the self-consistent potential $V$ of that nanoparticle. Assuming $V$ to be a spherical hard-wall potential, the coupling between plasmonic and single electronic degrees of freedom is Brandstetter-Kunc et al. (2015)

[TABLE]

with

[TABLE]

III Quasistatic plasmonic band structure

In this section, we start by focusing on the plasmonic degrees of freedom alone. We first present the diagonalization procedure of the Hamiltonian (II), from which we obtain the quasistatic plasmonic band structure. The latter is then analyzed in detail for both the cases of Bravais (Sec. III.1) and non-Bravais lattices (Sec. III.2). The discussion about how the band structure is influenced by retardation effects in the dipole–dipole interaction is postponed to Sec. IV.

The plasmonic Hamiltonian $H_{\mathrm{pl}}$ given in Eq. (II) is quadratic, and can thus be diagonalized exactly by means of a bosonic Bogoliubov transformation. In the large metasurface limit, going to Fourier space, we then introduce a set of bosonic operators (see Appendix A for details)

[TABLE]

annihilating a collective plasmon with wave vector $\mathbf{q}$ in the band $\tau$ with polarization $\varepsilon$. Notice that, in general, $\varepsilon=\varepsilon_{\tau}(\mathbf{q})$ is a $\mathbf{q}$- and $\tau$-dependent quantity but we drop in the remaining of this paper both indexes for notational simplicity. In Eq. (15), $b_{s}^{\sigma}(\mathbf{q})$ corresponds to the Fourier transform of the bosonic ladder operator (5) [cf. Eq. (61)], while $u_{\tau s}^{\varepsilon\sigma}(\mathbf{q})$ and $v_{\tau s}^{\varepsilon\sigma}(\mathbf{q})$ are complex coefficients which are determined by imposing that $H_{\mathrm{pl}}$ is diagonal in this new basis, i.e.,

[TABLE]

where $\omega_{\tau}^{\varepsilon}(\mathbf{q})$ is the quasistatic collective plasmon dispersion relation. The Bogoliubov operators $\beta_{\tau}^{\varepsilon}(\mathbf{q})$ and ${\beta_{\tau}^{\varepsilon}}^{\dagger}(\mathbf{q})$ act on an eigenstate $|n_{\tau}^{\varepsilon}(\mathbf{q})\rangle$ of the Hamiltonian $H_{\tau}^{\varepsilon}(\mathbf{q})$ corresponding to a collective plasmon in the band $\tau$ with wave vector $\mathbf{q}$ and polarization $\varepsilon$ as $\beta_{\tau}^{\varepsilon}(\mathbf{q})|n_{\tau}^{\varepsilon}(\mathbf{q})\rangle=\sqrt{n_{\tau}^{\varepsilon}(\mathbf{q})}|n_{\tau}^{\varepsilon}(\mathbf{q})-1\rangle$ and ${\beta_{\tau}^{\varepsilon}}^{\dagger}(\mathbf{q})|n_{\tau}^{\varepsilon}(\mathbf{q})\rangle=\sqrt{n_{\tau}^{\varepsilon}(\mathbf{q})+1}|n_{\tau}^{\varepsilon}(\mathbf{q})+1\rangle$, respectively. Here, $n_{\tau}^{\varepsilon}(\mathbf{q})$ is a non-negative integer.

The dispersion relation $\omega_{\tau}^{\varepsilon}(\mathbf{q})$, as well as the coefficients of the Bogoliubov transformation (15), are obtained from the Heisenberg equation of motion [cf. Eq. (66)], which yields the system of equations

[TABLE]

In Eq. (17) the lattice sum

[TABLE]

takes into account the quasistatic dipolar interaction between each pair of nanoparticles composing the metasurface. In the remainder of the paper, the lattice sum (18) is calculated for a specific metasurface numerically, until satisfactory convergence is obtained. In practice, we perform the summation in Eq. (18) up to $\rho_{ss^{\prime}}^{\mathrm{max}}=300d$, which yields a relative error of the order of $10^{-9}$.

The system of equations (17) needs to be satisfied for all integer $s\in[1,\mathcal{S}]$ and for all polarizations $\sigma=x,y,z$, yielding a $6\mathcal{S}\times 6\mathcal{S}$ eigensystem. Due to the structure of the lattice sum (18), such an eigensystem decouples into a block-diagonal matrix composed of a $4\mathcal{S}\times 4\mathcal{S}$ and a $2\mathcal{S}\times 2\mathcal{S}$ block, corresponding to the in-plane (IP, $\sigma=x,y$) and out-of-plane (OP, $\sigma=z$) polarized modes, respectively. Each block then yields a secular equation of order $2\mathcal{S}$ and $\mathcal{S}$ in $[\omega_{\tau}^{\varepsilon}(\mathbf{q})]^{2}$, respectively, which then gives access to the quasistatic band structure.

Additional insight about the nature of the quasistatic collective plasmons can be obtained from their corresponding eigenstates, from which we can deduce the polarization $\varepsilon$ of the collective modes. Introducing the vector $\mathbf{u}_{\tau}^{\varepsilon}(\mathbf{q})=(\sum_{s}u_{\tau s}^{\varepsilon x}(\mathbf{q}),\sum_{s}u_{\tau s}^{\varepsilon y}(\mathbf{q}),\sum_{s}u_{\tau s}^{\varepsilon z}(\mathbf{q}))$, we define the polarization angle as Lamowski et al. (2018)

[TABLE]

With such a definition, purely longitudinal (transverse) collective plasmons correspond to $\phi_{\tau}^{\varepsilon}(\mathbf{q})=0$ ($\pi/2$).

III.1 Bravais lattices

In the case of an arbitrary Bravais lattice (i.e., $\mathcal{S}=1$), the sublattice indexes $s$ and $s^{\prime}$, as well as the band index $\tau$ are irrelevant and can then be dropped from the system of equations (17). Such a system can be fully analytically solved, yielding for the OP plasmonic modes polarized in the $\varepsilon=z$ direction the dispersion relation

[TABLE]

and

[TABLE]

for the IP modes with polarizations $\varepsilon=\varepsilon_{\parallel,\pm}$. Explicit expressions for the corresponding Bogoliubov coefficients can be found in Appendix A.

We notice that the quasistatic dispersion relations (20) and (21) have been obtained by Zhen et al. Zhen et al. (2008) using a classical point-dipole model [cf. their Eq. (3.3)]. It is not surprising that, within the quasistatic approximation, our quantum approach and the classical one of Ref. Zhen et al. (2008) (which however is limited to simple Bravais lattices) give the same results since the Hamiltonian at hand [cf. Eq. (2)] is quadratic.

As an example of application of our general method for obtaining the quasistatic band structure of plasmonic modes in generic Bravais lattices, we consider in the following the simple square lattice sketched in Fig. 2(a), whose corresponding first Brillouin zone (1BZ) is depicted in Fig. 2(b). The plasmonic dispersion relation (20) for the square lattice is plotted in Fig. 2(c) for the OP polarization as a solid line. For comparison, we also show (dashed line) the plasmonic band structure considering only dipolar interactions between nearest neighbors (nn) in the lattice, for which the lattice sum (18) reduces to $f_{\mathrm{nn}}^{zz}(\mathbf{q})=2[\cos(q_{x}d)+\cos(q_{y}d)]$. As can be seen from Fig. 2(c), the nearest-neighbor approximation qualitatively reproduces the full band structure in most of the 1BZ, except for wave numbers close to the $\Gamma$ point. There, the long-range nature of the quasistatic dipolar interaction leads to a pronounced cusp of the dispersion relation. Such a behavior, which is purely an artifact of the quasistatic approximation and which points to the importance of taking into account retardation effects, as we shall do in Sec. IV, can be understood by a mean-field treatment of the dipolar interactions beyond those between nearest neighbors (see Appendix B).

We now turn to the discussion of the plasmonic modes polarized within the plane formed by the square lattice. The band structure (21) is plotted in Fig. 2(d) as solid lines. The color code corresponds to the polarization angle defined in Eq. (19). While for the high-symmetry axes $\Gamma M$ or $\Gamma X$, the IP collective plasmons are purely longitudinal or transverse, for less-symmetric axes such as the $MX$ direction in the 1BZ, such modes can be of a mixed type. For comparison, we further plot the dispersion relation taking into account nearest-neighbor couplings only, for which the lattice sums in Eq. (21) are replaced by $f^{xx}_{\mathrm{nn}}(\mathbf{q})=-4\cos{(q_{x}d)}+2\cos{(q_{y}d)}$, $f^{yy}_{\mathrm{nn}}(\mathbf{q})=2\cos{(q_{x}d)}-4\cos{(q_{y}d)}$, and $f^{xy}_{\mathrm{nn}}(\mathbf{q})=0$. In contrast to the OP modes [Fig. 2(c)], the long-range nature of the dipolar interaction has a more pronounced effect on the plasmonic band structure for IP polarized modes. For instance, the dipolar interaction lifts the degeneracy of the dispersion induced by the symmetry of the square lattice within the nearest-neighbor approximation along the $\Gamma M$ direction of the 1BZ. In addition, the long-range dipolar interaction leads to a cusp of the upper plasmonic band at the $\Gamma$ point, while the lower band does not show such a singularity in the derivative of $\omega^{\varepsilon_{\parallel,\pm}}(\mathbf{q})$. This behavior can be explained along the lines of the mean-field approximation discussed in Appendix B.

III.2 Bravais lattices with a basis

Our general method for obtaining quasistatic plasmonic band structures further applies to arbitrary Bravais lattices with a basis. In the following, we start by considering Bravais lattices with a basis of two.

III.2.1 Bipartite lattices

In the case of a bipartite lattice ($\mathcal{S}=2$), the $4\times 4$ matrix resulting from the system of equations (17) for the OP polarization $\sigma=z$ can be straightforwardly solved, yielding the two bands with dispersion relations

[TABLE]

The corresponding Bogoliubov coefficients are given in Eqs. (71) and (72). For IP polarization ($\sigma=x,y$), the $8\times 8$ eigenvalue problem can in principle be solved analytically, but provides cumbersome expressions. For practical purposes we therefore solve for the eigenproblem numerically.

To illustrate our method in the special case of bipartite lattices of near-field coupled metallic nanoparticles, we consider the celebrated honeycomb lattice sketched in Fig. 3(a). The corresponding 1BZ is shown in Fig. 3(b). Such a metasurface has been predicted Weick et al. (2013); Sturges et al. (2015) to exhibit Dirac-like collective plasmonic modes at the $K$ and $K^{\prime}$ points of the 1BZ, with appealing topological properties such as a nontrivial Berry phase (and its related topologically protected edge states Delplace et al. (2011)) or the absence of backscattering. Importantly, the results put forward in Refs. Weick et al. (2013); Sturges et al. (2015) rely on short-range dipolar couplings between nearest neighbors alone. Moreover, Refs. Weick et al. (2013); Sturges et al. (2015) considered the case of orientated dipoles, relevant, e.g., for elongated metallic rods, while we consider here the case of spherical nanoparticles.

The two plasmonic bands (22) for $\sigma=z$ are plotted in Fig. 3(c) as solid lines. For comparison, we also show by dashed lines the plasmonic band structure with dipolar interactions between nearest neighbors only Weick et al. (2013), given by Eq. (22) and replacing $f_{11}^{zz}(\mathbf{q})$ and $f_{12}^{zz}(\mathbf{q})$ by $f_{\mathrm{nn},11}^{zz}(\mathbf{q})=0$ and $f_{\mathrm{nn},12}^{zz}(\mathbf{q})=\sum_{j=1}^{3}\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{e}_{j}}$, respectively. Here, the vectors $\mathbf{e}_{1}=\mathbf{d}_{2}-\mathbf{t}_{2}$, $\mathbf{e}_{2}=\mathbf{d}_{2}$, and $\mathbf{e}_{3}=\mathbf{d}_{2}-\mathbf{t}_{1}$ connect a lattice site belonging to the red sublattice in Fig. 3(a) to its three (blue) nearest neighbors. We observe in Fig. 3(c) the presence of a cusp for the upper ($\tau=+$) band when all quasistatic interactions are taken into account, while no cusp appears for the lower ($\tau=-$) band. Notice that the upper (lower) band corresponds to bright (dark) modes, where the two dipolar LSPs are in-phase (out-of-phase) within each unit cell.

As can be seen from Fig. 3(c), the presence of a Dirac point at the $K$ point of the 1BZ located at $\mathbf{K}=\frac{4\pi}{3\sqrt{3}d}(1,0)$ is not ruled out by long-range interactions. Indeed, in the vicinity of the $K$ point, expanding the lattice sums to linear order in $|\mathbf{k}|$, where $\mathbf{q}=\mathbf{K}+\mathbf{k}$ with $|\mathbf{k}|\ll|\mathbf{K}|$, yields $f_{11}^{zz}(\mathbf{q})\simeq f_{11}^{zz}(\mathbf{K})\simeq-0.449$ and $f_{12}^{zz}(\mathbf{q})\simeq-1.16(k_{x}+\mathrm{i}k_{y})d$. Therefore, in the weak-coupling regime $\Omega\ll\omega_{0}$, the dispersion (22) is conical and reads

[TABLE]

with the group velocity $v^{z}=1.16\Omega d$. Comparing the dispersion above with the nearest-neighbor result Weick et al. (2013) $\omega_{\mathrm{nn},\tau}^{z}(\mathbf{k})=\omega_{0}+\tau v_{\mathrm{nn}}^{z}|\mathbf{k}|$ with $v_{\mathrm{nn}}^{z}=3\Omega d/2$, we see that the intrasublattice coupling $f_{11}^{zz}$ leads to an inconsequential redshift of the Dirac point frequency, while the intersublattice coupling $f_{12}^{zz}$ renormalizes the precise value of the group velocity.

Since the Bogoliubov coefficients (72) are negligible as compared to the coefficients (71) close to the $K$ point, we can safely disregard the former, which amounts to performing the rotating wave approximation (RWA) SM . Within this limit, the associated effective Hamiltonian reads in terms of the spinor operator $\hat{\Psi}_{\mathbf{k}}=(b_{1}^{z}(\mathbf{k}),b_{2}^{z}(\mathbf{k}))$ as $H_{\mathrm{pl}}^{\mathrm{eff}}=\sum_{\mathbf{k}}\hat{\Psi}_{\mathbf{k}}^{\dagger}\mathcal{H}^{\mathrm{eff}}_{\mathbf{k}}\hat{\Psi}_{\mathbf{k}}$, with the massless Dirac Hamiltonian

[TABLE]

where $\mathbb{1}_{n}$ is the $n\times n$ identity matrix and $\boldsymbol{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ is the vector of Pauli matrices

[TABLE]

acting on the sublattice pseudospin $1/2$. The corresponding plasmonic eigenstates [with eigenfrequencies given by Eq. (23)] read

[TABLE]

where $\theta(\mathbf{k})$ is the angle between the wave vector $\mathbf{k}$ and the $x$ axis, such that $\tan{\theta(\mathbf{k})}=k_{y}/k_{x}$. The spinors (26) are also eigenstates of the chirality (or helicity) operator Allain and Fuchs (2011), $\mathcal{C}=\boldsymbol{\sigma}\cdot\hat{k}$, which corresponds to the projection of the pseudospin operator onto the wave vector direction, with eigenvalues $C=-\tau$. Since $[\mathcal{H}^{\mathrm{eff}}_{\mathbf{k}},\mathcal{C}]=0$, chirality is conserved. As a consequence, collective plasmonic modes close to the $K$ point show similar effects as electrons in graphene Weick et al. (2013), such as the absence of backscattering off smooth inhomogeneities (since $\langle\psi_{\tau}^{z}(-\mathbf{k})|\psi_{\tau}^{z}(\mathbf{k})\rangle=0$ Allain and Fuchs (2011)), the related Klein tunneling phenomenon, or a Berry flux of $-\pi$. Note that the spinor eigenstates (26) do not depend on the details of the Hamiltonian (24), such as, e.g., the precise value of the group velocity $v^{z}$, but only on its matrix structure. Henceforth, the properties mentioned above do not depend on the level of approximation one may use (nearest-neighbor coupling vs long-range quasistatic interactions). In addition, we will demonstrate in Sec. IV that the retarded part of the dipolar coupling does not modify the above physics. Therefore, the long-range character of the dipole–dipole interaction does not rule out the massless Dirac nature of the plasmonic quasiparticles in the vicinity of the $K$ point, and the nearest-neighbor approximation is sufficient in catching the relevant physics. The same conclusion applies to the inequivalent Dirac point located at $K^{\prime}$.

We now turn to the description of the IP polarized plasmonic modes. We show in Fig. 3(d) the plasmonic band structure obtained numerically for $\sigma=x,y$ as solid lines. The color code corresponding to the polarization angle (19) reveals that two bands correspond to purely transverse plasmons, and two other bands to purely longitudinal plasmons along the high-symmetry axes $\Gamma K$ and $\Gamma M$. We also plot for comparison the dispersion relations with nearest-neighbor coupling only. In the latter case, the intrasublattice sums $f^{\sigma\sigma^{\prime}}_{\mathrm{nn},11}(\mathbf{q})$ vanish, so that the $8\times 8$ system given by Eq. (17) simplifies and can be straightforwardly solved analytically. The resulting four dispersion relations read

[TABLE]

with

[TABLE]

and

[TABLE]

In the expressions above, the intersublattice sums within the nearest-neighbor approximation read $f_{\mathrm{nn},12}^{xx}(\mathbf{q})=\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{e}_{1}}-\frac{5}{4}(\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{e}_{2}}+\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{e}_{3}})$, $f_{\mathrm{nn},12}^{yy}(\mathbf{q})=-2\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{e}_{1}}+\frac{1}{4}(\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{e}_{2}}+\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{e}_{3}})$, and $f_{\mathrm{nn},12}^{xy}(\mathbf{q})=\frac{3\sqrt{3}}{4}(\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{e}_{3}}-\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{e}_{2}})$. As can be seen from Fig. 3(d), the nearest-neighbor approximation is sufficient to qualitatively describe the plasmonic dispersion relations, apart from the second less energetic band for wave numbers close to the $\Gamma$ point, where a cusp appears in the full quasistatic band structure. Importantly, we note the presence of two inequivalent conical intersections (where the band degeneracy point occurs at the frequency $\omega_{0}$) at the $K$ and $K^{\prime}$ points of the 1BZ for IP polarized modes. In the vicinity of the $K$ point, we find for the second and third bands $\omega_{\mathrm{nn},\tau}^{\varepsilon_{\parallel}}(\mathbf{k})\simeq\omega_{0}\pm v_{\mathrm{nn}}^{\varepsilon_{\parallel}}|\mathbf{k}|$, with $v_{\mathrm{nn}}^{\varepsilon_{\parallel}}=9\Omega d/4$. The presence of conical intersections for IP polarized modes has been reported by Han et al. Han et al. (2009) using a numerical solution to Maxwell’s equations. Our transparent method allows us to analytically describe such a complex band structure hosting Dirac-like bosonic modes.

III.2.2 Tripartite lattices

For tripartite lattices ($\mathcal{S}=3$), the $6\times 6$ system obtained from Eq. (17) for $\sigma=z$ can still be solved analytically, yielding the three plasmonic bands

[TABLE]

where

[TABLE]

with

[TABLE]

and

[TABLE]

For $\sigma=x,y$, the $12\times 12$ eigensystem is solved numerically.

As an illustration, we consider in the remainder of this section the Lieb lattice sketched in Fig. 4(a). Its corresponding 1BZ is shown in Fig. 4(b). Such a lattice, together with the kagome one SM , is known to display a flat (i.e., nondispersive) band together with conical dispersions in tight-binding models with nearest-neighbor interactions (see, e.g., Refs. Bergman et al. (2008); Guo and Franz (2009); Weeks and Franz (2010); Niţă et al. (2013); Palumbo and Meichanetzidis (2015)). It is therefore of interest to study if these features survive in the case of plasmonic metasurfaces, where the long-range nature of the dipolar interaction might qualitatively change the band structure.

We plot in Fig. 4(c) the three dispersion relations (30) for the Lieb lattice and for OP polarization as solid lines. For comparison, we also show by dashed lines the plasmonic band structure within the nearest-neighbor approximation, for which we have $f_{\mathrm{nn},11}^{zz}(\mathbf{q})=f_{\mathrm{nn},22}^{zz}(\mathbf{q})=f_{\mathrm{nn},33}^{zz}(\mathbf{q})=f_{\mathrm{nn},23}^{zz}(\mathbf{q})=0$, $f_{\mathrm{nn},12}^{zz}(\mathbf{q})=2\cos{(q_{x}d)}$, and $f_{\mathrm{nn},13}^{zz}(\mathbf{q})=2\cos{(q_{y}d)}$. Within this approximation, the resulting band structure is given by $\omega_{\mathrm{nn},\tau}^{z}(\mathbf{q})=\omega_{0}[1+\tau\frac{4\Omega}{\omega_{0}}\sqrt{\cos^{2}{(q_{x}d)}+\cos^{2}{(q_{y}d)}}]^{1/2}$ and is characterized by the presence of a flat band at $\omega_{0}$ ($\tau=0$) and two dispersive bands ($\tau=\pm 1$). For wave vectors $\mathbf{k}$ in the vicinity of the $M$ point of the 1BZ located at $\mathbf{M}=\frac{\pi}{2d}(1,1)$, these latter form a conical dispersion $\omega_{\mathrm{nn},\pm 1}^{z}(\mathbf{k})\simeq\omega_{0}\pm v_{\mathrm{nn}}^{z}|\mathbf{k}|$, with group velocity $v^{z}_{\mathrm{nn}}=2\Omega d$.

As shown in Fig. 4(c), the long-range dipolar interactions affect differently the topological features of the plasmonic band structure described above (compare solid and dashed lines). Indeed, while the presence of a conical dispersion in the vicinity of the $M$ point is not ruled out by long-range interactions (as is the case for the honeycomb array, see Sec. III.2.1), the band which is flat in the whole 1BZ within the nearest-neighbor approximation becomes dispersive, and is only locally flat close to the $M$ point. One can understand these features by expanding the band structure (30) in the vicinity of the $M$ point. With $\mathbf{q}=\mathbf{M}+\mathbf{k}$, to linear order in $|\mathbf{k}|\ll|\mathbf{M}|$ we find $f_{ss}^{zz}(\mathbf{q})\simeq f_{ss}^{zz}(\mathbf{M})\simeq-0.331$, $f_{12}^{zz}(\mathbf{q})\simeq-1.65k_{x}d$, $f_{13}^{zz}(\mathbf{q})\simeq-1.65k_{y}d$, and $f_{23}^{zz}(\mathbf{q})\simeq 0$, resulting in

[TABLE]

i.e., a band $\tau=0$ which is locally flat, and two conical dispersions $\tau=\pm 1$ with renormalized group velocity $v^{z}=1.65\Omega d$.

Interestingly, the plasmonic Hamiltonian within the RWA for modes polarized in the $z$ direction takes the form $H_{\mathrm{pl}}^{\mathrm{eff}}=\sum_{\mathbf{k}}\hat{\Psi}_{\mathbf{k}}^{\dagger}\mathcal{H}^{\mathrm{eff}}_{\mathbf{k}}\hat{\Psi}_{\mathbf{k}}$, with $\hat{\Psi}_{\mathbf{k}}=(b_{1}^{z}(\mathbf{k}),b_{2}^{z}(\mathbf{k}),b_{3}^{z}(\mathbf{k}))$ and where

[TABLE]

with $\mathbf{S}=(S_{x},S_{y},S_{z})$. Here, the pseudospin-$1$ matrices (corresponding to the three sublattices of the Lieb lattice) are defined as

[TABLE]

and fulfill the angular momentum algebra $[S_{i},S_{j}]=\mathrm{i}\epsilon_{ijk}S_{k}$, with $\epsilon_{ijk}$ the Levi-Civita symbol. The matrices $S_{x}$, $S_{y}$, and $S_{z}$ (which correspond, respectively, to the Gell-Mann matrices $\lambda_{1}$, $\lambda_{4}$, and $\lambda_{7}$ Gell-Mann (1962); Arfken et al. (2013)) therefore correspond to a three-dimensional representation of the special unitary group $SU(2)$. However, contrary to the Pauli matrices, they do not form a Clifford algebra (i.e., $\{S_{i},S_{j}\}\neq 2\delta_{ij}\mathbb{1}_{3}$), so that Eq. (35) does not correspond to a massless Dirac Hamiltonian Goldman et al. (2011), despite presenting a conical spectrum. The eigenspinors corresponding to the Hamiltonian (35) with eigenfrequencies (34), and characterized by a vanishing Berry phase, read

[TABLE]

for $\tau=0,\pm 1$ and are nevertheless eigenstates of the helicity operator $\mathbf{S}\cdot\hat{k}$ with eigenvalues $-\tau$, so that backscattering is suppressed for the two dispersive bands with $\tau=\pm 1$ ($\langle\psi_{\pm 1}^{z}(-\mathbf{k})|\psi_{\pm 1}^{z}(\mathbf{k})\rangle=0$). The conclusion above is not modified when we take into account the retarded part of the dipolar interaction (see Sec. IV).

We show in Fig. 4(d) the plasmonic dispersion relations calculated numerically for $\sigma=x,y$ as solid lines. We further plot, for comparison, the numerical results with nearest-neighbor couplings only. We observe from the figure that for IP polarization, the long-range dipolar interactions completely reconstruct the plasmonic band structure. Notably, the topological features (flat bands and conical dispersions) occurring in the nearest-neighbors coupling approximation are not preserved when long-range interactions are included.

IV Perturbative treatment of the retardation effects onto the plasmonic band structure

We now consider the effects of the photonic environment alone [encapsulated in the Hamiltonian (8)] onto the collective plasmonic excitations supported by our generic metasurface of ordered metallic nanoparticles. The plasmon–photon interaction [cf. Eq. (9)] leads to two effects crucially affecting the quasistatic band structure discussed above: (i) the photon-induced frequency shifts which we unveil in the present section and resulting from retardation effects in the dipole–dipole interaction renormalize the quasistatic plasmonic dispersion; (ii) moreover, the spontaneous decay of plasmons into free photons leads to a finite radiative lifetime of the collective excitations, and consequently to a broadening of the plasmon lines (see Sec. V).

We start our analysis of the effects of the photonic environment onto the quasistatic plasmonic dispersion by considering the radiative frequency shifts induced by the light–matter coupling. Along the lines of Refs. Downing et al. (2017, 2018); Downing and Weick (2018), we treat the plasmon–photon coupling Hamiltonian (9) to second order in standard nondegenerate perturbation theory, yielding for the collective plasmon energy levels the result $E_{n_{\tau}^{\varepsilon}(\mathbf{q})}=n_{\tau}^{\varepsilon}(\mathbf{q})\hbar\omega_{\tau}^{\varepsilon}(\mathbf{q})+E_{n_{\tau}^{\varepsilon}(\mathbf{q})}^{(1)}+E_{n_{\tau}^{\varepsilon}(\mathbf{q})}^{(2)}$. The first term in the right-hand side of the equation above corresponds to the energy levels of the unperturbed Hamiltonian $H_{\tau}^{\varepsilon}(\mathbf{q})$ [cf. Eq. (16)]. The first-order (in $H_{\mathrm{pl\textrm{-}ph}}$) correction stems from the diamagnetic term in Eq. (9) (proportional to the vector potential squared) and reads $E_{n_{\tau}^{\varepsilon}(\mathbf{q})}^{(1)}=2\pi\mathcal{SN}\hbar\omega_{0}^{2}a^{3}\mathcal{V}^{-1}\sum_{\mathbf{k}}\omega_{\mathbf{k}}^{-1}$. Since the latter expression does not depend on the quantum number $n_{\tau}^{\varepsilon}(\mathbf{q})$, it does not participate in the renormalization of the collective plasmon frequency, and merely represents an irrelevant global energy shift. The second-order correction to $n_{\tau}^{\varepsilon}(\mathbf{q})\hbar\omega_{\tau}^{\varepsilon}(\mathbf{q})$ arising from the first (paramagnetic) term on the right-hand side of Eq. (9) reads

[TABLE]

where the summation over $\mathbf{k}$ excludes the singular term for which $\omega_{\mathbf{k}}=\omega_{\tau}^{\varepsilon}(\mathbf{q})$. In Eq. (IV), we have defined the array factor

[TABLE]

and

[TABLE]

In the large-metasurface limit (number of unit cells $\mathcal{N}\gg 1$), the modulus squared of the array factor above entering Eq. (IV) takes the simpler form

[TABLE]

so that the radiative frequency shift, defined through the relation $\delta_{\tau}^{\varepsilon}(\mathbf{q})=[E_{n_{\tau}^{\varepsilon}(\mathbf{q})+1}-E_{n_{\tau}^{\varepsilon}(\mathbf{q})}]/\hbar-\omega_{\tau}^{\varepsilon}(\mathbf{q})$, is given by

[TABLE]

and scales with the nanoparticle sizes as $a^{3}$, to leading order in the coupling constant $\Omega/\omega_{0}$ [cf. Eq. (7)].

The analytical result above is valid for both IP and OP plasmon polarizations, and depends on the quasistatic plasmon dispersion $\omega_{\tau}^{\varepsilon}(\mathbf{q})$ and the Bogoliubov coefficients [cf. Eq. (40)]. Equation (IV) can be expicitly evaluated SM , and we find

[TABLE]

and

[TABLE]

for the radiative shifts of the OP and IP plasmonic modes, respectively. Notably, the calculation leading to the above results does not need the introduction of an ultraviolet frequency cutoff for the photonic degrees of freedom [which prevent to take into account photons with a wavelength smaller than the nanoparticle size, cf. the dipolar approximation in Eq. (9)], as is the case for nanoparticle dimers Downing et al. (2017) and linear chains Downing et al. (2018); Downing and Weick (2018). We have checked that the inclusion of such a cutoff, of the order of $k_{\mathrm{c}}\simeq 1/a$, does not significantly modify the results (IV) and (IV). Notice also that Eqs. (IV) and (IV) are not periodic in reciprocal space since we consider only the interaction between the collective plasmons and photons for which the associated light cone belongs to the 1BZ. Finally, we note from Eqs. (IV) and (IV) that $\delta_{\tau}^{\varepsilon}(\mathbf{0})=0$, so that the dispersion relation at the $\Gamma$ point is not renormalized by the light–matter interaction.

We plot in Figs. 5(a)–5(c) the plasmonic band structure for OP-polarized modes $\omega_{\tau}^{z}(\mathbf{q})+\delta_{\tau}^{z}(\mathbf{q})$ including Eq. (IV) for (a) the square, (b) the honeycomb, and (c) the Lieb lattices along the high symmetry lines of their respective 1BZ. As shown in Fig. 5(a) for the square lattice, the radiative frequency shift $\delta_{\tau}^{z}(\mathbf{q})$ induces an important renormalization of the quasistatic band structure, which is of purely transverse nature [see Fig. 2(c)]. A similar effect is also observed for other simple Bravais lattices, such as the rectangular and hexagonal ones SM . Comparing the retarded dispersion relation (red line) to the quasistatic one (thin gray line), we first observe that the retarded one diverges at the intersection of the quasistatic band structure with the light cone (dashed lines). A polaritonic singularity arising from the strong light–matter interaction was already reported in the literature for the transverse modes in 1D plasmonic chains by means of involved numerical calculations based on the fully retarded solutions to Maxwell’s equations Weber and Ford (2004); Citrin (2004); Simovski et al. (2005); Citrin (2006); Koenderink and Polman (2006); Markel and Sarychev (2007); Fung and Chan (2007); Petrov (2015). Recently, studies Downing et al. (2018); Downing and Weick (2018) using the open quantum system approach which we use in this work have shown a similar divergence occurring in the dispersion relations of transverse plasmonic modes in chains. Importantly, the results presented in Ref. Downing et al. (2018) show a good qualititative agreement with numerical electromagnetic calculations for regular nanoparticle chains. We here point out that the singularities observed in Fig. 5 arise from calculations based on a second-order perturbative treatment of the light–matter coupling. Consequently, important variations from the LSP resonance frequency $\omega_{0}$ should be treated carefully. Notice also that the renormalized dispersion relation diverges as the inverse of a square root [see Eq. (IV)] instead of logarithmically as is the case for 1D arrays Downing et al. (2018); Downing and Weick (2018).

As can be seen from Fig. 5(a), taking into account the radiative shift implies that the cusp appearing at the $\Gamma$ point of the 1BZ within the quasistatic approximation disappears (compare the red and gray lines in the figure). In the vicinity of the $\Gamma$ point ($|\mathbf{q}|d\ll 1$), Eq. (IV) applied to the square lattice reduces (considering $\Omega/\omega_{0}\ll 1)$ to $\delta^{z}(\mathbf{q})\simeq 2\pi\Omega|\mathbf{q}|d+\mathcal{O}(|\mathbf{q}|d)^{3}$. Such a linear $|\mathbf{q}|$-dependence cancels out exactly the one of $\omega^{z}(\mathbf{q})$ in this regime of parameters [see Eq. (75)]. We hence find $\omega^{z}(\mathbf{q})+\delta^{z}(\mathbf{q})\simeq\omega_{0}+\Omega[4+\sqrt{2}\pi+(|\mathbf{q}|d)^{2}(\pi/\sqrt{2}-1)]$ close to the $\Gamma$ point, leading to a quadratic dependence of the dispersion relation. The results presented in Fig. 5(a) thus demonstrate that considering retardation effects is crucial for the study of the collective plasmonic modes in metasurfaces of near-field coupled nanoparticles, since the dispersion relations are qualitatively affected by the interactions with free photons. Such renormalization effects are not that prominent in 1D plasmonic systems, where, apart from the pronounced renormalization of the dispersion relation at the crossing of the light cone, the band structure is qualitatively unaffected by retardation effects Downing et al. (2018); Downing and Weick (2018); Weber and Ford (2004); Citrin (2004); Simovski et al. (2005); Citrin (2006); Koenderink and Polman (2006); Markel and Sarychev (2007); Fung and Chan (2007); Petrov (2015).

We show in Fig. 5(b) the band structure including retardation effects (colored lines) of the OP plasmonic modes for the honeycomb lattice. For comparison, we also reproduce as gray solid lines the transverse-polarized quasistatic band structure shown in Fig. 3(c). The upper band ($\tau=+1$) shows a similar profile as that in Fig. 5(a). It displays a divergence at the intersection between the light cone (dashed lines) and the quasistatic dispersion relation. In addition, it does not present a cusp at the $\Gamma$ point. In contrast, the lower band ($\tau=-1$) does not experience a noticeable (on the scale of the figure) renormalization induced by the light–matter coupling. Such low-energy plasmonic modes are thus coined “dark” modes since they only weakly couple to light. They correspond to LSPs within each sublattices that are out of phase. Conversely, plasmonic modes which show a significant coupling to light are called “bright” modes. Such modes correspond to LSPs within each sublattices that are in phase. Notably, the different nature of the bright and dark modes has also important consequences on their respective radiative lifetimes, as we will discuss in Sec. V.

Importantly, the Dirac cone exhibited by the quasistatic band structure remains unaffected by the light–matter coupling since the Dirac point lays outside of the light cone within the regime of parameters which we consider in this work. In Fig. 5(b), we nevertheless observe a slight mismatch between the two bands in the vicinity of the $K$ point. Such an artifact stems from the fact that we only consider the light cone belonging to the 1BZ in the evaluation of Eq. (IV). Full polaritonic numerical calculations Mann et al. (2018), where the light–matter interaction is taken into account exactly (and not perturbatively), show however that Dirac cones are unaffected by retardation effects.

We show in Fig. 5(c) the dispersion relation including retardation effects of the OP plasmonic modes for the Lieb lattice. Along the lines of the above discussion on the honeycomb lattice [cf. Fig. 5(b)], the plasmonic dispersion relation of the Lieb lattice present bright and dark modes. The most energetic band (red line) corresponds to bright transverse collective modes and thus displays a singularity at the intersection between the quasistatic dispersion (gray solid lines) and the light cone (dashed lines) as well as the absence of a cusp at the $\Gamma$ point. Conversely, the two low-energy bands correspond to dark modes.

We now focus on the effects of retardation in the light–matter interaction onto the collective plasmonic modes polarized within the plane of the metasurface [cf. Eq. (IV)]. We display in Fig. 5(d) the dispersion relation of the IP modes in the square lattice (see red and orange lines). The quasistatic band structure from Fig. 2(d) is reproduced here by gray lines for comparison. Three important features appear from the comparison of these two results. The low-energy band (cf. red line), which corresponds essentially to transverse modes, present a singularity at the crossing of the quasistatic dispersion relation with the light cone. Conversely, the high-energy band (cf. orange line), which corresponds essentially to longitudinal modes, does not present such a singularity at the crossing, as is the case for plasmonic chains Weber and Ford (2004); Citrin (2004); Simovski et al. (2005); Citrin (2006); Koenderink and Polman (2006); Markel and Sarychev (2007); Fung and Chan (2007); Petrov (2015); Downing et al. (2018); Downing and Weick (2018). In addition, the cusp that present the quasistatic high-energy band in the vicinity of the $\Gamma$ point is washed away by retardation effects. Similar conclusions can be drawn for the honeycomb and Lieb lattices, see Figs. 5(e) and 5(f), respectively, as well as for the kagome lattice SM . In addition, in the case of Bravais lattices with a basis, some of the IP bands are only weakly modified by retardation effects [see the cyan lines in Figs. 5(e) and 5(f)] and thus correspond to dark modes.

In Appendix C, we compare our perturbative approach to numerical calculations based on the solution to the fully retarded Maxwell’s equations performed in Refs. Han et al. (2009); Zhen et al. (2008), and find a rather satisfying agreement. Contrarily to three-dimensional plasmonic systems Lamowski et al. (2018) or metamaterials embedded in a photonic cavity Downing et al. (2019); Mann et al. (2018, ), where matter and photonic degrees of freedom are strongly coupled and which require an exact diagonalization of the light-matter coupling Hamiltonian, here the collective plasmons supported by the metasurface couple to a continuum of photonic modes as $\mathcal{V}\rightarrow\infty$ in Eq. (11). Therefore, the coupling is not as strong as in the previously mentioned systems, and leading-order perturbation theory is able to catch the main features of the polaritonic band structure.

V Radiative damping

We now concentrate on the evaluation of the radiative decay rate of the collective plasmons. To this end, we treat the plasmon–photon coupling Hamiltonian (9) as a weak perturbation to the plasmonic subsystem. In such a regime, the radiative decay rate of the plasmonic eigenmode $|1_{\tau}^{\varepsilon}(\mathbf{q})\rangle$ with band index $\tau$, polarization $\varepsilon$, and wave vector $\mathbf{q}$ is given by the Fermi golden rule expression

[TABLE]

where $F_{\mathbf{k},\mathbf{q}}^{-}$ and $P_{\tau s}^{\varepsilon\sigma}(\mathbf{q})$ are defined in Eqs. (39) and (40), respectively. In the large metasurface limit ($\mathcal{N}\gg 1$), using Eq. (41) yields the result

[TABLE]

which is valid for arbitrary polarizations. We can then show that SM

[TABLE]

for the OP modes, and

[TABLE]

for the IP modes.

We plot in Fig. 6 the radiative damping rate (V) for OP plasmonic modes for the square [Fig. 6(a)], the honeycomb [Fig. 6(b)], and the Lieb lattices [Fig. 6(c)] along the high symmetry lines originating from the $\Gamma$ point of their respective 1BZ. Note that we do not show in the figure the results along the entire 1BZ since the damping rate vanishes for wave vectors outside of the light cone, as can be easily inferred from the Heaviside step function in Eq. (V). In Fig. 6 the displayed results are scaled by the radiative decay rate of a single isolated nanoparticle,

[TABLE]

For the square lattice [Fig. 6(a)], the OP plasmonic modes present a highly superradiant profile ($\gamma^{z}(\mathbf{q})\gg\gamma_{0}$) inside of the light cone for wave vectors not too close to the $\Gamma$ point, while displaying a vanishing rate outside of the light cone. From the figure we observe that the radiative decay rate $\gamma^{z}(\mathbf{q})$ increases rapidly as $\mathbf{q}$ moves away from the center of the 1BZ and diverges at wave vectors corresponding to the intersection of the quasistatic band structure with the light cone. Such singularities are related to those observed in the fully retarded dispersion relation [cf. Fig. 5(a)]. The superradiant behavior of the radiative damping rate observed in Fig. 6(a) is reminiscent to the one reported for 1D plasmonic chains Weber and Ford (2004); Citrin (2004); Markel and Sarychev (2007); Citrin (2006); Koenderink and Polman (2006); Petrov (2015); Brandstetter-Kunc et al. (2016); Downing and Weick (2017); Downing et al. (2018); Downing and Weick (2018). Such a behavior is nevertheless much more prominent in 2D metamaterials, due to the enhanced constructive interferences between the dipolar electric fields produced by each LSP.

We plot in Figs. 6(b) and 6(c) the radiative decay rates of the OP plasmonic modes for the honeycomb and Lieb lattices, respectively. In both cases, one of the plasmonic band shows a similar profile as that in Fig. 6(a). Indeed, in both Figs. 6(b) and 6(c), the red lines display a superradiant behavior which diverges at the intersection of the quasistatic band structure and the light cone. Conversely, the blue and cyan lines show the existence of subradiant modes (i.e., $\gamma_{\tau}^{\varepsilon}(\mathbf{q})\ll\gamma_{0}$) corresponding to dark, nonradiative modes, for which retardation effects on the plasmonic band structure are essentially negligible [cf. blue and cyan lines in Figs. 5(b) and 5(c)].

We show in Fig. 6(d) the radiative damping rate (V) for IP polarized modes in the square lattice. In the figure, the red line corresponds to the lower transverse plasmonic band in Fig. 5(d) [see also Fig. 2(d)] and presents singularities coinciding with those in Fig. 6(a). Conversely, the orange line in Fig. 6(d) which corresponds to the upper longitudinal band [see Fig. 5(d)] displays an opposite trend, as the radiative decay rate decreases for wave vectors moving away from the center of the 1BZ. This is reminiscent to what has been previously reported for longitudinal plasmonic modes in 1D chains Weber and Ford (2004); Citrin (2004); Markel and Sarychev (2007); Citrin (2006); Koenderink and Polman (2006); Petrov (2015); Brandstetter-Kunc et al. (2016); Downing and Weick (2017); Downing et al. (2018); Downing and Weick (2018). We draw similar conclusions for the honeycomb and Lieb lattices, see Figs. 6(e) and 6(f), respectively. Additionally, some of the bands [cyan lines in Figs. 6(e) and 6(f)] display an almost vanishing radiative decay rate as they correspond to dark, out-of-phase modes.

VI Effects of the electronic environment on the collective plasmon excitations

In this section we now focus on the effects induced by the second environment the collective plasmons are subject to [cf. Eq. (1)], i.e., the electronic environment, which is represented by the Hamiltonian (12). Similarly to the photonic bath which we considered in the preceding Secs. IV and V, the coupling between plasmonic and single-particle electronic degrees of freedom, encapsulated in the Hamiltonian (13), leads to two distinct effects. First, the collective plasmons dissipate their energy by producing electron-hole pairs inside each nanoparticles composing the metasurface, corresponding to the well-known Landau damping (Sec. VI.1), and yielding a second (nonradiative) decay channel which adds up to the radiative one. Second, the electronic environment induces an additional renormalization of the quasistatic dispersion relation, which comes on top of the one induced by free photons (Sec. VI.2).

VI.1 Landau damping

We start this section by first evaluating the Landau damping of the collective plasmonic modes. To this end, we treat the coupling (13) between plasmonic and electronic degrees of freedom perturbatively. Since the typical Fermi temperature of ordinary metals is of the order of $10^{4}\,\mathrm{K}$, we employ the zero-temperature limit which is a very good approximation. Within this regime, the Landau damping linewidth $\Gamma_{\tau}^{\varepsilon}(\mathbf{q})$ of the plasmonic eigenmode $|1_{\tau}^{\varepsilon}(\mathbf{q})\rangle$ with band index $\tau$, polarization $\varepsilon$, and wave vector $\mathbf{q}$ is given by the Fermi golden rule

[TABLE]

where

[TABLE]

is given in terms of the Bogoliubov coefficients entering Eq. (15). In Eq. (VI.1), $\Lambda$ is the coupling constant defined in Eq. (14), $\omega_{eh}=(\epsilon_{e}-\epsilon_{h})/\hbar$ corresponds to the frequency associated to an electron-hole pair, where $\epsilon_{e}$ ($\epsilon_{h}$) is the energy of a single-particle electron (hole) state in the self-consistent hard-wall potential associated to each nanoparticle, and $\langle e|\sigma|h\rangle$ is the corresponding dipolar matrix element SM .

Equation (VI.1) can be explicitly calculated SM , yielding for the Laudau damping rate of the collective plasmons the analytical expression

[TABLE]

Here, $v_{\mathrm{F}}$ and $E_{\mathrm{F}}$ are the Fermi velocity and energy of the considered metal, respectively, while the function $g$ is a is a monotonically decreasing function of the parameter $\hbar\omega_{\tau}^{\varepsilon}(\mathbf{q})/E_{\mathrm{F}}$ SM ; Yannouleas and Broglia (1992). To leading order in the coupling constant (7), the Landau damping decay rate of the collective plasmonic modes thus scales with the nanoparticle size as $1/a$, as it is the case for the single-particle result SM ; Yannouleas and Broglia (1992); Kawabata and Kubo (1966); Barma and Subrahmanyam (1989). This is in stark contrast to the radiative linewidth, which increases with the nanoparticle radius as $a^{3}$ [cf. Eq. (V)].

In Fig. 7, we illustrate the result of Eq. (52) for the special case of the OP modes in the honeycomb lattice. In the figure, the displayed results are scaled by the single-particle Landau damping linewidth $\Gamma_{0}=(3v_{\mathrm{F}}/4a)g(\hbar\omega_{0}/E_{\mathrm{F}})$ SM , and the blue (red) curve corresponds to the lower (higher) energy band in Fig. 3(c). We therefore conclude that the higher the energy of the mode, the lower its nonradiative Landau decay rate, as is the case for single nanoparticles Seoanez et al. (2007). Importantly, in contrast to the radiative decay rate analyzed in Sec. V, the Landau damping is nonvanishing over the whole 1BZ, and is of the order of $\Gamma_{0}$. It is therefore crucial to take into account such a decay mechanism, especially for small nanoparticles where it dominates over radiative losses SM .

VI.2 Electronic-induced frequency shift

We now calculate analytically the frequency shift induced by the electronic environment on the plasmonic band structure. Treating the plasmon-electron coupling Hamiltonian (13) to second order in perturbation theory yields for the collective plasmon energy levels $n_{\tau}^{\varepsilon}(\mathbf{q})\hbar\omega_{\tau}^{\varepsilon}(\mathbf{q})+\mathcal{E}_{n_{\tau}^{\varepsilon}(\mathbf{q})}^{(1)}+\mathcal{E}_{n_{\tau}^{\varepsilon}(\mathbf{q})}^{(2)}$. While the first-order correction $\mathcal{E}_{n_{\tau}^{\varepsilon}(\mathbf{q})}^{(1)}$ vanishes due to the selection rules contained in the coupling Hamiltonian (13), the zero-temperature second-order correction reads

[TABLE]

The electronic-induced frequency shift, defined as $\Delta_{\tau}^{\varepsilon}(\mathbf{q})=[\mathcal{E}_{n_{\tau}^{\varepsilon}(\mathbf{q})+1}^{(2)}-\mathcal{E}_{n_{\tau}^{\varepsilon}(\mathbf{q})}^{(2)}]/\hbar$, is therefore given by

[TABLE]

The above expression can be evaluated analytically SM , leading to the result

[TABLE]

where

[TABLE]

Here, $\eta$ is a constant of order $1$ SM . Note that the frequency shift (55) scales with the size of the nanoparticles as $1/a$, up to a logarithmic correction. As is the case for the Landau damping linewidth (52), such a frequency renormalization is therefore of relevance for the smallest nanoparticles. Moreover, the electronic shift of the collective plasmons involves only the plasmonic band structure $\omega_{\tau}^{\varepsilon}(\mathbf{q})$ in contrast to the radiative frequency shifts [cf. Eqs. (IV) and (IV)] which depend on the eigenvectors as well. Notice that substituting $\omega_{0}$ with $\omega_{\tau}^{\varepsilon}(\mathbf{q})$ in Eq. (55) allows us to recover the electronic-induced frequency redshift of an isolated nanoparticle Weick et al. (2006),

[TABLE]

We plot in Fig. 8 the result of Eq. (55) applied to the special case of the OP collective plasmonic modes in the honeycomb lattice. As in Fig. 7, the blue (red) line in Fig. 8 corresponds to the lower (higher) energy band in Fig. 3(c). The displayed results are scaled by the absolute value of the frequency shift (57) corresponding to a single nanoparticle. Immediately noticeable from the figure is that the electronic shift (55) is negative in the entire 1BZ, thus corresponding to a redshift, and is of the order of $\Delta_{0}$. This is in contrast to the radiative frequency shift (IV), whose sign is both depending on the wave vector and polarization of the collective mode (see Sec. IV). Finally, the higher the energy of the mode, the lower is its associated electronic shift. Such a conclusion is reminiscent of what occurs in isolated nanoparticles Seoanez et al. (2007).

VII Observability of the collective plasmonic modes

Experimentally, the ability to observe the plasmonic dispersion relations presented, e.g., in Fig. 5, is governed by the resolution of the separation between the bands with respect to their respective linewidths. The spectral function $\mathcal{A}(\omega,\mathbf{q})$, which characterizes the response of the system to an external perturbation at a given frequency $\omega$ and in-plane wave vector $\mathbf{q}$, is the key quantity which is usually determined in a spectroscopy experiment (using, e.g., photons or hot electrons). Assuming a Breit–Wigner form for $\mathcal{A}(\omega,\mathbf{q})$, we have

[TABLE]

for the IP modes, respectively. In Eq. (58), the renormalized resonance frequency

[TABLE]

takes into account both the radiative shift $\delta_{\tau}^{\varepsilon}(\mathbf{q})$ [cf. Eq. (IV)] and the electronic one $\Delta_{\tau}^{\varepsilon}(\mathbf{q})$ [cf. Eq. (55)]. Additionally, the quantity

[TABLE]

is the total linewidth of the plasmonic modes and includes three distinct contributions: (i) the radiative losses $\gamma_{\tau}^{\varepsilon}(\mathbf{q})$ [Eq. (V)], (ii) the Landau damping $\Gamma_{\tau}^{\varepsilon}(\mathbf{q})$ [Eq. (VI.1)], and (iii) the Ohmic losses, inherent to any (bulk) metal, characterized by the damping rate $\gamma_{\mathrm{O}}$, and which we assume to be mode and size independent.

We show in Fig. 9 the spectral function (58) for both the OP and IP modes and for the square, the honeycomb, and the Lieb lattices along the high symmetry lines of their respective 1BZ as a function of $\omega$. The parameters used in the figure are $\omega_{0}=2.6\,\mathrm{eV/\hbar}$ and $E_{\mathrm{F}}=5.5\,\mathrm{eV}$, corresponding to Ag nanoparticles placed in an effective medium with a refractive index of $2$ Kreibig and Vollmer (1995). The interparticle distance is $d=3a$ and the (reduced) nanoparticle radius is $k_{0}a=0.15$ (corresponding to $a=11\,\mathrm{nm}$). In Fig. 9 the Ohmic losses entering Eq. (60) are neglected.

As can be seen in Fig. 9(a), which displays $\mathcal{A}(\omega,\mathbf{q})$ for the OP modes in the square lattice, the spectral function is reminiscent of the dispersion relation including retardation effects shown in Fig. 5(a), since the electronic shift (55) only induces a finite $\mathbf{q}$-dependent redshift of the band structure. In the figure, we clearly distinguish two different profiles corresponding to wave vectors inside ($|\mathbf{q}|\lesssim k_{0}$) and outside ($|\mathbf{q}|\gtrsim k_{0}$) of the light cone. Within the light cone, the total linewidth (60) is dominated by the radiative damping [see Fig. 6(a)], so that it is difficult to resolve the plasmonic dispersion relation. However, outside of the light cone, only the Landau damping (52) contributes to the total linewidth of the spectral function, allowing for a clear resolution of the plasmonic band. Note that, since the plasmonic modes outside of the light cone are essentially dark, nonoptical techniques are required to excite such modes. Dark modes in nanoparticle dimers and chains have been observed using electron energy loss spectroscopy experiments Chu et al. (2009); Koh et al. (2009); Barrow et al. (2014), and such a technique may be transposed to study plasmonic metasurfaces.

We plot in Figs. 9(b) and 9(c) the spectral function (58a) corresponding to the OP modes for the honeycomb and Lieb lattices, respectively. In both cases, we observe a similar trend as that in Fig. 9(a). Indeed, the entire band structure is clearly visible outside of the light cone in both Figs. 9(b) and 9(c), while the upper bands inside the light cone display a rather large linewidth. However, we see that the lower bands for both the honeycomb and the Lieb lattices are still well resolved inside of the light cone. These dispersion relations correspond to dark plasmonic modes, so that their radiative linewidths for wave vectors $|\mathbf{q}|\lesssim k_{0}$ are nearly vanishing [see the blue line in Fig. 6(b) and the blue and cyan lines in Fig. 6(c)], and only the Landau damping contributes to the observed linewidth.

In Figs. 9(d)–9(f) we display the spectral function (58b) for the IP plasmonic modes for (d) the square, (e) the honeycomb, and (f) the Lieb lattices, respectively. Similar conclusions as that drawn above for the OP modes can be put forward: while the nonradiative bands outside of the light cone are clearly visible, the bright ones inside of the light cone are essentially almost not resolvable.

We conclude this section by discussing the effect of Ohmic damping on the above results. As is evident from Eqs. (58) and (60), the inclusion of a finite $\gamma_{\mathrm{O}}$ leads to a global, size- and mode-independent increase of the width of the spectral lines shown in Fig. 9. From the experiments of Ref. Charlé et al. (1989) on silver clusters, we can extract the value $\gamma_{\mathrm{O}}=70\,\mathrm{meV/\hbar}$, which corresponds to $0.027\omega_{0}$. The increased linewidth induced by Ohmic losses results in a global resolution of the spectral function which is then significantly lower than the one shown in Fig. 9. As a consequence, several plasmonic bands cannot be distinguished properly as it is the case, e.g., for the OP modes in the honeycomb and Lieb lattices in the vicinity of the corners of their respective 1BZ SM . Nevertheless, within the regime of parameters used in Fig. 9, the linewidth of a majority of the plasmonic modes is still sufficiently small to allow for an experimental detection (except for the bright modes within the light cone). Embedding the metasurface in a gain media material Lawandy (2004); Citrin (2006); Maier (2007); Noginov et al. (2007, 2009) should help diminish the effects induced by Ohmic losses on the spectral function and further improve the experimental observability of the collective plasmons.

VIII Conclusion

We have considered the plasmonic properties of metasurfaces constituted by an ordered arbitrary two-dimensional array of spherical metallic nanoparticles. We have focused on the case where the interparticle distance is much smaller than the wavelength associated with the dipolar localized surface plasmon resonance frequency of single nanoparticles, where the near-field, quasistatic dipole–dipole interaction dominates. We have developed a comprehensive open quantum system framework to analyze in full analytical detail the dispersion relations and the lifetimes of the resulting collective plasmonic modes supported by the various metasurfaces which we have studied, including, e.g., the honeycomb and Lieb lattices. Such metasurfaces present appealing topological features in their band structures, such as massless Dirac-like conical dispersions, as well as nearly flat bands, and may be a possible experimental platform to explore new states of hybrid light–matter waves.

Our model enabled us to unveil analytical expressions for the retarded dispersion relations of the plasmonic collective modes, which are in good agreement with existing numerical solutions to Maxwell’s equations. Our model also includes the effects of the particle-hole environment to which such modes are coupled to, and that are of particular relevance for nanoparticles of only a few nanometers in size. Our theory further allowed us to provide analytical expressions for the radiative and nonradiative (Landau) damping rates of the plasmonic modes, which enabled us to critically examine their experimental observability. While Ohmic losses, inherent to the metallic nature of plasmonic particles, may make the detection of the collective modes elusive, the use of gain materials should give scope to their experimental observation.

Acknowledgements.

We thank C. A. Downing, D. Felbacq, K. R. Fratus, F. Gautier, Y. Henry, E. Mariani, P. I. Tamborenea, and D. Weinmann for helpful discussions. We acknowledge financial support from Agence Nationale de la Recherche (Grant No. ANR-14-CE26-0005 Q-MetaMat). This research project has been supported by the University of Strasbourg IdEx program.

Appendix A Diagonalization of the plasmonic Hamiltonian

In this Appendix, we present the details of the diagonalization procedure of the plasmonic Hamiltonian (II). Since we consider large metasurfaces, we use periodic boundary conditions and move to wave vector space using the Fourier transform

[TABLE]

of the bosonic ladder operator $b_{s}^{\sigma}(\mathbf{R}_{s})$ defined in Eq. (5), with $\mathcal{N}=\mathcal{N}_{1}\mathcal{N}_{2}\gg 1$ the total number of unit cells in the lattice. The Hamiltonian (II) then reads

[TABLE]

where $f_{ss^{\prime}}^{\sigma\sigma^{\prime}}(\mathbf{q})$ is the lattice sum defined in Eq. (18). Note that for $s=s^{\prime}$, $f_{ss^{\prime}}^{\sigma\sigma^{\prime}}(\mathbf{q})$ is real due to the inversion symmetry of the Bravais lattice. Moreover, $f_{ss^{\prime}}^{\sigma\sigma^{\prime}}(\mathbf{q})=[f_{s^{\prime}s}^{\sigma\sigma^{\prime}}(\mathbf{q})]^{*}$ since $\boldsymbol{\rho}_{s^{\prime}s}=-\boldsymbol{\rho}_{ss^{\prime}}$.

The plasmonic Hamiltonian $H_{\mathrm{pl}}$ given in Eq. (62) is quadratic, and can thus be diagonalized exactly by means of the bosonic Bogoliubov transformation given in Eq. (15). The inverse transformation reads

[TABLE]

Note that the bosonic commutation relations $[\beta_{\tau}^{\varepsilon}(\mathbf{q}),{\beta_{\tau^{\prime}}^{\varepsilon^{\prime}}}^{\dagger}(\mathbf{q}^{\prime})]=\delta_{\tau\tau^{\prime}}\delta_{\varepsilon\varepsilon^{\prime}}\delta_{\mathbf{q}\mathbf{q}^{\prime}}$ and $[\beta_{\tau}^{\varepsilon}(\mathbf{q}),\beta_{\tau^{\prime}}^{\varepsilon^{\prime}}(\mathbf{q}^{\prime})]=0$ impose that the Bogoliubov coefficients in Eq. (15) fulfill the relations

[TABLE]

and

[TABLE]

where $u_{\tau s}^{\varepsilon\sigma*}(-\mathbf{q})=u_{\tau s}^{\varepsilon\sigma}(\mathbf{q})$ and $v_{\tau s}^{\varepsilon\sigma*}(-\mathbf{q})=v_{\tau s}^{\varepsilon\sigma}(\mathbf{q})$. The dispersion relation $\omega_{\tau}^{\varepsilon}(\mathbf{q})$, as well as the coefficients of the Bogoliubov transformation (15), are then obtained from the Heisenberg equation of motion

[TABLE]

which yields the system of equations (17).

In the case of a generic Bravais lattice (see Sec. III.1), the diagonalization procedure above yields the dispersion relations (20) and (21) for the OP and IP modes, respectively. We find for the corresponding coefficients of the Bogoliubov transformation (15)

[TABLE]

and

[TABLE]

in the case of the OP polarized modes. For the IP modes, the condition $\zeta_{\pm}^{x}(\mathbf{q})\zeta_{\pm}^{y}(\mathbf{q})=\left[2\omega_{0}\Omega f^{xy}(\mathbf{q})\right]^{2}$ must be fulfilled, with $\zeta_{\pm}^{\sigma}(\mathbf{q})=[\omega^{\varepsilon_{\parallel,\pm}}(\mathbf{q})]^{2}-\omega_{0}^{2}-2\omega_{0}\Omega f^{\sigma\sigma}(\mathbf{q})$ ($\sigma=x,y$), so that we obtain

[TABLE]

and

[TABLE]

We point out that neglecting the nonresonant terms in Eq. (62), which corresponds to performing the rotating wave approximation (RWA), yields the same dispersion relations (20) and (21) to first order in $\Omega\ll\omega_{0}$ and the same Bogoliubov coefficients (67) and (69), while $v^{\varepsilon\sigma}(\mathbf{q})=0$ within such an approximation SM .

In the case of a generic bipartite lattice (see Sec. III.2.1), the Bogoliubov coefficients corresponding to the OP modes are given by

[TABLE]

and

[TABLE]

Alternatively to the exact diagonalization procedure presented in this Appendix, it may be useful to treat the nonresonant terms in Eq. (62) perturbatively, since, for all practical purposes, the coupling constant $\Omega\ll\omega_{0}$ [cf. Eq. (7)]. Such a procedure has the advantage of dividing by two the dimension of the system of equations leading to the plasmonic band structure, which may be helpful in deriving the spectrum and the associated Bogoliubov coefficients analytically SM .

Appendix B Mean-field approximation for the quasistatic band structure of the square lattice

Here, we provide an analytical understanding of the cusp presented by the quasistatic collective plasmon dispersion close to the $\Gamma$ point in the case of the OP and IP polarized modes in the case of the square lattice. As demonstrated in Sec. IV, such a cusp is actually an artifact of the quasistatic approximation, and retardation effects strikingly renormalize this behavior. It is however important to get a full description of the quasistatic band structure in order to understand the fully retarded one. The analysis presented below, which relies on the mean-field treatment of the long-ranged quasistatic dipolar interactions, can be straightforwardly extended to other Bravais and non-Bravais lattices.

The behavior of the quasistatic plasmonic dispersion relation of the square lattice for OP polarization close to the $\Gamma$ point mentioned in the main text (see Sec. III.1) can be understood by treating the nearest neighbors in the lattice sum (18) exactly, while averaging the interactions beyond nearest neighbors in the spirit of the mean-field (mf) approximation, leading to $f^{zz}(\mathbf{q})\simeq f_{\mathrm{nn}}^{zz}(\mathbf{q})+f_{\mathrm{mf}}^{zz}(\mathbf{q})$. Here,

[TABLE]

with

[TABLE]

where $J_{n}(z)$ and $\mathbf{H}_{n}(z)$ are the Bessel functions of the first kind and the Struve functions, respectively. Consequently, in the vicinity of the $\Gamma$ point ($|\mathbf{q}|d\ll 1$) and in the weak-coupling regime of interest ($\Omega\ll\omega_{0}$), we find

[TABLE]

so that a cusp appears in the dispersion relation (20). In Fig. 2(c) we show by a dotted line the plasmonic band structure within the mean-field approximation detailed above. As can be seen from the figure, the mean-field approximation accurately describes the cusp of the full band structure in the vicinity of the $\Gamma$ point, while it tends toward the nearest-neighbor approximation away from the $\Gamma$ point.

In the case of the IP-polarized modes, the lattice sums in Eq. (21) are replaced within the mean-field approximation presented above by $f^{\sigma\sigma^{\prime}}(\mathbf{q})\simeq f^{\sigma\sigma^{\prime}}_{\mathrm{nn}}(\mathbf{q})+f^{\sigma\sigma^{\prime}}_{\mathrm{mf}}(\mathbf{q})$ ($\sigma,\sigma^{\prime}=x,y$), where

[TABLE]

and

[TABLE]

The resulting band structure is represented in Fig. 2(d) by a dotted line, and reproduces quite well the full quasistatic dispersion in the vicinity of the $\Gamma$ point. In particular, to first order in $|\mathbf{q}|d\ll 1$, we find $\omega^{\varepsilon_{\parallel,+}}(\mathbf{q})\simeq\omega_{0}+\Omega(-2-\pi/\sqrt{2}+2\pi|\mathbf{q}|d)$ and $\omega^{\varepsilon_{\parallel,-}}(\mathbf{q})\simeq\omega_{0}-\Omega(2+\pi/\sqrt{2})$, demonstrating the presence (absence) of a cusp for the high- (low-) energy plasmonic branch.

Appendix C Comparison of the open quantum system approach to classical electrodynamic calculations

In this Appendix, we compare our analytical results for the retarded plasmonic band structure derived from the perturbative open quantum system approach (see Sec. IV) to numerical calculations based on the solution to Maxwell’s equations Han et al. (2009); Zhen et al. (2008).

In Ref. Han et al. (2009), Han et al. used the multiple scattering theory developed in Refs. Fung and Chan (2007); Zhen et al. (2008) which gives a solution to Maxwell’s equations including retardation effects in order to numerically extract the plasmonic dispersion of a honeycomb array of silver nanoparticles in vacuum, with resonance frequency $\omega_{0}=3.5\,\mathrm{eV/\hbar}$, radius $a=10\,\mathrm{nm}$, and lattice constant $60\,\mathrm{nm}$, corresponding to an interparticle distance $d=35\,\mathrm{nm}$. The classical calculation by Han et al. requires as an input the permittivity of the considered metal, which they chose to be of the Drude-type, without Ohmic losses. We note that the use of a nonlocal dielectric function is not a well-defined concept for a finite-size nanoscale system Raza et al. (2015), while the quantum model developed here is inherently nonlocal. Moreover, the method of Refs. Fung and Chan (2007); Zhen et al. (2008); Han et al. (2009) is fully numerical and hence is limited to specific choices of material parameters and geometry of the metasurface, and involves nonconverging sums that need to be regularized. In contrast, our perturbative approach is straightforward to implement and universal, as it allows to consider any metallic nanoparticles arranged in an arbitrary Bravais or non-Bravais lattice.

We reproduce in Fig. 10 the results of Ref. Han et al. (2009) (green dots), which we compare to our analytical results based on the open quantum system approach (solid lines). Note that the method of Ref. Han et al. (2009) cannot access the plasmonic modes inside of the light cone (i.e., for $|\mathbf{q}|\lesssim k_{0}$). As can be seen from the figure, the quantitative agreement between both theories is rather good for collective modes outside of the light cone, thus validating our approximate treatment of the light–matter coupling, at least in this region of the 1BZ.

We also compared our results to the fully retarded classical calculations by Zhen et al. Zhen et al. (2008) performed on a square lattice of silver particles with radius $a=25\,\mathrm{nm}$ spaced by an interparticle distance $d=75\,\mathrm{nm}$. There, we obtain a qualitative agreement with the classical simulations. The lack of a quantitative agreement in this case can be understood by the rather large spacing between nanoparticles, for which $k_{0}d\simeq 1.4$, which is clearly out of the range of validity of our approach ($k_{0}d\ll 1$).