lopscience

¡opscience.iop.org

Home Search Collections Journals About Contact us My IOPscience

Optical response of C60 fullerene from a time dependent Thomas Fermi approach

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: 131.91.169.193

This content was downloaded on 06/08/2015 at 12:24 Please note that terms and conditions apply.

J. Phys. B: At. Mol. Opt. Phys. 48 (2015) 185102 (8pp) doi:10.1088/0953-4075/48/18/185102

Optical response of C60 fullerene from a time dependent Thomas Fermi approach

D I Palade1,3 and V Baran2

1 National Institute of Laser, Plasma and Radiation Physics, PO Box MG 36, RO-077125 Magurele, Bucharest, Romania

2 Faculty of Physics, University of Bucharest, Romania

E-mail: dragos.i.palade@gmail.com and baran@nipne.ro

Received 9 March 2015, revised 11 June 2015 Accepted for publication 14 July 2015 Published 5 August 2015

CrossMark

Abstract

We study the collective electron dynamics in C60 clusters within the time dependent Thomas Fermi method in the frame of jellium model. The results regarding the optical spectrum are in good agreement with the experimental data, our simulations being able to reproduce both resonances from 20 eV and 40 eV. We compare also, the results with those from other theoretical approaches and investigate the implications of quantum effects including exchange-correlation corrections, or gradient corrections from a Weizsacker term. The nature of the second resonance is studied using transition densities.

Keywords: Thomas Fermi, Fullerene, Optical, C60

(Some figures may appear in colour only in the online journal)

1. Introduction

The discovery of C60 fullerene in 1985 [1] confirmed Buck-minster's Fuller mental construction [2] of a structure with high symmetries and outstanding stability properties. Since then, many investigations were focused on the properties of fullerene including optical, electronic, thermal and mechanical [3] behavior. The interest is not only from a fundamental research point of view; the possibilities of using such clusters in biological [4], chemical and even cancer therapy [5, 6] applications are being actively discussed. Consequently a solid understanding of the dynamical (optical) response of C60 to external excitations such as lasers, radiation, or particles is required. With the discovery of other similar systems as C20, C180, C240, the fullerenes became an important domain of study while C60 cluster still remaining the main figure in the class.

One of the most intriguing and important properties of C60, also from the point of view of possible applications, is related to the presence of a collective giant resonance situated at an unusually high energy in the optical spectrum. The first theoretical calculations within the linear response theory [7] predicted a position of the centroid around 20 eV. Later,

Author to whom any correspondence should be addressed.

photo-ionization experiments [8] confirmed the dominance in the photo-ionization yield of a resonance observed at 20 eV with a width around 10 eV. It can be interpreted accordingly with the Mie's theory [9] as a surface plasmon, which is usually present in metal nano-structures as a dipolar vibration of delocalized electrons.

This aspect can be explained from the electronic structure 1s22s22p2 of carbon atoms, considering that the valence electrons 2s22p2 are quasi-delocalized in the cluster structure and capable of being collectively excited. Even more peculiar was the presence of a second, smaller peak in the optical spectrum at around 40 eV revealed in photo-ionization experiments [10, 11] both for neutral or ionized fullerenes. The recent studies provided different interpretations of this resonance in terms of surface and volume plasmons. For the optical response of C60 several investigations were devoted to the single particle and collective aspects as well to the role of quantum effects in the dynamics [12-16].

The considered theoretical approaches, as Hartree-Fock (HF) [17, 18], Random Phase Approximation (RPA) [7], Density Functional Theory (DFT) [11, 12, 19], can be quite involved from a computation point of view.

We intend in this paper to provide an investigation of the photo-absorption spectrum based on a semi-classical treatment, namely Time Dependent Thomas Fermi (TDTF). This

0953-4075/15/185102+08$33.00

© 2015 IOP Publishing Ltd Printed in the UK

method or other variations of the same approximation, have been employed to explore specific features of systems as the atomic nuclei [20], the atoms [21], metallic clusters [22], thin metal films [23] or the laser interaction with the clusters [24] and even transport processes in semiconductors [25].

