Topological plasmons in dimerized chains of nanoparticles: robustness against long-range quasistatic interactions and retardation effects

TL;DR

This paper models topological plasmons in dimerized nanoparticle chains, demonstrating their robustness against long-range interactions and retardation effects, and analyzing their spectral shifts and topological properties.

Contribution

It provides an exact analytical bandstructure including long-range interactions and shows topological robustness despite broken chiral symmetry.

Findings

Topological edge states are robust against long-range Coulomb interactions.

Radiative shifts depend on nanoparticle size and polarization.

Topologically protected states remain unaffected by plasmon-photon coupling.

Abstract

We present a simple model of collective plasmons in a dimerized chain of spherical metallic nanoparticles, an elementary example of a topologically nontrivial nanoplasmonic system. Taking into account long-range quasistatic dipolar interactions throughout the chain, we provide an exact analytical expression for the full quasistatic bandstructure of the collective plasmons. An explicit calculation of the Zak phase proves the robustness of the topological physics of the system against the inclusion of long-range Coulomb interactions, despite the broken chiral symmetry. Using an open quantum systems approach, which includes retardation through the plasmon-photon coupling, we go on to analytically evaluate the resulting radiative frequency shifts of the plasmonic spectrum. The bright plasmonic bands experience size-dependent radiative shifts, while the dark bands are essentially unaffected by the light-matter coupling. Notably, the upper transverse-polarized band presents a logarithmic singularity where the quasistatic spectrum intersects the light cone. At wavevectors away from this intersection and for subwavelength nanoparticles, the plasmon-photon coupling only leads to a quantitative reconstruction of the bandstructure and the topologically-protected states at the edge of the first Brillouin zone are essentially unaffected.

Figures

Equations

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

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

$$H_{\mathrm{pl}}^0=\hslash {\omega }_0\sum _{\sigma =x,y,z}\sum _{n=1}^{\mathcal{N}}\left({a_n^{\sigma }}^{†}{a}_n^{\sigma }+{b_n^{\sigma }}^{†}{b}_n^{\sigma }\right)$$

$$H_{\mathrm{pl}}^{\mathrm{int}}=H_{\mathrm{pl}}^{A\leftrightarrow A,B\leftrightarrow B}+H_{\mathrm{pl}}^{A\leftrightarrow B}.$$

$$H_{\mathrm{pl}}^{A\leftrightarrow A,B\leftrightarrow B}=$$

$$\times [\left({a_n^{\sigma }}^{†}+a_n^{\sigma }\right)\left({a_m^{\sigma }}^{†}+a_m^{\sigma }\right)$$

$$+\left({b_n^{\sigma }}^{†}+b_n^{\sigma }\right)\left({b_m^{\sigma }}^{†}+b_m^{\sigma }\right)].$$

$$H_{\mathrm{pl}}^{A\leftrightarrow B}=$$

$$\times \left({a_n^{\sigma }}^{†}+a_n^{\sigma }\right)\left({b_m^{\sigma }}^{†}+b_m^{\sigma }\right)$$

$$+\sum _{m=n+1}^{\mathcal{N}}\frac{{\eta }_{\sigma }}{(m-n-1+d_2/d{)}^3}$$

$$\displaystyle\times\left({b_{n}^{\sigma}}^{\dagger}+b_{n}^{\sigma}\right)\left({a_{m}^{\sigma}}^{\dagger}+a_{m}^{\sigma}\right)\Bigg{\}}.$$

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

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

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

$$\displaystyle+\frac{N_{\mathrm{e}}e^{2}}{2m_{\mathrm{e}}c^{2}}\mathbf{A}^{2}(\mathbf{d}_{n,s})\bigg{\}},$$

$${\mathbf{\Pi }}_{n,A}=\mathrm{i}\sqrt{\frac{N_{\mathrm{e}}{m}_{\mathrm{e}}\hslash {\omega }_0}{2}}\sum _{\sigma =x,y,z}\hat{\sigma }({a_n^{\sigma }}^{†}-a_n^{\sigma }),$$

$${\mathbf{\Pi }}_{n,B}=\mathrm{i}\sqrt{\frac{N_{\mathrm{e}}{m}_{\mathrm{e}}\hslash {\omega }_0}{2}}\sum _{\sigma =x,y,z}\hat{\sigma }({b_n^{\sigma }}^{†}-b_n^{\sigma }).$$

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

$$a_n^{\sigma }$$

$$b_n^{\sigma }$$

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

$$+\frac{\hslash \Omega }{2}{f}_q^{\sigma }\left({a_q^{\sigma }}^{†}{a_{-q}^{\sigma }}^{†}+a_{-q}^{\sigma }{a}_q^{\sigma }+{b_q^{\sigma }}^{†}{b_{-q}^{\sigma }}^{†}+b_{-q}^{\sigma }{b}_q^{\sigma }\right)$$

$$\displaystyle+\hbar\Omega\left[g_{q}^{\sigma}{a_{q}^{\sigma}}^{\dagger}\left(b_{q}^{\sigma}+{b_{-q}^{\sigma}}^{\dagger}\right)+{g_{q}^{\sigma}}^{\ast}{a_{q}^{\sigma}}\left({b_{q}^{\sigma}}^{\dagger}+{b_{-q}^{\sigma}}\right)\right]\bigg{\}},$$

$$f_q^{\sigma }=2{\eta }_{\sigma }\sum _{n=1}^{\infty }\frac{\cos (nqd)}{n^3}$$

$$g_q^{\sigma }={\eta }_{\sigma }\sum _{n=0}^{\infty }\left[\frac{{\mathrm{e}}^{\mathrm{i}nqd}}{(n+d_1/d{)}^3}+\frac{{\mathrm{e}}^{-\mathrm{i}(n+1)qd}}{(n+d_2/d{)}^3}\right],$$

$$f_q^{\sigma }={\eta }_{\sigma }\left[{\mathrm{Li}}_3\left({\mathrm{e}}^{\mathrm{i}qd}\right)+{\mathrm{Li}}_3\left({\mathrm{e}}^{-\mathrm{i}qd}\right)\right]$$

$$g_q^{\sigma }={\eta }_{\sigma }\left[\Phi \left({\mathrm{e}}^{\mathrm{i}qd},3,d_1/d\right)+{\mathrm{e}}^{-\mathrm{i}qd}\Phi \left({\mathrm{e}}^{-\mathrm{i}qd},3,d_2/d\right)\right].$$

$$H_{q\tau }^{\sigma }=\hslash {\omega }_{q\tau }^{\sigma }{{\beta }_{q\tau }^{\sigma }}^{†}{\beta }_{q\tau }^{\sigma },$$

$${\omega }_{q\tau }^{\sigma }={\omega }_0\sqrt{1+2\frac{\Omega }{{\omega }_0}\left(f_q^{\sigma }+\tau |g_q^{\sigma }|\right)},$$

