Quantum Plasmonic Nanoantennas

TL;DR

This paper investigates quantum effects in plasmonic excitations of sodium chains, revealing a quantum-classical transition, field enhancement potential at the molecular level, and the impact of gap size on charge-transfer states in dimers, advancing ultra-small plasmonic device design.

Contribution

It introduces a detailed analysis of quantum plasmonic excitations in atomic-scale sodium chains, bridging classical and quantum models, and explores their potential for ultra-small tunable plasmonic devices.

Findings

Quantum-classical transition occurs around 10 atoms in longitudinal modes.

Large electric field enhancements are achievable at the molecular scale.

Charge-transfer states depend on gap and chain size in sodium dimers.

Abstract

We study plasmonic excitations in the limit of few electrons, in one-atom thick sodium chains, and characterize them based on collectivity. We also compare the excitations to classical localised plasmon modes and find for the longitudinal mode a quantum-classical transition around 10 atoms. The transverse mode appears at much higher energies than predicted classically for all chain lengths. The electric field enhancement is also considered which is made possible by considering the effects of electron-phonon coupling on the broadening of the electronic spectra. Large field enhancements are possible on the molecular level allowing us to consider the validity of using molecules as the ultimate small size limit of plasmonic antennas. Additionally, we consider the case of a dimer system of two sodium chains, where the gap can be considered as a picocavity, and we analyse the charge-transfer states and their dependence on the gap size as well as chain size. Our results and methods are useful for understanding and developing ultra-small, tunable and novel plasmonic devices that utilise quantum effects that could have applications in quantum optics, quantum metamaterials, cavity-quantum electrodynamics and controlling chemical reactions, as well as deepening our understanding of localised plasmons in low dimensional molecular systems.

Figures

c in the main text.

c in the main text.

Equations

$$\delta n^I(\mathbf{r},{\omega }_I)=-\frac{\delta {\rho }^I(\mathbf{r})}{i\eta }{\mathbf{\mu }}^I\cdot \mathbf{E}.$$

$$FE=\frac{|\delta \mathbf{E}(\mathbf{r},\omega )+{\mathbf{E}}_{ext}|}{|{\mathbf{E}}_{ext}|}$$

$$\left(\begin{array}{cc}\mathbf{A}& \mathbf{B}\\ {\mathbf{B}}^{\ast }& {\mathbf{A}}^{\ast }\end{array}\right)\left(\begin{array}{c}\mathbf{X}\\ \mathbf{Y}\end{array}\right)=\omega \left(\begin{array}{c}\mathbf{X}\\ \mathbf{Y}\end{array}\right)$$

$$\begin{array}{r}{A}_{ia,jb}=({ϵ}_a-{ϵ}_i){\delta }_{ij}{\delta }_{ab}+2K_{ia,jb}\\ B_{ia,jb}=2K_{ia,bj}\end{array}$$

$$\begin{array}{r}{K}_{ia,jb}=\int d\mathbf{r}d{\mathbf{r}}^{\prime }{\varphi }_i^{\ast }(\mathbf{r}){\varphi }_a(\mathbf{r})\\ \times \left(\frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}+f_{XC}(\mathbf{r},{\mathbf{r}}^{\prime })\right){\varphi }_b^{\ast }({\mathbf{r}}^{\prime }){\varphi }_j({\mathbf{r}}^{\prime })\end{array}$$

$$[{\mathbf{U}}^{†}\mathbf{T}\mathbf{V}{]}_{ij}=\sqrt{{\lambda }_i}{\delta }_{ij},$$

$$R=\frac{1}{{\sum }_i{\lambda }_i^2}$$

$${\mathbf{C}}^{+}{\mathbf{Z}}^{+}={\omega }^2{\mathbf{Z}}^{+}$$

$$\begin{array}{r}{\mathbf{C}}^{+}=(\mathbf{A}-\mathbf{B}{)}^{\frac{1}{2}}(\mathbf{A}+\mathbf{B})(\mathbf{A}-\mathbf{B}{)}^{\frac{1}{2}}\\ {\mathbf{Z}}^{+}=\sqrt{\omega }(\mathbf{A}-\mathbf{B}{)}^{-\frac{1}{2}}(\mathbf{X}+\mathbf{Y})\end{array}$$

$${\mathbf{C}}^{-}{\mathbf{Z}}^{-}={\omega }^2{\mathbf{Z}}^{-}$$

$$\begin{array}{r}{\mathbf{C}}^{-}=(\mathbf{A}+\mathbf{B}{)}^{\frac{1}{2}}(\mathbf{A}-\mathbf{B})(\mathbf{A}+\mathbf{B}{)}^{\frac{1}{2}}\\ {\mathbf{Z}}^{-}=\sqrt{\omega }(\mathbf{A}+\mathbf{B}{)}^{-\frac{1}{2}}(\mathbf{X}-\mathbf{Y}).\end{array}$$

$${|{\mathbf{Z}}^{\pm }|}^2=1$$

$$|\mathbf{X}{|}^2-|\mathbf{Y}{|}^2=1.$$

$${\mathbf{Z}}^{-}=\frac{(\mathbf{A}+\mathbf{B}{)}^{\frac{1}{2}}(\mathbf{A}-\mathbf{B}{)}^{\frac{1}{2}}}{\omega }{\mathbf{Z}}^{+}$$

$$\begin{array}{r}\mathbf{X}={\mathbf{D}}^{-}{\mathbf{Z}}^{+}+{\mathbf{D}}^{+}{\mathbf{Z}}^{-}\\ \mathbf{Y}={\mathbf{D}}^{-}{\mathbf{Z}}^{+}-{\mathbf{D}}^{+}{\mathbf{Z}}^{-}\end{array}$$

$$\begin{array}{r}{\mathbf{D}}^{+}=\frac{(\mathbf{A}+\mathbf{B})}{2\sqrt{\omega }}\\ {\mathbf{D}}^{-}=\frac{(\mathbf{A}-\mathbf{B})}{2\sqrt{\omega }}\end{array}$$

