Scholarly article on topic 'Spatially resolved collective excitations of nano-plasmas via molecular dynamics simulations and fluid dynamics'

Spatially resolved collective excitations of nano-plasmas via molecular dynamics simulations and fluid dynamics Academic research paper on "Nano-technology"

Share paper
Academic journal
New Journal of Physics
OECD Field of science

Academic research paper on topic "Spatially resolved collective excitations of nano-plasmas via molecular dynamics simulations and fluid dynamics"



Home Search Collections Journals About Contact us My IOPscience

Spatially resolved collective excitations of nano-plasmas via molecular dynamics simulations and fluid dynamics

This content has been downloaded from IOPscience. Please scroll down to see the full text. View the table of contents for this issue, or go to the journal homepage for more

Download details: IP Address:

This content was downloaded on 12/03/2015 at 12:10 Please note that terms and conditions apply.

New Journal of Physics

The open access journal for physics

Spatially resolved collective excitations of nano-plasmas via molecular dynamics simulations and fluid dynamics

T Raitza1'2'6, H Reinholz3'4'6, P-G Reinhard5 6, G Röpke3 and I Broda1

1 Johannes Kepler Universität Linz, Altenberger Strasse 69, 4040 Linz, Austria 2Laboratoire de Physique Theorique, IRSAMC, Universite Paul Sabatier, 118 Route de Narbonne, 31062 Toulouse, France

3 Institut für Physik, Universitat Rostock, 18051 Rostock, Germany

4 School of Physics, University of Western Australia, WA 6009 Crawley, Australia

5 Institut für Theoretische Physik, Universitat Erlangen, Staudtstrasse 7, 91058 Erlangen, Germany

E-mail:, and

New Journal of Physics 14 (2012) 115016 (14pp)

Received 21 August 2012 Published 21 November 2012 Online at


Abstract. Collective excitations in nano-plasmas are described by dynamical bi-local auto-correlation functions. These excitations, which are related to the plasmon excitations in bulk plasmas, arise in the classical as well as the quantum regime. Instead of the wave-vector-dependent dynamical structure factor, which is not well defined in finite systems, two different signatures are considered to characterize collective excitations: the bi-local particle density correlation function and the bi-local current density correlation function. The relation between both signatures is not as trivial as in the homogeneous case and is given here. Exemplary calculations are performed for expanding nearly spherical clusters of sodium atoms after excitation by a high-intensity short pulse laser beam. The lowest collective excitations obtained in the classical regime using molecular dynamics simulations agree well with the lowest collective excitations

6 Authors to whom any correspondence should be addressed.

Content from this work may be used under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

New Journal of Physics 14 (2012) 115016 1367-2630/12/115016+14$33.00

© IOP Publishing Ltd and Deutsche Physikalische Gesellschaft

obtained from quantum calculations using fluid dynamics. The energy, damping and structure of the lowest collective modes are given. The dynamical bi-local correlation functions are of relevance for the optical properties, in particular the determination of photo absorption coefficients of nano-plasmas.


1. Introduction 2

2. Definition of correlation functions 4

3. Molecular dynamics simulations of collective excitations 5

4. Collective modes via fluid dynamics 7

5. Resonance structures of the local current density and local density fluctuation 10

6. Conclusions and outlook 12 Acknowledgments 13 References 13

1. Introduction

The investigation of small clusters with respect to collective plasma excitations has been of increasing interest in recent years [1, 2]. In particular, nano-plasmas are formed in metallic clusters after irradiation with intense short pulse laser beams [3]. In previous works [4, 5], a complex excitation spectrum of hot weakly bound electrons inside clusters was obtained experimentally; however, there is still a lack of theoretical understanding. Experimental investigations on the surface plasmon amplification of nano particles [6] suggest applications such as a nanoscale quantum generator and an ultrafast amplifier. Unusual resonances in nano-plasmonic structures [7] were found as a non-local response after the excitation of electrons inside nano cylinders.