$${\beta }_{q\tau }^{\sigma }=\cosh {\theta }_{q\tau }^{\sigma }{\alpha }_{q\tau }^{\sigma }+\sinh {\theta }_{q\tau }^{\sigma }{{\alpha }_{-q\tau }^{\sigma }}^{†},$$

$${\alpha }_{q\tau }^{\sigma }=\frac{1}{\sqrt{2}}\left(a_q^{\sigma }+\tau \frac{g_q^{\sigma }}{|g_q^{\sigma }|}{b}_q^{\sigma }\right),$$

$$\cosh {\theta }_{q\tau }^{\sigma }=\frac{{\omega }_{q\tau }^{\sigma }+{\omega }_0}{2\sqrt{{\omega }_0{\omega }_{q\tau }^{\sigma }}},$$

$$\sinh {\theta }_{q\tau }^{\sigma }=\frac{{\omega }_{q\tau }^{\sigma }-{\omega }_0}{2\sqrt{{\omega }_0{\omega }_{q\tau }^{\sigma }}}.$$

$$ϵ=\frac{d_1-d_2}{d}.$$

$${\omega }_{\mathrm{nn},q\tau }^{\sigma }={\omega }_0\sqrt{1+2\tau \frac{\Omega }{{\omega }_0}|g_{\mathrm{nn},q}^{\sigma }|}.$$

$$g_{\mathrm{nn},q}^{\sigma }={\eta }_{\sigma }{d}^3\left(\frac{1}{d_1^3}+\frac{{\mathrm{e}}^{-\mathrm{i}qd}}{d_2^3}\right).$$

Full text

1 Introduction

The simple system of a chain of regularly-spaced sites coupled by short-range interactions is an insightful model across many areas of physics. Dimerized chains are of only a minor additional complexity, but present significant differences in the physics manifested in the system. For example, Peierls made the profound discovery in the 1930’s that the distortion of a one-dimensional atomic chain may be energetically favorable, leading to the celebrated Peierls transition Peierls1955 . While modeling electrons hopping along the alternating single and double C–C bonds of polyacetylene, Su, Schrieffer and Heeger wrote down a physically rich Hamiltonian which encompasses solitons and charge fractionalization Su1979 ; Su1980 ; Heeger1988 . More recently, Kitaev proposed an influential toy model which supports unpaired Majorana modes at the end of a quantum wire Kitaev2001 .

The field of plasmonics Maier2007 ; Stockman2011 aims to control and manipulate plasmons, the collective oscillations of the electron gas, in order to facilitate technological applications at the nanoscale, including, e.g. sensing and waveguiding Halas2007 . The localized surface plasmon (LSP) is a particularly promising quasiparticle to exploit due to advances in the controlled synthesis of metallic nanoparticles, which host such a collective excitation Bertsch ; Kreibig .

Various arrays of nanoparticles have been investigated Meinzer2014 ; Wang2018 , where Coulomb interactions between the nanoparticles give rise to collective plasmons spread out over the whole metamaterial. The simplest array hosting collective plasmons is the linear chain of metal nanoparticles, which is the prototypical plasmonic waveguide Quinten1998 . The regular plasmonic chain has been extensively studied both theoretically Brongersma2000 ; Maier2003b ; Park2004 ; Weber2004 ; Citrin2004 ; Simovski2005 ; Citrin2006 ; Koenderink2006 ; Markel2007 ; Fung2007 ; Lee2012 ; Compaijen2013 ; Pino2014 ; Compaijen2015 ; Petrov2015 ; Brandstetter2016 ; Downing2018 ; Compaijen2018 and experimentally Krenn1999 ; Maier2002 ; Maier2003a ; Koendrick2007 ; Crozier2007 ; Apuzzo2013 ; Barrow2014 ; Gur2017 . Crucially, it was noticed theoretically that retardation in the dipole-dipole interaction amongst the LSPs can have a dramatic influence on the plasmonic bandstructure Weber2004 ; Citrin2004 ; Simovski2005 ; Citrin2006 ; Koenderink2006 ; Markel2007 ; Fung2007 ; Compaijen2015 ; Petrov2015 ; Downing2018 ; Compaijen2018 . In particular, it was demonstrated that plasmonic modes polarized perpendicular to the chain, and with wavenumbers corresponding to the intersection of the quasistatic plasmonic dispersion and the light cone, experience a significant radiative redshift.

The inevitable marriage of dimerized chains and nanoplasmonics was studied theoretically in recent years, primarily at the level of nearest-neighbor interactions and neglecting retardation effects Ling2015 ; DowningBiPartite2017 ; Gomez2017 . Within the two above-mentioned approximations, it was shown that a bipartite chain of plasmonic nanoparticles should host midgap topological edge states which are robust against disorder. Zigzag chains of plasmonic particles offer an alternative to the dimerized geometry and were studied both theoretically Poddubny2014 and experimentally Sinev2015 .

Nonplasmonic one-dimensional photonic systems such as photonic lattices Kanshu2012 ; Schomerus2013 , dielectric resonators Poli2015 ; Slobo2015 , silicon waveguides Blanco2016 , polariton cavities Solny2016 ; Solny2016b , metallic nanowires liu16_preprint , and microring resonators Zhao2018 ; Parto2018 were also demonstrated as being the hosts of nontrivial topological states of light. Higher-dimensional topological artificial systems are also at present intensively studied, leading to many possible exciting applications in photonics (for recent reviews, see references Lu2014 ; Ozawa2018 ).

Recently, Pocock et al. Pocock2018 legitimately challenged the findings of references Ling2015 ; DowningBiPartite2017 ; Gomez2017 against the long-ranged nature of the quasistatic dipole-dipole interaction, as well as retardation effects and radiation damping stemming from the finite velocity of light. Based on numerical solutions to Maxwell’s equations, it was argued that a nontrivial Zak phase Berry1984 ; Zak1989 ; Asboth2016 , which dictates the existence of topologically-protected edge states via the bulk-edge correspondence Hasan2010 ; Delplace2011 ; Rhim2017 , persists even in the presence of such long-ranged, retarded interactions.

In this work, we extend the theory of the dimerized plasmonic chain with a transparent model, built within an open quantum system approach Brandstetter2015 ; Brandstetter2016 . In our model, we account for long-range quasistatic interactions exactly, which enables us to obtain an analytical expression for the plasmonic bandstructure. By analyzing the properties of the Bloch Hamiltonian, we show that despite the broken chiral (or sublattice) symmetry, time-reversal and inversion symmetries are preserved, guaranteeing a quantized Zak phase vanMiert2017 which we calculate exactly.

