# Quantum theory of plasmon polaritons in chains of metallic nanoparticles: From near- to far-field coupling regime

Thomas Allard, Guillaume Weick

► To cite this version:

Thomas Allard, Guillaume Weick. Quantum theory of plasmon polaritons in chains of metal- lic nanoparticles: From near- to far-field coupling regime. Physical Review B, 2021, 104 (12), 10.1103/PhysRevB.104.125434. [hal-04075133]

HAL Id: hal-04075133

https://hal.science/hal-04075133v1

Submitted on 26 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Quantum theory of plasmon polaritons in chains of metallic nanoparticles: From near- to far-field coupling regime

Thomas F. Allard and Guillaume Weick

Université de Strasbourg, CNRS, Institut de Physique et Chimie des Matériaux de Strasbourg, UMR 7504, F-67000 Strasbourg, France

We develop a quantum theory of plasmon polaritons in chains of metallic nanoparticles, describing both near- and far-field interparticle distances, by including plasmon–photon umklapp processes. Taking into account the retardation effects of the long-range dipole–dipole interaction between the nanoparticles, which are induced by the coupling of the plasmonic degrees of freedom to the photonic continuum, we reveal the polaritonic nature of the normal modes of the system. We compute the dispersion relation and radiative linewidth, as well as the group velocities of the eigenmodes, and compare our numerical results to classical electrodynamic calculations within the point-dipole approximation. Interestingly, the group velocities of the polaritonic excitations present an almost periodic sign change and are found to be highly tunable by modifying the spacing between the nanoparticles. We show that, away from the intersection of the plasmonic eigenfrequencies with the free photon dispersion, an analytical perturbative treatment of the light–matter interaction is in excellent agreement with our fully retarded numerical calculations. We further study quantitatively the hybridization of light and matter excitations, through an analysis of Hopfield’s coefficients. Finally, we consider the limit of infinitely spaced nanoparticles and discuss some recent results on single nanoparticles that can be found in the literature.

## I. INTRODUCTION

Following the rise of miniaturization of physical systems, the issue of energy and information transport through light within structures with nanoscale dimensions has become increasingly important. Notably, a way to go beyond the diffraction limit encountered in usual dielectric waveguides, such as optical fibers, is highly sought-after, since it is expected to have applications in various ranges of modern physics, from integrated optical circuits or microscopy to biophotonics or data storage [1, 2].

An interesting proposal made more than 20 years ago is to use electrodynamic interparticle interactions in a linear chain of equidistant metallic nanoparticles to be utilized as an effective optical waveguide [3]. It is indeed well known that metallic nanoparticles present a strong absorption peak at the Mie frequency, due to the presence of a collective excitation, the localized surface plasmon (LSP), which corresponds to a dipolar collective oscillation of the conduction electrons of the nanoparticle [4]. Systems based upon such collective excitations have been shown to achieve the confinement and control of light at the nanoscale [5, 6].

On the one hand, linear chains of nanoparticles have been extensively studied in the past through numerous experimental [7–13] and theoretical [14–31] works, especially in the near-field regime where propagation is most promising, that is, when the interparticle distance $d$ is such that $d \ll \lambda_0$, $\lambda_0$ being the wavelength associated with the dipolar LSP excitation. On the other hand, the far-field regime, with $d \gtrsim \lambda_0$, has been much less studied [32], even if it also has interesting properties. Notably, a topological phase transition was recently theoretically unveiled in the far-field regime of a dimerized plasmonic chain [33], which represents a direct topological analog [34] of the regularly spaced chain studied here.

Both the near- and far-field regimes still remain challenging to realize experimentally since it requires cutting-edge technologies. Recently, however, impressive progress has been made in producing highly regular samples of near-field cou-

pled gold nanoparticles with long-range order [35, 36].

After early theoretical studies focused on the quasistatic limit [3, 14, 15], the authors of Refs. [16, 17] noticed that retardation effects of the dipolar interaction along the chain are crucial and imply strong modifications of the energy spectrum. These works, followed by numerous others [19–24, 27, 30, 31], are based on the numerical resolution of the coupled-dipole equations in the point dipole approximation, which arise in the framework of Mie scattering theory. These are therefore fundamentally classical studies, which originate from Maxwell’s equations.

In recent years, however, quantum nanoplasmonics has attracted a lot of attention in the community [37]. The description of nanoscale optical properties of plasmonic metamaterials can exceed the limits of a classical description and requires a more general, quantum-mechanical treatment. For this purpose, approaches using quantum electrodynamics have been developed to tackle the properties of LSPs [38] and linear chains of metallic nanoparticles [26, 28]. A quantum approach turns out to be essential when the size of the particles are sufficiently small to present quantum-size effects, such as Landau damping [39] or electron spill-out [40, 41], namely, for radii smaller than approximately 10–20 nm.

Here, we build upon a Hamiltonian quantum theory of a chain of metallic nanoparticles [28, 29] which has the benefits to be an analytically tractable approach and to be able to readily incorporate quantum effects, such as those previously mentioned. It also permits the use of tools borrowed from quantum optics as well as nonrelativistic quantum electrodynamics, which are familiar to a growing part of the condensed matter community.

Several plasmonic systems have been studied with this quantum theory both in [42–44] and out [29, 45–47] of a finite photonic cavity. While for systems inside a cavity the plasmonic degrees of freedom couple mostly to a finite number of photonic modes, in vacuum the interaction with the whole continuum must be considered. This makes the complete diagonalization of the Hamiltonian more difficult, especially for low-dimensional systems for which crystal momentum is not
conserved for directions transverse to the plasmonic lattice, where translational symmetry is absent. For this reason, until now, a perturbative treatment of the light-matter interaction has been used for these low-dimensional systems, which does not take into account the polaritonic nature of the excitations but presents the advantage of remaining analytically tractable [29]. However, the approach of Ref. [29] remains limited since it only considers the near-field regime. To extend the model of Ref. [29] to nanoparticles with interparticle distances in the far field, the effects of plasmon-photon umklapp processes must be included. To fully take into account the light-matter coupling, we rely on a diagonalization method elaborated by Hopfield [48] in the context of exciton-polaritons. Such a method has already been employed for describing the plasmonic properties of three-dimensional lattices [47], and has been recently successfully compared with experimental results [36].

The system which we study in the following consists of a one-dimensional chain of equidistant resonant dipolar meta-atoms interacting through a long-range retarded dipole-dipole interaction. A typical platform corresponding to our theory is spherical metallic nanoparticles hosting LSPs in the point-dipole approximation, but any other system where dipolar interactions are dominant could also be used, such as, e.g., plasmonic nanorods or microwave helical resonators [42]. In this paper, we focus on the consequences of the light-matter coupling on such a chain and hence on retardation effects of the Coulomb interaction.

Here, we provide a full quantum description of the plasmon-polariton excitations propagating along the chain. Within our approach, we work in the Coulomb gauge such that the scalar and vector potentials describe, respectively, the longitudinal and transverse components of the electromagnetic field. Within this choice of gauge, the scalar potential accounts for purely instantaneous, long-range Coulomb interactions between the LSPs, leading to a collective plasmonic excitation extended over the whole array. The retardation effects arise from the light-matter coupling of the plasmonic modes with the transverse vector potential, which contains the photonic degrees of freedom, present in the vacuum electromagnetic field surrounding the chain. This leads to a hybridization of plasmons and transverse photons into plasmon-polaritons. Notably, the light-matter coupling results in a complex band structure which takes into account radiative losses.

The paper is organized as follows: Section II presents the microscopic quantum model of the plasmonic chain coupled to vacuum electromagnetic modes. In Sec. III, we describe how retardation effects are taken into account, either exactly with the numerical solution of a transcendental dispersion relation or perturbatively, using second-order perturbation theory to find approximate analytical expressions. This allows us to provide both the band structure and the radiative decay rates, as well as the group velocities associated with the normal modes in Sec. IV. In the latter section, we also take advantage of our quantum formalism to examine the hybridization of light and matter degrees of freedom, and we study the limit of infinitely spaced nanoparticles. Finally, we summarize our findings and draw conclusions in Sec. V. Four appendixes

<figure><img src="image_2.png" /><figcaption>FIG. 1. Sketch of a linear chain of identical spherical metallic nanoparticles of radius $a$, separated by a center-to-center distance $d$, arranged along the $z$ direction.</figcaption></figure>

complement the discussion presented in the main text.

## II. MODEL

In this section, we present a microscopic Hamiltonian model of a chain of interacting spherical metallic nanoparticles which host LSPs in the point-dipole approximation. This approach relies on the quantum theory developed in Refs. [28, 29] but with two substantial improvements: First, we consider here plasmon-photon umklapp scattering processes, which will allow us to study both the near- and far-field regimes of interaction. Second, we treat the retardation effects exactly and therefore reveal the polaritonic behavior of the collective excitations along the array. In Appendix A, we present a brief overview of the analogous, widely used macroscopic classical model based on the solution of the fully retarded Maxwell equations to compare it to our quantum approach.

The chain, orientated along the $z$ direction, is embedded in vacuum and comprises $N \gg 1$ nanoparticles<fnref n="1" /> with radius $a$, separated by a center-to-center distance $d$ (see Fig. 1). We consider the particular regime of nanoparticles with radius $a$ much smaller than the inverse wave vector $k_0^{-1}$ associated with the single LSP mode ($k_0a \ll 1$), such that each LSP can be considered as an oscillating point-dipole. We also consider center-to-center distances $d \gtrsim 3a$. For such separations, one can neglect multipolar interactions [15] and tunneling effects between nanoparticles [49]. Hence, each LSP corresponds to a harmonic dipole oscillation of the electronic center of mass at an individual resonance frequency $\omega_0 = k_0c$, with $c$ the speed of light in vacuum. For the case of alkaline nanoparticles in vacuum, $\omega_0 = \omega_p/\sqrt{3}$, where $\omega_p = \sqrt{4\pi n_e e^2/m_e}$ is the plasma frequency, $-e(<0)$ is the electron charge, $m_e$ its mass, and $n_e$ the electronic density.<fnref n="2" />

The LSPs in the chain couple through a long-range dipolar interaction, leading to collective plasmons extended over the whole array. The resulting band structure is then significantly modified by retardation effects. Indeed, due to the large size $(N-1)d + 2a \simeq Nd$ of the chain, one must deal with a

<footnote n="1">¹ Such a limit has been shown to be a good approximation for chains with $N \ge 20$ [16, 28].</footnote>

<footnote n="2">² Throughout this paper, we use cgs units.</footnote>
retarded dipole–dipole Coulomb interaction. In our quantum theory, the latter is encapsulated in the light–matter coupling of the plasmonic modes to a three-dimensional photonic environment [50].

Within our quantum formalism, quantum effects such as Landau damping [51–53] and the associated shift it induces in the plasmonic resonance frequency [40] can be readily incorporated in the theory by considering the interaction of the plasmon with electronic degrees of freedom (specifically, particle–hole excitations inside each nanoparticle). Landau damping has already been computed for a chain with nearest-neighbor interactions in Ref. [28], while the associated electronic-induced redshift can be found in Ref. [54]. Throughout this paper, however, we neglect the effects of an internal electronic environment and hence of Landau damping, since we are primarily focused on radiative effects, which are dominant as long as the nanoparticles are not too small (i.e., their radius should be larger than 10–20 nm [39]).

Using the Coulomb gauge [50, 55], the fully retarded Hamiltonian of the chain coupled to vacuum electromagnetic modes in a volume $V = NdL^2$ of linear size $L \to \infty$ reads [29]

$$H^{\sigma} = H_{\text{pl}}^{\sigma} + H_{\text{ph}} + H_{\text{pl-ph}}^{\sigma}, \quad (1)$$

where $\sigma$ is a fixed parameter accounting for either one of the transverse ($\sigma = x, y$) or the longitudinal ($\sigma = z$) polarization of the plasmonic modes, according to the chosen experimental configuration. We emphasize here the fact that the volume $V$ of the photonic cavity is infinite, the photonic environment representing a continuum of states, in contrast with situations with a finite cavity which could induce strong coupling [44]. We also note that, recently, techniques have been developed to quantize the electromagnetic field not through the Fourier components of the vector potential as we shall be doing here, but directly in position space [56, 57].

The purely plasmonic Hamiltonian describing the LSPs coupled through the long-ranged quasistatic dipole–dipole interaction is

$$H_{\text{pl}}^{\sigma} = \hbar\omega_0 \sum_{n=1}^{N} b_n^{\sigma\dagger} b_n^{\sigma} + \frac{\hbar\Omega}{2} \sum_{\substack{n,m=1 \\ (n \neq m)}}^{N} \frac{\eta^{\sigma}}{|n-m|^3} (b_n^{\sigma} + b_n^{\sigma\dagger}) (b_m^{\sigma} + b_m^{\sigma\dagger}). \quad (2)$$

Here, the bosonic operator $b_n^\sigma (b_n^{\sigma\dagger})$ annihilates (creates) an LSP with polarization $\sigma$ on nanoparticle $n \in [1, N]$. The long-range dipolar coupling strength between LSPs scales as the inverse interparticle distance cubed and is given by

$$\Omega = \frac{\omega_0}{2} \left(\frac{a}{d}\right)^3, \quad (3)$$

and the polarization-dependent factor $\eta^{x,y} = 1$ ($\eta^z = -2$) for the transverse (longitudinal) mode accounts for the anisotropy of the quasistatic dipole–dipole interaction.

The photonic environment in Eq. (1) is described by the Hamiltonian

$$H_{\text{ph}} = \sum_{l} \sum_{\mathbf{k}, \lambda_{\mathbf{k}}^{l}} \hbar \nu_{\mathbf{k}}^{l} c_{\mathbf{k}}^{l, \lambda_{\mathbf{k}}^{l \dagger}} c_{\mathbf{k}}^{l, \lambda_{\mathbf{k}}^{l}}, \quad (4)$$

where the upper index $l$ on $\mathbf{k}$-dependent variables is a short form signifying that the quantity depends on $\mathbf{k} - \mathbf{G}^l$, with $\mathbf{G}^l = 2\pi l \hat{z}/d$ ($l \in \mathbb{Z}$) representing the set of reciprocal lattice vectors. Hence the ladder operator $c_{\mathbf{k}}^{l, \lambda_{\mathbf{k}}^{l}} (c_{\mathbf{k}}^{l, \lambda_{\mathbf{k}}^{l \dagger}})$ annihilates (creates) a photon with wave vector $\mathbf{k} - \mathbf{G}^l$ in the first Brillouin zone, transverse polarization $\lambda_{\mathbf{k}}^{l} = \lambda_{\mathbf{k}-\mathbf{G}^l}$, i.e., $(\mathbf{k} - \mathbf{G}^l) \cdot \lambda_{\mathbf{k}-\mathbf{G}^l} = 0$, and angular frequency

