Peter Koval
email: koval.peter@gmail.com

Federico Marchesin

Dietrich Foerster

Daniel Sánchez-Portal

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

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

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.

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 clusters. The prime interest in silver clusters originates from the reduced chemical activity of silver [START_REF] Newton | Silver[END_REF], the best electron conducting 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 [START_REF] Brack | [END_REF][4][START_REF] Madison L R | Frontiers in Quantum Methods and Applications in Chemistry and Physics[END_REF][6], and in the past were also widely theoretically studied [START_REF] Kreibig | Optical Properties of Metal Clusters[END_REF][8][9][10][11]. It is experimentally established that silver nano-particles (NPs) with diameters 10 to 80 nm exhibit a surface-plasmon resonance 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 frequency [12][13][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 chemical sensing [18][19][20] and in surface-enchanced 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 [START_REF] Kreibig | Optical Properties of Metal Clusters[END_REF]. 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][25][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][10][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], densityfunctional theory (DFT) [27,[34][35][36][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.

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][START_REF] Koval | Proc. of the 19th Int. Workshop on Inelastic Ion-Surface Collisions[END_REF][START_REF] Iida | [END_REF]. 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-thanlinear 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 3 ), 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 [START_REF] Dronskowski | Computational Chemistry of Solid State Materials[END_REF]) to the realm of the optical response, as we will describe below.

Response functions formalism

The basic quantity of the linear response TDDFT is the density response function ( )

χ ω ′ r r
, , , which is a kernel of an integral operator delivering the density change ( ) δ ω r n , in response to a small external perturbation

