On the inclusion of dissipation on top of mean-field approaches

We discuss extensions of time-dependent mean-field theories such as time-dependent local density approximation (TDLDA) in order to include incoherent dynamical correlations, which are known to play a key role in far-off equilibrium dynamics. We focus here on the case of irradiation dynamics in clusters and molecules. The field, still largely unexplored, requires quantum approaches which represents a major formal and computational effort. We present several approaches we have investigated to address such an issue. We start with time-dependent current-density functional theory (TDCDFT), known to provide damping in the linear regime and explore its capability far-off equilibrium. We observe difficulties with the scaling of relaxation times with deposited energy. We next briefly discuss semi-classical approaches which deliver kinetic equations applicable at sufficiently large excitation energies. We then consider a first quantum kinetic equation at the level of a simplified, though rather elaborate in its content, relaxation time approximation (RTA). Thanks to its sophistication, the method allows us to address numerous realistic irradiation scenarios beyond the usual domain of reliability of such theories. We demonstrate in particular the key role played by dense spectral regions in the impact of dissipation in the response of the irradiated system. RTA nevertheless remains a phenomenological approach which calls for more fundamental descriptions. This is achieved by a stochastic extension of mean field theory, coined stochastic time dependent Hartree–Fock (STDHF), which provides an ensemble description of far-off equilibrium dynamics. The method is equivalent to a quantum kinetic equation complemented by a stochastic collision term. STDHF clearly leads to proper thermalization behaviors in 1D test systems considered here. It remains limited by its ensemble nature which requires possibly huge ensembles to properly sample small transition rates. An alternative approach, coined average STDHF (ASTDHF), consists in overlooking mean field fluctuations of STDHF. ASTDHF provides a robust tool, properly matching STDHF when possible and allowing extension to realistic dynamical scenarios in full 3D. It can also be used in open systems to explore, as done in RTA, the competition between ionization and dissipation.


Introduction
The theoretical description of far-off equilibrium dynamics in quantum many-body systems remains since numerous decades a widely open and unresolved problem. First attempts can be traced back to Bohr's early pioneering work on charged-particle penetration and stopping in matter, not mentioning the seminal work on nuclear collision dynamics [1,2]. The nuclear physics case, namely dynamics of heavy-ion collisions with rich outcomes at a high level of experimental detail, has been extensively explored during the past four decades, but mostly in terms of classical or semi-classical modelings with only few excursions into approximate quantum treatments [3][4][5][6][7][8][9][10].
More recently, far-off equilibrium dynamics has become a key experimental and theoretical issue in transport processes in solids [11], in ultracold gases ("quenches") [12][13][14][15] and in laser irradiation of atoms, molecules and condensed matter [16][17][18][19][20][21][22][23][24][25][26][27][28]. These elaborate experiments pose a considerable challenge to theory. Quantum manybody systems far-off equilibrium indeed raise substantial difficulties due to the huge accessible phase space. Manybody correlations become increasingly important while a fully coherent description quickly gets out of reach. The way is to go for a reduced description (usually at one-body level) which introduces statistical concepts of dissipation and fluctuations [29]. And even this is a demanding task, the more so as the basic processes of excitation, dissipation through collisional redistribution, relaxation and coupling between very different species (electrons and ions) occur on vastly different time scales. This contribution discusses a couple of recent attempts to describe dissipative dynamics of fermion systems. Among the many fields of applications, we focus here on the generic case of thermally isolated, finite many-electron systems (atoms, molecules, clusters) excited by high intensity lasers or swift ions, although the developed theories should be applicable (with proper adaptation) to other scientific fields.
In such laser irradiation scenarios, electrons immediately react to the external electro-magnetic field driving the system far-off equilibrium (thus well beyond the Born-Oppenheimer approximation, in molecular terms) but still fully in the quantum regime, as clearly pointed out by recent experiments, which find clear quantum and thermal pattern in photo-electron spectra [30][31][32]. The excited electron cloud then rapidly thermalizes and will, on a slower pace, couple to the ionic degrees of freedom. Describing quantum mechanical far-off equilibrium dynamics of electrons is thus the first key issue to address while coupling to the ions comes second. We shall thus mostly focus here on electronic degrees of freedom, except for Section 5.4 where we briefly comment on the impact of thermal ionic motion.

Mean-field dynamics and its limitations
Current theoretical approaches in complex electronic systems are often based on time-dependent density-functional theories (TDDFT), in practice typically at the level of the time-dependent (adiabatic) local density approximation (TDLDA) [25,33,34]. TDLDA and similar approximate mean-field theories usually provide a correct picture of the dynamical evolution up to moderate excitations. But, by construction, they cannot describe correlation effects at higher energies, the related dissipation and particularly not the fluctuations. This is illustrated in Figure 1 for photo-electron spectra (PES) and/or photo-angular distributions (PAD) [35,36] emerging after laser irradiation of C 60 . The figure compares TDLDA calculations and experimental measurements on PES (panel (c)) and PAD (panels (b), (d) and (e)) following irradiation by various lasers [37], both in the multiphoton regime (panels (a) and (b)) or in the monophoton regime (panels (c)-(e)). The comparison shows both a qualitative agreement and discrepancies in quantitative detail, especially what concerns PAD, whose more or less isotropic character is hardly recovered within pure mean field dynamics. This indicates the importance of dynamical correlations beyond mean field which thermalize excitation energy and so deliver more isotropic emission in the multiphoton regime (panel (b)). Finally, panel (a) is a pure experimental result displaying PES slopes as a function of laser fluence in the multiphoton regime. The exponential slope of PES is here fully attributed to thermal effects (although this point can be theoretically debated [38]), delivering temperatures as high as typically 20% of the ionization potential of the system, which definitely corresponds to a very hot electron system. All in all, the figure points out the need of more elaborate approaches than the mere effective mean-field approach used here, with a proper account of dynamical correlations beyond time-dependent mean field. Even though we focus on extensions beyond meanfield, TDLDA remains the very basis of all further-going approaches. It will thus be reviewed briefly in Section 2.
The above calculations were performed at the level of the so-called adiabatic local density approximation (ALDA) of the exchange-correlation (xc) potential (Sect. 2), where the xc potential is an instantaneous and local functional of the density (r, t). This approximation may become questionable in violently dynamical processes. One way to go beyond ALDA, but still remaining in the framework of a local density-functional, is timedependent current-density-functional theory (TDCDFT), with the current density j(r, t) as a further crucial ingredient besides the local density (r, t) [39,40]. A practical local functional of the current density, the so-called Vignale and Kohn (VK) functional, can be derived from dynamical response in a weakly perturbed electron gas [40] and has been used since then with some successes for finite and infinite electronic systems within the linearresponse regime [41][42][43]. TDCDFT has also been used beyond linear response in a real-time description for models or simple systems [44][45][46], where symmetries allow the complexity of the equations to be reduced to that of TDDFT. These studies indicate that the VK functional introduces damping mechanisms into electron dynamics, which was also observed in the linear response of atoms [47]. TDCDFT within this VK approach could thus provide a way to simulate dissipative dynamics at affordable expense. This will be tested in Section 3.