We explicitly include retardation through the light-matter minimal coupling Hamiltonian, which we account for perturbatively. Our analytical treatment is advantageous as compared to a numerical classical approach Pocock2018 , which requires the introduction of a radiative correction to the polarizability via a term which can lead to acausal behavior Weber2004 ; Jackson1998 . Our modus operandi here is underpinned by the fluctuation-dissipation theorem Kubo1966 : the photonic environment leads to both a finite radiative lifetime and a shift in the plasmonic energy levels Downing2017 , in direct analogy with the famous Lamb shift of atomic physics Lamb1947 ; Bethe1947 ; Milonni1994 . The radiative frequency shifts calculated here are a manifestation of the cooperative Lamb shift, well known from many-atom systems Friedberg1973 ; Scully2009 , and give rise to insignificant modifications to the plasmonic bandstructure for small-enough nanoparticles, thus preserving the existence of the topologically-protected edge states.

The paper is organized as follows: We present our model of a dimerized chain of plasmonic nanoparticles coupled to vacuum photonic modes in section 2, and we derive a simple expression for the exact quasistatic plasmonic bandstructure in section 3. In section 4, we comment on the robustness of the topologically nontrivial physics of the system. We unveil analytical expressions for the radiative shifts of the collective plasmonic bandstructure in section 5. Finally, we draw some conclusions in section 6. Details of our calculations, as well as complementary analytical results on the radiative damping decay rates of the system, are presented in the appendices.

2 Model

We consider a dimerized one-dimensional array of spherical metallic nanoparticles, each of radius $a$, characterized by two inequivalent sublattices which we denote $A$ and $B$ (see figure 1). The formed bipartite chain is characterized by two interparticle separations, $d_{1}$ and $d_{2}$, which alternate along the array and hence give rise to the periodicity $d=d_{1}+d_{2}$. Every nanoparticle harbors three degenerate LSPs, each polarized along one of the two transverse directions $x$ and $y$, or along the longitudinal direction $z$. The LSP corresponds to an oscillation of the electronic center of mass at the Mie frequency $\omega_{0}$ Bertsch ; Kreibig . For nanoparticles in vacuum, and neglecting the screening of the valence electrons by the core electrons, the Mie frequency corresponds to $\omega_{\mathrm{p}}/\sqrt{3}$, where $\omega_{\mathrm{p}}=\sqrt{4\pi n_{\mathrm{e}}e^{2}/m_{\mathrm{e}}}$ is the plasma frequency. Here, $n_{\mathrm{e}}$ is the electronic density of the considered metal, $-e<0$ the electronic charge, and $m_{\mathrm{e}}$ its mass.111Throughout the paper, we use cgs units.

Crucially, the LSPs in the chain couple to form collective plasmons spread out over the whole chain due to Coulomb interactions, which are essentially dipolar for the interparticle separations $d_{1,2}\gtrsim 3a$ we consider here Park2004 . Retardation effects, accounted for perturbatively in our model by the coupling of the plasmonic modes to the photonic environment, lead to a reconstruction of the resulting collective plasmonic bandstructure.

Within the Coulomb gauge cohen ; craig , the fully-retarded Hamiltonian of the plasmonic chain coupled to vacuum electromagnetic modes is

[TABLE]

where $H_{\mathrm{pl}}$ and $H_{\mathrm{ph}}$ respectively describe the electronic and photonic subsystems coupled through the light-matter interaction $H_{\mathrm{pl}\textrm{-}\mathrm{ph}}$.222In contrast to reference DowningBiPartite2017 , in this work we focus on radiative effects and neglect Landau damping Bertsch ; Kreibig ; Kawabata1966 ; Weick2005 , and hence the associated shift it induces in the plasmonic resonance frequency Weick2006 . This is justified as long as the nanoparticles are not too tiny (i.e. those with a radius $a\lesssim 5\,\mathrm{nm}$), rendering quantum-size effects irrelevant.

The plasmonic Hamiltonian, describing the LSPs coupled via the long-ranged quasistatic dipolar interaction, reads

[TABLE]

where the noninteracting part

[TABLE]

corresponds to the Hamiltonian for uncoupled harmonic LSPs. Here, the index $n\in[1,\mathcal{N}]$ denotes the dimer number in the chain of $2\mathcal{N}$ nanoparticles333Without loss of generality, we consider chains with an even number of nanoparticles. and $\sigma$ labels the two transverse ($x,y$) and the single longitudinal ($z$) polarizations of the plasmonic modes (see figure 1). In equation (3), the bosonic operators $a_{n}^{\sigma}$ (${a_{n}^{\sigma}}^{\dagger}$) and $b_{n}^{\sigma}$ (${b_{n}^{\sigma}}^{\dagger}$) annihilate (create) an LSP with polarization $\sigma$ in dimer $n$ in a nanoparticle belonging to the $A$ and $B$ sublattices, respectively.

In equation (2), the long-range quasistatic dipolar interaction among the LSPs is encapsulated in the term $H_{\mathrm{pl}}^{\mathrm{int}}$, which can be decomposed into two contributions:

[TABLE]

The first of these contributions corresponds to the coupling between pairs of nanoparticles belonging to the same sublattices ($A\leftrightarrow A$ and $B\leftrightarrow B$) and reads

[TABLE]

The second contribution $H_{\mathrm{pl}}^{A\leftrightarrow B}$ to equation (4) is due to the coupling between pairs of particles belonging to the two inequivalent sublattices ($A\leftrightarrow B$). It can be written as

[TABLE]

In equations (2) and (2), the coupling constant $\Omega$, defined by

[TABLE]

is much smaller than $\omega_{0}$ for all realistic values of the ratio $d/a$,444Note that touching nanoparticles would correspond to $d/a=4$, i.e. $\Omega/\omega_{0}=1/128$. while the polarization-dependent factor $\eta_{x,y}=1$ ($\eta_{z}=-2$) for the transverse (longitudinal) modes arises from the anisotropy of the dipole-dipole interaction.

We note that the nearest-neighbor approximation adopted in reference DowningBiPartite2017 amounts to two simplifications in the interacting Hamiltonian (4). Firstly, one disregards the intrasublattice Hamiltonian (2). Secondly, one retains in the summation over the dimer index $m$ only the contributions $m=n$ and $m=n+1$ in the first and second terms on the right-hand side of the intersublattice Hamiltonian (2), respectively (cf. equation (2) in reference DowningBiPartite2017 ).

In equation (1), the photonic environment is described by the Hamiltonian

[TABLE]

where $c_{\mathbf{k}}^{\hat{\lambda}_{\mathbf{k}}}$ (${c_{\mathbf{k}}^{\hat{\lambda}_{\mathbf{k}}}}^{\dagger}$) annihilates (creates) a photon with wavevector $\mathbf{k}$ and transverse polarization $\hat{\lambda}_{\mathbf{k}}$, so that $\mathbf{k}\cdot\hat{\lambda}_{\mathbf{k}}=0$ (here and in what follows, hats designate unit vectors). The photonic dispersion is $\omega_{\mathbf{k}}=c|\mathbf{k}|$, where $c$ is the speed of light in vacuum.

In the long-wavelength approximation $|\mathbf{k}|a\ll 1$, the plasmon-photon coupling Hamiltonian in equation (1) is

[TABLE]