$$R=R_X^{({\sum }_{ia}{X}_{ia}^2)}{R}_Y^{({\sum }_{ia}{Y}_{ia}^2)}.$$

$$\Re \left(\delta P_{ia\sigma }(\omega )\right)=\sum _{jb\tau }\frac{1}{{\omega }^2-C_{ia\sigma ,jb\tau }}{v}_{jb\tau }^{appl}(\omega )$$

$$\begin{array}{r}\Re \left(\delta P_{ia\sigma }(\omega )\right)=\sum _I\frac{1}{{\omega }^2-{\omega }_I^2}\\ \times \sum _{jb\tau }\sqrt{{\omega }_{ai\sigma }}{Z}_{ia\sigma }^I{\left({Z_{jb\tau }}^I\right)}^{\ast }\sqrt{{\omega }_{bj\tau }}{v}_{jb\tau }^{appl}(\omega ).\end{array}$$

$$\delta n(\mathbf{r},\omega )=\sum _{ia\sigma }{\varphi }_{i\sigma }(\mathbf{r}){\varphi }_{a\sigma }(\mathbf{r})\delta P_{ia\sigma }(\omega )$$

$$\begin{array}{r}\delta n(\mathbf{r},\omega )=\sum _I\frac{2}{{\omega }^2-{\omega }_I^2}\\ \times \sum _{ia\sigma ,jb\tau }{\varphi }_{i\sigma }(\mathbf{r}){\varphi }_{a\sigma }(\mathbf{r})\sqrt{{\omega }_{ai\sigma }}{Z}_{ia\sigma }^I{\left({Z_{jb\tau }}^I\right)}^{\ast }\sqrt{{\omega }_{bj\tau }}{\mathbf{r}}_{jb\tau }\cdot \mathbf{E}(\omega ).\end{array}$$

$$\begin{array}{r}\delta n(\mathbf{r},\omega )=\sum _I\frac{1}{{\omega }_I}\left(\frac{1}{\omega -{\omega }_I+i\eta /2}-\frac{1}{\omega +{\omega }_I+i\eta /2}\right)\\ \times \sum _{ia\sigma ,jb\tau }{\varphi }_{i\sigma }(\mathbf{r}){\varphi }_{a\sigma }(\mathbf{r})\sqrt{{\omega }_{ai\sigma }}{Z}_{ia\sigma }^I{\left({Z_{jb\tau }}^I\right)}^{\ast }\sqrt{{\omega }_{bj\tau }}{\mathbf{r}}_{jb\tau }\cdot \mathbf{E}(\omega ).\end{array}$$

$$\delta n(\mathbf{r},\omega \approx {\omega }_I)\approx -\frac{1}{i\eta }\delta {\rho }^I(\mathbf{r}){\mathbf{\mu }}^I\cdot \mathbf{E}.$$

$$\begin{array}{r}\delta {\rho }^I(\mathbf{r})=⟨{\Psi }_0⟩\sum _i^N\delta (\mathbf{r}-{\mathbf{r}}_i){\Psi }_I\\ =\sqrt{2}\sum _{ia}{\varphi }_i(\mathbf{r}){\varphi }_a(\mathbf{r})\left(X_{ia}^I+Y_{ia}^I\right)\end{array}$$

$$\delta {\mathbf{\mu }}^I(\mathbf{r})=⟨{\Psi }_0⟩\hat{\mathbf{\mu }}{\Psi }_I=\sqrt{2}\sum _{ia}{\hat{\mathbf{\mu }}}_{ia}\left(X_{ia}^I+Y_{ia}^I\right)$$

$${ϵ}_m+L_i\left(ϵ(\omega )-{ϵ}_m\right)=0$$

$$\begin{array}{r}{L}_a=\frac{1-e^2}{e^2}\left(\frac{1}{2e}\log (\frac{1+e}{1-e})-1\right)\\ e^2=1-\frac{b^2}{a^2}=1-\frac{c^2}{a^2}\end{array}$$

$${\omega }_p^2=\frac{n{e}^2}{{ϵ}_{0}m}=\frac{N{e}^2}{V{ϵ}_{0}m}=\frac{3N{e}^2}{4\pi abc{ϵ}_{0}m}=\frac{3N{e}^2}{4\pi (\Delta +(N-1)d/2){\Delta }^2{ϵ}_{0}m}$$

$$L_a+L_b+L_c=L_a+2L_b=1.$$

$${\alpha }_i=4\pi abc\frac{ϵ(\omega )-{ϵ}_m}{{ϵ}_m+L_i\left(ϵ(\omega )-{ϵ}_m\right)}$$

Full text

I Introduction

Localised surface plasmons (LSPs) are collective oscillations in finite free-carrier systems driven by electromagnetic fields which dominate the optical response of metals near the plasmon frequency Maier (2007). Classically, they can be understood as a resonance arising from the restoring force of induced surface charges and, as such, are highly sensitive to the geometry and surrounding environment and lead to strong electric field confinement and enhancement which results in a wealth of potential applications such as sensing Homola et al. (1999), metamaterials/metasurfaces Boltasseva and Atwater (2011), colour generation Kristensen et al. (2016), controlling chemical reactions Yan et al. (2016) and energy conversion Atwater and Polman (2010).

In quantum plasmonics Fitzgerald et al. (2016), one aims to describe how the constituent electrons behave under a perturbation to form a collective mode. This is an important task in modern plasmonics as recent state-of-the-art experiments are exploring size regimes, such as extremely small metallic nanoparticles (MNPs) Scholl et al. (2012) and sub-nanometer gaps, Ciracì et al. (2012); Savage et al. (2012); Scholl et al. (2013, 2015) where classical laws, which have been used successfully in the vast majority of plasmonic research done to date, will fail. In fact, even quasi-quantum mechanical theories, like the hydrodynamic model Raza et al. (2011); Mortensen et al. (2014), may not be adequate for quantitative predictions Stella et al. (2013). Instead the electrons must be treated using ab initio methods such as time dependent density functional theory (TDDFT). This will automatically take into account modification of the optical response from electron spill out, energy level discreteness and non-locality: which arises as a consequence of the finite extent of the electronic wavefunction. It can also be important to explicitly include the underlying atomic structure rather than relying on simplifications such as the jellium model Zhang et al. (2014). In this picture, MNPs behave as large molecules and the plasmon, which captures much of the spectral weight in the extinction cross-section, is built up out of transitions between discrete energy levels mediated by the Coulomb interaction.