We are interested in studying nano-plasmas in clusters considered after short pulse laser irradiation with intensities of the order of 1012 Wcm-2. In earlier investigations [8], the role of the additional ionization and heating due to collective effects inside clusters has been discussed. We address the collective effects of the nano-plasma in thermal equilibrium after the excitation process. The plasma properties of finite systems are mainly determined by the dynamics of the electrons that are bound to the cluster but delocalized from the original atoms. They are comparable to the conduction electrons in bulk systems. As was already shown in earlier publications, see [9], the electrons relax to local thermodynamic equilibrium (LTE) within time scales of a few fs after the laser is switched off. Therefore, the temperature T and the electron density p (r) as used in bulk plasmas can serve as the characteristic parameters of the electron system (nano-plasma) after LTE is established. Additionally, the cluster size (the number of ions N{) and the net charge Z have to be considered in finite systems.

Scattering, absorption of radiation and corresponding phenomena are determined by collective excitations of plasmas. The volume plasmon is well known as a collective excitation in an infinite, homogeneous plasma that occurs at eigen-frequencies w2(k) = (1 + 3k2/k2 + • • •), see [10], with the wave number k, the plasmon frequency = ne2/s0me and the Debye screening parameter k2 = ne2/s0kB T. In the case of homogeneously charged spheres the strongest signal is the dipole mode with the Mie resonance wMie = wpl/V3 [1,2].

So far [11-13], we have obtained excitation spectra for bulk and cluster materials via correlation functions. In the case of clusters, the correlation functions were analyzed using spherical harmonics. In order to bridge from finite systems to bulk plasmas, first results with respect to size effects have been reported in [9, 14]. In particular, we are interested in the dynamical structure factor and the response function for such finite nano-plasmas, which are related specifically to correlation functions of the spatially resolved electron density fluctuations. After some general remarks on the relations of different correlation functions in section 2, two different theoretical approaches—molecular dynamics (MD) simulations and fluid dynamics calculations—will be presented here to discuss collective excitation modes.

In section 3, the results for frequency-dependent bi-local current density correlation functions calculated via MD simulations are presented. For this approach, excited electrons in clusters at temperatures larger than the Fermi energy are considered, where the nano-plasma can be treated classically. Strong correlations are taken into account via collisions of all particles. Concepts such as the dynamical collision frequency that have been well established for infinite bulk systems near thermodynamic equilibrium [15] have to be modified for applications to finite systems, e.g. the clusters considered here. Therefore, a restricted MD scheme [9, 16] was used, where the ion configuration is frozen representing thermodynamic parameters at a particular instant of time within the cluster expansion. The dynamics of electrons inside the cluster can then be calculated separately. Using the ergodic theorem, the ensemble averaging over the classical trajectory of multiple simulations of the same physical process can be replaced by temporal averaging over a single trajectory. From the spatially resolved time-dependent current density j(r, t), the bi-local current density correlation function was calculated in the time domain. The corresponding frequency-dependent bi-local correlation functions are accessible after Laplace transformation. Finally, we are able to treat the optical response of a thermalized nano-plasma via bi-local current density correlations (j(r); j(r7)>^ in the frequency domain. For arbitrary degeneracy, the Laplace transform of the correlation function,

is defined with the equilibrium statistical operator p0, and 5 = 1/ kB T.

Alternatively in section 4, we consider a fully quantum mechanical T = 0 description in terms of time-dependent density functional theory (TDDFT), see, e.g., [2]. The dynamics of excitations is treated in the linear domain, which allows a straightforward sorting in terms of eigenmodes and associated eigenfrequencies. The excitation spectrum thus obtained will be used to derive the various correlation functions as were considered in the analysis of the MD simulation. This then allows a direct comparison of MD and TDDFT. For the quantum case, we consider the solution spectrum in random phase approximation (RPA) and fluid dynamical approximation [13, 17]. The latter relies on the full quantum mechanical ground state but includes approximations of the excitations in terms of local flow pattern, in this way establishing a bridge between full quantum mechanical treatment and classical MD. This will be discussed in section 5.

We would like to stress that collective excitations are not specific quantum phenomena but appear also in the classical case, as is well known from the bulk plasmons. Similar collective resonances in both approaches are expected. We therefore attempt a comparison of the results of quantum description and those of the classical approach.

Within MD simulations, the interaction of electrons with cluster ions is described by the error-function potential

Vei (r) =

where X = 6.02a0 was used in order to reproduce the ionization energy of Vei(0) = -5.1 eV for solid sodium. The ion positions are randomly distributed within the cluster sphere at the solidstate density. The resulting ionic potential is not sensitive to the actual position of the ions due to the smooth feature of the error function potential.

For TDDFT and fluid dynamics, the jellium-like potential

j (r) =

d ri Pjei(|r - ri |) — ri

was used. The Coulomb interaction is convoluted with a soft surface [18] using the Wood-Saxon-type radially symmetric ionic density distribution [19]

Pjel (r) =

1+exp [(r - rsN1/3)/a]'

where p0 has to be determined from the normalization Ni = f d3rpjel(r). For the jellium potential Vjel(r) to be in accordance with the error function potential used for MD simulations, the surface width is taken as a = 1.6a0.

2. Definition of correlation functions

As correlation functions of fluctuations around thermodynamic equilibrium will play a key role in the following discussions, we briefly state the major relevant properties. We treat the near-equilibrium system in linear response in terms of the generalized Zubarev formalism [20]. Within this approach, induced quantities in non-equilibrium due to external fields are determined via equilibrium correlation functions (A(t); B) and the frequency spectrum lime^0 (A; B)w+ie, after Laplace transformation. We introduce the external Hamilton operator

Hxt = / d3 r pWext (r, t), (5)

due to an external scalar field Uext(r, t). This induces a density fluctuation (¿p(r))t, ¿p(r) = p(r) — (p(r)), which reads in the frequency domain

<«p(r)>* = j d3r7 x(r, r, w) Uext(r7, w), (6)

where x(r, r7, w) = ^ (5p(r); (r7))w is the response function as also introduced within the Kubo formalism. Using partial integration, one is able to express the dynamic response function

X(r, r7, w) = (Sp(r); ¿p(r7))|t=0 + iw (Sp(r); ¿p(r7))w (7)

via the respective bi-local density correlation function; for details see [11, 20]. Alternatively, one can describe the excitation by considering the current density-dependent Hamiltonian

Hxt = / d3 rj(r)A(r, t), (8)

with the time derivative of the vector potential A(r, t) = c2grad Uext(r, t) = —c2E(r7, t). Accordingly, a current density will be induced,

<8j(r))o = / d3r7 a(r, r7, o) E(r7, o), (9)

where a(r, r7, o) = ^ <j(r); j(r7))o is known as the conductivity tensor. Using the continuity equation, the bi-local current density correlation function can also be related to the response function equation (7)