$$\nu_{\mathbf{k}}^{l} = c|\mathbf{k} - \mathbf{G}^{l}|, \quad (5)$$

where hats designate unit vectors. The two photon polarizations are parameterized as $\lambda_{\mathbf{k}}^{l,(1)} = [\hat{z} \times (\mathbf{k} - \mathbf{G}^l)]/|\hat{z} \times (\mathbf{k} - \mathbf{G}^l)|$ and $\lambda_{\mathbf{k}}^{l,(2)} = [(\mathbf{k} - \mathbf{G}^l) \times \lambda_{\mathbf{k}}^{l,(1)}]/|(|\mathbf{k} - \mathbf{G}^l) \times \lambda_{\mathbf{k}}^{l,(1)}|$ for the case $\mathbf{k} - \mathbf{G}^l \neq \hat{z}$, while we choose $\lambda_{\mathbf{k}}^{l,(1)} = \hat{x}$ and $\lambda_{\mathbf{k}}^{l,(2)} = \hat{y}$ for the case $\mathbf{k} - \mathbf{G}^l = \hat{z}$. We note that in previous works [28, 29], only the photonic band $l = 0$ was considered. This amounts to neglect umklapp processes and limits the theory to $k_0 d \ll 1$, i.e., to short center-to-center distances $d$. Taking into account umklapp processes is therefore essential to describe interacting nanoparticles in the far-field region.

The plasmon–photon coupling Hamiltonian in Eq. (1) reads in the long-wavelength point-dipole approximation ($|\mathbf{k}|a \ll 1$) [50]

$$H_{\text{pl-ph}}^{\sigma} = \frac{e}{m_e c} \sum_{n=1}^{N} \Pi_{n}^{\sigma} \cdot \mathbf{A}(\mathbf{d}_{n}), \quad (6)$$

where $\mathbf{d}_n = d(n-1)\hat{z}$ corresponds to the location of the center of nanoparticle $n$. The momentum associated with the LSP in nanoparticle $n$ reads

$$\Pi_{n}^{\sigma} = i \sqrt{\frac{N_e m_e \hbar \omega_0}{2}} (b_{n}^{\sigma\dagger} - b_{n}^{\sigma}) \hat{\sigma}, \quad (7)$$

where $N_e$ is the number of electrons in a given nanoparticle, while the quantized vector potential is given by

$$\mathbf{A}(\mathbf{d}_n) = \sum_l \sum_{\mathbf{k}, \lambda_{\mathbf{k}}^l} \hat{\lambda}_{\mathbf{k}}^l \sqrt{\frac{2\pi\hbar c^2}{\mathcal{V}\nu_{\mathbf{k}}^l}} \times (c_{\mathbf{k}}^{l, \lambda_{\mathbf{k}}^l} e^{i(\mathbf{k} - \mathbf{G}^l) \cdot \mathbf{d}_n} + \text{H.c.}) \Theta(\omega_c - \nu_{\mathbf{k}}^l). \quad (8)$$

In the equation above, we have introduced a Heaviside step function $\Theta(x)$ to take into account the finite size of the nanoparticles and to switch off the light–matter interaction for photonic modes with a frequency $\nu_{\mathbf{k}}^l$ larger than the cutoff frequency $\omega_c = c/a$, for which the point-dipole approximation breaks down.<fnref n="3" />

<footnote n="3">³ In our situation, the point-dipole approximation amounts to considering</footnote>
A significant simplification we are making in writing the coupling Hamiltonian (6) is to disregard the quadratic $A^2$ term, known as the diamagnetic term, which accounts for a photon self-interaction energy and does not contribute to retardation effects. In the light-matter coupling regime considered in this paper, the diamagnetic term can be safely neglected, as it would only lead to very slight changes to the results presented here. Details about the $A^2$ term are given in Appendix B, as well as justifications of the above-mentioned simplification.

### III. EXACT AND PERTURBATIVE TREATMENTS OF THE LIGHT-MATTER COUPLING

#### A. Full diagonalization of the polaritonic Hamiltonian

The light-matter coupling Hamiltonian (6) includes the retardation effects of the dipole-dipole interaction along the chain, which modify the nature of the normal modes which become plasmon-polaritons. Taking such a coupling into account is therefore essential and will be the subject of this section.

Considering the infinite chain limit with $N \gg 1$, it is convenient to choose, without loss of generality, periodic boundary conditions. We then move into wave vector space using the Fourier transform $b_n^\sigma = N^{-1/2} \sum_q e^{inqd} b_q^\sigma$, where the plasmonic wave number is $q = 2\pi p/Nd$, with the integer $p \in [-N/2, +N/2]$. It is important to note that the periodicity of the chain in the $z$ direction implies that the longitudinal component of the photonic wave vector $\mathbf{k}$ is conserved with the plasmonic one. We thus choose to write $k_z = q$ and, for commodity, we denote the $\mathbf{k}$-dependency as $\mathbf{k} = (\mathbf{\kappa}, q)$ where $\mathbf{\kappa} = (k_x, k_y)$.

The Hamiltonian (1) can hence be rewritten in Fourier space as $H^\sigma = \sum_q H_q^\sigma$, where the $q$-dependent polaritonic Hamiltonian is

$$H_q^\sigma = H_{\text{pl},q}^\sigma + H_{\text{ph},q} + H_{\text{pl-ph},q}^\sigma. \quad (9)$$

The plasmonic part can be written as

$$H_{\text{pl},q}^{\sigma} = \hbar\omega_0 b_q^{\sigma\dagger} b_q^{\sigma} + \frac{\hbar\Omega}{2} [f_q^{\sigma} b_q^{\sigma\dagger} (b_q^{\sigma} + b_{-q}^{\sigma\dagger}) + \text{H.c.}], \quad (10)$$

the photonic one as

$$H_{\text{ph},q} = \sum_l \sum_{\mathbf{\kappa}, \hat{\lambda}_{\mathbf{\kappa}q}^l} \hbar \nu_{\mathbf{\kappa}q}^l c_{\mathbf{\kappa}q}^{l, \hat{\lambda}_{\mathbf{\kappa}q}^l \dagger} c_{\mathbf{\kappa}q}^{l, \hat{\lambda}_{\mathbf{\kappa}q}^l}, \quad (11)$$

that exponentials of the type $e^{ik \cdot r}$ equal unity, where $r$ is the position of the electrons in the nanoparticle. Since at most $r \simeq a$, we require that $k \ll 1/a$. Therefore, we consider that the point-dipole approximation is accurate up to the wave number $k_c = \Lambda/a$, where $\Lambda \lesssim 1$ is a cutoff parameter which must be close to unity. In what follows, we set $\Lambda = 1$.

and the coupling as

$$H_{\text{pl-ph},q}^{\sigma} = i\hbar\omega_0 \sum_l \sum_{\mathbf{\kappa}, \hat{\lambda}_{\mathbf{\kappa}q}^l} \xi_{\mathbf{\kappa}q}^l (\hat{\sigma} \cdot \hat{\lambda}_{\mathbf{\kappa}q}^l) \times \left[ e^{-iqd} c_{\mathbf{\kappa}q}^{l, \hat{\lambda}_{\mathbf{\kappa}q}^l} (b_q^{\sigma\dagger} - b_{-q}^{\sigma\dagger}) - \text{H.c.} \right]. \quad (12)$$

The lattice sum $f_q^\sigma$ in Eq. (10) can be expressed in closed form in terms of the polylogarithm function $Li_5(z) = \sum_{n=1}^\infty z^n/n^s$ as

$$f_q^\sigma = \eta^\sigma [\text{Li}_3(e^{iqd}) + \text{Li}_3(e^{-iqd})], \quad (13)$$

and the light-matter coupling strength in Eq. (12) is encapsulated in

$$\xi_{\mathbf{\kappa}q}^l = \sqrt{\frac{\omega_0 \pi a^3}{\nu_{\mathbf{\kappa}q}^l d L^2}} \Theta(\omega_c - \nu_{\mathbf{\kappa}q}^l). \quad (14)$$

We start by diagonalizing exactly the Hamiltonian (9). Toward this goal, we define the bosonic plasmon-polariton annihilation operator $\mu_q^\sigma$ as a linear combination of plasmon and photon ladder operators,

$$\mu_q^\sigma = W_q^\sigma b_q^\sigma + X_q^\sigma b_{-q}^{\sigma\dagger} + \sum_l \sum_{\mathbf{\kappa}, \hat{\lambda}_{\mathbf{\kappa}q}^l} \left( Y_{\mathbf{\kappa}q}^{l, \sigma, \hat{\lambda}_{\mathbf{\kappa}q}^l} c_{\mathbf{\kappa}q}^{l, \hat{\lambda}_{\mathbf{\kappa}q}^l} + Z_{\mathbf{\kappa}q}^{l, \sigma, \hat{\lambda}_{\mathbf{\kappa}q}^l} c_{\mathbf{\kappa}, -q}^{l, \hat{\lambda}_{\mathbf{\kappa}, -q}^l} \right), \quad (15)$$

where the Hopfield coefficients $W_q^\sigma$, $X_q^\sigma$, $Y_{\mathbf{\kappa}q}^{l, \sigma, \hat{\lambda}_{\mathbf{\kappa}q}^l}$, and $Z_{\mathbf{\kappa}q}^{l, \sigma, \hat{\lambda}_{\mathbf{\kappa}q}^l}$ are a priori complex. Due to the bosonic nature of plasmon-polaritons, the coefficients are normalized following the condition

$$|W_q^\sigma|^2 - |X_q^\sigma|^2 + \sum_l \sum_{\mathbf{\kappa}, \hat{\lambda}_{\mathbf{\kappa}q}^l} \left( |Y_{\mathbf{\kappa}q}^{l, \sigma, \hat{\lambda}_{\mathbf{\kappa}q}^l}|^2 - |Z_{\mathbf{\kappa}q}^{l, \sigma, \hat{\lambda}_{\mathbf{\kappa}q}^l}|^2 \right) = 1. \quad (16)$$

We note that in the one-dimensional geometry which we consider here, the polariton polarization is aligned with the plasmonic one, $\sigma$. In higher-dimensional metamaterials, this is generally not the case [39, 47].

We now use the ansatz

$$\mu_q^\sigma(t) = \mu_q^\sigma e^{-i\Omega_q^\sigma t}, \quad (17)$$

where, in the following, $\Omega_q^\sigma$ are referred to as the polaritonic eigenfrequencies of the system, to solve the Heisenberg equation of motion

$$\dot{\mu}_q^\sigma(t) = \frac{i}{\hbar} [H^\sigma, \mu_q^\sigma(t)], \quad (18)$$