where $\mathbf{d}_{n,A}=d(n-1)\hat{z}$ and $\mathbf{d}_{n,B}=\mathbf{d}_{n,A}+d_{1}\hat{z}$ correspond, respectively, to the location of the center of the nanoparticle belonging to the $A$ and $B$ sublattices in the dimer $n$ (see figure 1). The momenta associated with the LSPs on the $A$ and $B$ sublattices are

[TABLE]

The vector potential appearing in equation (2) is given by

[TABLE]

where we have quantized the electromagnetic modes in a box of volume $\mathcal{V}$. Notably, the Hamiltonian (2) fully takes into account retardation effects in the dipole-dipole interaction between the LSPs craig ; Downing2018 ; Lamowski2018 .555Note that, since we consider interparticle separation distances much smaller than the wavelength associated with the LSP resonances, we neglect Umklapp processes in equations (8) and (2).

3 Quasistatic plasmonic bandstructure

We begin our analysis of the model presented in section 2 by considering the purely plasmonic Hamiltonian (2) and deriving its associated quasistatic spectrum, thus gaining an understanding of the bandstructure without the influence of retardation effects.

Since we are interested in long chains where $\mathcal{N}\gg 1$,666Such a limit has been shown to be a good approximation for chains comprising $20$ or more nanoparticles Weber2004 ; Brandstetter2016 ; Downing2018 . it is expedient to employ periodic boundary conditions. We move into wavevector space via the pair of Fourier transforms

[TABLE]

where the plasmonic wavevector $q=2\pi p/\mathcal{N}d$, with the integer $p\in[-\mathcal{N}/2,+\mathcal{N}/2]$. Hence the plasmonic Hamiltonian (2) becomes

[TABLE]

where we have introduced the lattice sums

[TABLE]

and

[TABLE]

which correspond to the intrasublattice ($A\leftrightarrow A$ and $B\leftrightarrow B$) and intersublattice ($A\leftrightarrow B$) interactions, respectively. The latter quantities can be expressed analytically in terms of the polylogarithm function $\text{Li}_{s}(z)=\sum_{n=1}^{\infty}z^{n}n^{-s}$ and the Lerch transcendent $\Phi(z,s,u)=\sum_{n=0}^{\infty}z^{n}(n+u)^{-s}$ as follows:

[TABLE]

and

[TABLE]

The behavior of the intra- and intersublattice functions, respectively equations (16) and (17), is analyzed in detail in appendix A.

Since the plasmonic Hamiltonian (3) is quadratic, it can be diagonalized exactly with the aid of a bosonic Bogoliubov transformation. We arrive at the form $H_{\mathrm{pl}}=\sum_{\sigma q\tau}H_{q\tau}^{\sigma}$, where

[TABLE]

with the index $\tau=+$ ($-$) corresponding to the upper (lower) collective plasmonic band. The eigenfrequencies $\omega_{q\tau}^{\sigma}$ of the collective plasmonic modes are

[TABLE]

which provides a fully analytical description of the quasistatic plasmonic bandstructure.

The bosonic Bogoliubov operators appearing in the Hamiltonian (18) are given by

[TABLE]

with

[TABLE]

where the coefficients are

[TABLE]

The Bogoliubov operators $\beta_{q\tau}^{\sigma}$ and ${\beta_{q\tau}^{\sigma}}^{\dagger}$ act on an eigenstate $|n_{q\tau}^{\sigma}\rangle$ of the Hamiltonian (18), which represents a collective plasmon with polarization $\sigma$, wavevector $q$, band index $\tau$, and eigenenergy $\hbar\omega_{q\tau}^{\sigma}$, as follows: $\beta_{q\tau}^{\sigma}|n_{q\tau}^{\sigma}\rangle=\sqrt{n_{q\tau}^{\sigma}}|n_{q\tau}^{\sigma}-1\rangle$ and ${\beta_{q\tau}^{\sigma}}^{\dagger}|n_{q\tau}^{\sigma}\rangle=\sqrt{n_{q\tau}^{\sigma}+1}|n_{q\tau}^{\sigma}+1\rangle$, respectively. Here $n_{q\tau}^{\sigma}$ is a non-negative integer.

The quasistatic plasmonic bandstructure of the system, unperturbed by the photonic environment, is analytically encapsulated in equation (19). It is plotted in figure 2 (see solid lines) for both the transverse (panel (a)) and longitudinal polarizations (panel (b)) for a chain with periodicity $d=8a$ and for different degrees of dimerization, characterized by the parameter $\epsilon$ defined through the relations $d_{1}=(1+\epsilon)d/2$ and $d_{2}=(1-\epsilon)d/2$, or equivalently777For interparticle distances $d_{1,2}\geqslant 3a$, the interactions amongst the LSPs are essentially dipolar Park2004 , which restrict the dimerization parameter in our model to $|\epsilon|\leqslant 1-6a/d$. Note that $\epsilon=0$ corresponds to the regular (undimerized) chain.

[TABLE]

In the figure, we also plot as dashed lines the plasmonic dispersion considering only the dipolar interaction between nearest neighbors (nn) in the chain, explicitly (see equation (15) in reference DowningBiPartite2017 )

[TABLE]

Here the quantity which plays the role of the intersublattice function (17), restricted to nearest neighbors only, is

[TABLE]

The weak coupling constant (7) implies that the nearest-neighbor spectrum (24) is essentially $\omega_{\mathrm{nn},q\tau}^{\sigma}\simeq\omega_{0}+\tau\Omega|g_{\mathrm{nn},q}^{\sigma}|$. Hence $\omega_{\mathrm{nn},\,q\tau}^{\sigma}$ possesses a band reversal symmetry around the resonance frequency $\omega_{0}$ (or, in the language of solid state physics, a particle-hole symmetry). This symmetry is broken when long-range interactions between LSPs belonging to the same sublattice are included in the model, due to the appearance of the intrasublattice function (16) in equation (19), which also breaks the chiral symmetry of the system Asboth2016 , such that $\omega_{q\tau}^{\sigma}\simeq\omega_{0}+\Omega(f_{q}^{\sigma}+\tau|g_{q}^{\sigma}|)$. This raises the question of the impact on the protection of the topological properties of the system. We discuss this important issue in section 4.

As can be seen from the solid lines in panels (c) and (e) of figure 2 (see also panels (a) and (b)), the bandgap $\Delta\omega^{\sigma}=\omega^{\sigma}_{\pi/d,+}-\omega^{\sigma}_{\pi/d,-}\simeq\omega_{0}(a/d)^{3}|g_{\pi/d}^{\sigma}|$ at the edge of the Brillouin zone increases as a function of the dimerization parameter $|\epsilon|$. Such a bandgap is almost unaffected by the long-range dipolar interactions, as it very well matches the weak coupling, nearest-neighbor bandgap

[TABLE]

represented by a dashed line in figures 2(c) and 2(e).