In this context it is important to have clear and robust methods to analyse the results from calculations which may return a large amount of data and many excitations over the frequency range of interest. The surface of the MNP breaks translational invariance and allows mixing of plasmons and single particle excitations (Landau damping) meaning identification of peaks can be difficult. Often any peaks with a large dipole strength are claimed to be plasmonic Gao and Yuan (2005); Yan et al. (2007); Yan and Gao (2008) using the argument that this indicates collectivity. This can be erroneous as the dipole strength function, as the name suggests, simply gives a measure of the dipole moment and hence how strongly the excitation interacts with light, a single-electron hole pair separated over a large distance can give a large dipole moment but it does not seem sensible to describe it as plasmonic! Thus to identify plasmons in small clusters more sophisticated methods of identification are needed. Examples in the literature include analysis of the electron population dynamics Townsend and Bryant (2014); Ma et al. (2015), the Coulomb scaling method Bernadotte et al. (2013) and the plasmonicity index Bursi et al. (2016). Collectivity is one of the most intuitive characteristics of a plasmon, in this work we focus on measuring collectivity to identify plasmonic excitations. Using excited state structural analysis to quantify the collectivity of an excitation, which can be obtained from the standard output of any electronic structure code (as long as the transition density matrix can be calculated), we define a new index which is the product of the collectivity index Martin (2003); Plasser et al. (2014) and the dipole strength. This provides a simple and convenient method to quantify how plasmon-like a molecular excitation is.

One of the most widely studied plasmonics systems are metallic nanoantennas Muskens et al. (2007); Biagioni et al. (2012); Giannini et al. (2011) where oscillation of free carriers in a confined space allows the reception or emission of electromagnetic radiation. To achieve a large optical response in the visible one needs to tune the resonant frequency by shrinking the antenna down to smaller sizes. Nanorods, in particular, have been studied in depth and, within the quasi-static limit ($\lambda\gg L,R$, where $L$ is the length of the nanorod and $R$ is the radius), the plasmon resonance depends only on the aspect ratio $L/R$ Bryant et al. (2008). They are also one of the simplest anisotropic structures which lead to strong and tunable resonances. Single atomic thick atomic chains, which have become a standard ‘workhorse’ of quantum plasmonics Gao and Yuan (2005); Yan et al. (2007); Yan and Gao (2008); Yasuike et al. (2011); Bernadotte et al. (2013); Piccini et al. (2013), represent the ultimate aspect ratio limit of nanorods. Their small size could lead to added tunability and novel quantum plasmonic behaviour as well as aid in our understanding of how plasmons arise in larger systems within a quantum picture. Nanorods have also been explored within the jellium approximation using TDDFT Zuloaga et al. (2010). Using linear response TDDFT (a.k.a. the Casida method) on the open source code OCTOPUS Andrade et al. (2015) we find that while the electronic excitations of such chains may share some characteristics with classical plasmons, the longitudinal mode in particular can be well described by a classical ellipsoid model above a threshold of about 10 atoms. In other respects, they differ considerably with the transverse mode appearing at much higher energies. This reveals that the atomic chain is a very quantum system no matter how long the chain is made.

One of the most crucial physical quantities to consider in plasmonics is the field enhancement. Unfortunately this is difficult to calculate in ab initio methods as the field enhancement depends on the damping and, in ultra-small quantum systems, the damping is not in general known whereas in bulk systems typically an experimental value is taken. This is particularly a problem when calculating the contributions to the Raman signal in surface-enhanced Raman scattering (SERS) for molecules close to metal clusters where both the field enhancement and chemical interactions (charge transfer between the molecule and metal) depend strongly on the linewidths of excited states Jensen et al. (2008). Therefore, to analyse the atomic chain in terms of a plasmonic component we need an estimate of how phonons couple to electronic excitations. In a separate calculation we have used density-functional perturbation theory (DFPT) Baroni et al. (2001) to calculate the normal modes for sodium chains and use the highest energy mode as an estimate of the upper bound for the line-broadening of the molecular excitations. We find that high field enhancement values, comparable to what is achievable in much larger plasmonic systems, are possible in atomic chains opening up the possibility of designing quantum antennas.

Another commonly explored plasmonic system is the dimer system which, for instance, has been explored for metallic nanospheres Hao and Schatz (2004), nanorings Tsai et al. (2012) and nanorods Jain et al. (2006); Biagioni et al. (2012). For separation distances above about $1\ nm$ the results can be understood classically using the plasmon hybridisation model Nordlander et al. (2004). The interference of the two LSPs leads to an extreme field enhancement and concentration in the gap region. For gap sizes below $1\ nm$ non-local effects start to play a role and limits the achievable field enhancement Ciracì et al. (2012). For even smaller gap sizes the electronic wavefunctions begin to overlap and electrons can tunnel across the junction forming a conductive bridge between the two MNPs and a new mode called the charge transfer plasmon (CTP) can appear Zuloaga et al. (2009), full quantum models become essential to explain this phenomenon Zhu et al. (2016). There have been a number of state-of-the-art experiments exploring nanogap dimer systems Ciracì et al. (2012); Savage et al. (2012); Scholl et al. (2013, 2015) which have been considered in the context of SERS Zhu and Crozier (2014), third harmonic generation Hajisalem et al. (2014), strong coupling Chikkaraddy et al. (2016), optical rectification Ward et al. (2010) and molecular tunnel junctions Tan et al. (2014). Remarkably, given the intense interest in this area of research, there has been no excited state structural analysis of gap plasmons in ultra-small systems that we are aware of. To rectify this, we have analysed the case of a sodium chain with a gap opened in the middle and have explored how this modifies the electronic excitations and consequent plasmonic behaviour. We find that the CTPs in atomic chains can be understood as combination of a few electron-hole pairs around the HOMO-LUMO and are not very plasmon-like. We also find, analogous to larger systems, bonding dipole plasmons (BDPs) which correspond to large field enhancements in the gap region. Such knowledge may aid understanding of tunnelling in larger systems which would be difficult to model using ab initio methods due to the large number of electrons. There are also applications in picocavity physics where the small gap region and high quality factor could create strongly mixed light-matter states Chikkaraddy et al. (2016).