( ) δ ω r V , ext ( ) ( ) ( ) ∫ δ ω χ ω δ ω = ′ ′ ′ r rr r n V r , , , , d . ext (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 ( ) ( ) ( )

δ δ = Ψ Ψ + r r r n t t t , , , 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

= + + r r r r V V V n V n , , eff e xt H x c
( ) ( ) ( ) ( ) on the electronic density ( ) r n via Hartree and exchange-correlation (xc) potentials, respectively. However, the so-called noninteracting response function ( )

( ) ( ) χ ω ≡ ′ δ ω δ ω ′ r r , , r r n V 0 , ,
eff 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 [START_REF] Nogueira | A Primer in Density Functional Theory[END_REF][START_REF] Gross | Basics of TDDFT[END_REF], we will repeat it here for the sake of completeness

( ) ( ) ( ) ( ) ( ) ( ) ( ) ∑ χ ω ω ε = - Ψ Ψ Ψ Ψ - - + ′ ′ ′ * * r r r r r r f f E E , , i . nm n m n m m n m n 0 (2)
Here, the occupations f n and energies E n of KS eigenstates ( ) Ψ r n do appear and ε is an infinitesimal number that accounts for the proper causality of the response. A finite value can be given to ε, 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 [START_REF] Petersilka | [END_REF][START_REF] Gross | Time-Dependent Density-Functional Formalism Benasque[END_REF] via the so-called interaction kernel ( )

( ) ( ) ≡ ′ δ δ ′ r r K , r r V n eff ∫ ″ ″ ″ ″ ″ ″ χ ω χ ω χ ω χ ω = + ′ ′ ′ ′ ′ ′ r r r r r r r r r r K r r , , , , , , , , , d d .
0 0 ( ) ( ) ( ) ( ) ( ) (3)
The TDDFT interaction kernel ′ r r K , ( ) is commonly separated into the Hartree and xc kernels [START_REF] Petersilka | [END_REF][START_REF] Gross | Time-Dependent Density-Functional Formalism Benasque[END_REF]

( ) ( ) ( ) δ δ = | -| + ′ ′ ′ r r r r r r K v n , 1 . xc (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

( ) δ ω ∝ r r V , ext
. Furthermore, the far-field response is connected to the polarisability tensor of the quantum system

( ) ( ) ∫ ω χ ω = ′ ′ ′ r r r r P r r , , d d ij i j (5)
which gives rise to the orientation-averaged optical cross section [59] ( ) ( )

∑ σ ω πω ω = c P 4 3
.

i ii opt (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 ( ) Ψ r m .

Product basis set

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

∑ Ψ = r r X f , n a a n a ( ) ( ) (7)
where the expansion coefficients X a n are determined in a diagonalisation procedure. The atomic orbitals ( ) r f a 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 . 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 3 .

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

( ) ( ) ( ) = µ µ r r r f f V F . a b ab (8)
In this equation, the complex conjugation does not appear because we use real-valued atomic orbitals [66]. The product 'vertex' coefficients µ V ab and the dominant products ( ) µ r F are determined in a diagonalisation-based procedure. Namely, we aim at identifying linear combinations of the original atomicorbital products

( ) ( ) r r f f a b ( ) ( ) ( ) = Λ µ µ r r r F f f . ab a b (9)
that are orthogonal to each other with respect to a Coulomb metric

∫ = | -| ′ ′ ′ ′ ′ ′ ′ ′ r r r r r r g f f f f r r 1 d d . ab a b a b a b , ( ) ( ) ( ) ( ) (10)
The linear combinations are built with eigenvectors of the metric

′ ′ g ab a b , λ Λ = Λ µ µ µ ′ ′ ′ ′ g ab a b a b ab , (11)
that guarantee the orthogonality requirement. Moreover, the eigenvalues of Coulomb metric λ µ are used as an indicator of the importance of a particular linear combination μ to the completeness of the basis { ( )} µ r F . Namely, the norm of dominant products ( ) µ r F (9) is proportional to the eigenvalue λ µ . Therefore, we can consistently limit the number of dominant product (9) by ignoring the eigenvector such that λ µ 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 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 ( ) µ r F 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 ( ) µ r A are sufficient in practice to express all orbital products ( ) ( )

r r f f a b
. 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 ( ) µ r F (i.e. dominant products generated for orbitals in the same atom, as opposed to bilocal pairs) as the atom-centered functions ( ) µ r A . The ansatz for atomic orbitals, analogous to equation ( 8) can be immediately written as

( ) ( ) ( ) = µ µ r r r f f P A , a b ab ( 12
)
where the atom-centered product vertex coefficients µ P ab must be still determined. In this work, we extend the ansatz (12) with a recipe for choosing these vertex coefficients µ P 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) . Notice that in order to termine ν P ab we do not need to invert the whole matrix of the Coulomb metric, but a smaller submatrix corresponding to atom-centered functions inside the contributing sphere. Thus, this step is not computationally prohibitive. However, the table µ P ab has a limited value for practical calculations. Namely, the table µ P ab can have an order of magnitude more non-zero elements that the product vertex coefficients µ V ab . This dramatic difference arises because of the distant bilocal atom pairs for which we have very few dominant products ( ) µ r F instead of many atom-centered functions ( ) µ r A contributing to such pairs. A more fruitful idea proves to be a re-expression of the bilocal dominant products ( ) µ r F in terms of atom-centered products ( ) ν r A (that are also chosen within the sphere of contributing centres)

( ) ( ) = µ ν µ ν r r F c A , ( 13
)
where the projection coefficients ν µ c can be also readily expressed as

= ν µ µν ν ν - ′ ′ c M v 1 ( ) (14)
in terms of the matrix elements of the Coulomb interaction

∫ ∫ = | -| = | -| ′ ′ ′ ′ ′ ′ µν µ ν µν µ ν r r r r r r r r M F A r r v A A r r d d , dd . ( ) ( ) ( ) ( ) (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 quanti ty 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.

Iterative method for computing the polarisability

By inserting the ansatzes [START_REF] Kreibig | Optical Properties of Metal Clusters[END_REF] and ( 8) into the non-interacting response function (2) we obtain

χ ω χ ω = ′ ′ µ µν ν r r r r F F , , , 0 0 ( ) ( ) ( ) ( ) (16)
where the tensor of non-interacting response is given by

( ) ( ) ( ) ( ) ( ) ∑ χ ω ω ε = - - - + µν µ ν f f X V X X V X E E i .
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

χ χ χ χ = + K , 0 0 (18)
where the product indices µν are dropped and the kernel matrix K is given by

( ) ( ) ( ) ( ) ( ) ∫ ∫ = | -| + ′ ′ ′ µν µ ν µ ν r r r r r r r K F F r r F K F r d d d .
xc (19) In this work, we use LDA and GGA kernels

( ) ( ) ( ) = δ δ δ r K r r E n n xc 2 xc
, which are computed in the atom-centered functions ( ) µ r A . Explicit expressions for the GGA kernel are discussed in the appendix B.

The whole response matrix ( ) χ ω µν is superfluous in the computation of the electronic polarisability [START_REF] Madison L R | Frontiers in Quantum Methods and Applications in Chemistry and Physics[END_REF]. Therefore, we further introduce the product basis set into the equation for polarisability [START_REF] Madison L R | Frontiers in Quantum Methods and Applications in Chemistry and Physics[END_REF] and using equation (18) obtain

ω δ χ ω χ ω = - µ ν µ µν ν µ µ ν ν - ′ ′ ′ ′ ′ P d K d, ij i j 0 1 0 ( ) ( ( ) ) ( ) (20)
where the dipole moments of the product functions ( )

∫ = µ µ r r d F r d i i
appear. The calculation of the polarisability ( 20) is split into a calculation of the density change in the product basis

( ) δ ω µ n [ ( ) ] ( ) ( ) δ χ ω δ ω χ ω - = K n d 0 0 ( 21
)
and a final trace with the dipole moments ( ) ( )

ω δ ω = µ µ P d n ij i j , .
The linear equation ( 21) is solved with an iterative method similar to the Arnoldi method and optimised for delivering the polarisability ( ) ω P ij with a given precision rather than the density change

( ) δ ω µ n j
, [42] that would be the target for generalpurpose iterative solvers [70].

Iterative linear equation solvers require only the action of a matrix onto given vectors to find the solution

( ) δ ω µ n j , .
In our case, the matrix reads

[ ( ) ] δ χ ω = - A K 0 .
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 χ 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 ( ) χ ω µν 0 with a vector ν z 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

χ ω χ ω = ′ ′ µ µ µ µν ν ν ν r r r r A c c A , , , 0 0 ( ) ( ) ( ) ( )
˜˜˜˜ (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

( ) χ ω µν ν z 0 .
There are several sequences of operations to organise the matrix-vector product ( ) χ ω µν ν z 0 . 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 ˜α = µ µ

V X an ab b n for occupied eigenstates states n.

The table αµ

an 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 [START_REF]SIESTA[END_REF]. The table αµ an is the major 'memory consumer' that, however, can be obviously eliminated when larger systems need to be treated. A vector µ z , on which the response function has to be applied, is converted to the basis of dominant products ˜= , 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 γ mn with the frequencyoccupation mask . 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 nearlydegenerate 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 Ag 561 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).

γ γ ω ε = - - - + - f f E E / i mn mn n m m n ˜(( ) ( ( ) ) ( ) ( ( ) )) ω ε - + - + f f E E / i n m m n

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 acc uracy 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 ω = 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).

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 [START_REF] Haas | [END_REF] 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 [START_REF] Haas | [END_REF], 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 contrib utions 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 147 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].

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 147 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)

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.

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.