It is further apparent from figure 2 (panels (d) and (f)) that the bandwidth of the high- ($\tau=+$) and low-energy ($\tau=-$) bands $\delta\omega_{\tau}^{\sigma}=|\omega_{0,\tau}^{\sigma}-\omega_{\pi/d,\tau}^{\sigma}|$, which decreases for increasing $|\epsilon|$, is noticably affected by the long-range interactions. Indeed, while in the weak-coupling, nearest-neighbor approximation, such a bandwidth is $\tau$-independent and given by

[TABLE]

represented by dashed lines in figures 2(d) and 2(f), the long-range interactions induce an asymmetry in the bandwith of the high- and low-energy branches (solid lines).

At the edge of the first Brillouin zone, we find that the collective plasmonic bandstructure (19) maps onto the form of a quasirelativistic spectrum, which indicates the presence of a nontrivial, topologically-protected edge state within the gap DowningBiPartite2017 . Performing a Taylor expansion at $qd=\pi+kd$ yields to leading order in $kd\ll 1$

[TABLE]

where $A^{\sigma}$ and $B^{\sigma}$ are (polarization-dependent) constants, given explicitly in equations (40) and (B) in appendix B, and the renormalized frequency $\omega_{0}^{\sigma}=\omega_{0}-3\zeta(3)\eta_{\sigma}\Omega/2$, where $\zeta(3)\simeq 1.20$ is Apéry’s constant. Notably, $\omega_{0}^{\sigma}$, which reduces in the nearest-neighbor approximation to $\omega_{0}$ DowningBiPartite2017 , is polarization dependent due to the long-range interactions. Up to the constant energy shift $\hbar\omega_{0}^{\sigma}$, the dispersion (28) recalls the spectrum of one-dimensional Dirac particles with mass $m$ and momentum $p$, $E_{\mathrm{D}}=\pm\sqrt{(mc^{2})^{2}+(pc)^{2}}$. Previously, this mapping was known to hold only in the nearest-neighbor approximation (see equation (16) in reference DowningBiPartite2017 ).

4 Robustness of the topological properties of the collective plasmons within the quasistatic limit

Here we investigate the robustness of the nontrivial topological nature of the dimerized plasmonic chain against long-range quasistatic interactions. To do so, we analyze the symmetries of the plasmonic Bloch Hamiltonian and calculate exactly the Berry phase Berry1984 specialized to one-dimensional systems, the so-called Zak phase Zak1989 ; Asboth2016 .

We start by writing the plasmonic Hamiltonian (3) in Bogoliubov–de Gennes form as $H_{\mathrm{pl}}=\frac{1}{2}\sum_{\sigma q}\hat{\Psi_{q}^{\sigma}}^{\dagger}\mathcal{H}_{q}^{\sigma}\hat{\Psi}_{q}^{\sigma}$, with the Bloch Hamiltonian

[TABLE]

and where $\hat{\Psi}_{q}^{\sigma}=(a_{q}^{\sigma},b_{q}^{\sigma},a_{-q}^{\sigma\dagger},b_{-q}^{\sigma\dagger})$. The $2\pi/d$-periodic part of the Bloch wavefunction, denoted $|\psi_{q\tau}^{\sigma}\rangle$, satisfies the eigenproblem $\mathcal{H}_{q}^{\sigma}\mathcal{J}|\psi_{q\tau}^{\sigma}\rangle=\hbar\omega_{q\tau}^{\sigma}|\psi_{q\tau}^{\sigma}\rangle$, where the plasmonic dispersion $\omega_{q\tau}^{\sigma}$ is given by equation (19). Here $\mathcal{J}=\sigma_{z}\otimes\mathbbm{1}_{2}$, in terms of the third Pauli matrix $\sigma_{z}$ and the $2\times 2$ identity matrix $\mathbbm{1}_{2}$. This complication is due to the requirement of obeying bosonic statistics Tsallis1978 ; Colpa1978 . Explicitly, in terms of the phase $\phi_{q}^{\sigma}$ defined via

[TABLE]

we obtain

[TABLE]

which may be found via a formal identification with the Bogoliubov coefficients (20). Notably, $|\psi_{q\tau}^{\sigma}\rangle$ is normalized with respect to the metric $\mathcal{J}$, $\langle\psi_{q\tau}^{\sigma}|\psi_{q\tau}^{\sigma}\rangle_{\mathcal{J}}=\langle\psi_{q\tau}^{\sigma}|\mathcal{J}|\psi_{q\tau}^{\sigma}\rangle=1$. We should also mention that the mathematical structure of the wavefunction (31) is reminiscent of the same quantity found in the nearest-neighbor approximation (see equation (21) of reference DowningBiPartite2017 ).

It was previously discussed in section 3 that the inclusion of both the resonance frequency $\omega_{0}$ and the intrasublattice term (2) in the full plasmonic Hamiltonian (2) leads to broken chiral symmetry (see the discussion after equation (25)). In mathematical terms, we indeed do not have $\mathcal{C}\mathcal{H}_{q}^{\sigma}\mathcal{C}=-\mathcal{H}_{q}^{\sigma}$ (chiral symmetry) due to $\omega_{0}\neq 0$ and $f_{q}^{\sigma}\neq 0$, where $\mathcal{C}=\mathbbm{1}_{2}\otimes\sigma_{z}$ is the $4\times 4$ generalized third Pauli matrix. Therefore, this raises the issue of the implications for the topological protection of the plasmonic edge states. In order to have a quantized Zak phase, it has been recently shown in reference vanMiert2017 that a 1D insulator is required to respect both time-reversal and inversion symmetries. The Bloch Hamiltonian (29) indeed satisfies both aforementioned criteria:

[TABLE]

where the $4\times 4$ exchange matrix is $\mathcal{I}=\sigma_{x}\otimes\sigma_{x}$, in terms of the first Pauli matrix $\sigma_{x}$. Hence the Zak phase is a meaningful quantity for the present problem, despite the long-range nature of the dipole-dipole interaction, the presence of nonresonant terms in the interacting part of the plasmonic Hamiltonian (4), and the intrinsic bosonic quasiparticle statistics.

Formally, the Zak phase is defined as the integral of the Berry connection $\mathcal{A}_{q}^{\sigma}=\mathrm{i}\langle\psi_{q\tau}^{\sigma}|\partial_{q}|\psi_{q\tau}^{\sigma}\rangle_{\mathcal{J}}$ over the first Brillouin zone,

[TABLE]

Here, the Berry connection should be calculated while taking into account the aforementioned metric $\mathcal{J}$. Together with equation (31), we find $\mathcal{A}_{q}^{\sigma}=\mathrm{i}\langle\psi_{q\tau}^{\sigma}|\mathcal{J}\partial_{q}|\psi_{q\tau}^{\sigma}\rangle=-\frac{1}{2}\partial_{q}\phi_{q}^{\sigma}$. Carrying out the integral in equation (33) with the above result, and taking care with the principal value of the phase $\phi_{q}^{\sigma}$, we find the exact result

[TABLE]

