Quantum Hydrodynamic Theory for Plasmonics: Impact of the Electron Density Tail

2 MB·

by  Cristian Ciracì

Quantum Hydrodynamic Theory for Plasmonics: Impact of the Electron Density Tail
Cristian Cirac`ı1 and Fabio Della Sala2,1
1Istituto Italiano di Tecnologia (IIT), Center for Biomolecular
Nanotechnologies@UNILE, Via Barsanti, 73010 Arnesano, Italy.∗
2Istituto Nanoscienze-CNR, Euromediterranean Center for Nanomaterial
Modelling and Technology, Via per Arnesano 16, 73100 Lecce, Italy.
(Dated: May 10, 2016)
Multiscale plasmonic systems (e.g. extended metallic nanostructures with sub-nanometer inter-
distances)playakeyroleinthedevelopmentofnext-generationnano-photonicdevices. Anaccurate
6 modelingoftheopticalinteractionsinthesesystemsrequiresanaccuratedescriptionofbothquan-
1
tum effects and far-field properties. Classical electromagnetism can only describe the latter, while
0
Time-Dependent Density Functional Theory (TD-DFT) can provide a full first-principles quan-
2
tum treatment. However, TD-DFT becomes computationally prohibitive for sizes that exceed few

y nanometers, which are instead very important for most applications. In this article, we introduce
a amethodbasedonthequantumhydrodynamictheory(QHT)thatincludesnonlocalcontributions
M of the kinetic energy and the correct asymptotic description of the electron density. We show that
our QHT method can predict both plasmon energy and spill-out effects in metal nanoparticles in

9 excellent agreement with TD-DFT predictions, thus allowing reliable and efficient calculations of
both quantum and far-field properties in multiscale plasmonic systems.

]
l
l
a I. INTRODUCTION Recently, TD-DFT has been applied to metallic wires
h
withadiameterupto20nm25,26 andtometallicspheres
-
s (and sphere dimers) with around one thousand valence
Plasmonic systems have received a renewed great deal
e electrons27–31. For larger systems TD-DFT rapidly be-
m of attention for their ability to localize electromagnetic
comes computationally prohibitive, as all single particle
radiation at visible frequencies well below the diffrac-
t. tion limit and enhance the local electric fields hun- orbitals need to be computed.
a
m dred times the incident radiation1–3. These properties An alternative approach is to use a simple lin-
make plasmonic structures valuable candidates for en- earized Thomas-Fermi hydrodynamic theory (TF-HT),
-
d hancing nonlinear optical phenomena4, controlling sur- also known simply as the hydrodynamic model, which
n face reflectance properties5, enhancing the far-field cou- takes into account the nonlocal behavior of the electron
o pling with nanometer-sized elements6–8, such as quan- response by including the electron pressure32–35. The
c tumemitters,andstudyingfundamentalphenomena9–11. introduction of an electron pressure term in the free-
[
In particular, nanostructures supporting gap-plasmon electron model accounts for the Pauli exclusion princi-

2 modes12 constitute an important platform. Advances in ple within the limit of the Thomas-Fermi (TF) theory36.
v nano-fabricationtechniques13,14 havemadeitpossibleto In contradistinction to the treatment of the electron
4 achieve separation between two metallic elements, i.e. response in classical electromagnetism, where induced
8
particlesornanowires, ofonlyafractionofananometer. charges are crushed into an infinitesimally thin layer at
5
At such distances nonlocal or quantum effects becomes the surface of the metal, the induced electron density in
1
0 non-negligible. It has been shown that the resonance of theTF-HTapproachratherspreadsoutfromthesurface
. nanoparticle dimers or film-coupled nanospheres can be into the bulk region34. In fact, the TF-HT method is
1
perturbed by such effects9,11,15. Form the theory stand usuallycombinedwiththeassumptionthattheelectrons
0
6 point, it may be really challenging, if not impossible, to cannot escape the metal boundaries (hard-wall bound-
1 accurately describe at once the entire multiscale physics ary conditions). The advantage of TF-HT with respect
: involved in such systems. On the one hand, one has a tofullquantummethodsisthatitcanbeeasilyemployed
v
macroscopic electromagnetic system constituted by the forstructuresoftheorderofseveralhundrednanometers
i
X whole plasmonic structure, on the other hand, as the in size.
r gap closes it is crucial to take into account the quantum The TF-HT model dates back to the ’70s37–39 and
a
nature of the electrons in the metal16,17. closed-form analytic solutions exist for a homogeneous
A time-dependent density functional theory (TD- sphere40–42. Nonetheless, its applicability in complex
DFT) approach allows the exact calculation of plasmon, plasmonics has been limited by the absence of experi-
aswellassingle-particle,excitationenergiesinbothfinite mentalconfirmationorthevalidationbyhigherlevelthe-
and extended systems18. It can be implemented both in oreticalmethods,whichareneededtoverifytheassump-
real-time propagation19,20, and frequency domain linear- tions and approximations used in constructing solutions.
response21, and can serve as reference for developing ap- In the last decade, however, the improvement of fabrica-
proximate schemes. In the context of plasmonics, TD- tion techniques and the proliferation of self-assembling
DFT has been largely applied to nanosystems with an colloidal plasmonic structures have provided a robust
atomisticdescription22oremployingajelliummodel23,24. platform43,44 for studying extremely sub-wavelength op-
2
tical phenomena, thus reinvigorating the interest in TF- hydrodynamicquantities77,78: theelectrondensityn(r,t)
HT45–53. In particular, Cirac`ı et al. applied the TF- andtheelectronvelocityfieldv(r,t). Undertheinfluence
HT method to plasmonic nanostructures consisting of of the electromagnetic fields E and B the electronic sys-
film-couplednanoparticlesandfoundthatthemodelcan tem can be described by the equation35,79:
provide predictions that are both in qualitative and in
quantitative agreement with experiments9. More re- (cid:18) (cid:19)
∂ δG[n]
cently,TF-HTresultshavebeencomparedwithTD-DFT m +v·∇+γ v=−e(E+v×B)−∇ ,
calculations25,54, showing some limitations. In fact, es- e ∂t δn
(1)
sentialeffectssuchaselectronspill-outandquantumtun-
withm ande,theelectronmassandtheelectroncharge
nelingarecompletelyneglected. Inordertoincludesuch e
(inabsolutevalue)respectively,andγthephenomenolog-
effects, other methods based on effective descriptions
have also been proposed26,27,51,55,56, as well as methods ical damping rate. The energy functional G[n] contains
based on the real-time orbital-free TD-DFT57,58. the sum of the interacting kinetic energy (T) and the
exchange-correlation (XC) potential energy (U ) of the
To include spill-out effects in TF-HT the first step is xc
electronic system. In DFT the XC potential energy is
to consider the spatial dependence of the electron den-
definedasU =E −(T−T )36,80 whereE istheXC
sity. This scheme traces back to the ’70s both for fi- xc xc s xc
nite systems (e.g. atoms)59–61 and surfaces38,62,63, and energyandTs isthenon-interactingkineticenergy: thus
we have G[n]=T[n]+U [n]=T [n]+E [n]. In this
hasbeenrecentlyreconsideredusingequilibriumelectron XC s XC
density from DFT calculations64. When spatial depen- work we employ the following approximation for G[n]:
dence of the electron density is included, however, one
should consider nonlocal contributions, namely the von (cid:18) 1 (cid:19)
Weizs¨acker term, to the free-electron gas kinetic energy, G[n]≈Gη[n]= TsTF[n]+ ηTsW[n] +EXLDCA[n], (2)
in place of the simple TF kinetic energy. This approach
isusuallynamedquantumhydrodynamictheory(QHT),
and it has been widely used photoabsorption of atoms65, where TTF[n] is the kinetic energy functional in the TF
s
metallicnanoparticles66,67,plasmaphysics68–71,andtwo- approximation, TW[n] is the von Weizs¨acker kinetic en-
s
dimensional magnetoplasmonics72,73. Very recently, the ergy functional36 and ELDA[n] is the local density ap-
XC
QHT method has been systematically investigated for proximation (LDA) for the XC energy functional. The
surfaces74 and a self-consistent version of QHT has been expressions for these functionals can be easily found in
presented and applied to plasmonic systems75,76. How- the literature (see e.g. Ref. 81); here we report the ex-
ever,theimpactoftheelectronicground-statedensityon pressionfortheirpotentialsobtainedbytakingthefunc-
the QHT optical response is yet unclear. tional derivative with respect to n:
In this paper, we will first study the influence of
ground-state electron density profile on the linear re-
sponse of metallic nanospheres described by using the δTsTF =(E a2)5c n2/3, (3a)
QHT method and compare our results with reference δn h 0 3 TF
TD-DFTcalculations. WefindthatQHTcanaccurately δTW 1(cid:18)∇n·∇n ∇2n(cid:19)
s =(E a2) −2 (3b)
describe both the plasmon resonance and the spill-out δn h 0 8 n2 n
effects only when it is combined with the exact DFT
δELDA (cid:18) 4 (cid:19)
ground-state density. XC =(E ) −a c n1/3+µ [n] =v (r) (3c)
Secondly, we will show that by using an analytical δn h 03 X C xc
model for the ground-state electronic density, it is possi-
bletoreproduceTD-DFTresultsandtoincluderetarda-
tion effects simultaneously, thus allowing the calculation where Eh = m(cid:126)e2a20 is the Hartree energy, a0 is the Bohr
of plasmonic systems exceeding hundred nanometers. radius, c = 3 (3π2)2/3, and c = 3(3)1/3. Eqs. (3),
TF 10 X 4 π
as well as other formulas in this paper are in S.I. units;
expressions in atomic units (a.u.) can also be easily ob-
II. QUANTUM HYDRODYNAMIC THEORY tained by considering that E = a = m = (cid:126) = 1. The
h 0 e
term in Eq. (3c) is the XC potential v (r); the corre-
xc
Withinthehydrodynamicmodel,themany-bodyelec- lation potential µ [n] (in atomic units) is obtained from
C
trondynamicofanelectronicsystemisdescribedbytwo the Perdew-Zunger LDA parametrization82:
(cid:40)ln(r )(cid:0)a+ 2cr (cid:1)+(cid:0)b− 1a(cid:1)+ 1(2d−c)r , r <1
µC[n]= α+((s761α+ββ11)√√rrs3s++β(2s43rsα)β22)rs, 3 3 s rss ≥1
3
witha r =( 3 )1/3 beingtheWigner-Seitzradius. The whereε isthevacuumpermittivity,cthespeedoflight,
0 s 4πn 0
coefficients are a = 0.0311, b = −0.048, c = 0.002, d = n (r)istheunperturbed(ground-state)electrondensity,
0
−0.0116, α=−0.1423, β1 =1.0529 and β2 =0.3334.82 and ωp(r)=(cid:112)e2n0(r)/(meε0) is the spatially dependent
ThekeyparameterinEq. (2)isη,whichisintherange plasmafrequency. Thefirstordertermsforthepotential
[1,∞] (usually in the literature λ = 1/η is used). While can be calculated as:
the TF approximation (η =∞) is exact only in the bulk
region where the electron density becomes uniform, the
von Weizs¨acker term adds a correction that depends on
the gradient (i.e. on the wavevector k in the reciprocal
space). In general choosing the parameter η = 9 gives a
goodapproximationforaslowlyvaryingelectrondensity
(kk.83(cid:28)In1)t,hwishwiloertkakwiengwηill=co1nsgiidveersbeoxtahctηr=esu1ltasnfdorηl=arg9e. (cid:18)δδGnη(cid:19) =(cid:90) δnδ(rG)ηδ[nn(]r(cid:48))(cid:12)(cid:12)(cid:12)(cid:12) n1(r(cid:48))dr(cid:48) (5)
The latter has been used in Ref. 75. 1 n=n0(r)
Equation(1)hastobecoupledtoMaxwell’sequations.
By linearizing the system with the usual perturbation
approach84, taking into account the continuity equation,
and the fact that ∂P/∂t = J = −nev, we obtain in the
frequency domain the following system of equations:
∇×∇×E− ωc22E=ω2µ0P, (4a) wpeitrhtunrb1a=tion1e∇. C·lPearblyeiEng,Pthaenedlenc1traornecdoemnspilteyxfiqrusatnotritdieesr
(cid:18) (cid:19) and depend on ω. Using the expressions (3) and Eq. (5)
en0∇ δGη +(cid:0)ω2+iγω(cid:1)P=−ε ω2E, (4b) the first order terms for the potentials are:
m δn 0 p
e 1
(cid:18)δTTF(cid:19) 10
s =(E a2) c n −1/3n (6a)
δn h 0 9 TF 0 1
1
(cid:18)δTW(cid:19) 1(cid:34)∇n ·∇n ∇2n |∇n |2 ∇2n (cid:35)
s =(E a2) 0 1 + 0n − 0 n − 1 , (6b)
δn h 0 4 n2 n2 1 n3 1 n
1 0 0 0 0
(cid:18)δELDA(cid:19) (cid:18) 4 (cid:19)
XC =(E ) −a c n−2/3n +a3µ(cid:48) [n ]n , (6c)
δn h 09 X 0 1 0 C 0 1
1
with
(cid:40)−4π (cid:2)a+ 1(1+2ln(r ))cr + 2dr (cid:3)r3, r <1
µ(cid:48)C[n]= α27π95β1√rs+3(7β12(1++8ββ12√)rrss++s2β12βr1sβ)s23rs3/23+16sβ22rs2srs3, rss ≥1
Inthefollowingthesystemofequations(4)willbenamed solution of the system of Eqs. (4) within a commercially
QHTη, i.e. QHT1 if η =1 or QHT9 if η =9. available software based on the finite-element method,
Comsol Multiphysics85 (see Appendix A). In particu-
It is useful to notice that Eq. (4b) reduces to the TF-
lar, we have implemented the method using the 2.5D
HT model if the XC and the von Weizs¨acker functionals
technique86, which allows to easily compute absorption
areneglected,andassumingthattheequilibriumelectron
spectra for spheres or, more generally, axis symmetric
density is uniform in space, n (r) ≡ n . In this case, in
0 0
structuresoftheorderoffewhundrednanometersinsize.
fact, the first term on the right-hand side of Eq. (4b)
becomes β2∇(∇·P) with β2 = 10cTFn2/3 = v2/334,54.
9 me 0 F
As already pointed-out in the introduction, TF-HT will
alwaysbeassociatedwithhard-wallboundaryconditions,
i.e., P·nˆ =0, atthemetalboundaries, withnˆ beingthe III. JELLIUM NANOSPHERES
unit vector normal to the surface.
The equations QHTη can be solved with a plane-wave The results of the QHT approach will directly depend
excitation for a range of frequencies ω; the solution ontheinputgroundstateelectronicdensityn (r),which
0
vectors E and P can then be used to compute the linear defines the system under consideration. This is very dif-
optical properties. We have implemented a numerical ferent from classical plasmonics, where the system is de-
4
fined by its local dielectric constant. TF-HT KS/QHT1
Ideal systems to test the QHT approach are repre-
sented by jellium nanospheres23, where Ne electrons are 1 a) 1 e)
confinedbytheelectrostaticpotentialgeneratedbyauni- n+
n
ftoivremlcyhacrhgaergdeednssipthyerne+of=ra(dri3u4sπR/3)=−1rsiNnse1i/d3e,waitnhdpzoesrio- nsity0.5 0 nsity 0.5
o2rsutt=osid46ea;a.hu.ue.,.rewirhnsicihrseratehlpemrWeesteiagnlntsse.rss-oSIdneiiuttzmhir.saIdwniuoosrr,kdreawrnetgoicnoegxnasfricdotemlyr Charge de1.015 RIme((nn11)) b) Charge de 00..024 f)
0.5
include all quantum effects, n (r) should be the exact
0 0
quantum-mechanicaldensityofthesystemunderconsid- 0
eration, obtained for example from a full ground-state 0 20 40 0 20 40
Distance (bohr) Distance (bohr)
Kohn-Sham (KS) DFT calculation. We have developed
an in-house code for the self-consistent solution of the c) g)
KS equations for jellium nanospheres, with the LDA XC
functional. Calculations are performed with finite differ-
ences on a linear numerical grid up to R+50 bohr. The n n
final ground-state electron density can be written as: 1 1
L(cid:88)max(cid:88)nl 2l+1 d) h)
nKS(r)= f R (r)2 (7)
0 4π nl
l=0 n=0
where R (r) is the solution of the radial Schr¨odinger
nl
equation and f = 2 (we consider only the spin- |E| |E|
restricted case). The electronic configuration of a jel-
lium nanosphere is characterized by L and a se-
max
quence of shell-number S = [n ,n ,...,n ] with FIG.1. Jelliumnanosphere(r =4)withN =508electrons
0 1 Lmax s e
n ≥ n ... ≥ n , i.e. there are n occupied or- (R= 31.9 a.u.) as obtained from the TF-HT (panel a,b,c,d)
0 1 Lmax l
bitals with angular momentum l. The total number of andtheKS/QHT1(panelse,f,g,h)approaches. Ground-state
electrons is then N = (cid:80)Lmaxn f(2l + 1), which can density n0(r), panel a) and e); Real and imaginary part of
e l=0 l the induced charge density n (r) at the plasmon resonance,
be called shell-closing numbers. The so called “magic- 1
panel b) and f); Imaginary part of n (r) at the cross-section
number” clusters not only have a shell-closing but also a 1
plane, panel c) and g); Norm of the induced electric field at
(large) positive KS energy-gap87–90.
the cross-section plane, panel d) and h).
Wehaveimplementedacodethatcomputesallmagic-
number jellium nanospheres up to an arbitrary number
of electrons. Starting from S = [1], i.e. a system where
there are only N = 2 electrons in the lowest 1s-shell,
e In Fig. 1 we compare the results of the TF-HT ap-
the program tries to fill other shells in order to keep the
proach for a jellium nanosphere with N = 508 elec-
e
KS energy-gap as large as possible. In the first step the
tronswiththesolutionoftheQHTequationswithη =1
program thus compares the KS energy-gap between jel-
(QHT1) using the exact KS ground-state electronic den-
lium nanospheres with S = [2] and S = [1,1], and ob-
sity; this approach will be referred to as KS/QHT1.
viously it finds that the latter is the next magic-number
While the TF-HT approach assumes that n (r) = n+,
0
jellium nanosphere. Then the same procedure is applied
see Fig. 1a, the KS ground-state density spreads out
toS =[1,1],comparingS =[2,1]andS =[1,1,1]andso
from the jellium boundary, see Fig. 1e). The resulting
on. Inthiswaythefirstjelliumnanospheresobtainedare
induceddensityn (r)fromtheTF-HTmodelisconfined
1
N = 2,8,18,20,34,40,58,68,90,92,106,132.., which
e insidethejelliumboundary,seeFig. 1b)andc),whereas
are well established in the literature91. However, jel-
there is a significant spill-out in the KS/QHT1 method,
liumnanosphereswithN =68,90,106,...arecharacter-
e see Fig. 1 g) and h), as recently discussed in Ref. 75.
ized by a negative KS energy-gap (i.e. there are occupa-
This difference will lead to a different description of the
tionholesbelowthehighest-occupiedmolecularorbitals):
electric field at the surface, which is the key quantity
this means that the so obtained electronic density is not
for plasmonic applications, such as enhancement of the
a ground-state density and these clusters are not magic-
spontaneous emission rates6, sensing14,92, and nonlinear
numberclusters. Moreover,wenotethatwhenN isvery
e optical effects93,94.
large, the electronic configuration cannot be established
using simple models89–91, due to the almost degeneracy In the following of the paper, we aim to verify if
ofthehigh-lyingKSorbitals. Inthiswork,weconsidered KS/QHT1yieldscorrectn (r)andplasmonenergies,and
1
all shell-closing magic-number jellium nanospheres up to to investigate alternative paths to compute the ground-
N =5032 (see Table S1 in the Supplemental Material). state density.
e
5
IV. INPUT GROUND-STATE DENSITIES V. ASYMPTOTIC ANALYSIS
The computation of the KS ground-state density is In the KS or OF approach, if we assume that vs(r)=
out-of-reach for all but the smallest systems (computa- vxc(r)−eφ(r)goesexponentiallytozero(thisisthecase
tional cost scales as O(N3)). An alternative, computa- for a neutral system and using LDA for the XC func-
e
tionally cheaper but less accurate, is to use OF-DFT to tional), then the density asymptotically decays as36
compute the ground-state density (nOF(r)). In this case
0 A
we have to solve the Euler equation36: n (r)→ exp(−κr). (12)
0 r2
δT
δn0(sr) +vxc(r)−eφ(r)=µOF (8) I(n1/tηhge)TOsWFawpephroaavceh:,iftheKEisapproximatedasTsTF+
(cid:18) (cid:19)
whereµOF isaconstantrepresentingthechemicalpoten- 1 (cid:112) √
κOF = √ 2 −2µOF η . (13)
tialandφ(r)isthetotal(i.e. frombothelectronsandthe a E g
0 h
bare positive background) electrostatic potential. The
Eulerequation(8)canberecastintoaneigenvalueequa- In the KS approach we have36
tion for the square root of the electron density95; if the
(cid:18) (cid:19)
kineticenergy(KE)isapproximatedasTsTF+(1/ηg)TsW κKS = √1 2(cid:112)−2(cid:15)HOMO (14)
it takes the form: a0 Eh
(cid:18) 1 (cid:126)2∇2 δTTF (cid:19)(cid:112) and (cid:15)HOMO is the eigenvalue of the highest occupied
η 2m + δns(r) +vxc(r)−eφ(r) n0(r) molecular orbital (HOMO). Note that (cid:15)HOMO < 0 for
g e 0 stable electronic systems and it coincides with the neg-
(cid:112)
=µOF n0(r) (9) ative of the ionization potential only for the exact XC-
functional98.
which we solved as a self-consistent KS equation (see The values of µOF1 (i.e., OF-DFT with η = 1), µOF9
above), considering only the lowest eigenvalue (with an- (i.e., OF-DFT with η = 9) and (cid:15)HOMO for all the jel-
gular momentum l = 0). The self-consistent calculation lium nanospheres considered are reported in Fig. 2. It
of nOF(r) can be, in principle, obtained for spherical is found that |µOF9| is only a factor 1.1-1.4 smaller than
0
nanoparticles of any size (computational cost is O(Ne)), µOF1. Thus,unlessηg =1,wehavethatthedensitycom-
even if we experienced very slow convergence, especially puted in OF-DFT is decaying faster than the exact one,
for η =9. as numerically shown in Fig. 3 for a jellium nanosphere
Athirdapproachistouseamodelexpressionthatap- with N =338 electrons. This is consistent with the fact
e
proximates the exact density. For a sphere, the approxi- that the von Weizs¨acker KE approximation is exact in
mated unperturbed electron density can be described by the asymptotic region36,99.
using the model66,67,96,97: We remark that Eq. (14) is valid only in the asymp-
totic region, i.e. where the density is dominated only by
f theHOMO.However, inthecaseofjelliumnanospheres,
nMod(r)= 0 (10)
0 1+exp(κMod(r−R)) there are several KS orbitals with energies very close to
the HOMO, so that the asymptotic limit will be reached
wherer isthedistancefromthecenterofthesphereand only very far from the jellium boundary, in a region that
R is the radius of the nanosphere. The expression (10) isnotrelevantfortotalenergies,norfortheopticalprop-
has to be normalized such that the total charge equals erties(seeFig. 4). Ifinthe“near”asymptoticregion(i.e.
the total number of electron: within the simulation domain) we assume that the den-
sity decays as in Eq. (12), with κ = κeff then we can
(cid:90) +∞ 4 define an effective energy:
4π nMod(r)r2dr = πR3n+ =N . (11)
0 3 e
0 (κeff)2
µeff =−(E a2) . (15)
h 0 8
This approach, if successful, is particularly useful to
compute the spectral response of arbitrary big systems, The values of µeff are also reported in Fig. 2, and they
since it provides the ground-state density without any areclearlylarger(inabsolutevalue)thanthe(cid:15)HOMO;the
computational cost. We underline that Eq. (10) is not difference increases with the number of electrons, due to
employed for a variational calculation of the ground- the increasing contribution of other (low-lying) orbitals.
state density66,67,96,97. Instead we will fix κMod, which Foraninfinitenumberofelectrons,alinearextrapolation
describes the asymptotic decay of the electronic density givesµeff,∞ ≈−3.75eV.Wethenusethisvaluetodefine
and is the only parameter in Eq. (10), as described in
(cid:18) (cid:19)
the next section. 1 (cid:112)
κMod = √ −8µeff,∞ ≈1.05. (16)
a E
0 h
6
N R R
0 e 100 sim
-2 532 338 92 34 8 a.u.)10-20 a)
y (
sit10-40
n
-2.5 De
eV) µOF9 10-60
gy ( -3 µOF1 V)3.5 b)
Ener-3.5 εµHeOffMO nergy (e 3 εHOMO
E
2.5
eff,∞
ε
-4 0 50 100 150
0 0.5 1 1.5 Distance (bohr)
-1
1/D (nm )
FIG. 4. (a) KS Ground-state electronic density n (r) for
0
FIG. 2. Eigenvalues versus the inverse of the jellium a jellium sphere with Ne = 338 electrons (R = 27.86
nanosphere(r =4a.u.) diameter;chemicalpotential(µOF9) bohr, indicated by the solid-blue vertical line) ; (b) Plot of
of orbital-freesDFT calculations with η=9 (purple squares); −(E a2)1(cid:16)dln(r2n0(r))(cid:17)2, which represents the ’local’ µeff in
chemical potential (µOF1) of orbital-free DFT calculations h 0 8 dr
the case of real density. Only in the far asymptotic region
with η = 1 (orange circles); HOMO eigenvalues ((cid:15)HOMO) (e.g. forr>150bohr)µeff approaches(cid:15)HOMO indicatedbya
ofKS-DFTcalculations(bluediamonds);effectiveeigenvalue
horizontalgreendashedline. Inthesimulationdomain, indi-
(µeff) from the KS-DFT electronic density decay (green tri- cated by the vertical red-line, (cid:15)eff is significantly larger than
angles), see text for details. (cid:15)HOMO.
0
10
KS hereby limiting our investigation to dipolar excitations.
OF1
-3 OF9 To proceed, we take the divergence of Eq. (4b), and
10
Model we use the quasistatic approximation (so that ε ∇·E=
0
u.) 10-6 ∇·P=en1), obtaining:
a. 5×10-3
y (
Densit1100-1-92 Density (a.u.)43××1100--33 ∇·emne0∇(cid:18)δδGn(cid:19)1+(cid:0)ω2(cid:1)en1 =−me2e (cid:18)εe0n0n1+∇n0·(1E8(cid:19)) .
2×10-30 10 20 30 In Eq. (18) we also assume no damping (i.e. γ =0) and
10-15 Distance (bohr) noexternalfield(i.e. weareconsideringonlyfreeoscilla-
tions). TheasymptoticsolutionofEq. (18)canbeeasily
0 10 20 30 40 50
Distance (bohr) found considering that the second term on the left-hand
side is proportional to n : thus all terms which decay
1
FIG. 3. Ground-state electronic density (in a log-scale) for exponentially faster than n can be neglected. These
1
a jellium nanosphere (rs = 4) with 338 electrons (R=27.8 are: the TF and XC contributions in the first term on
bohr)computedbyKS-DFT,OF-DFTwithηg =1(OF1)and the left-hand side, which are proportional to n2/3n and
with η = 9 (OF9), and the model density. The inset shows 0 1
g n1/3n , respectively (see Eq. (3a) and (3c)), and the
theground-stateelectronicdensityinalinearscaleinsidethe 0 1
first term to the right-hand side (proportional to n n ).
nanosphere. 0 1
The second term on the right-hand side requires spe-
cial attention. Asymptotically it decays proportionally
to (n d)/r3, where d is the dipole moment of n . Thus
Figure 3 shows that very good agreement is obtained in 0 1
Eq. (18) has an asymptotic solution only and only if n
the asymptotic region, between the model and the KS 0
decays faster than n , i.e. if
density. 1
We now move to consider the asymptotic solution of
β <κ. (19)
theQHTη equationsforsphericalsystems,extendingthe
work in Ref. 74, where only slabs have been considered,
and the early one in Ref. 65. If we assume the ground- Using Eqs. (12) and (17) in Eq. (18), we obtain (after
state density decay in Eq. (12) then we want to verify if some algebra) that Eq. (18) is asymptotically satisfied if
Eq. (4) has solutions of the type
(cid:18)E a2(cid:19) 1 (cid:18)κ2β2 β4 κβ3(cid:19)
h 0 + − =ω2 . (20)
n1(r)→Bexp(−βr)cos(θ) (17) me η 4 4 2
7
1.5
VI. PLASMON RESONANCE AND SPILL-OUT
EFFECTS.
β
4
1
β
3 In Fig. 6 we plot the absorption cross-section, σ, nor-
malized to the geometrical area σ = πR2, for a Na
0
κ 0.5 jellium sphere (r = 4 a.u.) with N = 338 and thus
β/ s e
β R = 27.86 a.u. (D = 2.94 nm), using different ap-
2
proaches.
0 Figure6a)reportsthereferenceTD-DFTresults. TD-
β1 DFT calculations (in the adiabatic LDA) have been
performed using an in-house developed code, following
-0.5
0 0.5 1 1.5 2 the literature24,102–104. Details of our TD-DFT numer-
ω/ω ical implementation, which allows calculations for large
c
nanospheres will be discussed elsewhere. In TD-DFT
FIG.5. GraphicalrepresentationofthesolutionsinEq. (21), (where no retardation is included) the absorption cross-
see text for details. section can be computed as:
ω
σ(ω)= Im[α (ω)] , (23)
This equation has four solutions of the type: cε zz
0
(cid:114)
κ κ ω where the frequency dependent polarizability is:
β = ± 1± (21)
2 2 ωc (cid:90)
α (ω)=−e2 drdr(cid:48)zχ(r,r(cid:48),ω)z(cid:48) , (24)
zz
where the critical energy is:
κ2(cid:115)(cid:18)E a2(cid:19) 1 η and χ(r,r(cid:48),ω) = δn(r)/δ(eVext(r(cid:48))) is the interacting
(cid:126)ω =(cid:126) h 0 =|µ|√g (22) density response function18.
c 8 m η η
e Panels b), c) and d) of Fig. 6 report the QHT absorp-
tion cross section computed as:
where we used Eq. (13). Note that it turned out that
these solutions are identical to the slab case74. 1 ω (cid:90)
The four solutions are shown in Fig. 5. Solution β σ(ω)= Im{E·P∗}dV (25)
1 I 2
is negative, i.e. it is asymptotically increasing, thus it 0
is excluded by the boundary conditions. Solution β4 is where Io is the energy flux of the incident plane-wave.
excluded by the condition in Eq. (19). Solution β2 and Panel6b)showsthespectrumobtainedbyapplyingthe
β3 arerealonlyforω <ωc. Forω >ωc,β2,3 arecomplex QHT method with η = 9 (QHT9) to the OF9 density;
withtherealpartfixedtoκ/2. Theaboveresultsarecon- this approach will be called self-consistent OF9/QHT9
sistent with the TD-DFT calculations of finite systems andcoincideswiththeapproachofToscanoet al.75,with
(η = ηg = 1), where (cid:126)ωc = (cid:15)HOMO can be interpreted theonlydifferencebeingthechoiceoftheXCfunctional.
as the ionization threshold100. In fact, in TD-DFT the Theenergypositionofthefirstpeak(≈3.2eV)isinvery
computation of excitation energies higher than (cid:15)HOMO good agreement with the TD-DFT result (≈ 3.15eV).
(i.e. the plasmon peak, too) can be challenging because Obviously TD-DFT results are much broadened due to
all the continuum of virtual orbitals must be accurately quantum size effects as N is quite small. However, the
e
described. Inthesamewaythespectracalculatedwithin decay of n is very different from the reference TD-DFT
1
QHT are well convergent up to energy (cid:126)ωc. results. In fact from Eq. (22) we obtain that (cid:126)ωc =
When the sphereis excited by photons with an energy 3|µOF9|. From Fig. 2 we see that µOF9 ≈ −2.4 eV and
larger than (cid:126)ωc, we experienced a large dependence thusthecritialfrequencyisartificiallymovedtoveryhigh
on the domain size. This is due to the fact that the energy ((cid:126)ω ≈7.2 eV). A good point of the OF9/QHT9
c
induced charge density acquires a propagating charac- approach is that the computation of all the spectrum
teristic typical of electrons in vacuum. The boundary (i.e. up to 7.2 eV) will be numerically stable. On the
condition we used (P = 0) is no longer valid since other hand, from Eqs. (21) and (13) we see that the
it produces an artificial scattering of the electrons at first solution will have a decay βOF9/QHT9 ≈0.87κOF9 ≈
the simulation boundary and appropriate boundary 2.61κKS,i.e. muchmoreconfinedthanthereferenceTD-
conditions should be developed65,101. For the jellium DFTresults,asnumericallyshowninFig. 7. Recallthat
nanospheres considered in this work we have that µeff ≈ in TD-DFT the induced density will also decay as in Eq.
3.5 eV (see Fig. 2) which is a bit above the Mie energy (17) with κKS/2≤β ≤κKS as discussed in Ref. 26. The
(cid:126)ω = (E )(cid:112)1/r3 =3.4 eV. Numerically we found so-called spill-out effects in computational plasmonics,
Mie h s
that only the calculation of the main (first) plasmon which indeed refer to the profile of induced density, are
peak is stable. thuslargelyunderestimatedintheOF9/QHT9approach.
Thus the good accuracy of the OF9/QHT9 resonance
8
energyseemsoriginatingfromerrorcancellationbetween 2
a) TD-DFT
between the too confined ground-state electron density
(from OF9) and the approximated kinetic energy-kernel 0
σ 1
(QHT9, which is valid only for slowly varying density). /
σ
In Fig. 6c) we report the results from the KS/QHT1
approach, already introduced in Fig. 1. In this case the 0
resonance peak (≈ 3.13 eV) is in even better agreement b) OF9/QHT9
4
with the TD-DFT results. Almost the same results are 0
σ
obtained by applying the QHT1 method to the model
σ/ 2
ground-statedensity(Mod/QHT1),asshowninFig. 6d).
More importantly Fig. 7 shows that the induced density
0
n from KS/QHT1 has almost the same decay of the
1 c) KS/QHT1
TD-DFT result, i.e. the KS/QHT1 approach correctly 4
0
describes the spill-out effects of the induced charge den- σ
sity. σ/ 2
It is useful to remark that, despite the fact the all
0
the spectra presented in Fig. 6b-d) result quite simi-
d) Mod/QHT1
lar, QHT is very sensible to the density tail. Using a 4
model density with a larger (smaller) κMod yields a red- 0
σ
shifted (blue-shifted) plasmon peak. In a similar way, σ/ 2
usingtheQHT1methodwiththeOF1ground-stateden-
sity yields a plasmon peak red-shifted by about 0.3 eV 0
(data not reported). This is not surprising considering 2 2.5 3 3.5 4 4.5 5
thattheOF1densityisdecayingmuchmoreslowlythan Energy (eV)
the effective one, see Fig. 2. As mentioned at the end of
Section V, spectral features appearing at energies higher FIG. 6. Absorption cross section (σ) normalized to the geo-
than the ionization energy ((cid:126)ω ) are not stable and will metricalarea(σ )forajelliumnanosphere(r =4a.u.) with
c 0 s
be investigated elsewhere. We also point out that QHT Ne = 338 electrons as obtained form TD-DFT, OF9/QHT9,
withη =9(whichyieldstheexactdielectricresponsefor KS/QHT1 and Mod/QHT1. All the spectra have been ob-
smallwavevectors71)cannotbeusedincombinationwith tained using an empirical broadening of 0.066 eV.
a ground-state density with the exact asymptotic decay.
This seems surprising, but it can be easily justified by
tl(cid:126)ohωorckeei=ntgim|µaeteffsE|s/qm3.a≈l(le2r12.)t1.h6aeInVf ,tηhgie.e=.Mi1teheafrnecdqruitηeicn=acly,9frseowqeutheonabtctytahiines y (a.u.)1100-01 TOKMDSFo-/9dQD//QQFHTHHTT1T91
whole absorption spectrum can be hardly computed. sit
We now move on to describe the shift of the plasmon en10-2
D
resonance as a function of the particle size. This prob- ed 10-3
lem has been extensively studied in the literature, both c
u
theoretically41,105 and experimentally106–108, and it rep- d
resentsarelevantbenchmarkforestimatingtheaccuracy In10-4
of QHT. In Fig. 8 we report the energy position of the -5
10
main resonance peak as a function of the inverse of the 30 40 50
Distance (bohr)
jellium nanosphere diameter, D. The exact energy posi-
tion of the main peak has been extracted from the com-
FIG. 7. Induced density (complex modulus of n ) at the
puted spectra (with an empirical broadening of (cid:126)γ =0.1 1
plasmon energy for a jellium nanosphere (r =4 a.u.) with
s
eV and using a spline interpolation). These procedure is
N =338electrons(R=27.86bohr),ascomputedfromTD-
e
not unique for some of the smallest clusters, where there DFT, OF9/QHT9, KS/QHT1 and Mod/QHT1.
are many peaks with similar intensity (for large cluster
there is always an unique main peak).
Thedot-dashedhorizontallinerepresentstheMieplas- systems33,109,110. For all the other cases the peak reso-
mon energy ((cid:126)ωMie =3.4 eV). Note that for the diame- nancesslidetolowerenergies(longerwavelengths)asthe
ters considered in Fig. 8, retardation effects can be ne- particle radius shrinks. While for noble metals like Au
glected. or Ag the plasmon resonance undergoes a blue shift as
The first thing to notice is the striking difference ob- the radius R decreases111,112, this is not the case for Na.
tained using TF-HT. It predicts in fact a resonance shift Theoriginoftheblueshiftfornoblemetalnanoparticles
toward higher energies (shorter wavelengths) as the par- is due to size dependent changes of the optical interband
ticle radius gets smaller, as previously observed in other transitions113,114.
9
N
e R=25.0 R = 25 nm
0
2 8
3 3 2 4
5 3 9 3 8 R=17.5
4.5
TF-HT
Mie n R=12.2
o
KS/QHT1 ti |E|
TD-DFT p
4 MOFod9//QQHHTT19 3.3 bsor R=8.5
a R=6.0