Contribution of different atoms to the polarisability.

In our framework, the interacting polarisability ( ) α ω is given by

( ) ( ) () ∫ α ω δ ω = = µ µ µ µ r r d n d A dr , , (23)
where we drop all the Cartesian indices for the sake of clarity. The index μ runs over all product functions µ r A ( ), 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:

( ) ( ) α ω α ω = ∑ L L
with a partial polarisability ( ) α ω L given by

( ) ( ) ( ) α ω δ ω = µ µ ν ν d P L n . L ( 24
)
Here a padding operator ( )

µ ν
P L is introduced. The padding operator is a diagonal matrix whose elements are equal to one ( ) = µ µ P L 1 if the product index μ belongs to an atom in the layer L and zero otherwise.

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 ( ) δ ω r n , in terms of electron-hole pairs

( ) ( ) ( ) ( ) ∑ δ ω δ ω = Ψ Ψ * r r r n n , , ij ij i j ( 25
)
where ( ) Ψ r j are the KS eigenstates. This kind of expansion naturally arises in the Casida formulation of TDDFT [START_REF] Mark | Time-dependent density functional response theory for molecules Recent Advances in Density Functional Methods[END_REF], 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 [START_REF] Li | [END_REF].

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

( ) ( ) ( ) ∫ δ ω χ ω δ ω = ′ ′ ′ r rr r n V r , , , , d , 0 eff (26)
where