II Results and discussions

II.1 The single sodium chain

To investigate plasmons in small molecular systems we have studied, what has become a prototypical plasmonic system to study with TDDFT, the sodium chain Gao and Yuan (2005); Yan et al. (2007); Yan and Gao (2008); Yasuike et al. (2011); Bernadotte et al. (2013). We now discuss the results for a single sodium chain of 10 atoms with an interatomic distance of $d=3.08\ \mbox{\AA}$ and aligned along the $x$ axis. The calculation techniques are discussed in appendix A and in appendix B we show how to measure the collectivity. In figure 1a we present the results showing the dipole strength function. The red line shows the full Casida spectra and the blue dashed line, for comparison, shows the spectra with interactions switched off (excitations are given by the difference of Kohn-Sham (KS) orbitals). From the transition densities (density plots in figure 1a) one can clearly identify the longitudinal (L), transverse end (TE) and transverse central (TC) modes that were first identified by Yan et al Yan et al. (2007) and which have also been studied for other metal chains such as gold Piccini et al. (2013). Similar results are also seen for graphene nanoribbons Cocchi et al. (2012). Both the L and TC modes have a dipolar character, along the long (x direction) and short lengths (y and z direction) of the chain respectively. For the widths of the peaks we use a broadening of $\eta=0.0228\ eV$ which is obtained by considering phonon-electron coupling and is discussed in the next section. The area under the curve must give the number of electrons according to the f-sum rule, we find that considering the spectrum up to $6\ eV$ as well as $105$ unoccupied states satisfies the sum rule to over $96\%$ which confirms the quality of the basis set truncation. The transverse modes are naturally higher in energy due to the greater electron confinement. The L mode is the lowest optically active mode, as is found for a classical nanorod Bryant et al. (2008). As chain length is increased the L and TC mode grow in size whilst the TE mode saturates (this mode seems to be a consequence of the atomic nature of the ends of the chain, we have found it does not appear in simplified models where the atomic structure is ignored. In fact, it can be identified as an end mode), this accumulation in oscillator strength does not necessarily mean an increase in electron collectivity. To explore the nature of these excitations we analyse the Casida eigenvectors to see which transitions between KS orbitals make up each excitation. We find that the L mode is predominately made up of the HOMO-LUMO gap transition (states $5\rightarrow 6$ in this case). We have checked this by repeating the Casida calculation but freezing all other orbitals apart from states $5$ and $6$ which recreates the L excitation (apart from a negligible depolarization shift).

The TC excitation has a much more complicated structure than the L mode and is made up of multiple KS orbital transitions. This already hints that the TC mode is much more collective and is confirmed by the collectivity index of $4.6$: this indicates that virtually all the electrons take part in this excitation. The explanation for this collectivity is the large number of transitions with a change in the quantum number associated with the y and z direction. These have the same symmetry and are close in energy (because the length is much longer than the radius, states with different quantum numbers associated in the x direction are close in energy) which has been identified as a condition for single-particle excitations to couple Guidez and Aikens (2014). This effect is increased by a degeneracy of two associated with the equivalence of the $y$ and $z$ direction. In figure 1b we illustrate the different KS transitions that make up the L and TC excitations 111Note that there are contributions from higher energy transitions to unbound states for the TC excitation not shown in figure 1b..

We find, in agreement with Yasuike et al Yasuike et al. (2011), that the L mode has a collectivity index of $1$ which led them to conclude that the L mode is a long-range charge transfer excitation rather than a plasmon. In fact, the L mode is very similar to high absorption excitations seen in conjugated polymers with high persistence lengths Vezie et al. (2016). What complicates matters is that the L mode has been identified as plasmonic in previous works because of it scaling with the Coulomb strength, Bernadotte et al. (2013) and a large plasmonicity index Bursi et al. (2016), furthermore the L mode is the second lowest energy excitation (rather than first which is a dark mode and so can’t be seen in figure 1a) which hints that the Coulomb interaction does play a role but is not able to couple different transitions Piccini et al. (2013). Thus the L mode should be identified as a ‘protoplasmon’ - a special type of mode, peculiar to the extreme geometry of the atomic chain, that exhibits a duality of plasmonic and low collectivity behaviour. Also shown in figure 1a is our plasmon index indicated by the vertical black lines. For convenience we divide by the number of electrons multiplied by $\pi\eta$ (which comes from the normalisation of the Lorentzian which is convolved with the oscillator strengths to obtain a continuous spectra). We find it identifies the TC mode as the most plasmon-like excitation due to its combined dipole strength and collectivity. The TE and L mode are identified as moderately plasmonic. Interestingly the second order L mode at $2.09\ eV$ has a much larger collectivity of $3.72$ although it is not identified as very plasmonic due to a small dipole moment.

In figure 2a we show the results of a sweep over different chain lengths, from 2 to 24 atoms in steps of 2 so all chains have an even number of valence electrons (so the system is always spin unpolarized). How the excitation energy depends on the chain length can help aid identification. The L mode redshifts with increasing length which can be understood either via a decrease in the restoring force on the plasmon mode or via a decrease in the HOMO-LUMO gap. The TE and TC modes, once formed for chain lengths above about 6 atoms, do not shift as the chain width is kept constant. We find that the L mode does not become more collective as the chain size is increased. The largest chain we checked was a \ceNa60 chain and again found a collectivity of 1 for the L mode. Interestingly it is now the 5th lowest energy excitation again indicating the role of the Coulomb interaction. The low energy absorption spectrum is shown in appendix D revealing many higher order longitudinal modes which have increasingly higher collectivities.

