# Optical response of silver clusters and their hollow shells from linear-response TDDFT

Peter Koval, Federico Marchesin, Dietrich Foerster, Daniel Sánchez-Portal

► To cite this version:

Peter Koval, Federico Marchesin, Dietrich Foerster, Daniel Sánchez-Portal. Optical response of silver clusters and their hollow shells from linear-response TDDFT. Journal of Physics: Condensed Matter, 2016, 28 (21), pp.214001 (1-21). 10.1088/0953-8984/28/21/214001. hal-01342926

HAL Id: hal-01342926

https://hal.science/hal-01342926

Submitted on 7 Jul 2016

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

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

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

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
Optical response of silver clusters and their
hollow shells from linear-response TDDFT

Peter Koval<fnref n="1" />, Federico Marchesin<fnref n="1" />, Dietrich Foerster<fnref n="2" />
and Daniel Sánchez-Portal<fnref n="1" />

E-mail: koval.peter@gmail.com and sqbsapod@ehu.es

Abstract

We present a study of the optical response of compact and hollow icosahedral clusters
containing up to 868 silver atoms by means of time-dependent density functional theory. We
have studied the dependence on size and morphology of both the sharp plasmonic resonance
at 3–4 eV (originated mainly from sp-electrons), and the less studied broader feature appearing
in the 6–7 eV range (interband transitions). An analysis of the effect of structural relaxations,
as well as the choice of exchange correlation functional (local density versus generalised
gradient approximations) both in the ground state and optical response calculations is also
presented. We have further analysed the role of the different atom layers (surface versus inner
layers) and the different orbital symmetries on the absorption cross-section for energies up to
8 eV. We have also studied the dependence on the number of atom layers in hollow structures.
Shells formed by a single layer of atoms show a pronounced red shift of the main plasmon
resonances that, however, rapidly converge to those of the compact structures as the number of
layers is increased. The methods used to obtain these results are also carefully discussed. Our
methodology is based on the use of localised basis (atomic orbitals, and atom-centered and
dominant-product functions), which bring several computational advantages related to their
relatively small size and the sparsity of the resulting matrices. Furthermore, the use of basis
sets of atomic orbitals also allows the possibility of extending some of the standard population
analysis tools (e.g. Mulliken population analysis) to the realm of optical excitations. Some
examples of these analyses are described in the present work.

Keywords: TDDFT, atomic orbitals, product basis, silver clusters, silver shells, GGA kernel, response function

1. Introduction

Scientific interest in metallic clusters persists after many
decades of extensive study [1]. Among the other metals, silver
is one of the most popular materials for the production of clus-
ters. The prime interest in silver clusters originates from the
reduced chemical activity of silver [2], the best electron con-
ducting properties, a sharp plasmonic resonance at a relatively
low frequency and an affordable price of the metal. Silver
clusters are widely used in various applications [3–6], and in

the past were also widely theoretically studied [7–11]. It is
experimentally established that silver nano-particles (NPs)
with diameters 10 to 80nm exhibit a surface-plasmon reso-
nance in a wide frequency range from blue to green visible
light depending on their size. Moreover, by controlling the
shape of silver NPs one can further lower the resonance fre-
quency [12–14]. Recently, subnanometric noble metal clusters
have also been produced and have attracted much attention
in the research community [15]. The use of subnanometric
silver clusters has been demonstrated in catalysis [16, 17], in

<footnote n="1">¹ Centro de Física de Materiales, Centro Mixto CSIC-UPV/EHU and Donostia International Physics Center (DIPC), Paseo Manuel de Lardizabal 5, 20018 Donostia-San Sebastián, Spain</footnote>

<footnote n="2">² Laboratoire Ondes et Matière d'Aquitaine, Université de Bordeaux 1, 351 Cours de la Liberation, 33405 Talence, France</footnote>
chemical sensing [18–20] and in surface-enhanced Raman scattering [21, 22], among other fields.

From a theoretical point of view, the plasmonic properties of large silver particles (20 nm and larger) can be satisfactorily described by classical Mie theory [7]. Further extensions of the classical description include the discrete-dipole approximation [8] and the finite-difference time-domain techniques [23] that are capable of describing the response of classical objects of any shape. For smaller clusters with an effective diameter less than 3 nm (~600 atoms), the classical theory cannot give a rigorous description [8, 24–26] because the atomistic details can significantly alter the classically averaged picture and it is necessary to give an accurate description of the electronic distribution and scattering at the surfaces. Therefore, quantum mechanical methods had been widely applied in studies of silver clusters [9–11, 27], silver shells [28, 29], and alloys [30, 31]. In this paper, we will focus on relatively large icosahedral clusters containing up to 561 atoms and shells of up to 868 atoms that were not addressed in the past.

From a methodological perspective, there are early studies of silver clusters of different sizes and shapes using Hückel models [32], tight-binding techniques [25, 33], density-functional theory (DFT) [27, 34–37], time-dependent DFT (TDDFT) [11, 36, 38, 39] and many-body perturbation theory [40, 41]. In this work, we apply linear-response TDDFT using a linear combination of atomic orbitals (LCAO) method to describe the electronic states of the clusters. Our method [42, 43] has been recently enhanced in several respects. Besides the generality of the geometries and chemical species characteristic of ab initio methods (in the present case our linear-response solver is coupled to the SIESTA method [44, 45]), the main advantage of the method is its computational efficiency that stems from the use of an iterative scheme to compute the optical response, as well as an efficient basis to express the products of atomic orbitals. Importantly, the number of iterations does not scale significantly with the size of the system studied. The frequency-by-frequency operation is a useful feature of the method, for instance in the calculation of the intensity of non-resonant Raman scattering, and in other situations when the target frequency range is small but arbitrarily-placed.

In this paper, after a description of the iterative TDDFT method in section 2, we apply the method to find the optical response of silver clusters and shells in section 3. We show that the surface-plasmon resonance of atomically-thin shells is substantially different from the response of compact clusters, although the difference quickly becomes marginal as the thickness of the shell is increased. A comparison with other theories and experiments is made in section 4, and section 5 summarises our results.

## 2. Methods

We focus on an ab initio atomistic description of small silver NPs with a diameter less than 3 nm. The NPs of that size typically contain several hundreds of atoms, which represents a

major computational difficulty for the quantum mechanical description of such systems with most current methodologies. To the best of our knowledge, only TDDFT can currently cope with such large systems practically within an ab initio framework. Moreover, due to the size of these systems, we need a method of low computational complexity for which one can envision several candidates. One widely known method is the wave-packet propagation [11, 36, 46–48]. If properly implemented the wave-packet propagation can be realised in O(N) operations where N is the number of atoms. A second type of method that is up to the task would be the recently developed stochastic methods [49, 50], for which a lower-than-linear computational complexity scaling has been claimed. A Sternheimer approach to the linear-response TDDFT also seems a viable alternative [51, 52].

The method that we are utilising in this work is an efficient iterative way of solving linear response equations in which we exploit the locality of the operators and use an LCAO expansion of the Kohn–Sham (KS) eigenstates [42, 43]. Although the method has a relatively high asymptotic scaling of the computational complexity O(N³), it has been demonstrated to be a useful alternative to other methods [26, 53]. Moreover, because of the description offered by LCAOs, we can adapt population analysis tools that allow the connection of the electronic structure of the system to the chemical intuition (e.g. the Mulliken population [54]) to the realm of the optical response, as we will describe below.

### 2.1. Response functions formalism

The basic quantity of the linear response TDDFT is the density response function $\chi(r, r', \omega)$, which is a kernel of an integral operator delivering the density change $\delta n(r, \omega)$ in response to a small external perturbation $\delta V_{\text{ext}}(r, \omega)$

$$ \delta n(r, \omega) = \int \chi(r, r', \omega) \delta V_{\text{ext}}(r', \omega) dr'. \quad (1) $$

In general, the density response function can be expressed in terms of eigenstates and eigenvalues of the Schrödinger equation, starting the derivation from a perturbative ansatz for the density change $\delta n(r, t) = \Psi(r, t)\delta\Psi(r, t) + \text{c.c.}$ In the case of the KS formulation of DFT, the density response function becomes more complicated due to the dependence of the KS effective potential $V_{\text{eff}}(r) = V_{\text{ext}}(r) + V_{\text{H}}(n, r) + V_{\text{xc}}(n, r)$ on the electronic density $n(r)$ via Hartree and exchange-correlation (xc) potentials, respectively. However, the so-called non-interacting response function $\chi_0(r, r', \omega) \equiv \frac{\delta n(r, \omega)}{\delta V_{\text{eff}}(r', \omega)}$ remains completely analogue to the Schrödinger’s response function. Although the explicit expression of the non-interacting response function is widely known in the literature [55, 56], we will repeat it here for the sake of completeness

$$ \chi_0(r, r', \omega) = \sum_{nm} (f_n - f_m) \frac{\Psi_n^*(r) \Psi_m(r) \Psi_m^*(r') \Psi_n(r')}{\omega - (E_m - E_n) + i\epsilon}. \quad (2) $$

Here, the occupations $f_n$ and energies $E_n$ of KS eigenstates $\Psi_n(r)$ do appear and $\epsilon$ is an infinitesimal number that accounts for the proper causality of the response. A finite value can be
given to $\epsilon$, in which case it can be thought to represent the finite lifetime of excitations. The interacting response function, defined by equation (1) can be related to the non-interacting response function [57, 58] via the so-called interaction

$$\text{kernel } K(\mathbf{r}, \mathbf{r}') \equiv \frac{\delta V_{\text{eff}}(\mathbf{r})}{\delta n(\mathbf{r}')}$$

