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

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


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 nanoplasmonic 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 10 12 W cm −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 ρ(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 i ) 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 ω 2 (k) = ω 2 pl (1 + 3k 2 /κ 2 + · · ·), see [10], with the wave number k, the plasmon frequency ω 2 pl = ne 2 /ε 0 m e and the Debye screening parameter κ 2 = ne 2 /ε 0 k B T . In the case of homogeneously charged spheres the strongest signal is the dipole mode with the Mie resonance ω Mie = ω pl / √ 3 [1,2].
So far [11][12][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(r ) ω in the frequency domain. For arbitrary degeneracy, the Laplace transform of the correlation function, is defined with the equilibrium statistical operator ρ 0 , and β = 1/k B 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 where λ = 6.02a 0 was used in order to reproduce the ionization energy of V ei (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 was used. The Coulomb interaction is convoluted with a soft surface [18] using the Wood-Saxon-type radially symmetric ionic density distribution [19] ρ jel (r ) = ρ 0 where ρ 0 has to be determined from the normalization N i = d 3 rρ jel (r ). For the jellium potential V jel (r ) to be in accordance with the error function potential used for MD simulations, the surface width is taken as σ = 1.6a 0 .

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 lim ε→0 A; B ω+iε , after Laplace transformation. We introduce the external Hamilton operator due to an external scalar field U ext (r, t). This induces a density fluctuation δρ(r) t , δρ(r) = ρ(r) − ρ(r) , which reads in the frequency domain where χ(r, r , ω) = β δρ(r); δρ(r ) ω is the response function as also introduced within the Kubo formalism. Using partial integration, one is able to express the dynamic response function χ (r, r , ω) = (δρ(r); δρ(r ))| t=0 + iω δρ(r); δρ(r ) ω 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 with the time derivative of the vector potentialȦ( Accordingly, a current density will be induced, where σ (r, r , ω) = β j(r); j(r ) ω 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) χ (r, r , ω) = − 1 iω div r div r j(r); j(r ) ω , where div r 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(r ) ω , on the one hand, and the bi-local density correlation function δρ(r); δρ(r ) ω , 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.

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 +19 1000 cluster immediately after laser irradiation at an average ion density n i = 2.8 × 10 22 cm −3 . This is the solid-state density of sodium and corresponds to a Wigner-Seitz radius of r s = 3.9a 0 . The initial electron distribution is thermalized at a temperature T = 1 eV. The electron dynamics has been calculated considering the trajectory from t MD = 10 · · · 100 ps using time steps of δt = 0.1 fs. The z component of the total current density of all electrons inside the cluster has been determined from the electron momenta p z l (t), with the cluster volume. We calculated the current density auto-correlation function Top: the total current density correlation function (13). Bottom: the leading modes of bi-local correlation matrix K a,b (ω) for spatially resolved current density; both for the Na 1000 cluster at electron temperature T = 1 eV, cluster charge Z = 19 and ionic density r s = 3.9a 0 .
The spectrum of the auto-correlation function 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 K j r,r (ω) = j z (r); j z (r ) ω or density fluctuations K δρ r,r (ω) = δρ(r); δρ(r ) ω . For the spatial resolution, the cluster volume was divided into a finite number of cells (a, b) using spherical coordinates, centered at (r, r ). The discrete correlation matrix K a,b (ω) was then analyzed by solving the following eigenvalue problem: for the leading eigenvectors µ,a (ω). They characterize the spatial structure of the individual modes ( µ,a (ω) → µ (r, ω)). The strength of each mode is given by the eigenvalues K µ (ω). Details are given in [11], where the corresponding mode analysis of the response function has been presented for a different cluster size (N i = 309).
The decomposition of the locally resolved current correlation matrix into eigenvalues K µ (ω), 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 into the spherical Bessel functions j l (k n,l r ) and spherical harmonics Y l,m (θ, φ) was performed according to where S n,l,m (ω) is the spherical Fourier component with ordinal numbers n, l, m. The normalization factor N n,l as well as the wave number k n,l 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 Y l,m (θ, φ). The leading dipole-like mode, represented via a solid line (black online), is characterized by an overlap of Y 0,0 (θ, φ) and Y 2,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 µ (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.

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 C † ν and associated eigenfrequencies ω ν . The operators are optimized by variation with respect to C ν , i.e.
where | 0 is the ground state. A special approximation is determined by the class of operators C † ν 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 . The variationally optimized field Q(r) represents a velocity potential such that the current associated with that excitation mode becomes j(r) = ρ(r)∇ Q. Therefore, it is an irrotational flow, ∇ × (j/ρ) = 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 where a † , a are fermion operators. The superposition coefficients are, again, to be determined variationally by equation (16) where variation with respect to C ν means, in practice, variation with respect to x (ν) * ph and y (ν) * ph . 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 † n L M . Key observables for an excitation are the excitation frequency ω n L M and the corresponding multipole strength. More detailed information on the structure of excitation is provided by the transition density δρ n L M and transition current j n L M defined as where ρ(r) and j(r) are the operators of local particle density and current density, respectively. Spherical symmetry leads to the separable form where Y L M are spherical harmonics, e q the three-dimensional (3D) unit vector in spherical representation and Y Ll M vector spherical harmonics [25]. The continuity equation ω ν L M δρ n L M = ∇ · j n L M admits only two independent radial distributions. Fluid dynamics, as mentioned above, produces laminar flow which ties together j (L+1) n L M and j (L−1) n L M and reduces the number of independent radial distributions to a single one. In any case, the transition density δρ n L M (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 δρ(r); δρ(r ) ω = n L (ω − ω n L ) 2 + 2 δρ * n L (r )δρ n L (r )P L (cos θ), with the Legendre polynomials P L (x); θ is the angle between r and r . The folding width is taken as = 0.01 Ry. For a better overview, we produce an ω 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 This form has the advantage that photo-absorption strengths and the corresponding sum rule can be obtained by integration over r and ω. Figure 3 shows the frequency-dependent transition densities equation (20a) for a test case Na 16+ 354 . Fluid dynamical L = 1 and 2 modes are visualized and a comparison with full RPA is given for L = 1. The r 2 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 ω pl √ L/(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.

Resonance structures of the local current density and local density fluctuation
Using the equation of continuity ωδρ(r) = div j(r), the current density mode j(r) is related to the density fluctuations mode δρ(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 Y * L ,L M (θ, φ) as they were introduced in [28], resulting in two radial profiles j − (r ) and j + (r ) for a given angular momentum number L, The density fluctuation (see equation (18a)) where δρ L (r ) gives the radial dependence of the respective oscillation mode, is related to the current density in the following way. The current density j L (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 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, 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 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 (ρ L (r) for L = 1) the complete solution of the current density is a superposition of the spherical harmonics Y 0,0 (θ, φ) and Y 2,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 δρ(r); δρ(r ) ω via fluid dynamics as well as from the current density correlations j(r); j(r ) ω 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. cluster at ionic density r s = 3.9a 0 and T = 0 calculated via fluid dynamics. Top: transition density modes δρ L (r ) (L = 1). Bottom: current density modes j z (r) calculated from the transition density (above) via equation (25). The spatial xand z-axes are scaled in units of the cluster radius R = 27.6a 0 .
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 = 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.

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 ε −1 (ω) = 1 − iβn C j z ; j z ω /(ε 0 ω), see [15]. Here, n C 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 nanoplasmas 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.