Dynamical correlations and dissipative dynamics
The ground-state of a many-electron system features correlations, defined as everything not embraced in a basic mean-field description, such as Hartree-Fock (HF). Numerous well established methods treat ground-state correlations at various levels. Density-functional theories incorporate them implicitly while maintaining formally a mean-field description [48]. A great variety of approaches exists for an explicit description of correlations, see [49]. A coherent superposition of mean-field states is often used in case of a few dominating correlation channels [50,51] or, otherwise for an overseeable amount of correlating states as in multi-configurational HF [52]. More involved situations can be attacked with advanced wave function-based methods such as complete-active space self-consistent field, multi-reference configuration interaction or coupledcluster methods [52,53], as well as with alternative routes through reduced density matrix functional theory, Green's functions and many-body perturbation theory [54][55][56][57].
Even more demanding and still less well developed is the description of correlations in the time evolution of interacting quantum many-body systems. In analogy to the static case, these correlations correspond to all deviations from time-dependent mean-field theories such as TDHF [57] and, again, a great background of correlations can be accounted for effectively by TDDFT approaches [25,33,34,[58][59][60][61][62][63][64]. These are very efficient and even allow Typical ionization observables, that is photo-electron spectra (PES) and photo-angular distributions (PAD), from C60 irradiated by fs laser pulses, in the multiphoton (panels (a) and (b)) or in the monophoton regime (panels (c) to (e)). In (a), the electronic temperature is extracted from an exponential fit of experimental PES [30]. The other panels are adapted from [37]. applications to large systems, but with only approximate accounts of exchange and correlations. There thus remain genuinely dynamical effects unaccounted for, as, e.g., particle-particle correlations in multi-electron emission or energy relaxation towards a thermal limit. These correlations are referred to as dynamical correlations (DC).
The description of DC in systems far-off equilibrium remains, to a large extent, open for development. The time-dependent extensions of ground state approaches as, e.g., the multi-configurational TDHF (MCTDHF) [65][66][67], can only access DC for short times, small systems and mostly low excitation energies, with little hope to extend them to the dissipation dominated regime [18,20,24]. The recently introduced time-dependent second-order reduced density matrix (TD-2RDM) theory [68] delivers a quantum kinetic equation with results almost equivalent to MCTDHF ones but with better computational performances. TD-2RDM can thus address larger systems (<10-20 electrons) on longer times (<10-20 fs) even if realistic applications remain very limited. Although potentially "exact", these schemes grow with increasing excitation (thus phase space) quickly untractable unless one escapes, again, to approximate schemes as reconstruction of many-body functionals, purification, etc. There also exist alternatives beyond standard TDDFT [69] namely works based on dynamical extensions of density matrix functional theory [34,70], or non-equilibrium Green's function approaches [71], even with a few attempts beyond the linear response [72,73]. But these highly elaborate approaches are again strongly limited in system size and time span, not mentioning their complexity.
Semi-classical (and even classical) methods provide a manageable, approximate treatment for large systems and high excitation energies where they can work very well [27,[74][75][76]. But they leave a gap between low and high energies thus being inapplicable to many irradiation scenarios with dominant quantum effects in the entrance channel [30][31][32]. To partially fill this gap, we have developed a quantum method, well suited to high excitation energies, by approximating the correlated wave function by an incoherent (stochastically explored) set of TDHF states (coined STDHF [77][78][79] following its nuclear physics origin). STDHF effectively treats DC via a quantum collision term evaluated in stochastic manner. In its general form, it allows one to cover dynamics with huge fluctuations of the mean field as it appears in electron transfer or multi-fragmentation processes. At this point, we have to remark that stochastic mean-field approaches on top of TDDFT are often used in a different context, namely for open quantum systems where stochastic elements are exploited to model the interaction with an external thermal bath. There is a considerable amount of literature on that, see e.g. [80][81][82][83][84][85]. In contrast to that, the approaches here deal with thermalization within an isolated (closed) system. STDHF can be reduced to an average version thereof (then equivalent to a quantum kinetic equation) by imposing a common mean field to all members of the stochastic ensemble (coined average STDHF, ASTDHF [86]). AST-DHF offers an efficient method for large systems on long times and for lower excitation energies (where fluctuations of the mean field are negligible). A further step of simplification is attained by approximating the quantum collision term by a relaxation time approximation (RTA). We have thus developed a phenomenological RTA approach [87,88] in parallel to ASTDHF, allowing realistic simulations in full 3D for long time dynamics. Both RTA and (A)STDHF are expected to deliver a quantum description of DC in the range where the latter are dominated by incoherent effects building up in the course of the dynamical evolution. These theories are therefore ideal candidates to provide a quantum description of dissipative effects in such dynamical scenarios.
The range of applicability of major available theories is schematically illustrated in Figure 2. Guiding line is the one-body entropy as a function of time from attoseconds to picoseconds for a typical cluster/molecule irradiation by a short (below a few tens of fs) laser pulse of moderate intensity. Entropy provides here a typical indicator of DC driving thermalization. The upper (red thick) curve represents DC attainable by coherent approaches such as MCTDHF or TD-2RDM (green curve). The latters include all correlation effects but are limited to (very) short times (attoseconds to femtoseconds). While probably much less accurate for short times, incoherent approximations based on incoherent electron-electron scattering (lower blue thin curve), such as STDHF, AST-DHF and RTA (and semi classical approximations for sufficiently energetic processes and/or long times) become increasingly accurate at later times (for an example in an exactly solvable model, see [89]). They provide a description of thermalization ("heating") until the system has reached a thermal stage, whose characteristics depend on total deposited excitation energy/irradiation time as schematized by the broken curves between short ( 10 fs) and intermediate (∼100 fs) times. This is precisely the aim of this paper to analyze the performances of such incoherent approximations. This is done in the following steps. After a brief summary of TDDFT in Section 2, we address the TDCDFT extension thereof (Sect. 3). Section 4 will then briefly recall semi-classical approaches in the field and provide and introduction to RTA which we shall introduce phenomenologically in Section 5. The following Sections 6 and 7 will be devoted to STDHF and ASTDHF respectively. Each section presents a brief review of the theoretical background together with typical examples of applications.

TDDFT
2.1 Basis of our description: time-dependent density functional theory (TDDFT) The state of the electronic system is described in terms of single-particle (s.p.) wave functions and their occupation probabilities. In the standard Kohn-Sham (KS) picture of DFT [90], occupied one electron states ϕ n (n = 1, . . . , N ) follow effective one-body Schrödinger equations : where U S = U ion + U ext + U H + U xc is the KS potential successsively composed of the ionic, the external (laser), the Hartree, and the xc potentials. The ionic background is described by soft local [91] or Goedecker-type [92] pseudopotentials.
We use the TDLDA [48,93,94] in which the exchange correlation potential U xc is approximated by a function of the density (r, t), (U xc → U ALDA xc ). TDLDA is augmented by an average density self-interaction correction (ADSIC) [95] to obtain a correct simulation of ionization properties [96]. The latter approximation nevertheless only involves the local density (r, t), constructed as: so that the exchange correlation potential keeps a simple form U xc (r, t) = U xc [ (r, t)], as a function of the density (not a functional). The associated one-body density operator ρ is defined by The equation of evolution of ρ reads as: For later use (in particular in Sect. 5), we summarize the mean-field propagation according to TDLDA formally as: where U (t, t ) is the unitary one-body time-evolution operator, T the time-ordering operator, h the KS meanfield operator, and U S the (density-dependent) actual KS potential [48].

Numerical details
The static and time-dependent LDA equations are solved in real time on a grid with standard techniques [59,97]. While TDLDA calculations can be routinely performed in full 3D [25,33,34], approaches beyond mean field rapidly become very demanding computationally speaking. This is the case of TDCDFT whose full 3D implementation beyond linear response domain is feasible but difficult (Sect. 3). Dissipative approaches built on top of TDLDA, as discussed here, are still in a developing stage and only few examples of applications have been done in full 3D, for isolated system, see [88] and for extensions of TDCDFT in open systems, see [82,84]. The computations presented in Section 5 were thus performed in the cylindrically averaged pseudo-potential scheme [98,99], which has proven to be an efficient and reliable approximation for metal clusters close to axial symmetry. Wave functions and fields are then represented on a 2D cylindrical grid in coordinate space [100]. The most recent and elaborate developments on the basis of stochastic extensions of mean-field dynamics (Sects. 6 and 7) were performed in 1D, although extensions to full 3D are straightforward, especially in the simplified version of the theory (Sect. 7). Finally, semi-classical calculations (Vlasov and VUU in Sect. 4) are performed using dedicated techniques, on the basis of the propagation of numerical particles, thus rather different from grid approaches [4].
To solve the (time-dependent) KS equations for the s.p. wave functions, we use time-splitting for time propagation [101] and accelerated gradient iterations for the stationary solution [102]. The Coulomb field is computed with successive over-relaxation [100] or Fast Fourier Transform techniques [103]. To access ionization dynamics, we use absorbing boundary conditions [97,104] which gently absorb outgoing electron flow at the boundaries of the numerical grid. A proper choice thereof allows us to avoid artifacts from erroneously back-scattered electrons and thus provides a relevant description of ionization with limited/controlled numerical artifacts [105].

Excitation mechanisms, lasers and boosts
The external, coherent laser field is handled as a classical field in the long wavelength limit, adding to the mean-field Hamiltonian the potential with θ the Heaviside function. The laser characteristics are: linear polarization e z along the symmetry axis, peak field strength E 0 related to laser intensity (I ∝ E 2 0 ), photon frequency ω las , and total pulse length T pulse . The full width at half maximum (of intensity) is given as F W HM T pulse /3.
Typical pulse durations in experiments with intense laser fields are in the range of a few fs to a few dozens of fs. This time scale lies in the range between electronic and ionic times. Shorter pulses are needed to discriminate electronic time scales. Thus we also employ extremely short pulses, idealized by an initial and instantaneous boost of the wave functions in the spatial z direction as |ϕ n (0 + ) = e ik b z |ϕ n (0 − ) with boost wavevector k b and corresponding excitation energy E * = 2 k 2 b /(2m).

Typical observables
The analysis of electron dynamics is performed through a set of well established observables [27,33,64,97] that we briefly recall here. We start with the dipole response which is the routine work of spectroscopic studies but which is also known to play a key role in numerous irradiation scenarios, even beyond the linear (spectroscopic) regime [27,33,59,106]. We compute the dipole moment of electrons (with respect to ionic background) as: from which, by time-frequency Fourier transform, D(t) → D(ω), one accesses photoabsorption cross-section and spectroscopic analysis, via the strength function Im(D(ω)) or the power spectrum | D(ω)| 2 . The latter is often called optical response and will be used at some places below.
In the far-off equilibrium regime, electron emission becomes an important mechanism. It can be analyzed at various levels of sophistication thanks to the use of absorbing boundary conditions [97,104], starting from the rough total ionization where N is the initial electron number and where integration runs over the computing box. Electron emission can be analyzed at a finer level of detail by computing energy-resolved (photo-electron spectra, PES) or angularresolved (photo-angular distribution) electron spectra. This also allows direct comparisons to experimental results when available [64]. The strategy consists in both cases to define a set of "measuring points" r M near the absorbing boundaries, and to measure the time evolution of s.p. wave functions ϕ n (r M , t) at these points. The time integration of |ϕ n (r M , t)| 2 provides ionization at each measuring point, thus immediately delivering the PAD, once properly accounting for the solid angle Ω r M associated with the direction of r M [38]. The PES in turn, is obtained by Fourier transforming the s.p. wave functions (augmented by an additional phase factor taking into account the external laser field) [107] from time to energy ϕ i (r M , t) −→ ϕ j (r M , E kin ) and defining kinetic energy as E kin = 2 k 2 /(2m) = ω. This delivers the PES at r M (within solid angle Ω r M ) by summing up contributions from each s.p. state The total yield is then simply obtained by integration over solid angle delivering Y(E kin ) = dΩ r M Y Ωr M (E kin ).

Theory
Standard TDDFT employs a scalar multiplicative KS potential U S (r, t) composed of the external field U ext , the Hartree potential U H and an xc potential U xc (Sect. 2.1). The most widely used practical approach is adiabatic LDA (ALDA) which expresses these potentials in terms of the local, instantaneous density (r, t). The next step and natural extension is to look for the impact of the local current j(r, t) with the hope to get hold of more aspects of manybody effects in dynamical response. The current j is a vector field and, accordingly, TDCDFT complements the KS equations by a vector potential A S to where q = −e is the electron charge and c the light velocity. The vector potential A S is composed of an external contribution A ext and the xc vector potential A xc . The density (r, t) given in equation (2) and the current density computed as: can be expressed in terms of the KS wave functions and are independent of the gauge chosen to represent the electromagnetic potentials. Widely used is the Vignale-Kohn (VK) approximation which derives the xc vector potentials from linear response to a dynamical perturbation [40]. Up to second order in spatial derivatives, under the basic assumption that the gradients of the density and the velocity are small (see [108]), and by choosing a gauge with U xc = 0, the VK functional in real time reads where the second term on the right-hand side is usually referred to as the "memory term". It is determined by the visco-elastic xc stress tensor σ xc,ij (r, t) which is a function of the time-dependent velocity field v(r, t) = j(r, t)/ (r, t). It can be expressed through coefficients of shear η(r, t, t ) and bulk viscosity ζ(r, t, t ) where the dependence on two times, t and t , indicates that the viscosities in general can embody memory effects. The η and ζ, in turn, depend directly on the longitudinal (L) and transverse (T) response kernels f L,T xc ( , ω) of the homogeneous electron gas. Memory effects are carried in its frequency dependence. Practicable approximations for new xc kernels are given in [109][110][111]. For the following example, we use the VK parametrization from [110,111], which was designed to satisfy all known limits and relations of the longitudinal and transverse xc kernels, and which allows a simple analytic evaluation of η(r, t, t ). Still, with its third order derivatives and double integrals over time, the VK approximation renders the KS equations highly involved. To simplify the original equations and to make computations feasible, while keeping the important aspects of the VK approximation, we make an instantaneous approximation for the basic ingredients, the viscosities η(r, t, t ) and ζ(r, t, t ). This is legitimated by the fact that their memory kernels contribute only during a time interval t ∈ [t − T, t] which is much shorter than the most relevant time scale, namely the plasma period of the electron cloud. This allows us to take the instantaneous η 0 (r, t) and ζ 0 (r, t) which are the η and ζ integrated over t . This yields a negligible longitudinal viscosity ζ 0 = 0, leading to the following expression of the instantaneous stress tensor (in 3D) [110,111]: For further details of the evaluation and the value of η 0 , see [112]. It is important to note that the above instantaneous approximation yields a purely real viscosity η 0 and it is precisely this real part of the complex viscosity kernels that is responsible for the dissipative effects described in TDCDFT [45]. This is a general feature also observed in theories dealing with explicit collision terms: a Markovian (instantaneous) approximation exclusively describes the dissipation [113,114].

Numerical difficulties
Even in its simplified form, the implementation of the VK functional in a real-time propagation of finite systems remains demanding because it involves third order derivatives of the wave functions and because the velocity v = j/ becomes an unsafe quantity for very small , e.g., in the tails of the electron cloud. Amplified by high-order derivatives, this leads in a few time steps to insurmountable numerical instabilities. We have recovered stability by using two essential measures: (i) for the derivatives in the current functionals, we use merely second-order finite differences which, in presence of fluctuations, are more robust than higher orders; (ii) in the computation of the velocity, we employ a numerical smoothing technique to eliminate the spurious signal in the low density area, which consists in driving progressively the velocity toward a limiting constant while ensuring its continuity along all axes.
In addition, TDLDA calculations are performed on valence electrons, the core electrons and the remaining atomic nucleus being modeled through a pseudo-potential. There are thus two kinds of density: the total one, , stemming from all electrons, and the (partial) density from valence electrons only. The latter is used in the (conventional) scalar KS potential, while we have to use the total density in the above expressions of the TDCDFT functional.

Some results
We show here results for two closed-shell systems, each with 2 active electrons, namely the Mg atom and the Na 2 dimer, both computed in full 3D. We excite the systems by an instantaneous initial boost k b (Sect. 2.2) of the valence electron cloud and then let the system evolve freely. Note that we do not use absorbing boundary conditions here. The system we explore is thus closed, so that we can focus on how the initial excitation energy redistributes toward "thermal" degrees of freedom, thus omitting the competing de-excitation channel through direct electron emission. Note finally that both systems possess one dominant dipole mode with only very little spectral fragmentation elsewhere. This avoids interference with damping of the dipole signal by distribution over many sub-modes (Landau damping [59,115]) and thus makes it ideal test cases for checking dissipative contributions from TDCDFT. Figure 3 compares the time evolution of the dipole moment, see equation (8), calculated using TDLDA (black dashed lines) and TDCDFT within the VK approximation (solid curves), after an initial boost of k b = 0.01 a −1 0 (dark blue) or 0.1 a −1 0 (light green). The dipole moments are rescaled proportionally to the inverse of the initial boost, so that they would become identical for an undamped motion.
We first discuss the difference between TDLDA and TDCDFT. The dipole moment calculated in mere TDLDA carries on to oscillate without visible damping. Instead, the TDCDFT functional produces a strong damping. This was already observed in previous works on model systems [44,46], although it had been shown to be overestimated in small systems while becoming correct in the thermodynamic limit [46]. The VK approximation, hence, can describe dissipation even in a real-time description and in the instantaneous approximation (14).
We also see that the position of the oscillation peaks remains the same as in TDLDA. This can be understood from the fact that the sequence of approximations used above to derive the visco-elastic coefficients η 0 and ζ 0 corresponds to a purely imaginary transverse kernel f T xc which can affect only the width of a mode. It requires the real part of viscosity (thus memory effects) to produce an effect on the position of the oscillation peaks [47].
We have also explored the time evolution of the total energy (not shown here). As it should be, the latter remains constant with TDLDA, while dissipation in the course of time is observed when using VK. It has been indeed demonstrated that the adiabatic KS energy decreases monotonically in the absence of an external field when the VK approach is used [116], which indicates that the system is irreversibly driven to equilibrium. This scenario has already been observed in previous TDCDFT works on 1D models [44,46], where it is argued that the energy loss is, in fact, distributed over configuration space (thermalization). However, the compensating energy increase by thermalization is not accounted for in the present form of TDCDFT.
We now compare the time evolution obtained for two different initial boosts. From Figure 3, we actually observe only little dependence on k b . This can be understood since the viscosities, combined with the instantaneous approximation (14), do not exhibit any entry for the actual excitation energy visible. However, this lack of energy-dependence is unphysical: damping does depend on excitation energy, as seen in detailed microscopic descriptions of electron-electron collisions, as in Fermi liquids in bulk [117] and in finite electron systems [74,75,87], and as in the dissipative theories discussed later on in this review. The VK approximation to TDCDFT does not reproduce this expected trend. It is therefore limited to the linear regime and should better not be used for non-linear excitations. This is consistent with the way in which the VK functional has been derived, i.e. for a weakly perturbed electron gas.

Semi-classical approximation and Vlasov-LDA
Semi-classical approximations to fully fledged quantum theories are expected to provide valid and valuable schemes at sufficiently high excitation and for sufficiently large systems. They have been explored for describing laser irradiation of simple metal clusters [27,74,75], but extension to other materials remain still an open question. A simple way to introduce such semi-classical descriptions is to start from TDLDA and to use the Wigner transform which provides the natural starting point for expansions in orders of [118]. (A more stable expansion is obtained, in fact, by using the Husimi transform [119] which, however, yields the same result in lowest order.) This delivers at lowest order in the Vlasov-LDA equation as the semi-classical analog of equation (4), in terms of the phase space distribution f (r, p, t), itself the semi-classical analog of ρ: The semi-classical Hamiltonian has the same density dependence as its quantum mechanical counterpart, using now the semi-classical estimate of the density of matter from the phase space distribution as (r, t) = d 3 pf (r, p, t).

From Vlasov-LDA to kinetic equations
The semi-classical approximation provides a simple framework to complement mean field dynamics (Vlasov) by dynamical correlations through a collision term "à la Boltzmann". The latter is to be augmented by a phasespace factor which properly accounts for Pauli blocking effects delivering together the Vlasov-Uehling-Uhlenbeck (VUU) equation [120] which reads: The collision term Photo-angular distribution of Na + 41 following a laser irradiation corresponding to an "on"-resonant (ω las = 2.7 eV, bottom panel) and an "off"-resonant (ω las = 2.3 eV, top panel) case. We compare mean fields approaches (quantum TDLDA and semi-classical Vlasov) to the semi-classical VUU and the dissipative quantum RTA.
where we use as compact notation f i = f (r, p i , t) and where dσ/dΩ is the elementary in-medium scattering cross section. The latter is associated to the residual interaction not contained in the mean field. The sharing between mean field and residual interaction is properly accounted for (i.e. no double counting) by using a screened Coulomb interaction [121][122][123]. An important property of the collision term is that each collision is computed at a given point in space and that it locally enforces conservation of momentum and s.p. energy (which reduces to kinetic energy for a local potential). The performance of VUU, particularly its dissipative aspect, is illustrated in Figure 4, where we plot the PAD of electrons emitted from a medium size metal cluster Na + 41 . following excitation by a laser of moderate intensity. We consider two laser frequencies, one in the dense part of the dipole response spectrum (ω las = 2.7 eV, bottom panel), corresponding to an "on"-resonant case, and one outside (ω las = 2.3 eV, top panel), that is safely "off"-resonant. The figure provides a comparison between various levels of theories. The starting point is TDLDA (Sect. 2) which, as expected, delivers a PAD with alignment of emission (thin red curves) along laser polarization axis, especially for the resonant case in the bottom panel. We first discuss this resonant case as it is the more telling one. The Vlasov-LDA, as a semi-classical approximation to TDLDA, delivers a very similar angular distribution (thin dashes), with preferential forward/backward emission along laser polarization axis. The VUU case (thick dashes) is notably different, producing a much more isotropic PAD. This is a direct consequence of the account of collisional correlations which have allowed some thermalization of the electron cloud. This leads to a balance between directed emission, well described at TDLDA/Vlasov-LDA level, and isotropic emission resulting from the impact of electron-electron collisions, and ultimately leading to fully isotropic emission. Mind that both processes compete in terms of time scales with direct emission at short times, followed by isotropic one (reflecting thermalization) at later times. The collision term I coll [f ] thus, as expected, leads to more isotropic electron emission. The situation is qualitatively similar but with less marked effects in the off-resonant case.

Semi-classical relaxation time approximation
As mentioned above, electron-electron collisions in the VUU collision term I coll [f ] are local, modifying for a given r only the momentum distribution at this specific point. This has as consequences the observation of local conservation laws [124] which preserve local density (r, t), local current j(r, t), and local kinetic energy E kin (r, t) (in the case of a local potential). A collision event thus drives the system towards a local and instantaneous equilibrium which can be characterized by local density, local current and local kinetic energy f eq (r, p; , j, E kin ). As usual in kinetic theory, global equilibrium is established later on by interplay with the long range transport resulting from mean-field (Vlasov) propagation. The VUU equation with its involved nine-fold integration in the collision term (17) can then be simplified in terms of a RTA which reads [125][126][127][128] where f eq is the thermal equilibrium at given (r, t), j(r, t), E kin (r, t). Note that these constraints are fully local in space and time. The crucial parameter τ relax will be discussed in more details in the next section around equation (20b). The RTA approximation is interesting because it has a simple form which can be rather easily generalized to a quantum formulation. This is precisely the topic of the forthcoming Section 5. Before addressing that case, let us mention that this quantum RTA approximation has also been shown in Figure 4 (thick full curves) and compared to VUU. Interestingly enough, quantum results (at RTA level) are very close to VUU ones. And again, one observes that, generally speaking, dissipative effects are less influential in the off-resonant case (top panel). This is not an accident but rather a general rule. Indeed, as we shall see in Section 5, resonance conditions play a key role on the impact of dissipation. While a semiclassical approximation is able to recover gross features of the excitation spectrum, like an average plasmon peak, it is unable by construction to address finer details so that its eigenfrequency spectra are in general different from the quantal one. As a consequence, the impact of dissipation will be different, which, once more, justifies a quantum account thereof.
5 Quantum relaxation-time approximation 5.1 Brief review on RTA Mere TDLDA could be described by pure states (Slater states) having s.p. occupations only 0 or 1. Dissipation produces mixed states described formally in a compact way by a one-body density operator, which reads, in natural orbitals representation: The sum runs up to Ω which is the size of the configuration space and which has to be larger than the actual electron number N . The weight ν n represents the occupation probability of s.p. state |ϕ n . Pure mean-field propagation (5) leaves the occupation weights ν n unchanged and propagates only the s.p. states which reads for the density Instead, dynamical correlations through instantaneous two-electron collisions change the occupation weights. RTA approximates this as relaxation towards a local, instantaneous equilibrium state ρ eq [ , j, E] [129] as: where h[ ] is the KS Hamiltonian (5c) in TDLDA (with ADSIC), depending on the actual local density distribution (r, t) = n ν n |ϕ n (r, t)| 2 . The relaxation time τ relax is estimated in semi-classical Fermi liquid theory. For metal clusters serving in the following as test cases, it becomes where E * intr is the intrinsic (thermal) energy of the system, N the actual number of electrons, σ ee the in-medium electron-electron cross-section and r s = (3/(4π )) 2/3 is the Wigner-Seitz radius of the electron cloud [129]. It employs an average density because τ relax is a global parameter. This approximation is appropriate for metallic systems with their rather homogeneous electron density. σ ee also depends on this average density. The actual σ ee is taken from the careful evaluation of [122,123] computing electron screening for homogeneous electron matter in Thomas-Fermi approximation. This yields σ ee = 6.5 a 2 0 for the case of small Na clusters whose Wigner-Seitz radius at zero temperature becomes typically r s ≈ 3.7 a 0 , somewhat smaller than the bulk value due to surface tension.
The instantaneous equilibrium ρ eq [ , j, E] entering equation (20a) is the thermal mean-field state of minimum energy under the constraints of given local density (r, t), local current j(r,t), and total energy E(t). It is computed by density constrained mean-field (DCMF) techniques as developed in [130], extended to account also for the constraint on current j(r). The weights ν (eq) n are determined according to thermal equilibrium where the temperature T is tuned to reproduce the desired total energy E, for details see [129].
The mean-field propagation (5b) of the s.p. wave functions is done on a fine time mesh with δt = 0.005 fs. The right-hand side in equation (20a) is evaluated at coarse time intervals ∆t, typically 0.2-0.5 fs. To that end, first, the actual , j, and E are computed. These are used to determine the local-instantaneous equilibrium state ρ eq . This is used to step to the new one-body density ρ(t+∆t) = ρ + (∆t/τ relax ) ρ − ρ eq [ , j, E] . In a final cleanup, this new state ρ(t+∆t) is mapped into natural orbitals representation, see equation (19), thus delivering the new s.p. wave functions ϕ α (t + ∆t) and occupation weights ν n (t+∆t), from which on the next step is performed. More details can be found in [129].

Time scales
A prominent observable characterizing the electron dynamics is the dipole moment equation (8) We consider here the component D z (t) ≡ D(t) in direction of laser polarization. The lower panel of Figure 5 shows the time evolution of the dipole moment for Na + 9 after an initial instantaneous boost depositing and excitation energy E * (Sect. 2.3). The dipole signal is heavily oscillating. Nonetheless, one can spot a global trend to shrinking dipole amplitude whereby RTA, not surprisingly, produces a faster attenuation than TDLDA.
For a further analysis, we define a smoother observable from the dipole signal to overcome the oscillations as the typical "dipole energy" which one can define as: The definition is applicable in connection with metal clusters where the dipole signal is dominated by the frequency of the Mie surface plasmon ω plas [131]. The prescription might fail for a system with well separate different eigenmodes in which case we should replace ω plas D → −D (or the envelope of the dipole signal). The upper panel of Figure 5 shows the time evolution of this dipole energy (21) (in logarithmic scale). It is obviously orders of magnitude smoother than the dipole signal as such and can serve as reasonable basis to deduce relaxation times. One can distinguish two clearly different phases of relaxation. There is a fast drop of the signal within the first 10-20 fs. This is caused by the spectral fragmentation of the dipole strength because the Mie plasmon couples to the 1-particle-1-hole states in its spectral vicinity. This drives the dipole signal quickly out of phase, a process known as Landau damping in bulk plasma [132]. After going to near extinction, the signal revives by sort of fluctuating recurrences and as overlaid with some long range damping. Relaxation in this late phase is caused by slowly  (21), for Na + 9 after an initial instantaneous boost, providing an excitation energy of E * = 2.7 eV, corresponding to the energy of one plasmon with ω plas = 2.7 eV. Compared are results with dissipation (RTA) and without (TDLDA). The dashed lines indicate fits to exponential relaxation, the black dashed line stands for a fit in the early time window t ∈ [0, 20] fs and the colored lines for a fit in the late window t ∈ [100, 400] fs. continuing electron emission and by conversion to intrinsic energy, the latter more pronounced in RTA.
We explore these relaxation times by fitting an exponential to the dipole energies (see dashed lines in Fig. 5). As already observed above, we distinguish two different phases of relaxation: An early phase up to the first deep minimum in the signal, i.e. in the range from t = 0 up to ∼10−20 fs, depending on boost strength k b , and a late phase from 50−100 fs on (again depending on k b ). Each time window provides a different exponential fit and therefore, a different relaxation time. They are all collected in Figure 6 which plots their dependence on the boost energy E * . The lower panel shows the damping times deduced in the early time window. They are practically the same for TDLDA and RTA because the early stages are dominated by Landau damping. Dissipative processes can take place only after Landau damping has weakened the Fermi surface and so opened the space for two-body collisions [54,97,117]. Although expected, the closeness of the two values is astonishing. Besides that, we see a significant dependence on boost energy, particularly for small E * . This is most probably due to the increasing contribution from direct electron emission which is a very fast process (within 1-3 fs). Instead, for the late phase (upper panel), we see large differences between LDA and RTA. This shows immediately the dissipation, modeled in RTA, at work. However, even then, the late relaxation times remain considerably larger than the time for initial There is also an effect on the spectral distribution of dipole strength [106]. A comparison between RTA and TDLDA is shown again for the case of Na + 9 in Figure 7. The boost energy was chosen such that it covers exactly one quantized plasmon state which provides a relevant picture of the intrinsic excitation related to a plasmon [133] and so calibrates the amount of dissipation in RTA properly. The RTA spectrum is, indeed, somewhat smoothed as compared to that from LDA. However even if the effect is visible, it is rather small, which complies very well with the late relaxation times in Figure 6. Indeed, at E * = 2.7 eV, one reads a relaxation time of about 200 fs, which relates to an extra amount of width of order of 0.01 eV, in accordance of the smoothing we see in Figure 7.

Energy balance
We now consider the impact of laser pulses with finite length T pulse of order of tens of fs. The system absorbs a Fig. 7. Spectral distribution of dipole strength calculated with TDLDA and RTA for Na + 9 obtained with spectral analysis after an initial instantaneous boost chosen such as E * = ω plas to cover one plasmon quantum.
certain amount of energy E abs from the laser field where E (mask) abs is a correction for the small amount of potential energy lost by electron emission at the absorbing bounds. The question is how this absorbed energy is distributed over the system. A first contribution comes from the emitted electrons, with some energy taken away, while another part is eaten up by an increasing potential energy from charging. Both together compose the charging energy E charging . The remaining energy delivered by the laser is shared between a collective kinetic energy E kin,coll and an "intrinsic" excitation energy E intr of the electron cloud itself consisting of a kinetic E intr,kin and a potential E intr,pot component. All terms sum up to absorbed energy [88]: The intrinsic kinetic energy is obtained from densityconstrained mean-field (DCMF) as (24) where E TDLDA (t) is the actual TDLDA (+ADSIC) energy and E DCMF ( , j, T = 0) is the DCMF energy at T = 0 (=ground state for fixed and j). A simpler estimate can be obtained from a Thomas-Fermi approximation to E DCMF ( , j, T = 0), for details see [129]. We use this estimate for counter checking. The intrinsic kinetic energy (24) accounts for thermal intrinsic energy. A further contribution also being an intrinsic kinetic energy stems from the collective flow, namely E kin,coll (t) = d 3 r j 2 (r, t) 2m (r, t) .
Since E kin,coll = E DCMF ( , j, T = 0) − E DCMF ( , j = 0, T = 0), this shows that E kin,coll is part of the intrinsic energy, but a coherent one. It turns out that E kin,coll is generally a small contribution and it fades away on the long run due to dissipation. Therefore, we can safely ignore it for the following considerations. We do not consider in detail the charging energy. But we will look at charging as such (i.e. ionization) which is quantified in terms of the number of escaped electrons N esc (t) defined in equation (9). For all observables, it is important to keep in mind that TDLDA and RTA produce averages. For example, N esc represents the ionization averaged over an ensemble of similar processes; a detailed distribution of (integer) ionization stages may be recovered approximately from the final wave functions [134], situations which we do not consider here. Figure 8 exemplifies trends of intrinsic energy E intr,kin (lower panels) and ionization N esc (upper panels) as functions of laser frequency ω las comparing TDLDA with RTA. The right panels show scans for fixed laser intensity, the strategy usually followed when measuring photoelectron cross-sections. In the upper right panel, the total ionization N esc is thus proportional to the photo-electron cross-section with the typical fluctuating pattern of resonant and non-resonant situations with huge differences in total yield, corresponding to largely differing energy absorption. Particularly at the resonant points with high yield, we see considerable differences between TDLDA and RTA whereby the former generally delivers more electron emission. The dissipation in RTA obviously manages to convert a larger fraction of invested laser energy into internal heating.
The energy E intr,kin as such undergoes similar hefty fluctuations. To elucidate the internal energy branching, we consider the intrinsic energy relative to the total absorbed energy, i.e. the ratio E intr,kin /E abs , see lower right panel of Figure 8. Generally, there is a global trend to decreasing E intr,kin with increasing ω las . This is plausible because the lower frequency photons excite more intrinsic states below emission threshold. Besides the global trend, E intr,kin /E abs fluctuate with laser frequency, reflecting the fluctuations of absorptions. An important feature arises in the narrow frequency region of the plasmon resonance around 2.7 eV, corresponding to the region of strongest electron emission: there is a remarkable difference in relative intrinsic energy in that RTA channels off more energy into E intr,kin than TDLDA. This feature is nicely correlated with the difference in emission (upper right panel) leaving less N esc with RTA. This systematic difference is totally absent for higher frequencies above ionization potential (IP). The reason is that in this region direct onephoton emission prevails which is a fast process leaving little time for dissipative processes.
The above analysis scanning frequency for fixed laser intensity has the disadvantage that it runs over very different dynamical regimes due to the huge changes in dynamical susceptibility. Therefore, we complement the analysis by a frequency scan where the laser intensities were tuned for each frequency such that the absorbed energy was constant at E abs = 8.2 eV. The results are shown on the left panels of Figure 8. Ionization (upper left panel) now remains of the same order of magnitude throughout, confirming that we are, indeed, dealing with comparable dynamical situations. It is gratifying to see that this different way of looking at frequency dependence comes to the same result which can be summarized as: dissipation, modeled with RTA, gets its grip in the multiphoton regime for ω las below IP (indicated by the vertical dashes), becoming particularly active at resonance while dissipation plays a negligible role in the regime of direct emission triggered by one photon (above IP).
The laser frequency scan done in two manners indicates that laser intensity may also be an important parameter determining the amount of dissipation. Figure 9 shows the trend of intrinsic energy and ionization for a fixed (resonant) frequency ( ω las = ω plas = 2.7 eV) as functions of laser intensity I. There is, indeed, a moderate, but steady, trend. Larger intensity I leaves less intrinsic energy. This is the corollary of what we had seen when comparing right and left lower panels of Figure 8. At this given frequency ω las = 2.7 eV, we are still in the regime of two-photon emission. Increasing intensity makes twophoton emission faster which, in turn, leaves less time for conversion to intrinsic energy. The effect is, of course, more pronounced for RTA such that there is an increasing difference to TDLDA when the intensity decreases. The difference in N esc (upper panel) reflects that trend. There is one more interesting feature in the trend of relative ionization N esc /E abs : it decreases with increasing I. That looks at first glance surprising because we stay in the regime of two-photon emission throughout. The point is that the energy distribution of the emitted electrons changes with I. It becomes less steep with increasing I. This means that each emitted electron carries away relatively more kinetic energy which, in turn, requires more energy for emitting an electron thus reducing the output.  9. Energy balance as fraction of total intrinsic energy to total absorbed energy Eintr,tot/E abs as a function of laser intensity for laser pulses of total length T pulse = 96 fs and frequency ω las = ω plas = 2.7 eV.

On ionic motion effects
In several of the above examples, we have considered rather long laser pulses, and thus followed time evolution of the systems up to several hundreds of fs, which is enough to see sizable effects of ionic motion. Ionic motion has not been included in the above examples, though, in order to focus on electronic effects. The inclusion of ionic motion is strictly straightforward in RTA and exactly follows the path used to include ionic motion in TDLDA [97]. Therefore we have also started making a few test computations including ionic motion together with RTA. We are still in the process of complementing earlier findings with ionic motion. But first computations indicate that there is no major impact of ionic motion on electronic dissipation, at least in the few test scenarios considered. There remains to systematically check such preliminary results by varying the laser parameters. Work along that line is in progress.
Still, with respect to experimental observables, it is interesting to note that electronic temperature effects are not the single ones to be accounted for. Indeed, as soon as ionic motion is accounted for, one should also consider potential impact of ionic temperature on measurements. This effect is well known in cluster physics where temperature is a key ingredient in cluster production [135], and the impact of ionic temperature has been clearly observed in the optical response of metal clusters [136] and is suspected to be large in PES and PAD, for example in C 60 irradiation scenarios [37].
The point is illustrated in the case of Na + 9 studied at two different ionic temperatures in comparison to the zero temperature case. To generate the thermal ensemble, we first run a TDLDA-MD calculation, starting from the ground-state configuration at T = 0 K and giving each ion an initial velocity distributed stochastically around the intended mean velocity. After a few cycles of ionic oscillations, we obtain a randomly fluctuating cluster with the initial kinetic energy equally distributed between ionic potential and kinetic energies. We deduce the actual temperature of the ensemble from the kinetic energy averaged over time, as E kin = 0.5 × (3N ion − 6) T where 3N ion − 6 is the number of ionic degrees-of-freedom. We finally stochastically pick a certain number of configurations from this fluctuating system and store their ionic positions and velocities, constituting at the end of the day our thermal ensemble.
In all cases, we used a long laser pulse clearly covering ionic motion typical times scales which, in that case, lie in the 100 fs range. The laser intensity is tuned so that the total ionization remains small (a few percents in both cases). Figure 10 compares both PES (left panels) and PAD (right panels) at zero and non-vanishing temperatures, that is T = 150 K (top) and 475 K (bottom). The impact of temperature is rather evident in both cases with a smearing of PES peaks (increasing with ionic temperature) and a more isotropic PAD. Still, one might conclude that these effects do not appear extremely large. In fact, the amplitude of the effect strongly depends on the laser frequency, in a way similar to what has been observed in RTA, namely in relation to the spectral distribution [87]. The case displayed in Figure 10 corresponds to a laser frequency of ω las = 6 eV, slightly above the ionization potential and the domain of high spectral density, where effects are observed to be larger in amplitude [137]. Clearly, ionic temperature effects are thus to be accounted for in order to understand numerous experimental results, in complement to electronic effects which are the focus of the present paper.

Formalism
The idea underlying STDHF is to approximate Multi-Configurational TDHF [138][139][140] and to simplify it by occasional reduction of the correlated many-body state to an ensemble of mean-field states (Slater states). The practical description aims at an ensemble of N Slater states {|Φ (α) , α = 1, . . . , N }. Each single state formally possesses a pure-state N -body density matrix D (α) = |Φ (α) Φ (α) |. The N -body density of the full ensemble can then be written compactly as the following incoherent superposition We sketch now one step in the STDHF propagation. We start at initial time from the ground state of our theory, namely a Slater state, which will develop in time into an ensemble of Slater states via incoherent dynamical correlations. At a certain time, we are thus given an incoherent ensemble of D (α) . Each state α of the ensemble is propagated with the same scheme. We explain it here for one particular α. Starting from the given D (α) we propagate the full Schrödinger equation via an expansion in terms of a basis |Φ (α) κ of n-particle-n-hole (nph) states κ about the Slater state |Φ (α) . This yields the correlated N -body density with c (α) κ some weights unknown at this stage. Propagation is halted at a certain time t = τ STDHF , long enough to justify the assumption that the fast phase oscillations in the off-diagonal terms κ = κ wipe out their contributions in the average, as often done in non-equilibrium statistical physics [29]. This allows us to simplify the coherent expansion (27) to an incoherent sum of density operators as: where the probabilities w (α) κ are given by Mind that the physical quantities in the equations above are all time-dependent and are here evaluated at t = τ STDHF . We employ in equation (28b) time-dependent perturbation theory for computing the w (α) κ and reducing them in terms of Fermi's Golden rule [77]. The energies appearing in the energy matching δ-function above are the ones associated to the mean field states α and κ, respectively. The scheme implies that the mean field provides the dominant component to validate the use of time-dependent perturbation theory for a time τ STDHF which is long enough to develop the energy matching δ(E (α) (t) − E (κ) (t)), see [77,141,142] for the subtleties of time matching. The term W therein stands first for the residual interaction beyond mean field. However, closer inspection of the dynamical correlations shows that one should not insert here the given two-body interaction but rather an effective in-medium interaction which incorporates properly screening effects [113,114,121]. In practice, one chooses a simple form, e.g. contact interaction, whose strength reproduces the in-medium low-energy scattering cross-section.
The incoherent correlated state (28a) is a new, weighted ensemble of Slater states, so to say another ensemble within the ensemble (26). Carrying that through will soon lead to an insurmountable proliferation of states. To keep the scheme manageable, we reduce that back again to an unweighted ensemble (26) with a constant number of samples N . This is achieved by picking for each α in Monte-Carlo fashion one Slater state κ with probability w (α) κ from the ensemble equation (28a) and renaming that as the new reference state D (α) κ −→ D (α) . This then produces the next STDHF states D(t + τ STDHF ) in the form (26), ready for the next step in STDHF propagation.
Note that the time evolution of each state in the ensemble is dominated by the mean-field associated with |Φ (α) . This implies that STDHF allows to cover (large) fluctuations of the mean field within its ensemble. This becomes important in large amplitude dynamics with, e.g., electron transfer reactions where the final states may carry totally different charge states. On the other hand, handling so many mean fields in parallel is inefficient in processes with small fluctuations. We will come to that in Section 7.

Practical handling
The starting states |Φ (α) of the STDHF ensemble are built from single-particle (s.p.) states {ϕ (α) i , i = 1, . . . , Ω}, where i = 1, . . . , N stands for the occupied (or hole, (h)) states. These are complemented by a sufficient amount of unoccupied (or particle, (p)) states i = N +1, . . . , Ω. The p states supply space for the stochastic jumps. The enlarged s.p. space provides a frame for the nph excitations. The first correlation contribution beyond mean field comes from the 2ph states because 1ph excitations are already accounted for in the self-consistent meanfield state [143,144]. The states κ into which correlations develop are then all of the form The practical procedure then proceeds as follows:  t −→ t = t + τ STDHF according to TDLDA where U (α) is the mean-field propagator (5b) for Slater state α. 3. At the new time t , we evaluate the probability w (α) pp hh for a jump to κ = pp hh by using equations (28b) and (28c). Note that in this expression, the infinitely sharp δ function in the energy conservation is replaced by a finite width δ Γ function taken as δ Γ ( ) = θ( − Γ )θ(Γ − )/2Γ with θ is the Heaviside function. This accounts for a finite sampling time τ STDHF and provides, at the same time, the necessary smoothness to cope with discrete spectra of finite systems. The probability of no jump, κ = 0, thus keeping the unmodified state Φ (α) , is w pp hh . 4. Finally, we pick one particular κ in standard Monte-Carlo fashion with probability w (α) κ . The search for good candidates in the step 3 is accelerated by requiring for the difference of s.p. energies ε p +ε p −ε h −ε h ≤ 2Γ as pre-selector. The energy resolution of the scheme is given by the finite width Γ . The formally correct choice Γ = /τ STDHF in connection with practical τ STDHF could easily produce poor energy conservation. In practice, we use Γ as a free numerical parameter to achieve a good compromise between energy conservation and sampling quality. The hope is that reality is often more forgiving than strictly derived bounds.
The computation of observables follows standard prescriptions with the full density matrix D finally obtained by equation (26). Generally, most commonly used observables are one-body observables. These can be evaluated from the one-body density which is obtained as a trace of the N -body density over all but one particle, i.e.
Of particular interest for measuring equilibration is the one-body entropy S which is computed as (with k B the Boltzmann constant): in the natural orbital basis where the one-body density matrix is diagonal, i.e. ρ nm = δ nm ρ n .

A simple 1D model
To test the performances of the STDHF propagation, we use a simple 1D model. The mean field Hamiltonian is given (in x representation) as: where (α) (x) is the local one-body density associated to the actual Slater state |Φ (α) . The external potential V ext (x) is a Woods-Saxon profile V ext (x) = v 0 /[1 + exp((x − x 0 )/a)] with v 0 = −68 eV, x 0 = 15 a 0 , a = 2 a 0 . With such parameters, V ext represents a typical ionic background of a simple atom/molecule. The system is furthermore embedded into a confining harmonic oscillator tailored to act only outside the potential well (|x| > x 0 + 2a). This confining potential ensures soft reflecting boundary conditions. Closing the system this way avoids the conceptual complications of non-equilibrium thermodynamics in open systems. Note that as a result, absolute values of energies are somewhat arbitrary. What has been kept "realistic" is the structure of s.p. spectrum. The last term in equation (33), with strength λ, serves to simulate a typical self-consistent mean-field (HF or LDA) and its impact (relative to V ext (x)) can be tuned with λ. The two-body interaction behind the (α) 2 is a simple zerorange potential δ(x − x ). Actually, we use a moderate self-consistent term with λ = 27.2 eV a 2 0 . The set of chosen mean-field parameters leads to a spectrum of s.p. energies dense enough to allow a sufficient number of 2ph transitions. As written above, the residual interaction W ought to describe the in-medium scattering cross-section. We use here W (x, x ) = W 0 δ(x − x ) with W 0 = 18.3 eV. This choice leads to realistic relaxation times for the entropy. It is to be noted that we use W only for evaluating the collision term in the Markovian approximation which accounts exclusively for dissipation and thus does not interfere with the correlations embodied in the TDLDA propagation [113,114].
The considered system has N = 9 physical electrons. Without loss of generality, we do not consider spin effects here. The working space thus contains 9 hole levels complemented by an arbitrary number of empty levels, which is usually 14 in the following results (basically independent of this number). This makes the total number of holes and particles Ω = 23 here. The lowest orbital has an energy of −58.8 eV, while the HOMO is at −43.0 eV. The HOMO-LUMO gap is equal to 2.4 eV and thus corresponds to the minimal value of the excitation energy accessible by a particle-hole (1ph) transition.
The stationary state is obtained by solving the static mean-field equations with the Hamiltonian (33) which is independent of α at initial time, as STDHF is by construction unable to build static correlations in the ground state. We then initiate the dynamics by an instantaneous nph excitation. This delivers various initial excitations at different excitation energies E * . To excite lower energies, we use 1ph excitations. For larger energies, we step up to 2ph states and even further if more energy is wanted. This way to excite the system is very efficient for testing because it leads quickly to a sufficient number of level crossings between occupied and unoccupied levels and with it a sufficient number of 2ph transitions which match the δ Γ function in the jump rate (28c).
As already discussed in Section 6.1, one has to find a good compromise on the time interval for the stochastic jumps τ STDHF which should be at the same time rather small, to account for the dynamical correlations that develop in the course of time, and large enough to justify the use of an incoherent ensemble. For the results presented in the next section, we used τ STDHF between 0.1 and 2 fs. The second crucial parameter is the size N of the stochastic ensemble. It has to be sufficiently large to achieve proper statistics. On the other hand, it represents the major cost factor for STDHF. By construction, STDHF implies at least the price of N times mean field (TDLDA) propagation, not counting the evaluation of jumps at each multiple of τ STDHF . Experience shows that N = 100 is a minimum value to ensure a proper convergence of the results for the excitation energies under consideration here. Larger values of N (up to 300 or even 1000) become necessary when reducing the excitation energy. Last but not least, the width Γ of the δ function for the energy conservation in the jump probability (28b) is taken as Γ = 1.36 eV.

Results
To give an insight of how STDHF works, we show in the left panel of Figure 11 the time evolution of s.p. energies for the case of an initial excitation at E * = 25.2 eV. We compare results from a pure TDLDA propagation (upper left panel) with those from one sample of the STDHF ensemble (lower left panel). Both cases start from the same configuration. One clearly sees the 1ph excitation at initial time. In case of pure TDLDA, the s.p. energies fluctuate but maintain their ordering. By construction, occupation numbers are frozen at their original n α = 1 or 0 values (with an occupied state above empty ones and one deeply lying empty state). Instead, STDHF occasionally exploits one of the many level crossings between occupied and unoccupied levels which develop in the course of time. This generates situations in which the two energies in the jump probability (28b) match, thus allowing a jump between the two configurations as indicated by the dashed vertical lines in the lower left panel. The density of crossings is moderate in this rather small test case, but suffices to produce a sufficient amount of jumps to see equilibration towards a thermalized state.
Averaging over the different sequences from all samples of the STDHF ensemble produces a mixed state. The right panels of Figure 11 shows results deduced from this ensemble averaged state, the one-body entropy S and, as inserts, snapshots of the distribution of occupation numbers ρ n , see equation (32). The entropy shows nearly exponential convergence to an equilibrium value. The distributions of occupation clearly show the initial 1ph state and a gradual transition to a final Fermi distribution, of course with some fluctuations due to finite sampling.
Having seen in detail how STDHF evolution works and how entropy is produced in the course of time, we show in Figure 12 the time evolution of entropy for various excitation energy E * . Shaded areas around the curves represent the impact of numerical parameters, namely from variations of τ STDHF (in their range of relevance), Ω and N (which should be increased up to convergence). We see the typical pattern of a thermal relaxation process with an early rapid growth turning later to steadily slower increase. All patterns resemble an exponential relaxation towards the maximum entropy S lim for a given excitation energy E * , with relaxation times roughly decreasing with increasing E * and of the order of a few fs, representing realistic values for such a system [31].
The cases with lowest excitation energy show large fluctuations which renders the fitted relaxation times rather uncertain. These cases suffer from too low statistics because the dynamical phase space is small at low energies which, in turn, reduces jump probabilities too much. A brute force solution is to enhance the number of samples N . But this has its practical limitations. On the other hand, a low excitation energy means low fluctuations of the mean field. This allows us to reduce the expense by ignoring the fluctuations of the mean field and to use propagation in one common mean field. This leads to Averaged STDHF (ASTDHF) which is outlined in the next section.

Principle
The idea here is to replace the many different mean fields h (α) by one common mean field Hamiltonian h. This strategy is therefore coined "Average STDHF" (ASTDHF) [86]. If it is certainly justified at moderate excitations, one may hope that it still approximately holds at larger excitations, at least what concern average properties.
The common mean field is associated to one common basis of s.p. states which is used to span all D (α) . The numerical expense is thus considerably reduced at the side of mean-field propagation since only one set of s.p. states has to be handled. The D (α) built on a common set {|ϕ n } are distinguished only by the choice of which |ϕ n are actually occupied. This reduces the variety of D (α) dramatically. As a consequence, many D (α) , which emerge during STDHF propagation, will be identical. It is then preferable to recollect the total density matrix into a sum of different states Here, the weights W β are linked to those of the STDHF ensemble (26) as W β = N β /N , with N β being the number of identical samples in the STDHF ensemble covering N states. Equation (34) can thus be seen as the corollary of equation (26) but now with the difference that there is only one single set of s.p. states involved. However, it is not very efficient to unfold first the full STDHF ensemble and then to condense it to a common mean field. It is much more efficient to perform the reduction to an orthonormal basis in each stochastic step immediately.
The natural way to do so is to refer to the natural orbital basis which provides a diagonal representation of the one-body density matrix, see equation (19), and thus allows a direct evaluation of the associated mean field, as in the case of "pure" Slater states with integer (0 or 1) occupation numbers. The key of ASTDHF is to use the STDHF dissipative step for evaluating dynamical correlations. This means to create, at each dissipative step, an instantaneous STDHF ensemble of Slater states with integer occupation numbers, for which, in particular, transitions between hole and particle levels are well defined. One thus needs to map the actual one-body density matrix with its fractional occupation numbers into a set of Slater states with integer occupation numbers.
The pathway from a given N -body density D to the corresponding one-body density ρ is straightforward. The reverse relation is however much more involved. It is still simple and unique for pure Slater states |Φ (β) which are distinguished by the occupation numbers ν (β) n ∈ {0, 1} for the s.p. states |ϕ n . It then reads formally where A stands for anti-symmetrization and ρ i stands for the one-body density in the space of the ith particle. The ν n 's in equation (19) emerge then as ν n = β ν (β) n W β . The task of reverse mapping is to build a set of appropriate occupations ν (β) n with weights W β . We restrict the reverse mapping to the D (β) which have (approximately) the same energy as ρ, selected within a reservoir of ℵ Slater states through an energy criterion, using once again a finite width δ function (this width denoted by Γ map ). This reverse mapping corresponds to a non-negative least square (NNLS) problem [145] and is performed by minimizing the following quantity [146]: where η is some small positive numerical parameter. The term ∝ η serves to prefer the solution with the most evenly distributed coefficients W β . In practice, the minimization scheme proceeds until χ 2 < 10 −6 .

Time propagation in ASTDHF
The aim is to propagate the one-body density ρ(t) in natural orbital representation (19) which amounts to propagate the ϕ n (t) together with their occupations ν n (t). This is done in combination of mean-field propagation for the set {ϕ n (t)} and two-body collisions for changing the ν n (t). The latter are evaluated in practice by referring to the well developed strategy of STDHF jumps (Sect. 6.1). As in STDHF, time stepping is done at two time scales. The mean-field propagation is done on a fine time grid because it has to resolve the possibly high Fourier components of the mean field. The rearrangement of occupations by two-body collisions can proceed at slower pace and thus the time stepping of the ν n is done at a coarse time grid t j = jτ STDHF with j = 0, 1, 2, . . . as in STDHF. Starting from time t j−1 , the s.p. states are propagated from t j−1 to t j by solving TDLDA in standard manner with the common mean field Hamiltonian. This formally means a propagation |ϕ n (t j−1 ) → | ϕ n (t j ) . We use on purpose a specific notation for this pure mean field propagation with a "∼" to distinguish the thus propagated states from the forthcoming natural orbitals |ϕ n (t j ) at time t j which will account for collisions.
The collisional change of ν n is then performed at time t j by recourse to STDHF. To that end, the preliminary one-body density at t j is unfolded to many-body densities with weights W β , formally amounting to find appropriate sets of ν (β) n with W β which represent the ν n (t j−1 ). As in STDHF, each Slater state D (β) is now evolved to a correlated state κ with the jump rates as given in equation (28c) where κ[β] labels the space of 2ph states about the Slater state |Φ (β) plus the state as such for κ = 0 (Sect. 6.2). This yields altogether the new manybody density at t j as Now we have to recall that the STDHF step includes a redefinition of the s.p. basis to minimize energy uncertainty. The optimized s.p. basis thus depends again on β as usual in STDHF. But this is no obstacle for the further proceeding. We have to map anyway the new state D(t j ) to its one-body density. This is done in straightforward manner with equation (31). We finally diagonalize the one-body density matrix and obtain the density operator at time t j in natural orbitals with the final new ϕ n (t j ) and ν n (t j ), ready for the next coarse time step.

ASTDHF and quantum kinetic equation
ASTDHF replaces the ensemble of pure mean field trajectories by one mixed (in the statistical physics sense) dynamical evolution with a time evolution of occupation numbers of the one-body density matrix. This reminds standard quantum kinetic equations. As demonstrated in [3,77], the key to derive kinetic equations from STDHF is to gather the total STDHF ensemble in sub-ensembles, each characterized by one single mean field. By doing so, STDHF can be reduced to a stochastic kinetic equation, the quantum Boltzmann-Langevin equation. Dealing with a set of sub-ensembles allows one to keep a stochastic component in the picture and is responsible for the appearance of the stochastic collision term of the Boltzmann-Langevin equation, acting as a complement to the standard quantum collision term. In ASTDHF, by merging the set of sub-ensembles into one common ensemble with one single mean field, one cancels the stochastic collision term from the Boltzmann-Langevin equation which then reduces to a standard quantum Boltzmann kinetic equation. Since ASTDHF provides a realization of this sub-ensemble/total ensemble reduction, it actually allows a practical approach to solve the quantum Boltzmann equation.

Numerical details
At the side of the numerical cost of each scheme, AST-DHF is roughly only a factor 3-5 more expensive than pure mean field (TDLDA) propagation. It is thus typically 2 orders of magnitude less expensive than STDHF. Because of the high computational expense of STDHF, to allow a comparison with ASTDHF, extensive performance tests have been carried in the 1D model introduced in Section 6.3 to mimick a molecular system. Note however that, because of the order of magnitude reduction in computational cost when stepping from STDHF to AST-DHF, realistic computations in full 3D become accessible to ASTDHF, as the effort in ASTDHF mostly concerns the number of Slater states to be included at each coarse time step t j = jτ STDHF , an aspect which presumably little depends on the dimensionality of the coordinate space. Several numerical parameters have been explored to decide on the robustness of the method and on the cost of computations. The first one is the coarse time step τ STDHF , which has to be kept in a reasonable range (see [142] for detail). We have explored values of τ STDHF from 0.1 to 1 fs in ASTDHF (and up to 2 fs in STDHF, see Sect. 6.3). We use here the value τ STDHF = 0.242 fs, which is a typical time scale in clusters and molecules and dynamics is propagated on rather long times for such systems (120 fs).
The size N of the STDHF ensemble finds its equivalent in ASTDHF through the size ℵ of the reservoir of Slater states to be used to reconstruct the one-body density, see equation (35). It is controlled by the the finite width Γ map of the energy selection for the ASTDHF set. The latter has been varied from 1.4 to 12.2 eV. If one chooses too small a Γ map , the ASTDHF set of states will be too small and will then hinder the mapping from the one-body density matrix ρ to the N -body density matrix D. On the other hand, too large a value of Γ map produces too large an energy violation. We have checked that a good compromise between the size of the ASTDHF set and energy conservation is achieved with Γ map = 8.2 eV.
Finally, for the width Γ of the energy-conserving δ Γ function in the jump probability (28b), common in STDHF and ASTDHF, we again use the typical value of Γ = 1.36 eV which delivers energy conservation within 0.3% in both schemes.

Results in closed systems
The system is excited at initial time by a 1ph transition chosen from the ground state s.p. spectrum, following again the strategy presented in [78,142], and which has proven efficient for validating STDHF. Figure 13 shows typical results obtained in STDHF and ASTDHF, for the case E * = 20.7 eV. The left panel compares the occupation number distributions at 150 fs. They both resemble a Fermi-Dirac distribution. A fit provides a Fermi energy and a temperature of ε F = −39.43 eV and T = 5.74 eV for STDHF, and ε F = −40.37 eV and T = 5.81 eV for ASTDHF. Therefore, even if the distributions slightly differ, their thermal characteristics compare very well. This is even more visible when comparing the time evolution of the one-body entropy (right panel) which also shows that ASTDHF and STDHF display very close relaxation pattern in terms of time scales.
We have pursued the comparison by studying the dependence of the asymptotic value of the entropy, S lim , as a function of the excitation energy E * . In a simple Fermi gas, one expects that the square of the equilibrium entropy S equil ≈ S lim scales as √ E * . It was shown in [78] that STDHF indeed fulfills this behavior. Figure 14 presents S lim 2 as a function of E * for STDHF and ASTDHF. Both follow a linear correlation (see dashed and dotted lines), with slightly different slopes (2.24 for ASTDHF and 2.14 for STDHF). In addition, the ASTDHF results almost perfectly superimpose the STDHF ones, up to the error bars accounting for variations of model parameters, similarly as discussed for STDHF (Sect. 6.4). We however observe that the ASTDHF asymptotic entropy is systematically above that obtained in STDHF. This is however accidental : we have checked that, by changing the shape of the external potential defined in equation (33), and therefore the level spacing, the agreement between STDHF and ASTDHF remains the same but their relative position does depend on the actual spectrum. It is thus not a significant feature at that level of detail.

Open systems
We have focused, up to now, all STDHF (Sect. 6.4) and ASTDHF (Sect. 7.5) applications to closed systems. This allowed us to avoid competition between ionization and thermalization and simplifies the analysis. This competition was, in turn, analyzed in the RTA framework (Sect. 5). We now want to examine how ASTDHF performs in the case of an open system. It should be noted here that at the core of STDHF and ASTDHF approaches lie well identified transitions between s.p. states.
Stepping to an open system means to use absorbing boundary conditions on top of usual dynamics of s.p. wave functions (Sect. 2.2). This has the severe consequence that the transitions may take place between s.p. wave functions not entirely contained in the computational box. This directly impacts the fundamental inputs entering STDHF and ASTDHF, in particular in the meaning of the notion of particle-hole transitions to be computed between non orthonormal states (within the computational box). We have thus been led to develop a dedicated strategy to correct such a technical issue, in order to preserve, within ASTDHF, the clear-cut definitions of particle-hole transitions with the associated energy conservation. This led us to split any wave function ϕ n into two parts, namely ϕ in n and ϕ out n , representing the part of the wave function inside and outside the computational box respectively.
The practical implementation of this framework is rather involved and uses several s.p. bases in parallel, following a strategy somewhat similar to that adopted to analyze TDLDA dynamics within the natural orbital basis [147]. Details of this approach will be found in [148]. We thus only illustrate here the capabilities of the method on an example in our 1D framework (Sect. 6.3). Results are shown in Figure 15. At variance with previous particle-hole excitations (Sects. 6.4 and 7.5), and in order to generate a simple dipole oscillation and associated ionization, we excite the system by an initial instantaneous boost, which mocks up a very short (possibly violent) laser pulse (Sect. 2.3). As expected, the AST-DHF approach (similar to RTA) provides a significant damping of the dipole moment, visible in the lower panel of Figure 15. This is associated in that case to a sizable reduction of ionization (upper panel). This clearly shows that ASTDHF can be used in realistic irradiation scenarios where there is a competition between direct ionization and thermalization.

Conclusions
We have presented in this paper several directions of research aiming at extending TDDFT at the level of the TDLDA by the inclusion of dissipative effects, which become important in electron dynamics at high excitations. The description of dissipative mechanisms in far-off equilibrium dynamics is a long standing question. It raises major formal and practical difficulties, especially in the regime of moderate excitations where quantum effects are still playing a key role. One thus needs to develop a truly quantum-mechanical theory of dissipation on top of a robust time-dependent mean-field description, as TDLDA has proven to be. We have investigated several possible strategies going into that direction.
We first considered TDCDFT as an extension which is designed to provide a proper description of damping in the linear regime. It turned out, at least within the limitations of the functional used in our tests, that TDCDFT is unable to simulate correctly far-off equilibrium situations because it does not reproduce the strong energy dependence of relaxation time as is typical for highly excited systems.
We then stepped to semi-classical approaches, here in particular the VUU equation. They provide a detailed description of two-body collisions including the basic quantum effect of Pauli blocking and they are very successful in the realm of violent fermion dynamics, particularly well suited for Fermi liquids as nuclei and electrons in metals. But the underlying Vlasov dynamics lacks all quantum shell effects which become crucial at moderate excitation energies and/or for small systems.
All further steps consider quantum mechanical versions of a kinetic equation. This still remains a formidable task in finite quantum systems. A still simple and manageable scheme is provided by combining the RTA with TDLDA for finite systems. A crucial input is here the local relaxation time for which we adopt the well tested semi-classical approximation. RTA makes it possible to consider realistic irradiation scenarios irrespective of the nature of the material considered. This allowed us to make systematic investigations of laser irradiation scenarios. We found, in particular, that dissipation is especially active if the laser frequency is in resonance with system's modes, thus, for metal clusters, particularly in the regime of the surface plasmon resonance. This complies with the findings of semi-classical approaches, but confirms exactly the urgent need for a quantum description because only that can guarantee a proper account of spectral density.
As a next step upward in refinement, we considered a stochastic extension of mean-field dynamics, coined STDHF. This corresponds to a full quantum kinetic equation complemented by a collision term which is evaluated with stochastic methods. This very elaborate approach treats the dynamics in terms of an ensemble of Slater states and can account for mean field fluctuations. It thus allows one to address even complex dynamical scenarios including bifurcation and multi-fragmentation. Nevertheless, it requires a heavy computational effort and we thus use in the present development stage a simplified 1D molecular model to test it. With that, we have shown that STDHF provides a proper account of thermalization. A practical limitation of the full STDHF treatment appears in the regime of low excitation. The collision rate shrinks which, in turn, requires even larger ensembles to collect sufficient statistics. On the other hand, we see that fluctuations of the mean field remain small in this regime. This opens the path to a simplification, which deals with one common mean field for all states in the ensemble, an approach coined ASTDHF. The common mean field provides a common s.p. basis which allows a much more complete book-keeping of collisions. It amounts eventually to a mere (fully quantal) kinetic equation. Thus ASTDHF provides a valuable alternative to STDHF in the regime of moderate excitations (small fluctuations) whose reduced expense allows more easily an extension to full 3D treatment. Finally, we have sketched the most recent results which focus, at the level of ASTDHF, on the interplay between dissipation and ionization, a line of development still underway.
What is still missing is a modeling of very short times where coherent correlations (ground state and dynamical ones) play a role. There are several theories offering such a possibility, but all limited to very short times.
These need yet to be interfaced with (A)STDHF to properly bridge the transition from the coherent to the incoherent regime. For instance, one can envision to start from a more correlated ground state, as that delivered by multi-configurational Hartree-Fock or reduced density matrix theory, and compare a fully coherent time evolution with the incoherent evolution delivered by (A)STDHF. Work along that line is in progress. This work was supported by the CNRS and the Midi-Pyrénées region (doctoral allocation number 13050239), and by the Institut Universitaire de France. This work was granted access to the HPC resources of IDRIS under the allocation 2014-095115 made by GENCI (Grand Equipement National de Calcul Intensif), of CalMiP (Calcul en Midi-Pyrénées) under the allocation P1238, and of RRZE (Regionales Rechenzentrum Erlangen).