We take the advantage of a great simplification from a computational point of view and provide a detailed analysis of the optical spectrum of C60, the agreement with experimental data being reasonable for both resonances. By including consistently additional quantum effects we also compare our predictions with those from more complex studies. In any case, to the best of our knowledge, such a TDTF model has never been used before in connection with C60 in its fully nonlinear form. In spite of the inherent limitations there are also some specific advantages of our approach:

and connected with the density of particles. This local form is used further in complex system with interactions to approximate, also locally, the true kinetic energy density. It is useful to see how TF emerges as a genuine Density Functional Theory. From Hohenberg-Kohn (HK) theorem [30] results, that given a system of particles with a two-body interaction v (r, r') placed in an external potential V (r), there is an unique energy functional of density E [p] which is minimized only by the true ground state density p0 (r) of the system:

E[p] = G[p(r)] + f3 p(r)V(7)d3r

+ 1 f 3 p(?)v(7, 7')p(r')d3r'd3r 2 *3 (1)

(i) TDTF requires the propagation of a single pseudo-wave function in time in comparison with TDDFT or TDHF which propagate 240 or even 360 electron orbitals, and therefore, the gain in computational complexity and time is obvious;

(ii) The method allows us from its first principles to include or not, different corrections which account for pure quantum effects and so, we can investigate when these effects become important in the dynamics;

(iii) The numerical scheme gives direct access to the local dynamics of electronic density, from which a quantitative description of the nature of the resonances in terms of the transition densities is possible.

Moreover, we can analyze numerically the capabilities of the semi-classical hydrodynamic theory in the case of complex systems as fullerenes with the goal to extend such investigations to similar structures, at least for aspects related to the gross properties of dynamics.

2. Theory

2.1. Jellium model and Thomas Fermi theory

In order to explore the optical response of C60 in the framework of TDTF we adopt for the core electrons and atomic nuclei a description based on jellium model. So, instead of nuclei and bound electrons we consider a homogeneous positive charge distribution which creates a fixed potential for the valence electrons. The self-consistent jellium model has proved to be an appropriate method predicting quantitative results in good agreement with the experimental data in cluster physics [26, 27]. In our case we use the classical model of a positive charged shell centered on the support surface for the carbon nuclei in C60 cluster which was considered also in other approaches [12, 28, 29].

Before presenting the TDTF method, we mention some useful features of the stationary Thomas Fermi (TF) theory. The approximation is derived considering an infinite system of non-interacting fermions. From the associated plane wave functions the local density of kinetic energy is constructed

Here G [p] is an universal (unknown) functional of the particles density. It should contain the kinetic energy of the system as well as contributions from pure quantum effects as exchange-correlation. The variational principle with the constraint to have a constant number of particles N leads to:

SE [p] Sp

where / is the chemical potential. The associated Euler-Lagrange equation becomes (with the notation SG/dp = g [p]):

g[p(r)] + v(r) + f3 v7')p(')d3r' = m (3)

The equation (3) has been derived from the H-K only and the result is exact in the frame of DFT. Now, within TF approximation the functional g [p] has the following local form: g [p] = gp2/3, with g = (3p2)2/3A2/2m. Being obtained within the homogeneous electron gas approximation, the functional works only in the case of systems with delocalized electrons and slow varying density (low gradients). The theory can be improved by adding gradient corrections [31] and exchange-correlation [32] terms in g [p]. Nevertheless the TF method provided already positive results in the study of metallic clusters [33, 34] and recently of C60 [35] and this justifies our extension to the time-dependent case.

2.2. Derivation of Time Dependent Thomas Fermi equation

The TDTF can be formulated in terms of a hydrodynamical-like system of equations which can be derived from the Vlasov semi-classical approximation of quantum Wigner equation [36]. Here we will obtain these equations from Time Dependent DFT principles by using the Lagrangian formalism. In this way we expect a more transparent picture of the chain of approximations. A short description of this approach with an application to the electron dynamics in metal clusters can be found in [22]. The main idea is to start, analogous with the stationary version (TF), from the quantum mechanical

action and the variational principle:

A(t) = f' (Y(r)\mdT — H\y(t))dt J 0

where now the state \y) is an unique functional of the local density: \y) ° \y [p (r, ')]) [37]. Here H is the Hamiltonian operator. From the variational principle the following set of hydrodinamical equations can be derived [37]:

d'p(r, ') + Vj (r, ') = 0 (4a)

dtj (r, ') = P [p]((, ') (4b)

P [p]((, ' ) = -i (y [p]\ [j , H ] \y [p]) (4c)

In the equation (4c) j is the current operator while P is obtained as the expectation value of the commutator of the current and the Hamiltonian. Unfortunately, the Hamiltonian operator is not known explicitly and we can make approximations just on its expectation values, so there is no easy way to compute the commutator.

Here we adopt, as in [38], a fluid description of the system of particles having an internal energy given by the energy functional provided by DFT. The dynamical part of the kinetic energy has the collective component defined by a velocity field u (r, '). The expectation value of the Hamiltonian operator (y\H\y) = H[p, u] is:

H [a u= L p ( t )

и (r, t)2 so [p(r, t)]

V(r, t) + J3 v(r, r')p(r', t)dr

dr3 (S)

дt p =--

дtS =

As an additional approximation we assume that the velocity field u is irrotational, i.e. u = VS. With this representation the equations of motion 6a and 6b become:

дtS + !vs2 + g [p] + v + u = 0

The above set of equations can be embedded in a single non-linear Schroedinger-like equation by introducing the Madelung transform [39], which defines the complex field F = p1/2 exp (iS/ ft ) in terms of {p, S} variables:

iftдt F =--v2f + wF

j 2 v2

= g [ IFI2 J + V + U +

ft2 v2 If I

2m IfI

In the equation (5), v(r, r') ° \r — r' \-1 is the Coulomb interaction between electrons and further we will use as notation for the electrostatic potential U (r) generated by the density distribution which is the solution of Poisson equation v2U(r, ') = \e \p (r, ')/e0. In order to obtain the equations of motion, a scalar S (r, ') as the conjugated variable to the density such th at u (r, ') ° u [S (r, ')] is introduced. Let us drop further (r, ') arguments for legibility, but remaining implied for quantities as p, S, u, etc. By performing a Legendre Transformation H [p, S] = fpdtS — £ the following Hamilton equations are derived:

Where w is an effective pseudo-potential w. Finally, by setting g [p] = gp2/3, we obtain the TDTF model for the dynamics of a system of identical particles in an external potential. From a time-dependent DFT perspective it can be observed that, within the indicated approximations, all Kohn-Sham orbitals f = \ fj \ e' in the kinetic energy

functional T [{fj}] = — y- ^ ^=i f f* V2fjd3r acquire a constant amplitude \ fj \ = cons' and have an isotropic distribution inside a sphere of radius proportional to the Fermi momentum. For the dynamic part, all the phases of fj are the same arg (fj) — arg (f) = S, Vj = 1, N and the correlation and exchange effects are neglected. Arriving to a single NLSE (8a), very efficient computational methods such as Cranck-Nicholson [40], which is known to be unconditionally stable and norm conserving, can be used to simplify the numerical task. Perhaps the most important aspect of the above presented derivation of TDTF is the fact that, even heavily simplified, it remains a genuine Density Functional Theory and so, we should expect at least a qualitative agreement with results of TDLDA and in general KS-like schemes.

2.3. Ground state and corrections

Let us first notice that the stationary solution of the equation (8a) is the same with the one obtained from the TF equation (3). Indeed, with the ansatz F (r, ') = F0 (r)e—'m'n the equation (8a) transforms into an eigenvalue problem:

— V2 + w (F0 ) 2m

which by some simple algebra can be proven to be equivalent to the Euler-Lagrange (tf) 3 equation:

m = g [i f0i2] + v + и

д{ p + V (p VS) = 0

Let us notice that starting from the Wigner-Kirkwood expansion [41] of Bloch density function, one derives TF equation of state just as the 0th order approximation which is valid only for very slow varying densities. The additional higher order (gradient) corrections could be added. This is done through a Weiszacker term A (v2^1/2)/p1/2 with A = 1/8 [42] which accounts for local oscillations of the density profile [31]. Moreover the universal functional G

contains not only the kinetic energy of the system but also the exchange-correlations term. Consequently we can correct g [p] by adding a supplementary contribution, present in the DFT potential of the Kohn-Sham equations:

-2 y2p1/2

g [p] = gp2'3 + A ---p— + vxc [p] (11)

2m p1'2

In our calculations for vx Lundqvist functional (12)

Vxc (p) = -

we choose the Gunnarson-

(( Г [

1 + 11.4—

Figure 1. The radial profile of the jellium density (step-like, dot-dashed) and of the electron density from TF (continuous) and from TFWXC (dashed) in C60.

In this expression rs depends on local density and is given by rs [p(r)] = (4pp (r)/3)-^3. Along the next sections we shall compare TF predictions with those obtained when both Wiezsacker and xc (which will be referred as TFWXC) terms are included, since from kinetic perspective, both gradient corrections and exchange-correlations effects have the same magnitude [43].

3. Results and discussion

3.1. C60 ground state

Before the analysis of the optical response of C60 we shall present the properties of the ground state obtained by solving the equation (3). The jellium structure is given by a homogeneous positively charged shell centered on the spherical surface that contains the geometric position of nuclei pjei (r) k Q [(r — r1)(r2 — r)]. This sphere has a radius of r0 = 3.54A and the thickness is taken as in [29, 44, 45] d » 1.5A. Therefore r1 » 2.8A and r2 » 4.2A.

The solution of the equation (3) can be obtained through iterations: we start with an initial guess on the ground state density p0 from which the Coulomb potential is constructed with the help of Poisson equation v2f = —ep/e0 and then from g functional, the new density and the chemical potential are calculated in such a way that the normalization condition Jp (r) d3r = 240 is fulfilled. The process continues until the expected convergence is reached.

Since in the case when the gradient and exchange-correlation corrections are included solving the equations for both p and p can be more difficult, we have chosen instead, to consider the eigenvalue problem (9). We start with a guess for f0 from which the potential w is constructed. Then, in a finite difference description on real space of the Laplacian operator, the spectrum of eigenvalues and eigenvectors is computed. The eigenvector with the lowest eigenvalue is taken as new solution and the loop is repeated until convergence. In this way, the g functional can be computed directly from the last solution and it is not necessary to solve complicated equation such as g(x) = a. In practice, to achieve convergence, a

Figure 2. The radial dependence of the total potential, w [pTFwxc ] (solid line) and of Bohm potential V2pl!2/p1/2 (dashed line).

weighted iteration has been used pnew = apold + (1 — a) pnew where a e (0, 1).

The results for the electron density distribution are plotted in figure 1, the differences between the simple TF approximation and TFWXC being quite small for this quantity. In figure 2 we also compared the effective w potential and Bohm potential corresponding to the density calculated in the TFWXC approach. Naturally, the Bohm potential is the term of the main contribution since it must balance the fake diffusion term.

We have also performed a DFT-LDA calculation solving the KS equations for the ground state orbitals in C60 based on a numerical method which is essentially the same as for the TF ground state calculations but for 240 electrons. Following [44] we have added a constant pseudo-potential V0 in the jellium region to account for the ionic structure effects. The Kohn-Sham potential can be decomposed in vKS = U + V + V0 + vxc + Vl, where vxc is the exchange correlation potential as the one used in TFWXC simulations, U is the Hartree potential given by the electrostatic interaction in the electron cloud, V is the external potential created by the jellium distribution and Vl = l (l + 1)/r2 corresponds to the centrifugal effects. For a detailed discussion about the jellium

Figure 3. We plot the differences between the radial averaged profile of electron density obtained in a present DFT-LDA calculation (pKS) and the ones obtained with TF (dashed) method (pTF) and TFWXC (solid line) version (pTFWXc).

model-LDA applicability to fullerene see [12, 46]. In figure 3(a) comparison between the electron density profile obtained from DFT-LDA calculation and the ones obtained with TF and TFWCX approximations is presented. As expected, the TFWXC is closer to DFT-LDA results, the differences being less than 10%. We extended the present comparison to other observables. The LDA obtained chemical potential is 7.54 eV while the TFWXC chemical potential obtained in this work is m » 2.37 eV. This large difference comes from the fact that in LDA the chemical potential is a single particle effect and our approximation is able to describe such quantities only on average. Concerning the static polarizability we have obtained a value a » 97A while DFT provided a value of 89a . Let us mention that the experimental value is a = 76.5 ± 8A [47]. In figure 4 we also show a comparison between the ground state KS potential for the I = 0 electrons and the same potential obtained with the PTFWXC density.

From these comparisons regarding the ground state properties we are confident that also in the dynamic regime the results can get close to those predicted by the the time dependent version of DFT, namely TDLDA.

3.2. Numerical details

From numerical point of view, we have chosen to solve both the stationary equation (8a) and the time dependent one (9) in real space based on the finite difference approximation. Since the ground state is taken as being spherically symmetric the only variable involved being the radial one, very refined grids can be adopted. For example, a radial domain of 0 < r < 3r0 divided in n = 1000-10000 equidistant points has been employed usually. A very accurate picture for both Coulomb potential and for more stiff terms as Bohm potential or the correlation potential was attained. Convergence has been achieved usually within 200-300 iterations when a criterion on the solution norm ||$0+1 — ^0II < 10-4 was imposed.

For the dynamical simulations a cylindrical coordinates system was adopted. Due to the angular symmetry only the 2D

Figure 4. The KS potential determined by the TFWXC density (dashed line) and the one obtained from DFT-LDA treatment for l = 0 orbitals (solid line).

dependence in (r, z) coordinates, where — 3r0 < z, r < 3ro, was kept. Such extended range is necessary in order to prevent spurious electron escaping from the cluster. Moreover is important to use such large distances as boundaries of the box since the Coulomb potential is obtained from Poisson equation with a set of boundary conditions extracted from a multipole expansion, valid only in the large distance approximation.

The radial coordinate has been taken to cover the whole diameter, as described in [48], in order to avoid the singularity from 0 and the implementation of boundary conditions for Poisson equation in r » 0. The usual discretization involved a number nr = nz = 100 — 200 of points on each axis. We have found that this value offers a good precision-computational cost ratio and treats properly the effects determined by the stiffness of the potential.

For the time propagation we have used a simulation interval between 10-50 fs. This interval has been discretized with a step of St = 10—3 — 10—2 fs. The numerical method adopted for the propagation of was the Crank-Nicholson (CN). This assures stability and norm conservation. Due to the self-consistency between the mean field potential w and <P, we use a mid-point rule to avoid numerical energy losses. We use w [f] to propagate a false solution f* then, we reconstruct a mid-point potential (w [f (t)] + w [f*])/2 which will propagate f (t) to f (t + St). The self consistency between the potential and the solution can induce spurious numerical variations of the norm due to non-local and nonlinear character. Still, these are negligble, having a rough magnitude of 10-3 from the total norm.

In order to study electron dynamics we induced the initial excitations through a dipolar shift in the real space which will correspond to an initial condition for <P, f (r, t = 0) = ehpF0 (r) = f0 (r — hez) , where n is a small displacement along the Oz axes and p is the corresponding boost momentum. Another possible type of excitation is through impact with fast particles (the so-called electron energy loss (EEL) experiments). This time, since the incoming particles or beams of particles have small sizes, the electric field created is no longer uniform. This can be interpreted as a local

О 0.0

1.0 t(fs)

LU 0.6

LU 0.4

Figure 5. The time evolution of the normalized dipole moment (see the text): TF (dashed lines) and TFWXC (solid line).

non-uniform displacement of electron density with effects in the volume of the cluster.

By solving the equation (8a) in the numerical frame described the time evolution of the density distribution p (r, t) = | F (r, t)|2 can be obtained. Then we can compute from our simulations the total dipole moment:

V(w) = Jeiwtd(t)dt = Jd3rdteiwtzp(7, t)

and the corresponding Fourier transform. Eventually we study the optical response through the power spectrum of the dipole moment ID(w)l2 [22, 49, 50].

The time evolution of the dipole moment is plotted in figure 5. The damping of the oscillations is larger in the presence of exchange-correlations and gradient corrections. In the present approach, the collective behavior is the only possible interpretation since the method has no single particle effects. This interpretation is actually supported by more sophisticated calculations as is TDLDA [11, 12, 15, 51].

3.3. Plasmons: existence and interpretation

In figure 6 we have plotted the power spectrum of the dipole moment and we notice the presence of two maxima at around 20 eV and 40 eV respectively which we associate with the two resonances identified experimentally [11]. Regarding the nature of the second peak, at 40 eV it was firstly assigned to a volume plasmon [11]. Other arguments [52, 53] suggested a surface plasmon also for the second resonance. Since then, more work [14-16, 51] has been devoted to achieve a clear understanding and interpretation of this collective oscillation. Based on the plasmon resonance approximation and a simple model for the fullerene shell in [15] the authors provided a decomposition of optical spectrum in three plasmons. They argued that a uniform electric field can excite only surface density variations that will determine two surface plasmons (symmetric or asymmetric coupled), while only a non-uniform electric field can induce volume density variation and consequently a volume plasmon. They connect the theoretical model with a set of electron energy loss experiments at different angles and show that indeed the spectrum is reproduced by two surface plasmons and a volume one. The 40 eV

Figure 6. The energy weighted power spectrum —w ID (—w)|2 from TF (black line), TFWXC (red dashed line) and from TDLDA (blue dot-dashed line). For TDLDA results, the single particle effects are visible in the irregular character of the spectrum. In the inset: the photo-ionization experimental spectrum (gray triangles) [11] deconvoluted in 2 lorentzians centered at 22 eV and 38 eV.

Figure 7. In radial profile, the density p (solid line), the induced Sp (t = 0) density and the limit of volume of the cluster (dotted-dashed vertical line).

plasmon observed for the first time in a photoionization experiment [11] must be interpreted as a surface one.

Still, one could argue that during the dynamics of the electron cloud the inner modulations of density might appear. In order to clarify this aspect we define as inner surface of the cluster the sphere with r » 0.75r0 and the outer surface the r » 1.25r0. This choice is somehow arbitrary but is motivated by the fact that inside this imaginary shell is contained around 85% of the total electronic charge and so can be interpreted as the volume of the cluster. Now we can see in figure 7 how the initial induced density by the dipolar shift have its peaks on what we have defined as "surface" consistent with the fact that the modulation of density by an uniform electric field is minimal in the volume of the system.

Next, to see if volume fluctuations appear over time we should look on Sp (r, t) = p (r, t) — p0 (r) at all moments. While this is possible, a straightforward image is hard to obtain and somehow ambiguous. Therefore, we intend to put Sp in connection with the spectral response, performing a local Fourier Transform and obtaining a spectral transition

1.6 1.4 1.2 1.0 0.8 0.6 0.4

0 10 20 30 40 50 ficj (eV)

Figure 8. The transition densities in the fuw-z (r0 ) space.

density:

x (r, w) = J 6p (r, t)

We plot | c (0, z, w) |2 in figure 8. As expected, the higher peak of the response is present around 20 eV and a smaller one around 40 eV, consistent with the total spectrum, located spatially on what we have defined as surfaces of the system. In the volume region, as defined above, the amplitude of density fluctuations is practically zero. We therefore conclude that the analysis based on the transition densities suggests that at 40 eV manifests a surface plasmon.

Conclusions

Starting from the fundamental theorems of DFT and TDDFT we have presented a derivation of the so-called TDTF method which uses as approximations the TF equation of state and the picture of single valued velocity field for a fermionic system. The semi-classical and exchange-correlation corrections are discussed and included for comparison in our simulations. The method is used in the case of dipolar oscillations in C60 cluster. Under spherical symmetry and jellium model approximation, the electronic ground state radial profile is calculated in good agreement with the experimental descriptions of the inner and outer surfaces.

Furthermore, in the dynamical regime, we excite the electrons by a small initial dipole shift and let the system evolve under the TDTF equation in its NLSE form. We found a good result in the optical spectrum where the resonance from 20 eV and that from 40 eV are semi-quantitatively obtained. The inclusion of corrections is found to produce

mainly more damping effects in comparison with the pure TF approximation.

In order to have a quantitative description of the second resonances we analyze the spectral transition densities from which we conclude that in C60 optical response we have two surface collective oscillations which account both for the 20 eV peak and also for the 40 eV one. The local oscillations in the volume of the system are minimal.

In summary, the TDTF(WXC) is a good choice over more complex approaches, as TDDFT is, from the perspective of computational costs, this being the main advantage of the method. On the downside, the approximations introduced through the kinetic energy functional and the assumption of irrotational velocity field reduce the accuracy of the results and make the theory capable of describing only collective phenomena. For smartly chosen systems, the disadvantages of the method are less important and the gross properties are well described. Therefore, it is advisable to use such a method for systems in which only the dynamics of the valence electrons is involved and the single particle effects do not play a major role. Finally, in contrast with DFT, TDTF(WXC) gives direct access to macroscopic quantities as density of particles or current density.

Acknowledgments

For this work V. Baran was supported by a grant of the Romanian National Authority for Scientific Research, CNCS —UEFISCDI, project number PN-II-ID-PCE-2011-3-0972.

References

Kroto H W, Heath J R, O'Brien S C, Curl R F and Smalley R E

1985 Nature 318 162-3 Richard B 1965 Laminar geodesic dome uS Patent 203 144

http://google.com/patents/US3203144 Dresselhaus M S, Dresselhaus G and Eklund P C 1996 Science of Fullerenes and Carbon Nanotubes: Their Properties and Applications (Academic Press) Bosi S, da Ros T, Spalluto G and Prato M 2003 European

Journal of Medicinal Chemistry 38 913-23 Mroz P, Pawlak A, Satti M, Lee H, Wharton T, Gali H, Sarna T and Hamblin M R 2007 Free Radical Biology and Medicine 43 711-9 Mroz P, Tegos G P, Gali H, Wharton T, Sarna T and Hamblin M R 2007 Photochemical & Photobiological Sciences 6 1139-49 Bertsch G F, Bulgac A, Tomanek D and Wang Y 1991 Phys.

Rev. Lett. 67 2690-3 Hertel I V, Steger H, de Vries J, Weisser B, Menzel C,

Kamke B and Kamke W 1992 Phys. Rev. Lett. 68 784-7 Mie G 1908 Annalen der physik 330 377-445 Reinköster A, Korica S, Prümper G, Viefhaus J, Godehusen K, Schwarzkopf O, Mast M and Becker U 2004 J. Phys. B: At. Mol. Opt. Phys. 37 2135 Scully S W J et al 2005 Phys. Rev. Lett. 94 065503

[12] Madjet M E, Chakraborty H S, Rost J M and Manson S T 2008

J. Phys. B: At. Mol. Opt. Phys. 41 105101

[13] Belyaev A, Tiukanov A, Toropkin A, Ivanov V,

Polozkov R and Solov'Yov A 2009 Phys. Scr. 80 048121

[2] [3]

[9] [10]

[14] Verkhovtsev A V, Korol A V, Solov'yov A V, Bolognesi P,

Ruocco A and Avaldi L 2012 J. Phys. B: At. Mol. Opt. Phys. 45 141002

[15] Bolognesi P, Avaldi L, Ruocco A, Verkhovtsev A,

Korol A and Solovyov A 2012 The European Physical Journal D- Atomic, Molecular, Optical and Plasma Physics 66 1-9

[16] Verkhovtsev A, Korol A and Solov'yov A 2013 Plasmon

excitations in photo-and electron impact ionization of fullerenes Journal of Physics: Conference Series vol 438 (IOP Publishing) p 012011

[17] Sheka E F 2007 Int. J. Quantum Chem. 107 2803-16

[18] Talapatra G B, Manickam N, Samoc M, Orczyk M E,

Karna S P and Prasad P N 1992 The Journal of Physical Chemistry 96 5206-8

[19] Bauernschmitt R, Ahlrichs R, Hennrich F H and Kappes M M

1998 J. Am. Chem. Soc. 120 5052-9

[20] Holzwarth G and Eckart G 1979 Nucl. Phys. A 325 1-30

[21] Horbatsch M and Dreizler R 1981 Zeitschrift für Physik A

Atoms and Nuclei 300 119-27

[22] Domps A, Reinhard P G and Suraud E 1998 Phys. Rev. Lett. 80

5520-3

[23] Crouseilles N, Hervieux P A and Manfredi G 2008 Phys. Rev.

B 78 155412

[24] Fennel T, Bertsch G and Meiwes-Broer K H 2004 The

European Physical Journal D-Atomic, Molecular, Optical and Plasma Physics 29 367-78

[25] Jüngel A 2001 Quasi-hydrodynamic Semiconductor Equations

vol 41 ((Springer)

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

[27] de Heer W A 1993 Rev. Mod. Phys. 65 611

[28] Yannouleas C and Landman U 1994 Chem. Phys. Lett. 217

175-85

[29] Rüdel A, Hentges R, Becker U, Chakraborty H S,

Madjet M E and Rost J M 2002 Phys. Rev. Lett. 89 125503

[30] Hohenberg P and Kohn W 1964 Phys. Rev. 136 B864

[31] Yang W 1986 Phys. Rev. A 34 4575

[32] Hodges C 1973 Can. J. Phys. 51 1428-37

[33] Kresin V 1988 Phys. Rev. B 38 3741

[34] Serra L, Garcías F, Barranco M, Navarro J, Balbás L,

Rubio A and Mañanes A 1989 J. Phys.: Condens. Matter. 1 10391

[35] Palade D and Baran V 2014 Romanian Journal of Physics (In

press) arXiv:1406.3826

[36] Manfredi G 2005 Fields Inst. Commun. 46 263-87

[37] Runge E and Gross E K 1984 Phys. Rev. Lett. 52 997

[38] Lundqvist S and March N H 1983 Theory of the

Inhomogeneous Electron Gas (New York: Plenum Press)

[39] Madelung E 1927 Zeitschrift für Physik A Hadrons and Nuclei

40 322-6

[40] Taha T R and Ablowitz M I 1984 J. Comput. Phys. 55 203-30

[41] Fujiwara Y, Osborn T and Wilk S 1982 Phys. Rev. A 25 14

[42] Weizsäcker C v 1935 Zeitschrift für Physik A Hadrons and

Nuclei 96 431-58

[43] Bonitz M and Dufty J 2004 Cond. Matt. Phys 7 483

[44] Puska M and Nieminen R 1993 Phys. Rev. A 47 1181

[45] Weaver J, Martins J L, Komeda T, Chen Y, Ohno T, Kroll G,

Troullier N, Haufler R and Smalley R 1991 Phys. Rev. Lett. 66 1741

[46] Verkhovtsev A, Polozkov R, Ivanov V, Korol A and

Solov'yov A 2012 Journal of Physics B: Atomic Molecular and Optical Physics 45 215101

[47] Antoine R, Rayane D, Benichou E, Dugourd P and Broyer M

2000 The European Physical Journal D-Atomic, Molecular, Optical and Plasma Physics 12 147-51

[48] Mohseni K and Colonius T 2000 J. Comput. Phys. 157 787-95

[49] Calvayrac F, Reinhard P and Suraud E 1997 Ann. Phys., NY

255 125-62

[50] Calvayrac F, Reinhard P G, Suraud E and Ullrich C 2000 Phys.

Rep. 337 493-578

[51] Verkhovtsev A V, Korol A V and Solov'yov A V 2013 Phys.

Rev. A 88 043201

[52] Korol A and SolovYov A 2007 Phys. Rev. Lett. 98 179601

[53] Scully S W J etal 2007 Phys. Rev. Lett. 98 179602