d
eV) gy (eV) alize R=4.2 Mie R = 1 nm
( er3.2 m
gy3.5 En or R=3.0 FT
ner N D-D
E R=2.0 T
0.1 0.2 0.3 R=1.5 |E|
3 1/D (nm-1) R=1.0
2.8 3 3.2 3.4 3.6
Energy (eV)
2.5 FIG. 9. Mod/QHT1 spectra for jellium nanospheres ranging
0 0.5 1 1.5 from R = 2 to 25 nm (red curves). The black (blue) curve
1/D (nm-1) shows the Mie (TD-DFT) resonance trajectory.
FIG. 8. Plasmon resonance for jellium nanospheres (r =4
s
a.u.) as a function of the inverse of the sphere diameter, as range of nanoparticle sizes. The comparison between
computed by different approaches. The behavior for large
TD-DFT and KS/QHT1 is important because both ap-
particles is shown in the inset.
proaches use the same KS density (in the former addi-
tionalinformationisusedfromtheKSorbitalandeigen-
values). The good accuracy in Fig. 8 means that for
We now compare the QHT models investigated in this large nanospheres the full TD-DFT linear response can
work, with respect to the TD-DFT results, which can be bewellapproximatedbythesimplerQHT1method. The
considered as a reference. As widely investigated in the verygoodresultsobtainedforMod/QHT1areevenmore
literature for jellium nanospheres105, the TD-DFT main important. In fact it means that QHT can be use with-
peak oscillates for small Ne, but it converges to (cid:126)ωMie out the need of calculating the KS ground-state density,
for large Ne. which is a bottleneck for large system (scales as O(Ne3)).
Results for the OF9/QHT9 approach are signifi- Moreover, the simple model density employed here can
cantly blue-shifted with respect to TD-DFT (note that be constructed at no cost for systems of any size, and
OF9/QHT9 predicts a red-shift with respect to TF-HT, can also be generalized to the non-spherical case. The
as also found in Ref. 75), and do not present quantum parameter κMod, is clearly material-dependent (e.g. will
oscillations. In fact, it is well known115 that orbital-free dependonrs)butcanbeparameterizedonceandforall.
(OF1orOF9)electronicdensitydoesnotshowquantum InFig.9weshowtheabsorptionspectrumforparticles
(i.e. Friedel) oscillations inside the nanosphere (see inset with a diameter going from 2 up to 50 nm. The solid
of Fig. 3). On the other hand when QHT1 is applied to black line shows the trajectory of the Mie resonance as
the KS density, quantum-oscillations are clearly visible the particle size increases: for particles with R > 10
for small nanospheres, even if TD-DFT features are not nm the Mie resonance peak undergoes a red-shift due
fully reproduced. For Ne ≥ 338, KS/QHT1 reproduces to retardation effects. The red curves represent the
TD-DFT plasmon energies with great accuracy (with a spectrum calculated within the Mod/QHT1 method.
maximum error of 20 meV, about a half of the error ob- For small nanoparticles the peaks follow the TD-DFT
tained with the OF9/QHT9 approach, see Table S1 in trajectory, moving toward higher energies. As the
the Supplemental Material). particle size grows the plasmon energy tends toward
Finally, we analyze results for Mod/QHT1. Also in the Mie trajectory, up to the big particle regime, where
this case no quantum-oscillations are present, and for retardation effects become predominant (see Fig. 8). It
N ≥ 338 TD-DFT results are reproduced almost ex- is striking how Mod/QHT1 can describe the full range
e
actly,withamaximumerrorofonly10meV.Thususing of effects going from the nonlocal/spill-out effects up to
a simple model density it is possible to match the whole retardation effects. That is, the resonance shift due to
10
z
N = 508 TD-DFT
) e a)
. KS/QHT1
u
. 0.04 Mod/QHT1
a
( E

]
α
[ 0 k
m
I
-0.04
N = 1074
e
)
.
u n n+ E E
a. 0.04 4 1 4 0 25
( b) c)
] 3 100 3
[α 0 2 2 20
m
50
I m) 1 m) 1 15
-0.04 n 0 0 n 0
N = 4570 z (−1 z (−1 10
) e −50
. −2 −2
u
a. 0.04 −3 −100 −3 5
(
] −4 −4 0
α
[ 0 −2 −1 0 1 2 −2 −1 0 1 2
m y (nm) y (nm)
I
-0.04
FIG. 11. Dimer of Na spheres constituted by N =398 elec-
0 20 40 60 80 e
tronseach,andseparatedbyadistanceg=1nm. Thedimer
Distance (bohr) is excited by a plane wave oscillating with energy (cid:126)ω=−2.8
eV, γ = γ +v /R has been used with (cid:126)γ = 0.066 eV. (a)
0 F 0
Geometry and incident field. (b) Real part of the induced
FIG. 10. Induced polarization charge density (normalized charge density normalized with respect to the bulk positive
by R3) at the plasmonic resonance ω0 for different jellium charge density. (c) Electric field norm distribution.
nanospheres.
the incident field E =ˆzE can be defined from:
0 0
microscopic and macroscopic effects are incorporated
1
in a single model, which makes the potential of QHT n (r,θ,φ)=−eE cos(θ)α (r,ω) (26)
1 0r2 zz
with respect to DFT approach very clear. Although this
advantage was already outlined by the authors of Ref. so that
75, we remark that their method was only qualitatively
4π (cid:90) +∞
verified for nanowires. α (ω)= α (r,ω)rdr. (27)
zz 3 zz
0
In Fig. 10 we plot the imaginary part of the induced
polarization charge density for different particles sizes in
correspondence with the plasmon resonance ω . For the
VII. THE INDUCED CHARGE DENSITY 0
smaller particles big oscillations of the density can be
seeninthecaseoftheTD-DFTcalculationsthatarenot
SofarwehaveseenthatKS/QHT1orMod/QHT1re- presentifnotinaverymodestfrominthecaseofQHT1.
produce with very good accuracy the reference TD-DFT These oscillations are in fact due to a purely quantum
results, both the energy position and the asymptotic de- size effect (Friedel oscillations). As the particle size in-
cay of the induced density. In this section we closely an- creases however, these oscillations diminish. The main
alyzethenear-fieldproperties. Accurateinduceddensity induced peak, however, is very well reproduced by the
translatesintoagooddescriptionofthelocalfieldsatthe QHT1approach,bothwiththeKSorthemodelground-
surface of the plasmonic system. Such knowledge is cru- state density.
cialforestimatingthemaximumfieldenhancements,and
hence nonlinear optical efficiencies, and more in general
light-matter interactions. VIII. APPLICATION TO THE SPHERE DIMER
In particular we compare the QHT1 induced polariza-
tionchargedensitytothefullTD-DFTcalculations. The Up to this point we have considered spheres. In this
polarization charge density α (r) of a sphere excited by section we are going to extend the applicability of QHT
zz