X(r, r7, o) = — - divr divr (j(r); j(r7)) , (10)

io x /o

where divr is the tensor divergence with respect to r. Simulation data of different approaches are now used to calculate the bi-local current density correlation function <j(r); j(r7))o, on the one hand, and the bi-local density correlation function <8p(r); 8p(r7))o, on the other hand. Our aim is to show the relation between the bi-local current density and the bi-local particle density correlation function using the expressions in equations (7) and (10).

The most prominent feature of electronic modes in metal clusters is collective excitations (such as, e.g., the Mie plasmon) which appear as sharp resonance peaks in the spectrum. The spatial properties of the corresponding modes are described by the electronic transition density or the transition current density. We will discuss the relation between these two spatially resolved distribution functions for the strongest collective modes.

3. Molecular dynamics simulations of collective excitations

Starting from MD simulations [21], a restricted MD scheme [11] was introduced for the cluster case, treating the electron dynamics within a fixed ion configuration due to adiabatic approximation. Exemplarily, we discuss the collective resonances for an Na+00o cluster immediately after laser irradiation at an average ion density ni = 2.8 x 1022cm—3. This is the solid-state density of sodium and corresponds to a Wigner-Seitz radius of rs = 3.9a0. The initial electron distribution is thermalized at a temperature T = 1 eV. The electron dynamics has been calculated considering the trajectory from tMD = 10 • • • 100 ps using time steps of 81 = 0.1 fs. The z component of the total current density of all electrons inside the cluster

« 1 ^ z

jz (t) = —- j^pz (t) (11)

has been determined from the electron momenta pf (t), with a the cluster volume. We calculated the current density auto-correlation function

1 P tMD

(jz; jz(t)) =-^ drjz(T) jz(t + T). (12)

tMD\ J 2) J0