( ) δ ω r V ,
eff is an effective (screened) perturbation. Using sum-over-states representation (2) of the non-interacting density response function we can immediately see the electronhole expansion coefficients

( ) δ ω n ij ( ) ( ) ( ) ( ) ( ) ∫ δ ω ω ε δ ω = - - - + Ψ Ψ ′ ′ ′ ′ * r r r n f f E E V r i , d . ij i j j i j i eff (27)
Now it remains only to derive an equation for the effective perturbation

( ) δ ω r V , eff .
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

[ ( ) ( ) ( ) ] ( ) ( ) ″ ″ δ χ ω δ ω δ ω -- = ′ ′ ′ r r r r r r r r K V V , , , , , .
0 eff e xt (28) In terms of product basis (8) and for optical absorption

( ) δ ω ≡ r r V , ext
, the latter equation transforms into a linear algebraic equation

[ ( )] ( ) δ χ ω δ ω - = K V d, 0 eff (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

( ) ( )( ) ( ) ( ) δ ω ω ε δ ω = - - - + µ µ n f f X V X E E V i . ij i j a j ab b i j i eff ( 30
)
Having the density change in terms of electron-hole pairs (30) we can define a transition-resolved optical polarisability

( ) ( ) ( ) ( ) ∫ α ω δ ω = Ψ Ψ * r r r n r d ij ij j i (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 [START_REF] Dronskowski | Computational Chemistry of Solid State Materials[END_REF]. 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 (32) where the expansion coefficients

( ) ( ) ( ) ( ) ∑ δ ω δ ω = Ψ Ψ ∈ ∈ * r r r n n , , i j ij j i occ, unocc
( ) δ ω n ij read ( ) ( )( ) ( ) ( )( ) ( ) ( ) δ ω ω ε ω ε δ ω = - - - + - - + - + µ µ µ ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ n f f X V X E E f f X V X E E V i i . ij i j a i ab b j j i i j a i ab b j j i 0 (33)
Now we can define a polarisability resolved in the angular momentum l of the occupied states

( ) ( ) ∑ α ω δ ω δ = ∈ ∈ n X d X . l i j ij a j ab l l a i occ, unocc , a (34)
Obviously, the partial polarisability ( ) α ω l add up to the total interacting polarisability ( ) ( )

α ω α ω = ∑ l l
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.

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).

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 -A ˚1 using the GGA functional after Wu and Cohen [START_REF] Haas | [END_REF]78], as we discussed in section 2.5. In other words, we 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 1 4d 10 ) 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 pseudopotential (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 optim ised 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.

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) = 10 l 2 -20 l + 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 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 12 is 19.68 and 19.40 Å for the ideal and the GGA geometries, respectively, while the distance d 34 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 relaxations, we note that relaxation leads to minor geometrical distortions 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.

Optical absorption of Ag N 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 performed TDDFT calculations of the absorption spectra for compact clusters and shells. The smallest cluster is composed of two layers (Ag 13 ), while the largest cluster consists of six layers (Ag 561 ). 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 application 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 ( ) µ r A and the calculation of interaction matrices µν K takes a relatively small amount of time. For instance, in the case of our largest Ag 561 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. Photoabsorption 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 309 and Ag 561 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 147 cluster at the maxima of the lowfrequency and high-frequency bands, respectively. The direction 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 ( ) δ ω r n , were plotted for a 10% of the corresponding maximal value. One can see that the low-frequency resonance 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.

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 309 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 162 ) exhibits a low-frequency resonance at 2.8 eV, which must be compared to the 3.37 eV in the case of the Ag 309 cluster. If we leave two atomic layers as in the silver shell S5L2 (Ag 254 ), 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 309 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 309 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. )