To further understand the nature of the excitations we try and identify them with classical modes, we fit the excitation energies with the two plasmon modes of a prolate spheroid which can be expressed analytically (see supplementary material). We treat the plasma frequency as dependent on the minor semiaxis $\Delta$ via the electron concentration, this reduces the number of fitting parameters and gives a more realistic test of the classical model 222Similar analysis was performed in reference Yan and Gao (2008) they found a surprisingly good agreement with the classical ellipsoid model by fitting both the shorter semiaxis (which will depend on the amount of electron spill-out and we will label as $\Delta$) and using the bulk plasma frequency of $3.83\ eV$. In fact, whilst the fitting to the L mode is indeed good they seem to have fitted to the TE rather than the TC mode.. We calculate the L mode energy for chains up to 40 atoms and for the largest chain size (which can be expected to behave the most classical) we fit for $\Delta$ and get a value of $1.45\times d$ (see figure 2b). For larger chains the classical model works well but below about 10 atoms the model starts to fail and predicts energies too low which can be attributed to quantum effects as well as a failure of approximating the system as an ellipsoid. Interestingly for the above mentioned fitting parameter the TC mode energies predicted are much too low (about an $eV$) highlighting that the classical model breaks down for such small length scales, the increased energy is probably down to the strong quantum confinement. Also the classical model predicts the wrong behaviour for the (small) shift in energy of the TC mode with chain length, the classical model predicts a slight blueshift with increasing chain length for small chain sizes before becoming constant whilst the quantum calculations show a redshift. In figure 2c we show the corresponding, classically calculated, absorption cross section (we show only up to \ceNa24). We find that the relative cross section amplitudes in the classical calculation are similar to the Casida calculation showing that the identification of the TC mode with the classical transverse mode is correct and they do show some qualitative similarities. The radius of the corresponding classical ellipsoid is quite large and the electron density is practically zero at this point, this reveals why the plasma frequency is found to be much smaller than the 3D bulk value. We have also tried other methods for fitting which are shown in appendix E .

II.2 Field enhancements

The local field enhancement of an excitation, as induced by an electric field $\bm{E}_{0}e^{-i\omega_{n}t}$, can be found from the induced density $\delta n$ (not to be confused with the transition density $\delta\rho$ - although they are closely connected) which can be obtained from the real part of the density matrix Casida (1995) and, close to an excitation, written in the following revealing form for the Ith excitation

[TABLE]

The response is purely imaginary as the electrons oscillate out of phase at the resonance. In appendix C we show a derivation of this equation, it has also been derived, in a different manner, by Cocchi et al Cocchi et al. (2012). The field enhancement

[TABLE]

can be obtained within the electrostatic approximation. Clearly the induced field depends crucially on the magnitude of the damping for the excitation. Using DFPT we calculate the vibrational normal modes of the \ceNa10 chain. There are a total of $3N=30$ modes with only the top $9$ modes corresponding to non-zero energies. In figure 3a we show the excitation energies and the corresponding motion of the sodium atoms which all turn out to correspond to longitudinal motion. We take the largest energy of $0.0228\ eV$ and use it as the linewidth in the Casida calculations. This is the simplest method to include electron-phonon coupling and should give an upper-bound on the true value. We note that this is a particularly simple method and assumes the phonon is a perturbation on the electronic spectra. For larger plasmon-phonon coupling there will be richer spectral features as recently shown for aromatic hydrocarbons Cui et al. (2016). For longer chains the L mode energy reduces and the phonon energy increases as shown in figure 3a which suggest for certain chain lengths there should be a strong-coupling regime. In figures 3b and 3c the field enhancement in $\log_{10}$ scale is shown for the L and TC mode respectively. We find the maximum field enhancement to be $80.5$ for the L mode and $78.3$ for the TC mode (which corresponds to intensity enhancements of $6400$), these values are comparable to the field enhancements in larger nanostructures. The electric field is qualitatively very similar to that seen for a classical antenna and shows that the modes will couple to the far field. From a practical point of view, whilst the maximum field enhancements are attractive, the sharp fall off of the induced field, due to the small number of electrons, may limit potential plasmonic applications. Fortunately this does lead to extremely large field gradients that could be useful in exciting higher order multipole excitations in small molecules and atoms, for example the TC mode has a maximum gradient enhancement $\nabla|\bm{\delta E}|/|\bm{E}_{ext}|$ of $17.1\ a_{0}^{-1}$ along the z direction for the slice $x=y=0$. If one makes do with more modest field enhancement then molecular antennas may be useful for inducing large local fields: for the TC mode field enhancements of around 5 are achievable at about $2.2$ bond lengths away which hints it might be possible to alter the optical density of states with minimal electronic coupling between the chain and a small absorbing molecule.

II.3 Double sodium chains