up to modulo $2\pi$, since the Berry connection $\mathcal{A}_{q}^{\sigma}$ is not gauge invariant. Thus, the dimerized chain has a trivial Zak phase ($\vartheta_{\mathrm{Z}}^{\sigma}=0$) when $d_{1}<d_{2}$ and a nontrivial one ($\vartheta_{\mathrm{Z}}^{\sigma}=\pi$) when $d_{1}>d_{2}$, irrespective of the collective mode polarization $\sigma$. There is a topological phase transition at the crossover point of the regularly-spaced chain $d_{1}=d_{2}$, where the bandgap closes (see figure 2). Equation (34) confirms the presence of plasmonic edge states in the nontrivial phase with $d_{1}>d_{2}$, corresponding to topologically-protected localized collective plasmons at the Brillouin zone edge DowningBiPartite2017 , as directly follows from the bulk-edge correspondence theorem Asboth2016 . Such behavior can be understood by a simple physical picture where plasmonic edge states are exponentially pinned at the unpaired nanoparticles located at the two edges of the chain Downing2017 .

In this section, we have analytically proved the robustness of bosonic edge states against long-range quasistatic interactions and nonresonant couplings. This general result may be important for studies of generalized Su-Schrieffer-Heeger (SSH) models, where previous works have focused on the inclusion of next- and third-nearest neighbor hopping terms only Li1990 ; Li1993 ; Li1995 ; Li2014 ; Perez2018 . Furthermore, there is a recent rise in interest in non-Hermitian SSH models Zhu2014 ; Lieu2018 ; Yin2018 ; Gong2018 ; Yao2018 ; Martinez2018 for which our results may have implications.

We should also mention that it has been shown recently via sophisticated numerical simulations on finite dimerized chains that plasmonic edge states are further robust against both retardation effects and moderate disorder Pocock2018 . We comment on how retardation affects the topological behavior of the system at the end of section 5.

5 Radiative frequency shifts of the collective plasmonic bandstructure

We now consider the effect of retardation in the dipole-dipole interaction amongst the nanoparticles in the bipartite chain on the plasmonic bandstructure. Along the lines of reference Downing2018 , we treat the light-matter coupling Hamiltonian (2) up to second order in perturbation theory, with the aim of obtaining an analytical understanding of the radiative frequency shift induced by the photonic environment.

In the long-chain limit ($\mathcal{N}\gg 1$), our perturbative calculation yields the renormalized collective plasmon dispersion relation (see appendix C for details)

[TABLE]

where the radiative shifts are given by

[TABLE]

Here, $\Theta(z)$ is the Heaviside step function and the ultraviolet cutoff $\omega_{\mathrm{c}}$ is of the order of $c/a$, corresponding to the frequency above which the dipolar approximation used in equation (2) breaks down. Notably, the frequency shift (5) is only logarithmically dependent on the cutoff frequency $\omega_{\mathrm{c}}$. This is reminiscent of the so-called Bethe logarithm in the equivalent expression for the Lamb shift in atomic physics Bethe1947 ; Milonni1994 . Also immediately apparent from equation (5) is the characteristic size dependence $\delta_{q\tau}^{\sigma}\propto a^{3}$, as is the case of the regular chain Downing2018 .

In equation (5), we have introduced the brightness factor

[TABLE]

where the phase $\phi_{q}^{\sigma}$ is defined in equation (30). This quantity measures the strength of the light-matter coupling, and arises from the first term on the right-hand side of the coupling Hamiltonian (2). Importantly, $\mathcal{B}_{q\tau}^{\sigma}$ depends on the polarization $\sigma$ and the band index $\tau$ only through the product $\tau\,\mathrm{sgn}\{\eta_{\sigma}\}$. This is exemplified in figure 3, where we show the brightness factor in the first Brillouin zone for several values of the dimerization parameter $|\epsilon|$. As can be seen from the figure, there are drastically different behaviors of $\mathcal{B}_{q\tau}^{\sigma}$ for $\tau\,\mathrm{sgn}\{\eta_{\sigma}\}=+1$ ($-1$) which corresponds to bright (dark) plasmonic bands. Therefore, for the transverse polarization ($\sigma=x,y$) the high-energy band ($\tau=+$) is bright (i.e. it significantly couples to light) and corresponds to an orientation where the dipole moments on each nanoparticle in a given dimer are in-phase ($\uparrow\uparrow$). Whereas, the low-energy band ($\tau=-$) refers to a situation where the dipole moments on each nanoparticle in a dimer are out-of-phase ($\uparrow\downarrow$), or a dark band. For the longitudinal polarization ($\sigma=z$), the aformentioned situation is inverted: $\tau=+$ ($\tau=-$) corresponds to a dark (bright) band.888Whether a band is bright or dark can in fact be easily envisioned, at least for the $q=0$ mode, from the very nature of the dipole-dipole interaction, which energetically favors an antiparallel (parallel) alignment of the dipole moments in the transverse (longitudinal) configuration. Notice that for plasmonic modes with wavevectors close to the Brillouin zone edge, this classification is less pronounced in the dimerized case ($|\epsilon|\neq 0$).

The strikingly different polarization-dependent behavior of $\mathcal{B}_{q\tau}^{\sigma}$ shown in figure 3 has important implications for the frequency shifts (5), as well as for the radiative decay rates of the system (see appendix D). In figure 4, we plot the renormalized plasmonic bandstructure (35) for both the transverse (panel (a)) and longitudinal polarizations (panel (b)) as colored lines, for nanoparticle radii $k_{0}a=0.1$ (green lines), $k_{0}a=0.2$ (blue lines), $k_{0}a=0.3$ (red lines) with $k_{0}=\omega_{0}/c$, and fixed periodicity $d=6.5a$. In the figure, the dimerization parameter is $|\epsilon|=1/13$ (cf. equation (23)) and the ultraviolet cutoff frequency is chosen to be $\omega_{\mathrm{c}}=c/a$. With increasing nanoparticle size, the radiative corrections (5) significantly modify the quasistatic result in figure 4 (see black lines) for the bright bands and are negligible for the dark bands. Thus, for the transverse polarization it is the upper-in-energy ($\tau=+$) band which experiences a sizable redshift, while for the longitudinal polarization it is the lower band ($\tau=-$).

Most noticeable from figure 4 are the presence of cusps in the dispersion relation of the upper band of the transverse modes (panel (a)). In the large-chain limit ($\mathcal{N}\gg 1$), such cusps corresponds to logarithmic singularities, see the last term in equation (5). These singularities occurs at the intersection between the plasmonic dispersion and the light cone ($\omega_{q+}^{x,y}=c|q|$), i.e. at $q\simeq\pm k_{0}$ to zero order in $\Omega/\omega_{0}\ll 1$.999Formally, there are also logarithmic singularities in the lower transverse band at $\omega_{q-}^{x,y}=c|q|$. However, this is a dark band so that these singularities are heavily suppressed by the near-zero behavior of the brightness factor $\mathcal{B}_{q\tau}^{\sigma}$ shown in figure 3. These striking features also manifest in the associated radiative damping, where they appear as discontinuities in the decay rate at the same points in wavevector space (see appendix D). Remarkably, the influence of the cusp singularities on the bandstructure leads to a band inversion for larger particles (red lines in figure 4(a)), since the $\tau=+$ band dips below the $\tau=-$ band in some regions of the first Brillouin zone. However, this band crossing phenomenon should be treated with some caution, since it arises at the limit of applicability of our dipolar approximation (cf. equation (2)).