ω (eV) Ag S5L3 Ag 309 cutout relaxed c) d)
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 S7L1 the relative intensity of the 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 943 , 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 )

ω (eV) Ag S6L3 Ag 561 cutout relaxed c) d)
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 doublelayered shells for smaller clusters, while in the case of the single-layered shells this discrepancy is visible only for the smallest shell. )

ω (eV) Ag S7L2 relaxed S7L2 S6L2 S5L2 S4L2 S3L2 S2L2 a) b) 0 0.5 1 1.5 2 2.5 3 3.5 4 1 2 3 4 6 7 8 σ opt (nm 2 ) ω (eV)
Ag S7L* relaxed

S7L4 S7L3 S7L2 S7L1 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 0 0.2 0.4 0.6 0.8 1
ω sp (eV) 1 / d (nm -1 ) S3L1 S4L1 S5L1 S6L1 S7L1 (S7L4) Clusters 2-layered 1-layered c) d)

Analysis of the optical absorption

In figure 8 we show the analysis of the atom-layer contributions 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 polarisabilities (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 contrib ution of the first 'layer', which is composed 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 resonance, the main contribution comes from the surface layer. For this mode, it is critical to describe accurately the polarisation 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 contributions, 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 contributions 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 309 cluster. We see that, on the absolute scale, the contribution 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 constructively to the absorption cross section, while the d channel contributes 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 electrons [95]. This is because the occupied states are sp-hybridised 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 trans itions 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.

Discussion

Because the clusters that we have considered are rather large, we expect that their properties will approach the properties of classical Mie spheres [START_REF] Kreibig | Optical Properties of Metal Clusters[END_REF]73] to the extent that the nonsphericity of icosahedra and charge spillage effects allow. The classical absorption cross section in the electrostatic limit depends on the dielectric function of the material ( ) ω ε and the dielectric function of the embedding medium ( ) ω ε m [START_REF] Kreibig | Optical Properties of Metal Clusters[END_REF]96] 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 309 (S5L5) cluster as a function of frequency.

-0.

( ) ( ) ( ) ( ) ( ) σ ω ω ω ω ω ω ∼ - + ε ε ε ε Im 2 . m m ( 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 ω sp from a linear dependence ω = + kd b sp 1

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][100][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 calcul ations based on Hedin's GW band structures [79]. In this paper, we showed that GGA functional by Wu and Cohen produces a strong lowfrequency 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 147 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 309 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 561 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,[START_REF] Yugang | [END_REF][111][START_REF] Geddes | Reviews in Plasmonics 2010[END_REF][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 inter ior 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]. 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 92 (S4L1 shell) we can find 3.85 (LB94) [29] and 2.8 (PBE) [28] eV, while we get 2.97 eV. For Ag 12 , Ag 42 , Ag 162 and Ag 252 (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 effectively reduced dimensionality: rods [36,37,114,[START_REF] Sönnichsen | [END_REF] and platelets [12,116] also exhibit red-shifted resonance frequencies of the plasmon excitations with respect to spherical or quasi-spherical clusters of similar sizes.

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 resonance of silver clusters depends on the size and morphology of the clusters. Namely, the frequency of the maximal absorption 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 negligible as the thickness of the shells increases. Both observations 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 realisation 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 frequency-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 extension of the number of treated atoms by an order of magnitude, and even further if full MPI parallelisation is implemented.

∫

= - - r R r r R d f f r d ab a a b b
. 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 ), but not generally for ss and pp symmetries.

( ) ( ) ( ) ( ) ∫ ∫ = + d r r r R r r f f r f f r d d ,
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 309 cluster. The contrib ution 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 309 . 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 Ag 309 s p Tot 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 ε explicitly depends on the charge density ( ) r n and its gradient In the Kohn-Sham formulation of the DFT we need the functional derivative of the energy with respect to the electron density ( ) ( )