Figure 1. Top: the total current density correlation function (13). Bottom: the leading modes of bi-local correlation matrix Ka,b for spatially resolved current density; both for the Na10oo cluster at electron temperature T = 1 eV, cluster charge Z = 19 and ionic density rs = 3.9a0.

The spectrum of the auto-correlation function

( jz; jzL = / dt eiwt (jz; jz (t)), (13)

see figure 1 (top), shows a main maximum at w = 0.23 Ry marked with a dashed vertical line. It is related to the Mie resonance. Additional structures are visible at higher frequencies.

The complete response function for finite systems, which resolves the resonances spatially,

can be determined from the bi-local correlation functions of the current density Kjr (&>) = (jz(r); jz(r7))w or density fluctuations K^= (¿p(r); ¿p(r7))^.

For the spatial resolution, the cluster volume was divided into a finite number of cells (a, b) using spherical coordinates, centered at (r, r7). The discrete correlation matrix Kab (&>) was then analyzed by solving the following eigenvalue problem:

^ Re KabM = KM(^M,aM, (14)

for the leading eigenvectors They characterize the spatial structure of the individual

modes (^M,a(&>) ^ ^M(r, &>)). The strength of each mode is given by the eigenvalues KM(&>). Details are given in [11], where the corresponding mode analysis of the response function has been presented for a different cluster size (Ni = 309).

The decomposition of the locally resolved current correlation matrix into eigenvalues KM(&>), see figure 1 (bottom), shows a fairly complex structure of resonance modes. To analyze the spatial structure within the cluster, a Fourier decomposition of the eigenvectors

lOP Institute of Physics <j)deutsche physikalische Gesellschaft a) mode nr. 1 o>=0.23 Ry b) mode nr. 2 o>=0.30 Ry c) mode nr. 3 o>=0.44 Ry

-0.5 -1

-0.5 -1

-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 x [R] x [R] x [R]

Figure 2. Selected eigenvectors ^M(r, of the dipole-like mode (black line in figure 1, bottom) of the Na|00o clusters with ionic density rs = 3.9a0 and temperature T = 1 eV. The spatial x- and z-axes are scaled in units of the cluster radius R = 39.0a0.

into the spherical Bessel functions jl (knlr) and spherical harmonics Yl,m (0, 0) was performed according to

Nn Ni l

= ^^ ^ Sn,i,mM Nn j(kn,ir)Yim(0,0), (15)

n=1 l =0 m=—l

where Snj,m is the spherical Fourier component with ordinal numbers n, l, m. The normalization factor Nnl as well as the wave number knl are chosen in the way the eigenvector has a node at the cluster surface.

In figure 1 (bottom), the four strongest eigenvalue modes are shown. They are characterized by pairs of ordinal numbers l, m, which determine the angular part of the eigenvector in terms of the spherical harmonics Yl,m (0, 0). The leading dipole-like mode, represented via a solid line (black online), is characterized by an overlap of Y0 0(0, 0) and Y2,0(0, 0). We find three resonance frequencies which are identical to the ones found in the total current density autocorrelation function. They are indicated by vertical dashed lines (blue online).

In a next step, we analyze selected eigenvectors ^M(r, &>) of the dipole-like mode, which is the strongest mode shown in figure 1 (bottom). The results are given in figure 2. At the resonance frequency ^ = 0.23 Ry, figure 2(a), all electrons are moving in the same direction. No nodes can be seen, as one would expect in the case of a Mie oscillation. The mode at the third resonance frequency, see figure 2(c), is similar to a plane wave oscillation of electrons, but trapped inside the cluster. To identify a wave number of the plane wave oscillation, a Fourier decomposition into plane waves has been done. A characteristic plane wave number in the z-direction of k = 1.6 nm-1 was found. The resonance structure in figure 2(b) appears as a mix of the first and the third resonance structure.

4. Collective modes via fluid dynamics

To describe collective modes at low temperatures, a quantum mechanical approach has to be used. In this section, we investigate the electron dynamics of the cluster at T = 0 using the time-dependent density functional theory [2]. The electron cloud confined by the potential of the ion distribution equation (4) is described in local-density approximation (LDA) [22] using

the density functional of Perdew and Wang [23]. The excitation spectrum is computed with linearized time-dependent LDA, also known as random-phase approximation (RPA), using the efficient operator techniques developed by Reinhard et al [13, 17]. It yields a spectrum of excitation operators Cj and associated eigenfrequencies &>v. The operators are optimized by variation with respect to Cv, i.e.

&cv<Oq| [Cv, [H, Cj]] |Oo)= 0, (16)