The behavior observed in figure 4 and encapsulated in equations (35) and (5) has been recently reported by means of sophisticated numerical calculations based on solutions of the fully-retarded Maxwell equations Pocock2018 . Our elementary open quantum system approach provides a simple and transparent analytical understanding of the plasmonic bandstructure, including the intriguing effects of retardation.

Radiative shifts and radiative decay rates are intrinsically linked by the fluctuation-dissipation theorem Kubo1966 . The reporting of the exact result of the quasistatic bandstructure (19) means that the accuracy of the radiation damping decay rate, previously derived within the nearest-neighbor approximation in reference DowningBiPartite2017 , can be improved upon. In appendix D, we provide the exact expression for the radiation damping decay rate including the effects of long-range interactions.

Since the radiation damping decay rate is vanishing outside the light cone, that is for wavevectors $|q|\gtrsim k_{0}$, the edge states discussed in section 4 are dark (due to being pinned at $|q|=\pi/d$). It then follows that these edge states weakly interact with light and so should not be destroyed by retardation effects. Indeed, while the inclusion of retardation acts to further break the band reversal (particle-hole) and chiral symmetries of the collective plasmons, these broken symmetries were already shown in section 4 to be inconsequential for the existence of topologically nontrivial states. Then, as long as there is a well-defined bandgap at the Brillouin zone edge, plasmonic edges states should be detectable. This physical argument is supported by the latest numerical calculations on finite plasmonic chains with retardation Pocock2018 , where it is also shown that plasmonic edge states are robust to moderate disorder.

Finally, we emphasize that the results presented in this section arise from a perturbative calculation. Consequently, the singularity in the transverse collective plasmon dispersion, leading to a large deviation from the natural LSP frequency $\omega_{0}$, should be treated with caution. An analysis of the strong coupling regime, where the quasiparticles are plasmon-polaritons Torma2015 , is beyond the scope of this work.

6 Conclusions

We have presented a simple Hamiltonian model of collective plasmons in a dimerized chain of spherical metallic nanoparticles. Taking into account long-range quasistatic interactions, we have derived an exact expression for the plasmonic bandstructure. We have shown that the Bloch Hamiltonian obeys time-reversal and inversion symmetries, such that the Zak phase is quantized, despite the broken particle-hole and chiral symmetries of the system. We have calculated the Zak phase explicitly, finding a topologically nontrivial regime of dimerization, which proves the existence of plasmonic edge states which are robust against long-range interactions, bosonic statistics, and nonresonant couplings.

We have also taken into account retardation effects and calculated analytical expressions for the resulting radiative frequency shifts. While the dark bands are essentially unaffected by retardation, the bright bands present a size-dependent frequency shift. The latter is most pronounced for the transverse-polarized modes in the vicinity of the intersection between the quasistatic spectrum and the light cone. Otherwise, the fully-retarded plasmonic bandstructure presents only a quantitative difference with respect to the quasistatic result for small nanoparticles. Therefore, the topologically nontrivial edge states are essentially unaffected by retardation.

Acknowledgments

We thank François Gautier and Dietmar Weinmann for useful discussions. We acknowledge financial support from the French Agence Nationale de la Recherche (Project ANR-14-CE26-0005 Q-MetaMat) and from the Spanish Ministerio de Economía y Competitividad (Project FIS2015-64951-R CLAQUE).

Author contribution statement

Both authors contributed equally to this paper.

Appendix A Lattice sums $f_{q}^{\sigma}$ and $g_{q}^{\sigma}$

In this appendix we discuss the behavior of the lattice sums (16) and (17), which are both $2\pi/d$-periodic as a function of the wavevector $q$.

The even intrasublattice function $f_{q}^{\sigma}$ is plotted as a dashed line in figure 5. It is bounded by its values at the center and edge of the first Brillouin zone, $f_{0}^{\sigma}=2\zeta(3)\eta_{\sigma}$ and $f_{\pi/d}^{\sigma}=-3\zeta(3)\eta_{\sigma}/2$, respectively, where Apéry’s constant $\zeta(3)\simeq 1.20$ is defined through the Riemann zeta function $\zeta(s)=\sum_{n=1}^{\infty}n^{-s}$. In the vicinity of the center of the first Brillouin zone ($|q|d\ll 1$), we have $f_{q}^{\sigma}/\eta_{\sigma}\simeq 2\zeta(3)-(3-2\ln|q|)q^{2}/2$, while close to the Brillouin zone edge, $f_{q}^{\sigma}/\eta_{\sigma}\simeq-3\zeta(3)/2+\ln{(2)}(|q|d-\pi)^{2}$.

The absolute value of the intersublattice function $g_{q}^{\sigma}$ which enters the quasistatic spectrum (19) is also even and is plotted as solid lines in figure 5 for increasing absolute values of the dimerization parameter $\epsilon$ defined in equation (23). Note that $|g_{q}^{\sigma}|$ is an even function of the parameter $\epsilon\in[-1,+1]$. The intersublattice function (17) has the boundary values

[TABLE]

and

[TABLE]

where the Hurwitz zeta function is defined as $\zeta(s,z)=\sum_{n=0}^{\infty}(n+z)^{-s}$.

Appendix B Mapping of the quasistatic plasmonic bandstructure at the Brillouin zone edge to a quasirelativistic spectrum

Here we provide details of the expansion of the quasistatic spectrum (19) close to the edge of the first Brillouin zone, showing that it maps to the quasirelativistic spectrum (28).

To leading order in $kd\ll 1$, with $qd=kd+\pi$, we find after a lengthy but straightforward Taylor expansion, that equation (19) reduces to equation (28), with

[TABLE]

where $g_{\pi/d}^{\sigma}$ is given in equation (A), and

[TABLE]

Here we have introduced the auxiliary functions

[TABLE]

where the Hurwitz zeta function $\zeta(s,z)$ is defined in appendix A. This completes the analytical description of the quasirelativistic plasmonic dispersion at the Brillouin zone edge (28).

Appendix C Photonic-induced frequency shifts of the collective plasmons

In this appendix, we provide details of the calculation of the renormalized plasmonic dispersion (35) and the radiative frequency shifts (5) induced by the photonic bath in which the nanoparticles are embedded.

Following reference Downing2018 , we treat the plasmon-photon coupling Hamiltonian (2) up to second order in perturbation theory. For a given mode of polarization $\sigma$, band index $\tau$, and wavevector $q$, the plasmonic energy levels are