δ δ ≡ = ∂ ∂ -∇ ∂ ∂∇ ε ε r r v E n n n . (B.2)
and in the linear response TDDFT we also need a second variational derivative of the energy with respect to the electronic density While the calculation of the KS potential in GGA is widely discussed in the literature [117][118][119], the TDDFT kernel (B.3) is less commonly discussed [START_REF] Gaiduk | Theory of model Kohn-Sham potentials and its applications[END_REF]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 between the atomic orbitals ( ) r f a . 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 ( ) σ ≡ ∇n 2 . After applying the chain rule, we transform the last formula to 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.

( ) ( ) ( )

. Exercising the same approach as for the potential in the previous section, we obtain for matrix elements of the GGA kernel

  matrix β an is a full rectangular matrix. The calculation continues with a matrix-matrix multiplication γ β

.

  The next O(N 3 ) matrixmatrix multiplication reads ˜β γ tracing over a and n indices in the product of β an with the precom-

Figure 1 .Figure 2 .

 12 Figure 1. Absorption cross section of the icosahedral Ag 147 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 0 xc starting from the WC results. Relaxation of geometries was done with PZ or WC functionals.

Figure 3 .

 3 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.

Figure 4 .

 4 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 ( ( )) δ ω r n Re , for the Ag 147 cluster at the frequencies 3.4 and 6.6 eV, respectively.

Figure 5 .

 5 Figure 5. The absorption cross sections of icosahedral silver shells constructed from the Ag 309 cluster. In panel (a) the cross section of the Ag 309 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.

Figure 6 .

 6 Figure 6. The absorption cross sections of icosahedral silver shells constructed from the Ag 561 cluster. In panel (a) the cross section of the Ag 561 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.

Figure 7 .

 7 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).

Figure 8 .

 8 Figure 8. Contribution of the different atom layers to the absorption cross section of the Ag 309 cluster (S5L5).

Figure 10 .

 10 Figure10. 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.

  matrix element d ab is zero if both terms in the last equation are zero. This is always the case for sd channels (

Figure A1 .

 A1 Figure A1. Analysis of the product-function angular-momentum contributions to the absorption cross section of the Ag 309 cluster. The total absorption cross section is also plotted for comparison.

4 )

 4 In order to 'exchange' the gradients of density variation δ the density variation δn, we represent the last equation in form of an integral with a δ-function ∫ and apply the integration by parts whenever necessary. Finally, we can get[122]

5 )

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

Figure A2 .B. 2 .

 A22 Figure A2. Analysis of atomic orbital angular momenta of the products of atomic-orbital contributions to the optical absorption cross section of the Ag 309 cluster. The total absorption cross section is also plotted for comparison.

  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 ( one more integration by parts, we can transform the expression (B.8) into a more suitable form

B. 3 .

 3 Computation of the GGA matrix elements of the kernelIn our implementation of the linear response TDDFT, we need the matrix elements of the xc kernel (B.6) defined by

10 )

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

  Namely, the vertex coefficients µ V 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 μ generates a matrix

	properties of the product vertex coefficients	µ V ab	= Λ µ ab in
	equation (8). = ∑ µ µ s V ab ab , which has
	the sparsity of the usual overlap matrix	S	ab	=	∫	f	a	( ) ( ) r r f d r b	,
	while the summation over one of the orbitals generates a rec-
	tangular sparse matrix	µ R a	= ∑	b	V	µ ab because, by construction,
	a dominant product index μ 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.
										F	( ) µ r	and also determines the sparsity

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.

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 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 ( ) µ r F

. 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 μ carries a well-defined angular momentum µ l . We can write the polarisability as a sum over the angular momenta of the product basis ( ) ( ) α ω α ω = ∑ l l with a product angular momentum-resolved polarisability ( ) α ω l given by

where the dipole moments µ d refer to a global origin of the coordinate system. The dipole moments are defined for product functions ( ) µ r F that are centered on atoms Yet another type of analysis involves the angular momentum symmetry of the atomic orbitals, the products of which expand the density change ( ) δ ω r n , . Namely, the product function ( ) µ r F is expressed in terms of a linear combination of the products of atomic orbitals (9). Therefore, we can immediately write the interacting polarisability ( ) α ω as a sum over atomic orbital angular momentum sub-sums ( ) ( )