Next we consider separating an atomic chain along the longitudinal direction. First we consider taking a \ceNa20 chain and separating it in the middle, figure 4a shows the low energy spectra as the gap length is changed as we separate the gap ($\delta$) up to $1.4d$, we resist going further as the LDA is not suitable to consider long range van der Waals interactions which will dominate. We find that as the gap is opened there is generation of a number of new modes and we find a qualitatively similar mode to the low energy charge-transfer plasmons that has been explored heavily in the context of nearly touching dimers Zhu et al. (2016). The non-interacting electronic spectra can be understood as a splitting of energy levels due to the gap and the creation of symmetric and anti-symmetric pairs of eigenstates. The transverse modes (for clarity not shown in figure 4a, the full spectra up to $4\ eV$ is shown in appendix F) are relatively unchanged via the presence of a gap as is the case for classical nanorod dimers Jain et al. (2006), we do though observe some interference effects between the TC and TE modes at small gap sizes. There is also the formation of a number of new longitudinal modes for gap sizes $0.4d\rightarrow 1.4d$, which we will label as ‘tunnelling’ or ‘gap’ modes as they involve the transfer of charge across the gap, due to electron tunnelling, every optical cycle. We find that they tend to have a low collectivity compared to the transverse modes, but crucially they have a larger collectivity than the L mode of a single chain indicating that they arise as a consequence of multiple electron-hole excitations. For larger separation distances ($>1.2d$) the low energy absorption is dominated by the longitudinal BDP which originates from the electrostatic coupling of the two chains and is red-shifted from the mode that would be present with no coupling. To identify the peaks from the complicated spectra we will take a practical definition of a CTP as involving each chain having a different sign of charge and the intensity should go to zero for large gap sizes. The BDP should have a dipole charge distribution for each chain and the intensity will not go to zero for large distances but rather merge with the L mode of a single \ceNa10 chain. Identification is easiest using the induced potential (see appendix G). For small gap sizes, where electron tunnelling is most important, two CTPs dominate the spectra. The lowest energy excitation we label CTP1 which redshifts and loses intensity as the gap increases and becomes negligible above a gap size of about $1d$. There is also another CTP which we label CTP2 which also loses intensity as the gap sizes increase. Both CTP modes merge with the \ceNa20 L mode as $\delta\rightarrow 0$. We identify a single BDP that starts to become apparent about $0.7d$ and grows in intensity, it is the highest energy tunnelling mode and is close in energy to the \ceNa10 L mode. In between the CTPs and the BDP there are 3 modes which we believe are best explained as hybrid modes: they have a similar charge distribution as a BDP although the change in the potential over the gap region is smaller which is a consequence of charge transfer neutralising the induced charge setup across the gap, they also die out as the gap length is increased. This is a pleasing example of how plasmonic concepts and analogies can be used to simplify analysis of an electronic excitation spectra, as compared to analysis of the orbitals, which in this case is complicated and not very illuminating.

To aid identification of the modes we have also considered holding the gap distance fixed and sweeping the chain length: see figure 4b. We consider only the lower energy L modes allowing us to consider chain sizes up to 40 atoms, we only look at atomic numbers that are multiples of 4 to make sure that the separated chains have an even number of valence electrons each to ensure the validity of unpolarized spin calculations. We find that the CTPs are present for very small chains highlighting that they do not involve many electrons and originate from the interaction of electronic wavefunctions near the HOMO-LUMO gap. Crucially the L modes do not become more plasmonic as the chain length is increased.

Next we look in detail at the case of a gap size of $1d$ (which is equivalent to removing an atom from the middle of a 21 atom chain) as shown in figure 5 where the low energy spectra and the field enhancement of 4 of the tunnelling modes is shown. We find that at this distance our plasmon index predicts the hybrid mode and BDP to be much more plasmonic excitations than CTP1 and CTP2, this is backed up by the calculated field enhancement which reaches maximum values of $38.2$ and $80.1$ for hybrid and BDP modes respectively compared to $18$ and $10.6$ for the CTP1 and CTP2 (see figure 5b). The large plasmonicity for the BDP comes from its large dipole moment rather than a large collectivity, although of the 4 modes analysed it has the largest collectivity of $3.29$ as compared to values close to $2$ for the other 3 modes. The BDP mode is also by far the most interesting in terms of potential applications as it exhibits the largest field enhancement in the gap region where potentially a small atom could be placed and its emission and absorption properties modified. Thus atomic chains could potentially be used as a picocavity and modify light-matter interactions on the smallest possible scale.

III Conclusions

In conclusion we have used TDDFT to study the optical response of single-atom-thick atomic chains. We have identified plasmon modes using a new measure of ‘plasmonicity’ based on excited state structural analysis and have found that there are two modes that correspond to the classical plasmon modes but the transverse mode is blueshifted because of the extreme thinness. For the \ceNa10 chain we have used DFPT, in a separate calculation, to calculate the phonon modes and have used the highest energy normal mode to estimate the upper bound of the linewidth and hence calculate the field enhancement to test the suitability of the atomic chain as a quantum antenna. We found that high field enhancements are achievable and the rapid spatial decay of the field away from the molecule, while it may limit some potential application, results in large field gradients. Finally, we have considered an atomic chain dimer system and have explored the spectra and excited-state structure of low energy tunnelling modes as a function of gap size and chain size. We find that they are not very collective but do bear some similarities to the charge transfer plasmons of larger dimer systems in their spectral behaviour and charge distribution. A large field enhancement can be achieved in the gap region in analogy to dimers in classical plasmonics. This work shows there is pleasing synergy between the ideas of plasmonics and electronic structure which will aid development of new quantum nano-optic devices.

Acknowledgements.

The work of J. M. Fitzgerald was supported under a studentship from the Imperial College London funded by the EPSRC Grant 1580548.

Appendix A Methods

We used the open source real-space code OCTOPUS Andrade et al. (2015). The sodium atoms are described by a norm-conserving Troullier-Martins pseudopotential Troullier and Martins (1991) with only the 3s electrons explicitly modelled in the calculation. The (modified) local density approximation (LDA) Perdew and Zunger (1981) was used for the exchange-correlation potential for both the ground-state and for excited-state calculations, this is known to work well for sodium. We use the simplest exchange-correlation functional as it is well known that the existence of plasmons depends only on the Coulomb interaction and that exchange-correlation effects only provide a small correction (a small red shift) to the plasmon energy. All fields, such as the KS wavefunctions and the density, are represented on equally spaced real-space grids, the simulation box is built via a union of spheres centred on each atom. We used a radius of $12\ \mbox{\AA}$ and a grid spacing of $0.3\ \mbox{\AA}$ which gives a good compromise between accuracy and speed. For the spacing of the sodium atoms we use the interatomic distance of \ceNa2 dimers of $d=3.08\ \mbox{\AA}$ Vasiliev et al. (2002), no attempt was made of structural optimisation of any of the structures.

