A quantum relaxation-time approximation for finite fermion systems
AA quantum relaxation-time approximation for finite fermion systems
P.-G. Reinhard ∗ and E. Suraud
2, 3, 4 Institut f¨ur Theoretische Physik, Universit¨at Erlangen, Staudtstraße 7, D-91058 Erlangen, Germany Universit´e de Toulouse; UPS; Laboratoire de Physique Th´eorique, IRSAMC; F-31062 Toulouse Cedex, France CNRS; UMR5152; F-31062 Toulouse Cedex, France Physics Department, University at Buffalo, The State University New York, Buffalo, NY 14260, USA (Dated: February 15, 2018)We propose a relaxation time approximation for the description of the dynamics of strongly excitedfermion systems. Our approach is based on time-dependent density functional theory at the levelof the local density approximation. This mean-field picture is augmented by collisional correlationshandled in relaxation time approximation which is inspired from the corresponding semi-classicalpicture. The method involves the estimate of microscopic relaxation rates/times which is presentlytaken from the well established semi-classical experience. The relaxation time approximation impliesevaluation of the instantaneous equilibrium state towards which the dynamical state is progressivelydriven at the pace of the microscopic relaxation time. As test case, we consider Na clusters of varioussizes excited either by a swift ion projectile or by a short and intense laser pulse, driven in variousdynamical regimes ranging from linear to strongly non-linear reactions. We observe a strong effectof dissipation on sensitive observables such as net ionization and angular distributions of emittedelectrons. The effect is especially large for moderate excitations where typical relaxation/dissipationtime scales efficiently compete with ionization for dissipating the available excitation energy. Tech-nical details on the actual procedure to implement a working recipe of such a quantum relaxationapproximation are given in appendices for completeness.
PACS numbers: 05.30.Fk,31.70.Hq,34.10.+x,36.40.Cg
I. INTRODUCTION
The analysis of the non linear response of finite fermionsystems subject to strong perturbations constitutes acentral issue in many areas of physics. Prominent ex-amples are low energy nuclear physics and fission [1, 2]as well as excitation of clusters and molecules by in-tense laser pulses [3, 4]. But similar situations are alsoencountered elsewhere, for example, in trapped Fermigases [5] or electron transport in nano systems [6]. Thecase of irradiated clusters and molecules has become par-ticularly interesting due to recent progress in experi-mental techniques at the side of photon sources [4, 7]as well as at the side of detection accessing more andmore detailed information from the decaying system.Particularly useful observables are from the distribu-tions of emitted electron (energy, angular distributions...)through elaborate imaging techniques such as VelocityMap Imaging (VMI) [8, 9]. These highly sophisticatedmeasurements call for elaborate theoretical modeling toreconstruct the underlying dynamics of both the irradi-ation and de-excitation process. The choice of modelsranges from a highly detailed quantum description tomacroscopic rate equations [4]. Each one of these ap-proaches is valid in a limited range of dynamical scenar-ios. The robust and rather versatile Time-DependentDensity Functional Theory (TDDFT), mostly realizedat the level of the Time-Dependent Local-Density Ap- ∗ Corresponding author : [email protected] proximation (TDLDA), [10–12] certainly provides one ofthe best compromises in the domain of quantum dynam-ics from the linear regime up to highly non-linear pro-cesses [4, 13, 14]. Still TDLDA basically remains a mean-field approach and misses by construction dissipative ef-fects from electron-electron collisions which are expectedto play a role in the course of violent excitation/de-excitation scenarios. This limits the range of applica-bility of TDLDA to moderate excitations and/or the en-trance phase of the dynamics. There is thus a strong in-terest to extend TDLDA by collisional dynamics in orderto push the limits farther out. It is the goal of this paperto investigate such an extension by including the effect ofcollisional correlations (or dynamical correlations). Thenotion of “collisional correlations” is taken from Fermiliquid theory [15] with incoherent reduction of two-bodycorrelations to two-Fermion collisions. As first example ofapplication, we discuss the case of irradiated metal clus-ters. But we emphasize that the strategy outlined here isapplicable to any many-fermion system where TDLDA isa good starting point, namely where it provides a gooddescription of low energy dynamics.The much celebrated Boltzmann equation constitutesthe prototype approach to collisional dynamics in classi-cal systems [16]. In quantum systems, one has to accountfor the Heisenberg uncertainty relation and the Pauliprinciple. Pauli blocking can be accounted for by ex-tending the Boltzmann collision term to the Boltzmann-Uehling-Uhlenbeck (BUU) form [17]. This semi-classicalBUU approach (also known as Vlasov-Uehling-Uhlenbeck(VUU) equation) provides an acceptable picture at suf-ficiently large excitations where quantum shell effectscan be ignored. It has been extensively used in nu- a r X i v : . [ phy s i c s . a t m - c l u s ] S e p clear physics [18, 19] and also explored in metal clusters[20, 21] in a high excitation domain. A further limita-tion of the BUU/VUU approach appears in the case ofclusters: it is not clear that it could be used in systemsothers than simple alkalines, where electron wave func-tions are sufficiently delocalised and smooth to allow asemi-classical treatment [22]. This hinders, e.g., an ap-plication to C which is one of the systems attractingthe most elaborate analysis of dissipative dynamics sofar [7, 23]. It should finally be noted that even in thehigh-excitation domain a major de-excitation channel isionization which may quickly take away large amounts ofexcitation energy cooling the system down into a regimewhere quantum effects cannot be neglected any more.This limits BUU/VUU often to the initial stages of a dy-namical process. All this shows that there is an urgentneed for a quantum description augmented by relaxationeffects.In spite of their limitations, semi-classical approachesnevertheless provide extremely useful guidance in threemajor aspects. First, they constitute a long standingtesting basis for strategies to implement approximate re-laxation pictures into a quantum framework [24]. Second,as we shall exploit below, they serve as source of inspira-tion for developing simple approximations for collisionalrelaxation in TDLDA. Last, but not least, they deliverapproved estimates for the key quantities (collision rate,intrinsic relaxation time) to be used in quantum approx-imations.The aim of this paper is to propose an extension of(real-time) TDLDA by collisional correlations. We shallexploit the experience from semi-classical approximationsbut keep the level of description quantum mechanicalthroughout. The paper is organized as follows. SectionII provides the basics of the mean field description andintroduces associated notations. Section III describes theproposed scheme for the Relaxation Time Approximation(RTA) and details of its handling. First results for thecase of metal clusters are then presented in section IVand conclusions and perspectives are drawn in section V.Finally, the text is completed by a series of appendiceswhich give some technical details on the RTA and itshandling in practice. II. MEAN-FIELD PROPAGATION
The starting point and dominant feature of the dynam-ics is the propagation at the level of the mean field. Inthis paper, we are dealing with the electron dynamics inmetal clusters and we describe it by time-dependent den-sity functional theory at the level of the Time-DependentLocal-Density Approximation (TDLDA) treated in thereal time domain [10, 11]. It is augmented by a self-interaction correction (SIC) approximated by average-density SIC (ADSIC) [25] in order to attain correct ion-ization properties [26] in the course of the dynamical sim-ulation. TDLDA is formulated within the usual Kohn- Sham picture in terms of a set of occupied single-particle(s.p.) wavefunctions {| φ α (cid:105) , α = 1 ...N } . Their dynamicsis described by the time-dependent Kohn-Sham equationi ∂ t | φ α (cid:105) = ˆ h [ (cid:37) ] | φ α (cid:105) (1)where ˆ h is the Kohn-Sham mean-field Hamiltonian whichis a functional of the instantaneous local density (cid:37) ( r , t ) = (cid:80) α | φ α ( r , t ) | [27, 28]. The time evolution delivered byEq. (1) can be expressed formally by the unitary one-body time-evolution operatorˆ U ( t, t (cid:48) ) = ˆ T exp (cid:32) − i (cid:90) t (cid:48) t ˆ h ( t (cid:48)(cid:48) ) dt (cid:48)(cid:48) (cid:33) (2a)where ˆ T is the time-ordering operator. This yields aclosed expression for the time-evolution of s.p. states | φ α ( t ) (cid:105) = ˆ U ( t, t (cid:48) ) | φ α ( t (cid:48) ) (cid:105) . (2b)So far, TDLDA propagates pure states. Dissipationwhich we will add later on leads inevitably to mixedstates. This requires to generalize the description fromfully occupied s.p. wavefunctions to a one-body densityoperator ˆ ρ . It is often denoted as one-body density ma-trix ρ ( r , r (cid:48) ) which is, in fact, the coordinate represen-tation of ˆ ρ . A (natural orbitals) representation of theone-body density operator in terms of s.p. wavefunctionsreads ˆ ρ = ∞ (cid:88) α =1 | φ α (cid:105) W α (cid:104) φ α | . (3)New are here the weights W α , the probability with whicha state | φ α (cid:105) is occupied. The mean-field propagation (1)then becomes i ∂ t ˆ ρ = (cid:104) ˆ h [ (cid:37) ] , ˆ ρ (cid:105) (4)where ˆ h [ (cid:37) ] is formally the same as before and the localdensity is now computed as (cid:37) ( r , t ) = (cid:88) α W α | φ α ( r , t ) | . (5)The pure mean-field propagation (4) leaves the occupa-tion weights W α unchanged and propagates only the s.p.states. The mean-field propagation of an initial state (3)then reads ˆ ρ ( t ) = ∞ (cid:88) α =1 | φ α ( t ) (cid:105) W α (cid:104) φ α ( t ) | = ˆ U ( t, ρ (0) ˆ U − ( t,
0) (6)where ˆ U is the mean-field evolution operator (2a). III. THE RELAXATION-TIMEAPPROXIMATION (RTA) FOR FINITEFERMION SYSTEMSA. Motivation: RTA in a semi-classical framework
In homogeneous fermion systems the phase space dis-tribution f only depends on momentum and dynamicalcorrelations can be described by the Uehling-Uhlenbeckcollision term I UU [ f ( p )] [17, 29]. It is a functional ofthe momentum space distribution f ( p ) which drives thedynamics steadily towards the thermal equilibrium distri-bution f eq . As the collision term I UU conserves particlenumber, mean momentum and energy, the f eq ( p ; (cid:37), j , E )represents the equilibrium for given density (cid:37) , current j ,and kinetic energy E . Sufficiently close to equilibrium,one can approximate the convergence as exponential re-laxation which allows to model the dynamical processsimply as ∂ t f ( p , t ) = − τ − ( f ( p , t ) − f eq ( p ; (cid:37), j , E )) (7)where τ relax is the relaxation time. The expectation val-ues (cid:37), j , E for density, current, and energy are the onesassociated to f ( p , t ) and computed as (cid:37) = (cid:82) d p f ( p , t ), j = (cid:82) d p p f ( p , t ), and E = (cid:82) d p p f ( p , t ) / (2 m ). Thisis called the Relaxation Time Approximation (RTA). Itwas introduced in [30] and it has been used used in thatform for a wide variety of homogeneous systems [29, 31].Finite fermion systems are spatially inhomogeneous.An obvious generalization to this case is to extend f ( p ) toa phase-space distribution f ( r , p ). This leads to a semi-classical description which is valid for high excitations (ortemperatures) where quantum shell effects are obsolete.The mean-field dynamics is then described by the Vlasovequation and the dynamical correlations by an additionalcollision term I UU [ f ( r , p )] yielding together the Vlasov-Uehling-Uhlenbeck (VUU) equation [32] ∂ t f − { h, f } = I UU [ f ( r , p )] (8)where h is the (classical) mean field Hamiltonian. In thissemi-classical approach, collisions are local, changing fora given r only the momentum distribution at this point.It thus establishes local conservation laws [33] such thatcollisional relaxation conserves local density (cid:37) ( r , t ), localcurrent j ( r , t ), and local kinetic energy E kin ( r , t ). Thecollision term thus drives towards a local and instanta-neous equilibrium f eq ( r , p ; (cid:37), j , E kin ). The global equili-bration is achieved at slower pace by interplay with thelong range transport described by the mean-field propa-gation (Vlasov part of the VUU equation). The RTA forthe VUU equation (8) reads [30, 34–36] ∂ t f − { h, f } = − τ relax ( f ( r , p , t ) − f eq ( r , p ; (cid:37), j , E kin ))(9)where the constraints (cid:37), j , E kin depend on position r andtime t . This is the model which we will now generalizeto the case of quantum mean-field theory. B. RTA in quantum-mechanical framework
The generalization of the one-body phase-space distri-bution f ( r , p ) to a quantum-mechanical mean-field the-ory is the one-body density operator ˆ ρ , or one-body den-sity matrix ρ ( r , r (cid:48) ) respectively. The equation of motionfor ˆ ρ including dynamical correlations reads in general[37, 38] i ∂ t ˆ ρ − (cid:2) ˆ h, ˆ ρ (cid:3) = ˆ I [ˆ ρ ] . (10)The left hand side embraces the mean-field propagation.It may be time-dependent Hartree-Fock or the widelyused LDA version of TDDFT. The right-hand side con-sists in the quantum-mechanical collision term. Moti-vated by the successful semi-classical RTA, we importEq. (9) for the quantum case as ∂ t ˆ ρ = − i (cid:2) ˆ h, ˆ ρ (cid:3) − τ relax (ˆ ρ − ˆ ρ eq [ (cid:37), j , E ]) , (11)where ˆ ρ eq is the density operator of the thermal equilib-rium for local density (cid:37) ( r , t ), current distribution j ( r , t )and total energy E ( t ) given at that instant of time t .These constraining conditions are, in fact, functionals ofthe actual state ˆ ρ , i.e. (cid:37) [ˆ ρ ], j [ˆ ρ ], and E [ˆ ρ ]. For the diago-nal representation Eq.(3) of the density operator ˆ ρ , theyread (cid:37) ( r ) = (cid:88) α | φ α ( r ) | W α , (12a) j ( r ) = (cid:88) α W α φ ∗ α ( r ) → ∇ − ← ∇ φ α ( r ) . (12b)The energy E ( t ) is taken as the total energy becausethe semi-classical concept of a local kinetic energy is am-biguous in a quantum system. This RTA equation (11)looks innocent, but is very involved because many entriesdepend in various ways on the actual state ˆ ρ ( t ). Theself-consistent mean field is a functional of the actuallocal density, i.e. ˆ h = ˆ h [ (cid:37) ]. The instantaneous equilib-rium density ˆ ρ eq is the solution of the stationary, thermalmean-field equations with constraint on the actual (cid:37) ( r ), j ( r ) and energy E , for details see Appendix B.The relaxation time τ relax is estimated in semi-classicalFermi liquid theory, for details see appendix A 2. For themetal clusters serving as test examples in the following,it becomes (cid:126) τ relax = 0 . σ ee r s E ∗ intr N , (13)where E ∗ intr is the intrinsic (thermal) energy of the sys-tem (appendix C), N the actual number of particles, σ ee the in-medium electron-electron cross section, and r s theeffective Wigner-Seitz radius of the electron cloud. C. Summary of the procedure
The solution of the RTA equations is rather involved.We explain the necessary steps here from a practical side starting point: ˆ ρ ( t ) = (cid:80) α | φ α ( t ) (cid:105) W α ( t ) (cid:104) φ α ( t ) | (cid:105) (cid:63) mean-field propagation: | φ (mf) α (cid:105) = ˆ U ( t + ∆ t, t ) | φ α ( t ) (cid:105) W α ( t ) = W ( mf ) α = const.ˆ ρ mf = ˆ ρ mf ( t + ∆ t ) = (cid:80) α | φ ( mf ) α (cid:105) W ( mf ) α (cid:104) φ ( mf ) α | (cid:45) (cid:37) mf ( r , t + ∆ t ) , j mf ( r , t + ∆ t ) , E mf (cid:105) (cid:105) (cid:63) density-constrained mean field (DCMF)1. ˆ ρ eq = ˆ ρ eq [ (cid:37) mf ( r ) , j mf ( r ) , E mf ]= (cid:80) α | φ (cid:48) α (cid:105) W (cid:48) α (cid:104) φ (cid:48) α |
2. intrinsic excitation energy E ∗ intr (cid:1)(cid:1)(cid:1)(cid:1)(cid:11) relaxation time: (cid:126) τ − = 0 . σ ee r − s E ∗ intr /N (cid:63) (cid:16)(cid:16)(cid:16)(cid:16)(cid:16)(cid:16)(cid:16)(cid:16)(cid:16)(cid:16)(cid:16)(cid:16)(cid:16) (cid:105) (cid:105) (cid:105) (cid:63) ˆ ρ ( t +∆ t ) = ˆ ρ mf − ∆ tτ relax [ˆ ρ mf − ˆ ρ eq ]diagonalize to natural orbitals:ˆ ρ ( t +∆ t ) = (cid:80) α | φ α ( t +∆ t ) (cid:105) ˜ W α (cid:104) φ α ( t +∆ t ) | final fine-tuning of W α to reproduce E mf ˆ ρ ( t +∆ t ) = (cid:80) α | φ α ( t +∆ t ) (cid:105) W α ( t +∆ t ) (cid:104) φ α ( t +∆ t ) | (cid:72)(cid:72)(cid:89)(cid:63) (cid:105) (cid:63) (cid:105) Figure 1: Sketch of the scheme for performing one large time step t −→ t +∆ t in solving the RTA equations. The numbers inopen circles indicate the steps as outlined in the text. and unfold details in the appendices. We briefly summa-rize the actual scheme for one step from t to t +∆ t . Notethat mean-field propagation (actually TDLDA) runs ata much faster pace than relaxation. We resolve it bystandard techniques [13, 28] on a time step δt which ismuch smaller (factor 10–100) than the RTA step ∆ t . Wesummarize this TDLDA propagation in the evolution op-erator ˆ U from Eq. (2a) and discuss only one RTA step.Its sub-steps are sketched in Figure 1 and explained inthe following whereby the label here correspond to theones in the Figure:1. We first propagate ˆ ρ by pure TDLDA. This meansthat the s.p. states in representation (3) evolve as | φ α ( t ) (cid:105) → | φ (mf) α (cid:105) = ˆ U ( t + ∆ t, t ) | φ α ( t ) (cid:105) , while theoccupation weights W α are kept frozen (pure mean-field propagation).2. We compute density (cid:37) ( r , t +∆ t ), current j ( r , t +∆ t ),and total energy E mf associated to the TDLDA-propagated density matrix ˆ ρ mf .3. We determine the thermal mean-field equilibriumstate ˆ ρ eq constrained to the given (cid:37) , j , and E mf .This is achieved by Density-Constrained MeanField (DCMF) iterations as outlined in AppendixB. The actual equilibrium state ˆ ρ eq is representedby new s.p. states {| φ (cid:48) α (cid:105)} and new occupation num-bers W (cid:48) α in diagonal form (3). 4. We compose the new density matrix from theTDLDA propagated state ˆ ρ mf and the equilibra-tion driving term ˆ ρ mf − ˆ ρ eq with the appropriateweight ∆ t/τ relax , as outlined in Appendix D. Therelaxation time Eq. (13) requires the actual intrin-sic excitation energy E ∗ intr which is also obtainedfrom DCMF, see appendix C.5. We diagonalize the state emerging from step 4 tonatural-orbital representation Eq. (3). This yieldsthe s.p. states {| φ α ( t +∆ t ) (cid:105)} for the next step andpreliminary new occupations ˜ W α .6. After all these steps, the initial energy E mf = E TDLDA ( t ) may not be exactly reproduced. Wemay remain with a small energy mismatch as com-pared to the goal E mf . We now apply a small iter-ative thermalization step to readjust the energy, asoutlined in Appendix E. This then yields the finaloccupation weights W α ( t +∆ t ) which comply withenergy conservation.The scheme can be used also in connection with absorb-ing boundary conditions [13, 39]. The particle loss willbe mapped automatically to loss of occupation weightsin step 4. A word is in order about the choice of the timesteps. The δt for propagation of TDLDA is limited by themaximal energy on the grid representation and thus verysmall (for Na clusters typically 0.005 fs). The steppingfor the relaxation term needs only to resolve the changesin the actual mean field which is achieved already with∆ t ≈ . t and find thesame results for all ∆ t ≤ . t = 0 . D. Numerical representation and computation ofrelevant observables
The numerical implementation of TDLDA is done instandard manner [13, 28]. The coupling to the ions ismediated by soft local pseudopotentials [40]. The Kohn-Sham potential is handled in the Cylindrically Aver-aged Pseudo-potential Scheme (CAPS) [41, 42], whichhas proven to be an efficient and reliable approximationfor metal clusters close to axial symmetry. Wavefunc-tions and fields are thus represented on a 2D cylindricalgrid in coordinate space [43]. For the typical example ofthe Na cluster, the numerical box extends up to 104a in radial direction and 208 a along the z -axis, whilethe grid spacing is 0.8 a . To solve the (time-dependent)Kohn-Sham equations (1) we use time-splitting for timepropagation [44] and accelerated gradient iterations forthe stationary solution [45]. The Coulomb field is com-puted with successive over-relaxation [43]. We use ab-sorbing boundary conditions [13, 39], which gently ab-sorb all outgoing electron flow reaching the bounds ofthe grid and thus prevent artifacts from reflection backinto the reaction zone. We take the exchange-correlationenergy functional from Perdew and Wang [46].A great manifold of observables can be deduced fromthe ˆ ρ ( t ) thus obtained. We will consider in the fol-lowing the dipole signal, dipole spectrum, ionization,angular distribution of emitted electrons, and entropy.We focus here on the dipole moment along symme-try axis z , which is obtained from the local density as (cid:104) ˆ d z (cid:105) ( t ) = (cid:82) d r d z (z) (cid:37) (r) where d z (z) = z is the (lo-cal) dipole operator. The dipole strength distribution iscomputed with the methods of spectral analysis [47]. Itis attained by an instantaneous dipole-boost excitation,collecting (cid:104) ˆ d z (cid:105) ( t ) during propagation, and finally Fouriertransforming (cid:104) ˆ d z (cid:105) ( t ) into frequency domain. The angulardistribution of emitted electrons is obtained from record-ing the absorbed electrons as in TDLDA [48, 49]. The an-gular distribution is characterized by the anisotropy pa-rameter β , the leading parameter in the photo-electronangular cross section dσ/d Ω ∝ (1 + β P ( cos ( θ ) + .... )[50, 51] where P is the second order Legendre polyno-mial and θ the direction with respect to laser polarizationaxis (here z -axis in 2D cylindrical geometry). A specificquantity to track relaxation processes is the one-body en-tropy which is computed in diagonal representation (3) by the standard expression [52] S = − (cid:88) α [ W α log W α + (1 − W α ) log(1 − W α )] (14)in units of Bolzmann constant. It serves as a direct indi-cator of thermalization and allows to read off the typicaltime scale of relaxation processes. IV. RESULTSA. The test cases
As test cases, we will consider the clusters, Na , Na +9 ,Na +41 , and Na +93 . All test cases have electronic shell clo-sures ( N el = 8 , ,
92 [53]) and are thus close to sphericalsymmetry. This is no principle restriction because com-putations with deformed systems show similar pattern.In fact, shell closures with their large HOMO-LUMOgaps are the most demanding situations (thus critical testcases) for the RTA scheme. The ionic ground-state con-figuration is optimized by iterative cooling in the spirit ofsimulated annealing [13, 28]. In the following we are in-terested exclusively in electronic dissipation and we areconsidering rather short time intervals. We thus keepthe ions frozen at their ground-state configurations. TheWigner-Seitz radius required in the estimate of local re-laxation time Eq. (13) is computed with the recipes ofAppendix A 2. It turns out to be almost the same for alltest cases mentioned above. We use in practice r s = 3 . . The TDLDA equations are propagated with a timestep of δt = 0 .
005 fs. The larger time step for evaluationof the dissipative term is ∆ t = 0 . B. Trends for boost excitations
In the first round, we use the simplest excitationmechanism to elucidate the basic effects of dissipativeterm. This is an instantaneous dipole boost φ α → exp( − i p ˆ d z ) φ α applied to all s.p. wavefunctions in thesame manner [13, 28]. The boost momentum p reg-ulates its strength. We are characterizing the booststrength henceforth in terms of the initial excitation en-ergy E ∗ = N p / (2 m ) brought into the cluster by theboost. Mind that E ∗ is the initial excitation energy de-posited into the system, not to be confused with the timedependent intrinsic excitation energy E ∗ intr used in es-timating the relaxation time (13). The instantaneousdipole boost models to a good approximation the time-dependent Coulomb field at the cluster site for collisionswith very fast ion passing by the cluster [54, 55].
1. Optical response as initial example
Optical response is the basic observable characterizingthe reaction of a cluster to an electromagnetic perturba- d i po l e s t r e ng t h [ a r b . u . ] excitation energy [eV]relaxno relax Figure 2: (Color online) Spectral distribution of dipolestrength for Na with ionic background in the CAPS, eval-uated for a boost strength E ∗ = 2 . tion and it serves as key to the analysis of a large varietyof dynamical scenarios [28]. It is thus of interest to checkthe impact of RTA on optical response. We take hereNa as test example. Its optical response is dominatedby a surface plasmon resonance in a rather narrow spec-tral range around 2.7 eV.Figure 2 shows the effect of dissipation on the spec-tral distribution of dipole strength for an excitation inthe one-plasmon regime, i.e. E ∗ = 2 . E ∗ .
2. Trends with varied excitation strength
Figure 3 shows the time evolution of relative ioniza-tion N esc /E ∗ , envelope of the dipole signal (cid:104) ˆ d z (cid:105) , entropy S , and anisotropy β for Na after boosts with differ-ent strengths. Results from pure TDLDA are shownwith dotted lines and those including dissipation with fulllines. The three cases for E ∗ represent three regimes: lin-ear response for E ∗ = 0 .
27 eV, the one-plasmon regimefor E ∗ = 2 . E ∗ = 8 . N esc / E * [ / e V ] time [fs]d) relax, E* =0.27 eVrelax, E* =2.7 eVrelax, E* =8.1 eVnorelax E* =0.27 eVnorelax E* =2.7 eVnorelax E* =8.1 eV 0.1 1 d i po l e a m p li t ud e [ a ] c) 0 0.2 0.4 0.6 0.8 1 e n t r op y S / S asy b) 1.1 1.2 1.3 1.4 1.5 1.6 1.7 a n i s o t r op y β a) Figure 3: (Color online) Time evolution of four key observ-ables after an instantaneous boost with three different boostenergies as indicated. Test case is Na with ionic backgroundin the CAPS. The observables are: d) ionization N esc /E ∗ rescaled to initial excitation energy E ∗ , c) envelope of thedipole signal (cid:104) ˆ d z (cid:105) , b) entropy S relative to the asymptoticvalue S asy , and a) anisotropy β of the angular distributionof emitted electrons. The cases with relaxation are drawnwith full lines and computed with the standard scatteringcross section σ ee = 6 . . TDLDA results are drawn withdotted lines. The dipole amplitudes (panel c) are shown in logarith-mic scale to visualize the long-time trend. They decaywith time, the faster the higher the excitation E ∗ . Eventhe envelope of the dipole signal is heavily fluctuatingsuch that it is hard to read off a global relaxation time.Moreover, already the TDLDA amplitudes are attenu-ated due to Landau damping and direct electron emission[13, 28]. The dissipative effect corresponds to the differ-ence between TDLDA and RTA propagation. Unfortu-nately, it cannot be extracted cleanly from the dipolesignal.The entropy, Eq. (14), stays at S = 0 for pure TDLDA.It is thus a selective signal for dissipative effects andonly results from the RTA calculations are relevant here.Panel b in Figure 3 shows the entropy relative to theasymptotic value to allow a direct comparison of the threecases. All three cases show a nice exponential conver-gence. The corresponding global relaxation time dependssensitively on the excitation energy E ∗ . There is prac-tically no dissipation for the faintest excitation. This isclear from the recipe Eq. (13) for the local relaxation ratewhich is ∝ E ∗ and thus disappears for E ∗ →
0. Dissipa-tion increases with increasing excitation amounting to aglobal relaxation time of about 50 fs for the one-plasmonregime E ∗ = 2 . E ∗ = 8 . N esc , accompanied by slowly decreasing dipole ampli-tude. Dissipation offers an alternative channel for damp-ing, namely internal excitation. Consequently, ionizationis much suppressed, the more the stronger dissipation isat work, i.e. practically not for E ∗ = 0 .
27 eV and in-creasing with E ∗ .Anisotropy is a well known signal of thermalization.Angular distributions become more and more isotropicthe more thermalized a system is. It is thus importantto compare TDLDA with RTA in that respect. This isdone in panel a of Figure 3. Dissipation clearly leadsto a reduction of β . It is, nevertheless, interesting tonote the counteracting effects pulling on β , one due toTDLDA, the other one to RTA. The trend is simple forpure TDLDA: the stronger the boost the larger the driveto forward/backward emission and thus the larger β .On the other hand, larger excitation enhances dissipationand reduces increasingly β . Eventually, the trend can goboth ways.Finally, remind that the analysis of angular distribu-tions is done here over finite times which means that the β shown here are not yet the asymptotic values. Theexcitation energy stored in intrinsic degrees of freedomby dissipation will be released slowly later, to a large ex-tent in terms of thermal electron emission. This adds anisotropic background of thermal electrons which will fur-ther reduce the anisotropy. A strategy to estimate thiseffect is proposed in [55, 56]. We skip this final clean-upin the present exploratory study as it is not crucial for N esc / E * [ / e V ] time [fs]d) σ ee =3.25 a σ ee =6.5 a σ ee =13 a σ ee =0 0.01 0.1 1 d i po l e a m p li t ud e [ a ] c) 1.1 1.2 1.3 1.4 1.5 1.6 a n i s o t r op y β b) 0 0.2 0.4 0.6 0.8 1 e n t r op y S / S asy a) Figure 4: (Color online) Same as Fig. 3, but for varied crosssections σ ee defining τ relax which correspond to the standard σ ee = 6 . and half as well as twice that value. The boostenergy is E ∗ = 2 . understanding RTA.
3. Trends with electron-electron cross section
The choice of the cross section σ ee for in-medium elec-tron scattering is, of course, a crucial ingredient for theestimate of the relaxation time Eq. (13). To explore thesensitivity on the choice of σ ee , we show in Figure 4 re-sults for three values of σ ee in steps of factors 2, namelyfor the standard value of σ ee = 6 . (see appendix A 2)used in most of the paper, and for 13 a (twice that crosssection), and for 3.25 a (half of it). The comparison isdone for the one-plasmon excitation regime E ∗ = 2 . σ ee are large. From the entropysignal (panel a), we can read off that the global relax-ation time shrinks here almost inversely proportional tothe cross section, from 100 fs over 50 fs down to about25 fs. One can spot the same trend in the attenuation ofthe dipole signal in panel c. Correspondingly, we see anincreasing suppression of ionization with increasing crosssection. In this case, however, the changes are more mod-erate, far slower than proportional. The same holds forthe anisotropy (panel b). Reduction of β increases with σ ee , but rather slowly.It is to be noted that the excitation of E ∗ = 2 . σ ee around the standard choice. Very little happens, ofcourse, deep in the linear regime where dissipation is neg-ligible anyway. Somewhat less sensitivity is also seen forheftier excitations (not shown here).
4. Trends with system size
The next question is how dissipative effects depend oncharge state and system size. Figure 5 shows results forthree cluster cations of different size, Na +9 , Na +41 , andNa +93 , all with the same excitation E ∗ = 2 . +41 here in Fig-ure 5 with the case of E ∗ = 2.7 eV for Na in Fig-ure 3 shows that one more charge changes very little inall respects (emission, dipole amplitude, relaxation time,anisotropy). However, system size makes a huge differ-ence. The smaller system Na +9 has much stronger dissi-pative effects (shorter relaxation time, more suppressionof ionization) while the heavier cluster Na +93 shows verysmall relaxation. A quick glance at recipe (13) for therelaxation rate explains this. The rate is proportional to E ∗ intr /N . Large electron number N thus reduces the ratewhile low N enhances it. The effect is obvious: the plas-mon energy depends only weakly on system size [57] butthe thermal effects are related to the energy per particle.The latter quantity shrinks with increasing size. As aconsequence, a one-plasmon excitation remains safely inthe linear regime for heavy clusters while small clustersexperience more thermal effects as, e.g., relaxation. N esc / E * [ / e V ] time [fs] 0.1 1 d i po l e a m p li t ud e [ a ] e n t r op y S / S asy a n i s o t r op y β Na Na Na Na norelaxNa norelaxNa norelax Figure 5: (Color online) Same as Fig. 3, but for varied clustersize comparing Na +9 , Na +41 , and Na +93 . Boost energy is in allcases E ∗ = 2 . C. Laser excitation
Having explored the basic features of RTA for in termsof the simple boost excitation, we now present a quickfirst exploration of laser excitations. The laser field isdescribed as a classical electro-magnetic wave handled inthe limit of long wavelengths. This amounts to add toTDLDA a time-dependent external dipole field U ext ( r , t ) = e r · e z E f ( t ) sin( ω las t ) , (15a) f ( t ) = sin (cid:18) π tT pulse (cid:19) θ ( t ) θ ( T pulse − t ) . (15b)The laser features therein are: the (linear) polarization e z along the symmetry axis, the peak field strength E related to laser intensity as I ∝ E , photon frequency ω las , and pulse length T pulse . Actually we use T pulse =96 fs corresponding to a Full-Width at Half Maximum(FWHM) of 48 fs for field strength.In order to track the energy balance in the process,we use here two energetic observables: first the energyabsorbed from the laser field E exc = (cid:90) t dt (cid:48) (cid:90) d r (cid:37) ( r , t (cid:48) ) ∂ t (cid:48) U ext ( r , t (cid:48) ) (16)where (cid:37) is the usual electron density, and second, theintrinsic kinetic energy E ∗ intr (Eq. (C1)), see appendix C.Figure 6 shows results for a laser intensity I =10 W/cm . Calculations were done up to time 250 fsfor a set of laser frequencies. Shown in panels b–e are theasymptotic values of the observables which are practicallyreached at this final time. Panel a shows for compari-son the dipole strength distribution obtained from spec-tral analysis after small boost (linear regime). The fre-quency dependence of ionization (panel e) maps roughlythe dipole strength distribution (compare with panel aand mind that the present scan of ω las has low resolu-tion). Relaxation suppresses ionization considerably atthe peaks of the distribution, i.e. for resonant excita-tion. Practically no dissipative effects are seen in theoff-resonant minima between the peaks and in the regionabove IP. Absence of dissipation is related to an emissionprocess faster than the relaxation time. This is obviousfor the energies above IP. These are direct emission pro-cesses which run for Na clusters at a time scale of fewfs [13, 28]. Very instructive is the different behavior foron- and off-resonant processes below IP. Resonant ex-citations are known to oscillate for a longer amount oftime [58] which, in turn, gives dissipation long time tointerfere. Off-resonant excitations, even if they involvemulti-photon processes, are confined to short times. Thephotons involved have to cooperate instantaneously.The intrinsic kinetic energy E ∗ intr shown in panel c ofFigure 6 complements the information from ionization.No differences between TDLDA and RTA are seen forthe off-resonant processes while resonant processes shifta larger part of the excitation energy to E ∗ intr , that partbeing thus lost for direct ionization. It is interestingto note that the total excitation energy E exc absorbedfrom the laser (panel d) is practically the same with andwithout relaxation term. It is the balance between in-trinsic excitation (thermalization) and direct ionizationwhich is regulated by the relaxation term. As expected,anisotropy β (panel b) is visibly lowered for all resonantprocesses where dissipation is at work. A similar effecthad been see in previous study using VUU [59]. N esc frequency ω [eV] I P e) with relaxno relax 0 1 2 3 4 5 6 7 8 E ecx [ e V ] I P d) with relaxno relax 0 1 2 3 4 5 E * i n t r [ e V ] I P c) with relaxno relax 0.8 1 1.2 1.4 1.6 1.8 2 a n i s o t r op y β I P b) with relaxno relax 0 5 10 15 20 25 30 35 d i po l e s t r e ng t h I P a) Figure 6: (Color online) Asymptotic values of four keyobservables after laser excitation with a laser pulse of in-tensity I = 10 W/cm , pulse with sin profile of duration T pulse = 96 fs (thus FWHM=48 fs), and varying frequency(drawn along x axis). Test case is Na with ionic backgroundin CAPS. The observables are: e) ionization N esc , d) excita-tion energy deposited by the laser E exc , c) intrinsic (thermal)kinetic energy E ∗ intr and b) anisotropy β of the angular dis-tribution of emitted electrons, and entropy. Compared areresults from RTA (using the standard scattering cross sec-tion 6.5 a ) and TDLDA without relaxation. The uppermostpanel a) shows for comparison the spectral distribution ofdipole strength computed from an instantaneous boost (seesection IV B 1). N esc frequency ω [eV] I P d) with relaxno relax 0 20 40 60 80 100 E ecx [ e V ] I P c) with relaxno relax 0 5 10 15 20 25 E * i n t r [ e V ] I P b) with relaxno relax 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 a n i s o t r op y β I P a) with relaxno relax Figure 7: (Color online) Asymptotic values of four keyobservables after laser excitation with a laser pulse of in-tensity I = 10 W/cm , pulse with sin profile of duration T pulse = 96 fs (thus FWHM=48 fs), and varying frequency(drawn along x axis). Test case is Na with ionic back-ground in CAPS. The observables are: d) ionization N esc ,c) excitation energy E exc absorbed from the laser, b) intrinsic(thermal) kinetic energy E ∗ intr , a) anisotropy β of the angulardistribution of emitted electrons. Compared are results fromcalculations with RTA (using the standard scattering crosssection 6.5 a ) and without. Figure 7 shows results for a higher laser intensity I = 10 W/cm . The pattern are much different fromthe previous case. Ionization (panel d) shows one broadresonance peak. This happens already for TDLDA be-cause ionization acts here as a strong dissipation mech-anism [60]. In fact, ionization is always fast in thisregime of violent excitations and thus always the domi-nant mechanism. Little is left to do for collisional relax-ation. The energy absorption (panel c) is huge and, ofcourse, shows no difference with and without relaxationterm. The effect of dissipation is seen more clearly in the intrinsic energy (panel b). The ratio between the fullcalculation and pure TDLDA amounts to about a factortwo. Mind, however, that this is a large difference on asmall quantity. In any case, only a small fraction of thetotal absorbed energy (panel c) is moved into intrinsicexcitation and 80% or 90% are eaten up for ionization.The E ∗ intr shows another feature which seems surprisingat first glance: a difference between RTA calculation andTDLDA persists in some region above the IP. This has asimple explanation. The strong ionization enhances theIP in the course of the dynamics. Thus the system seesin the average a larger IP than in the ground state andthis shifts the regime of purely direct processes to higherlaser intensity. The fact that somewhat more of the ab-sorbed energy is moved into intrinsic excitation is alsoreflected by the anisotropy (panel a) which is reduced ascompared to TDLDA up the point where E ∗ intr differs.After all, these two short examples indicate that theinterplay of laser excitation and collisional relaxation car-ries a world of interesting effects waiting to be uncovered.This calls for further studies. V. CONCLUSIONS
In this paper, we have proposed a practical way toinclude collisional correlations in finite fermion systemsat a quantum mechanical level. The issue is crucial forenergetic processes in many systems from nuclei to clus-ters and molecules. Our approach is inspired from thesemi-classical picture, but remains strictly quantum me-chanical. It relies on a Relaxation Time Approximation(RTA) of the quantum collision term. The key ingredi-ents of our RTA are the instantaneous equilibrium den-sity matrix ˆ ρ eq ( t ) and the instantaneous relaxation time τ relax ( t ). The scheme is applicable to any finite fermionsystem. It turns out to be efficient and robust in variousdynamical scenarios.As typical application we have investigated simplemetal clusters subject to a possibly strong electromag-netic perturbation (collision with by-passing ion or laserirradiation). Inclusion of dissipation through RTA pro-vides the expected behaviors: enhanced damping of os-cillations, reduced ionization, more energy transfer tointrinsic degrees of freedom (electronic thermalization),and more isotropic emission of electrons. The effectsstrongly depend on excitation conditions and strengths(boost amplitude, laser frequency or intensity) which isvery plausible. For the given collision rates used in ourRTA, we find a strong competition between collisionalrelaxation and damping through direct electron emissionwhose outcome depends sensitively on the actual dynam-ical conditions. RTA becomes particularly beneficial forresonant laser excitations. In this case, pure TDLDA isplagued by long lasting oscillations of the system whichare unphysical. RTA yields a realistic attenuation of thedipole signal.The present RTA provides a valuable extension of1mean field theories such as the LDA version of DFT forenergetic dynamical scenarios requiring a proper accountof dissipation. Although RTA is restricted to situationswhich can be modeled by one global (but instantaneous)relaxation rate it certainly provides a valuable step inthe right direction. It surely deserves further explorationend extension to other observables as, e.g., the widelyexplored photo-electron spectra. Work along that line isin progress. Acknowledgments
We thank Thomas Fennel for helpful discussions. Wethank the RRZE (Regionales Rechenzentrum Erlangen)for supplying the computing resources for the studiespresented here. We also thank Institut Universitaire deFrance, ITN network CORINF and french ANR Musesprogram for financial support.
Appendix A: Solution of the RTA equations1. Interlaced mean-field and RTA time steps
For completeness we remind the RTA equation (11) ∂ t ˆ ρ = − i (cid:2) ˆ h, ˆ ρ (cid:3) − τ relax (ˆ ρ − ˆ ρ eq [ˆ ρ ]) , Its r.h.s. contains first the term (1) driving mean-fieldpropagation and second the dissipative term. Mean-fieldpropagation covers all s.p. oscillations and runs at amuch faster pace than relaxation. We thus treat the twocontributions at different time scales. The dissipativeterm is evaluated in time steps of ∆ t (thus at discretetimes t n = n ∆ t ), while the mean-field propagation (12a)is resolved on a finer mesh δt . To make the dependenciesmore transparent, we express the instantaneous equilib-rium density more generally as ˆ ρ eq [ˆ ρ ], i.e. as functional ofgiven density operator ˆ ρ which is communicated throughthe local density, current, and energy of ˆ ρ .Formally, we express the higher resolution for themean-field propagation in terms of the evolution operator(2a) and use that to transform to the interaction picture.During the step from t n = n ∆ t to t n +1 = ( n + 1)∆ t wethus map ˜ ρ ( t ) = ˆ U − ( t, t n )ˆ ρ ( t ) ˆ U ( t, t n ) , (A1a)ˆ ρ ( t ) = ˆ U ( t, t n )˜ ρ ( t ) ˆ U − ( t, t n ) , (A1b)which turns Eq. (11) into ∂ t ˜ ρ = − τ relax (˜ ρ − ˜ ρ eq [ˆ ρ ]) . Integrating time over the interval [ t n , t n +1 ] yields˜ ρ ( t n +1 ) = ˜ ρ ( t n ) − (cid:90) t n +1 t n dt (cid:48) ˜ ρ ( t (cid:48) ) − ˜ ρ eq [ˆ ρ ( t (cid:48) )] τ relax ≈ ˜ ρ ( t n ) − ∆ tτ relax [˜ ρ ( t n +1 ) − ˜ ρ eq [ˆ ρ ( t n +1 )]]where the last step represents a perturbative evaluationof the dissipation in the time interval [ t n , t n +1 ]. Using Eq.(A1b), we transform back into the Schr¨odinger picture toˆ ρ ( t n +1 ) = ˆ ρ mf − ∆ tτ relax [ˆ ρ mf − ˆ ρ eq [ˆ ρ mf ]] , (A2a)ˆ ρ mf = ˆ U ( t n +1 , t n )ˆ ρ ( t n ) ˆ U − ( t n +1 , t n ) (A2b)with all terms now expressed at time t n +1 . This equationdelivers the stepping scheme. We first perform a mean-field propagation in standard manner from t n to t n +1 .This yields the propagated ˆ ρ mf . We then compute thecorresponding ˆ ρ eq [ˆ ρ mf ] ≡ ˆ ρ eq [ (cid:37) mf , j mf , E mf ], and use thatfinally to compose the relaxation step (A2a).
2. Estimate of the relaxation time
We refer to a semi-classical estimate of relaxationrates developed in [61] for homogeneous nuclear matter.Rescaling that from the nuclear spin-isospin degeneracyfactor g = 4 to the electronic spin degeneracy g = 2, weobtain as starting point (cid:126) τ relax = 3 . (cid:126) m σ ee k F ρ T ε (A3)where m is the electron mass, k F the Fermi momentum, ρ the matter density, ε F the Fermi energy, and σ ee the ef-fective in-medium cross section for electron-electron colli-sions. Using the standard relations for a degenerate elec-tron gas, r s = 1 . /k F , ε F = 3 . /r s , and ρ = 3 / (4 πr s )(see [62]), we obtain (cid:126) τ relax = 0 . (cid:126) m σ ee T . (A4)It is preferable to express the rate in term of the intrinsicthermal excitation energy [61, 62] E ∗ intr N = π T ε F (A5)where N is the actual particle number. A word of cau-tion is necessary here. We are working in a system withabsorbing boundary conditions to account for ionization.This means that the number of particles is a time de-pendent quantity. The value N is the one characteriz-ing the system at a given instant and thus time depen-dent. We resolve the above expression to T and expressagain ε F through r s . This yields finally Eq. (13), i.e.2 (cid:126) /τ relax = 0 . E ∗ intr σ ee / ( N r s ). The intrinsic excitationenergy E ∗ intr is a dynamical observable which has to bedetermined anew at each time step, for details see Ap-pendix C.The Wigner-Seitz radius r s is a crucial parameter inthe estimate of the relaxation time. It characterizes theaverage electron density and it is naturally given in thecase of the jellium model. For a cluster with detailed ionicbackground, there are two ways to deduce r s from thegiven cluster. One can take either the electronic diffrac-tion radius R el , diffr (box equivalent radius) as defined in[63] and identify r s = R el , diffr N − / (A6a)or one can take the ionic r.m.s. radius r ion which is oftensimpler to evaluate. The ionic structure has practicallyno surface zone and thus the ionic diffraction radius is R el , diffr = (cid:112) / r ion . We know that the electron distri-bution in metals follows closely the ionic background dueto the strong, attractive Coulomb interaction. This al-lows to identify alternatively r s = (cid:114) r ion N − / . (A6b)Both definitions yield in practice very similar results. Weuse the form (A6b) which is simpler to handle, particu-larly in case where the ionic structure is frozen. In prac-tice, it is sufficiently well described by a constant r s = 3 . .It remains to determine the effective cross section σ ee .We use here the careful evaluation of [64, 65]. They com-pute electron screening for homogeneous electron matterin Thomas-Fermi approximation, compute from that thescattering cross-section and apply a Pauli correction offactor 1/2. This yields σ ee = 6 . for the case of Naclusters at r s ≈ . . It is this value which was used asreference throughout this paper. Appendix B: Density Constrained Mean-Field(DCMF)
A key task is to determine the instantaneous equilib-rium density-operator ˆ ρ eq [ (cid:37), j , E ] in the RTA equation(11). It is the mean-field state of minimum energy underthe constraints of given local density (cid:37) , current j , andenergy E . A scheme for density constrained mean-field(DCMF) calculations was developed in [66]. We use ithere in a version extended to account also for the con-straint on current j ( r ). The density constrained mean-field Hamiltonian then readsˆ h dens . co . [ (cid:37) ] = ˆ h mf [ (cid:37) ] − (cid:90) d r λ (cid:37) ( r )ˆ (cid:37) ( r ) − (cid:90) d r λ j ( r )ˆ j ( r ) (B1) where ˆ h mf [ (cid:37) ] is the standard mean-field Hamiltonian forthe given local density (cid:37) ( r ), ˆ (cid:37) ( r ) is the operator of localdensity, ˆ j ( r ) the operator of local current, while λ (cid:37) and λ j stand for the associated Lagrange parameters. Theseare determined iteratively such that the solution of thecorresponding Kohn-Sham equations yields the wanteddensity (cid:37) ( r ) and current j ( r ) [66].As a further constraint, we want to adjust the equi-librium state to given energy E goal = E mf = E [ˆ ρ ] andparticle number N goal . N goal represents the actual num-ber of particles in the system, also denoted N in the coreof the text of the paper. Its computation is simple. Itis just given by the integral over the computational boxof the density (cid:37) ( r ) at a given RTA time step t n . Theadjustment of E goal and N goal is achieved by consideringa thermalized mean-field state ˆ ρ eq of the form Eq. (3)where each s.p. state φ α is augmented with an occupa-tion weight W (eq) α = 11 + exp (( ε α − µ ) /T ) , (B2a)( µ, T ) ↔ tr { ˆ ρ eq } = N goal , E [ˆ ρ eq ] = E goal . (B2b)Temperature T and chemical potential µ are to be ad-justed such that the wanted total particle number N goal and energy E goal is reproduced. To simplify the computa-tions, the occupation numbers are adjusted to the Fermiform (B2a) once before DCMF and once after the DCMFiterations. The final reoccupation after DCMF slightlyspoils the reproduction of density and current. But thedeviations remain small along all time propagation. Appendix C: The intrinsic excitation energy
DCMF is also used to compute the intrinsic excitationenergy E ∗ intr which is, in fact, the intrinsic kinetic energy.The total kinetic energy is the difference between the ac-tual energy E and the minimal energy at given local den-sity (cid:37) ( r ) which is the energy of the (stationary) DCMFstate for fixed (cid:37) ( r ), j = 0 and T = 0. Part of the total ki-netic energy goes into the kinetic energy of the collectiveflow represented by j ( r ). The other part is the intrinsickinetic energy E ∗ intr . We evaluate it by computing theDCMF state for given (cid:37) ( r ), j ( r ) (cid:54) = 0, but now for T = 0to find the minimum energy E DCMF [ (cid:37) ( r ) , j ( r ) , T = 0] un-der the given constraints. The difference between totalenergy and this energy is then the intrinsic kinetic energy E ∗ intr = E − E DCMF [ (cid:37) ( r ) , j ( r ) , T = 0] . (C1)This fully quantum-mechanical definition is used for esti-mating the relaxation time in Eq. (13) which is a crucialparameter in the propagation. For analyzing purposes(see section IV C), we can use the less expensive semi-classical Thomas-Fermi approximation for E kin , DCMF asdone in previous works [13, 28].3
Appendix D: Mixing of two one-body densitymatrices
The dissipative step delivers a density matrix (A2a)which is a mixˆ ρ ( t n +1 ) = (1 − η ) ˆ ρ mf + η ˆ ρ eq , η = ∆ tτ relax , ˆ ρ mf = (cid:88) α | φ (mf) α (cid:105) W α (cid:104) φ (mf) α | , ˆ ρ eq = (cid:88) α | φ (cid:48) α (cid:105) W (cid:48) α (cid:104) φ (cid:48) α | , where the | φ (mf) α (cid:105) = ˆ U ( t n +1 , t n ) | φ α ( t n ) (cid:105) constitute thebasis of TDLDA-propagated states while the | φ (cid:48) α (cid:105) are newstates from DCMF which also delivers new occupancies W (cid:48) α . We expand the composed state ˆ ρ ( t n +1 ) with respectto the set {| φ (mf) α (cid:105)} because it represents the majoritycontribution. It readsˆ ρ ( t n +1 ) = (cid:88) αβ | φ (mf) α (cid:105) ρ αβ (cid:104) φ (mf) β | ,ρ αβ = (1 − η ) δ αβ W α + η (cid:88) γ (cid:104) φ (mf) α | φ (cid:48) γ (cid:105) W (cid:48) γ (cid:104) φ (cid:48) γ | φ (mf) β (cid:105) . For the further propagation, we want to use again thediagonal representation (3). To this end, we diagonal-ize ρ αβ which finally delivers the new set | φ α ( t n +1 ) (cid:105) and W α ( t n +1 ) which will remain constant over the next RTAstep from t n +1 to t n +2 . Appendix E: Iterative correction of particle numberand energy
After mixing (A2a), one may end up with a slight de-viation from the wanted energy E goal although we havetaken care that both entries as such meet this energy (seesection B). But minimal energy violations may accumu-late during long propagation. The task is thus to correctfor a possibly occuring, small energy error δE = E (ˆ ρ mix ) − E goal (E1)where E (ˆ ρ mix ) is the actual total energy after the mixingstep. We need a scheme to correct the total energy withas little modification as possible. To this end, we mayemploy an iterative strategy. We avoid direct use of the equilibrium distribution and express a corrective step interms of the given occupation numbers W α . To this end,we take over from the equilibrium distribution the crucialproperty δW (equi) α = − (cid:18) δT ε α − µT + δµ T (cid:19) W (equi) α (cid:0) − W (equi) α (cid:1) . (E2)We take this as motivation to postulate for a generalchange δW α = ( δ ε α + δ ) W α (1 − W α ) . (E3)The δ i are tuned to desired change in particle number δN and in energy δE (where δE is determined from theenergy loss in the relaxation step). The conditions tofulfill are (cid:88) α δW α = 0 , (cid:88) α ε α δW α = δE . The energy conservation is formulated in terms of thes.p. energies. This is valid in the linear regime becauseof E (ˆ ρ + δ ˆ ρ ) = E (ˆ ρ ) + tr { ˆ hδ ˆ ρ } = E (ˆ ρ ) + (cid:88) α δW α ε α . Thus we want to fulfill (cid:88) α δW α = δN = δ δ ε , (E4a) (cid:88) α δW α ε α = δE = δ ε + δ ε , (E4b)with A = (cid:80) α w α (1 − w α ) A α . This linear system is re-solved by δ = ε δN − εδEε − ( ε ) , (E5a) δ = − εδN + 1 δEε − ( ε ) . (E5b)Thus we obtain the δ i and subsequently the δW α accord-ing to eq. (E3). Finally, it has to be checked whetherthe readjusted density matrix ˆ ρ mix , cor = (cid:80) α | φ α (cid:105) ( W α + δW α ) (cid:104) φ α | fulfills E (ˆ ρ mix , cor ) = E goal . In case, the energymatching is not good enough, the above stepping has tobe repeated. [1] T. Wada, N. Carjan, and Y. Abe, Nucl. Phys. A ,283 (1992).[2] Y. Abe, S. Ayik, P.-G. Reinhard, and E. Suraud, Phys.Rep. , 49 (1996).[3] U. Saalmann, C. Siedschlag, and J. M. Rost, Journal ofPhysics B (2006). [4] T. Fennel, , K.-H. Meiwes-Broer, J. Tiggesb¨aumker, P.-G. Reinhard, P. M. Dinh, and E. Suraud, Rev. Mod.Phys. , 1793 (2010).[5] J. Dalibard, in Proceedings of the International School ofPhysics-Enrico Fermi, Course CXL , edited by M. Ingus-cio, S. Stringari, and C. Wieman (IOS Press, Amsterdam, Nanoscale Energy Transport and Conversion:A Parallel Treatment of Electrons, Molecules, Phonons,and Photons (Oxford University Press, New York, 2005).[7] M. Kjellberg, O. Johansson, F. Jonsson, A. V. Bulgakov,C. Bordas, E. E. B. Campbell, and K. Hansen, Phys.Rev. A , 023202 (2010).[8] P. Wopperer, C. Z. Gao, T. Barillot, C. Cauchy,A. Marciniak, V. Despr´e, V. Loriot, G. Celep, C. Bor-das, F. L´epine, et al. (2014), submitted.[9] J. Pinar´e, B. Baguenard, C. Bordas, and M. Broyer, Eur.Phys. J. D , 21 (1999).[10] E. K. U. Gross and W. Kohn, Adv. Quant. Chem. ,255 (1990).[11] E. K. U. Gross, J. F. Dobson, and M. Petersilka, Top.Curr. Chem. , 81 (1996).[12] M. A. L. Marques, N. T. Maitra, F. M. S. Nogueira,E. K. U. Gross, and A. Rubio, Fundamentals of Time-Dependent Density Functional Theory (Lect. Notes inPhys. vol 837,Springer-Verlag, Berlin, 2012).[13] F. Calvayrac, P.-G. Reinhard, E. Suraud, and C. A. Ull-rich, Phys. Rep. , 493 (2000).[14] M. A. Marques, C. A. Ullrich, F. Nogueira, A. Rubio,K. Burke, and E. K. Gross,
Time Dependent DensityFunctional Theory (Springer, Berlin, 2006).[15] L. P. Kadanoff and G. Baym,
Quantum Statistical Me-chanics: Green’s Function Methods in Equilibrium andNonequilibrium Problems (Frontiers in physics, Ben-jamin, New York, 1962).[16] C. Cercignani,
The Boltzmann equation and its applica-tions (Applied Mathematical Sciences 67, Springer, NewYork, 1988).[17] E.-A. Uehling and G.-E. Uhlenbeck, Phys. Rev. , 552(1933).[18] G. F. Bertsch and S. Das Gupta, Phys. Rep. , 190(1988).[19] D. Durand, E. Suraud, and B. Tamain, Nuclear Dynam-ics in the Nucleonic Regime (Institute of Physics, Lon-don, 2000).[20] A. Domps, P.-G. Reinhard, and E. Suraud, Phys. Rev.Lett. , 5524 (1998).[21] T. Fennel, G. F. Bertsch, and K.-H. Meiwes-Broer, Eur.Phys. J. D , 367 (2004).[22] P. Ring and P. Schuck, The Nuclear Many-Body Problem (Springer, Berlin, 1980).[23] K. Hansen,
Statistical physics of nanoparticles in the gasphase (Springer Netherlands, Amsterdam, 2013).[24] F. Calvayrac, A. Domps, S. El-Gammal, C. Kohl, P.-G.Reinhard, and E. Suraud, in
Proceedings of Intern. Conf.NuCl Prag 97 (Czech. Journ. Phys, 1998), vol. 48, p. 715.[25] C. Legrand, E. Suraud, and P.-G. Reinhard, J. Phys. B , 1115 (2002).[26] P. Kl¨upfel, P. M. Dinh, P.-G. Reinhard, and E. Suraud,Phys. Rev. A , 052501 (2013).[27] R. M. Dreizler and E. K. U. Gross, Density FunctionalTheory: An Approach to the Quantum Many-Body Prob-lem (Springer-Verlag, Berlin, 1990).[28] P.-G. Reinhard and E. Suraud,
Introduction to ClusterDynamics (Wiley, New York, 2004).[29] D. Pines and P. Nozi`eres,
The Theory of Quantum Liq-uids (W A Benjamin, New York, 1966).[30] P. L. Bhatnagar, E. P. Gross, and M. Krock, Phys. Rev. , 1954 (511).[31] N. W. Ashcroft and N. D. Mermin, Solid State Physics (Saunders College, Philadelphia, 1976).[32] E. M. Lifschitz and L. P. Pitajewski,
PhysikalischeKinetik , vol. X of
Lehrbuch der Theoretischen Physik (Mir, Moscow, 1988).[33] K. G¨utter, P.-G. Reinhard, and C. Toepffer, Phys. Rev.A , 1641 (1988).[34] S. Kohler, Nucl. Phys. A , 315 (1980).[35] S. Kohler, Nucl. Phys. A , 181 (1982).[36] C. Y. Wong and K. T. R. Davies, Phys. Rev. C , 240(1983).[37] H. Reinhardt, P.-G. Reinhard, and K. Goeke, Phys. Lett.B , 177 (1985).[38] K. Goeke, P.-G. Reinhard, and H. Reinhardt, Ann. Phys.(N.Y.) , 257 (1986).[39] P.-G. Reinhard, P. D. Stevenson, D. Almehed, J. A.Maruhn, and M. R. Strayer, Phys. Rev. E , 036709(2006).[40] S. K¨ummel, M. Brack, and P.-G. Reinhard, Eur. Phys.J. D , 149 (1999).[41] B. Montag and P.-G. Reinhard, Phys. Lett. A , 380(1994).[42] B. Montag and P.-G. Reinhard, Z. Phys. D , 265(1995).[43] K. T. R. Davies and S. E. Koonin, Phys. Rev. C23 , 2042(1981).[44] M. D. Feit, J. A. Fleck, and A. Steiger, J. Comp. Phys. , 412 (1982).[45] V. Blum, G. Lauritsch, J. A. Maruhn, and P.-G. Rein-hard, J. Comp. Phys , 364 (1992).[46] J. P. Perdew and Y. Wang, Phys. Rev. B , 13244(1992).[47] F. Calvayrac, P.-G. Reinhard, and E. Suraud, Ann. Phys.(N.Y.) , 125 (1997).[48] A. Pohl, P.-G. Reinhard, and E. Suraud, Phys. Rev. A , 023202 (2004).[49] P.-G. Reinhard and E. Suraud, in Time-dependent den-sity functional theory , edited by M. A. L. Marques, C. A.Ullrich, and F. Nogueira (Springer, Berlin, 2006), vol.706 of
Lecture Notes in Physics , p. 391.[50] P. Wopperer, B. Faber, P. M. Dinh, P.-G. Reinhard, andE. Suraud, Phys. Lett. A , 39 (2010).[51] P. Wopperer, B. Faber, P. M. Dinh, P.-G. Reinhard, andE. Suraud, Phys. Rev. A , 063416 (2010).[52] L. E. Reichl, A Modern Course in Statistical Physics (Wiley, New York, 1998).[53] M. Brack, Rev. Mod. Phys. , 677 (1993).[54] M. B¨ar, B. Jakob, P.-G. Reinhard, and C. Toepffer, Phys.Rev. A , 022719 (2006).[55] P. Wopperer, P. M. Dinh, P.-G. Reinhard, and E. Suraud,Phys. Rep. to appear (2014), arXiv:1407.4965.[56] P. Wopperer, C. Z. Gao, T. Barillot, C. Cauchy,A. Marciniak, V. Despr´e, V. Loriot, G. Celep, C. Bor-das, F. L´epine, et al., preprint (2014).[57] P.-G. Reinhard, S. Weisgerber, O. Genzken, andM. Brack, Lecture Notes in Physics , 254 (1992).[58] P.-G. Reinhard, F. Calvayrac, C. Kohl, S. K¨ummel,E. Suraud, C. A. Ullrich, and M. Brack, Eur. Phys. J.D , 111 (1999).[59] E. Giglio, P.-G. Reinhard, and E. Suraud, Phys. Rev. A , 043202 (2003).[60] C. A. Ullrich, P.-G. Reinhard, and E. Suraud, Phys. Rev.A , 1938 (1998).[61] G. Bertsch, Z. Phys. A , 103 (1978).[62] J. Maruhn, P.-G. Reinhard, and E. Suraud, Simple mod- els of many-fermions systems (Springer, Berlin, 2010).[63] J. Friedrich and N. V¨ogler, Nucl. Phys. A , 192(1982).[64] J. K¨ohn, R. Redmer, K.-H. Meiwes-Broer, and T. Fennel,Phys. Rev. A , 033202 (2008). [65] J. K¨ohn, R. Redmer, and T. Fennel, New J. Phys. ,055011 (2012).[66] R. Cusson, P.-G. Reinhard, J. Maruhn, W. Greiner, andM. Strayer, Z. Phys. A320