[TABLE]

where $n_{q\tau}^{\sigma}$ is a non-negative integer. Here, the zeroth-order contribution $E_{n_{q\tau}^{\sigma}}^{(0)}=n_{q\tau}^{\sigma}\hbar\omega_{q\tau}^{\sigma}$ is simply the unperturbed quantum harmonic oscillator levels corresponding to the Hamiltonian (18), where the frequency $\omega_{q\tau}^{\sigma}$ is given in equation (19).

The first-order contribution $E_{n_{q\tau}^{\sigma}}^{(1)}$ arises from the second, “diamagnetic” term on the right-hand side of equation (2). It corresponds to the process whereby a virtual photon, with wavevector $\mathbf{k}$ and transverse polarization $\hat{\lambda}_{\mathbf{k}}$, is emitted and absorbed. Notably, this self-energy term does not involve any plasmonic operators and so the plasmonic eigenstate $|n_{q\tau}^{\sigma}\rangle$ remains unchanged by the process. Explicitly, one finds

[TABLE]

which corresponds to a global energy shift, independent of the quantum number $n_{q\tau}^{\sigma}$. Since only interlevel energy differences are observable through plasmonic transitions, this seemingly divergent contribution does not lead to a renormalization of the bare collective mode resonance frequency $\omega_{q\tau}^{\sigma}$ and may be discarded.

The second-order contribution $E_{n_{q\tau}^{\sigma}}^{(2)}$ in equation (43) arises from the first, “paramagnetic” term on the right-hand side of the coupling Hamiltonian (2). It corresponds to the emission and reabsorption of virtual photons via the intermediate plasmonic states $|n_{q\tau}^{\sigma}\pm 1\rangle$. One obtains

[TABLE]

where the summation over $\mathbf{k}$ excludes the singular term for which $\omega_{\mathbf{k}}=\omega_{q\tau}^{\sigma}$. In the expression above, the array factor reads as

[TABLE]

with $k_{z}=\mathbf{k}\cdot\hat{z}$, and the phase $\phi_{q}^{\sigma}$ is defined in equation (30). In the continuum limit, where $\sum_{\mathbf{k}}\to\mathcal{P}\int\mathrm{d}^{3}\mathbf{k}\;\mathcal{V}/{(2\pi)^{3}}$ (here, $\mathcal{P}$ denotes the Cauchy principal value), the second-order correction (C) is divergent. This divergence can be regularized by introducing a physically-motivated ultraviolet cutoff $k_{\mathrm{c}}\sim 1/a$, corresponding to the wavelength below which the dipolar approximation used in equation (2) breaks down.

The renormalized frequency difference between successive plasmonic energy levels $\tilde{\omega}_{q\tau}^{\sigma}=(E_{n_{q\tau}^{\sigma}+1}-E_{n_{q\tau}^{\sigma}})/\hbar$ is then independent of the quantum number $n_{q\tau}^{\sigma}$, up to second order in perturbation theory. Hence, the radiative frequency shift (see equation (35)) is given by

[TABLE]

The identity

[TABLE]

allows one to perform the summation over photon polarization in equation (C). Transforming the wavevector summation into a principal-value integral, we obtain in spherical coordinates $(k,\theta,\varphi)$

[TABLE]

where the cutoff $k_{\mathrm{c}}>\omega_{q\tau}^{\sigma}/c$. The integral over the azimuthal angle $\varphi$ is readily evaluated with the identity

[TABLE]

Since we are working in the long-chain regime ($\mathcal{N}\gg 1$), the subsequent integral over the polar angle $\theta$ is easily obtained using the result

[TABLE]

where $\delta(z)$ is the Dirac delta function. The final integral over the radial coordinate $k$ then yields the desired result (5).

Appendix D Radiation damping decay rates

Here we present a detailed calculation of the radiation damping decay rates of the collective plasmons, including the long-ranged quasistatic dipole-dipole interaction, thereby extending the results of reference DowningBiPartite2017 which was restricted to nearest-neighbor interactions only.

The radiative damping of a $\sigma$-polarized collective mode with wavevector $q$ in the band $\tau$ arises from the decay of such a state, resulting in the spontaneous emission of photons with momentum $\mathbf{k}$ and transverse polarization $\hat{\lambda}_{\mathbf{k}}$ in the far field. The associated radiation damping decay rates $\gamma_{q\tau}^{\sigma}$ of the system follow directly from a Fermi golden rule calculation with the coupling Hamiltonian (2), yielding

[TABLE]

where the array factor $F_{\mathbf{k},q}^{-}$ is given in equation (46) and the phase $\phi_{q}^{\sigma}$ is defined in equation (30). In the continuum limit, we replace the summation over photonic wavevectors by an integral, and using equation (48), we obtain

[TABLE]

It should be noted that the $k$-integral in the above equation does not require the introduction of the ultraviolet cutoff $k_{\mathrm{c}}$, contrary to the case of the radiative frequency shifts (see equation (C)). This is because the Dirac delta function in equation (D) imposes conservation of energy, which is compatible with the dipolar approximation used in equation (2). With the help of equations (50) and (51), we then obtain in the long-chain limit ($\mathcal{N}\gg 1$) the final result

[TABLE]

where the brightness factor $\mathcal{B}_{q\tau}^{\sigma}$ is defined in equation (37). The expression above presents an $a^{3}$-dependence, as is the case for the radiative damping of a single nanoparticle $\gamma_{0}=2\omega_{0}^{4}a^{3}/3c^{3}$. Immediately apparent in the above expression is the influence of the light cone at $|q|\simeq k_{0}$ (through the Heaviside step function $\Theta(z)$) as was noticed for the associated radiative shifts (5). Clearly, the decay rate is vanishing outside the light cone ($|q|\gtrsim k_{0}$).

In figure 6, we plot the radiation damping decay rate (D) for the bright bands (solid lines) and the same quantity in the nearest-neighbor approximation given in equation (33) of reference DowningBiPartite2017 (dashed lines). Due to the brightness factor $\mathcal{B}_{q\tau}^{\sigma}$ in equation (D), the results for the dark bands are negligible in comparison to the bright bands (which display a superradiant behavior), and so we do not plot them. In the figure, we show results for increasing nanoparticle sizes $k_{0}a=0.1$ (green lines), $k_{0}a=0.2$ (blue lines), $k_{0}a=0.3$ (red lines), and fixed periodicity $d=6.5a$. The dimerization of the chain is taken to be $|\epsilon|=1/13$. Clearly, for the transverse polarization (panel (a)) the two expressions are in excellent agreement, with the dashed lines being barely noticeable. The longitudinal polarization (panel (b)) displays a very good agreement between the exact and approximate expressions, with small deviations only visible at the center of the first Brillouin zone. We can therefore conclude that the long-ranged nature of the quasistatic dipole-dipole interaction has almost no quantitative effects on the radiative damping rates of the collective plasmons.