$$\chi(\mathbf{r}, \mathbf{r}', \omega) = \chi_0(\mathbf{r}, \mathbf{r}', \omega) + \int \chi_0(\mathbf{r}, \mathbf{r}'', \omega) K(\mathbf{r}'', \mathbf{r}'') \chi(\mathbf{r}'', \mathbf{r}', \omega) d\mathbf{r}'' d\mathbf{r}'''. \quad (3)$$

The TDDFT interaction kernel $K(\mathbf{r}, \mathbf{r}')$ is commonly separated into the Hartree and xc kernels [57, 58]

$$K(\mathbf{r}, \mathbf{r}') = \frac{1}{|\mathbf{r} - \mathbf{r}'|} + \frac{\delta V_{\text{xc}}(\mathbf{r})}{\delta n(\mathbf{r}')}. \quad (4)$$

In this work, we mainly use the generalised gradients approximation (GGA) kernel (see appendix B). The effect of the TDDFT kernel and the DFT potential on the optical response is evaluated below in section 2.5 by comparing the GGA results with those of the local density approximation (LDA).

Because we are interested in optical perturbations, the wavelength of which exceeds 150 nm, which is much larger than the characteristic length of our systems, the coupling to the external electromagnetic stimuli can be correctly taken into account via a simple dipole operator $\delta V_{\text{ext}}(\mathbf{r}, \omega) \propto \mathbf{r}$. Furthermore, the far-field response is connected to the polarisability tensor of the quantum system

$$P_{ij}(\omega) = \int \mathbf{r}_i \chi(\mathbf{r}, \mathbf{r}', \omega) \mathbf{r}_j' d\mathbf{r} d\mathbf{r}' \quad (5)$$

which gives rise to the orientation-averaged optical cross section [59]

$$\sigma_{\text{opt}}(\omega) = \frac{4\pi\omega}{3c} \sum_{i} P_{ii}(\omega). \quad (6)$$

In this paper we are concerned with the calculation of the optical cross section (6) using LDA and GGA DFT functionals and a basis set of local orbitals to expand the KS orbitals $\Psi_m(\mathbf{r})$.

## 2.2. Product basis set

The eigenstates entering the response function (2) are sought within LCAO

$$\Psi_n(\mathbf{r}) = \sum_a X_a^n f^a(\mathbf{r}), \quad (7)$$

where the expansion coefficients $X_a^n$ are determined in a diagonalisation procedure. The atomic orbitals $f^a(\mathbf{r})$ are given by a product of radial functions and spherical harmonics. The LCAO solution is set up and solved within the DFT package SIESTA [44, 45].

The products of eigenstates $\Psi_n^*(\mathbf{r})\Psi_n(\mathbf{r})$ in equation (2) give rise to products of atomic orbitals $f^{a*}(\mathbf{r})f^b(\mathbf{r})$. Furthermore, we aim at solving the integral equation (3) for the interacting response function. In order to turn this equation into an algebraic equation that is easily solved, one has to use a basis set of functions capable of spanning the space of atomic-orbital

products. This set of basis functions should be as small as possible and contain preferably localised functions. There are several options to construct such a set of functions, hereafter product basis. The most widely known is probably the auxiliary functions for Gaussian basis sets [60, 61]. However, the SIESTA method is based on so-called numerical atomic orbitals (NAOs), which can be more flexible and economic than Gaussian basis sets. There are methods to construct the product basis sets for numerical orbitals [62, 63] and also our method of so-called dominant products [64, 65]. The dominant product basis can be very accurate but requires a large number of functions that becomes prohibitive for large systems. Thus, in this work we project our dominant products onto a basis set of atom-centered functions in order to reduce the basis set size. As we will see below, the use of this more economical basis increases the range of applicability and efficiency of the iterative scheme without a significant loss of accuracy<fnref n="3" />.

In the method of dominant products, we utilise a simple ansatz for the products of atomic orbitals

$$f^a(\mathbf{r})f^b(\mathbf{r}) = V_{\mu}^{ab}F^{\mu}(\mathbf{r}). \quad (8)$$

In this equation, the complex conjugation does not appear because we use real-valued atomic orbitals [66]. The product 'vertex' coefficients $V_{\mu}^{ab}$ and the dominant products $F^{\mu}(\mathbf{r})$ are determined in a diagonalisation-based procedure. Namely, we aim at identifying linear combinations of the original atomic-orbital products $f^a(\mathbf{r})f^b(\mathbf{r})$

$$F^{\mu}(\mathbf{r}) = \Lambda_{ab}^{\mu} f^{a}(\mathbf{r}) f^{b}(\mathbf{r}). \quad (9)$$

that are orthogonal to each other with respect to a Coulomb metric

$$g^{ab,a'b'} = \int f^a(\mathbf{r}) f^b(\mathbf{r}) \frac{1}{|\mathbf{r} - \mathbf{r}'|} f^{a'}(\mathbf{r}') f^{b'}(\mathbf{r}') d\mathbf{r} d\mathbf{r}'. \quad (10)$$

The linear combinations are built with eigenvectors of the metric $g^{ab,a'b'}$

$$g^{ab,a'b'} \Lambda_{a'b'}^{\mu} = \lambda^{\mu} \Lambda_{ab}^{\mu} \quad (11)$$

that guarantee the orthogonality requirement. Moreover, the eigenvalues of Coulomb metric $\lambda^{\mu}$ are used as an indicator of the importance of a particular linear combination $\mu$ to the completeness of the basis $\{F^{\mu}(\mathbf{r})\}$. Namely, the norm of dominant products $F^{\mu}(\mathbf{r})$ (9) is proportional to the eigenvalue $\lambda^{\mu}$. Therefore, we can consistently limit the number of dominant product (9) by ignoring the eigenvector such that $\lambda^{\mu}$ is lower than a certain eigenvalue threshold.

For our purposes here, it is only necessary to add that the procedure is applied to each atom-pair individually. This keeps the operation count at $O(N)$ scaling, generates localised dominant products $F^{\mu}(\mathbf{r})$ and also determines the sparsity

<footnote n="3">³ We performed many test calculations of the optical polarisability comparing the atom-centered functions and the dominant product basis sets. For double-zeta polarized and larger LCAO basis sets, the discrepancies are negligible and also decrease with system size. For smaller basis sets containing only s and p angular momentum symmetries, the description worsens, but can be easily recovered by adding the higher angular momentum orbitals while generating the atom-centered product basis.</footnote>
properties of the product vertex coefficients $V_{\mu}^{ab} = \Lambda_{ab}^{\mu}$ in equation (8). Namely, the vertex coefficients $V_{\mu}^{ab}$ form a double-sparse table, which needs asymptotically only $O(N)$ stored numbers. The term 'double-sparse table' means that only summation over two indices of this table generates a full object (vector), while the summation over one of the indices generates a sparse matrix. For instance, a summation over the product index $\mu$ generates a matrix $s^{ab} = \sum_{\mu} V_{\mu}^{ab}$, which has the sparsity of the usual overlap matrix $S^{ab} = \int f^{a}(\mathbf{r})f^{b}(\mathbf{r})d\mathbf{r}$, while the summation over one of the orbitals generates a rectangular sparse matrix $R_{\mu}^{a} = \sum_{b} V_{\mu}^{ab}$ because, by construction, a dominant product index $\mu$ is connected to orbital indices a and b of one atom (local pairs) or two atoms with overlapping orbitals (bilocal pair) rather than to all the orbital indices in the molecule.

The dominant products described above have been used in TDDFT, Hedin's GW approximation and for solving a Bethe-Salpeter equation [42, 43, 67, 68]. However, the construction of dominant products, although mathematically rigorous and sparsity-preserving, has the important disadvantage of generating a large number of functions. This disadvantage stems from the construction procedure, which is repeated independently for each atom pair. It is easy to see that the dominant products $F^{\mu}(\mathbf{r})$ can strongly overlap because different atom pairs can have the same or close centers at which the products have their maximal values. This fact results in a redundant description of the orbital products by the dominant product basis when looking from a perspective of the whole system. In order to correct for this, we use an ansatz for the auxiliary basis set that is widely known in quantum chemistry, [60, 61] and also in more 'physics-oriented' proposals [62, 63, 69]. Namely, the cited works assume that solely atom-centered product functions $A^{\mu}(\mathbf{r})$ are sufficient in practice to express all orbital products $f^{a}(\mathbf{r})f^{b}(\mathbf{r})$. This very statement, although only justified a posteriori, is a useful piece of advice that allows the reduction of the linear dependencies in a product basis set because the atom centers are separated from each other at least with a bonding distance which prevents strong overlaps of the resulting functions.

Here we take the local dominant products $F^{\mu}(\mathbf{r})$ (i.e. dominant products generated for orbitals in the same atom, as opposed to bilocal pairs) as the atom-centered functions $A^{\mu}(\mathbf{r})$. The ansatz for atomic orbitals, analogous to equation (8) can be immediately written as

$$f^{a}(\mathbf{r})f^{b}(\mathbf{r}) = P_{\mu}^{ab}A^{\mu}(\mathbf{r}), \qquad (12)$$

where the atom-centered product vertex coefficients $P_{\mu}^{ab}$ must be still determined. In this work, we extend the ansatz (12) with a recipe for choosing these vertex coefficients $P_{\mu}^{ab}$. We propose to draw a sphere around a given atom pair with a radius corresponding to the maximal spatial extension of its orbital products (here defined approximately as the maximum radius of their atomic orbitals), and consider the atom centers within that sphere as contributing to the atom pair. The ansatz (12) can be easily resolved to obtain $P_{\nu}^{ab} = T^{ab,\mu}(v^{\mu\nu})^{-1}$,

where $v^{\mu\nu}$ is the Coulomb matrix element between functions $A^{\mu}(\mathbf{r})$ and $A^{\nu}(\mathbf{r})$, and $T^{ab,\mu}$ is the corresponding matrix element between $A^{\mu}(\mathbf{r})$ and the product of orbitals $f^{a}(\mathbf{r})f^{b}(\mathbf{r})$. Notice that in order to termine $P_{\mu}^{ab}$ we do not need to invert the whole matrix of the Coulomb metric, but a smaller sub-matrix corresponding to atom-centered functions inside the contributing sphere. Thus, this step is not computationally prohibitive. However, the table $P_{\mu}^{ab}$ has a limited value for practical calculations. Namely, the table $P_{\mu}^{ab}$ can have an order of magnitude more non-zero elements that the product vertex coefficients $V_{\mu}^{ab}$. This dramatic difference arises because of the distant bilocal atom pairs for which we have very few dominant products $F^{\mu}(\mathbf{r})$ instead of many atom-centered functions $A^{\mu}(\mathbf{r})$ contributing to such pairs.

A more fruitful idea proves to be a re-expression of the bilocal dominant products $F^{\mu}(\mathbf{r})$ in terms of atom-centered products $A^{\nu}(\mathbf{r})$ (that are also chosen within the sphere of contributing centres)

$$F^{\mu}(\mathbf{r}) = c_{\nu}^{\mu}A^{\nu}(\mathbf{r}), \qquad (13)$$

where the projection coefficients $c_{\nu}^{\mu}$ can be also readily expressed as

$$c_{\nu}^{\mu} = M^{\mu\nu'}(v^{\mu\nu'})^{-1} \qquad (14)$$

in terms of the matrix elements of the Coulomb interaction

$$M^{\mu\nu} = \int \frac{F^{\mu}(\mathbf{r})A^{\nu}(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r} d\mathbf{r}', \quad v^{\mu\nu} = \int \frac{A^{\mu}(\mathbf{r})A^{\nu}(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r} d\mathbf{r}'. \qquad (15)$$

The projection ansatz (13) is useful because it is computationally very fast to turn a vector in the atom-centered basis to the dominant product basis and back, depending on the quantity that needs to be treated on the product basis. Namely, we found that it is faster to apply the non-interacting response on the basis of dominant products and, on the other hand, it is faster to compute and easier to store the TDDFT kernel in the basis of atom-centered functions.

## 2.3. Iterative method for computing the polarisability

By inserting the ansatzes (7) and (8) into the non-interacting response function (2) we obtain

$$\chi_{0}(\mathbf{r}, \mathbf{r}', \omega) = F^{\mu}(\mathbf{r})\chi_{\mu\nu}^{0}(\omega)F^{\nu}(\mathbf{r}'), \qquad (16)$$

where the tensor of non-interacting response is given by

$$\chi_{\mu\nu}^{0}(\omega) = \sum_{nm} (f_{n} - f_{m}) \frac{(X_{a}^{n} V_{\mu}^{ab} X_{b}^{m})(X_{c}^{m} V_{\nu}^{cd} X_{d}^{n})}{\omega - (E_{m} - E_{n}) + i\epsilon}. \qquad (17)$$

Furthermore, we assume for the interacting response function (3) an expression similar to equation (16). This turns the Petersilka-Gossmann-Gross equation (3) into a matrix equation

$$\chi = \chi^{0} + \chi^{0}K\chi, \qquad (18)$$

where the product indices $\mu\nu$ are dropped and the kernel matrix K is given by
$$K^{\mu\nu} = \int \frac{F^{\mu}(\mathbf{r})F^{\nu}(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r} d\mathbf{r}' + \int F^{\mu}(\mathbf{r})K_{\text{xc}}(\mathbf{r})F^{\nu}(\mathbf{r}')d\mathbf{r}'. \quad (19)$$

In this work, we use LDA and GGA kernels $K_{\text{xc}}(\mathbf{r}) = \frac{\delta^2 E_{\text{xc}}}{\delta n(\mathbf{r})\delta n(\mathbf{r}')}$, which are computed in the atom-centered functions $A^{\mu}(\mathbf{r})$. Explicit expressions for the GGA kernel are discussed in the appendix B.

The whole response matrix $\chi_{\mu\nu}(\omega)$ is superfluous in the computation of the electronic polarisability (5). Therefore, we further introduce the product basis set into the equation for polarisability (5) and using equation (18) obtain

$$P_{ij}(\omega) = d_i^{\mu}(\delta_{\nu'}^{\mu} - \chi_{\mu\nu'}^{0}(\omega)K_{\nu'\mu'}^{0})^{-1}\chi_{\mu'\nu'}^{0}(\omega)d_j^{\nu}, \quad (20)$$

where the dipole moments of the product functions $d_i^\mu = \int F^\mu(\mathbf{r})r_i d\mathbf{r}$ appear. The calculation of the polarisability (20) is split into a calculation of the density change in the product basis $\delta n_\mu(\omega)$

$$[\delta - \chi^0(\omega)K]\delta n(\omega) = \chi^0(\omega)d \quad (21)$$

and a final trace with the dipole moments $P_{ij}(\omega) = d_i^\mu \delta n_{\mu,j}(\omega)$. The linear equation (21) is solved with an iterative method similar to the Arnoldi method and optimised for delivering the polarisability $P_{ij}(\omega)$ with a given precision rather than the density change $\delta n_{\mu,j}(\omega)$ [42] that would be the target for general-purpose iterative solvers [70].

Iterative linear equation solvers require only the action of a matrix onto given vectors to find the solution $\delta n_{\mu,j}(\omega)$. In our case, the matrix reads $A = [\delta - \chi^0(\omega)K]$. The product of this matrix with a vector $z$ can be computed in terms of subsequent matrix-vector products of the TDDFT kernel $K$ with a vector and of the non-interacting response function $\chi^0$ with another vector. The former product is easy to organise because we represent the kernel as a full matrix between the atom-centered product functions. The latter product is more involved and is explained below.

The matrix-vector product of the non-interacting response function $\chi_{\mu\nu}^0(\omega)$ with a vector $z^\nu$ is also split into a sequence of matrix-vector and matrix-matrix operations. First, we use the projection ansatz (13), insert that expression into the equation (16) and get

$$\chi_0(\mathbf{r}, \mathbf{r}', \omega) = A^{\mu}(\mathbf{r})c_{\mu}^{\bar{\mu}}\chi_{\bar{\mu}\nu}^{0}(\omega)c_{\nu}^{\bar{\nu}}A^{\nu}(\mathbf{r}'), \quad (22)$$

where the tilde indices run over the dominant products and the simple indices run over the atom-centered product functions. Equations (17) and (22) define the matrix expression for the matrix-vector product $\chi_{\mu\nu}^0(\omega)z^\nu$. There are several sequences of operations to organise the matrix-vector product $\chi_{\mu\nu}^0(\omega)z^\nu$. However, our original algorithm suggested in [42] proved to be the fastest alternative, with moderate memory requirements. The algorithm starts with precomputing of a quantity $\alpha_{\bar{\mu}}^{\bar{\nu}} = V_{\bar{\mu}}^{ab}X_b^n$ for occupied eigenstates states $n$. The table $\alpha_{\bar{\mu}}^{\bar{\nu}}$ is stored in a block-sparse storage that uses $O(N^2)$ elements of the random access memory (RAM) and enables matrix operations with ordinary basic linear algebra

subroutines [71]. The table $\alpha_{\bar{\mu}}^{\bar{\nu}}$ is the major ‘memory consumer’ that, however, can be obviously eliminated when larger systems need to be treated. A vector $z^\mu$, on which the response function has to be applied, is converted to the basis of dominant products $z^\bar{\mu} = c_\mu^{\bar{\mu}}z^\mu$, then the product indices are summed to produce a matrix $\beta^{\bar{\nu}} = \alpha_{\bar{\mu}}^{\bar{\nu}}z^\bar{\mu}$. The matrix $\beta^{\bar{\nu}}$ is a full rectangular matrix. The calculation continues with a matrix-matrix multiplication $\gamma^{\bar{\mu}\bar{\nu}} = X_a^m\beta^{\bar{\mu}\bar{\nu}}$, where the index $m$ now runs over the unoccupied KS orbitals. The latter multiplication determines a maximal computational complexity of the whole algorithm, which is $O(N^3)$. The calculation continues with an update of the matrix $\gamma^{\bar{\mu}\bar{\nu}}$ with the frequency-occupation mask $\tilde{\gamma}^{\bar{\mu}\bar{\nu}} = \gamma^{\bar{\mu}\bar{\nu}}((f_n - f_m)/(\omega - (E_m - E_n) + i\epsilon) - (f_n - f_m)/(\omega + (E_m - E_n) + i\epsilon))$. The next $O(N^3)$ matrix-matrix multiplication reads $\tilde{\beta}^{\bar{\mu}\bar{\nu}} = X_a^m\tilde{\gamma}^{\bar{\mu}\bar{\nu}}$. Finally the non-interacting density change $\delta n_{\bar{\mu}}^0 = \chi_{\bar{\mu}\bar{\nu}}^0 z^\bar{\nu}$ is obtained by tracing over $a$ and $n$ indices in the product of $\tilde{\beta}^{\bar{\mu}\bar{\nu}}$ with the precomputed quantity $\alpha_{\bar{\mu}}^{\bar{\nu}}$, i.e. $\delta n_{\bar{\mu}}^0 = \tilde{\beta}^{\bar{\mu}\bar{\nu}}\alpha_{\bar{\mu}}^{\bar{\nu}}$.

In summary, we described the iterative algorithm of $O(N^3)$ computational complexity that uses $O(N^2)$ memory and enables a relatively fast calculation of interacting polarisability in plasmonic systems, i.e. in systems that have many nearly-degenerate transitions. Although the presented algorithm possesses a relatively high asymptotic computational complexity, the algorithm is relatively inexpensive in terms of computational resources. This allowed us to perform calculations for system sizes containing hundreds of atoms, despite the fact that our implementation uses only OpenMP parallelisation. For example, the calculation of the Ag561 icosahedral cluster was done on a 12-core node (Intel Xeon CPU X5550 2.67 GHz, release date 2009) in 25 h of walltime. The iterative procedure took most of the walltime (19.7 h) while it performed for 200 frequencies and only for the xx-component of polarisability tensor (20).

## 2.4. Accuracy of the methods

In the algorithm presented above, several approximations are involved. The approximations originating from the input DFT calculation, including the choice of the xc functional and the applied basis set of atomic orbitals, are discussed below in sections 2.5 and 2.6. The approximations originating from the implementation of TDDFT are related to the usage of the product basis sets and an iterative procedure to compute the induced density change for a given perturbation. Both of the latter approximations were carefully analysed in our previous works, in which we identified the corresponding accuracy indicators. The accuracy indicators include the difference between the overlap and dipole-matrix elements computed directly from the atomic orbitals and via the moments of the product basis functions, and a convergence test of the iterative procedure. The accuracy indicators were routinely controlled in the calculations we present in this work. As a result of this control, we found that the spectra presented in this work are unaffected by the usage of product basis sets until a high frequency approximately $\omega = 50$ eV for all cluster sizes;
the iterative procedure provides the same results for polarisability, within a given, previously specified small tolerance, as computed via the Casida formulation (possible only for small clusters with fewer than about 20 atoms in our realisation).

## 2.5. Choice of the exchange-correlation functional

All correlation effects should be captured in DFT through a single xc functional. Previous studies revealed that both the geometry [72] and the electronic structure [11, 73, 74] of silver-containing compounds are affected by the choice of xc functional. In the extensive comparative study in [72], it was shown that the lattice parameter predicted by the LDA is 1.6% shorter than the experimental value. This situation is improved by some GGA functionals. In particular, the functionals by Wu and Cohen (WC) [75], the Perdew-Burke-Ernzerhof adapted to solids [76], and that due to Armiento and Mattsson [77] provided the best results (deviations of the lattice parameter below 0.5%). Because the WC functional shows the best performance for bulk silver and also performs well for finite systems [78], we have chosen this functional for geometry optimisation.

With respect to electronic structure, it is well known from the literature that LDA and GGA predict a too low onset of the *d*-bands in the solid, and this reduces the intensity of the low-frequency plasmonic resonance produced mainly by the *s*-electrons [79]. This deficiency of LDA and GGA can be corrected by using so-called long-range-corrected (LRC) xc functionals, which contain a portion of Fock exchange [37, 41, 74], the van-Leeuwen-Baerends explicit ansatz with the correct asymptotic behaviour [80], or orbital-dependent functionals [11, 81]. Unfortunately, these functionals are not available within the publicly-available version of the SIESTA package [44]. Moreover, the LRC functionals referenced above give rise to a non-local two-point TDDFT kernel, which is either not known or computationally too expensive to treat the large systems addressed here. From this point of view, the Sternheimer approach and wave-packet propagation approach have an advantage over our iterative TDDFT: the former approaches only need the xc potential, and do not involve the TDDFT kernel. Fortunately, local and semi-local functionals correctly capture trends of the plasmonic response in nanoparticles, as was well documented in the past [73], and are still widely used [82].

For these reasons, we have chosen to use the WC functional for both the DFT and TDDFT steps of our calculation. In order to assess the effect of the xc functional, we compare the LDA and GGA spectra and also analyse the contributions from *d*-electrons to the total absorption cross section (in sections 2.7 and 3.5). In figure 1 we compare the optical absorption cross sections computed with the LDA and GGA functionals for the Ag<sub>147</sub> icosahedral cluster. The geometries of the cluster were relaxed with Perdew-Zunger (PZ) LDA or WC functionals, although the relaxations themselves did not affect the spectra significantly. We can immediately confirm that GGA increases the intensity of the plasmon peak as compared with LDA. However, it does not significantly shift the

resonant frequencies (3.5 eV LDA, 3.54 GGA). Moreover, it is interesting to note that the effect of the gradient corrections in the kernel $K_{xc}$ is marginal. Namely, in figure 1 we can hardly distinguish the spectra in which the GGA kernel is substituted by the LDA kernel, and the other computational parameters are kept the same (WC+WC versus WC+PZ curves). As expected, however, the effect of the Hartree kernel is crucial. The non-interacting response (not shown in the plots) is dominated at low energies by a peak at approximately 1 eV, much too low as compared to experiment. Therefore, we see a minor influence of xc kernel on the optical response of our Ag nanoparticles. It is important to note, however, that the improved approximation to the quasi-particle spectrum provided by the KS eigenvalues computed with the GGA functional is reflected in an improved description of the optical properties. This is important information from a methodological point of view because the calculation of real-space integrals of the xc kernel is much more time consuming (two orders of magnitude) than the calculation of the Hartree kernel. The Hartree kernel is calculated with the help of fast Bessel transforms [83, 84] and the multipole expansions [42, 85]. Moreover, the GGA xc kernel is much more cumbersome than the LDA kernel (see appendix B) and for many of the most sophisticated functionals the explicit expressions of the kernel are still lacking. The analysis presented above with respect to the influence of the xc kernel (RPA versus GGA) is consistent with that previously published for smaller clusters [39].

## 2.6. Choice of atomic-orbital basis set

The choice of the atomic-orbital basis set determines the quality of the LCAO calculations to a large extent. The numerical orbitals used in the SIESTA package are capable of approaching the results of plane-wave calculations for bulk systems [86, 87], at least for ground state properties and the description of the low lying unoccupied states. At the same time, it was found that semi-infinite systems (surface properties) need special care when being described using confined NAOs [88]. Namely, it was found that adding a single layer of floating orbitals can considerably improve the surface properties with only a small impact on the computational performance. Because the surface-to-volume ratio of our target clusters is rather large [1], we assess the quality of the default SIESTA basis set by augmenting it with an extra layer of floating orbitals.

In figure 2, we plot the absorption spectra of the Ag<sub>147</sub> cluster computed with a double-zeta polarised (DZP) basis (Energy Cutoff = 50 meV) and with the same DZP basis augmented with an extra layer of *s*- and *p*-floating orbitals. The floating orbitals are placed at the positions of the next (fifth) Mackay layer (the geometry of the clusters is discussed below in section 3) and are excluded from the relaxation procedure. The WC functional was used in both calculations. One can see in the figure that the extra layer of floating orbitals slightly red-shifts the frequency of the main plasmon resonance (by 0.02 eV) and increases the absorption cross section around the resonance. A direct analysis of the absorption cross section in terms of the different cluster layers (see section 2.7)
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Absorption cross section of the icosahedral Ag<sub>147</sub> cluster computed with LDA and GGA functionals. The curves are labelled according to the xc potential (first label) and xc kernel (second label) used. The Perdew-Zunger (PZ) LDA functional and Wu-Cohen (WC) GGA functionals are compared. The blue solid line represent the RPA K<sub>xc</sub> = 0 starting from the WC results. Relaxation of geometries was done with PZ or WC functionals.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Absorption cross section of the icosahedral Ag<sub>147</sub> cluster computed with a DZP basis set and with a DZP basis augmented with an extra layer of floating orbitals (FLTs). Relaxations have been done for all atoms except ‘ghost atoms’ that carry floating orbitals.</figcaption></figure>

shows that the enhancement is due to the extra layer of ‘ghost atoms’ which carry the floating orbitals (see sections 2.7 and 3.5), and that the main contribution to the cross section is still due to the outer layer of real atoms. At the same time, the layer of floating orbitals does not affect the computation cost dramatically because it does not unnecessarily decrease the sparsity of the matrices involved in the calculation, which is the case when spatially extended diffuse orbitals are added to the basis set.

## 2.7. Analysis of the interacting density change.

Here we focus on the analysis of the interacting polarisability in terms of the spatial distribution of atoms and in terms of the non-interacting electron-hole pairs. Other types of analysis

are too technical for a general presentation and we reserved them for appendix A.

2.7.1. Contribution of different atoms to the polarisability. In our framework, the interacting polarisability α(ω) is given by

$$ \alpha(\omega) = d^{\mu}\delta n_{\mu}(\omega), \quad d^{\mu} = \int A^{\mu}(\mathbf{r})rdr, \qquad (23) $$

where we drop all the Cartesian indices for the sake of clarity. The index μ runs over all product functions A<sup>μ</sup>(r), which are centered on real and ghost atoms. The sum over product functions can be split into sub-sums according to a given criterion. For instance, we can split the sum over the product indices into sub-sums over atomic layers in the cluster L: α(ω) = Σ<sub>L</sub> α<sub>L</sub>(ω) with a partial polarisability α<sub>L</sub>(ω) given by
$$ \alpha_L(\omega) = d^\mu P_\mu^\nu(L) \delta n_\nu(\omega). \qquad (24) $$

Here a padding operator $P_\mu^\nu(L)$ is introduced. The padding operator is a diagonal matrix whose elements are equal to one $P_\mu^\mu(L) = 1$ if the product index $\mu$ belongs to an atom in the layer L and zero otherwise.

### 2.72. Electron-hole expansion of the interacting induced density.
The analysis presented above for the interacting polarisability can be completed using an expansion of the interacting density change $\delta n(r, \omega)$ in terms of electron-hole pairs

$$ \delta n(\mathbf{r}, \omega) = \sum_{ij} \delta n_{ij}(\omega) \Psi_i^*(\mathbf{r}) \Psi_j(\mathbf{r}), \qquad (25) $$

where $\Psi_j(\mathbf{r})$ are the KS eigenstates.

This kind of expansion naturally arises in the Casida formulation of TDDFT [89], and this is perhaps at the root of the popularity of Casida's formulation, because the electron-hole expansion (25) allows the classification of a given excitation in terms of the character of the non-interacting transitions contributing to it, as was recently elaborated with an alternative method [90].

Obtaining the expansion (25) is relatively straightforward in the iterative formulation of the TDDFT presented above. The interacting density change $\delta n(r, \omega)$ is given by (compare with equation (1))

$$ \delta n(\mathbf{r}, \omega) = \int \chi_0(\mathbf{r}, \mathbf{r}', \omega) \delta V_{\text{eff}}(\mathbf{r}', \omega) d\mathbf{r}', \qquad (26) $$

where $\delta V_{\text{eff}}(\mathbf{r}, \omega)$ is an effective (screened) perturbation. Using sum-over-states representation (2) of the non-interacting density response function we can immediately see the electron-hole expansion coefficients $\delta n_{ij}(\omega)$

$$ \delta n_{ij}(\omega) = \frac{f_i - f_j}{\omega - (E_j - E_i) + i\epsilon} \int \Psi_j^*(\mathbf{r}') \Psi_i(\mathbf{r}') \delta V_{\text{eff}}(\mathbf{r}', \omega) d\mathbf{r}'. \qquad (27) $$

Now it remains only to derive an equation for the effective perturbation $\delta V_{\text{eff}}(\mathbf{r}, \omega)$. Using the Petersilka-Gossmann-Gross equation (3) and the two alternative expressions (26) and (1) for the interacting density change, it is straightforward to derive an equation for the effective perturbation

$$ [\delta(\mathbf{r} - \mathbf{r}') - K(\mathbf{r}, \mathbf{r}'')\chi_0(\mathbf{r}'', \mathbf{r}', \omega)] \delta V_{\text{eff}}(\mathbf{r}', \omega) = \delta V_{\text{ext}}(\mathbf{r}, \omega). \qquad (28) $$

In terms of product basis (8) and for optical absorption $\delta V_{\text{ext}}(\mathbf{r}, \omega) \equiv \mathbf{r}$, the latter equation transforms into a linear algebraic equation

$$ [\delta - K\chi_0(\omega)] \delta V_{\text{eff}}(\omega) = d, \qquad (29) $$

where we dropped the Cartesian indices and product indices for clarity. This equation can be solved iteratively with a generalised minimal residue solver [70]. Using the product basis again in equation (27) we get

$$ \delta n_{ij}(\omega) = \frac{(f_i - f_j)(X_a^j V_\mu^{ab} X_b^i)}{\omega - (E_j - E_i) + i\epsilon} \delta V_{\text{eff}}^\mu(\omega). \qquad (30) $$

Having the density change in terms of electron-hole pairs (30) we can define a transition-resolved optical polarisability

$$ \alpha_{ij}(\omega) = \delta n_{ij}(\omega) \int \Psi_j^*(\mathbf{r}) \Psi_i(\mathbf{r}) d\mathbf{r} \qquad (31) $$

and thus assess the contribution of each non-interacting pair of states to the true, interacting polarisability. Moreover, it is now possible to answer questions related to the symmetry of the charge density with a strongest contribution at a given frequency [38] and perform other types of analysis analogous to the crystal orbital overlap and Hamiltonian populations analysis of the density matrix [54]. Here, we will focus on the analysis of the polarisability in terms of the dominant atomic angular-momentum contributions in the initial state. For this, we will explicitly separate the occupied and virtual states and look for the electron-hole expansion in the form

$$ \delta n(\mathbf{r}, \omega) = \sum_{i \in \text{occ}, j \in \text{unocc}} \delta n_{ij}(\omega) \Psi_j^*(\mathbf{r}) \Psi_i(\mathbf{r}), \qquad (32) $$

where the expansion coefficients $\delta n_{ij}(\omega)$ read

$$ \delta n_{ij}(\omega) = \left( \frac{(f_i - f_j)(X_a^i V_\mu^{ab} X_b^j)}{\omega - (E_j - E_i) + i\epsilon} - \frac{(f_i - f_j)(X_a^i V_\mu^{ab} X_b^j)}{\omega + (E_j - E_i) + i\epsilon} \right) \delta V_0^\mu(\omega). \qquad (33) $$

Now we can define a polarisability resolved in the angular momentum $l$ of the occupied states

$$ \alpha_l(\omega) = \sum_{i \in \text{occ}, j \in \text{unocc}} \delta n_{ij}(\omega) X_a^j d^{ab} \delta_{l_a, l} X_a^i. \qquad (34) $$

Obviously, the partial polarisability $\alpha_l(\omega)$ add up to the total interacting polarisability $\alpha(\omega) = \sum_l \alpha_l(\omega)$ and each of the partial polarisabilities gives an idea of the contribution of a given symmetry in the occupied states to the total polarisability. The result of this analysis for silver clusters is presented below in section 3.5.

## 3. Results

In this work, we focus our attention on silver clusters and shells, which are well-known plasmonic systems. We will address the dependence of the plasmonic resonances on the system size, the thickness of the shells and the details of the cluster geometry (relaxation method).

### 3.1. Calculation parameters

Clusters of icosahedral shapes were constructed using the atomic simulation environment [91] according to a Mackay motif [92]. The initial atomic positions were obtained using a 4.0 Å lattice constant that is close to the GGA-relaxed geometry and is smaller than the experimental lattice constant for bulk Ag (4.09 Å). The ‘ideal’ geometries were relaxed by minimising the forces acting on atoms below 0.02 eV Å⁻¹ using the GGA functional after Wu and Cohen [72, 78], as we discussed in section 2.5. In other words, we
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Relaxed geometries of icosahedral silver clusters. Labels S3, S4 and S5 refer to the number of atom layers present in the clusters. Labels L1 and L2 refer to the number of layers that are kept in a given cluster shell. The structures in the lowest row represent compact clusters. The atoms whose relative distances are discussed in text are marked by coloured numbers.</figcaption></figure>

used similar parameters in all the calculations. The spatial extension of the orbitals was set in a default procedure by an EnergyShift parameter of 50 meV and a double-zeta polarised set of atomic orbitals was used. The real-space mesh was set via a MeshCutoff parameter of 150 Ry. Only valence electrons (5s¹⁴d¹⁰) are represented in the atomic orbitals set, which give rise to 15 atomic orbitals per atom in the DZP basis set. Moreover, we added a layer of floating orbitals (see section 2.6) in the cluster calculations and two layers of floating orbitals, inner and outer, in the case of shell geometries. The coordinates of these 'ghost atoms' were kept fixed during relaxations.

The core electrons are removed by using the pseudo-potential (PP) of the Troullier-Martins type. The PPs were generated with the ATOM program, part of the SIESTA distribution. The parameters for PP generation were taken from the SIESTA database [93] except for the use of Wu and Cohen functional. It is interesting to note that the optimised PP described in [94], which is supposed to provide a band structure in better agreement with all-electron calculations, failed to satisfactorily describe the optical absorption cross section. Indeed, the description is severely worsened, and the main plasmonic resonance disappears.

## 3.2. Relaxed geometries

A set of representative clusters and cluster shells are shown in figure 3. We characterise icosahedral clusters by the number of atom layers present in the cluster and refer to this number as size of the cluster. Cluster shells are constructed starting from a cluster and keeping only several outer atomic layers; similar to the approach adopted by other groups [28, 29]. The number of atoms in a particular layer $l > 1$ is given by $N(l) = 10l^2 - 20l + 12$. The first 'layer' is composed of one atom. For example, the largest cluster we considered is composed of six atom layers (denoted S6L6), which makes in total $1 + 12 + 42 + 92 + 162 + 252 = 561$ atoms, while the icosahedral shell S7L4 will contain $92 + 162 + 252 + 362 = 868$ atoms.

As already mentioned, we optimise the geometries of the clusters and shells in order to account for the rearrangement effects caused by surface stresses. Geometrical relaxations were done by minimising the total forces, as implemented in SIESTA. The DFT relaxation only slightly compressed the ideal Mackay structures. For instance, the distance between extreme atoms 1 and 2 marked in figure 3 for the S4L4 cluster is 14.76 and 14.52 Å for the ideal geometry (in which the
<figure><img src="image_5.png" /><figcaption>**Figure 4.** The absorption cross sections of silver clusters of icosahedral shape are shown in panel (a). One can recognise resonances around 3–4 eV (sharp) and 6–7 eV (very broad). Panels (b) and (c) show the isosurfaces of density change Re(δn(r, ω)) for the Ag<sub>147</sub> cluster at the frequencies 3.4 and 6.6 eV, respectively.</figcaption></figure>

experimental lattice constant 4.09 Å is used) and the Wu and
Cohen geometry, respectively.

The surface stress also tends to distort (round up) the ideal
geometries. For instance, the distance between atoms 1 and
2 in the S5L5 cluster is slightly different from the distance
between extreme atoms 3 and 4 (see figure 3). The distance
d<sub>12</sub> is 19.68 and 19.40 Å for the ideal and the GGA geom-
etries, respectively, while the distance d<sub>34</sub> is 19.23 Å for the
GGA geometry.

Both geometry distortions (overall compression and rounding) are present in the silver shells. For instance, the GGA-relaxed length $d_{12}$ in the S4L1 shell is 13.79 Å. Distances $d_{12}$ and $d_{34}$ in the S5L1 shell are 18.18 and 18.14 Å, respectively. Note that these distances are smaller than those in the corresponding compact structures, indicating a larger compression in the case of mono-layered shells.

Summarising the outcome of the GGA geometric relaxa-
tions, we note that relaxation leads to minor geometrical dist-
ortions of the ideal icosahedral clusters. However, this leads
to relatively large compressions for thin silver shells. For
instance, if we characterise the compression with an averaged
bond length at the edge of the clusters (chosen as a simple
measure and representative for the other nearest neighbour
distances in the cluster), then the mono-layered shells have
that bond length compressed by 6% (2.80 Å) as compared to
the average bond length in the compact clusters (2.97 Å). The
cluster averaged bond length is also compressed by 2.4%, as
compared to the ideal bond length (3.04 Å) calculated with the
experimental lattice constant of bulk silver.

3.3. Optical absorption of Ag<sub>N</sub> icosahedral clusters

Qualitatively, one expects that the plasmon resonance is
affected by the size and morphology of the clusters. In
order to assess the magnitude of this dependence, we per-
formed TDDFT calculations of the absorption spectra for
compact clusters and shells. The smallest cluster is com-
posed of two layers (Ag<sub>13</sub>), while the largest cluster con-
sists of six layers (Ag<sub>561</sub>). In order to assess the effect of
removing internal atoms, we also performed calculations of

hollow structures, keeping up to four outer atomic layers in
a given system.

In order to expand the response and induced density we
have chosen an atom-centered product basis, as described in
section 2 with 77 functions per atom. The dominant product
basis [42, 64, 65] is used as an intermediate basis in the appli-
cation of the non-interacting response, as explained above.
Depending on the system, the size of the dominant product
basis is 2.4–4.1 times larger than that of the atom-centered
product basis. The generation of the product basis A<sup>μ</sup>(r) and
the calculation of interaction matrices K<sup>μν</sup> takes a relatively
small amount of time. For instance, in the case of our largest
Ag<sub>561</sub> cluster, a 12-core Intel machine spends (approximately)
45, 1 and 260 min, respectively, for the basis set generation,
and the calculations of Hartree and GGA kernels. The iterative
procedure normally takes most of the walltime. We decided
to compute the absorption spectra in a range 0–10 eV, with
a frequency step Δω = 0.05 eV, and a broadening constant
ε = 0.08 eV (i.e. the full-width at half maximum is 0.16 eV).
This consistent choice of frequency step and broadening
ensures that we do not ‘overlook’ any feature present in
the computed data and also produces data that can be well
interpolated.

We focus our attention on the optical absorption properties
of compact icosahedral geometries and then will move on to
a comparison with silver shells in the next section 3.4. Photo-
absorption cross sections are shown in figure 4(a). The cross
section possesses two maxima: a sharp peak at 3–4 eV and
a broad maximum at around 6–7 eV. The frequencies of the
resonances decrease for larger clusters. For instance, the two
largest clusters Ag<sub>309</sub> and Ag<sub>561</sub> have the absorption maxima
at 3.37 and 3.25 eV in the low-frequency band and at 6.7 and
6.2 eV in the high-frequency band, respectively. Panels (b)
and (c) in figure 4 show the induced density change (solution
of equation (21)) in the Ag<sub>147</sub> cluster at the maxima of the low-
frequency and high-frequency bands, respectively. The direc-
tion of the external field is set along the x-axis, i.e. collinear
with the plot plane and horizontal. The isosurfaces of the real
part of δn(r, ω) were plotted for a 10% of the corresponding
maximal value. One can see that the low-frequency resonance
<figure><img src="image_6.png" /><figcaption>**Figure 5.** The absorption cross sections of icosahedral silver shells constructed from the Ag<sub>309</sub> cluster. In panel (a) the cross section of the Ag<sub>309</sub> cluster is shown together with the cross sections of the one- to four-layered shells. In panel (b) the optical absorption cross sections of the ‘cutout’ and relaxed single-layered silver shells are compared (see text for detail). Panels (c) and (d) show the cross sections of the cutout and relaxed two and three-layered silver shells, respectively.</figcaption></figure>

is caused by an oscillation of charge with a pronounced dipole character, while the high-frequency band is supported by more localised, homogeneously-spread oscillations. The lower energy resonance obviously corresponds to the dipole Mie plasmon of the particle. From a quantum mechanical point of view, both bands consist of many nearly-degenerate transitions. The number of the transitions makes it impractical to analyse each of them in detail. However, one can provide some analysis of the electronic transitions and atom-layer contributions, as demonstrated below in section 3.5.

### 3.4. Optical absorption of silver shells

Silver shells can be constructed from the initial geometry of the corresponding cluster in several ways. We suggest two useful ways to construct the silver shells. In the first way, we simply delete atoms of several inner atomic layers from the compact cluster that was previously relaxed within the DFT. In the second way, we additionally relax the positions of the atoms remaining in the shell. The former way is less computationally demanding and is also eventually more useful to approximate the response of the whole cluster by the response of its outer shell. The latter way should give results that are closer to the corresponding experimental values for the shells.

In figure 5, we analyse the optical absorption cross section of the Ag<sub>309</sub> cluster and of the silver shells derived from that cluster. In panel (a) we show the cross sections of the relaxed shells. It is worth noting that a single-layered shell (S5L1, Ag<sub>162</sub>) exhibits a low-frequency resonance at 2.8 eV, which must be compared to the 3.37 eV in the case of the Ag<sub>309</sub> cluster. If we leave two atomic layers as in the silver shell S5L2 (Ag<sub>254</sub>), then the low-frequency band shifts back (3.28 eV) almost to the frequency of the full cluster. The broad high-frequency resonance of the silver shells qualitatively follows the same behaviour. Namely, the maximum frequency of the single-layered shell is red-shifted to 5.8 eV, which must be compared to 6.7 eV in the case of the Ag<sub>309</sub> cluster. However, already a two-layered shell almost recovers (6.5 eV) the position of the maximum of the high-frequency band for the compact cluster. The optical cross sections of the three and four-layered shells approach that of the cluster steadily.

The effect of geometry relaxation in the silver shells is also quantified in figure 5. Panels (b)–(d) show the optical absorption cross sections of the one, two and three-layered silver shells, the geometry of which was either ‘cut out’ from the geometry of the Ag<sub>309</sub> cluster or optimised (relaxed). The comparison shows that the effect of geometry relaxations is much larger for the single-layered shell than for the other shells.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** The absorption cross sections of icosahedral silver shells constructed from the Ag<sub>561</sub> cluster. In panel (a) the cross section of the Ag<sub>561</sub> cluster is shown together with the cross sections of the one to three-layered shells. In panel (b) the optical absorption cross sections of cutout and relaxed single-layered silver shells are compared. Panels (c) and (d) show the cross sections of the cutout and relaxed two and three-layered silver shells, respectively.</figcaption></figure>

The relaxation leads to a slight blue shift of both bands. This is due to the additional compression of the structure when relaxed. Moreover, the relaxations result in a broadening of the low-frequency resonance and in a relative increase in the cross section in the high-frequency band in the single-layered shell (panel (a)). Panels (c) and (d) show that the effect of geometry relaxation is less important for thicker shells. Cutout structures, created by simply removing the internal layers of atoms, give results very similar to those of the relaxed shells.

In figure 6 we show the optical absorption cross sections of the six-layered icosahedral silver cluster and shells constructed from that cluster. The cross sections show qualitatively the same behaviour as those of the five-layered cluster and shells (see figure 6 panels (a)–(d)). Namely, the cutout and relaxed one-layered shells S6L1 have the red-shifted low-frequency resonance at 2.6 and 2.62 eV, correspondingly, while the cross sections of the two and three-layered shells are much closer to that of the full cluster (3.25 eV). However, the low-frequency band of the two-layered shells S6L2 differs more from that of the compact cluster (3.15 eV for cutout and relaxed geometries). It is also interesting to note that the geometry relaxations of the one and two-layered shells lead to a small blue shift relatively to the unrelaxed calculations. In the case of three-layered shell S6L3, however, the geometry relaxations

lead to a slight red shift (3.25 eV) of the low-frequency maximum in comparison to the ‘cut out’ geometry (3.26 eV).

In figure 7 we collected the optical absorption of all the single (panel (a)) and double-layered (panel (b)) shells computed in this work. The cross section of the shells is qualitatively similar to the cross section of the clusters: there are low-frequency sharp and high-frequency broad resonances, the maxima of these resonances steadily decreases with increasing the cluster size. The ratio of the low-frequency to high-frequency intensities is nearly the same in compact structures and shells of sizes up to 6, while for the largest single-layered shell S7L1 the relative intensity of the low-frequency peak increases. The low-frequency peak in the absorption cross section of the S7L1 shell is more intense than for the S7L2 shell (figure 7, panel (c), in contrast to what is observed in smaller shells (figures 5 and 6). The absorption cross section of the thicker largest shells (S7L3 and S7L4) should approach the response of the compact cluster of the same size. Unfortunately, we could not compute the absorption of the cluster Ag<sub>943</sub>, because of the large memory requirements, but we can be reasonably sure that the thickest shell S7L4 represents well the absorption of the compact cluster of size 7. The frequency of maximal absorption shows an approximately linear dependency on
<figure><img src="image_8.png" /><figcaption>**Figure 7.** The absorption cross sections of the single (panel (a)) and double-layered (panel (b)) icosahedral silver shells. In panel (c) the absorption cross sections of shells of size 7 are compared. The plasmon frequency versus the size of the clusters and shells in panel (d). The lines in panel (d) are drawn through the points corresponding to largest structures (S6 and S7).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Contribution of the different atom layers to the absorption cross section of the Ag<sub>309</sub> cluster (S5L5).</figcaption></figure>

the inverse diameter of the clusters, as shown in panel (d) both for the clusters and shells. The effective diameter of the clusters and shells has been computed from the spatial cross section of the cluster, while the diameter of the cluster is defined as the diameter of a circle with the same area. The

data shown in figure 7 (d) indicate small discrepancies from the linear trend both for compact geometries and double-layered shells for smaller clusters, while in the case of the single-layered shells this discrepancy is visible only for the smallest shell.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Analysis of the different angular momenta of the basis orbitals involved in the description of the occupied states contributing to the absorption cross section of the Ag<sub>309</sub> (S5L5) cluster as a function of frequency.</figcaption></figure>

**3.5. Analysis of the optical absorption**

In figure 8 we show the analysis of the atom-layer contrib-
utions to the optical absorption cross section according to
the method presented above in section 2.7. The partial cross
sections corresponding to different atom-layer polarisabili-
ties (24) are plotted together with the total absorption cross
section. The outer layer (5th layer, 162 atoms) gives rise to
the largest partial cross section, which is even larger than
the total cross section at resonance. The 6th layer of ‘ghost
atoms’ contributes constructively in the whole frequency
range and with a rather large magnitude that is approximately
equal to that of the second inner layer (4th layer) of ‘true’
atoms. The inner atomic layers give rise to negative partial
cross sections, at least in some frequency ranges. At the main
resonance 3.37 eV, all the inner layers contribute destructively.
The second outer layer (4th layer, 92 atoms) gives rise to a
negative partial cross section in most of the frequency range
considered. The contribution of the first ‘layer’, which is com-
posed of one atom is much lower than that of the next layer
(2nd layer, 12 atoms) and is not shown.

The interpretation of these results in the low-frequency
range is quite clear. As expected from the plot in figure 4(b),
and consistent with the Mie plasmon character of the reso-
nance, the main contribution comes from the surface layer.
For this mode, it is critical to describe accurately the polarisa-
tion of the surface, which explains the large contribution of
the layer of ‘ghost atoms’ above the cluster surface. Those
basis orbitals control the extension of the electronic states
towards a vacuum and, therefore, are instrumental to correctly
account for the polarisability of the surface. The change of
signs is related to the relative phases of the different contrib-
utions, reflecting to some extent the global nodal structure of
the contributing electronic states.

At higher energies, the smaller dielectric function of bulk
silver makes the charge screening less efficient, the modes
loose their predominant surface character, and the contrib-
utions of different layers become more similar in intensity.

In figure 9 the partial cross sections (corresponding to the
partial polarisabilities (34)) for the atomic orbitals of different
symmetry expanding the *occupied* states are shown for the
Ag<sub>309</sub> cluster. We see that, on the absolute scale, the contrib-
ution of each angular momentum channel is important even
at low frequencies. However, there is a striking difference
between the *s*, *p* and *d* contributions to the absorption cross
section. Namely, the *s*- and *p*-channels contribute construc-
tively to the absorption cross section, while the *d* channel con-
tributes destructively in the low frequency range (3–4 eV) and
only starts to contribute constructively at high frequencies,
starting approximately at 5.5 eV. This conclusion is similar
to that presented in a previous LDA study for smaller silver
clusters [73]. This theoretical analysis supports the view of the
low-frequency excitation as only partially produced by 5s elec-
trons [95]. This is because the occupied states are *sp*-hybrid-
ised close to Fermi energy, and there is also a non-negligible
contribution of *s*-symmetrical density at higher frequencies
(6–7 eV). From the other side, the strong over-screening of the
*s* plasmon due to *sp*-d interband transitions has been discussed
at length in the literature [74, 80]. The comparison of our data
with those in [11, 74] suggests that GGA could underestimate
the intensity of *s* plasmon by a factor about 1.5.

In order to comprehend better the outcome of *ab initio*
modelling we compare our results with the available relevant
literature in the following section.

**4. Discussion**

Because the clusters that we have considered are rather large,
we expect that their properties will approach the properties
of classical Mie spheres [7, 73] to the extent that the non-
sphericity of icosahedra and charge spillage effects allow.
The classical absorption cross section in the electrostatic limit
depends on the dielectric function of the material $\epsilon(\omega)$ and the
dielectric function of the embedding medium $\epsilon_m(\omega)$ [7, 96]
<figure><img src="image_11.png" /><figcaption>Figure 10. The absorption cross section in a classical quasi-static limit for a silver sphere in vacuum. The experimentally-determined dielectric function of silver is used.</figcaption></figure>

$$ \sigma(\omega) \sim \omega \operatorname{Im} \frac{\epsilon(\omega) - \epsilon_m(\omega)}{\epsilon(\omega) + 2\epsilon_m(\omega)}. \qquad (35) $$

If we plot the cross section (35) with the experimentally determined dielectric function of silver [97, 98], then we will get a sharp low-frequency resonance and a less intense and broad high-frequency resonance, as shown in figure 10 panel (a). The position of the resonances slightly varies with the set of experimental data used. The more recent experiment [98] gives resonance maxima at 3.53 and 6.8 eV, while the older and widely cited experimental results [97] give them at 3.47 and 6.1 eV, correspondingly. The response in the low-frequency part of the spectra is detailed in panel (b) in figure 10.

Comparing the quasi-static absorption cross section (35) with our results, we see similarities and discrepancies. First of all, the resonance frequencies from Mie theory do not depend on the size of the object, but only on its shape. Classical estimations of the shape influence [11] lead to a conclusion that icosahedral particles have to resonate at 2.9% lower frequencies than perfectly spherical clusters.

TDDFT accounts for charge spillage at surfaces (i.e. the fact that surfaces are not perfectly sharp, as assumed usually in classical electrodynamics) which leads to a divergence of the plasmon frequency $\omega_{sp}$ from a linear dependence $\omega_{sp} = kd^{-1} + b$ with respect to the inverse diameter for small clusters (see figure 7(d)). Measurements of the optical absorption of small silver clusters support this conclusion, at least qualitatively [24, 99–101]. The divergence from the linear dependence in the experimental data must be due to charge spillage, because the other effect that leads to red shifts, higher multipole contributions to the electron–photon coupling, only affects relatively large clusters of diameters 50 nm and more [102].

In the case of the smallest cluster presented here Ag$_{13}$, there are experimental data [103, 104] for argon embedded clusters that show a maximum absorption strength at 3.4–3.5 eV, while we get 3.63 eV. Because the effect of argon embedding can lead to red shifts as large as 0.3–0.5 eV [101, 105], we estimate that TDDFT with a GGA functional might be delivering red-shifted frequencies by about 0.2 eV. In the case of larger clusters, we can compare with the experimental data [24] measured on free silver clusters. For instance, for a cluster

diameter of 1.2 nm, which corresponds approximately to our Ag$_{55}$ cluster, the experimental peak position is 3.51 eV, while we get 3.7 eV. For a cluster diameter of 2.2 nm, which corresponds to our Ag$_{309}$ cluster, the experimental peak position is 3.47, while we get 3.37 eV. For a cluster diameter of 2.8 nm, which corresponds to our Ag$_{561}$ cluster, the experimental peak position is 3.45, while we get 3.25 eV. These comparisons indicate that our calculations deviate by no more than 0.2 eV from the experiment, presenting a somewhat stronger size dependency than that observed experimentally. The computed plasmon frequencies are red-shifted for large clusters.

Besides the experimental results, there are many calculations of silver optical absorption properties available in the literature. The first observation is that very little can be found on the broad resonance in the frequency range 5.5–7.5 eV, although this feature must be available in the calculations done with atomistic codes. This is probably due to little practical relevance of the frequency range and the much stronger low-frequency plasmonic response.

The position and strength of the low-frequency plasmonic peak depends on the used functional and the basis set: the agreement with experiment improves while using more sophisticated functionals and larger basis sets. For instance, the use of the simplest LDA functional results in a too low energy onset of the d-bands in the electronic structure, [79] which leads to a reduced strength of the low-frequency resonance and its red shift as compared to the experiment and calculations based on Hedin's GW band structures [79]. In this paper, we showed that GGA functional by Wu and Cohen produces a strong low-frequency peak for icosahedral clusters. GGA functionals have been extensively used in the past for silver clusters, [36, 104] as well as the more sophisticated long-range-corrected functionals: van Leeuwen–Baerends (LB94) [30], Gritsenko–van Leeuwen–van Lenthe–Baerends (GLLB) [11] and long-range corrected PBE (LC-ωPBE) [37, 106, 107].

Quantitatively assessing the position of the low-frequency resonance, we can compare our results with several theoretical calculations. For the Ag$_{13}$ cluster, we can find the values (in eV) 3.5 (PBE) [104], 3.5 (LC-ωPBE) [37], 3.2 (LDA) [38], 3.7 (LDA) [39] and 3.2 (PBE) [108], while we get the maximum at 3.63 eV. For the Ag$_{55}$ icosahedral cluster, we can find values
3.5 (LDA) [39] and 4.2 (GLLB) [11] eV, while we get the maximum at 3.71 eV. For the Ag<sub>147</sub> icosahedral cluster, we can find values 3.8 (GLLB) [11], 3.2 (PBE) [36], 4.5 (LB94) [30], while we get the maximum at 3.5 eV. For Ag<sub>309</sub> icosahedral cluster, we can find values 3.7 (GLLB) [11], 3.50 (polarisability interaction model) [109], while we get the maximum at 3.37 eV. Finally, for the Ag<sub>561</sub> icosahedral cluster, there is a value 3.65 eV from [11], while we get 3.25 eV. Summarising, we see that the deviations between our calculations and other calculations that utilise the LDA or GGA functionals do not exceed 0.2 eV and might be due to differences in the used basis sets or pseudo-potentials. The larger deviations of about 0.4 eV, we have with LRC functionals.

Hollow structures have been produced and characterised [14, 110–113]. The experimental evidence supports our finding that the hollow structures have lower frequencies of the main plasmonic resonance. However, the measurements were performed with relatively large and thick shells (at least one order of magnitude larger than those considered here) and a quantitative comparison of our results with experiments is hardly possible.

From the theoretical side, there are several estimations for the hollow metallic particles available. First, the solution of the Maxwell equations in the quasi-static limit for a hollow spherical particle is known [96] and the effect of removing the interior part of the sphere is a red shift of the plasmon frequency. As in the case of compact spheres, the resonance frequencies do not depend on the size of the system, but only on the ratio of the inner and outer radius of the sphere. Second, there are atomistic TDDFT calculations for icosahedral shells of up to a size of six layers available [28, 29]. These calculations show red-shifted plasmonic resonances for hollow clusters compared to the filled ones. Quantitatively comparing the position of the low-frequency resonance, we get a fair agreement. For instance, for Ag<sub>92</sub> (S4L1 shell) we can find 3.85 (LB94) [29] and 2.8 (PBE) [28] eV, while we get 2.97 eV. For Ag<sub>12</sub>, Ag<sub>42</sub>, Ag<sub>162</sub> and Ag<sub>252</sub> (S2L1, S3L1, S5L1 and S6L1 shells) we extract from [29] 4.0, 4.35, 3.5 and 3.2 eV, while we get 4.17, 3.13, 2.8 and 2.63 eV, respectively. The discrepancies are sizeable, but can be explained by the differences of the functional (LB94 versus WC) and the geometry relaxations (ideal symmetric with the nearest-neighbour distance fixed at a equilibrium bulk value 2.89 Å versus fully relaxed geometries with no symmetry imposed in our case). The effect of geometry relaxations is estimated in section 3.4: it has a minor importance compared to the influence of the DFT functional and blue-shifts the resonances of the ideal structures with respect to those of the relaxed ones.

The response of other silver-made structures with effec-
tively reduced dimensionality: rods [36, 37, 114, 115] and
platelets [12, 116] also exhibit red-shifted resonance frequen-
cies of the plasmon excitations with respect to spherical or
quasi-spherical clusters of similar sizes.

5. Conclusion

We studied the optical response of silver clusters of icosahedral symmetry and their hollow counterparts by means of quantum mechanical, atomistic methods (linear response TDDFT within LCAO). The applied iterative methods allowed for comparatively fast calculations (about 1 day of walltime)

of compact clusters containing up to 561 atoms and hollow
clusters of up to 868 atoms. We found that the plasmonic reso-
nance of silver clusters depends on the size and morphology
of the clusters. Namely, the frequency of the maximal absorp-
tion of the icosahedral clusters that contain three and more
atom layers is inversely proportional to the cluster diameter.
Moreover, the single-layered shells show a sizeable red shift
of the resonance frequencies, which quickly becomes neg-
ligible as the thickness of the shells increases. Both obser-
vations are compatible with the experimental findings and
previous calculations, and can be partially understood within
classical electrodynamics. Furthermore, both observations are
valid also for the high-frequency, interband plasmon, which
has not been widely studied so far.

From a methodological point of view, we presented recent
developments of our iterative technique including a realisa-
tion of an atom-centered product basis and various analysis
tools. The iterative method used here [42], particularly with
the speedup allowed by the use of an atom-centered auxiliary
basis to express the orbital products, is advantageous in many
respects, which are discussed in this paper. In particular, the fre-
quency-range selectivity is useful in calculations of the Raman
response, for which one needs the response in a very narrow
spectral range. Moreover, the current implementation of the
response times vector operation can be easily optimised to be
less memory demanding. This optimisation will allow an exten-
sion of the number of treated atoms by an order of magnitude,
and even further if full MPI parallelisation is implemented.

Acknowledgments

This work is supported, in part, by the ORGAVOLT (ORGANIC solar cell VOLTAGE by numerical computation) Grant ANR-12-MONU-0014-02 of the French Agence Nationale de la Recherche (ANR) 2012 Programme Modèles Numériques. F Marchesin, P Koval and D Sánchez-Portal acknowledge support from the Deutsche Forschungsgemeinschaft (DFG) through the sFB1083 project, the Spanish MINECO MAT2013-46593-C6-2-P project, the Euroregion Aquitaine-Euskadi program and from the Basque Departamento de Educación, UPV/EHU (Grant No. IT-756-13). Peter Koval acknowledges financial support from the Fellows Gipuzkoa program of the Gipuzkoako Foru Aldundia through the FEDER funding scheme of the European Union. Research is conducted in the scope of the Transnational Common Laboratory "QuantumChemPhys Theoretical Chemistry and Physics at the Quantum Scale".

Appendix A. Dominant basis set symmetries
contributing to the optical polarisability

The prior analysis in section 2.7 has been formulated for a
simple, physically motivated splitting of the interacting
polarisability in terms of the angular momentum of the
occupied states and in terms of the atomic contributions to
the optical polarisability. In this section, we focus on a more
technical analysis of the interacting polarisability in terms of
<figure><img src="image_12.png" /><figcaption>**Figure A1.** Analysis of the product-function angular-momentum contributions to the absorption cross section of the Ag<sub>309</sub> cluster. The total absorption cross section is also plotted for comparison.</figcaption></figure>

the products of atomic-orbital functions, and in terms of the product basis functions contributing to the resulting induced screened density change in the cluster.

Similar to the analysis tools considered in section 2.7, the total interacting polarisability can be split into sums with fixed angular momentum of the product functions $F^{\mu}(\mathbf{r})$. This is so because the product functions are constructed as a linear combination of the products of atomic orbitals separately for each possible angular momentum of the product [65] and, therefore, the product function $\mu$ carries a well-defined angular momentum $l_{\mu}$. We can write the polarisability as a sum over the angular momenta of the product basis $\alpha(\omega) = \sum_{l} \alpha_{l}(\omega)$ with a product angular momentum-resolved polarisability $\alpha_{l}(\omega)$ given by

$$ \alpha_l(\omega) = d^\mu \delta_{l_\mu, l} \delta n_\mu(\omega), \quad (A.1) $$

where the dipole moments $d^\mu$ refer to a global origin of the coordinate system. The dipole moments are defined for product functions $F^{\mu}(\mathbf{r})$ that are centered on atoms

$$ \mathbf{d}^{\mu} = \int F^{\mu}(\mathbf{r} - \mathbf{R}_{\mu})\mathbf{r} d^3\mathbf{r} = \int F^{\mu}(\mathbf{r})\mathbf{r} d^3\mathbf{r} + \mathbf{R}_{\mu} \int F^{\mu}(\mathbf{r}) d^3\mathbf{r}. \quad (A.2) $$

The last equation makes it apparent that only angular momentum $l_{\mu} = 0$ and $l_{\mu} = 1$ can contribute to the dipole polarisability.

Yet another type of analysis involves the angular momentum symmetry of the atomic orbitals, the products of which expand the density change $\delta n(\omega, \mathbf{r})$. Namely, the product function $F^{\mu}(\mathbf{r})$ is expressed in terms of a linear combination of the products of atomic orbitals (9). Therefore, we can immediately write the interacting polarisability $\alpha(\omega)$ as a sum over atomic orbital angular momentum sub-sums $\alpha(\omega) = \sum_{l_1,l_2} \alpha_{l_1,l_2}(\omega)$ with the partial polarisability $\alpha_{l_1,l_2}(\omega)$ given by

$$ \alpha_{l_1,l_2}(\omega) = \delta_{l_1,l_2} d^{ab} \delta_{l_2,l_1} A_{ab}^{\mu} \delta n_{\mu}(\omega). \quad (A.3) $$

Here the dipole matrix elements between atomic orbitals is used $d^{ab} = \int f^a(\mathbf{r} - \mathbf{R}_a)rf^b(\mathbf{r} - \mathbf{R}_b)d\mathbf{r}$. Thus, recalling that

the atom-centered functions are constructed from local on-site products, we see that the total polarisability can be expressed in terms of the on-site atomic orbitals. From the one side, this reveals a certain arbitrariness of the proposed separation, but from the other side, it might be useful for developing simplified models. Because the global origin and the atom centres do not generally coincide, there are no dipole selection rules in the matrix-elements $d^{ab}$. Namely, analogously to the case of dipole moments (A.2), the matrix elements $d^{ab}$ will depend on the overlap between orbitals

$$ \mathbf{d}^{ab} = \int f^a(\mathbf{r})rf^b(\mathbf{r})d\mathbf{r} + \mathbf{R} \int f^a(\mathbf{r})f^b(\mathbf{r})d\mathbf{r}, \quad (A.4) $$

where $\mathbf{R} = \mathbf{R}_a = \mathbf{R}_b$. The matrix element $d^{ab}$ is zero if both terms in the last equation are zero. This is always the case for sd channels ($l_a = 0, l_b = 2$), but not generally for ss and pp symmetries.

In figure A1, we present the analysis of the product-function angular momentum contributions to the optical absorption cross section. The partial cross sections corresponding to the partial polarisability (A.1) are plotted together with the total absorption cross section for the Ag<sub>309</sub> cluster. The contribution of the angular momentum higher than p is strictly zero. The contribution of the s-symmetric product functions dominates in the cross section, while the contribution of p-symmetric functions is negative for the low-frequency resonance in the range 3-4 eV and positive in the frequency range 5-8 eV.

In figure A2, we plot the partial cross sections computed from the orbital angular-momentum-resolved polarisability (A.3) for the 5-layered icosahedral cluster Ag<sub>309</sub>. The sum of the partial cross sections is also shown for comparison. We can clearly appreciate the contribution of the p-orbitals to the total absorption cross section. Namely, the contribution of the pp angular-momentum channel is the most significant. The contributions have different signs and, for example, pp-products give rise to a partial cross section that is about
<figure><img src="image_13.png" /><figcaption>**Figure A2.** Analysis of atomic orbital angular momenta of the products of atomic-orbital contributions to the optical absorption cross section of the Ag<sub>309</sub> cluster. The total absorption cross section is also plotted for comparison.</figcaption></figure>

three times larger than the total cross section. The second most important type of product is of ss type. The partial cross section of ss type is negative in the whole frequency range we computed. The absolute value of the ss partial cross section is about 1.7 times larger than the total cross section. The other non-zero angular-momentum channels do not contribute significantly to the total absorption cross section. Therefore, a linear combination of the pp and ss products of atomic orbitals seems to be able to provide a reasonably accurate description of the total absorption cross section in the frequency range 0–10 eV and even better so in the frequency range 0–4 eV.

## Appendix B. GGA interaction kernel

In the generalised gradient approximation (GGA) of the DFT functional, the xc energy density $\epsilon$ explicitly depends on the charge density $n(\mathbf{r})$ and its gradient $\nabla n(\mathbf{r})$

$$E = \int \epsilon(n(\mathbf{r}), \nabla n(\mathbf{r})) d^3r. \qquad (B.1)$$

In the Kohn–Sham formulation of the DFT we need the functional derivative of the energy with respect to the electron density

$$v(\mathbf{r}) \equiv \frac{\delta E}{\delta n(\mathbf{r})} = \frac{\partial \epsilon}{\partial n} - \nabla \frac{\partial \epsilon}{\partial \nabla n}. \qquad (B.2)$$

and in the linear response TDDFT we also need a second variational derivative of the energy with respect to the electronic density

$$f(\mathbf{r}, \mathbf{r}') \equiv \frac{\delta v(\mathbf{r})}{\delta n(\mathbf{r}')}. \qquad (B.3)$$

While the calculation of the KS potential in the GGA is widely discussed in the literature [117–119], the TDDFT kernel (B.3) is less commonly discussed [120, 121] and we find it

convenient to state explicitly the equations necessary for the implementation of our iterative TDDFT method.

### B.1. Variational derivative of potential

In linear response TDDFT one needs a first variational derivative of the xc potential with respect to density (B.3). The potential (B.2) is in fact a functional of the density and its gradient. Computing the variation of the potential one gets

$$\delta v(\mathbf{r}) = \frac{\partial^2 \epsilon}{\partial n \partial n} \delta n + \frac{\partial^2 \epsilon}{\partial n \partial \nabla_i n} \nabla_i \delta n - \nabla_i \left( \frac{\partial^2 \epsilon}{\partial n \partial \nabla_i n} \delta n \right) - \nabla_i \left( \frac{\partial^2 \epsilon}{\partial \nabla_k n \partial \nabla_i n} \nabla_k \delta n \right). \qquad (B.4)$$

In order to ‘exchange’ the gradients of density variation $\nabla_i \delta n$ and $\nabla_k \nabla_i \delta n$ for the density variation $\delta n$, we represent the last equation in form of an integral with a $\delta$-function $\delta v(\mathbf{r}') = \int \delta(\mathbf{r}' - \mathbf{r}) \delta v(\mathbf{r}) d^3r$ and apply the integration by parts whenever necessary. Finally, we can get [122]

$$\delta v(\mathbf{r}') = \int \delta(\mathbf{r}' - \mathbf{r}) \delta n \frac{\partial^2 \epsilon}{\partial n \partial n} - \delta(\mathbf{r}' - \mathbf{r}) \delta n \nabla_i \left[ \frac{\partial^2 \epsilon}{\partial n \partial \nabla_i n} \right] - \delta n \nabla_k \left[ \nabla_i [\delta(\mathbf{r}' - \mathbf{r})] \frac{\partial^2 \epsilon}{\partial \nabla_k n \partial \nabla_i n} \right] d^3r. \qquad (B.5)$$

By virtue of the kernel definition (B.3), the kernel GGA reads [123]

$$f(\mathbf{r}, \mathbf{r}') = \delta(\mathbf{r} - \mathbf{r}') \frac{\partial^2 \epsilon}{\partial n \partial n} - \delta(\mathbf{r} - \mathbf{r}') \nabla_i \left[ \frac{\partial^2 \epsilon}{\partial n \partial \nabla_i n} \right] - \nabla_k \left[ \nabla_i [\delta(\mathbf{r} - \mathbf{r}')] \frac{\partial^2 \epsilon}{\partial \nabla_k n \partial \nabla_i n} \right]. \qquad (B.6)$$
## B.2. Computation of GGA matrix elements of potential

In DFT, we need matrix elements of the xc potential (B.2)

$$v^{ab} = \int f^a(\mathbf{r}) \frac{\partial \epsilon}{\partial n} f^b(\mathbf{r}) d^3\mathbf{r} - \int f^a(\mathbf{r}) \nabla \frac{\partial \epsilon}{\partial \nabla n} f^b(\mathbf{r}) d^3\mathbf{r} \quad (B.7)$$

between the atomic orbitals $f^a(\mathbf{r})$. Because we use LIBXC library [124], we need to work with derivatives of the xc energy density with respect to the square of the density gradient, i.e. in terms of $\sigma \equiv (\nabla n)^2$. After applying the chain rule, we transform the last formula to

$$v^{ab} = \int f^a(\mathbf{r}) \frac{\partial \epsilon}{\partial n} f^b(\mathbf{r}) d^3\mathbf{r} - 2 \int f^a(\mathbf{r}) \nabla_i \frac{\partial \epsilon}{\partial \sigma} \nabla_i n [f^b(\mathbf{r})] d^3\mathbf{r}, \quad (B.8)$$

Here, and in the following, the summation over repeated Cartesian indices ($i$ and $k$) is understood. The expression (B.8), similar to that in equation (B.7), is inconvenient in a numerical calculation because one must compute derivatives of the quantities delivered by LIBXC library ($\frac{\partial \epsilon}{\partial \nabla n}$ or $\frac{\partial \epsilon}{\partial \sigma} \nabla n$). However, with one more integration by parts, we can transform the expression (B.8) into a more suitable form

$$v^{ab} = \int f^a(\mathbf{r}) \frac{\partial \epsilon}{\partial n} f^b(\mathbf{r}) d^3\mathbf{r} \\ + 2 \int \nabla_i [f^a(\mathbf{r})] \nabla_i n \frac{\partial \epsilon}{\partial \sigma} f^b(\mathbf{r}) d^3\mathbf{r} \\ + 2 \int f^a(\mathbf{r}) \frac{\partial \epsilon}{\partial \sigma} \nabla_i [f^b(\mathbf{r})] \nabla_i n d^3\mathbf{r}. \quad (B.9)$$

This form is advantageous because it is easier to calculate the gradients of the basis functions, the form of which is known and does not change with the type of xc functional. Moreover, the last expression uses the same ingredients that are necessary to compute the matrix elements of the GGA kernel.

## B.3. Computation of the GGA matrix elements of the kernel

In our implementation of the linear response TDDFT, we need the matrix elements of the xc kernel (B.6) defined by $K^{\mu\nu} = \int F^{\mu}(\mathbf{r})f(\mathbf{r}, \mathbf{r}')F^{\nu}(\mathbf{r}')d^3\mathbf{r}'d^3\mathbf{r}$. Exercising the same approach as for the potential in the previous section, we obtain for matrix elements of the GGA kernel

$$K^{\mu\nu} = \int F^{\mu}(\mathbf{r}) \frac{\partial^2 \epsilon}{\partial n^2} F^{\nu}(\mathbf{r}) d^3\mathbf{r} \\ + 2 \int \nabla_i [F^{\mu}(\mathbf{r})] \nabla_i [n] \frac{\partial^2 \epsilon}{\partial n \partial \sigma} F^{\nu}(\mathbf{r}) d^3\mathbf{r} \\ + 2 \int F^{\mu}(\mathbf{r}) \frac{\partial^2 \epsilon}{\partial n \partial \sigma} \nabla_i [F^{\nu}(\mathbf{r})] \nabla_i [n] d^3\mathbf{r} \\ + 4 \int \nabla_k [F^{\mu}(\mathbf{r})] \nabla_k [n] \frac{\partial^2 \epsilon}{\partial \sigma^2} \nabla_i [F^{\nu}(\mathbf{r})] \nabla_i [n] d^3\mathbf{r} \\ + 2 \int \nabla_k [F^{\mu}(\mathbf{r})] \frac{\partial \epsilon}{\partial \sigma} \delta_{ik} \nabla_i [F^{\nu}(\mathbf{r})] d^3\mathbf{r}. \quad (B.10)$$

These matrix elements are computed with standard integration methods in quantum chemistry [125, 126].

## References

[1] de Heer W A 1993 *Rev. Mod. Phys.* **65** 611–76

[2] Newton D E 2015 *Silver. Chemical Elements: From Carbon to Krypton*. (Gale group)

[3] Brack M 1993 *Rev. Mod. Phys.* **65** 677–732

[4] Thämer M, Kartouzian A, Heister P, Lünskens T, Gerlach S and Heiz U 2014 *Small* **10** 2340–4

[5] Madison L R, Ratner M A and Schatz G C 2015 *Frontiers in Quantum Methods and Applications in Chemistry and Physics* (Progress in Theoretical Chemistry and Physics vol 29) ed M Nascimento *et al* (Cham: Springer) pp 37–52

[6] Cushing S K and Wu N 2016 *J. Phys. Chem. Lett.* **7** 666–75

[7] Kreibig U and Vollmer M 1995 *Optical Properties of Metal Clusters* (Springer Series in Materials Science vol 25) (Berlin: Springer)

[8] Hao E and Schatz G C 2004 *J. Chem. Phys.* **120** 357–66

[9] Koponen L, Tunturivuori L O, Puska M J and Hancock Y 2010 *J. Chem. Phys.* **132** 214301

[10] Bousquet B, Cherif M, Huang K and Rabilloud F 2015 *J. Phys. Chem. C* **119** 4268–77

[11] Kuisma M, Sakko A, Rossi T P, Larsen A H, Enkovaara J, Lehtovaara L and Rantala T T 2015 *Phys. Rev. B* **91** 115431

[12] Wiley B, Sun Y, Mayers B and Xia Y 2005 *Chem.: Eur. J.* **11** 454–63

[13] Chang R and Leung P T 2006 *Phys. Rev. B* **73** 125438

[14] Oldenburg S J 2015 wwwaldrich.com/materials-science/nanomaterials/silver-nanoparticles.html

[15] Mathew A and Pradeep T 2014 *Part. Part. Syst. Charact.* **31** 1017–53

[16] Shimizu K I and Satsuma A 2006 *Phys. Chem. Chem. Phys.* **8** 2677–95

[17] Leelavathi A, Bhaskara Rao T and Pradeep T 2011 *Nanoscale Res. Lett.* **6** 123

[18] Guável X L, Hötzer B, Jung G, Hollemeyer K, Trouillet V and Schneider M 2011 *J. Phys. Chem. C* **115** 10955–63

[19] Chakraborty I, Udayabhaskararao T and Pradeep T 2012 *Journal of Hazardous Materials* **211–2** 396–403 (nanotechnologies for the treatment of Water, Air and Soil www.sciencedirect.com/science/article/pii/S030438941101524X)

[20] Qu F, Li N B and Luo H Q 2012 *Anal. Chem.* **84** 10373–9

[21] Chakraborty I, Bag S, Landman U and Pradeep T 2013 *J. Phys. Chem. Lett.* **4** 2769–73

[22] Birke R L and Lombardi J R 2015 *J. Opt.* **17** 114004

[23] Hao F and Nordlander P 2007 *Chem. Phys. Lett.* **446** 115–8 (www.sciencedirect.com/science/article/pii/S0009261407010986)

[24] Charlé K P, König L, Nepijko S, Rabin I and Schulze W 1998 *Cryst. Res. Technol.* **33** 1085–96

[25] Negre C F A, Perassi E M, Coronado E A and Sánchez C G 2013 *J. Phys.: Condens. Matter* **25** 125304

[26] Barbry M, Koval P, Marchesin F, Esteban R, Borisov A G, Aizpurua J and Sánchez-Portal D 2015 *Nano Lett.* **15** 3410–9

[27] Aikens C M, Li S and Schatz G C 2008 *J. Phys. Chem. C* **112** 11272–9

[28] Weissker H C, Whetten R L and Lopez-Lozano X 2014 *Phys. Chem. Chem. Phys.* **16** 12495–502

[29] Barcaro G, Sementa L, Fortunelli A and Stener M 2014 *J. Phys. Chem. C* **118** 12450–8

[30] Barcaro G, Broyer M, Durante N, Fortunelli A and Stener M 2011 *J. Phys. Chem. C* **115** 24085–91

[31] Barcaro G, Sementa L, Fortunelli A and Stener M 2014 *J. Phys. Chem. C* **118** 28101–8

[32] Baetzold R C 1978 *J. Chem. Phys.* **68** 555–61

[33] Zhao J, Luo Y and Wang G 2001 *Eur. Phys. J. D: At. Mol. Opt. Plasma Phys.* **14** 309–16
[34] Bonačić-Koutecký V, Burda J, Mitrić R, Ge M, Zampella G and Fantucci P 2002 *J. Chem. Phys.* **117** 3120–31

[35] Johnson H E and Aikens C M 2009 *J. Phys. Chem. A* **113** 4445–50

[36] Lozano X L, Mottet C and Weissker H C 2013 *J. Phys. Chem. C* **117** 3062–8

[37] Heard C J and Johnston R L 2014 *Phys. Chem. Chem. Phys.* **16** 21039–48

[38] Baishya K, Idrobo J C, Öğüt S, Yang M, Jackson K and Jellinek J 2008 *Phys. Rev. B* **78** 075439

[39] Idrobo J C and Pantelides S T 2010 *Phys. Rev. B* **82** 085420

[40] Huda M N and Ray A K 2003 *Phys. Rev. A* **67** 013201

[41] Bae G T and Aikens C M 2012 *J. Phys. Chem. A* **116** 8260–9

[42] Koval P, Foerster D and Coulaud O 2010 *J. Chem. Theory Comput.* **6** 2654–68

[43] Koval P, Foerster D and Coulaud O 2010 *physica status solidi* b **247** 1841–8

[44] Soler J M, Artacho E, Gale J D, García A, Junquera J, Ordejón P and Sánchez-Portal D 2002 *J. Phys.: Condens. Matter* **14** 2745

[45] Artacho E, Anglada E, Diéguez O, Gale J D, García A, Junquera J, Martin R M, Ordejón P, Pruneda J M, Sánchez-Portal D and Soler J M 2008 *J. Phys.: Condens. Matter* **20** 064208

[46] Tsolakidis A, Sánchez-Portal D and Martin R M 2002 *Phys. Rev. B* **66** 235416

[47] Koval N E, Sánchez-Portal D, Borisov A G and Muiño R D M 2013 *Nucl. Instrum. Methods Phys. Res.* B **317** 56–60 (Proc. of the 19th Int. Workshop on Inelastic Ion-Surface Collisions (Frauenchiemsee, Germany, 16–21 September 2012))

[48] Iida K, Noda M, Ishimura K and Nobusada K 2014 *J. Phys. Chem. A* **118** 11317–22

[49] Xiang H, Zhang X, Neuhauser D and Lu G 2014 *J. Phys. Chem. Lett.* **5** 1163–9

[50] Gao Y, Neuhauser D, Baer R and Rabani E 2015 *J. Chem. Phys.* **142** 034106

[51] Andrade X, Botti S, Marques M A L and Rubio A 2007 *J. Chem. Phys.* **126** 184106

[52] Hübener H and Giustino F 2014 *Phys. Rev. B* **89** 085129

[53] Manjavacas A, Marchesin F, Thongrattanasiri S, Koval P, Nordlander P, Sánchez-Portal D and García de Abajo F J 2013 *ACS Nano* **7** 3635–43

[54] Dronskowski R 2005 *Computational Chemistry of Solid State Materials* (Weinheim: Wiley)

[55] Nogueira F, Marques M A L and Fiolhais C 2003 *A Primer in Density Functional Theory* (Lecture Notes in Physics) (Berlin: Springer)

[56] Gross E K U and Burke K 2006 *Basics of TDDFT* vol 706 (Berlin: Springer) chapter 1, p 1

[57] Petersilka M, Gossmann U J and Gross E 1996 *Phys. Rev. Lett.* **76** 1212

[58] Gross E 2010 *Time-Dependent Density-Functional Formalism Benasque 2010* (Summer School Lecture) (Benasque: European Theoretical Spectroscopy Facility) www.tddft.org/TDDFT2010/school/gross3.pdf

[59] Onida G, Reining L and Rubio A 2002 *Rev. Mod. Phys.* **74** 601–59

[60] Bauernschmitt R, Häser M, Treutler O and Ahlrichs R 1997 *Chem. Phys. Lett.* **264** 573–8

[61] Hamel S, Casida M E and Salahub D R 2001 *J. Chem. Phys.* **114** 7342–50

[62] te Velde G, Bickelhaupt F M, Baerends E J, Fonseca Guerra C, van Gisbergen S J A, Snijders J G and Ziegler T 2001 *J. Comput. Chem.* **22** 931–67

[63] Blase X and Ordejón P 2004 *Phys. Rev. B* **69** 085111

[64] Foerster D 2008 *J. Chem. Phys.* **128** 034108

[65] Foerster D and Koval P 2009 *J. Chem. Phys.* **131** 044103

[66] Blanco M A, Flórez M and Bermejo M 1997 *J. Mol. Struct.: THEOCHEM* **419** 19–27

[67] Koval P, Foerster D and Sánchez-Portal D 2014 *Phys. Rev. B* **89** 155417

[68] Ljungberg M P, Koval P, Ferrari F, Foerster D and Sánchez-Portal D 2015 *Phys. Rev. B* **92** 075422

[69] Ren X, Rinke P, Blum V, Wieferink J, Tkatchenko A, Sanfilippo A, Reuter K and Scheffler M 2012 *New J. Phys.* **14** 053020

[70] Frayssé V, Giraud L, Gratton S and Langou J 2005 *ACM Trans. Math. Softw.* **31** 228–38

[71] SIESTA www.netlib.org/blas/

[72] Haas P, Tran F and Blaha P 2009 *Phys. Rev. B* **79** 085104

[73] Yabana K and Bertsch G F 1999 *Phys. Rev. A* **60** 3809–14

[74] Rabilloud F 2014 *J. Chem. Phys.* **141** 144302

[75] Wu Z and Cohen R E 2006 *Phys. Rev. B* **73** 235116

[76] Perdew J P, Ruzsinszky A, Csonka G I, Vydrov O A, Scuseria G E, Constantin L A, Zhou X and Burke K 2008 *Phys. Rev. Lett.* **100** 136406

[77] Mattsson A E and Armiento R 2009 *Phys. Rev. B* **79** 155101

[78] Tran F, Laskowski R, Blaha P and Schwarz K 2007 *Phys. Rev. B* **75** 115131

[79] Marini A, Del Sole R and Onida G 2002 *Phys. Rev. B* **66** 115101

[80] Durante N, Fortunelli A, Broyer M and Stener M 2011 *J. Phys. Chem. C* **115** 6277–82

[81] Kuisma M, Ojanen J, Enkovaara J and Rantala T T 2010 *Phys. Rev. B* **82** 115106

[82] Weissker H C and Lopez-Lozano X 2015 *Phys. Chem. Chem. Phys.* **17** 28379–86

[83] Talman J D 2009 *Comput. Phys. Commun.* **180** 332–8

[84] Koval P and Talman J D 2010 *Comput. Phys. Commun.* **181** 2212–3

[85] Foerster D 2006 *arXiv Physics e-prints* (preprint arXiv: 0612187)

[86] Junquera J, Paz O, Sánchez-Portal D and Artacho E 2001 *Phys. Rev. B* **64** 235111

[87] Anglada E, Soler M, Junquera J and Artacho E 2002 *Phys. Rev. B* **66** 205101

[88] García-Gil S, García A, Lorente N and Ordejón P 2009 *Phys. Rev. B* **79** 075441

[89] Mark E C 1995 *Time-dependent density functional response theory for molecules Recent Advances in Density Functional Methods (Recent Advances in Computational Chemistry* vol 1) ed P C Delano (Singapore: World Scientific) pp 155–92

[90] Li Y and Ullrich C A 2015 *J. Chem. Theory Comput.* **11** 5838–52

[91] Bahn S R and Jacobsen K W 2002 *Comput. Sci. Eng.* **4** 56–66

[92] Kuo K 2002 *Struct. Chem.* **13** 221–30

[93] http://departments.icmab.es/leem/siesta/Databases/

[94] Rivero P, García-Suárez V M, Pereniguez D, Utt K, Yang Y, Bellaiche L, Park K, Ferrer J and Barraza-Lopez S 2015 *Comput. Mater. Sci.* **98** 372–89

[95] Otto A and Petri E 1976 *Solid State Commun.* **20** 823–6

[96] Lucas A A, Henrard L and Lambin P 1994 *Phys. Rev. B* **49** 2888–96

[97] Johnson P B and Christy R W 1972 *Phys. Rev. B* **6** 4370–9

[98] Babar S and Weaver J H 2015 *Appl. Opt.* **54** 477–81

[99] Tiggesbäumker J, Köller L, Meiwes-Broer K H and Liebsch A 1993 *Phys. Rev. A* **48** R1749–52

[100] Mochizuki S and Ruppin R 1993 *J. Phys.: Condens. Matter* **5** 135

[101] Haberland H 2013 *Nature* **494** E1–2

[102] Kreibig U and Zacharias P 1970 *Z. Phys.* **231** 128–43

[103] Fedrigo S, Harbich W and Buttet J 1993 *Phys. Rev. B* **47** 10706–15
[104] Harb M, Rabilloud F, Simon D, Rydlo A, Lecoultre S, Conus F, Rodrigues V and Félix C 2008 *J. Chem. Phys.* **129** 194108

[105] Conus F, Rodrigues V, Lecoultre S, Rydlo A and Félix C 2006 *J. Chem. Phys.* **125** 024511

[106] Rohrdanz M A, Martins K M and Herbert J M 2009 *J. Chem. Phys.* **130** 054112

[107] Silverstein D W and Jensen L 2010 *J. Chem. Phys.* **132** 194302

[108] Rao Y, Lei Y, Cui X, Liu Z and Chen F 2013 *J. Alloys Compd.* **565** 50–5

[109] Ewusi-Annan E 2010 Modeling the optical properties of metal nanoparticles and metal-molecule systems *PhD Thesis* The Pennsylvania State University, Lasse Jensen, Thesis Advisor Master of Science, ewusi@psu.edu

[110] Yugang S, Brian T M and Younan X 2002 *Nano Lett.* **2** 481–5

[111] Sun Y, Mayers B and Xia Y 2003 *Adv. Mater.* **15** 641–6

[112] Geddes C D (ed) 2010 *Reviews in Plasmonics 2010* (*Reviews in Plasmonics* vol 2010) (New York: Springer)

[113] Pattanayak S, Priyam A and Paik P 2013 *Dalton Trans.* **42** 10597–607

[114] Sönnichsen C 2001 Plasmons in metal nanostructures *PhD Thesis* Ludwig-Maximilians-Universität München

(https://edoc.ub.uni-muenchen.de/2367/1/Soennichsen_Carsten.pdf)

[115] Sönnichsen C, Franzl T, Wilk T, von Plessen G, Feldmann J, Wilson O and Mulvaney P 2002 *Phys. Rev. Lett.* **88** 077402

[116] Grigoryan V, Springborg M, Minassian H and Melikyan A 2013 *Comput. Theor. Chem.* **1021** 197–205 (clusters: from Dimers to Nanoparticles)

[117] White J A and Bird D M 1994 *Phys. Rev. B* **50** 4954–7

[118] Bylander D and Kleinman L 1997 *J. Comput. Phys.* **136** 599–602

[119] Balbás L C, Martins J L and Soler J M 2001 *Phys. Rev. B* **64** 165110

[120] Gaiduk A P 2013 Theory of model Kohn–Sham potentials and its applications *PhD Thesis* University of Western Ontario Electronic Thesis and Dissertation Repository paper 1099

[121] Caillie C V and Amos R D 2000 *Chem. Phys. Lett.* **317** 159–64

[122] Nazarov V U 2015 personal communications

[123] Nazarov V U and Vignale G 2011 *Phys. Rev. Lett.* **107** 216402

[124] Marques M A, Oliveira M J and Burnus T 2012 *Comput. Phys. Commun.* **183** 2272–81

[125] Delley B 1990 *J. Chem. Phys.* **92** 508–17

[126] Krack M and Köster A M 1998 *J. Chem. Phys.* **108** 3226–34