The electronic spectra is calculated using the Casida method Casida (1995); Jamorski et al. (1996) which is a representation of linear response theory in a configuration (particle-hole) space and is similar in spirit to the matrix formulation of RPA Bertsch and Broglia (1994). It involves solving a (pseudo) eigenvalue equation to find the electronic eigenmodes of the system

[TABLE]

where the elements of the matrices are

[TABLE]

and $K$ is the two-body Coulomb and exchange-correlation integral constructed with the KS eigenstates within the adiabatic approximation

[TABLE]

(we use the standard convention that $i,j$ and $a,b$ represent occupied and unoccupied orbitals respectively. The number of occupied and unoccupied orbitals will be given by $N_{o}$ and $N_{u}$). The $N_{o}N_{u}$ eigenvalues give the excitation energies and the eigenvectors can calculate physically relevant quantities like the transition density $\delta\rho$, dipole moments $\delta\bm{\mu}$ and dipole strength function.

It is well known that the one-particle transition density matrix $T_{ia}$ encodes, in a convenient and compact fashion, information on electronic transitions from the ground state to the excited states. It is in general non-symmetric so cannot be diagonalised, instead one can perform a singular value decomposition to express it in natural transition orbitals (NTOs) in which it will be diagonal Martin (2003)

[TABLE]

where $\bm{U}$ and $\bm{V}$ are unitary transformations from the ground state orbitals to the NTOs for the occupied and virtual orbitals respectively. The eigenvalue $\lambda_{i}$ represents the contribution of the ith electron-hole pair to a particular excitation and has the normalisation condition $\sum_{i}\lambda_{i}=1$. In the context of plasmonics it is useful to define a collectivity index (also called the inverse participation ratio) Martin (2003); Plasser et al. (2014)

[TABLE]

The collectivity index of an excitation measures the number of electron-hole pairs that contribute: the higher the value the more collective the excitation. If there is no electronic correlation, then each excitation represents a single transition and the index is equal to 1. Thus we have a quantitative measure of an excitations collectivity which will clearly be useful in identifying which excitations are plasmon-like. But on its own it is not a useful measure as many collective excitations are dipole forbidden and will not be excited by light, a useful plasmonic mode will couple strongly to light. With this in mind we introduce a new plasmon index given by the product of the collectivity and the dipole strength (of course one is not always only interested in dipole modes, for instance quadrupole modes could be of interest - in which case our index could be adjusted accordingly).

One can calculate the KS transition density matrix from a Casida calculation allowing the collectivity of an excitation to be calculated. Any excitation can be written as a superposition of p-h states with the eigenvectors of equation 3 acting as coefficients describing the contribution of each possible electron-hole transition. This clearly suggests that the Casida eigenvectors must contain information on the collectivity of an excitation. For practical considerations the Casida calculations are not solved in the form of equation 3 and direct access to the vectors $X$ and $Y$ is not possible in OCTOPUS. Also the presence of de-excitations (described by the $\bm{Y}$ eigenvector) complicates calculating the transition density matrix (within the Tamm-Dancoff approximation, where there is no coupling between the hole-particle and hole-particle correlations, it is possible to simply identify that $T_{ia}=X_{ia}$ Etienne (2015)). In appendix B the necessary post-analysis of typical TDDFT results is shown.

Phonon spectra and normal modes are obtained using DFPT Baroni et al. (2001) as implemented in Quantum Espresso code Giannozzi et al. (2009). We extensively employed DFPT for calculation of lattice dynamics and phonon energies of crystals Azadi and Foulkes (2013); Singh et al. (2014); Azadi et al. (2014) and molecules Azadi and Cohen (2016). As for the TDDFT calculations, we used the Perdew-Zunger Perdew and Zunger (1981) parametrization of the local-density-approximation (LDA) for the exchange-correlation functional. In our density functional calculations, a non-relativistic norm-conserving pseudopotential is used on an which 3s orbital with an occupation of 1 is determined as the non-local channel. Electronic structure results are calculated using the energy cutoff $60\ Ry$. For calculations on isolated one-dimensional Na-chain, we employ open boundary conditions. We performed this by locating the Na-chain in a large simulation box and truncating the Coulomb interaction at long distances to eliminate electrostatic interactions between periodic images.

Appendix B Calculation of the collectivity

Typically a slight variation (valid for real wavefunctions) of the Casida equation (equation 3) is solved in electronic structure codes

[TABLE]

where

[TABLE]

and $\bm{X}$, $\bm{Y}$, $\bm{A}$ and $\bm{B}$ were defined in the main text. This is a convenient transformation as it reduces the dimensionality of the problem by a half and one needs only deal with one set of eigenvectors when calculating physical quantities such as the transition density and dipole oscillator strength. It is also now an Hermitian eigenvalue problem. What is often overlooked is that there is an alternative and equivalent form of these equations Luzanov and Zhikol (2010)

[TABLE]

where

[TABLE]

Both eigenvectors are normalised to one

[TABLE]

which is consequence of the more unusual non-euclidean normalisation condition of the Casida eigenvectors

[TABLE]

The two sets of eigenvectors are related by the equation

[TABLE]

where the square roots should be understood as principal square roots of matrices. Thus if one has access to the two-body matrix elements $K_{ia,jb}$ then it is simple to construct both matrices $\bm{C}^{\pm}$ and then solving a simple Hermitian eigenvalue equation to obtain $\bm{Z}^{\pm}$. It is also possible to calculate the eigenvectors $\bm{X}$ and $\bm{Y}$ using the equations

[TABLE]

where

[TABLE]

In our experience for the systems considered in this work, taking the transition density matrix to be either $\bm{Z}^{+}$ or $\bm{Z}^{-}$ provides reasonable results but the non-uniqueness is not satisfactory. There has been some work on uniquely determining the transition density matrix including ph-hp correlation Etienne (2015); Luzanov and Zhikol (2010); Yasuike et al. (2011). Given that the different methods return similar values we use the simplest method from Yasuike et al. (2011) where two pseudo transition matrices are defined and then we use the following equation for the collectivity