where |O0) is the ground state. A special approximation is determined by the class of operators CJ considered for the variational problem.

A simple fluid dynamical picture [13, 24] is realized by choosing the excitation operators from a purely local distribution Q (r) as Cj + Cv e {Q (r)}, CJ — Cv a i[H, Q]. The variationally optimized field Q (r) represents a velocity potential such that the current associated with that excitation mode becomes j(r) = p(r)V Q. Therefore, it is an irrotational flow, V x (j/p) = 0. Fluid dynamics is a collective model in the sense that only the global flow determines the system while details of single-particle motion are ignored.

The full single-particle picture is maintained in RPA, which is realized by allowing a composition of all conceivable one-particle-one-hole (1ph) states Cj = J]oh (xp/*apah +

ahap), where a', a are fermion operators. The superposition coefficients are, again, to be determined variationally by equation (16) where variation with respect to Cv means, in practice, variation with respect to x^h)* and y^*.

The system which we consider here is modeled with spherical jellium background and has an electron number which yields a closed shell configuration. The many-electron ground state is spherically symmetric. As a consequence, the excitations can be sorted according to angular momentum (L, M). We will indicate that by the notation C\LM. Key observables for an excitation are the excitation frequency o>nLM and the corresponding multipole strength. More detailed information on the structure of excitation is provided by the transition density 8pnLM and transition current jnLM defined as

SPnLM(r) = <Oo|[p(r), C¡LM]|®0>, jnLM(r) = |[j(r), C^m]|®0>, (17)

where p(r) and j(r) are the operators of local particle density and current density, respectively. Spherical symmetry leads to the separable form

$PnLM (r) = SpnLM (r )Ylm (ft), (18a)

jnLM (r) = ^ jmi(r )Y LIM (ft), (18b)

I G{ L-1, L+1}

Y LIM (ft) = £ YLm (ft)eq V2Z+1 (-)L-1+M(m q

11 M . (18c)

m q - M) v y

where YLM are spherical harmonics, eq the three-dimensional (3D) unit vector in spherical representation and YUM vector spherical harmonics [25]. The continuity equation &>vLM8pnLM = V- jnLM admits only two independent radial distributions. Fluid dynamics, as mentioned above, produces laminar flow which ties together j(Land j^LM) and reduces the number of independent radial distributions to a single one. In any case, the transition density 8pnLM (r),

being a scalar and 1D function, is the simplest object for a detailed visualization of the transition. Finally, this is related to the bi-local density correlation as

(<5p(r); ¿p(r')}m = £ _ ^)2 + r2Sp*nL(r)SpnL(r')PL(cos0), (19)

with the Legendre polynomials PL (x); 0 is the angle between r and r7. The folding width T is taken as T = 0.01 Ry. For a better overview, we produce an w weighted and L-dependent radial density by superposing the various spectral states with a weight given by the corresponding transition matrix element and some smoothing

8pL (r; w) = V WnL--2 Sp*L (r) MnL, (20a)

^ (w _ WnL)2 + T2

MnL = / dr7 r72+L 8pnL (r7). (20b)

This form has the advantage that photo-absorption strengths and the corresponding sum rule can be obtained by integration over r and w.

Figure 3 shows the frequency-dependent transition densities equation (20a) for a test case Na36+. Fluid dynamical L = 1 and 2 modes are visualized and a comparison with full RPA is given for L = 1. The r2 weight has been applied to show directly the contributions to the corresponding multipole momenta. Two angular momentum modes are shown as indicated: dipole (L = 1) on the left panels and quadrupole (L = 2) on the right panel.

The mode at the lowest energy in the fluid dynamical picture for L = 1 (lower left panel) is clearly the Mie surface plasmon. Above that energy a second surface mode is found. Higher modes acquire more and more volume components [26] (a feature which is, however, not well visible at the present color scale). The transition densities of full RPA (upper left panel) show the same two dominant surface modes as the fluid dynamical approach in spite of the fact that the pattern indicates, of course, a much richer spectrum behind. Each one of these many detailed particle-hole states carries a small amount of strength. It is impressive that the entity of these states forms together a pattern similar to the collective model. This effective collectivity holds even for the volume modes found at higher frequencies. One can spot it for the mode at 0.41 Ry which is just visible in full RPA at a given color scale. However, one encounters a growing fragmentation with increasing energy such that separate volume modes become hard to uncover in a fully fledged RPA spectrum.

The right panel of figure 3 shows the fluid dynamics spectrum for the quadrupole mode. It looks very similar to the dipole modes. But all eigenfrequencies for L = 2 are systematically slightly higher than for L = 1. Such an up-shift can already be estimated from the simple Mie picture which predicts the frequency of the leading surface plasmon as wplVL /(2L+1) where L is the angular momentum of the mode [26]. This predicts an up-shift by about 10% when going from L = 1 to 2. The actual up-shift seen in figure 3 is smaller due to the soft jellium surface profile.

After all, we see that the fluid dynamics picture provides a pertinent model of the dominant flow pattern in the dynamics of a large cluster as being studied here. Fluid dynamics is, in fact, a classical concept and so we find a classical description validated by the results. It is then to be expected that fully classical MD simulations should display the same dominant collective pattern in the excitation spectrum.

to 0.4

<5 0.3

0.2 0.1 0

y] [R 0.5

ti a t 0.2

X e 0.1

IOP Institute of Physics <j)deutsche physikalische Gesellschaft full RPA L=1

10 20 30 fluid dynamics L=1

10 20 30 radial distance [a0]

transition density r2*5pL(o>) for

Mo 16+

soft jellium rs =3.9 a0, a=1.6 a0

fluid dynamics L=2

40 0 10 20 30 radial distance [a0]

Figure 3. Radially weighted, frequency-dependent transition density r28pL (r, according to equation (20a) for Na16+ as a color-scale plot: dipole modes (L = 1) on the left panels and quadrupole mode (L = 2) on the right panel, fluid dynamics results in the lower panel and full RPA in the left upper panel. The system radius is indicated by a faint vertical line. The frequencies of the first three fluid dynamical L = 1 eigenmodes are indicated by a faint horizontal line.

5. Resonance structures of the local current density and local density fluctuation

Using the equation of continuity ^8p(r) = div j(r), the current density mode j(r) is related to the density fluctuations mode 8p (r). A similar problem was considered earlier for very small cluster sizes in [27]. Here, we consider a more general description relevant for arbitrary excitation modes. Because the current density represents a vector field, its spatial dependence must be decomposed into vector spherical harmonics YL, LM(0,0) as they were introduced in [28], resulting in two radial profiles j— (r) and j+ (r) for a given angular momentum number L,

jL (r) = j— (r) YL — 1,LM(0, 0) + j+ (r) YL+1,LM(0, 0). (21)

The density fluctuation (see equation (18a))

8pL(r) = 8pL(r) YlM(0,0), (22)

where 8pL (r) gives the radial dependence of the respective oscillation mode, is related to the current density in the following way. The current density jL (r) is composed of angular parts of

vector spherical harmonics L + 1 and L — 1. We calculate the radial dependence of the current density from the density fluctuation via

L +1 / d L + 2\ / L id L —1\

'"pl (r > = V 2T+T (¡7 + —) J+(r > — V 2TT^(;J7 — —)'—(r >■ (23)

Note that rot j(r> is not determined by the density distribution. Assuming an irrotational current density field, one is able to find solutions for the radial functions of the current density separately,

j+ (r > = (1 — a >A/ dr1 rL+2 8pl (r1 >,

2 L + 1 m /,<x> L + 1 r

. 2L + 1 l—1 r a ¿PL (r1 >

j—(r > = aj—— m r J dr1 rL—1

where 0 ^ a ^ 1 appears as a free mixing parameter since the current density profile is not unambiguously determined by the density profile. Using equation (21), the complete solution of the current density in terms of the density fluctuations then reads

Y L+1, LM (g, —(

V L + 1

j L (r) = V2 L + 1

/1 xYL+1, LM _(L +2) A L +2 O / \

(1 — a>-—r (L+2> / drx rL+2 ¿pl (rx >

V L + 1 Jr

YL — 1,LM(g,0) rL—1 r dr I VL Jr 1 rL—1 J

+a L—— rL—M dr^^^^r^ ■ (25)

The spatial resolution of the density fluctuation correlations resulting from our MD simulations is statistically not good enough for further analysis. Subsequently, the direct comparison of the current density from MD simulations with that obtained via the electron density eigenvector is viable. We will therefore focus on a dipole mode of the density fluctuation oscillation, which consists purely of spherical harmonics L = 1. Thus, in our case (pL (r) for L = 1) the complete solution of the current density is a superposition of the spherical harmonics 70,o (0, 0) and y2,o (0, 0) as it already appeared in the results for the strongest collective mode of the MD simulations. In figure 3 (bottom, left), the spatial dependence of such a density fluctuation mode is shown in dependence on the excitation energy. This was calculated via fluid dynamics.

Using the density fluctuation mode shown in figure 3 (bottom, left), the current density mode is determined from equation (25). The z component is shown in figure 4 for the resonance frequencies found in the fluid dynamics calculations. The resulting current density compares well with the eigenvectors from the bi-local current density correlation function for the resonant excitation energies calculated in MD simulations, shown in figure 2. Thus we conclude that the Mie-like character of the strongest collective resonance observed from density fluctuation correlations (5p(r); ¿p(r7via fluid dynamics as well as from the current density correlations (j(r); j(r7))w via MD simulations is in remarkable agreement.

The most important difference between the two approaches to obtain the resonant oscillation structure is due to the ground-state distribution of electrons considered inside the nano-plasma. In the case of MD simulations, electrons thermalized at a finite temperature T = 1 eV are considered. In contrast, the fluid dynamics ground state of the nano-plasma electrons corresponds to T = 0. Obviously, the oscillation mode structures and resonance frequencies are not strongly affected by thermal effects.

12 lOP Institute of Physics <J>deutsche physikalische Gesellschaft

a) mode nr. 1 œ=0.21 Ry b) mode nr. 2 œ=0.32 Ry c) mode nr. 3 œ=0.41 Ry

-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1

x [R] x [R] x [R]

Figure 4. Mode structure in the z-x-plane for resonance frequencies of the Na+56 cluster at ionic density rs = 3.9a0 and T = 0 calculated via fluid dynamics. Top: transition density modes 8pL(r) (L = 1). Bottom: current density modes jz(r) calculated from the transition density (above) via equation (25). The spatial x-and z-axes are scaled in units of the cluster radius R = 27.6a0.

The statistical fluctuations of the MD simulations are considerably lower for the calculation of the current density correlations than for the electron density correlation. This is also clear from figure 4. The averages of the current density vanish whereas the fluctuating average electron densities have to be subtracted to get the electron density fluctuations. Now, we are able to solve this problem by transferring the density modes into current density modes.

Beyond that, collisions between electrons are included via MD simulations which lead to a damping of collective electron oscillations. This is not included in the fluid dynamics scheme. An additionally introduced damping width of T = 0.01 Ry is relatively small compared to the MD simulations case. Nevertheless, at resonance frequencies, the oscillation mode structure is in agreement for both schemes. However, in comparison to 1D chains, see [29], where the spatial mode structures were simply allocatable by their wave number k, the current density mode structures of 3D clusters are more complex.

6. Conclusions and outlook

We have been successful in using the continuity equation to relate the frequency-dependent bi-local correlation functions of the current and charge densities. This allows us to compare different approaches to investigate collective excitations in finite systems and to analyze their spatial excitation structures. From MD simulations as well as from a time-dependent density functional approach, we were able to calculate resonance frequencies at similar frequencies. Our approach offers a generalization of the dynamic structure factor relevant in bulk systems that are homogeneous in space, to spatially resolved correlation functions applicable to finite

inhomogeneous systems. Collective excitations in highly excited nano-materials are relevant for the response to external perturbation, in particular considering the emission and absorption of light.

A systematic treatment, as well as the proof of general properties such as the Thomas Reiche-Kuhn sum rules, is left for our future work. The relation of the total current density autocorrelation function to the dielectric function is given by s-1 (a) = 1 — ißnC (jz; jz)a /(s0a), see [15]. Here, nC is the density of clusters. Also the bi-local response function that results from the bi-local current density correlation function is now accessible with MD as well as fluid dynamics simulations. Hence, we are able to discuss the bi-local conductivity of nano-plasmas which is directly related to the bi-local current density correlation function. Finally, these material properties allow us to discuss the role of collective excitations and the damping of these excitations in nano-plasmas, which has not been completely understood until now. For example, this is important for discussing experimental results which suggest still unknown effects of nano-plasmonics. Experimental consequences, in particular absorption spectra as well as scattering of light and charged particles, are therefore of interest.


We thank Igor Morozov (Moscow) and Eric Suraud (Toulouse) for discussions. We would like to acknowledge financial support from SFB 652, which is funded by the DFG. GR is grateful for a JKU research fellowship that provided funding during his visits to Linz.


[1] Kreibig U and Vollmer M 1995 Optical Properties of Metal Clusters (Berlin: Springer)

[2] Reinhard P G and Suraud E 2003 Introduction to Cluster Dynamics (New York: Wiley)

[3] Hüfner S 1995 Photoelectron Spectroscopy (Berlin: Springer)

[4] Xia C, Yin C and Kresin V V 2009 Phys. Rev. Lett. 102 156802

[5] Andersson T et al 2011 J. Chem. Phys. 134 094511

[6] Stockman MI 2010 J. Opt. 12 024004

[7] Raza S, Toscano G, Jauho A-P, Wubs M and Mortensen N A 2011 Phys. Rev. B 84 121412

[8] Bauer D and Macchi A 2003 Phys. Rev. A 68 033201

[9] Raitza T, Reinholz H, Röpke G, Morozov I and Suraud E 2009 Contrib. Plasma Phys. 49 498

[10] Thiele R, Bornath T, Fortmann C, Holl A, Redmer R, Reinholz H, Ropke G, Wierling A, Glenzer S H and

Gregori G 2008 Phys. Rev. E 78 026411

[11] Raitza T, Reinholz H, Morozov I and Ropke G 2011 Phys. Rev. E 84 036406

[12] Morozov I, Reinholz H, Ropke G, Wierling A and Zwicknagel G 2005 Phys. Rev. E 71 066408

[13] Reinhard P-G, Genzken O and Brack M 1996 Ann. Phys., Lpz. 5 576

[14] Reinholz H, Raitza T, Ropke G and Morozov I 2008 Int. J. Mod. Phys. B 22 4627

[15] Reinholz H 2005 Ann. Phys., Fr. 30 4-5

[16] Raitza T, Reinholz H, Ropke G and Morozov I 2009 J. Phys. A: Math. Theor. 42 214048

[17] Reinhard P-G and Gambhir Y K 1992 Ann. Phys., Lpz. 1 598

[18] Ekardt W 1984 Phys. Rev. Lett. 52 1925

[19] Montag B, Hirschmann T, Meyer J, Reinhard P-G and Brack M 1995 Phys. Rev. B 52 4775

[20] Zubarev D, Morozov V and Ropke G 1997 Statistical Mechanics of Nonequilibrium Processes vol 2 (Berlin:


[21] Calvayrac F, Reinhard P G, Suraud E and Ullrich C A 2000 Phys. Rep. 337 493

[22] Dreizler R M and Gross E K U 1990 Density Functional Theory: An Approach to the Quantum Many-Body

Problem (Berlin: Springer)

[23] Perdew J P and Wang Y 1992 Phys. Rev. B 45 13244

[24] Brack M 1989 Phys. Rev. B 39 3533

[25] Edmonds A R 1957 Angular Momentum in Quantum Mechanics (Princeton, NJ: Princeton University Press)

[26] Brack M 1993 Rev. Mod. Phys. 65 677

[27] Kümmel S, Andrae K and Reinhard P-G 2001 Appl. Phys. B 73 293

[28] Serr F E, Dimitrescu T S, Suzuki T and Dasso C H 1983 Nucl. Phys. A 404 359

[29] Raitza T, Reinholz H and Röpke G 2010 Int. J. Mod. Phys. B 24 4961