that we can rewrite as $\hbar\Omega_q^\sigma \mu_q^\sigma = [\mu_q^\sigma, H^\sigma]$. We then compute the above commutator using Eq. (15) and the bosonic commutation relations for the plasmonic and photonic ladder operators, respectively, $[b_q^\sigma, b_{q'}^{\sigma\dagger}] = \delta_{q,q'} \delta_{\sigma,\sigma'}$ and
$[c_{\kappa q}^{l, \hat{\lambda}_{\kappa q}^{l,(i)}}, c_{\kappa' q'}^{l', \hat{\lambda}_{\kappa' q'}^{l', (j)}}]^{\dagger} = \delta_{q, q'} \delta_{\kappa, \kappa'} \delta_{l, l'} \delta_{i, j}$. Since the ladder op-
erators are independent from each other, we use them to fac-

torize the equation and separate it in a system of four coupled
equations with five unknown variables: the four Hopfield co-
efficients and the desired eigenfrequency Ω<sub>q</sub><sup>σ</sup> for a given wave
number q,

$$
W_q^\sigma = - \frac{1}{\Omega_q^\sigma - \omega_0 - \Omega f_q^\sigma} \left\{ \Omega f_q^\sigma X_q^\sigma + i\omega_0 e^{iqd} \sum_l \sum_{\kappa, \hat{\lambda}_{\kappa q}^l} \left[ \xi_{\kappa q}^l (\hat{\sigma} \cdot \hat{\lambda}_{\kappa q}^l) Y_{\kappa q}^{l, \sigma, \hat{\lambda}_{\kappa q}^l} - \xi_{\kappa, -q}^l (\hat{\sigma} \cdot \hat{\lambda}_{\kappa, -q}^l) Z_{\kappa q}^{l, \sigma, \hat{\lambda}_{\kappa q}^l} \right] \right\}, \quad (19a)
$$

$$
X_q^\sigma = \frac{1}{\Omega_q^\sigma + \omega_0 + \Omega f_q^\sigma} \left\{ \Omega f_q^\sigma W_q^\sigma i\omega_0 e^{iqd} \sum_l \sum_{\kappa, \hat{\lambda}_{\kappa q}^l} \left[ \xi_{\kappa q}^l (\hat{\sigma} \cdot \hat{\lambda}_{\kappa q}^l) Y_{\kappa q}^{l, \sigma, \hat{\lambda}_{\kappa q}^l} - \xi_{\kappa, -q}^l (\hat{\sigma} \cdot \hat{\lambda}_{\kappa, -q}^l) Z_{\kappa q}^{l, \sigma, \hat{\lambda}_{\kappa q}^l} \right] \right\}, \quad (19b)
$$

$$
Y_{\kappa q}^{l, \sigma, \hat{\lambda}_{\kappa q}^{l}} = i\omega_0 \frac{e^{-iqd} \xi_{\kappa q}^{l} (\hat{\sigma} \cdot \hat{\lambda}_{\kappa q}^{l})}{\Omega_{q}^{\sigma} - \nu_{\kappa q}^{l}} (W_{q}^{\sigma} + X_{q}^{\sigma}), \quad (19c)
$$

$$
Z_{\kappa q}^{l, \sigma, \hat{\lambda}_{\kappa q}^{l}} = i\omega_0 \frac{e^{-iqd} \xi_{\kappa, -q}^{l} (\hat{\sigma} \cdot \hat{\lambda}_{\kappa, -q}^{l})}{\Omega_{q}^{\sigma} + \nu_{\kappa, -q}^{l}} (W_{q}^{\sigma} + X_{q}^{\sigma}). \quad (19\text{d})
$$

Solving the system of equations (19) for Ω<sub>q</sub><sup>σ</sup> leads to the fol-
lowing transcendental equation,

$$
\frac{\Omega_q^{\sigma 2} - \omega_q^{\sigma 2}}{4\omega_0 \omega_q^{\sigma 2}} = \sum_l \sum_{\kappa, \hat{\lambda}_{\kappa q}^l} \frac{\nu_{\kappa q}^l (\hat{\sigma} \cdot \hat{\lambda}_{\kappa q}^l)^2 (\xi_{\kappa q}^l)^2}{\Omega_q^{\sigma 2} - \nu_{\kappa q}^{l 2}}, \quad (20)
$$

where the summation over $\kappa$ excludes singular terms. Note
that $\Omega_q^\sigma$ is a *priori* complex, its imaginary part taking into ac-
count the radiative decay of the collective excitation inside the
three-dimensional photonic continuum, leading to a radiative
contribution to the linewidth [cf. Eq. (17)].

Setting the light–matter coupling (14) to zero in Eq. (20)
naturally leads to $\Omega_q^\sigma = \omega_q^\sigma$, where

$$
\omega_q^\sigma = \omega_0 \sqrt{1 + 2 \frac{\Omega}{\omega_0} f_q^\sigma} \quad (21)
$$

is the quasistatic spectrum of the collective plasmonic
modes. Such a quasistatic dispersion can also be obtained
[29] through a direct diagonalization of the Hamiltonian
(2) by means of a bosonic Bogoliubov transformation (see
Sec. III B).

In Eq. (20), we can carry out the summation over photon polarizations using the relation

$$
\sum_{\hat{\lambda}_{\kappa q}^l} \left( \hat{\sigma} \cdot \hat{\lambda}_{\kappa q}^l \right)^2 = 1 - \left( \hat{\sigma} \cdot \frac{\mathbf{k} - \mathbf{G}^l}{|\mathbf{k} - \mathbf{G}^l|} \right)^2 . \quad (22)
$$

To then compute the summation over $\kappa$, we recall that we
consider a quantization volume $V \to \infty$ such that we can take
the continuum limit

$$
\frac{(2\pi)^2}{L^2} \sum_{\kappa} \rightarrow \mathcal{P} \int_{0}^{\infty} \kappa d\kappa \int_{0}^{2\pi} d\theta, \quad (23)
$$

where we used polar coordinates and $\mathcal{P}$ represents a principal value. The angular integral is trivial, however, in the radial one, the integrand suffers from an ultraviolet logarithmic divergence. The latter is naturally regularized with the help of the UV frequency cutoff $\omega_c$ appearing in the light–matter coupling (14) and which takes into account the finite size of the nanoparticles. The cutoff also selects how many photonic bands $l$ interact with the plasmonic chain, leaving us with a finite summation over the index $l$, from $-l_{\max}$ to $+l_{\max}$. In the first Brillouin zone, one can show that using a cutoff $\omega_c = c/a$ leads to

$$
l_{\max} = \left\lfloor \frac{d}{2\pi a} + \frac{1}{2} \right\rfloor, \qquad (24)
$$

where $\lfloor x \rfloor$ denotes the floor function. After a long but straight-
forward calculation, Eq. (20) translates into

$$
\begin{equation}
\begin{split}
\Omega_q^{\sigma 2} - \omega_q^{\sigma 2} ={}& \eta^{\sigma} \omega_0^2 \omega_q^{\sigma 2} \frac{a^3}{dc^2} \sum_{l=-l_{\text{max}}}^{l_{\text{max}}} \left( \frac{cq^l}{\Omega_q^{\sigma}} \right)^2 \left\{ \ln \left( \frac{\omega_c}{c|q^l|} \right) \right. \\
& \left. + \frac{1}{2} \left[ 1 + \operatorname{sgn} \{\eta^{\sigma}\} \left( \frac{\Omega_q^{\sigma}}{cq^l} \right)^2 \right] \right. \\
& \left. \times \log^{\phi} \left( \frac{(cq^l)^2 - \Omega_q^{\sigma 2}}{\omega_c^2 - \Omega_q^{\sigma 2}} \right) \right\} \Theta(\omega_c - c|q^l|),
\end{split}
\tag{25}
\end{equation}
$$

where we used the short form $q^l = q - 2\pi l/d$. The use of a
simple root-finding algorithm on Eq. (25) finally leads to the
complex eigenfrequencies $\Omega_q^\sigma$.

Importantly, the logarithms log<sup>φ</sup>(z) present in the third line
of the above equation are complex logarithms, and hence a
choice φ of branch cut must be made. We recall that one
can define any complex logarithm as log<sup>φ</sup>(z) = ln(|z|) +
i arg(ze⁻¹iφ) + iφ, with φ = φ(|z|) being any real function of the modulus of z, ln(x) the natural logarithm and Arg(z) the principal value of the argument. In our case, the use of the usual principal value logarithm, corresponding to the choice φ = 0, hides one of the roots behind the branch cut. To ensure that our algorithm can readily find the complex roots, and that they have a physical meaning, i.e., a positive (negative) real (imaginary) part [cf. Eq. (17)], we have to modify the branch cut and we choose φ = -π/2, leading to a logarithm which is continuous on the open set C \ {iR⁺}. Within such a choice, the root is no longer hidden and the algorithm can readily converge.

## B. Second-order perturbation theory

The approach developed in the previous subsection allows us to take into account the retardation effects in an exact manner, diagonalizing the open quantum system Hamiltonian (1) to find the complex eigenfrequencies of the system. However, it requires the numerical solution of the transcendental dispersion equation (25), in a similar way to the resolution using the classical approach (see Appendix A). An advantage of our quantum approach is that we can use a perturbative treatment of the light-matter interaction, considering the interaction of the photonic continuum as a weak perturbation to the plasmonic system, to find adequate analytical approximations for the eigenfrequencies and the radiative decay rates of the system.

In this subsection, we rely on the work presented in Ref. [29], with the substantial improvement of taking into account umklapp processes, which enables us to study both near- and far-field regimes.

### 1. Radiative frequency shift

We treat the coupling Hamiltonian (6) up to second order in perturbation theory. For this purpose, we begin by diagonalizing the plasmonic Hamiltonian (2) by means of a Bogoliubov transformation. As detailed in Ref. [29], one can rewrite the latter as $H_{pl}^{\sigma} = \sum_{q} \hbar \omega_{q} B_{q}^{\sigma \dagger} B_{q}^{\sigma}$, where the quasistatic dispersion $\omega_{q}^{\sigma}$ is given in Eq. (21), and the bosonic Bogoliubov ladder operator $B_{q}^{\sigma}$ ($B_{q}^{\sigma \dagger}$) acts on an eigenstate $|n_{q}^{\sigma}\rangle$ of the Hamiltonian (10) representing $n_{q}^{\sigma}$ quanta occupying the collective plasmon mode with polarization $\sigma$, wave vector $q$, and eigenfrequency $\omega_{q}^{\sigma}$, as $B_{q}^{\sigma}|n_{q}^{\sigma}\rangle = \sqrt{n_{q}^{\sigma}}|n_{q}^{\sigma} - 1\rangle$ ($B_{q}^{\sigma \dagger}|n_{q}^{\sigma}\rangle = \sqrt{n_{q}^{\sigma} + 1}|n_{q}^{\sigma} + 1\rangle$). Explicitly, one has

$$B_{q}^{\sigma} = \frac{\omega_{q}^{\sigma} + \omega_{0}}{2\sqrt{\omega_{0}\omega_{q}^{\sigma}}} b_{q}^{\sigma} + \frac{\omega_{q}^{\sigma} - \omega_{0}}{2\sqrt{\omega_{0}\omega_{q}^{\sigma}}} b_{-q}^{\sigma \dagger}. \quad (26)$$

For a given polarization $\sigma$ and mode $q$, the perturbed plasmonic energy levels are $E_{n_{q}^{\sigma}} = E_{n_{q}^{(0)}} + E_{n_{q}^{(1)}} + E_{n_{q}^{(2)}}$. The renormalized plasmonic eigenfrequency within perturbation theory (pt) $\Omega_{q}^{\sigma, \text{pt}}$ is defined as the difference between successive plasmonic energy levels, $\Omega_{q}^{\sigma, \text{pt}} = (E_{n_{q}^{\sigma}+1} - E_{n_{q}^{\sigma}})/\hbar$, and we write

it as

$$\Omega_{q}^{\sigma, \text{pt}} = \omega_{q}^{\sigma} + \delta_{q}^{\sigma}. \quad (27)$$

Since the unperturbed energy is $E_{n_{q}^{(0)}} = n_{q}^{\sigma} \hbar \omega_{q}^{\sigma}$ and the first-order contribution $E_{n_{q}^{(1)}} = 0$, the radiative frequency correction reads $\delta_{q}^{\sigma} = (E_{n_{q}^{(2)}} - E_{n_{q}^{(2)}})/\hbar$ and can be written as

$$\delta_{q}^{\sigma} = 2\omega_{0}^{2}\omega_{q}^{\sigma} \frac{\pi a^{3}}{dL^{2}} \sum_{l} \sum_{\kappa, \hat{\lambda}_{\kappa q}^{l}} \frac{(\hat{\sigma} \cdot \hat{\lambda}_{\kappa q}^{l})^{2} \Theta(\omega_{c} - \nu_{\kappa q}^{l})}{\omega_{q}^{\sigma 2} - \nu_{\kappa q}^{l 2}}, \quad (28)$$

where the summation over $\kappa$ excludes singular terms. We note that the only correction to the plasmonic energy levels which contributes is the second-order one, which corresponds to the emission and reabsorption of virtual photons by the plasmonic eigenstate $|n_{q}^{\sigma}\rangle$.

Similarly as in the previous subsection, one can carry out the summation over photon polarization using the relation (22) and taking the continuum limit (23) to write the summation over $\kappa$ as a principal value integral using polar coordinates. We then obtain

$$\delta_{q}^{\sigma} = \eta^{\sigma} \frac{\omega_{0}^{2} \omega_{q}^{\sigma}}{2} \frac{a^{3}}{dc^{2}} \sum_{l} \Theta(\omega_{c} - c|q^{l}|) \left(\frac{cq^{l}}{\omega_{q}^{\sigma}}\right)^{2} \left\{ \ln\left(\frac{\omega_{c}}{c|q^{l}|}\right) + \frac{1}{2} \left[ 1 + \operatorname{sgn}\{\eta^{\sigma}\} \left(\frac{\omega_{q}^{\sigma}}{cq^{l}}\right)^{2}\right] \ln\left(\left|\frac{(cq^{l})^{2} - \omega_{q}^{\sigma 2}}{\omega_{c}^{2} - \omega_{q}^{\sigma 2}}\right|\right) \right\}. \quad (29)$$

Note that the above radiative frequency shift, aside from a factor of $2\omega_{q}^{\sigma}$, has the same expression as the right-hand side of Eq. (25), except that the complex polaritonic eigenfrequency $\Omega_{q}^{\sigma}$ has been replaced by the quasistatic plasmonic frequency $\omega_{q}^{\sigma}$, and hence the complex logarithms have become natural logarithms. We further note that due to the plasmon-photon umklapp processes, and unlike what has been done in previous works [28, 29], the $2\pi/d$-periodicity of the eigenfrequencies $\Omega_{q}^{\sigma, \text{pt}}$ is conserved.

### 2. Radiative decay rate

Treating the light-matter coupling Hamiltonian (6) as a weak perturbation, one can use Fermi's golden rule to obtain the radiative decay rates $\gamma_{q}^{\sigma, \text{pt}}$ of the collective plasmonic excitation into the surrounding photonic environment. Using the expression (12) of the coupling Hamiltonian, we find that these decay rates, which also represent the radiative contribution to the inverse lifetime of a given mode $q$, read

$$\gamma_{q}^{\sigma, \text{pt}} = 2\pi^2 \omega_0^2 \omega_q^\sigma \frac{a^3}{dL^2} \sum_l \sum_{\kappa, \hat{\lambda}_{\kappa q}^l} \frac{(\hat{\sigma} \cdot \hat{\lambda}_{\kappa q}^l)^2}{\nu_{\kappa q}^l} \delta(\omega_q^\sigma - \nu_{\kappa q}^l) \\ \times \Theta(\omega_c - \nu_{\kappa q}^l). \quad (30)$$
After taking the continuum limit (23) and completing the integration, the above equation can be written as

$$ \gamma_q^{\sigma,pt} = \frac{\pi \eta^{\sigma} \omega_0^2 a^3}{2 \omega_q^{\sigma} d} \sum_l (q^l)^2 \left[ 1 + \operatorname{sgn}\{\eta^{\sigma}\} \left( \frac{\omega_q^{\sigma}}{cq^l} \right)^2 \right] \times \Theta(\omega_q^{\sigma} - c|q^l|). \quad (31) $$

## IV. RESULTS

We now present the results obtained through our numerical quantum approach presented in Sec. III A, which we compare to the analytical perturbative results presented in Sec. III B and to the classical model widely used in the literature and outlined in Appendix A. Specifically, the classical model which we employ uses the nanoparticle polarizability of Eq. (A3) that is obtained from the first (electric dipole) Mie coefficient [58–60], see Eq. (A2).

We begin this section by displaying our results for the polaritonic band structures and the radiative decay rates (Sec. IV A) and then discuss the corresponding group velocities (Sec. IV B). We then quantitatively demonstrate the hybridization of the light and matter degrees of freedom through a detailed discussion of the Hopfield coefficients (Sec. IV C), and, finally, we discuss the limit of infinitely spaced nanoparticles, i.e., of single nanoparticles (Sec. IV D).

### A. Polaritonic band structure and radiative decay rate

In this subsection, we present the results obtained from the numerical resolution of Eq. (25) for the band structure Re{$\Omega_q^\sigma$} and the radiative decay rate $\gamma_q^\sigma = -2\operatorname{Im}\{\Omega_q^\sigma\}$, for the longitudinal ($\sigma = z$) and transverse ($\sigma = x, y$) polarizations. The factor -2 in the above expression of the decay rate is such that $\gamma_q^\sigma$ corresponds to the radiative contribution to the inverse lifetime of a given mode of the system [see Eq. (17)].<fnref n="4" /> In Figs. 2–4 below, we compare these results with the perturbative ones, given by Eqs. (27) and (31), and to the results from the classical model obtained through the numerical resolution of Eq. (A5).

We choose a value of $k_0a = 0.3$ which corresponds to, e.g., nanoparticles with a Mie frequency $\omega_0 = 3$ eV/ħ and a radius $a = 20$ nm. We display the eigenfrequencies and decay rates in units of the bare frequency $\omega_0$ and as a function of the reduced longitudinal wave number $qd$. In this way, we consider $\omega_0$ and hence the plasma frequency $\omega_p$ as (essentially material-dependent) parameters of our theory. Only the positive half of the first Brillouin zone is shown, since all of the plotted quantities are symmetric around $q = 0$ and $2\pi/d$-periodic.

### 1. Near-field coupled nanoparticles

We first discuss in Fig. 2 the case of a center-to-center distance $d = 3a$. For such a separation, the LSPs supported by the nanoparticles are near-field coupled.

Compared to the quasistatic dispersion from Eq. (21) [black solid lines in Figs. 2(a) and 2(b)], the fully retarded one (red solid lines) experiences a redshift, especially in the area inside the first light cone, i.e., for frequencies larger than the first ($l = 0$) light line $c|q|$ (depicted by a black dotted line). In this area of the Brillouin zone, the light-matter interaction is significant and the plasmons are radiating. Thus, a finite radiative linewidth appears [see Figs. 2(c) and 2(d)], which vanishes at the wave number where the dispersion relation intersects the light line. We note that at this precise point, the dispersion curve for the longitudinal polarization [see Fig. 2(a)] presents a discontinuity, i.e., no roots are found. Taking into account Ohmic losses, however, leads to a definite root at the intersection [30, 61], where the slope of the curve changes drastically to be equal to the slope of the light line, namely, $c$. For larger values of $qd$, the plasmonic excitation is interacting with photonic modes of higher energy and the normal modes correspond to guided ones, which are immune to radiation damping.

From the results displayed in Figs. 2(a) and 2(c), we see that for the longitudinal polarization there are almost no differences between our fully retarded eigenfrequencies obtained from Eq. (25) and the analytical second-order perturbation result (27) (dashed blue lines). Likewise, we only have a slight increase of the decay rate $\gamma_q^z$ in comparison with the Fermi golden rule result [Fig. 2(c)]. The longitudinally polarized plasmonic chain is thus very well described by our analytical theory. This is due to the fact that in the Coulomb gauge, only transverse photons are present and hence they interact weakly with the longitudinally polarized plasmons.

For the transverse polarizations [Figs. 2(b) and 2(d)], the full treatment of the light-matter interaction implies the appearance of two distinct bands in the dispersion, anticrossing at the light line. This is a typical avoided crossing dispersion, which is a signature of the presence of polaritonic excitations [48]. Here and in the following, we denote as the upper band the one which is above the light line anticrossing the dispersion (here $c|q|$), and the lower band the one below. We see in Figs. 2(b) and 2(d) that for $d = 3a$, the upper band represents modes with radiating polaritons inside the light cone [i.e., with a finite decay rate, cf. the upper red line in Fig. 2(d)] and stops at the anticrossing, while the lower band is free of radiation damping in the entire Brillouin zone [red line corresponding to $\gamma_q^{x,y} = 0$ in Fig. 2(d)] and is therefore a real solution of Eq. (25). This polaritonic behavior is not described by the result of perturbation theory, which is nevertheless very close to the fully-retarded result aside from the region around the anticrossing with the light line, which appears at a wave number $q \simeq k_0$.

By comparing the above quantum results to the classical ones (green solid lines in Fig. 2), one notices that the redshift in the dispersion [Figs. 2(a) and 2(b)] is stronger (1–3 %) in the case of the classical computation. A significant differ-

<footnote n="4">⁴ Explicitly, with the decomposition of $\Omega_q^\sigma$ into real and imaginary parts mentioned above, the polaritonic annihilation operator (17) takes the form $\mu_q^\sigma(t) = \mu_q^\sigma e^{-i\operatorname{Re}\{\Omega_q^\sigma\}t} e^{-\gamma_q^\sigma t/2}$, such that the decay rate of the collective excitation corresponds to $\gamma_q^\sigma$.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 2. (a), (b) Band structure and (c), (d) radiative decay rates in units of the bare frequency $\omega_0$ and as a function of the reduced wave number $qd$ in half of the Brillouin zone. The left (right) panels represent longitudinal (transverse) polarizations. The red curves correspond to the results obtained within our fully retarded quantum formalism while the green ones are the results obtained within the classical model presented in Appendix A. The blue dashed lines represent the second-order perturbation theory and Fermi's golden rule analytical results. In the upper panels, we also display the quasistatic band structure (21) by solid black lines and the first ($l=0$) light line, $c|q|$, as a dotted black line. The center-to-center distance is $d = 3a$ and $k_0a = 0.3$.</figcaption></figure>

ence between the two models is present for the radiative decay rate only, where the classical model predicts a rate up to 25% lower than the quantum model. However, if the relative difference is noticeable, the absolute one is very small with a maximum deviation of $0.03 \omega_0$.

As experienced by other authors [20, 22], we do not find data points of the retarded band structure for the upper band after the anticrossing with the light line, in contrast with the usual picture of the avoided crossing where it extends upward. It has been argued that the large losses in this region prevent one from finding well-defined roots [20]. We also notice that in the transverse case the real roots, corresponding to the lower band, are difficult to find for values of $qd \ll 1$. The problem is even more stringent using the classical model, where no solutions can be found [20, 30].

Finally, we remark that for values of $qd \ge 3$, the fully retarded and quasistatic dispersion relations become indistinguishable. This is due to the ultraviolet cutoff $\omega_c$ we chose, which limits our theory to wave numbers $|q - 2\pi l/d| \le k_c = 1/a$ and is consistent with the fact that at the edge of the Brillouin zone, the plasmonic excitations are interacting with photonic modes with energy three times larger, which implies

weak coupling. Within this choice of cutoff, only the first light line $l=0$ is of sufficiently low energy to interact with the plasmonic excitations when the center-to-center distance is set to a value of $d$ such that $d \le \pi a$. This explains why the perturbative study of the case $d = 3a$ can be carried out without taking into account the umklapp processes [29].

## 2. Far-field coupled nanoparticles

Considering now a larger center-to-center distance $d = 13a$, far-field effects take place and umklapp processes become of great importance. We show in Fig. 3 the same quantities as previously, i.e., the dispersion relation for the longitudinal [Fig. 3(a)] and transverse [Fig. 3(b)] polarizations, and the corresponding radiative linewidths in Fig. 3(c) and Fig. 3(d), respectively. We see that the same conclusions as for the case $d = 3a$ of Fig. 2 can be drawn from the comparison between the different models, namely, that the perturbation theory reproduces, at least qualitatively in the transverse direction, the exact quantum result for the polaritonic band structure, which itself is close to the results obtained through a purely classical
<figure><img src="image_4.png" /><figcaption>FIG. 3. Same quantities as shown in Fig. 2, except for the black dotted line in (a) and (b), which is now the light line corresponding to the photonic band index $l = 1$, i.e., $c|q - 2\pi/d|$. In the figure, $k_0a = 0.3$ and $d = 13a$.</figcaption></figure>

cal model. As before, however, the perturbation theory misses the avoided crossing feature of the transverse polarization and hence is not reliable around the anticrossing. We also note that since the dipolar coupling strength $\Omega$ scales as $(a/d)^3$ [cf. Eq. (3)], the bandwidth is significantly reduced when we increase the spacing between the nanoparticles. This leads the quasistatic dispersion (21) to being reduced to almost $\omega_0$ in the entire Brillouin zone, so we no longer show it in Figs. 3(a) and 3(b).

What is also immediately remarkable by comparing Figs. 2 and 3 is the flipping of the dispersion curves, the slope of the longitudinal dispersion curve being now negative, as well as the upper band in the transverse case. This is due to the fact that here, the longitudinal (transverse) dispersion intersects (anticrosses) the light line corresponding to the photonic band index $l = 1$, i.e., the curve $c|q - 2\pi/d|$, which is also decreasing in the positive half of the first Brillouin zone, as $qd$ increases. The slope of the closest light line hence dictates the one of the longitudinal dispersion, and as we increase the center-to-center distance $d$, the curve continuously evolves from an increasing [Fig. 2(a)] to a decreasing [Fig. 3(a)] function, and similarly for the transverse polarization [Figs. 2(b) and 3(b)]. We note that we still observe a discontinuity in the longitudinal dispersion at the intersection with the light line, and that in the transverse case, no roots can be found along a

light line $l \neq 0$, neither upward nor downward, for both the quantum or classical approaches.

Since all the dispersion curve is included into the first light cone, i.e., is above the first light line $c|q|$, all the modes are now radiating, as can be seen in Figs. 3(c) and 3(d). This is in stark contrast with the case $d = 3a$ where guided modes immune to radiation damping are present [Figs. 2(c) and 2(d)]. The latter guided modes disappear completely, starting from a center-to-center distance of around $d \approx 11a$, for longitudinal (transverse) polarization and in both classical and quantum models, when the dispersion curve intersects (anticrosses) the first light line $c|q|$ at the edge of the Brillouin zone. Interestingly, this means that for a specific distance $d \approx 11a$, the dispersion curve will intersect both the first ($l = 0$) and second ($l = 1$) light lines. For this precise interparticle distance, one observes a double avoided-crossing scheme in the transverse case, leading to three distinct eigenfrequencies for a given mode $q$.

The second light line $l = 1$ also forms a light cone, containing the modes $q$ whose frequencies are larger than $c|q - 2\pi/d|$, that is to say, the modes on the right of the black dotted line in Figs. 3(a) and 3(b). For the radiative linewidth in the longitudinal case [Fig. 3(c)], one can observe a cusp at the point where the dispersion intersects the second light line, followed by a slight increase of the linewidth. More interestingly, in the
<figure><img src="image_5.png" /><figcaption>FIG. 4. Same quantities as shown in Figs. 2 and 3, except for the black dotted line in panels (a) and (b), which is now the light line corresponding to the photonic band index $l = -1$, i.e., $c|q + 2\pi/d|$. In the figure, $k_0a = 0.3$ and $d = 23a$.</figcaption></figure>

transverse case [Fig. 3(d)], we see that the effect of the second light cone is to enhance the radiative linewidth. The two bands hence radiate very differently, the one inside the second light cone having a radiative linewidth almost four times larger than the one outside. We note that for the modes around the anticrossing, the perturbation theory misses the decrease (increase) of the radiative linewidth of the band outside (inside) the second light cone, and approximates the behavior to an almost step-like function.

Finally, we consider the case of $d = 23a$ in Fig. 4. Now the intersecting light line is the one with band index $l = -1$, i.e., the curve $c|q + 2\pi/d|$ [see the black dotted lines in Figs. 4(a) and 4(b)]. Therefore, the dispersion curves are once again flipped and we again see a different radiating behavior [Figs. 4(c) and 4(d)], depending on whether the polaritonic band is inside or outside the light cone formed by the third photonic band $l = -1$.

We also notice from Fig. 4 that both the dispersion relations and radiative decay rates tend to become flatter. This is consistent with the fact that when the center-to-center distance tends to infinity, the nanoparticles become isolated from each other and there is no more collective effects along the chain. Therefore, all the eigenfrequencies must become degenerate to the value corresponding to the excitation of a single nanoparticle. We will come back in detail to this limiting case in Sec. IV D.

## B. Group velocity of the polaritonic modes

To insight into the propagation of the polaritonic modes, we define their group velocity as the derivative of the dispersion, $v_q^\sigma = \partial\text{Re}\{\Omega_q^\sigma\}/\partial q$. Figure 5 shows the results obtained using a center-to-center distance $d = 3a$, i.e., by taking the derivatives of the curves displayed in Figs. 2(a) and 2(b). We first note that as for the dispersion relations, the group velocities are in very good agreement between the classical and quantum models, for both polarizations.

Unlike what could be expected from its dispersion, the longitudinal case admits quite large group velocities, with a fair amount of modes propagating at almost $0.1c$ or faster, providing a notable window for propagation just after the intersection with the light line, where guided modes present no radiative damping and a rather large group velocity.

At the precise intersection with the light line, as discussed before, no definite value is found in our model free of Ohmic losses, and we observe huge variations of the derivative around the intersection [see Fig. 5(a)], implying difficulties to obtain a good numerical accuracy for values of $v_q^z > 0.2c$. Taking into account Ohmic losses, however, removes the discontinuity and leads to a group velocity of exactly $c$ at the intersection, as demonstrated in Refs. [30, 61] using two different approaches. The analytical second-order
<figure><img src="image_6.png" /><figcaption>FIG. 5. Group velocities $v_q^\sigma$ in units of the speed of light in vacuum $c$, for the (a) longitudinal and (b) transverse polarizations as a function of the reduced wave number $qd$ in half of the first Brillouin zone. The parameters used in the figure are the same as in Fig. 2. The inset shows the longitudinal group velocity in units of $c$ corresponding to the reduced wave number $qd = \pi/2$ as a function of the ratio $d/a$, computed using the analytical perturbative results.</figcaption></figure>

computation only qualitatively accounts for this effect with the appearance of a logarithmic singularity in the group velocity at $\omega_q^\sigma = c|q - 2\pi l/d|$, as can be seen taking the derivative of Eq. (29). This unphysical infinite group velocity at the crossing with the light line shows the limits of the analytical computation in the longitudinal case.

For the transversely polarized excitations [see Fig. 5(b)], we see that the upper band in the dispersion, i.e., the radiating one, shows a rather large negative group velocity with a maximum of about $-0.15c$. The lower band, i.e., the guided modes, which follows the photonic dispersion $c|q|$ inside the light cone [see Fig. 2(b)], propagates at the speed of light $c$ and then rapidly decreases around the intersection with the light line, where a very slow but notably slightly negative group velocity is found. The perturbative treatment, once again, reaches its limitations and presents an unphysical singularity at the crossing.

When we increase the center-to-center distance $d$, as observed in the previous subsection, the dispersion curves maintain a similar shape but get flatter, hence the group velocity globally decreases. However, the longitudinal modes will keep a large group velocity around the intersection with the

light line, which has been found to be of exactly $c$ at the precise intersection when Ohmic losses are taken into account [30, 61]. Interestingly, as discussed previously, increasing the distance between the nanoparticles also leads to a continuous reversing of the dispersion curves, implying an inversion of the sign of the group velocities.

In the inset of Fig. 5, we show the longitudinal group velocity computed from the perturbative results and associated to the mode with wave number $q = \pi/2d$, in units of the speed of light $c$ and as a function of the reduced center-to-center distance $d/a$. Remarkably, we observe an almost periodic sign change of the group velocity. One can show that when we increase the ratio $d/a$, the pseudoperiod of the sign change tends very quickly to $\pi/k_0 a$ in units of $d/a$, i.e., to $\sim 10.5$ within our choice of $k_0 a = 0.3$. In units of $d$, this amounts to an approximate period of $\lambda_0/2$, where $\lambda_0 = 2\pi/k_0$ is the resonance wavelength associated with the isolated single nanoparticle. In the limit of $d \gg a$, the group velocity is hence almost periodic with a period $\lambda_0$. This occurs for all modes $q$, except for $q = 0$ and $q = \pm \pi/d$, respectively the center and edges of the Brillouin zone, for which the group velocity remains zero regardless of the interparticle distance $d$. Interestingly, even if the first sign change appears at a different distance $d$ depending on the mode $q$, when $d \gg a$ the interparticle distance $r_n$ for which a sign change appears becomes approximately $r_n = n\lambda_0/2$, where $n$ is a positive integer. More precisely, the roots of the longitudinal group velocity are given by $r_{q,n}(d) = n\lambda_0/2 + s_q(d)$, where $s_q(d) \ll \lambda_0/2$ is a rapidly decreasing shift that tends to 0. We also note that for the particular mode $q = \pi/2d$ shown in the inset, the absolute value of the group velocity is maximal (and could be equal to $\pm c$ in a model where Ohmic losses would be taken into account [30, 61]) exactly for $r_n = (2n + 1)\lambda_0/2 \mp \lambda_0/4$.

Importantly, the above-discussed unusual behavior of the group velocity allows for a high tunability, since all given modes $q$ except the ones at the center and edges of the Brillouin zone are associated with a wide range of different group velocities, including negative and vanishing ones, depending on the center-to-center distance chosen experimentally. This phenomenon can be compared with the behavior of surface lattice resonances, which admit a strong narrowing of the plasmon resonance due to diffraction and interference effects when the distance between nanoparticle equals the resonance wavelength [62, 63].

### C. Hybridization of light and matter degrees of freedom

As previously discussed, the plasmonic excitations hybridize with the photonic modes and form plasmon-polaritons, which are directly apparent through the avoided-crossing feature of the transverse dispersion relations. Our quantum theory gives access to the share of plasmonic and photonic excitations within a given collective mode with wave number $q$. This distribution can be found through the Hopfield coefficients introduced in Eq. (15), where each coefficient represents the weighting of the plasmonic and photonic excitations.
<figure><img src="image_7.png" /><figcaption>FIG. 6. Modulus squared of the plasmonic Hopfield coefficients $W_q^\sigma$ (orange solid line) and sum of all the moduli squared of the photonic ones $Y_{\kappa q}^{l,\sigma,\hat{\lambda}_{\kappa q}^l}$ (purple solid line), as a function of the reduced wave number $qd$ in half of the Brillouin zone. The vertical dashed line represents the mode for which the polaritonic dispersion crosses a light line. The upper and lower panels show, respectively, the longitudinal and transverse polarizations, while the three columns are associated, respectively, to center-to-center distances $d = 3a$, $d = 13a$, and $d = 23a$. In the figure, $k_0a = 0.3$.</figcaption></figure>

Using the normalization condition (16), one can solve the system of equations (19), respectively, for the moduli squared of the two plasmonic Hopfield coefficients and for the sum over all band indexes $l$, wave vectors $\kappa$ and polarizations $\hat{\lambda}_{\kappa q}^l$ of the moduli squared of the two photonic coefficients. Transforming the sums over $\kappa$ into finite integrals, one obtains functions of the complex eigenfrequencies $\Omega_q^\sigma$. We give details on the computation and provide the analytical expressions of the Hopfield coefficients in Appendix C.

Using the eigenfrequencies found via Eq. (25), we show in Fig. 6 the results for the longitudinal (upper panels) and transverse (lower panels) polarizations, for center-to-center distances $d = 3a$, $d = 13a$, and $d = 23a$. We note that the two Hopfield coefficients $X_q^\sigma$ and $Z_{\kappa q}^{l,\sigma,\hat{\lambda}_{\kappa q}^l}$, associated with the counter-rotating terms in the polaritonic ladder operator (15) are not plotted in the figure, since their largest contributions are less than 0.006. This suggests that we could have employed the rotating wave approximation (RWA) in our derivations. However, since the RWA does not lead to any computational simplification, we chose to use the full Hamiltonian (9).

Remarkably, we notice that for the longitudinal polarization [Figs. 6(a), 6(c), and 6(e)], and hence without having any avoided-crossing dispersion, the normal modes inside a light cone are an equitable mix of light and matter degrees of freedom, proving that we do indeed have plasmon–polaritons, even in this case. In the near-field case [Fig. 6(a), $d = 3a$] at the end of the light cone, the excitation becomes almost purely

plasmonic, which is consistent with the suppression of radiative losses. For center-to-center distances such that the dispersion curve intersects a photonic band of index $l \neq 0$ [Figs. 6(c) and 6(e)], the entire Brillouin zone is located within the light cone $l=0$ and all modes are radiating. Thus in this case all modes are polaritonic. Approaching the mode for which the dispersion curve crosses the light line, represented as a vertical dashed line in the figure, the excitation becomes mostly photonic. At the exact crossing, the Hopfield coefficients are not rigorously defined due to the discontinuity and avoided crossing in the dispersions, but one can argue that at this precise point, the excitation tends to be purely photonic, consistent with the fact that the associated group velocity tends to $\pm c$.

For the transverse polarization, let us first discuss the near-field case, Fig. 6(b). The upper plasmon–polariton band, associated with a nonzero radiative damping, presents a strong part of photon excitation which, quite surprisingly, decreases when approaching the intersection with the light line. We note that the polaritonic nature is marked for all modes inside the light cone, and not only those close to the intersection. This can be understood by the fact that plasmons interact with all photonic modes $k$ and not only with the light line, which corresponds to $k = (0, 0, q)$. We further notice that the lower band of the transverse dispersion is almost exclusively composed of photons inside the light cone, which is in adequacy with the fact that it is not subject to radiative losses and has a dispersion relation almost equal to the light line. Outside
the light cone, the modes are nearly purely plasmonic, with a well-balanced mix only around the intersection.

In the far field, when all modes are radiating [Figs. 6(d) and 6(f)], we observe an inversion in the predominance between photonic and plasmonic excitations at the crossing with a light line. The band which radiates the most, i.e., the one located inside the light cone formed by light line *l* that intersects the dispersion curve, is mainly plasmonic. The other band, which in the near-field case corresponded to the nonradiating one, shows an equitable share between photonic and plasmonic excitations away from a light line. When it approaches the intersection, it becomes predominantly photonic, as it can be inferred from the dispersion curves where this band falls into the light line.

Finally, we remark that for both polarizations, increasing the center-to-center distance leads to an increasingly equitable mix between photonic and plasmonic excitations, with both contribution tending to 0.5. This is consistent with the fact that, as will be discussed in detail in the following subsection, when $d \gg a$, we end up with single noninteracting nanoparticles but which interact with the photonic continuum, leading to a polaritonic eigenstate with equal plasmonic and photonic weightings.

### D. Single nanoparticle limit

As mentioned above, we can use the limit of infinitely spaced nanoparticles, $d/a \to \infty$, to study the case of isolated single nanoparticles. These noninteracting nanoparticles are still coupled to the surrounding photonic continuum, leading to a radiative shift $\delta_0$ of the resonance frequency $\omega_0$, and to a radiative decay rate $\gamma_0$.

As discussed previously, when the nanoparticles are moved apart in the chain, the dispersion curves flatten and all the eigenfrequencies become degenerate at the eigenfrequency of a single nanoparticle. All the quantities hence must become polarization and mode-independent and the analysis of the Hopfield coefficients, in the previous subsection, showed that one tends, indeed, toward a single degenerate mode of polaritonic nature, consisting of an equal mix of light and matter excitations.

The single nanoparticle case has already been studied perturbatively in Ref. [46], using a similar quantum formalism as in the present paper. Using second-order perturbation theory (as was done for the chain, cf. Sec. III B) the radiative shift $\delta_0$ is given by $\delta_0 = (E_{n^\sigma+1}^{(2)} - E_{n^\sigma}^{(2)})/\hbar$, where the second-order energy correction reads

$$E_{n^\sigma}^{(2)} = \pi \hbar \omega_0^3 \frac{a^3}{V} \sum_{\mathbf{k}, \hat{\lambda}_{\mathbf{k}}} \frac{(\hat{\sigma} \cdot \hat{\lambda}_{\mathbf{k}})^2}{\nu_{\mathbf{k}}} \frac{(2n^\sigma + 1)\nu_{\mathbf{k}} - \omega_0}{\omega_0^2 - \nu_{\mathbf{k}}^2}. \quad (32)$$

In the above expression, $n^\sigma$ is the quantum number corresponding to the number of plasmonic excitations with polarization $\sigma$, and $\nu_k = \nu_k^{l=0}$ is the photon dispersion, where, since the single nanoparticle is not anymore a periodic system, we do not have to take into account umklapp plasmon-photon processes.

The correction (32) appears to be linearly divergent. The authors of Ref. [46] hence regularized such a divergence using a renormalization scheme analogous to that used by Bethe in his analysis of the Lamb shift in atomic physics [64]. Bethe's mass renormalization amounts to subtract from Eq. (32) the energy shift corresponding to free electrons coupled to the photonic continuum, which can be obtained taking the limit of vanishing transition frequency $\omega_0 \to 0$. Notably, this leads to a divergence which is only logarithmic.

However, as detailed in Sec. III B in the case of the chain, the second-order correction is logarithmically divergent. Trying to use the same renormalization procedure for the chain leads to unphysical results, with a singularity at the center of the Brillouin zone. Furthermore, a perturbative treatment of the light-matter interaction in a similar plasmonic model but in two-dimensional lattices with arbitrary geometries leads to finite, nondivergent corrections, which agree well with classical electrodynamic calculations [39]. This leads us to question the validity of Bethe's mass renormalization procedure in our specific case of a plasmonic system, for which the energy scales are quite different from the case of atomic physics studied originally by Bethe. In Appendix D, we also discuss the results obtained through this mass renormalization procedure for the case of nanoparticle dimers.

Keeping the linearly divergent second-order correction (32), one can go to the continuum limit and transform the sum over wave vectors $\mathbf{k}$ into principal value integrals using spherical coordinates. After calculating the integrals, the single nanoparticle radiative frequency shift reads

$$\delta_0 = \frac{\omega_0}{3\pi} (k_0 a)^3 \left[ \ln \left( \frac{\omega_c / \omega_0 + 1}{\omega_c / \omega_0 - 1} \right) - 2 \frac{\omega_c}{\omega_0} \right], \quad (33)$$

where, as expected, the expression now depends linearly on the cutoff $\omega_c$ and hence also on the cutoff parameter $\Lambda$ (see footnote 3).

The upper panel of Fig. 7 shows the different results for the radiative shift $\delta_0$, according to the hypothesis used. We see that the results from Ref. [46] are significantly different both from those obtained by the classical model (green solid line) and by the quantum formalism where Bethe's mass renormalization procedure was not used (red solid line and blue dashed line) for the chain via the limit $d/a \to \infty$, and from the quantum model of a single nanoparticle via Eq. (33). Reference [46] predicts a minute blueshift, while other methods predict a redshift which is one order of magnitude larger (in absolute value). We see that, as expected, the results from the chain using our quantum formalism are in perfect agreement with the one from the single nanoparticle computation (33), and the three curves are almost perfectly overlapping. The classical computation also leads to a redshift of the same order of magnitude but larger. We note that in the classical formalism, using the radiative correction to the classical polarizability (which is defined in Appendix A) leads to almost no frequency shift, while using the modified long wavelength approximation (MLWA, see Appendix A) leads to results which are almost the same as the one presented here for which we used the exact Mie polarizability.

The radiative decay rate $\gamma_0$ of the single nanoparticle can
<figure><img src="image_8.png" /><figcaption>FIG. 7. Single nanoparticle radiative shift $\delta_0$ (upper panel) and radiative decay rate $\gamma_0$ (lower panel) in units of the single nanoparticle resonance frequency $\omega_0$ and as a function of the reduced nanoparticle radius $k_0a$. The green solid line represents the result obtained from the classical model, as presented in Appendix A, while the thick grey solid line (almost perfectly overlapped by both the fully retarded quantum results in red and the perturbative results in dashed blue) represents the single nanoparticle results of Eqs. (33) and (34). In the upper panel, the orange line represents the result for the radiative shift presented in Ref. [46], and the grey dotted line is a guide for the eye.</figcaption></figure>

also be obtained perturbatively through our quantum formalism. Using Fermi's golden rule similarly as for the chain in Sec. III B, one obtains [28]

$$ \gamma_0 = \frac{2\omega_0^4 a^3}{3c^3}. \qquad (34) $$

This familiar expression, representing the rate at which the plasmons dissipate their energy by emitting photons, can be obtained classically through the direct evaluation of the radiated power of an oscillating electric dipole [65].

In the lower panel of Fig. 7, we show the comparison between Eq. (34) and the results obtained by taking the limit of infinitely spaced nanoparticles in the chain. We see that both the fully retarded and perturbative quantum results are in perfect agreement with the expression (34), as the three curves are almost indistinguishable. This indicates that taking the limit of infinite center-to-center distance $d$ is an accurate method to study the case of single nanoparticles.

We also observe in the lower panel of Fig. 7 that the classical chain gives slightly different results and predicts a smaller damping rate for values of $k_0a \gtrsim 0.2$. Notably, this differ-

ence depends on the approximation we adopt for the classical polarizability of the nanoparticle. Indeed, using the radiative correction gives results which agree with Eq. (34), while using the MLWA, which takes into account possible depolarization effects [see the discussion in Appendix A after Eq. (A3)] gives results similar to the one presented in Fig. 7 (see green solid line), where we used the Mie polarizability.<fnref n="5" /> The small difference in the obtained results for the single nanoparticle radiative decay rate are thus due to depolarization effects, which increase with the (reduced) nanoparticle radius $k_0a$. However, this cannot fully explain the small discrepancy in the decay rates of the chain discussed in Sec. IV A, as the classical approach slightly underestimates the rate compared to the quantum one, regardless of the approximation performed for the classical polarizability.

As can be seen from the lower panel of Fig. 7, for a reduced radius $k_0a = 0.3$, one obtains $\gamma_0 = 0.018 \omega_0$. This means that in the case of the chain, several modes are superradiant. In particular, in the near-field case presented in Fig. 2, all transverse radiating modes are superradiant, with a damping from 3 to 9 times larger than that of the single nanoparticle. Interestingly, when we increase the center-to-center distance $d$ in Figs. 3 and 4, only the modes inside the light cone formed by the light line $l$ which intersects the dispersion curve (respectively, $l = +1$ and $l = -1$) present superradiance, while the ones outside this light cone are subradiant [see Figs. 3(c), 3(d), 4(c), and 4(d)].

V. CONCLUSION

We have provided a quantum theory of a linear chain of metallic nanoparticles hosting plasmon-polaritons. In particular, the retardation effects of the long-range dipole-dipole interaction along the chain have been considered exactly. These retardation effects, arising from the light-matter interaction between the plasmonic degrees of freedom and a three-dimensional photonic continuum, lead to a polaritonic behavior of the eigenmodes of the system. By taking into account plasmon-photon umklapp processes, we have been able to describe nanoparticles at arbitrary center-to-center distances $d \gtrsim 3a$, that is, in both near- and far-field coupling regimes.

We have shown, by computing the band structure, the radiative linewidth, as well as the group velocities of the eigenmodes, that our quantum description is in good agreement with the usual classical approach used for plasmonic chains. The latter approach has also been shown to be in satisfying agreement with simulations using finite-difference time-domain techniques [66]. We have provided approximate analytical expressions for the above-mentioned quantities which reproduce the exact results, at least at a qualitative level, although they do not describe the polaritonic nature of the excitations. These analytical expressions are valid in both near-

<footnote n="5">⁵ Note that the MLWA can be obtained from the Mie polarizability (A3) through an expansion in the small parameter $\omega a/c$.</footnote>
and far-field situations and allow for a much more transparent formalism than in the classical case, where a numerical resolution is needed.

In the far-field regime and thanks to our analytical developments, we have observed a particularly unusual behavior of the group velocity, which changes signs and vanishes every time the center-to-center distance *d* approximately equals a positive multiple of half of the resonance wavelength associated with the single isolated nanoparticle. The group velocity of any polaritonic mode can hence be tuned to a variety of different values, including zero and large positive or negative group velocities, depending on the distance *d*.

We have studied quantitatively the hybridization of plasmonic and photonic degrees of freedom using Hopfiled's coefficients, showing a well-balanced hybridization inside the light cone even in the longitudinal polarization case, which leads us to associate the presence of plasmon-polariton excitations to the presence of radiative losses rather than to the presence of an avoided crossing dispersion.

Finally, we have studied the regime of infinitely spaced nanoparticles and found eigenfrequencies and radiative decay rates in perfect agreement with the single nanoparticle case. Notably, we have found results that are qualitatively different from the ones of Ref. [46]. We showed that these differences arise from a mass renormalization scheme that may not be adequate to the present situation and that we did not follow here.

Our detailed analysis supports the development of a quantum formalism, which is required to have a self-contained theory of arrays of interacting nanoparticles, where quantum effects, such as Landau damping or electronic spill-out, can be readily incorporated. Our work also opens the perspective of studying other systems where light-matter interaction is of major importance, such as cavity-embedded systems, arrays of larger dimensionality and/or with other geometries, and plasmonic systems composed of nonspherical nanoparticles.

## ACKNOWLEDGMENTS

We thank Charles A. Downing, Charlie-Ray Mann, and Eros Mariani for insightful discussions. We acknowledge financial support from the French National Research Agency (ANR) through Grant No. ANR-14-CE26-0005 Q-MetaMat. This work of the Interdisciplinary Thematic Institute QMat, as part of the ITI 2021-2028 program of the University of Strasbourg, CNRS, and Inserm, was supported by IdEx Unistra (ANR 10 IDEX 0002), and by SFRI STRAT'US Projects No. ANR-20-SFRI-0012 and No. ANR-17-EURE-0024 under the framework of the French Investments for the Future Program.

## Appendix A: Classical electromagnetic modeling of the nanoparticle chain

The standard theoretical approach for describing the plasmonic chain introduced in the main text (see Fig. 1) is based on macroscopic classical electrodynamics and benefits from extensive literature [3, 15–17, 19–25], which is still active

nowadays [30, 31]. In this appendix, we describe the dielectric properties of the metallic nanoparticles using the Drude model. Neglecting Ohmic losses, the nanoparticle local permittivity is given by $\epsilon(\omega) = 1 - \omega_p^2/\omega^2$.

We first start by looking for the polarizability of a single nanoparticle. Using Mie's theory [58, 59] and considering only the electric dipole term, one obtains

$$ \mathbf{E}_s = \frac{3ic^3 E_0 a_1(\omega)}{2\omega^3} \frac{e^{i\omega r/c}}{r} \left\{ \frac{\omega^2}{c^2}(\hat{r} \times \hat{u}) \times \hat{r} + [3\hat{r}(\hat{r} \cdot \hat{u}) - \hat{u}] \left( \frac{1}{r^2} - \frac{i\omega}{cr} \right) \right\} \quad (A1) $$

for the expression of the electric field scattered by the sphere. Here, $E_0$ is the amplitude of the incident field $\mathbf{E}_i = E_0 e^{i\omega r \cos\theta/c} \hat{u}$ in a given direction $\hat{u}$, and

$$ a_1(\omega) = \frac{n\psi_1(nka)\psi_1'(ka) - \psi_1(ka)\psi_1'(nka)}{n\psi_1(nka)\xi_1'(ka) - \xi_1(ka)\psi_1'(nka)} \quad (A2) $$

is the first (electric dipole) Mie coefficient [60], where $n = \sqrt{\epsilon(\omega)}$ is the refractive index, $k = \omega/c$ the wave number in the metal, and $a$ the nanoparticle radius. $\psi_1(x) = \sin(x)/x - \cos(x)$ and $\xi_1(x) = (-i/x - 1)e^{ix}$ are, respectively, the first Ricatti-Bessel functions of the first and second kinds. The scattered field (A1) corresponds to the radiation field emitted by an oscillating electric dipole, therefore taking fully into account retardation effects, with a dipolar moment $\mathbf{p} = (3i/2)(c/\omega)^3 E_0 a_1(\omega) \hat{u}$. Using that the dipolar moment is also equal to $\mathbf{p} = \alpha(\omega)\mathbf{E}_i|_{r=0}$ leads to an expression of the dynamic, frequency-dependent polarizability of the metallic nanoparticle [22, 25, 67],

$$ \alpha(\omega) = \left(-i \frac{2\omega^3}{3c^3}\right)^{-1} a_1(\omega). \quad (A3) $$

Many authors [16, 20, 21, 24, 31] use an *ad hoc* correction to the quasistatic polarizability of a spherical nanoparticle

$$ \alpha_{qs}(\omega) = a^3 \frac{\epsilon(\omega) - 1}{\epsilon(\omega) + 2}, \quad (A4) $$

instead of the exact Mie expression (A3). To take into account the radiative decay of the plasmons and preserve energy conservation, the choice $\alpha_{adhoc}^{-1} = \alpha_{qs}^{-1} - i2\omega^3/3c^3$ is usually made. Such a radiative correction can be found by adding a radiation reaction field to the quasistatic system [65, 68, 69]. Other authors [30] also include a contribution proportional to $\omega^2$ to the inverse polarizability, stemming from the introduction of a depolarization field in the system. This is known as the modified long wavelength approximation (MLWA), and it leads to a redshift of the polaritonic dispersion curves which has been interpreted as arising from a dephasing between radiation emitted by different parts of the sphere [70].

As already discussed in Refs. [60, 71, 72], an expansion of the exact polarizability (A3) for small particles to third order in $\omega a/c$ leads, up to a numerical factor, to the same *ad hoc* corrections mentioned above.

Since the radiative correction and the MLWA are approximations of Eq. (A3), and since using the latter equation does
not increase substantially the computation time, in this paper we have chosen to describe the nanoparticles within the classical model using the exact Mie polarizability (A3). Following Ref. [20], one can then obtain an implicit dispersion relation by writing the induced dipole moment $\mathbf{p}_n$ of the nanoparticle $n$ as $\mathbf{p}_n = \alpha(\omega)\mathbf{E}_{\text{loc}}$, where $\mathbf{E}_{\text{loc}}$ is the total electric field generated by all the others nanoparticles $m \neq n$. One finds

$$1 + \frac{\alpha(\omega)}{d^3} \Sigma_{\sigma}(\omega, q) = 0, \qquad (A5)$$

where $q$ is the wave number of the normal mode along the chain, belonging to the first Brillouin zone. Equation (A5) can be seen as the classical analog of Eq. (25). As it has been shown in Ref. [19], the function $\Sigma_{\sigma}(\omega, q)$ can be rewritten, using analytical continuation, in terms of a sum of polylogarithms,

$$\Sigma_z(\omega, q) = 2i \frac{\omega d}{c} [\text{Li}_2(\varphi^+) + \text{Li}_2(\varphi^-)] - 2 [\text{Li}_3(\varphi^+) + \text{Li}_3(\varphi^-)], \qquad (A6a)$$

$$\Sigma_{x,y}(\omega, q) = - \left( \frac{\omega d}{c} \right)^2 [\text{Li}_1(\varphi^+) + \text{Li}_1(\varphi^-)] - i \frac{\omega d}{c} [\text{Li}_2(\varphi^+) + \text{Li}_2(\varphi^-)] + [\text{Li}_3(\varphi^+) + \text{Li}_3(\varphi^-)], \qquad (A6b)$$

where $\varphi^{\pm} = e^{i(\omega/c \pm q)d}$. Here, $\sigma = x, y$ or $z$ characterizes the polarization of the collective excitation, depending on whether the induced dipole moment points along the array ($\sigma = z$) or orthogonal to it ($\sigma = x, y$). Since collective and individual LSP polarizations are aligned, it is equivalent to consider a given polarization $\sigma$ on each individual LSP.

Numerically solving Eq. (A5) for complex $\omega$ leads to the dispersion and to the radiative linewidth of the polaritonic normal modes discussed in the main text (see green solid lines in Figs. 2–4).

## Appendix B: The diamagnetic $A^2$ term

In our treatment of the photonic environment, we discarded the $A^2$ term, also known as the diamagnetic term. In this appendix, we discuss the origin of such a term and give justifications for this approximation.

Physically, the diamagnetic term represents a photon self-interaction energy which is not involved in the retardation effects of the light-matter coupling. It comes from the minimal-coupling substitution of the momentum $\mathbf{P} \rightarrow \mathbf{P} + (e/c)\mathbf{A}(\mathbf{r})$ in the kinetic Hamiltonian $H_{\text{kin}} = \mathbf{P}^2/2m$, with $m$ the considered mass. We recall that this coupling originates fundamentally in a relativistic theory by the replacement in the free Lagrangian of the usual derivatives with the covariant derivatives, $\partial_{\mu} \rightarrow D_{\mu} = \partial_{\mu} - (ie/c)A_{\mu}$, where $A_{\mu}$ is the four-vector gauge potential. This replacement is required to preserve the $U(1)$ gauge symmetry of electromagnetism, therefore, the light-matter coupling is fully determined by the gauge invariance of the theory itself [73].

Neglecting the diamagnetic term can thus lead to conceptual problems such as the loss of gauge invariance but also unphysical ground states [74, 75]. It is, however, frequently discarded due to its quadratic nature in the light-matter coupling and several well-known quantum optics models such as the Jaynes-Cummings or Dicke model do not consider it. Perceptible differences due to the diamagnetic term such as gauge ambiguities were found in the context of the ultra strong coupling (USC) or deep strong coupling (DSC) regime only [76].

Including the $A^2$-term in our model amounts to replacing Eq. (6) by the quadratic plasmon-photon coupling Hamiltonian

$$H_{\text{pl-ph}}^{\sigma, A^2} = \frac{e}{m_e c} \sum_{n=1}^{N} \Pi_n^{\sigma} \cdot \mathbf{A}(\mathbf{d}_n) + \frac{N_e e^2}{2m_e c^2} \sum_{n=1}^{N} \mathbf{A}^2(\mathbf{d}_n). \qquad (B1)$$

Since the quadratic $A^2$ term does not involve plasmonic degrees of freedom, it has no effect when one applies perturbation theory on the plasmonic Hamiltonian (2), except for a global energy shift [29]. Therefore, the results found through the perturbative treatment in Sec. III B are the same, whether we include such a term or not. However, it can renormalize the photonic dispersion by inducing what is usually referred to as a diamagnetic shift, hence slightly modifying the polaritonic dispersion.

To gain insight into the effects of such a term in our system, one can diagonalize the Hopfield model, corresponding to a model close to the one studied here but with a three-dimensional dipolar system instead of a one-dimensional one, with and without the quadratic term and show that different behaviors occur only for the USC and DSC regimes [77]. We recall that one enters in the USC regime when the counter-rotating terms of the Hamiltonian become sizable [78]. From the analysis of the Hopfield coefficients in Sec. IV C, one can then deduce that in the model studied here, we are far from the USC regime and therefore that the diamagnetic term can be safely neglected.

## Appendix C: Analytical expressions for the Hopfield coefficients

In Sec. IV C, we discussed the behavior of the Hopfield coefficients, namely the plasmonic and photonic weightings of the polaritonic excitation. Here, we provide the analytical expressions for the latter weightings, which are displayed in Fig. 6. Using the system of equations (19), one can write the photonic part of the weightings of the excitation as

$$\sum_l \sum_{\kappa, \lambda_{\kappa q}} \left| Y_{\kappa q}^{l, \sigma, \lambda_{\kappa q}} \right|^2 = \frac{4\omega_0^2\omega_q^{\sigma 4} |X_q^{\sigma}|^2}{|\omega_q^{\sigma 2} - \omega_0\Omega_q^{\sigma}|^2} \sum_l \mathcal{I}_{q,-l}^{l,\sigma} \qquad (C1)$$

and

$$\sum_l \sum_{\kappa, \lambda_{\kappa q}} \left| Z_{\kappa q}^{l, \sigma, \lambda_{\kappa q}} \right|^2 = \frac{4\omega_0^2\omega_q^{\sigma 4} |X_q^{\sigma}|^2}{|\omega_q^{\sigma 2} - \omega_0\Omega_q^{\sigma}|^2} \sum_l \mathcal{I}_{q,+l}^{l,\sigma} \qquad (C2)$$
and the plasmonic weights as

$$|W_q^\sigma|^2 = |X_q^\sigma|^2 \left| \frac{\omega_q^{\sigma 2} + \omega_0 \Omega_q^\sigma}{\omega_q^{\sigma 2} - \omega_0 \Omega_q^\sigma} \right|^2 \quad (C3)$$

and

$$|X_q^\sigma|^2 = \left( \frac{4\omega_0^2\omega_q^{\sigma 4}}{|\omega_q^{\sigma 2} - \omega_0\Omega_q^\sigma|^2} \sum_l \left[ \mathcal{I}_{q,-}^{l,\sigma} - \mathcal{I}_{q,+}^{l,\sigma} \right] + \left| \frac{\omega_q^{\sigma 2} + \omega_0\Omega_q^\sigma}{\omega_q^{\sigma 2} - \omega_0\Omega_q^\sigma} \right|^2 - 1 \right)^{-1} \quad (C4)$$

In the above expressions, $\mathcal{I}_{q,\pm}^{l,\sigma}$ are functions of the complex eigenfrequencies $\Omega_q^\sigma$, and read

$$\mathcal{I}_{q, \pm}^{l, \sigma} = \sum_{\kappa, \hat{\lambda}_{\kappa q}^l} \frac{(\xi_{\kappa q}^l)^2 (\hat{\sigma} \cdot \hat{\lambda}_{\kappa q}^l)^2}{|\Omega_q^\sigma \pm \nu_{\kappa q}^l|^2} \quad (C5)$$

Summing over the photon polarizations using Eq. (22), and transforming the sum over $\kappa$ to integrals via Eq. (23), leads to the following result:

$$\mathcal{I}_{q, \pm}^{l, \sigma} = \pm \frac{\eta^{\sigma} \omega_0 a^3}{4dc^2} \left| \frac{cq^l}{\Omega_q^{\sigma}} \right|^2 \Theta(\omega_c - c|q^l|) \left\{ \mp \left( \frac{1}{\omega_c} - \frac{1}{c|q^l|} \right) \right. \\ - \frac{\Omega_q^{r, \sigma}}{|\Omega_q^{\sigma}|^2} \ln \left( \frac{(\Omega_q^{i, \sigma})^2 + (\Omega_q^{r, \sigma} \pm c|q^l|)^2}{(\Omega_q^{i, \sigma})^2 + (\Omega_q^{r, \sigma} \pm \omega_c)^2} \frac{\omega_c^2}{(c|q^l|)^2} \right) \\ + \frac{\tan^{-1} \left( \frac{\Omega_q^{r, \sigma} \pm \omega_c}{\Omega_q^{i, \sigma}} \right) - \tan^{-1} \left( \frac{\Omega_q^{r, \sigma} \pm c|q^l|}{\Omega_q^{i, \sigma}} \right)}{\Omega_q^{i, \sigma}} \\ \left. \times \left[ \operatorname{sgn}\{\eta^{\sigma}\} \left| \frac{\Omega_q^{\sigma}}{cq^l} \right|^2 + \frac{(\Omega_q^{r, \sigma})^2 - (\Omega_q^{i, \sigma})^2}{|\Omega_q^{\sigma}|^2} \right] \right\} \quad (C6)$$

where $q^l = q - 2\pi l/d$, and where, for clarity, we have written the real and imaginary parts of the polaritonic eigenfrequencies $\Omega_q^\sigma$ as $\Omega_q^{r,\sigma}$ and $\Omega_q^{i,\sigma}$, respectively. We note that the term in the third line of the above equation is well-defined due to the fact that the imaginary part $\Omega_q^{i,\sigma} \to 0$ only when the real part $\Omega_q^{r,\sigma} < c|q^l|$, i.e., guided modes which are immune to radiation damping exist only outside the first light cone, as discussed in Sec. IV.

## Appendix D: Discussion of Bethe's mass renormalization procedure for nanoparticle dimers

As shown in Sec. IV D, we obtain drastically different results when we do not apply Bethe's mass renormalization procedure to our quantum model of a single metallic nanoparticle, and this led us to question the use of Bethe's approach in the specific case of a plasmonic system. The radiative shift induced by the vacuum electromagnetic modes on a dimer of interacting metallic nanoparticles was computed in Ref. [46]

using Bethe's scheme. In this appendix, we will reproduce the same quantities as presented in Ref. [46], but without using the mass renormalization procedure, and comment on the differences obtained.

The dimer of interacting metallic nanoparticles is modeled by the plasmonic Hamiltonian of the chain (2) for which we set the number of nanoparticles $N = 2$. The dimer interacts with a photonic environment, and since the system is not periodic, we do not have to take into account umklapp plasmon-photon processes. The Hamiltonian describing the photonic environment is thus given by Eq. (4) for which we only keep the term with $l = 0$. The coupling Hamiltonian is therefore given by Eq. (6) with $N = 2$ and $l = 0$.

As detailed in Ref. [46], such a system hosts bright and dark hybridized collective plasmonic modes, and one can easily diagonalize the plasmonic Hamiltonian to find the bare (quasistatic) eigenfrequencies,

$$\omega_{\tau}^{\sigma} = \omega_{0} \sqrt{1 + 2\tau |\eta^{\sigma}| \frac{\Omega}{\omega_{0}}}. \quad (D1)$$

The label $\tau = \pm$ distinguishes the high- and low-energy coupled plasmonic modes. Importantly, the high-energy transverse (↑↓) and low-energy longitudinal (→→) excitations correspond to symmetric bright modes coupled to the photonic continuum, while the low-energy transverse (↑↓) and high-energy longitudinal (→←) excitations are antisymmetric dark modes, hence weakly coupled to the electromagnetic continuum. Just as in the single-nanoparticle case or for the chain, these plasmonic modes will hybridize with the electromagnetic continuum, leading to dressed eigenfrequencies $\tilde{\omega}_{\tau}^{\sigma} = \omega_{\tau}^{\sigma} + \delta_{\tau}^{\sigma}$, where $\delta_{\tau}^{\sigma}$ are the radiative frequency shifts.

To compute the dressed eigenfrequencies, we rely on Ref. [46], which uses a similar methodology as the one we detailed in Sec. III B. Without making use of the mass renormalization procedure, we obtain $\delta_{\tau}^{\sigma} = (E_{n_{\tau}^{\sigma}+1}^{(2)} - E_{n_{\tau}^{\sigma}}^{(2)})/\hbar$, where the second-order correction reads

$$E_{n_{\tau}^{\sigma}}^{(2)} = \pi \hbar \omega_0^2 \omega_{\tau}^{\sigma} \frac{a^3}{V} \sum_{\mathbf{k}, \hat{\mathbf{k}}} \frac{(\hat{\sigma} \cdot \hat{\mathbf{k}})^2 (2n_{\tau}^{\sigma} + 1) \nu_{\mathbf{k}} - \omega_{\tau}^{\sigma}}{\nu_{\mathbf{k}} \omega_{\tau}^{\sigma 2} - \nu_{\mathbf{k}}^2} \times [1 + \tau \operatorname{sgn}\{\eta^{\sigma}\} \cos(k_z d)]. \quad (D2)$$

Computing the summation over photon polarizations through the relation (22) and transforming the sum over wave vectors $k$ into a principal value integral in spherical coordinates, we get

$$\begin{aligned} \delta_{\tau}^{\sigma} = 2\pi\omega_0^2\omega_{\tau}^{\sigma} \frac{a^3}{8\pi^3} \mathcal{P} \int_0^{k_c} \frac{k^2 dk}{\omega_{\tau}^{\sigma 2} - c^2 k^2} \\ \times \int_0^{\pi} d\theta \sin\theta [1 + \tau \operatorname{sgn}\{\eta^{\sigma}\} \cos(kd \cos\theta)] \\ \times \int_0^{2\pi} d\varphi [1 - (\hat{\sigma} \cdot \hat{k})^2]. \end{aligned} \quad (D3)$$

After a long but straightforward calculation, we obtain the ra-
<figure><img src="image_9.png" /><figcaption>FIG. 8. Radiative frequency shifts in units of the bare frequency $\omega_0$ as a function of the reduced nanoparticle radius $k_0a$. The same approach as for the chain and two-dimensional systems has been used, i.e., without making use of Bethe's mass renormalization. The grey dotted line is a guide for the eye. The interparticle distance is $d = 3a$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 9. Ratio $\Delta\tilde{\omega}^z/\Delta\tilde{\omega}^{x,y}$ of longitudinal and transverse frequency splittings in a nanoparticle dimer as a function of the (reduced) interparticle distance $k_0d$, for increasing (reduced) nanoparticle radius $k_0a$, from 0.1 (green line), 0.15 (yellow line), to 0.2 (red line). The grey dotted line shows the bare ratio $\Delta\omega^z/\Delta\omega^{x,y}$ in the absence of the coupling with the photonic environment. The UV frequency cutoff is chosen as $\omega_c = c/a$.</figcaption></figure>

radiative shifts

$$ \delta_{\tau}^{\sigma} = \frac{\omega_{\tau}^{\sigma 2}}{3\pi\omega_{0}} (k_{0}a)^{3} \left[ \ln\left(\frac{\omega_{c}/\omega_{\tau}^{\sigma} + 1}{\omega_{c}/\omega_{\tau}^{\sigma} - 1}\right) - 2\frac{\omega_{c}}{\omega_{\tau}^{\sigma}} \right] \\ + \tau|\eta^{\sigma}|\frac{\omega_{\tau}^{\sigma}\Omega}{\pi\omega_{0}}(k_{0}d)^{2}\left[\frac{1 + \operatorname{sgn}\{\eta^{\sigma}\}}{2}\right]g_{\tau}^{\sigma} \\ + \frac{h_{\tau}^{\sigma}}{k_{\tau}^{\sigma}d} - \frac{2\operatorname{Si}(k_{c}d) + g_{\tau}^{\sigma}}{(k_{\tau}^{\sigma}d)^{2}}, \qquad (\text{D4}) $$

where $k_{\tau}^{\sigma} = \omega_{\tau}^{\sigma}/c$,

$$ g_{\tau}^{\sigma} = \sum_{\zeta=\pm} \left[ \zeta \sin(k_{\tau}^{\sigma} d) \mathrm{Ci}(k_{c}d + \zeta k_{\tau}^{\sigma} d) - \cos(k_{\tau}^{\sigma} d) \mathrm{Si}(k_{c}d + \zeta k_{\tau}^{\sigma} d) \right], \qquad (\text{D5}) $$

$$ h_{\tau}^{\sigma} = \sum_{\zeta=\pm} \left[ \zeta \cos(k_{\tau}^{\sigma} d) \mathrm{Ci}(k_{c}d + \zeta k_{\tau}^{\sigma} d) + \sin(k_{\tau}^{\sigma} d) \mathrm{Si}(k_{c}d + \zeta k_{\tau}^{\sigma} d) \right], \qquad (\text{D6}) $$

and with $\mathrm{Si}(x)$ and $\mathrm{Ci}(x)$ denoting the sine and cosine integrals, respectively. As expected and as for the single nanoparticle case [see Eq. (33)], the expression (D4) depends linearly on the cutoff $\omega_c$.

In Fig. 8, we show the resulting radiative shifts for the dimer, together with the shift (33) obtained for the single nanoparticle. In Ref. [46], a blueshift has been obtained for the bright modes, whereas only the dark modes were redshifted. This resulted in an increase (decrease) of the splitting between hybridized plasmonic modes for the transverse (longitudinal) polarization. Here, we see that both the bright and dark modes are redshifted, coherently with the radiative redshifts observed for the chain and also for two-dimensional lattices [39]. The bright modes are significantly more altered

than the dark ones, which is consistent with the fact that they interact more strongly with light. Therefore, the splittings we predict are reversed from those presented in Ref. [46], since here we show an increase (decrease) of the splitting between hybridized plasmonic modes for the longitudinal (transverse) polarization. Importantly, the frequency shifts we predict are also more than one order of magnitude larger than the ones presented in Ref. [46]. Hence, the experimental protocol proposed in Ref. [46] to detect the frequency splitting between bright and dark modes could be more accessible than originally thought.

Another quantity that has been proposed to be observed experimentally is the ratio $\Delta\tilde{\omega}^z/\Delta\tilde{\omega}^{x,y}$ of the longitudinal and transverse frequency splittings $\Delta\tilde{\omega}^\sigma = \tilde{\omega}_+^\sigma - \tilde{\omega}_-^\sigma$ between bright and dark modes. In the absence of light-matter interaction, the bare ratio $\Delta\omega^z/\Delta\omega^{x,y} = 2$ is independent of the center-to-center distance $d$ up to quadratic corrections in $\Omega/\omega_0 \ll 1$. In Ref. [46], it was found that the ratio $\Delta\tilde{\omega}^z/\Delta\tilde{\omega}^{x,y}$ presents a universal scaling with the distance $d$, being independent of the nanoparticle radius $a$, unlike the radiative frequency shifts for individual dimer levels.

In Fig. 9, we show the result obtained for the dimensionless ratio $\Delta\tilde{\omega}^z/\Delta\tilde{\omega}^{x,y}$ using our approach that does not use Bethe's mass renormalization. We immediately see that the universal scaling put forward in Ref. [46] is no longer present, implying that the dimensionless ratio $\Delta\tilde{\omega}^z/\Delta\tilde{\omega}^{x,y}$ depends on the nanoparticle radius. Actually, this makes the experimental realization of the protocol proposed in Ref. [46] to detect the level splitting and this dimensionless ratio even more interesting, since it could allow one to gauge the relevance of Bethe's renormalization procedure for dipolar systems such as the one studied in this appendix.

To conclude, this appendix demonstrated that not using Bethe's mass renormalization procedure leads to drastically different conclusions about the frequency shifts induced by the light-matter interaction in metallic nanoparticle dimers.
Here, we predict radiative frequency shifts of more than one order of magnitude larger than the one predicted in Ref. [46], an increase of the splitting between longitudinal bright and dark modes, and a dependence of the dimensionless ratio $\Delta\tilde{\omega}^z/\Delta\tilde{\omega}^{x,y}$ on nanoparticle radius. Importantly, within this approach, both the single nanoparticle and dimer radiative frequency shifts depend linearly on the choice of cutoff $\omega_c$, as opposed to the chain radiative frequency shift which depends

on it only logarithmically, and to two-dimensional lattices radiative shifts which are cutoff independent [39]. In the three-dimensional case, an exact diagonalization which does not require any cutoff is even possible [47]. This suggests that our quantum formalism of light-matter-induced radiative frequency shift may be more adapted to the study of periodic systems such as lattices rather than to single nanoparticles or dimers.

[1] W. L. Barnes, A. Dereux, and T. W. Ebbesen, Surface plasmon subwavelength optics, Nature **424**, 824 (2003).

[2] M. I. Stockman, Nanoplasmonics: past, present, and glimpse into future, Opt. Express **19**, 22029 (2011).

[3] M. Quinten, A. Leitner, J. R. Krenn, and F. R. Aussenegg, Electromagnetic energy transport via linear chains of silver nanoparticles, Opt. Lett. **23**, 1331 (1998).

[4] S. A. Meier, Plasmonics: Fundamentals and Applications (Springer, New York, 2007).

[5] E. Ozbay, Plasmonics: Merging Photonics and Electronics at Nanoscale Dimensions, Science **311**, 189 (2006).

[6] D. K. Gramotnev and S. I. Bozhevolnyi, Plasmonics beyond the diffraction limit, Nat. Photonics **4**, 83 (2010).

[7] J. R. Krenn, A. Dereux, J. C. Weeber, E. Bourillot, Y. Lacroute, J. P. Goudonnet, G. Schider, W. Gotschy, A. Leitner, F. R. Aussenegg, and C. Girard, Squeezing the optical near-field zone by plasmon coupling of metallic nanoparticles, Phys. Rev. Lett. **82**, 2590 (1999).

[8] S. A. Maier, M. L. Brongersma, P. G. Kik, and H. A. Atwater, Observation of near-field coupling in metal nanoparticle chains using far-field polarization spectroscopy, Phys. Rev. B **65**, 193408 (2002).

[9] S. A. Maier, P. G. Kik, H. A. Atwater, S. Meltzer, E. Harel, B. E. Koel, and A. A. G. Requicha, Local detection of electromagnetic energy transport below the diffraction limit in metal nanoparticle plasmon waveguides, Nat. Mater. **2**, 229 (2003).

[10] A. F. Koenderink, R. de Maele, J. C. Prangsma, and A. Polman, Experimental evidence for large dynamic effects on the plasmon dispersion of subwavelength metal nanoparticle waveguides, Phys. Rev. B **76**, 201403(R) (2007).

[11] K. B. Crozier, E. Togan, E. Simsek, and T. Yang, Experimental measurement of the dispersion relations of the surface plasmon modes of metal nanoparticle chains, Opt. Express **15**, 17482 (2007).

[12] A. Apuzzo, M. Février, R. Salas-Montiel, A. Bruyant, A. Chelnokov, G. Lérondel, B. Dagens, and S. Blaize, Observation of near-field dipolar interactions involved in a metal nanoparticle chain waveguide, Nano Lett. **13**, 1000 (2013).

[13] S. J. Barrow, D. Rossouw, A. M. Funston, G. A. Botton, and P. Mulvaney, Mapping bright and dark modes in gold nanoparticle chains using electron energy loss spectroscopy, Nano Lett. **14**, 3799 (2014).

[14] M. L. Brongersma, J. W. Hartman, and H. A. Atwater, Electromagnetic energy transfer and switching in nanoparticle chain arrays below the diffraction limit, Phys. Rev. B **62**, R16356 (2000).

[15] S. Y. Park and D. Stroud, Surface-plasmon dispersion relations in chains of metallic nanoparticles: An exact quasistatic calculation, Phys. Rev. B **69**, 125418 (2004).

[16] W. H. Weber and G. W. Ford, Propagation of optical excitations by dipolar interactions in metal nanoparticle chains, Phys. Rev.

B **70**, 125429 (2004).

[17] D. S. Citrin, Coherent excitation transport in metal-nanoparticle chains, Nano Lett. **4**, 1561 (2004).

[18] C. R. Simovski, A. J. Viitanen, and S. A. Tretyakov, Resonator mode in chains of silver spheres and its possible application, Phys. Rev. E **72**, 066606 (2005).

[19] D. S. Citrin, Plasmon-polariton transport in metal-nanoparticle chains embedded in a gain medium, Opt. Lett. **31**, 98 (2006).

[20] A. F. Koenderink and A. Polman, Complex response and polariton-like dispersion splitting in periodic metal nanoparticle chains, Phys. Rev. B **74**, 033402 (2006).

[21] V. A. Markel and A. K. Sarychev, Propagation of surface plasmons in ordered and disordered chains of metal nanospheres, Phys. Rev. B **75**, 085426 (2007).

[22] K. H. Fung and C. T. Chan, Plasmonic modes in periodic metal nanoparticle chains: a direct dynamic eigenmode analysis, Opt. Lett. **32**, 973 (2007).

[23] M. Conforti and M. Guasoni, Dispersive properties of linear chains of lossy metal nanoparticles, J. Opt. Soc. Am. B **27**, 1576 (2010).

[24] I. B. Udagedara, I. D. Rukhlenko, and M. Premaratne, Complex-$\omega$ approach versus complex-$k$ approach in description of gain-assisted surface plasmon-polariton propagation along linear chains of metallic nanospheres, Phys. Rev. B **83**, 115451 (2011).

[25] B. Rolly, N. Bonod, and B. Stout, Dispersion relations in metal nanoparticle chains: necessity of the multipole approach, J. Opt. Soc. Am. B **29**, 1012 (2012).

[26] C. Lee, M. Tame, J. Lim, and J. Lee, Quantum plasmonics with a metal nanoparticle array, Phys. Rev. A **85**, 063823 (2012).

[27] M. Petrov, Disorder-induced Purcell enhancement in nanoparticle chains, Phys. Rev. A **91**, 023821 (2015).

[28] A. Brandstetter-Kunc, G. Weick, C. A. Downing, D. Weinmann, and R. A. Jalabert, Nonradiative limitations to plasmon propagation in chains of metallic nanoparticles, Phys. Rev. B **94**, 205432 (2016).

[29] C. A. Downing, E. Mariani, and G. Weick, Retardation effects on the dispersion and propagation of plasmons in metallic nanoparticle chains, J. Phys. Condens. Matter **30**, 025301 (2018).

[30] P. J. Compaijen, V. A. Malyshev, and J. Knoester, Time-dependent transport of a localized surface plasmon through a linear array of metal nanoparticles: Precursor and normal mode contributions, Phys. Rev. B **97**, 085428 (2018).

[31] A. M. Pikalov, A. V. Dorofeenko, and Y. E. Lozovik, Dispersion relations for plasmons in complex-shaped nanoparticle chains, Phys. Rev. B **98**, 085134 (2018).

[32] A. O. Pinchuk and G. C. Schatz, Nanoparticle optical properties: Far- and near-field electrodynamic coupling in a chain of silver spherical nanoparticles, Mater. Sci. Eng. B **149**, 251 (2008).
[33] S. R. Pocock, P. A. Huidobro, and V. Giannini, Bulk-edge correspondence and long-range hopping in the topological plasmonic chain, *Nanophotonics* **8**, 1337 (2019).

[34] C. A. Downing and G. Weick, Topological collective plasmons in bipartite chains of metallic nanoparticles, *Phys. Rev. B* **95**, 125426 (2017).

[35] F. Schulz, O. Pavelka, F. Lehmkühler, F. Westermeier, Y. Okamura, N. S. Mueller, S. Reich, and H. Lange, Structural order in plasmonic superlattices, *Nat. Comm.* **11**, 3821 (2020).

[36] N. S. Mueller, Y. Okamura, B. G. M. Vieira, S. Juergensen, H. Lange, E. B. Barros, F. Schulz, and S. Reich, Deep strong light-matter coupling in plasmonic nanoparticle crystals, *Nature* **583**, 780 (2020).

[37] M. Tame, K. McEnery, S. Özdemir, J. Lee, S. A. Maier, and M. S. Kim, Quantum plasmonics, *Nat. Phys.* **9**, 329 (2013).

[38] K. Miwa and G. C. Schatz, Quantum electrodynamics description of localized surface plasmons at a metal nanosphere, *Phys. Rev. A* **103**, L041501 (2021).

[39] F. Fernique and G. Weick, Plasmons in two-dimensional lattices of near-field coupled nanoparticles, *Phys. Rev. B* **102**, 045420 (2020).

[40] G. Weick, G.-L. Ingold, R. A. Jalabert, and D. Weinmann, Surface plasmon in metallic nanoparticles: Renormalization effects due to electron-hole excitations, *Phys. Rev. B* **74**, 165421 (2006).

[41] C. A. Downing and G. Weick, Plasmonic modes in cylindrical nanoparticles and dimers, *Proc. R. Soc. A.* **476**, 20200530 (2020).

[42] C. R. Mann, T. J. Sturges, G. Weick, W. L. Barnes, and E. Mariani, Manipulating type-I and type-II Dirac polaritons in cavity-embedded honeycomb metasurfaces, *Nat. Commun.* **9**, 2194 (2018).

[43] C. A. Downing, T. J. Sturges, G. Weick, M. Stobińska, and L. Martín-Moreno, Topological phases of polaritons in a cavity waveguide, *Phys. Rev. Lett.* **123**, 217401 (2019).

[44] C. A. Downing and L. Martín-Moreno, Polaritonic Tamm states induced by cavity photons, *Nanophotonics* **10**, 513 (2021).

[45] A. Brandstetter-Kunc, G. Weick, D. Weinmann, and R. A. Jalabert, Decay of dark and bright plasmonic modes in a metallic nanoparticle dimer, *Phys. Rev. B* **91**, 035431 (2015), **92**, 199906(E) (2015).

[46] C. A. Downing, E. Mariani, and G. Weick, Radiative frequency shifts in nanoplasmonic dimers, *Phys. Rev. B* **96**, 155421 (2017).

[47] S. Lamowski, C. R. Mann, F. Hellbach, E. Mariani, G. Weick, and F. Pauly, Plasmon polaritons in cubic lattices of spherical metallic nanoparticles, *Phys. Rev. B* **97**, 125409 (2018).

[48] J. J. Hopfield, Theory of the contribution of excitons to the complex dielectric constant of crystals, *Phys. Rev.* **112**, 1555 (1958).

[49] J. A. Scholl, A. García-Etxarri, A. L. Koh, and J. A. Dionne, Observation of quantum tunneling between two plasmonic nanoparticles, *Nano Lett.* **13**, 564 (2013).

[50] D. P. Craig and T. Thirunamachandran, *Molecular Quantum Electrodynamics* (Academic Press, London, 1984).

[51] A. Kawabata and R. Kubo, Electronic Properties of Fine Metallic Particles. II. Plasma Resonance Absorption, *J. Phys. Soc. Jpn.* **21**, 1765 (1966).

[52] G. F. Bertsch and R. A. Broglia, *Oscillations in Finite Quantum Systems* (Cambridge University Press, Cambridge, 1994).

[53] U. Kreibig and M. Vollmer, *Optical Properties of Metal Clusters* (Springer-Verlag, Berlin, 1995).

[54] F. Fernique, *Collective plasmonic excitations in two-dimensional metamaterials based on near-field coupled metallic nanoparticles*, PhD thesis, Université de Strasbourg (2019).

[55] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, *Atom-Photon Interactions: Basic Processes and Applications* (Wiley, New York, 1992).

[56] J. Southall, D. Hodgson, R. Purdy, and A. Beige, Locally acting mirror Hamiltonians, *J. Mod. Opt.* **68**, 647 (2021).

[57] D. Hodgson, J. Southall, R. Purdy, and A. Beige, Quantising the electromagnetic field in position space, arXiv:2104.04499.

[58] G. Mie, Beiträge zur optik trüber medien, speziell kolloidaler metallösungen, *Ann. Phys.* (Leipzig) **330**, 377 (1908).

[59] C. F. Bohren and D. R. Huffman, *Absorption and Scattering of Light by Small Particles* (Wiley-VCH, Weinheim, 2004).

[60] F. Capolino, *Theory and Phenomena of Metamaterials* (Taylor & Francis, Boca Raton, 2009).

[61] W. A. Jacak, Exact solution for velocity of plasmon-polariton in metallic nano-chain, *Opt. Express* **22**, 18958 (2014).

[62] B. Auguié and W. L. Barnes, Collective resonances in gold nanoparticle arrays, *Phys. Rev. Lett.* **101**, 143902 (2008).

[63] V. G. Kravets, A. V. Kabashin, W. L. Barnes, and A. N. Grigorenko, Plasmonic Surface Lattice Resonances: A Review of Properties and Applications, *Chem. Rev.* **118**, 5912 (2018).

[64] H. A. Bethe, The electromagnetic shift of energy levels, *Phys. Rev.* **72**, 339 (1947).

[65] J. D. Jackson, *Classical Electrodynamics*, 3rd ed. (Wiley, New York, 2007).

[66] S. R. Pocock, X. Xiao, P. A. Huidobro, and V. Giannini, Topological plasmonic chain with retardation and radiative effects, *ACS Photonics* **5**, 2271 (2018).

[67] W. T. Doyle, Optical properties of a suspension of metal spheres, *Phys. Rev. B* **39**, 9852 (1989).

[68] A. Wokaun, J. P. Gordon, and P. F. Liao, Radiation damping in surface-enhanced Raman scattering, *Phys. Rev. Lett.* **48**, 957 (1982).

[69] B. T. Draine, The discrete-dipole approximation and its application to interstellar graphite grains, *Astrophys. J.* **333**, 848 (1988).

[70] A. Moroz, Depolarization field of spheroidal particles, *J. Opt. Soc. Am. B* **26**, 517 (2009).

[71] M. Meier and A. Wokaun, Enhanced fields on large metal particles: dynamic depolarization, *Opt. Lett.* **8**, 581 (1983).

[72] V. Grigoriev, N. Bonod, and B. Stout, Optimizing nanoparticle designs for ideal absorption of light, *ACS Photonics* **2**, 263 (2015).

[73] M. D. Schwartz, *Quantum Field Theory and the Standard Model* (Cambridge University Press, Cambridge, 2013).

[74] V. Rokaj, D. M. Welakuh, M. Ruggenthaler, and A. Rubio, Light-matter interaction in the long-wavelength limit: no ground-state without dipole self-energy, *J. Phys. B* **51**, 034005 (2018).

[75] C. Schäfer, M. Ruggenthaler, V. Rokaj, and A. Rubio, Relevance of the quadratic diamagnetic and self-polarization terms in cavity quantum electrodynamics, *ACS Photonics* **7**, 975 (2020).

[76] O. D. Stefano, A. Settineri, V. Macri, L. Garziano, R. Stassi, S. Savasta, and F. Nori, Resolution of gauge ambiguities in ultrastrong-coupling cavity quantum electrodynamics, *Nat. Phys.* **15**, 803 (2019).

[77] A. F. Kockum, A. Miranowicz, S. D. Liberato, S. Savasta, and F. Nori, Ultrastrong coupling between light and matter, *Nat. Rev. Phys.* **1**, 19 (2019).

[78] P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, Ultrastrong coupling regimes of light-matter interaction, *Rev. Mod. Phys.* **91**, 025005 (2019).