[TABLE]

Appendix C Calculation of the field enhancement

To calculate the field enhancement using the results of a Casida simulation we begin by considering the real part of the density matrix Casida (1995)

[TABLE]

where $v_{jb\tau}=\int d^{3}r\phi_{j\tau}(\bm{r})v^{appl}(\bm{r},\omega)\phi_{b\tau}(\bm{r})$ and $v^{appl}$ is the perturbation. This equation can be rewritten using, a spectral expansion, in terms of the Casida eigenvalues and eigenvectors

[TABLE]

The induced density can be defined as Casida (1995)

[TABLE]

and, if we consider a perturbation of the form $v^{appl}(\bm{r},\omega)=\bm{r}\cdot\bm{E}$, then

[TABLE]

In its present form the pole structure of the equations means the induced density diverges at the excitation energies, thus we add a finite life time correction $\eta$

[TABLE]

which agrees with the result of Rossi Rossi (2013). Usually one is most interested in frequencies at the Casida excitation energies, therefore we can save some computational effort by expanding about the point $\omega=\omega_{I}$ and reordering the summations to get our final results

[TABLE]

Where we have used the following definitions for the transition density and transition dipole moment for the transition between the many body states $\Psi_{0}$ to $\Psi_{I}$ in the Casida formalism Casida (1995); Yasuike et al. (2011)

[TABLE]

and

[TABLE]

where we assume real KS orbitals and $\bm{\mu}_{ia}=-\int\ \phi_{i}(\bm{r})\phi_{a}(\bm{r})\ \bm{r}\ d^{3}r$. Equation 23 is an intuitive result and has been derived via a different method by Cocchi et al Cocchi et al. (2012). The induced density for a particular excitation has a spatial profile given by the transition density and its magnitude is given by the projection of the electric field on the transition dipole moment, which describes how strongly the excitation couples to the perturbation. It is also inversely proportional to the linewidth of the excitation. The induced density can then be used to calculate the induced electric field and the consequential field enhancement.

Appendix D Na60 chain

In figure 6 results for a Sodium Na60 chain with the spectrum calculated up to $2\ eV$ in figure 6. The chain has $60$ occupied states and a total of 80 unoccupied states were considered. The log of the dipole strength function is taken to clearly show the higher order longitudinal plasmon modes. As the chain size is increased these modes bunch together and in the infinite limit give the dispersion curve of a plasmon in a 1D electron gas Bernadotte et al. (2013).

Appendix E Fitting to the classical model

Here we fit our results to a classical ellipsoid model within the electrostatic approximation Bohren and Huffman (2008). The conditions for a localised resonance of an ellipsoid, with semiaxes $a,b,c$ is given by poles of the polarizability

[TABLE]

where $i=\{a,b,c\}$, $\epsilon_{m}$ is the permittivity of the surrounding medium (we take it to be one). We are interested in the particular case of a prolate spheroid where $b=c$ and $a$ is the largest semiaxes, in this case $L_{a}$ can be found analytically

[TABLE]

where $e$ is the eccentricity of the ellipsoid. Following Yan and Gao Yan and Gao (2008) we define $a=\Delta+(N-1)\frac{d}{2}$ and $b=c=\Delta$ where $d=3.08\mbox{\AA}$ is the bonding length used in the simulation and $\Delta$ is a measure of the electron spill out at the edges of the chain and is used as a fitting parameter. Yan and Gao Yan and Gao (2008) used a fixed bulk value for the plasma frequency in the Drude dielectric function used for $\epsilon(\omega)$, effectively adding another fitting parameter. Instead we define a chain size dependent plasma frequency

[TABLE]

meaning we need only fit $\Delta$. For large $N$ the plasma frequency becomes independent of $N$ as one would expect. Once we have fitted for $\Delta$ we can then calculate $L_{b}$ via the relation

[TABLE]

Applying the above procedure, we obtain the results show in the main text where $\Delta=1.45d$. The polarizability can be calculated using

[TABLE]

from which the (averaged over the 3 spatial dimensions) absorption cross section (which dominates over scattering for small particles) can be calculated

[TABLE]

which was used to obtain the results in figure 2c in the main text.

If we instead we take the plasma frequency to be fixed at the bulk value (we use $\omega_{p}=3.83\ eV$ as used by Yan and Gao Yan and Gao (2008)) then we obtain a best fit of $\Delta=1.13d$. The absorption cross section and the results of the fit for different chain sizes are shown in figure 7a and b respectively.

Another potential method is to consider only the ratio of the T and L mode energies $\frac{\omega_{T}}{\omega_{L}}$ which is convenient as the plasma frequency does not need to be considered. One can find that

[TABLE]

which can be compared to the linear response DFT results. The fitting parameter is obtained for the largest chain and comes out to $\Delta=0.79d$ and the results are shown in see figure 7c. In figure 7d we show a cross section view of the ground state density, along the z direction, calculated by DFT and the radii for the various fits discussed are shown in comparison.

Appendix F Transverse gap plasmons

In figure 8 we show the electronic spectra for a \ceNa20 chain with a gap opened up in the center. We observe some interference effects between the TC and TE modes at the gap sizes considered and consequently some new modes are setup in the energy gap between the TC and TE modes (which are approximately at the same energy for the \ceNa10 and \ceNa20 chain).

Appendix G Induced potential of gap plasmons

The induced potential can be calculated from the induced density by solving Poisson’s equation $\nabla^{2}(\delta\phi(\bm{r},\omega))=e\frac{\delta n(\bm{r},\omega)}{\epsilon_{0}}$, the results of this calculation for the 6 major low-energy tunnelling modes is shown in figure 9. As well as allowing the induced electric field to be calculated, it is also useful for identification of excitations. We find that the CTP1 mode is most CTP-like and the modes become increasingly BDP-like with increasing energy.