Phuong Mai Dinh
email: dinh@irsamc.ups-tlse.fr

Marc Vincendon

François Coppens

Eric Suraud

Paul-Gerhard Reinhard

P.-G Reinhard

Quantum Dissipative Dynamics (QDD): A real-time real-space approach to far-off-equilibrium dynamics in finite electron systems

Keywords: time-dependent density functional theory, electronic dissipation, photoelectron spectra, electron emission, beyond mean field

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The theoretical description of dynamics far off equilibrium in quantum many-body systems has been for decades a widely open problem in numerous fields of science. Nonequilibrium and dissipative quantum many-particle dynamics has a long history which can be traced back all the way to Bohr's early pioneering work on charged-particle penetration and stopping in matter [1,2]. The nuclear physics domain has been explored over the last four decades in connection with the rich empirical material from heavy-ion collisions. This was done mostly with classical or semi-classical approaches [3,4,5,[START_REF] Durand | Nuclear Dynamics in the Nucleonic Regime[END_REF][START_REF] Feldmeier | [END_REF]8,9]. Non-equilibrium dynamics has also become a key issue in transport processes in solids [10]. It is an equally important issue for ultracold bosonic and fermionic gases ("quenches") [START_REF] Bonitz | Quantum Kinetic Theory[END_REF][START_REF] Streltsov | [END_REF]13,14] as well as for the electronic dynamics in atoms, molecules, and condensed matter driven by ultrashort and strong laser fields [15,16,17,18,19,20,21,22,23,[START_REF] Haken | Time-Dependent Density-Functional Theory: Concepts and Applications[END_REF][START_REF] Wilken | [END_REF]. Multi-differential observables such as vectorial single-or multipleelectron emission patterns or density-density correlation functions and phase shifts of many-particle wave functions [14] can now be experimentally probed in unprecedented detail and provide benchmarks for state-of-the-art theories.

The electronic case is especially challenging in finite systems with its numerous irradiation scenarios by lasers or charged particles. Early stages of irradiation dynamics proceed at short time scale, predominantly at electronic level. For example, proper analysis of ultrafast processes in clusters and molecules is crucial for understanding microscopic mechanisms underlying radiation damage as manifested in tumor cell destruction in oncology [START_REF]Nanoscale Insights into Ion-Beam Cancer Therapies[END_REF] or in shielding electronic devices [START_REF] Autran | Soft Errors: From Particles to Circuits[END_REF]. Irradiation by charged particles (swift ion, electron) can deliver a violent electromagnetic kick to the system driving it far off equilibrium. Laser irradiation with the high versatility to design pulse profiles adds new possibilities. It is now possible to explore sub-fs (1 fs = 10 -15 s) dynamics down to the attosecond (1 as = 10 -18 s) domain, hence opening the door to fully time-resolved analysis of irradiation effects at the electronic level [START_REF] Krausz | [END_REF]. Photo-induced ultrafast processes are also involved in many other mechanisms such as vision, photo-synthesis or solar cells. Up to now, most investigations focused on ionization of atoms [START_REF] Krausz | [END_REF] and only recently on small molecules [29]. Forthcoming experiments will address ultrafast processes in complex systems such as e.g. C 60 . A challenging aspect lies here in the expected impact of electronic correlations [30] possibly driving relaxation pathways prior to decoherence and ionic relaxation. This implies both the development of new dedicated (coincidence [31]) experiments using elaborate observables and new theories properly including correlations.

Understanding violent irradiation processes requires to describe far-off-equilibrium dynamics, involving ionization, electron transport and strong electron dynamical correlation effects leading to dissipation and thermalization. Note that the term "dynamical correlations" is widely used and, unfortunately, with different meanings in different contexts. Here we associate with it specifically those correlations which are genuinely driven by the dynamics of the system and cannot be incorporated into standard approximations of TDDFT propagation (even not effectively). More specifically, we confine them here to that part of time-dependent correlations which can be described in Markov approximation which pick exactly the dissipative part of correlations [32] and which is often imagined as result of particle-particle collisions [4]. The role of dynamical correlations unfolds successively in the course of dynamical evolution. The immediate electron response to irradiation kicks the system far off equilibrium but leaves it still fully quantum. Initially coherent quantum correlations dominate at very short times. But sooner or later, the phase space for electrons collisions progressively opens up leading to a dominance of the incoherent dynamical correlations tractable with tools of statistical physics. On even longer time scales, electrons couple to ionic degrees of freedom (molecular vibrations, dissociation). The intermediate time window (from coherent excitation to statistical electron correlations) is the pivotal link from microscopic excitation to long-time evolution, which also gives access to larger systems then unfolding to macroscopic scale.

Today, the most widely used microscopic theoretical approaches addressing such intermediate time scales in realistic systems rely on Time-Dependent Density Functional Theory (TDDFT) [START_REF] Haken | Time-Dependent Density-Functional Theory: Concepts and Applications[END_REF]. They are efficient and reliable at early times and/or low excitation. However, they do not include dynamical correlations and associated dissipative effects which come into play at intermediate times. It is the aim of this paper to describe the first open source package based on electronic TDDFT augmented by dynamical correlations and dissipative effects for applications in finite electronic systems such as atoms, molecules and clusters. To indicate the combination of quantum treatment and dissipative dynamics, the code package is named Quantum Dissipative Dynamics (QDD). As often done with TDDFT, electrons are coupled to ionic motion described by classical molecular dynamics [33,[START_REF] Reinhard | Introduction to Cluster Dynamics[END_REF]. A word of caution is in order. This way of coupling to ionic dynamics, coined Ehrenfest dynamics, misses proper energy transfer at even longer times scales. This can only be cured by ionic correlations [START_REF][END_REF] which, however, is far beyond the scope of the present project. 4 41 and C 60 . We show ionic structure and electronic cloud, sequence of single particle energies (with the HOMO emphasized in blue and the LUMO in red dashes) and the optical response in some specific spatial directions as indicated.

Typical examples of systems attainable in this package are demonstrated in figure 1. We show there the structure properties of four typical examples of clusters and molecules which correspond to various system's sizes and kinds of binding [36]. We show for all of them the ionic structure dressed by the associated electron cloud. More details on the ground state properties are accessible via the sequence of single electron states on the left side of each panel. Two specific levels are especially important for dynamical applications, the Highest Occupied Molecular Orbital (HOMO) and the Lowest Unoccupied Molecular Orbital (LUMO) one. They are respectively plotted in blue (HOMO) and red (LUMO) while other occupied levels are black. We also show the optical response of each system in the bottom part of each panel. This is a crucial system property (optical response) which characterizes the coupling of the electrons cloud to photon fields. We compute it directly in real time from the electronic dipole response Fourier transformed to frequency domain (section 3.2.4). Note that the computed ground state properties of the four presented cases reproduce rather well the experimentally known data.

The paper is organized as follows. Section 2 gathers formal theoretical aspects underlying the theory. We describe numerical methods and tools in section 3 with details on computed observables and outputs. Section 4 provides a few typical examples of applications in various systems (clusters, molecules). Conclusions and perspectives are given in section 5. To complement this document, we provide a supplementary material that should be taken as a user manual. At several places, cross-references to the user manual can be found; they are indicated by a heading letter "M" to avoid any confusion with a section or an equation of the present paper.

2. Formal background: electronic DFT coupled to ionic motion

Brief review: mean-field theory and correlations

In the ground state of a system, correlations are defined as deviations from mean field, such as Hartree-Fock (HF). In electronic systems, numerous well established methods treat ground-state correlations at various levels such as Density-Functional Theories (DFT) [START_REF] Dreizler | Density Functional Theory: An Approach to the Quantum Many-Body Problem[END_REF], wave function-based methods as Configuration Interaction or Coupled-Cluster methods [START_REF] Schmidt | [END_REF][START_REF] Jensen | Introduction to Computational Chemistry[END_REF], Green's functions [START_REF] Martin | Interacting electrons: theory and computational approaches[END_REF], etc. Much less developed are descriptions of correlations in the time domain, which range from particle-particle correlations in electron emission to energy relaxation towards thermal equilibrium. As with DFT for the ground state, Time-Dependent DFT (TDDFT), mostly treated at the level of the Time-Dependent Local-Density Approximation (TDLDA), approximately accounts for some correlations and so provides a reliable and efficient starting point for the real-time simulation of electronic dynamics in large systems [START_REF] Marques | Time Dependent Density Functional Theory[END_REF][START_REF] Marques | Fundamentals of Time-Dependent Density Functional Theory[END_REF][START_REF] Shinohara | [END_REF][START_REF] Haken | Time-Dependent Density-Functional Theory: Concepts and Applications[END_REF]44,45,46]. TDLDA practically carries forth the ground-state correlations from LDA assumed to follow instantaneously the changing situations. This is why this approach is also called Adiabatic LDA [START_REF] Gross | Density Functional Theory[END_REF]. Still, at high excitation energies and long times, new correlations from electron-electron collisions build up. These "Dynamical Correlations" (DC) correspond to deviations from time-dependent mean field, here described by TDDFT at the level of TDLDA.

The theoretical treatment of DC in finite Fermion systems far off equilibrium is still under development. Fully coherent descriptions of DC are often developed from timedependent extensions of correlated ground state approaches. Wave function-based techniques such as Multi-Configurational Time-Dependent Hartree-Fock (MCTDHF) [START_REF] Kitzler | [END_REF]15] or Time-Dependent Complete Active Space Self-Consistent Field (TD-CASSCF) [23], can mostly access DC for short times (before dissipation sets in) and small systems [17,[START_REF] Stefanucci | Nonequilibrium many-body theory of quantum systems[END_REF]19], giving little hope to extend them directly to large systems and/or long times, because of their too demanding computational scaling properties. The Time-Dependent 2-particle Reduced Density Matrix (TD2RDM) method delivers a close-to-exact solution to the problem [START_REF] Lackner | [END_REF], with better computational performances, although again very limited (small systems ∼ 10 electrons, short times ∼ 10 fs). On the other side are semiclassical and even fully classical approaches. These simplify the rather involved DC by dismissing coherent correlations and quantum features, thus delivering a robust approach applicable to large systems and high excitation energies which often works reasonably well [51,52,53,54,55]. Of course, they leave a large gap between the true quantum regime of low to moderate excitations and the regime of high excitation density by collisions leading to relaxation and thermalization. They thus miss many irradiation scenarios at intermediate energies/early stages dominated by quantum effects [56].

To fill this gap, we have recently developed a set of quantum theories, coined QDD (Quantum Dissipative Dynamics) in which DC are treated incoherently while maintaining the quantum structure throughout the process [57]. These theories can deal with large systems (> 100 active electrons) on intermediate to long times (> 100 fs). They are based on Stochastic TDHF (STDHF) [58] and different levels of approximation to it. STDHF describes DC in terms of a stochastic treatment of electron-electron collisions leading to an incoherent ensemble of time-dependent mean-field states. The quantum collision term is evaluated in Markov approximation, thus neglecting memory effects. The latter depend on the spectral density and usually decrease with excitation energy [59]. STDHF has been validated (at high enough excitation energies) on an exactly solvable model [60]. At moderate excitation, STDHF can be reduced to an average version ignoring the fluctuations of the mean field [61]. A further simplification consists in approximating the collision term by a Relaxation-Time Approximation (RTA) [62]. Whatever approach, dissipative dynamics thus neglects coherent DC. Coherence may play a role at the very first stages of building up DC, the longer the lower the excitation energy. STDHF, ASTDHF, and particularly RTA, are thus best suited to sufficiently strong perturbations (as delivered by a swift charged particle or a fs laser) and analysis of long time spans with a bias on relaxation toward equilibrium. A detailed description of pre-equilibrium processes stays outside the scope of these approaches.

In the present paper, we present our QDD code, an open source implementation of the RTA on top of a 3D coordinate-space code for TDLDA augmented by a various versions of Self-Interaction Correction (SIC) [63]. TDLDA is formulated at standard Kohn-Sham level [64] using single particle wave functions as building blocks. Dissipative features at RTA level are formulated in terms of the one-body density matrix built from a large set of single particle wave functions with fractional occupation numbers. Ionic degrees of freedom are treated classically in a standard manner, as usually done in realtime implementations of TDDFT [START_REF] Marques | Time Dependent Density Functional Theory[END_REF][START_REF] Marques | Fundamentals of Time-Dependent Density Functional Theory[END_REF]. We formulate TDLDA starting from the total energy (section 2.2) from which stationarity of action delivers i) standard Kohn-Sham TDLDA equations for single electron wave functions and ii) classical Molecular Dynamics (MD) equations for ions (section 2.3). This altogether leads to standard combined TDLDA-MD. RTA is built on top of the TDLDA evolution by making occupation numbers time-dependent according to the RTA kinetic equation (section 2.3.2). RTA is thus coupled to ionic MD in the same standard manner as TDLDA.

The total energy of the model

The basic dynamical variables in mean-field theory and in QDD are the set of single particle (s.p.) wave functions with their occupation numbers at the side of valence electrons, and classical coordinates and momenta for the ions: s.p. wave functions: ϕ α , α = 1...Ω , s.p. occupation numbers:

w α , α = 1...Ω , ionic coordinates: R I , I = 1...N ion , ionic momenta:

P I , I = 1...N ion . (1)
Note that α also accounts for spin (up, ↑ and down, ↓) degrees of freedom. Ω denotes the total number of s.p. wave functions. It must be Ω ≥ N el for a system with N el valence electrons. Later on, we will denote by Z I the charge of ion I and by M I its mass. The term valence electrons needs explanation: Electrons in a cluster or a molecule have very different energetic properties. For the kind of excitation that can be dealt with in TDLDA and QDD, we only need to handle the least bound electrons which are commonly coined "valence" electrons. The "core" electrons remain practically inert because they are energetically too far away. They form together with the corresponding atomic nucleus an ionic core. The coupling of the core electrons to the valence electrons is described in terms of pseudopotentials, see section 2.2.3.

A key quantity in connection with density functionals is the electronic local density which covers, in fact, two separate densities for spins up and spins down

↑ (r, t) = α∈↑ w α ϕ * α (r, t)ϕ α (r, t) , ↓ (r, t) = α∈↓ w α ϕ * α (r, t)ϕ α (r, t) , (2a)
which sum up to the total density (r, t) = ↑ (r, t)

+ ↓ (r, t) = α w α ϕ * α (r, t)ϕ α (r, t). (2b)
2.2.1. The energy in Local-Density Approximation (LDA) Starting point is an expression for the total energy of the coupled electronic and ionic system:

E total = E kin + E C + E xc ELDA +E ext + E el,ion + E kin.ion + E pot,ion , (3a) E kin [{ϕ α }] = dr α w α ϕ + α p2 2m e ϕ α , (3b)
E C [ ] = e 2 2 dr dr (r, t) (r , t) |r -r | (3c) E xc [ ↑ , ↓ ] = dr ↑ (r, t) xc ( ↑ (r, t)) + dr ↓ (r, t) xc ( ↓ (r, t)) ,(3d) E ext [ , t] = dr (r, t)V ext (r, t) , (3e)
E el,ion [{ϕ α }, {R I }] = I dr α w α ϕ + α VPsP (r -R I )ϕ α , (3f)
E kin,ion (P I ) =

I P 2 I 2M I , (3g)
E pot,ion (R I ) = 1 2 J =I e 2 Z I Z J |R I -R J | + V ext,ion (R I , t) , (3h)
where functionals of density are indicated by square brackets, and functions of coordinates by round brackets. E C is the direct part of the electronic Coulomb energy. E xc is the density functional for electronic exchange and correlations computed from the exchange correlation functional xc . The sum of E C and E xc constitute the LDA energy E LDA . We use in the code two options: the Perdew-Wang functional of [65] or the older Gunnarsson-Lundqvist functional from [66]. These two options correspond to a Local Density Approximation (LDA). We furthermore use the stationary expression of the functionals in dynamical simulations (adiabatic approximation). The LDA exchange-correlation term then becomes a function of the density: E xc [ ] → E xc ( ); and the resulting DFT approach becomes an effective mean field theory.

E ext stands for the excitation mechanisms by external sources, either from a laser pulse or from the Coulomb field of a fast bypassing ion, see section 2.2.5.1. This part is, of course, absent in static calculations of the ground state. E el,ion carries the interaction of the electrons with the ions which is usually described by a pseudopotential VPsP , see section 2.2.3, or may be simplified in terms of the jellium model, see section 2.2.4. Ionic kinetic and potential energy are described by the obvious classical expressions. V ext,ion in the ionic potential energy describes the action of an external field on the ions.

Self-Interaction correction (SIC)

The LDA approximation introduces a self-interaction (SI) error because the total energy employs the total density which also includes the electron on which the interaction acts. This causes a couple of problems and correspondingly there were many attempts to mellow the SI errors, as the generalized gradient approximation [67], a kinetic extension thereof [68], or hybrid xc functionals see, e.g., [69]. All these attempts are of limited use for dynamics of electron emission where the SI error is particularly harmful. Its effect for the long-range Coulomb term is to shift the s.p. energy spectrum. This leads to a wrong ionization potential if estimated from the energy of the HOMO. This is a severe problem in a real-time simulation in which electron emission is evaluated on the fly. One must remove the SI error exactly. A safe way to achieve that is to augment the LDA density functional an SI Correction (SIC) [70]

E LDA -→ E LDA [ ↑ , ↓ ] - α E LDA [ α ] ( 4
)
where α labels the density associated to s.p.

α, that is α (r, t) = |ϕ α (r, t)| 2 .
In this formulation of SIC, the price to pay is to deal with a non-hermitean and state-dependent Hamiltonian (see section 2.3). The difficulty is especially severe in the time domain as the non-hermiticity of the Hamiltonian requires extra measures to guarantee unitary propagation, for tractable strategies see [71,72].

Thus, a solution of full dynamical SIC equations remains computationally expensive and there is still no universally accepted solution to the SIC problem, especially in the time domain. Fortunately, there are many situations in which simplified implementations of SIC are possible and justified. The Average Density SIC (ADSIC), which was proposed already in the 1930s [73] and applied since the 2000s in clusters [63], provides a simple and robust SIC. The idea underlying ADSIC is to assume that all valence electrons contribute about equally to the SI error. This amounts to replace Eq. ( 4) by

E LDA -→ E LDA [ ↑ , ↓ ] - σ∈{↑ ↓} N σ E LDA [ σ /N σ ] . ( 5
)
where N σ is the number of electrons with spin σ. ADSIC delivers again a functional of the local density only and thus can be treated in the same manner as any LDA scheme. It turns out that ADSIC works remarkably well in a wide class of compact atomic/molecular systems [74] and especially well for metal clusters. The major difficulty with ADSIC lies in the total number N el of electrons which explicitly enters the functional. It is thus not applicable to situations with fragmented electron density such as in dissociation and to bulk systems (the latter because N el grows infinite). Ionization is manageable as long as it remains moderate as compared to N el . And this is typically the physical situations we are interested in. There is also another case for which ADSIC is probably not suited, namely if the system combines various types of binding. This is the situation for instance in Na(H 2 O) n complexes in which metallic and covalent bonding are at play [75]. Such cases require a full SIC treatment which is feasible with some formal and algorithmic effort [71,72]. However, this is much more elaborate and more involved than ADSIC. The present release of QDD does not contain full SIC. It is, nevertheless, our plan to include it in a future release.

There exist more elaborate approaches to full SIC. One option is to treat Coulomb exchange exactly, but then ignoring correlations. The other options belong to the class of optimized effective potentials (OEP) which aim at best possible local potentials simulating exchange-correlation effects [76,77,78]. The next simple after ADSIC is the SIC-Slater approximation [79,80] which accounts for the local dominance of particular s.p. states. The method of Krieger-Li-Ifrate (KLI) puts a feedback loop on top of SIC-Slater to cure deficiencies of the mere SIC-Slater approach [81]. These three methods, SIC-Slater, KLI, and exact exchange, are also implemented in QDD, however with reduced functionality: they do not work for fractional occupation numbers. A safe application of SIC-Slater and KLI is confined to static calculation because time-dependent simulations can easily run into insurmountable difficulties [82,83]. Nonetheless, we find it useful to have these options available for occasional checks. For the formal details, see the brief summary in [63] or the extensive review on OEP [78].

After all, full TDDFT dynamics which is free of the SI error and applicable to all systems is still a matter of research and development. To say it positively, the present code, limited to ADSIC in that respect, provides a useful tool for many applications and with a bit of coding for further developments of a time-dependent SIC.

Pseudopotentials

There is a great variety of pseudopotentials available. The code employs two variants. Particularly efficient and simple to use are local pseudopotentials integrated from a pseudo-density which is represented as a sum of Gaussians. The corresponding potential is then a sum of error functions

V PsP (r) = Ze 2 2 i=1 c i erf(|r|/( √ 2σ i )) |r| , (6a) erf
(x) = 2 π x 0 dy e -y 2 , (6b)
where the σ i are widths and the strength parameters c i have to line up to the total charge of the ionic core c 1 + c 2 = Z ion . This pseudopotential is well suited for alkaline atoms for which it was originally developed [84]. More involved, but also more versatile in the applicability, are the Goedecker-like pseudopotentials from [85] which are composed from a local and a non-local part as

VPsP = V loc (r -R) + Vnloc , (7a)
V loc (r -R) = - Z ion x erf (x) + e -x 2 3 n=0 C n+1 2 n x 2n , x = |r -R| r loc , (7b)
Vnloc (r, r , R) = i,j l l m=-l G * ilm (r -R) h l ij G jlm (r -R) , ( 7c
)
where R is the position of the ionic core with respect to which the pseudopotential is defined. The non-local part (7c) serves to project out the electronic states which are occupied in the ionic core. It involves a summation over the angular momentum l, its component m = -l...+l, and i, j. The projector functions G are defined as:

G ilm (r) = √ 2 |r| l+2i-2 exp - |r| 2 2r 2 l r l+(4i-1)/2 l Γ l + 4i -1 2 Y lm (Ω r ) , ( 7d
)
where Γ is the gamma function and Y lm is the spherical harmonics, evaluated at the solid angle Ω r about the position vector r. For the chemical elements we deal with, we use orders i, j = 1, 2 and l runs from 0 up to 2. The parameters related to a given ionic core are thus: r loc , C 1...4 , r l , h l ij with i, j ∈ {1, 2}. The projection looks expensive at first glance. However, one can exploit the fact that the Gaussian projectors cover only a small region of space, of order of a few multiples of the non-local radii r l . Thus one needs to evaluate the projectors only on a small sub-grid which reduces the expense dramatically.

The soft jellium model for the ionic background

The electronic wave functions of bulk metals and of metal clusters are spread softly over the whole system and hardly resolve the spatial structure of the ionic cores. This allows one to approximate the detailed ionic background by a smooth, positive background density. Such a jellium approximation is a standard approach in the theory of bulk metals [START_REF] Ashcroft | Solid State Physics[END_REF]. It is also very popular in connection with metal clusters because it allows one to describe the gross features of cluster electronic shell structure and dynamics without bothering about the fine details of ionic structure. The adaptation of the jellium model to a finite cluster is straightforward. One carves from bulk jellium a finite element of constant positive charge density corresponding to the average bulk density. The volume is chosen such that its total charge coincides with the given ionic charge. For finite clusters, it is advantageous to use jellium with a soft surface profile. This is more suited for numerical handling and it improves the quality of the model, e.g., by producing a correct peak energy for the optical response of metal clusters [START_REF] Rubio | [END_REF]. Versatile and easy to handle in this respect is a Woods-Saxon profile for the jellium density

jel (r) = jel,0 1 + exp |r| -R(θ, φ) σ jel -1 , ( 8a
) with R(θ, φ) = R jel 1 + lm α lm Y lm (θ, φ) , (8b)
dr jel = N ion . (8c)
The jellium radius R jel is determined by Eq. (8c), the normalization to the desired number of ions. The central density is determined by the bulk density jel,0 = 3/(4πr 3 s ) and the Wigner-Seitz radius r s is a genuine material parameter [START_REF] Ashcroft | Solid State Physics[END_REF]. σ jel parametrizes the surface width and the transition from 90% to 10% bulk density is achieved within 4σ jel .

The model also allows one to describe deformed clusters by angular dependence R(θ, φ) of the extension. This is achieved through the deformation coefficients α lm weighting the impact of the spherical harmonics. Actually, we use in the code only the quadrupole deformations α 20 and α 2±2 as well as the hexadecapole deformation α 40 . Axially symmetric shapes are obtained if α 2±2 = 0 where then positive α 20 produce prolate shapes and negative α 20 oblate ones. The hexadecapole moment plays a role for fine-tuning the shape [88]. The cluster radius R jel is fixed by the other parameters through Eq. (8c) which is to be solved numerically by a root finding procedure, see the Supplemental material, section M.3.2.2. After all, the leading parameters of the soft jellium model ( 8) are the Wigner-Seitz radius r s and the surface thickness σ jel . They are universal for a given material. Typical values are r s ∼ 4 a 0 with σ jel ∼ 0.9 a 0 for Na clusters, r s ∼ 2.66 a 0 and σ jel ∼ 0.76 a 0 for Mg clusters, or r s ∼ 3 a 0 and σ jel ∼ 0.78 a 0 for Ag clusters. The deformation parameters α lm depend on the actual cluster and strongly vary with system and size.

The jellium approximation then consists in discarding the ionic contribution to the total energy (3), i.e. the terms E kin,ion and E pot,ion , and to replace the pseudopotential background in E el,ion by the Coulomb potential of the jellium density (8). There is no dynamics associated with the jellium. The model thus applies to situations where ionic motion can be ignored. Note that QDD offers the possibility to consider an even simpler modeling option by means of phenomenological electronic shell models of harmonicoscillator or Woods-Saxon type [89,[START_REF] Reinhard | Introduction to Cluster Dynamics[END_REF]. This can be done in very flexible manner by supplying the shell model choice in form of an explicit 3D field and reading it as background potential, for details, see the Supplemental material, section M.3.3.

External excitation fields and excitation mechanisms

The total energy (3) contains a contribution E ext from (time-dependent) external electromagnetic fields. These are supposed to come from external sources and they serve as dynamical excitation mechanisms. The code contains explicit entries for laser pulses or the Coulomb field of a by-passing charged ion. The excitation reduces to an instantaneous dipole boost in the limit of extremely short pulses or a very fast and sufficiently far by-passing ion. This constitutes an alternative way of excitation which we will also discuss at the end of this subsection, together with two further instantaneous excitation schemes for more subtle modes (spin dipole, scissors). This type of excitation, which changes instantaneously the electronic center-of-mass momentum, is also often called a kick excitation.

Laser fields

Lasers are the most important and very flexible means for a dedicated, well tuned excitation of electronic systems. They produce a strong coherent electromagnetic field which can be well approximated by a classical time-dependent electromagnetic field. Typical wavelengths are in the range of several hundreds of nm. This is a huge distance as compared to the spatial extension of atoms, molecules, and (most) clusters. One can thus treat the laser field in the limit of long wavelengths (k → 0). This amounts to deal with a spatially homogeneous electrical field E at the cluster site and we can also neglect the effect of the magnetic field for the laser intensities of relevance here. The coupling Hamiltonian leaves the freedom of gauge transformation [START_REF] Jackson | Classical Electrodynamics[END_REF]. The external field operator in velocity gauge reads:

V ext = eE 0 F (t) • p , F (t) = t 0 dt f (t ) exp (-iω las t ) , (9)
where f (t ) is the envelope of the laser pulse. The same pulse in length gauge becomes

V ext = eE 0 f (t)•r exp (-iω las t) (10)
which is simpler to handle because the laser field acts here simply as a time-dependent local operator. QDD thus uses the external field in length gauge. Transformation to velocity gauge, if needed, can be performed a posteriori by standard rules of gauge transformation [START_REF] Messiah | Quantum Mechanics[END_REF].

The laser pulse is characterized by frequency ω las , peak field strength E 0 = |E 0 |, polarization E 0 /E 0 , and time profile f (t). The field strength is usually parametrized in terms of the laser intensity I as

E 0 = c EI I 1/2 , c EI = 1.07 × 10 -8 eV a 0 W cm 2 -1/2 . ( 11
)
Even if numerous works on the characterization of the time profile of the laser pulse do exist [START_REF] Trebino | [END_REF]93,94,95], this profile is very often not precisely known experimentally, and is usually taken as a Gaussian with a certain Full Width at Half Maximum (FWHM).

Gaussians require long simulation times to cover the outer wings of the pulse. To avoid that, we shape a pulse with finite support by using a sin 2 envelop instead:

f (t) =    sin 2 π t 2T pulse for t ∈ {0, 2T pulse } 0 else . ( 12
)
It is close to a Gaussian pulse in the vicinity of peak field strength and combines high spectral selectivity with finite bounds. Note that the form ( 12) is scaled such that the pulse parameter T pulse is identical with the FWHM of a Gaussian pulse.

Charged projectiles

Probing electronic systems by beams of charged particles is a standard tool in atomic and molecular physics [START_REF] Bransden | Charge Exchange and Theory of Ion-Atom Collisions[END_REF] and is also used in cluster physics. Fortunately enough there exists a class of cases which can be addressed in a simplified, more robust way: distant collisions with highly charged ions, protons and, to some extent, also electrons can be considered as being structureless. What only counts is their Coulomb field. Charged ions are heavy and can be treated with classical trajectories R ext (t). For sufficiently heavy and fast projectiles, one can approximate these trajectories by straight lines and that is what is done in QDD. In any case, the effect of the ion, of charge Z ext , on the system can be described by a time-dependent external field [START_REF] Wang | [END_REF]:

V ext (r, t) = Z ext e 2 |r -R ext (t)| , (13a)
R ext (t) = R ext (0) + Ṙext (0)t . (13b)
Magnetic effects are neglected. They may play a role only for extremely fast ions in the relativistic domain. The ionic trajectories are characterized by the ion velocity v ion and the impact parameter b ion [98]. The velocity is more or less well defined by the experimental setup. But the impinging ion beam will cover a broad range of impact parameters. From the theoretical side, one has then to run several calculations with systematically varied impact parameters. Reaction cross-sections are then computed by integration of the reaction probability over impact parameters, for examples see [99,100]. Note that the mere Coulomb field (13) can only be used for distant collisions where the Coulomb singularity stay away from the system and that is the regime we intend with the charged projectile. Close collisions lead to significant electron transfer which requires careful modeling of a pseudo-potential for the colliding ion. More involved projectiles come into play when considering atomic collisions. Such scenarios could, in principle, also be treated in QDD and it had been investigated in the past in connection with collisions on surfaces [101]. However, the preparation of the initial condition is much handwork because it has too many degrees of freedom and has to be done in several steps. We thus do not discuss this option further in this public version of the QDD code.

Instantaneous boost

If the time of impact from a laser pulse or a charged projectile is shorter than the typical response time of occupied electron states, the pulse just instantaneously imprints a certain phase profile onto the s.p. wave functions which then unfolds in the subsequent electron dynamics. The leading mechanism here is a dipole boost which initializes dynamics as ϕ α (r, t = 0) = e ip0•r ϕ α,g.s. (r) (14) where p 0 is the boost momentum related to the boost energy

E boost = N el 2m e |p 0 | 2 . Such
an excitation mocks up realistic irradiation either through a very short laser pulse or a collision with a bypassing swift charged projectile. In both cases, the shortness of the excitation practically means that the initial boost delivers a perturbation covering all frequencies. As a consequence, the response of the system is dominated by the coupling to its eigenfrequencies and only depends on the amplitude of the boost. This considerably reduces the possibly large set of parameters in the case of a laser pulse or a charged projectile to explore various scenarios. Actually, very small boosts are typically used for evaluating the optical absorption spectrum, see section 3.2.4. The code also provides the option of an instantaneous displacement of the electron cloud which is a zero-time excitation as the boost, but with bias to higher energies.

Further instantaneous excitations

There are two more instantaneous excitations schemes built into the QDD code which trigger less prominent modes, namely spin-dipole and scissors oscillations. The spindipole mode is initiated by a spin-dipole boost (Eq. (M.8) in the Supplemental material) similar as the dipole boost in Eq. ( 14), however spin-up and spin-down boosted in opposite directions. The mode is tracked by spin-dipole moment, see section 3.2.5.

The scissors mode is an orbital J = 1 -excitation representing a rotational oscillation of one principle axis of the electron cloud against the corresponding ionic axis [102,103]. It is initialized by an instantaneous rotation (M.9) of the ionic background relative to the electron cloud. This mode is tracked by the angular momentum of the electron cloud, see section 3.2.5.

Equations of motion

The degrees of freedom of the theory are the s.p. wave functions ϕ α , their occupation number ω α , and ionic positions and momenta R I , P I . We derive in this section equations of motion for each of them. Those for the ϕ α 's as well as for the R I 's and P I 's are obtained from a variational principle on the total energy E total . For fixed ω α , the resulting equations of motion finally deliver coupled TDLDA-MD equations in which electrons are propagated according to TDLDA equations and ions according to classical Molecular Dynamics (MD). Such coupled dynamics is often quoted Ehrenfest dynamics [START_REF][END_REF]. The dynamical equations for the ω α 's requires a dedicated sequence of approximations well documented in the literature [104]. We shall not enter those details here and just recall briefly the path towards the RTA we use to effectively propagate the ω α 's.

Equations of motion for single particle electronic wave functions : TDLDA

Electronic Kohn-Sham equations and TDLDA

Static and dynamical equations determining structure and dynamics of the electron cloud are determined by variation of the total energy with respect to the s.p. wave functions ϕ α . This yields what is called the Kohn-Sham (KS) equations [64]. We use density functionals [65,66] i.e. functionals in Local Density Approximation (LDA). The dynamical application is called Time-Dependent LDA (TDLDA) and we use this acronym also in cases where we add SIC (see section 2.2.2) to the treatment.

The time-dependent KS equations determining the propagation of electronic wave functions are:

ĥσα ϕ α = i ∂ t ϕ α for α ∈ {1, ..., Ω} , (15a) ĥKS,σ
= p2 2m + V C (r, t) + Vxc [ρ σ ] + Vel,ion + V ext (r, t) VKS , (15b)
V C (r, t) = e 2 dr (r , t) |r -r | , (15c)
Vxc,σ = xc [ σ (r, t)] + σ (r, t) δ xc [ρ] δρ ρ= σ (r,t) , ( 15d
)
where ĥ is the KS mean-field Hamiltonian and σ ∈ {↑, ↓}. The exchange-correlation potential Vxc involves the construction δ xc [ρ]/δρ which is the functional derivative of the exchange-correlation energy per particle xc with respect to the local density ρ considered as an independent variable. This functional derivative is then evaluated at ρ = σ (r, t).

Vel,ion is the potential describing the interaction of the electrons with the ions, and is either a pseudopotential VPsP , see section 2.2.3, or stems from the jellium model, see section 2.2.4. The stationary KS equations analogously read

ĥσα ϕ α = ε α ϕ α for α ∈ {1, ..., Ω} , (15e)
where ĥ is composed in the same manner as in the time-dependent case (instantaneous approximation for V xc ). The KS equations (15e) pose a stationary eigenvalue problem. They provide the electronic ground state of a system. The time-dependent KS equations (15a) pose an initial value problem. The natural starting point is the ground state as obtained from the stationary KS equation. The numerical solution of the KS equations is explained in section 3. The time evolution delivered by the dynamical KS equations (15a) can be expressed formally by the unitary one-body time-evolution operator

Û (t, t ) = T exp - i t t ĥKS (t )dt ( 16a
)
where T is the time-ordering operator. This yields a closed expression for the timeevolution of s.p. states

|ϕ α (t) = Û (t, t )|ϕ α (t ) . (16b)
This compact form will be used later on in connection with the extension of the propagation by dissipation, see section 2.3.2.

LDA and TDLDA with mixed states

So far, we dealt with the equations for pure states. In the case of stationary states at finite temperature and in the case of dissipative dynamics, we encounter mixed states. This is described by associating an occupation number w α with each s.p. state ϕ α as provided already in the initial set-up in section 2.2. Even if the w α 's become crucial in the dissipative dynamics, see section 2.3.2, they are by definition kept frozen during a TDLDA propagation. They matter again for stationary states at finite temperature, in which case they are determined by:

w α = 1 1 + exp[(ε α -µ)/T ] (17)
where ε α is the s.p. energy of state α obtained from the stationary KS equations (15e), and µ is the chemical potential tuned such that the total electron number is reproduced by α w α = N el .

In case of a mixed state, it is advantageous to express the entity {ϕ α , w α } in compact manner by the one-body density operator:

ρ = ∞ α=1 |ϕ α w α ϕ α | Ω α=1 |ϕ α w α ϕ α | . ( 18
)
The solution of the thermal KS equations, i.e. Eqs. (15e) together with (17), provides the density operator immediately in this form which is called natural orbital representation (diagonal in the s.p. states). In general, the one-body density operator can be nondiagonal with respect to given s.p. states. This will play a role later on. As said above, the occupation numbers w α are kept frozen during TDLDA. The KS equations (15a) then can be written compactly as

i ∂ t ρ = ĥ, ρ (19)
where ĥ[ ↑ , ↓ ] is the KS Hamiltonian as above. This pure mean-field propagation (19) leaves the occupation numbers w α 's unchanged and propagates only the s.p. states. The mean-field propagation of an initial state (18) then reads:

ρ(t) = ∞ α=1 |ϕ α (t) w α ϕ α (t)| Ω α=1 |ϕ α w α ϕ α | = Û (t, 0)ρ(0) Û -1 (t, 0) ( 20
)
where Û is the mean-field evolution operator (16a).

Relaxation-Time Approximation (RTA)

The Relaxation-Time Approximation (RTA) is a robust and well manageable method to describe microscopically dissipation from two-particle collisions in a many-body system. Dissipation emerges generally from a reduction of the full dynamics of a large, iteracting system to a small subsystem. Here, in particular, we deal with the reduction of fully correlated many-body dynamics to a one-body dynamics within mean-field theory. A useful framework for that is the BBGKY hierarchy of density matrices which allows one to derive quantum kinetic equations [START_REF] Bonitz | Quantum Kinetic Theory[END_REF], the quantum generalization of the much celebrated Boltzmann equation [START_REF] Cercignani | The Boltzmann equation and its applications[END_REF]. These kinetic equations are very bulky and still extremely costly for finite fermion systems. The RTA, developed long ago for bulk fermionic matter [106,[START_REF] Ashcroft | Solid State Physics[END_REF], provides a robust and pertinent low-level approach to the (quantum) kinetic equations. We recently adapted them to finite fermion systems [62] which opens the door to affordable systematic investigations of dissipation processes in molecules [107].

In the following, we will present a brief summary of RTA and its pratical implementation and refer to [62] for a detailed derivation. To simplify the rather involved notations, we drop the spin index throughout this RTA part.

Formal background

The quantum Boltzmann equation is the quantum mechanical counterpart of the semiclassical Vlasov-Uehling-Uhlenbeck equation [4,3]. It complements the self-consistent TDLDA propagation of the one-body density operator ρ by dynamical correlations through a collision term. It reads in general [108,109] i ∂ t ρ -ĥ, ρ = Î[ρ] where the left-hand side contains the mean-field propagation. Î in the right-hand side stands for the quantum-mechanical collision term which, however, is extremely hard to handle for finite fermion systems. A great simplification can be achieved by the Relaxation-Time Approximation (RTA) which was used successfully in a wide variety of homogeneous systems [START_REF] Pines | The Theory of Quantum Liquids[END_REF][START_REF] Ashcroft | Solid State Physics[END_REF]. The RTA equations for the present case of finite fermion systems read [START_REF] Reinhard | [END_REF]:

∂ t ρ = - i ĥ, ρ - 1 τ relax (ρ -ρeq [ , j, E sp ]) , (21a)
(r, t) = α |ϕ α (r, t)| 2 w α , (21b)
j(r, t) = m e α w α ϕ * α (r, t) → ∇ - ← ∇ 2i ϕ α (r) , (21c)
τ relax = 0.40 σ ee r 2 s E * intr N el , r (el) s = 3¯ 4π -1/3 , σ ee = σ ee (¯ ) , ( 21d
)
where ρeq is the density operator of the thermal equilibrium for local density (r, t) current distribution j(r, t), and total energy E. The computation of E is often simplified by replacing it through a computation of the total s.p. energy E sp . This is justified for changes of energy at frozen local density (r, t). To see that, we decompose the energy as [START_REF] Maruhn | Simple models of many-fermions systems[END_REF]:

E = E sp + E rearr [ ] , E sp = α w α ε α = α w α ϕ α | ĥ|ϕ α . (22)
The rearrangement energy E rearr is a functional of local density and remains unchanged. Thus δE = δE sp , in other words, changes in total energy E at frozen are identical to changes in s.p. energy E sp . Thus we will use E sp in many of the subsequent equations whenever evaluating changes of E, but not always write it explicitely.

A crucial parameter is the relaxation time τ relax which is taken over from semi-classical Fermi liquid theory, for details see [START_REF] Reinhard | [END_REF]. Key entries are: the intrinsic (thermal) energy of the system E * intr (defined later on in Eq. ( 24), the actual number of valence electrons N el , the in-medium electron-electron cross-section σ ee , the effective Wigner-Seitz radius r (el) s of the electron cloud, and the average electron density ¯ . Note that r (el) s and σ ee depend on the average density ¯ because a spatially varying τ relax would be very cumbersome to implement in a quantum mechanical expression. The average density is deduced from the r.m.s. radius r rms of the actual electron cloud as ¯ = 3N el / 4π( 5/3 r rms ) 3 .

The RTA equation ( 21a) is rather involved because its entries depend in various ways on the actual state ρ(t). The most expensive piece is the instantaneous equilibrium density operator:

ρeq [ , j, E] = |ϕ (eq) α w (eq) α ϕ (eq) α | (23a)
which minimizes LDA energy with constraint on the actual (r, t), j(r, t) and energy E(t). It is determined by the Density-Constrained mean-Field (DCMF) equation

ĥDCMF [ , j]ϕ (eq) α = ε (eq) α ϕ (eq) α (23b) ĥDCMF [ , j] = ĥ[ , j, E] -drλ(r)ˆ (r) -dr λ j (r) • ĵ(r) -µ dr(ˆ (r) -(r, t)) 2 -µ j dr ( ĵ(r) -j(r, t)) 2 (23c)
where ĥ[ , j] is the KS mean-field Hamiltonian at the given instant. The ˆ (r) is the operator of local density at the space point r and ĵ(r) the operator of local current. The wanted energy E (or E sp (t) equivalently) and electron number N el are tuned by adjusting chemical potential µ (eq) and temperature T (eq) in the Fermi-Dirac distribution:

w (eq) α = 1 + exp ϕ (eq) α | ĥ|ϕ (eq) α -µ (eq) T (eq) -1 (23d) µ (eq) ↔ α w (eq) α = N el , T (eq) ↔ α w (eq) α ϕ (eq) α | ĥ|ϕ (eq) α = E sp (t) (23e)
Although cumbersome to evaluate, it is important to use exactly this local, instantaneous equilibrium in the relaxation term. This guarantees that the dissipative step conserves local density, current, and energy as it is mandatory for a correct collision term [START_REF] Gütter | [END_REF].

starting point:

ρ(t0) = α |ϕα(t0) wα(t0) ϕα(t0)| 1 ? mean-field propagation: |ϕ (mf) α (t1) = Û (t1, t0)|ϕα(t0) , w (mf) α = wα(t0) = const. ρmf = ρmf (t1) = α |ϕ (mf) α (t1) w (mf) α ϕ (mf) α (t1)| 2 ?
express ϕ and w through natural orbitals

ρmf (t1) = α |ϕ (nat) α w (nat) α ϕ (nat) α | X X X X X X X z 3 (r, t1) , j(r, t1), E (t1)
@ @ @ @ R ? Starting point: ϕ

(nat) α , w (nat) α
, εα (after RTA step 2) (r, t1) , j(r, t1) , E(t1) (after RTA step 3)

?

Fermi-Dirac distribution (23d): εα, N el , E -→ w (eq) α , µ (eq) , T (eq) ? Accelerated gradient step with constrained Hamiltonian: ϕ

(new) α ←→ O ϕα -D ĥDCMF -ϕα| ĥDCMF|ϕα ϕα ? ∆ 2 ĥDCMF ? < 1 ? yes Fermi-Dirac distribution (23d): εα = ϕα| ĥ|ϕα (new) , N el , E -→ w (eq) α , µ (eq) , T (eq)
no w (eq) α =const.

?

∆ 2ĥ DCMF ? < 0 & |E (new) sp -E (n) sp | ? < 2 ν |ρDCMF(rν ) -ρ(rν , t)|/ ν ? < ρ ?
yes compute E * intr , Eq. ( 24), and τ relax , Eq. (21d)

? back to RTA loop no 6 λ (n) + 2µ( (new) -) = λ (n+1) λ (n) j + 2µj(j (new) -j) = λ (n+1) j

Summary of the RTA procedure

We briefly summarize the solution scheme for one step from t ≡ t 0 to t 0 +∆t ≡ t 1 , for more details, see [START_REF] Reinhard | [END_REF]. The TDLDA propagation runs at a much faster pace than the relaxation. We resolve it by standard techniques [33,[START_REF] Reinhard | Introduction to Cluster Dynamics[END_REF] on a time step δt which is much smaller (by a factor 50-200) than the RTA step ∆t. We summarize this TDLDA propagation in the evolution operator Û from Eq. (16a) and only discuss the RTA step t 0 → t 1 :

1. We first propagate ρ by pure TDLDA. The s.p. states in diagonal representation (18) evolve as |ϕ α (t 0 ) → |ϕ

(mf) α (t 1 ) = Û (t 1 , t 0 )|ϕ α (t 0 )
, while the occupation numbers w α (t 1 ) = w α (t 0 ) are kept frozen (pure mean-field propagation). 2. Absorbing bounds (see section 3.1.6) may have removed parts from the s.p. wave functions and so destroy orthonormalization. We transform the propagated density operator to a representation in terms of natural orbitals which has the diagonal representation (18) with an orthonormal set of s.p. wave functions together with corresponding occupation numbers {ϕ

(nat) α , w (nat) α
}, for details on the procedure see appendix C.1. This step can be overridden if reflecting (or periodic) boundaries are used, in which case TDLDA preserves orthonormalization. 3. We compute density (r, t 1 ), current j(r, t 1 ), and total energy E, associated to the TDLDA-propagated density matrix ρmf . 4. We determine the thermal mean-field equilibrium state ρeq constrained to the given , j, and E from step 3. This is achieved by the DCMF equations ( 23) with an iterative algorithm sketched in figure 3. The equilibrium state ρeq is represented by new s.p. states {|ϕ (eq) α } and new occupation numbers w (eq) α in diagonal form (18). Having these, we determine finally the excitation energy as the energy relative to the zero-temperature state

E * intr = E sp - α w (T =0) α ϕ (eq) α | ĥ|ϕ (eq) α ( 24
)
where w

(T =0) α
are the occupation numbers from a Fermi distribution (23d) at temperature zero, i.e. a step function. The chemical potential µ (T=0) is tuned to match α w (T=0) α = N el using the given s.p. energies ϕ

(eq) α | ĥ|ϕ (eq) α
. This E * intr thus measures the amount of thermal excitation energy in the system, for details see [114]. 5. We compose the new density operator as mixture of TDLDA propagated state ρmf and equilibration driving term ρmfρeq with weight ∆t/τ relax as

ρmix = ρmf - ∆t τ relax [ρ mf -ρeq ]
where the relaxation time τ relax requires the actual intrinsic excitation energy E * intr , see Eq. ( 24) also obtained from DCMF. While evaluating the mixing, the new state is expressed in natural orbital representation (18). This yields the final s.p. states {|ϕ α (t 1 ) } for this step and preliminary new occupations wα . Technical details of this mixing step are found in appendix C.2.

6. The mixing in step 5 may have slightly changed the energy such that we remain with a small energy mismatch as compared to the aimed E. We now apply a small iterative thermalization step to readjust the energy, as outlined in Appendix A. This then yields the final occupation numbers w α (t 1 ) which comply with energy conservation.

The above steps are sketched in figure 2 whereby the step numbers here correspond to those encircled in the figure. The most involved part is step 4 for the determination of the solution of the DCMF equations. It is expanded in detail in figure 3. There are three termination criteria in DCMF iteration: convergence of the variance of s.p. energies 0 , energy convergence 2 , and convergence of average deviation of local density ρ . These are given as input parameters, for details see the Supplemental material. Note, that the actual thresholds for the values of these criteria may change with the system. As a rough guide, energy criteria scale with Fermi energy and density with average density of a system. An explicit example of convergence is discussed in section 3.1.7. As said above, the time step δt for propagation of TDLDA is very small because it is limited from above by the maximal energy on the grid representation. The stepping ∆t for the relaxation term needs only to resolve the changes in the actual mean field which allows much larger values. Typically, we perform 50-200 TDLDA steps before calling one RTA step. For detailed values, see the examples discussed in the Supplemental material, see section M.6.5.1 therein.

A word is in order about the system for which the present form of RTA can be used. The relaxation time τ relax is one global number chosen according to the average electron density ¯ , see Eq. (21d). This requires systems which can be characterized by such an average density, i.e., systems having only small density variations in the bulk as it holds typically for metallic bonds. The RTA rate is insensitive to many details of the microscopic collision term as energy-and angle-dependent scattering cross-sections [115] or a broad spectrum of relaxation rates. However, these details are usually resolved only (if at all) for fast and energetic processes which are anyway deep in the regime of semiclassical VUU. The grossly averaged treatment of RTA is acceptable for not too fast and not too energetic processes, preferably in compact systems.

Ionic dynamics

We finally add the ionic dynamics to complete the picture. Ions are described by classical Molecular Dynamics (MD), i.e. classical equations of motion, under the influence of their mutual Coulomb force, the forces experienced from the electrons, and possibly external forces. We start from the total energy (3) which serves as classical Hamiltonian in terms of R I and P I . Standard variation yields the classical equations of motion for the ions of positions R I and momenta P I

∂ t P I = -∇ R I E pot,ion (R I ) + N el α=1 ϕ α |V PsP (r -R I )|ϕ α , (25a)
∂ t R I = P I M I . ( 25b
)
They are to be propagated simultaneously with TDLDA for the electrons, here represented by wave functions ϕ α . The simultaneous propagation scheme is called TDLDA-MD. It applies to all dynamical situations including those that are far from the adiabatic limit and embraces truly diabatic scenarios. The practical realization adds the (simple) classical propagation to the evolution of the electronic states. For solution schemes, see section 3. All wave functions and fields are defined on a three-dimensional (3D) equidistant Cartesian grid of N x × N y × N z grid points. The spacing between the points is given as δx, δy, and δz in units of a 0 . For reasons of equilibrated accuracy, it is highly recommended to give the same value to all of them. Typical values depend on the atoms/ions involved. The pseudopotentials (see section 2.2.3) set the pace. The δx, δy, and δz must not be larger than √ 2 ln 2 times the smaller radius in the pseudopotential min(r loc , r l ) in Eqs. (7b) or (7d). In case of the jellium model (see section 2.2.4), the spacing should be of order of the surface width σ jel .

The grid is automatically arranged in such a way that in each direction the same number of grid points are located on both sides of the origin. The coordinate values for e. g. the x direction are thus:

1 - N x 2 δx, 2 - N x 2 δx, . . . , N x 2 -1)δx, N x 2 δx . (26)

Derivatives

The computation of currents and the action of the kinetic energy operator need first and second derivatives. We have two options for that in the code. Standard is the definition of derivatives via Fourier transform. This delivers high precision at affordable expense.

For simplicity, we explain here the Fourier strategy for one dimension. Given are N x discrete grid points x ν , ν ∈ {1, . . . , N x }, in coordinate space. They are mapped to the same number of grid points k n in Fourier space (physically equivalent to momentum space) as:

x ν = - N x 2 + ν δx , ν = 1, . . . , N x , (27a)
k n = (n -1)δk, n = 1, . . . , N x /2 , (27b)
k n = (n -N x -1) δk, n = N x /2 + 1, . . . , N x , δk = 2π N x δx . ( 27c
)
Note the particular indexing for the k-values. In principle, the values k n = (n -1)δk for all n are equivalent for the Fourier transform, but for the second half of this range the negative k-values should be chosen because of their smaller magnitude. For the Fourier expansion, k = -δk and k = (N x -1)δk are equivalent because of periodicity in k-space.

A function f (x ν ) in coordinate space is mapped to a function f (k n ) in Fourier space by:

f (k n ) = Nx ν=1 exp (-ik n x ν )f (x ν ) , (27d)
f (x ν ) = 1 N x Nx n=1 exp (ik n x ν ) f (k n ) (27e)
This complex Fourier representation implies that the function f is periodic:

f (x + N x δx) = f (x).
The appropriate integration scheme is the trapezoidal rule which complies with the above summations, adding up all terms with equal weight. The derivatives of the exponential basis functions are:

d m dx m exp (ik n x) = (ik n ) m exp (ik n x) . (28)
Computation of the m-th derivative thus becomes a trivial multiplication by (ik n ) m in Fourier space. Time critical derivatives are best evaluated in Fourier space using the Fast Fourier Transform (FFT). To that end, a forward transform (27d) is performed, then the values f (k n ) are multiplied by (ik n ) m as given in Eq. ( 28), and finally we transform backward by applying (27e) to (ik n ) m f (k n ) to come back to coordinate space. A word is in order about the first derivative. The upper point in the k-grid, δk N x /2, is ambiguous. Exploiting periodicity, it could be equally well -δk N x /2. Both choices introduce an unwanted bias. We circumvent the problem by setting k Nx/2 = -k Nx/2 = 0.

Handling of the pseudopotentials

As exposed in section 3.1.1.1, the maximal value of the grid spacing is related to the smallest length scale entering the pseudopotential. Local pseudopotentials for metals, see Eq. ( 6), are quite forgiving in the sense that their length scales σ i 's are relatively large. For instance, in the case of Na clusters, one can safely use δx = 0.8 a 0 . For Goedecker-like pseudopotentials, the values of r loc and r l 's are often smaller, particularly for elements with covalent binding which can require a large amount of grid points, thus rendering the calculations more expensive. The original values of r loc and r l 's for some chemical elements in the first, second and third rows [85,116] are implemented by default in the source code of QDD and can be used as such, with the drawbacks mentioned just above. To circumvent them, we have for some elements refitted pseudopotential parameters with a unique and possibly large length scale for r loc and r l . These are given with the input files of the examples of QDD application in the supplemental material.

Electronic ground state

The quality of ground state convergence is essential for the time propagation as a poor static convergence will generate spurious dynamical behaviors. A typical example is provided by the box size which needs to be chosen as large as possible in a balanced way between accuracy and computational cost. Maximum mesh size is basically given by the extension of the frozen core electrons as practically delivered by radial extension of a pseudopotential or surface thickness in the case of a jellium background. Box size itself is basically a free parameter once the box extension covers the ionic configuration and some electron skin around. However static convergence can be reached even with too small computational boxes and the defect may even not be clearly visible on static properties. Indeed the static convergence procedure involves reorthonormalization of s.p. wave functions at each static iteration which may blur the biasing effect of too small a box size. The ultimate test allowing to ensure that the static convergence is properly reached in a sufficiently large box is a dynamical one in which the computed ground state is propagated in time with absorbing boundary conditions and without any excitation (see section 3.1.6). Too small a box becomes then immediately visible as the statically constrained electrons become free to fly away and can be absorbed at the boundaries of the computational box.

The computation of the ground state amounts to diagonalize the KS Hamiltonian ĥ (15a) under the constraint of orthonormality of s.p. wave functions ϕ α , i.e. ϕ α |ϕ β = δ αβ . The static KS equations (15e), optionally combined with thermal occupation (17) of s.p. states, are solved iteratively. The wave functions are iterated with a gradient step (also known as steepest descent) which is accelerated by pre-conditioning with kineticenergy damping [117,118]. The static iteration from step n to step n + 1 can be written in the following compact form:

ϕ (n+1) α = O ϕ (n) α -D ĥ(n) KS -ϕ (n) α | ĥ(n) KS |ϕ (n) α ϕ (n) α (29a) D = δ T + E 0,damp (29b)
where T = p2 /(2m) is the operator of kinetic energy, O means orthonormalization of the whole set of new s.p. wave functions and ε

(n) α = ϕ (n) α | ĥ(n) KS |ϕ (n) α
is the expectation value of the KS Hamiltonian on ϕ (n) α at static iteration n. Because of self-consistency, the KS Hamiltonian is itself expressed as ĥ(n

) KS = ĥKS [ϕ (n) α ]
. This sort of kinetic-energy damping is particularly suited for the fast Fourier techniques that we use in QDD. The damped gradient step has two numerical parameters: the step size δ and the damping regulator E 0,damp . The latter should be chosen typically of the order of the depth of the local potential V KS (E 0,damp |ε 1 | where ε 1 is the lowest s.p. energy). Typical values are E 0,damp = 10, . . . , 50 eV, depending on the material under consideration. After proper choice of E 0,damp , the step size is of order of δ = 0.1, . . . , 0.5. Larger values yield faster iteration but can run more easily into pathological conditions.

In some cases, an initial finite electronic temperature, associated with Fermi distribution (17), can help the convergence of the statics, see e.g. section M.8.5. In that case, we have to realize that the wave function step has produced a new set of s.p. energies ε

(n+1) α = ϕ (n+1) α | ĥ(n+1) |ϕ (n+1) α
, which requires that we determine preliminarily new occupation numbers wα with Eq. (17). Using the wα immediately can occasionally lead to unstable iteration. It is safer to mix old and new occupation numbers gently as

w (n+1) α = η occ wα + (1 -η occ )w (n) α . ( 30
)
Typical values of mixing lie around η occ = 1/2. Critical cases may be stabilized with smaller values and forgiving cases can cope with larger values. Having the new s.p. wave functions and occupation weights, the new local electron densities (2a) are computed and then the new KS Hamiltonian (15). This provides the starting point for the next iteration. The process is continued until sufficient convergence is achieved. We consider as the convergence criterion the averaged energy variance of the s.p. states at step n:

∆ε (n) = α w α (∆ε 2 α ) (n) N el , (31a)
(∆ε 2 α ) (n) = ψ (n) α |( ĥ(n) KS ) 2 |ψ (n) α -(ε (n) α ) 2 , ( 31b
)
ε (n) α = ψ (n) α | ĥ(n) KS |ψ (n) α . ( 31c
)
Vanishing total variance ∆ε (n) signals that we have reached a stable energy minimum, i.e. a solution of the KS equations. However, this may be only a local minimum (isomeric state). It requires experience to judge whether one has found the absolute energy minimum. In case of doubt, one should redo a couple of static iterations from very different initial configurations. It turns out that ∆ε provides a stringent test of convergence, much superior to the mere decrease of total energy as often used. Alternative convergence criteria, also used in other codes, can be derived from checking the convergence of the local density distribution ρ(r). We explore here two options

δ (n) ρ = dr ρ (n+1) (r) -ρ (n) (r) dr , (32a) δ (n
) 2 ρ = dr ρ (n+1) (r) -ρ (n) (r) 2 dr . (32b)
This yields the density change per grid point. As reference density, one should take the average density ρ = 3N el 4π( 5/3 r rms ) 3 with r rms being the r.m.s. radius of the actual electron density. The dimensionless ratios δ (n) ρ/ρ or δ

(n)
2 ρ/ρ allow the comparison of the convergence between different systems. Figure 4 compares the three criteria using H 2 O as a test case. All three lines have exactly the same slope, although differing a bit in their offset. This means that all three criteria are equally well suited to measure the degree of convergence. The present code uses only the variance of s.p. energies as termination criterion (but prints the other criteria on demand).

The initial guess for the s.p. wave functions can be done in two ways. One option is to start with the wave functions of the deformed harmonic oscillator, for details see e.g. [START_REF] Maruhn | Simple models of many-fermions systems[END_REF]. These are characterized by n = (n x , n y , n z ), the number of nodes in each direction. We sort the states with increasing oscillator energy ε

(0) α = ω x n α
x + ω y n α y + ω z n α z and stop if the desired number of states is reached. The deformation of the initializing oscillator influences the initial state in two ways: first, through the deformation of the basis wave functions as such, and second, through the energy ordering of the ε (0) α and corresponding sequence of levels built. Variation of initial conditions means basically a variation of the oscillator radius and deformation. In particular, the initial deformation decides in which local minimum the KS iteration will terminate. This initialization by harmonic oscillators is well suited for metallic bonding where the wave functions spread over the whole system. Covalent bonding produces localized states and here it is more appropriate to start also from localized states. In that case, we place Gaussians at each ionic site and, if more is needed, higher harmonic oscillator wave functions. Bookkeeping is more involved then and will be explained in connection with the input to the code in the supplemental material. A word of caution is in order about open-shell systems, i.e. systems where the HOMO-LUMO gap is extremely small or zero. There can occur situations which have no DFT solution. The C atom is an emblematic example in that sense. It has four valence electrons. Two of them occupy the 1s state in the pseudo atom (frozen core electrons). There remain two spin up electrons (total spin has to be one) to be placed into the 1p shell with three available spin-up states. However, any choice of placing the two electrons leads to an unstable KS iteration. Thus one should check not only the numerics but also the physics if a system refuses to converge. A finite initial electronic temperature can help sometimes but not systematically. An instructive example is given in the Supplemental material, see section M.8.5.

(n) /ε F δ (n) 2 ρ δ (n) ρ

Ionic ground state

Three strategies are implemented for the optimization of the ionic configuration: steepest descent, dynamical cooling, and simulated annealing. The simplest method is steepest descent. For a given electronic configuration, one computes the Hellmann-Feynman forces on the ions and follows their direction for a short step. One re-iterates to the new electronic ground state and repeats the steps until convergence. This is the fastest method but it is prone to get stuck in local minima. One should use it only if one knows a reliable starting configuration.

Dynamical cooling allows one to better explore the energy landscape and thus to avoid distraction by unimportant side minima. From a given starting configuration, one runs full TDLDA-MD and keeps a protocol of ionic kinetic energy. Starting from a nonminimal energy configuration, the kinetic energy will first increase, which is accompanied by a corresponding decrease of potential energy. As soon as the kinetic energy turns to decrease, we stop propagation and reset all ionic velocities to zero. This defines the new starting point for the next round. The procedure is stopped if gain in kinetic energy falls below a given level of precision. This method is more forgiving than steepest descent. Still, it is also often kept in isomeric minima. One has to rerun it from different initial configurations to explore the landscape of minima.

The most elaborate, expensive, and reliable method is simulated annealing, for a detailed description see [119]. Simulated annealing explores in Monte-Carlo fashion the energy landscape with a supposed thermal distribution of configurations thereby reducing the temperature successively which leads at the end to a ground state minimum for temperature zero. The method has several parameters which need to be tuned carefully to a given situation. It requires some experience to use it efficiently [119], for details see section M.6.1 in the Supplemental material.

Electronic propagation

Electronic dynamics at the level of TDLDA is governed by the time-dependent KS equation (15a). QDD offers two different ways of determining the solution. Both methods starts directly from the formulation of the propagator Û (16a). The first method, called exponential evolution, simply employs a Taylor expansion of the propagator:

|ϕ α (t+δt) = T exp - i t+δt t dt ĥKS (t ) |ϕ α (t) ≈ m n=0 1 n! - iδt n ĥKS n |ϕ α (t) , (33)
where ĥKS in the expansion is taken at fixed time. At time t, only ĥKS (t) is known. However, using that in the approximate propagator (33) generates undue bias on t with disastrous consequences for energy conservation. The solution is to use a predictorcorrector scheme. For the predictor, we perform a half-time step, i.e. using Eq. (33) with ĥKS (t) and δt/2. This produces an intermediate set of wave functions ϕ α (t + δt/2) with subsequent densities and KS Hamiltonian ĥKS (t + δt/2). For the corrector step, we perform a full time step (33) within using the intermediate KS Hamiltonian ĥKS (t + δt/2) in the expansion. Appropriate values for the order of expansion are m = 4 for the predictor and 8 for the corrector. Below that, conservation laws are at stake, and above that, we will not gain much because the time step is also limited by the speed of change of ĥ(t). Altogether, this exponential evolution provides a reliable propagation with satisfying norm and energy conservation. An alternative propagation scheme is time-splitting method [120], often also coined split-step or split-operator technique and for the Schrödinger equation more precisely T Vsplitting technique because it relies on splitting propagation into kinetic energy T and potential V parts. We use in the following the term "time-splitting method". In addition, we extend it here for use in connection with non-local pseudopotentials presented in section 2.2.3. The KS Hamiltonian is split into three pieces: ĥ = Tkin + Vnloc + V KS,loc (r, t)

where Vnloc is the piece stemming from the non-local part of the pseudopotentials while V KS,loc (r, t) collects all parts which form together a local potential. The latter is strongly time-dependent due to the self-consistent electronic contributions. The non-local part depends only on ionic motion which is snail-slow at electronic scale such that we can ignore time dependence during an electronic step. For going from t 0 to t 1 = t 0 + δt, we define a propagator for each one of the three parts and factorize the full propagator as:

|ϕ α (t 1 ) = Û (t 1 , t 0 )|ϕ α (t 0 ) ≈ exp - i δt 2 V KS,loc (r, t 1 ) Ûnloc exp - i δt Tkin Ûnloc × × exp - i δt 2 V KS,loc (r, t 0 ) |ϕ α (t 0 ) , (34b) Ûnloc
= m n=0 1 n! - iδt 2 n Vnloc n . (34c)
The sequence (34b) is built in symmetric manner to minimize the separation error. The quality of the time-splitting method depends on the size of the commutators between the three pieces of ĥ. The smaller the commutators, the better it is. The advantage of the method is that the two crucial propagators, kinetic and potential energy, are performed exactly. The local potential energy operator can be easily exponentiated in coordinate space, while the kinetic energy operator is exponentiated in Fourier space for which we employ forward and backward FFT as for the evaluation of derivatives (see section 3.1.1). The propagator of the non-local potential could also be resolved exactly, however with considerable bookkeeping expense. It is thus handled by the Taylor expansion (34c). The time-splitting has a particular advantage concerning evaluations of the KS potential. Note that the separation (34b) employs the potential at t 0 in the first potential step and at t 1 in the second. This avoids any bias by giving the same weight to both times. The key saving now comes with the fact that the local propagator exp -i δt 2 V KS,loc (r, t 1 ) only changes the complex phase of the wave function and does not change the local density. Thus we can evaluate the final new density for t 1 already at the stage before applying this local propagator and no predictor-corrector strategy is needed any longer. We exemplify and compare the two propagation schemes for the case of the Na + 9 cluster with explicit ionic structure described by local pseudopotentials (see section 2.2.3). We aim at testing the stability of the steppings and do that by propagating the electronic ground state of Na + 9 with frozen ions. Let us start with the exponential evolution (33). Besides time step δt, it has the order of Taylor expansion m which appears twice: once in the predictor and once in the final step. A fast scheme is obtained with m = 2 for the predictor and m = 4 for the corrector. Even if it can work for a while, it develops instabilities rather early (in the present test case, after a few 100 fs). A sensitive check is the conservation of the norm of the s.p. wave functions because the exponential step at finite m is not exactly unitary, whereas an exact propagation is. Insufficiencies then show up first in a slight drift of the norm that we observe for the cheap choice m = 4 rather early, long before true instabilities pop up. A much better norm conservation is obtained with m = 4 for the predictor and m = 8 for the step. This is this choice that is shown in figure 5 for two different time steps. The case δt = 3.53 as (1 as = 10 -18 s) is rock stable for very long, the norm is conserved up to 12 digits, and energy as well as dipole stay solidly at with explicit ionic background. Lower panels: time evolution of dipole moments. Upper panels: time evolution of the total energy. Right panels: View over a long time interval up to 2.4 ps. Left panels: zoom on a shorter time interval usually sufficient for simulating electronic effects. The acronym "T V " stands for the time-splitting method (34b) and "exp" for exponential evolution (33) of order m = 8. The numbers in the legend indicate the time step in units of attoseconds.

the ground state value. The slightly larger δt = 3.54 as, on the other hand, explodes after 500 fs. Such a sudden transition is typical for the exponential propagation and has to do with the convergence radius of the expansion. The limit of stability is related to the maximal energy in the numerical representation. It may change in truly dynamical situations such that one must choose a smaller step size. With δt = 2.4 as, we are on the safe side for all dynamical scenarios in Na + 9 . Larger values are possible if carefully checked.

The situation is much different for the time-splitting (34b). This propagation is manifestly unitary (in that case without non local pseudo potential) such that the norm of the s.p. wave functions is strictly conserved. But there is a slight error in dipole amplitude and energy from the onset (better visible in the left panels) which is due to the fact that we do not use the exact propagator U but an approximate factorized expression of it. The splitting error depends on the time step δt as can be seen in the figure (mind the scales, though). Once we have accepted the initial error, the propagation time-splitting method proceeds in stable manner for long. The larger δt = 4.8 as produces gradually increasing deviations after 1 ps (see right panels). They are still rather small and may be acceptable when simulating hefty excitations with dipole amplitudes and energies larger than the deviations here and with the smaller δt = 2.4 as one remains at the safe side for much longer. It is also to be noted that the time-splitting scheme allows a time step of δt = 4.8 as, larger than the maximum value allowed for the exponential evolution (3.53 as). This indicates once more that the time-splitting is more forgiving. Most importantly, the time-splitting scheme is an order of magnitude faster than exponential evolution. This therefore makes it by far the preferred option, even at the price of some small initial deviations.

In any case, this example nicely demonstrates the commonly known fact that any numerical propagation scheme has a limited time span of reliability which depends on the numerical representation and, in case of self-consistent mean fields, on the actual dynamical scenario. It is, therefore, strongly recommended to launch a few test runs for any new application before starting large scale production.

Absorbing boundary conditions

A grid representation naturally leads to reflecting or periodic boundary conditions. Reflection emerges for finite difference schemes. A representation of the kinetic energy by complex Fourier transform is associated with periodic boundary conditions where flow leaving the box at one side is re-fed at the opposite side. Both can lead to artifacts if a sizable fraction of electronic flow hits the boundaries. Proper handling of electron emission requires absorbing boundary conditions which hinder outgoing electrons from coming back into the simulation box. There are several ways to solve the problem, for a detailed discussion see [121]. The QDD code employs the simple strategy of applying a mask function during time evolution [122]. This technique is particularly easy to implement and has been widely used in the past. Its robustness and efficiency allow one to develop advanced analyzing techniques on the grid as, e.g., the computation of Photo-Electron Spectra (PES) and Photo-Angular Distributions (PAD) [36,[START_REF] Reinhard | Timedependent density functional theory[END_REF]. We give below a short explanation.

The left panel of figure 6 sketches the implementation of absorbing boundary conditions with computation of PES and PAD on a coordinate space grid (as 2D cut of a 3D grid). The absorbing boundary conditions are indicated by the ring area in the figure in between the radii R in and R out , covering here 3 grid points in each direction (actual calculations typically use 6 and more points). The absorption is performed in each time step as:

ϕ α (r, t) -→ φα (r) = Û (t+δt, t) ϕ α (r, t) , (35a)
ϕ α (r, t+δt) = M(r) φα (r) , (35b)
M(r) =        1 for |r| < R in , cos γ M π 2 |r| -R in R out -R in for R in < |r| < R out , 0 for |r| > R out . (35c)
First comes one standard KS step (see section 3.1.5) expressed here in terms of the TDLDA propagator Û defined in (16a), which yields the intermediate wave function φα (r). This is followed by the action (35b) of the mask function M defined in (35c) and shown in the right panel of figure 6, which steadily reduces the norm of the wave functions from the inner mask radius R in to the outer one R out . We use here a spherically symmetric mask. The spherical profile is helpful to minimize griding artifacts when computing angular distributions [124] (simpler rectangular masks may be used if PAD and PES are not of interest).

R in < l a t e x i t s h a 1 _ b a s e 6 4 = " q u D n t 3 A q X H D M 6 n f o 8 R U V a A H 1 D u c = " > A A A C 0 H i c j V H L S s N A F D 3 G V 3 1 X X b o J F s F V S V T Q p e j G p Y p t B V v K Z J z W w b y c T E Q p R d z 6 A 2 7 1 q 8 Q / 0 L / w z p i C W k Q n J D l z 7 j 1 n 5 t 4 b p K H M t O e 9 j j i j Y + M T k 6 W p 6 Z n Z u f m F 8 u J S P U t y x U W N J 2 G i T g O W i V D G o q a l D s V p q g S L g l A 0 g s t 9 E 2 9 c C 5 X J J D 7 R t 6 l o R a w b y 4 7 k T B P V O m 4 3 I 6 Y v V N S T c b 9 d r n h V z y 5 3 G P g F q K B Y h 0 n 5 B U 2 c I w F H j g g C M T T h E A w Z P W f w 4 S E l r o U e c Y q Q t H G B P q Z J m 1 O W o A x G 7 C V 9 u 7 Q 7 K 9 i Y 9 s Y z s 2 p O p 4 T 0 K l K 6 W C N N Q n m K s D n N t f H c O h v 2 N + + e 9 T R 3 u 6 V / U H h F x G p c E P u X b p D 5 X 5 2 p R a O D H V u D p J p S y 5 j q e O G S 2 6 6 Y m 7 t f q t L k k B J n 8 D n F F W F u l Y M + u 1 a T 2 d p N b 5 m N v 9 l M w 5 o 9 L 3 J z v J t b 0 o D 9 n + M c B v W N q r 9 Z 3 T j a q u z u F a M u Y Q W r W K d 5 b m M X B z h E j b y v 8 I g n P D v H z o 1 z 5 9 x / p j o j h W Y Z 3 5 b z 8 A E E P Z T m < / l a t e x i t > R out < l a t e x i t s h a 1 _ b a s e 6 4 = " y P Y C U G m + L 7 u 2 O 7 w B I u 1 Q A S H 4 Y o A = " > A A A C 0 X i c j V H L S s N A F D 2 N r 1 p f V Z d u g k V w V d I q 6 L L o x m V 9 9 A F t L c l 0 2 g b z Y j I R S i m I W 3 / A r f 6 U + A f 6 F 9 4 Z U 1 C L 6 I Q k Z 8 6 9 5 8 z c e 5 3 I c 2 N p W a 8 Z Y 2 5 + Y X E p u 5 x b W V 1 b 3 8 h v b t X j M B G M 1 1 j o h a L p 2 D H 3 3 I D X p C s 9 3 o w E t 3 3 H 4 w 3 n 5 l T F G 7 d c x G 4 Y X M l R x D u + P Q j c v s t s S d T 1 R b f t 2 3 I o / H G Y y E k 3 X 7 C K l l 7 m L C i l o I B 0 V c P 8 C 9 r o I Q R D A h 8 c A S R h D z Z i e l o o w U J E X A d j 4 g Q h V 8 c 5 J s i R N q E s T h k 2 s T f 0 H d C u l b I B 7 Z V n r N W M T v H o F a Q 0 s U e a k P I E Y X W a q e O J d
The limiting radii need yet to be specified. The outer radius is the smallest distance from the origin of the box to the bounds, formally R out = min rν ∈bounds (|r ν )|) where the index ν indicates that r ν is a grid point. The inner radius is then

R in = R out -N M δx , ( 36
)
where δx is the grid spacing. The N M characterizes the number of absorbing points in one direction. It is a crucial numerical parameters, see table M.10 in the manual. This looks simple and straightforward. However, the mask technique is not perfect. One will always encounter a small amount of reflected flow, particularly for electrons with low kinetic energy. One can minimize the back-flow by proper choice of the exponent γ M entering the mask profile, see Eq. (35c). This depends, however, on the actual numerics (number of absorbing points, size of time step). Typical values of γ M are of order 1/8 or lower. A detailed description and discussion of this approach and its proper choice of numerical parameters is found in [121].

Absorbing boundary conditions introduce a subtle difficulty in the time propagation. Indeed, while electronic propagation is built such that it preserves orthonormality of s.p. wave functions (see section 3.1.5), orthonormality of s.p. wave functions is gradually lost as soon as absorbing boundary conditions are active. There are two ways to deal with that. The first is to assume that the wave functions in full space remain orthonormal and simply to carry on. Ionization is deduced in this case from the loss of norms of the s.p. wave functions. The second is to diagonalize the one-body density matrix in the simulation box which delivers the natural orbitals, for a detailed discussion see [125]. In the natural orbital basis, the s.p. basis remains orthonormal by construction. The change of electron content caused by the absorbing boundary conditions then leads to a change of occupation numbers in the natural basis. Both schemes to TDLDA yield exactly the same results (up to numerical accuracy), see e.g. [125], but count ionization in different manner, from loss of norms for straightforward TDLDA propagation, or from changing occupation numbers in the basis of natural orbitals. RTA does not leave a choice. It uses natural orbitals by construction and it produces a genuine change of occupation numbers from relaxation. We conclude from this discussion that the change of occupation numbers in RTA with absorbing bounds has two sources: relaxation processes modeled by RTA and ionization induced by absorbing boundary conditions.

DCMF

The solution of the DCMF equations ( 23) was already sketched in figure 3. We add here a few more explanations. The basic ingredient is the solution of the mean-field equation (23b) for which we use again the accelerated gradient step (29) with its own damping parameters however, with different damping parameters E 0,damp → E (RTA) 0,damp and δ → δ (RTA) . New is the iteration of the Lagrangian parameters λ(r) and λ j (r) for density and current constraint in the DCMF mean-field Hamiltonain (23c). These are driven by the mismatch of density and current as

λ (n+1) = λ (n) + 2µ( (new) -) λ (n+1) j = λ (n) j + 2µ j (j (new) -j) (37)
via the quadratic term whose coefficients µ and µ j are numerical parameters kept constant throughout the calculations. They should be large enough to drive efficiently. However too large values slow down iterations again. A proper compromise for them has to be found and it depends on the actual system. For more details and an example, see the supplemental material, section M.6.5. The constraint iteration adds one more termination criterion. As in the accelerated gradient step, we check for convergence of energy in terms of energy change

∆E sp = |E (n+1) sp -E (n) sp | and of variance of s.p. energies √ ∆ 2 h DCMF = (1/N e )( α w (n) α ϕ (n) α |( ĥDCMF ) 2 |ϕ (n) α -w (n) α ϕ (n) α | ĥDCMF |ϕ (n) α 2
), see also figure 3 in section 2.3.2.2. New is a check for convergence of density in terms of deviation from the goal

δρ DCMF = N el dr|ρ DCMF (r) -ρ(r, t)| dr , ( 38a
)
δj DCMF = N el dr|j DCMF (r) -j(r, t)| dr . ( 38b
)
In practice, we found that it suffices to check convergence of δρ DCMF .

The DCMF scheme also requires to tune the Fermi distribution to the given particle number and energy, see Eq. (23e). This is solved by interlaced bracketing [START_REF] Press | Numerical Recipes[END_REF] which is reliable and robust, although not very efficient. But this is not an issue for this inexpensive part of the RTA scheme. 34

The most critical and time consuming part in the RTA step is the solution of the DCMF equations ( 23) to obtain the local instantaneous equilibrium state. Figure 7 shows an example of typical convergence pattern for DCMF during RTA evolution. The left lower panel shows the dipole response to the external laser field to demonstrate the dynamical regime in which the example works. The difference between RTA and TDLDA indicates the times at which dissipation is most active, namely in the early stages where the two curves start to deviate and where the impact of the laser amplitude (e.g. in terms of energy input) is large. The left upper panel shows the number of DCMF iterations which are needed until the required precision demands are fulfilled. This number changes dramatically along the dynamical path and it is, indeed, largest where dissipation has most work. It is to be noted that the needed number of iterations sensitively depends on the system and strength of excitation. In the present case, we could afford to set a limit at 1000 iterations. There are situations where one needs considerably more. The right column of the figure shows detailed convergence patterns in the top panel, variance of s.p. energies √ ∆ 2 h DCMF and, in the bottom panel, deviation of density δρ DCMF /ρ, see Eq. (38a), at three different times. The convergence criteria are not simply going downward, particularly in the slowest case. The interplay of energy optimization and density matching in DCMF can lead to significant fluctuations in the iteration process. After all, with diligent tuning of the convergence parameters and a bit of patience, we obtain the nice result in the left lower panel which demonstrates that proper dissipation eliminates the unphysical long-term oscillations of a resonantly excited system.

Ionic propagation

Time propagation for ionic molecular dynamics (MD) is done with the robust and reliable Verlet algorithm [START_REF] Verlet | [END_REF][START_REF] Press | Numerical Recipes[END_REF], also known as leap-frog algorithm. We present it here for one ionic core described by the two vectors R and P, and thus drop the subscript I that labels the ionic core (we correspondingly label the ionic mass M ). We denote by δt I the ionic time step which is a multiple of the electronic time step δt.

Propagating the ionic equations of motion [START_REF] Wilken | [END_REF] raises the same problem with update of the forces as we had experienced in quantum mechanical propagation. The Verlet algorithm solves this by taking coordinate R and momentum P at different time grids, shifted by half a time step δt I /2. More precisely, this proceeds as:

R(t+δt I ) = R(t) + δt I P t+ 1 2 δt I /M , (39a)
P t+ 3 2 δt I = P t+ 1 2 δt I + δt I F R(t+δt I ), t+δt I , (39b)
where F denotes the total force acting on the ionic core. The advantage of this algorithm is that this force is needed only for propagating P and it uses the coordinate R evaluated on the shifted time grid, thus readily available at the moment when it is needed, because it was computed using the previous momentum P.

The disadvantage of the method appears at initialization and for computation of observables. In both cases, one has to move one of the variables R or P by half an ionic time step to have both at the same time. This problem is avoided by slight modification of the Verlet algorithm, called velocity Verlet [START_REF] Press | Numerical Recipes[END_REF]. Here, both variables are handled at the same time grid. The momentum P(t+ 1 2 δt I ) at intermediate time is estimated by a predictor and the force at intermediate time is taken as an average of the forces at t and t+δt I . This amounts to:

R(t+δt I ) = R(t) + δt I P(t+δt I )/M + 1 2 (δt I ) 2 F (R(t), t) /M (40a)
P(t+δt I ) = P(t) + δt I 1 2 F(R(t), t) + F(R(t+δt I ), t+δt I ) . (40b)
Note that the force at the new time F(R(t+δt I ), t+δt I ) is available for the P propagation (40b) because the R propagation (40a) has been completed before. Typical values of δt I are a couple of δt, the latter being the electronic time step. There is no general rule for how large δt I can be. This strongly depends on the system under study. For explicit examples, the reader can refer to the supplemental material. In any case, both schemes, Verlet and velocity Verlet, are equivalent in terms of precision, energy conservation, and stability.

A typical example of ionic motion is given in figure 8. Test case is the H 2 O molecule where the electrons are excited initially by an instantaneous boost (see section 2.2.5.3). The electronic response (see left column) is very fast. The dominant part of the excitation energy is used up for ionization, leaving a charged molecule out of Coulomb equilibrium and still with a considerable amount of electronic excitation energy as visible from the 36 Propagation was done with δt = 0.6 as and δt I = 10δt = 6 as. Plotted as functions of time are: left column, ionization (top) and electronic dipole moment in each spatial direction (bottom); middle panel, ionic kinetic energy E kin , total and its three spatial components; right column, ionic velocities (bottom) and positions (top) in x direction.

ongoing remaining dipole oscillations. Coulomb pressure and dipole fluctuations drive, at a slower pace, ionic motion. In spite of the considerable excitation, the ionic amplitudes remain small, reflecting a small transferred energy, at least in the short initial time interval shown here (full thermal equilibration takes much longer). Indeed, the maximum total ionic kinetic energy is less than 1 eV, while the initial boost brings almost 30 eV in the system. Along with this small transferred energy, we plot in the right column of the figure the ionic velocities and positions in x direction. We observe small amplitude oscillations, rather regular and at a pace given by the leading ionic frequencies, much smaller than the electronic ones, visible in the bottom left panel which shows the electronic dipole response with very fast oscillations. Through this example, we demonstrated the non-trivial cross-talk between electronic and ionic dynamics, made possible by the non-adiabatic coupling between electrons and ions in our TDLDA-MD at the level of the Ehrenfest approximation.

Observables

Energies

We denote by E(N el , R (Nion) ) the total binding energy of a system consisting of N el TDLDA electrons and N ion ionic cores. Energy is a most prominent observable and it naturally results from any calculation with density functionals. Comparison with measurements is usually done in terms of differences of energies, e.g., the monomer separation energy as the adiabatic energy difference

E mono = E(N el , R (Nion) ) -E(N el - 1, R (Nion-1)
) where both energies are to be taken from fully relaxed ionic configurations, or the vertical ionization potential (IP)

E IP = E(N el , R (Nion) )-E(N el -1, R (Nion)
) where the new electronic state in the N el -1 system has relaxed but the ions are kept in their original configuration.

As a byproduct of mean-field calculations, one also obtains the series of s.p. energies ε α . But it is known that the ε α 's from (TD)LDA are spoiled by the self-interaction error. This defect can be cured by a SIC, see section 2.2.2, after which the set of ε α provides a fair map of electron separation energies, particularly of the IP [63,74]. Experimental data on s.p. energies are mainly the IP and the sequence of peaks in a PES from onephoton processes with laser pulses of weak intensity. Note that the QDD code also offers the possibility to compute fully dynamically a PES in case of pure TDLDA, see section 3.2.6.4. Besides comparison with data, the s.p. energies are very instructive observables for theoretical interpretation as, e.g., analyzing electronic shell structure [128,89].

Other "theoretical" observables are the separate energy contributions which are helpful to analyze the energy balance in dissipative dynamics [107]. We will address later on in this paper the intrinsic excitation energy E * intr and the energy absorbed from the laser field E abs . The E * intr characterizes the amount of internal, so to say thermal, kinetic excitation of the electron cloud. This is a quantity which plays a key role in the model for dissipation, as was explained in section 2.3.2.2. The absorbed energy is computed as

E abs (t) = t 0 dt dr E(t ) • j(r, t ) (41)
with an additional small correction for particle loss at the absorbing bounds (for details see [107]). It provides the total excitation energy delivered to the system.

Densities and shapes

Density functionals also produce the electronic local density (r, t) as a natural outcome. The electronic density together with the ionic configuration {R I } constitutes the shape of a system in all detail. But this can become an overwhelming amount of information. The gross structure of a shape is more simply sorted in terms of its multipole moments. Leading quantity is the root-mean-square (r.m.s.) radius which reads for the electrons:

r rms (t) = dr (r + r cm , t) r 2 N el (42a)
where r cm is the electronic center of mass. Next come the quadrupole moments which quantify the deformation. The QDD code computes and prints the Cartesian quadrupole moments

Q ij = dr (r + r cm , t) r i r j - δ ij 3 r 2 . ( 42b
)
This is provided not only for the electronic density but also for the ionic configuration by replacing dr ... through I .... The quadrupole moments (42b) as such change scale with system size. This inhibits simple comparison of shapes between systems of different size. To that end, it is advisable to work with dimensionless quadrupole moments which can easily be deduced from the above radius and quadrupole moments. This is outlined in appendix B.

Polarizability

The static polarizability is a key observable of atoms, molecules, and clusters. We discuss here the most important case of dipole polarizability α D . One applies a static external dipole field V ext (r) = eE (stat) 0

•r and performs static calculations for a couple of E (stat) 0

. This delivers an electric dipole moment e D = e r = e dr r (r, t) as a function of E (stat) 0

. The tensor of static dipole polarizability is then

(α D ) ij = e ∂ Di ∂E (stat) 0,j E (stat) 0 =0 for i, j ∈ {x, y, z} . (43)
The polarizability matrix simplifies if the system has spatial symmetries. For example, an axially symmetric system aligned with the z-axis has (α

D ) xz = (α D ) yz = (α D ) xy = 0 and (α D ) xx = (α D ) yy . A spherically symmetric system has additionally (α D ) xx = (α D ) zz .

Optical response

Optical response is a key observable in cluster physics. In the limit of long wavelengths, the laser field at the system site is related to the dipole moment D = r where r is taken with respect to the center-of-mass of the electron cloud. Within QDD, the dipole excitation strength can be computed by spectral analysis of the dipole response to an instantaneous dipole excitation of the system. To that end, one starts from a well relaxed ground state and applies the instantaneous initial excitation by a small dipole boost ϕ α (r, t = 0 + ) ← exp(ip 0 • r)ϕ α (r, t = 0), see section 2.2.5.3. One then propagates electrons with TDLDA and samples a protocol of the dipole moment

D(t) = dr r (r, t) . (44a)
After a sufficient time, say T max , one multiplies the dipole signal D(t) with an appropriate window function W(t) [START_REF] Press | Numerical Recipes[END_REF], Fourier transforms it into D(ω), and finally obtains the spectral strength S Di (ω) for i ∈ {x, y, z} and corresponding spectral power P Di as:

S Di (ω) = { D i (ω)} p 0,i , D i (ω) = dt W(t) e iωt D i (t) , (44b)
P Di (ω) = | D i (ω)| 2 p 2 0,i , (44c)
where p 0 is the boost momentum initially applied. We here ignore the posssibility of (small) cross-talk between different channels i and j and concentrate on the diagonal elements of the strength tensor. The maximum possible spectral resolution is given by δω = 2π/T max . The window function W(t) serves to attenuate the dipole signal toward the end point T max and to avoid artifacts from non-zero D(T max ) [START_REF] Press | Numerical Recipes[END_REF]. Useful windows are W(t) = cos 2n tπ/(2T max ) where n is an integer number. The choice n = 1 produces often still a bit rough { Di (ω)} while n = 2 usually performs satisfyingly well. It is interesting to note that this treatment in connection with absorbing boundary conditions (see section 3.1.6) allows one to compute correctly the escape width of spectral states lying in the electron continuum. For details of spectral analysis and variants thereof, see [129]. Spectral analysis can equally well be performed with other observables, as e.g. higher multipoles, spin modes, etc. The full TDLDA furthermore allows one to go beyond the linear regime. The more appropriate observable is then the power spectrum (44c).

The same principles of spectral analysis can also be applied for computing the vibration spectra of clusters and molecules. Here we concentrate on the MD part of TDLDA-MD. The analyzing times are then to be taken much longer to supply sufficient spectral resolution for excitation in the meV range, characteristic of ionic vibrational states, and the ionic multipole moments ought to be used for the spectral analysis of ionic motion [130].

Momenta probing exotic modes

As explained in section 2.2.6, QDD also allows one to excite more subtle modes as spin-dipole and scissors oscillations. These modes require, of course, specific observables for proper analysis. The spin-dipole mode is tracked by the spin-dipole moment D S (t) = dr r ↑ (r, t)dr r ↓ (r, t) .

(

) 45
The observable related to the scissors mode is the orbital angular momentum

L = α w α ϕ α |r × p|ϕ α ( 46
)
where r is to be taken with respect to the center-of-mass of the electron cloud.

Ionization

Absorbing boundary conditions as explained in section 3.1.6 provide a pertinent picture of electron emission. There are several observables associated with emission from (photo-)excited systems: total ionization, ionization probabilities, photo-electron angular distribution (PAD), and photo-electron spectra (PES). We discuss them in the next four paragraphs.

Total ionization

The first observable is the total ionization, i.e. the number of escaped electrons denoted by N esc . The, now decreasing, norm of the s.p. wave function ϕ α defines the depletion N α of state α, from which one computes the ionization from this state, N esc,α , as the complement of N α :

N α (t) = ϕ α (t)|ϕ α (t) and N esc,α (t) = 1 -N α (t) . (47)
The total ionization then reads as the sum of the level ionizations:

N esc (t) = Ω α=1 ω α N esc,α (t) (48)
Both, total ionization and detailed level depletion as functions of time, are very instructive observables [131].

Ionization probabilities

Total ionization can be further analyzed in terms of detailed ionization probabilities. The probability to find the system at time t in one of the possible (positive) charge states k to which the system can ionize is denoted by P (+k) (t). A scheme to estimate the P (+k) (t) has been derived for pure mean-field states with integer occupation probabilities by combinatorial analysis [132]. We summarize it here briefly. The physical change of charge state is related to ionization which is realized numerically by absorbing boundary conditions. From the time evolution of the level depletions and ionizations per state, N α (t) and N esc,α (t) respectively, defined in Eq. ( 47), one starts with the following relation

1 = Ω k=0 P (+k) (t) = Ω α=1 N α (t) + N esc,α (t)) . ( 49
)
The product in the right-hand side is then expanded and the resulting terms are rearranged, collecting terms containing k factors N k α (N esc,α ) Ω-k which allows one to identify them to the ionization probability P (+k) [132]. For examples of application and details see, e.g., [99,33].

The irradiation of the C 3 chain by an intense laser pulse is taken as an example and various time evolutions are shown in figure 9. The dipole moment (lower left panel) gives an impression of the dynamical response of the molecule. It shows the typical pattern for an off-resonant excitation. The dipole signal follows the laser profile and some small asymptotic oscillations remain, due to a small excitation of the system's eigenmodes. The total ionization (top red curve in upper left panel) increases in the region of large dipole amplitudes and then levels off. We also show the ionizations from the HOMO, HOMO-1 and two degenerate HOMO-2, 3 levels. These are the key ingredients in evaluating the charge-state probabilities P (+k) . The emerging P (+k) are shown in the right panel. Their time evolution corresponds to the pattern of total ionization. The neutral state loses probability with increasing ionization while all other take increasingly their share. Finally, the average charge state of Q = 0.8, read off from the total ionization, is mainly distributed with 40% of charge state Q = 0 and 40% of charge state Q = 1. Charge state +2 gather appreciable 14% and the higher charges share the leftovers with rapidly decreasing amounts (only visible on logarithmic scale). The distribution visualizes a subtle feature of TDLDA in connection with absorbing boundary conditions: it generates inevitably a mixed state (which becomes apparent in a naturalorbital representation [133]) and this state embraces very different charge states when decomposed properly. The different charges produce different Coulomb fields. TDLDA can describe only the average charge and an average mean field.

So far, we have discussed the case of pure mean-field states with integer occupation. RTA inevitably produces mixed states with fractional occupations w α . The extension to this situation is formally trivial and actually implemented in QDD, provided that we merely generalize the level depletion of state α as N α (t) = w α ϕ α (t)|ϕ α (t) and apply the combinatorial decomposition as before. Although technically simple, the physical interpretation of the resulting P (+k) is involved. Fractional occupations w α as generated by RTA produces also a spread in charge states which has nothing to do with ionization, but reflects a general weakness of a thermal mean-field state, namely that it represents unavoidably an ensemble covering different electron numbers [START_REF] Reichl | A Modern Course in Statistical Physics[END_REF]. The resulting P (+k) of RTA with absorbing boundary conditions thus mixes effects from ionization with effects from thermalization. It is a task for future research to disentangle these two contributions.

Photo-electron angular distributions

Photo-angular distributions (PAD) are evaluated in angular segments, see Fig. 6. An angular segment can be described by an azimuthal angle θ and a polar angle φ. The reference frame for these two angles is usually the z axis identical with the laser polarization axis. We collect all probabilities which were removed by the absorption step (35b) and accumulate it. The grid points in a given angular segment and in the absorbing zone are denoted r n in Fig. 6, with n an index vector that differs from one segment to the other.

A straightforward collection of grid points in a segment n tends to produce noisy results because the number of grid points per segment fluctuates. We therefore associate with each grid point a smoothing function S centered at r n which distributes the strength over a vicinity of order of grid spacing. This suffices to produce acceptably smooth distributions. The total PAD A(θ, φ) and the state-specific PAD A α (θ, φ) are computed one direction allows one to identify uniquely energy and momentum as k = e k √ 2m e E kin . The energy is read off from the phase oscillations of ϕ(r M , t) and the direction e k ⇒ Ω M from the position of the measuring point. The practical evaluation of PES is done in a two-step process. QDD prints out the phase information at the measuring point to a file which is afterward analyzed by an auxiliary program, for details, see the Supplemental material, section M.7.3.1 therein.

The formula (51) applies for weak and for strong fields, but still fails for extremely strong fields. The limits of validity depend on system and time structure of the pulse. To give an order of magnitude, an intensity for long laser pulses impinging on Na clusters is I ≈ 10 15 W/cm 2 . For details, see [136].

This evaluation of PES through the above explained phase sampling is extremely efficient. However, it is not yet applicable to dissipative dynamics as RTA, because that causes jumps in the wave functions and thus spoils phase information. A proper extension of the scheme has yet to be developed.

An example of a PES is given in figure 10. Test case is the H 2 O molecule, irradiated by an off-resonant laser pulse. In this off-resonant regime, the total ionization remains small (2 × 10 -3 ). The values of the s.p. energies (see inset of the figure) thus remain almost constant in time. The PES clearly exhibits sharp peaks (mind the vertical log scale) above the background signal. These peaks correspond to a map of the s.p. energy spectrum, shifted by a multiple of the laser frequency ω las . A PES is thus a standard tool to access the single particle energy spectrum of a system. For instance, the 1st, 3rd and 4th peaks, respectively at 2.6, 5.5 and 7.7 eV, come from a 2-photon absorption of states 2, 3 and 4. Analogously, the 2nd, 7th, 8th and 9th peaks, respectively at 1.7, 13.1, 16.9 and 19.1 eV, come from a 3-photon absorption of the four states. As expected, the higher the multi-photon process, the smaller the peak height. A final word is in order on the peak widths. The longer the laser pulse, the better the resolution of the peaks. Still, there is physics contained in the width of the peak. In particular, this way to obtain such a PES, that is from the calculation of the ionization dynamics, allows the encoding in the peak width of the dynamical coupling between electrons but also between the electrons and the ionic background. In other words, a PES contains at the same time static information but also valuable dynamical ones at the side of the electron emission. This technique is therefore superior to the simple production of a PES by shifts of the energy spectrum and an artificial widening of the peaks by, e.g. convolution of a Gaussian with a constant width. Note also that this allows one to compute a PES also in the on-resonant regime (not shown here), then bearing all the impacts of a large electron emission can have on the irradiated system, in particular on the time evolution of the s.p. energy spectrum.

Observables specific to relaxation

Most of the observables computed with RTA are exactly the same as for TDLDA, e.g., energy, density, excitation spectra, ionization, or PAD. New are observables related to the mixed character of the one-body operator which is characterized by the occupation numbers w α . A specific quantity in that respect is the one-body entropy which is computed in diagonal representation (18) by the standard expression [START_REF] Reichl | A Modern Course in Statistical Physics[END_REF]

S = - α [w α log w α + (1-w α ) log(1-w α )] (52)
in units of Boltzmann constant. It serves as a direct indicator of thermalization and allows one to read off the typical time scale of relaxation processes. An example will be discussed along with Fig. 12, when a comparison of a RTA and a TDLDA dynamics is performed. Note that S remains strictly zero at TDLDA level with constant values of w α = 0 or 1.

The structure of the TDLDA and RTA packages

The QDD package, providing TDLDA as a basis with the more elaborate RTA on top, is a rather complex collection of routines. A rough schematic overview is given in Fig. 11. It shows the basic structure and switches in a self-explaining manner. A more detailed tree structure for the code and its subroutines at a deeper technical level is given in the Supplemental material.

Examples of RTA dynamics

A distinctive feature of RTA dynamics is that it accounts for electronic dissipation which manifests itself by a gradual transformation of the available excitation energy into heat. We will illustrate this process with two of the initial excitations available in QDD, namely an initial instantaneous boost on the one hand and a femtosecond laser pulse on the other hand. In addition to the standard observables accessible in TDLDA, we will study the time evolution of the one-body entropy S defined in Eq. ( 52) as a specific quantity characterizing dissipative features. Indeed, S is strictly zero in a TDLDA evolution. Any non-vanishing value of S thus provides a direct indicator of dissipation.

Electronic response to an initial boost

A simple way to visualize RTA dynamics is to consider an instantaneous excitation and to follow in time the dynamics of relaxation of the electrons. Such a perturbation can be delivered by an instantaneous initial boost of the electron cloud with respect to the ionic background, see section 2.2.5.3.

We take here such an initial boost along the z direction applied on Na + 11 as an illustrative example. The ionic background is here described by a soft jellium model, body entropy, the total ionization and the dipole moment along z for the three E * as indicated, on the one hand obtained in TDLDA (dashed curves) and on the other hand in RTA (full curves). Note that we chose to plot the envelop of the dipole amplitude rather than the full signal for a better readability and that we plot it in logarithmic scale to point out the attenuation of the signal. Without entering details, which can be found elsewhere [62,57], we briefly discuss the general outcomes. The dipole signal (upper panel) shows a fast, near exponential decay in the early phase up to 10 fs. This is the same for RTA and TDLDA thus having nothing to do with dissipation. It stems from Landau damping which is already present in TDLDA [33]. Dissipation becomes effective in the second phase and here we see clear differences. The TDLDA signal goes on undamped while RTA leads to further, slow decrease. The effect of dissipation is also seen in the ionization signal (middle panel). It reduces electron emission because part of the excitation energy is converted into intrinsic energy thus being not available anymore for direct electron emission. The most striking difference occurs in terms of the entropy S. It vanishes for TDLDA and acquires significant values for RTA. One observes a clear increase of S as a function of time, relaxing sort of exponentially toward an equilibrium value. This asymptotic value of S increases with initially deposited excitation energy E * and the typical relaxation time needed to reach the asymptotics in turn decreases with E * , features which are expected in such a dissipative theory.

Electronic response to a laser

We here illustrate how RTA practically works in case of a realistic laser irradiation. We take the H 2 O molecule as test case and consider the electronic response to laser pulses at two different frequencies [137]. The water molecule has an optical response with a well identified strong transition around 10 eV, well below its IP (around 15 eV) followed by a highly fragmented response spectrum mostly in the continuum [137]. We compare RTA and TDLDA for two laser frequencies: one on resonance with ω las = 11.4 eV, and one off-resonance with ω las = 10.2 eV, the other laser characteristics being maintained. The laser polarization axis goes along the z direction, while the H 2 O molecule lies in the xy plane, see Fig. 8 for the orientation of the molecule in space.

Figure 13 compares the time evolution in terms of three basic observables: electronic dipole, absorbed energy E abs as defined in Eq. ( 41), and ionization. Let us first discuss the resonant case (right panels) corresponding to a laser frequency slightly above the lowest excitation energy mode (10.2 eV). The early stages of irradiation progressively ionize the molecule which produces a slight blue-shift of the resonance and thus places the system, still during the pulse, exactly on top of the resonance with all consequences of resonant excitation. These are particularly well visible in TDLDA, namely enhanced dipole response whose envelope deviates soon from the pulse envelope and which continues to oscillate long after the pulse is over leading also to ongoing electron emission. The mechanism is well known and was used to analyze expansion of clusters irradiation by strong lasers [138,139]. Even if the mechanism is well understood, it has mostly been described theoretically using TDLDA-like approaches, see e.g. [140]. The major interest of the case here is the striking difference between RTA and TDLDA dynamics during the second phase of the response. While TDLDA dynamics leaves long lasting dipole oscillations, dissipation in RTA manages to curb down the dipole signal by turning energy away from the dipole channel and distributing it over the great manifold of other excitation modes. A remarkable feature is seen in the energy absorption from the laser field (right middle panel). TDLDA shows strong oscillations between absorption and (induced) emission of energy very similar to Rabi oscillations [141]. RTA reduces these oscillations because attenuation of the dipole reduces emission. This, in turn, gives way to absorb more energy from the laser field.

The situation is clearly different in the off-resonance case illustrated in the left column of figure 13. The electron dipole follows the external laser field closely and dies out once the latter is over. Emission is reduced as compared to the resonant case. Most striking is the comparison between TDLDA and RTA. At variance with the on-resonance case, differences between the approaches are much reduced. Dipole signals in particular are very similar. A larger difference is seen for energy absorption where RTA allows the system to soak more energy from the laser, similar as in the resonant case. Only small differences can be spotted on ionization, with RTA emitting slightly more because it absorbs much more energy (most of that used up for internal heating). Altogether, dissipation effects are much more pronounced near resonance and less so off resonance. This example, varying laser frequency, shows that dissipation depends sensitively on the details of the dynamical scenario. The great versatility of laser pulses offers many more variations of properties and a rich field for explorations [107].

Impact of dissipation on PAD

There is one more observable which we have not yet looked at in the context of RTA. This is the PAD introduced in section 3.2.6.3. To that end, we take up the test case from section 4.1, namely Na + 11 described with a soft jellium model, excited by an initial boost along the z direction. The Na + 11 cluster is a cylindrically deformed system and the laser is aligned with the symmetry axis. Thus the PAD does not depend on the angle φ and it suffices to look at the ϑ dependence. (less symmetric systems usually require some angular averaging [142]). Figure 14 and TDLDA (dashed curves), obtained after initial boost excitation in Na + 11 . Results from three different boost momenta, and consequently three different initial excitation energies E * as indicated, are compared. The stronger the excitation, the more forward peaked are the PAD. This is a well known and plausible feature because the stronger the field, the more it pulls electrons along its direction [36]. Comparing RTA with TDLDA results, we see that dissipation modeled in RTA reduces the anisotropy (trend to forward peaking). But this is a rather weak effect for the extremely short excitation here. Two aspects are to be remarked here: First, dissipation requires time to unfold and the longer a process takes, the more effective it is. Second, we have not added the isotropic background from thermal emission coming much later. This contribution reduces the anisotropy substantially [143].

Conclusions and perspectives

We have presented in this paper the first, to the best of our knowledge, open source package QDD (Quantum Dissipative Dynamics) allowing the account of dissipative electronic dynamics in finite systems, from atoms and molecules to clusters. The present implementation of dissipation is realized at the simplest level of a quantum Relaxation-Time Approximation which nevertheless includes essential features of dissipative dynamics. QDD has been constructed as an extension of a real-time real-space TDLDA package developed over the past two decades, in order to study off-equilibrium dynamics in clusters and molecules. As such, and to comply with major experimental signals such as PES and PAD, the TDLDA package already contained these key observables characterizing electronic emission. One objective of this paper is to present in a practical manner how these observables can be computed at the TDLDA level. These observables can, for most of them, be exported to RTA dynamics. The evaluation of PES in RTA is still a pending question, the difficulty being due to the technique used for computing the PES. Work along that line will be pursued.

In terms of physical properties, we have seen that the RTA dynamics displays the expected features of a dissipative approach. It has revealed quite different from TDLDA in many situations, the larger the excitation energy, the larger the effect. This is especially true when the system is driven towards one of its eigenfrequencies. One then observes both quantitative and qualitative differences between RTA and TDLDA. Differences shrink in fully off-resonance scenarios, although simulations remain quantitatively different on some observables. The implementation of dissipative features via RTA on top of TDLDA thus definitely appears as a crucial step forward in the theoretical modeling of far-off equilibrium dynamics such as encountered in numerous irradiation scenarios. This is true for cases in which one observes significant differences between TDLDA and RTA, as is obvious. But it is also true for cases in which TDLDA and RTA deliver qualitatively comparable results. Indeed, if one can conclude that TDLDA is an acceptable approximation in such cases, it is impossible to draw this conclusion without the input of a more elaborate approach such as RTA.

At variance with other open source real-time TDLDA packages, our QDD approach has not been focused on a large variety of otherwise well established parameters. By this, we mean a wide set of different DFT functionals or various kinds of pseudopotentials. This is justified in many respects. First, the whole spirit of the QDD project has always been to focus on far-off equilibrium scenarios for which structural details can be energetically overlooked. The set of parameters used in QDD have of course been checked against experimental or theoretical results when available, and allowed us to reproduce them in a perfectly acceptable manner. Such a reduced set of parameters thus allows us to focus on dynamical questions. It should also be noted that QDD routinely includes a simple strategy for the Self-Interaction Correction (SIC) strategy, namely Average-Density SIC (ADSIC). It restores key structural properties such as the IP which is crucial to obtain correct electron emission properties in real-time dynamics of irradiation scenarios.

This first release of QDD contains the full capabilities of TDLDA, particularly concerning all observables of electron emission. In the regime beyond TDLDA, it concentrates on simple and robust extensions, ADSIC for the self-interaction correction and RTA for dissipative dynamics. Extensions are in preparation for both approximations concerning SIC as well as dissipation. For the latter, we work on an implementation of Stochastic Time-Dependent Hartree-Fock (STDHF) which deals with an ensemble of mean-field trajectories to represent incoherent dynamical correlations. It includes by construction mean-field fluctuations, which opens the door to elaborate dynamical scenarios involving dissociative processes and, generally speaking, bifurcation dynamics. At the side of technical developments, the next steps are implementation of an MPI parallel version and a GPU version, the latter being already under way. At the formal side, we will prepare in near future a release which includes also a full time-dependent SIC as a further option, following our prototype implementation [144].

A simple signature of shape can already be read off here: Axially symmetric systems aligned along z-axis are distinguished by α lm =0 = 0.

However, the example of axial systems indicates a next problem, namely the orientation of the system. The appearance of α lm =0 = 0 can have two causes: first, the system is not aligned along z-axis, and second, there is a true breaking of axial symmetry. The first action is then to determine the principle axes of the system and to rotate the system such that the z-axis is identical with the dominant principal axis. This leads by definition to α 2 ±1 = 0 and α 2 +2 = α 2 -2 . Axial symmetry is truly broken if we then still find α 22 =0 = 0, signaling triaxial shapes. One often regroups the quadrupole deformation parameters into total deformation β 2 and triaxiality γ as This convention has been originally introduced to characterize shapes of nuclei [145] and has been taken over for clusters at several places, e.g. [START_REF] Lauritsch | [END_REF]149,150]. It is to be noted that γ = 0 • as well as γ = 60 • represent axially symmetric shapes. The case γ = 0 • corresponds to prolate shapes and γ = 60 • to oblate ones. The same definitions apply to ionic shapes if we replace dr (r) r i r j by I R I,i R I,j in Eqs. [START_REF] Marques | Fundamentals of Time-Dependent Density Functional Theory[END_REF]. We first determine eigenvalues and eigenvectors of N through whereby we take care to order the states to a sequence of decreasing occupations W (nat) A . This one-body density is still given in a space which is twice as large as the s.p. space before mixing. But the occupations for large A are extremely small. This allows to reduce R back to standard size and in natural orbitals as (C. 16)

N αβ V βl = V αl n l (C.

C 3 ,

 3 12 e - Optical response (a.u.) Frequency hω (eV) xy z s.p. energies (eV) Frequency hω (eV) xyz Optical response (a.u.) Frequency hω (eV) x y z s.p. energies (eV) Frequency hω (eV) xy (×20) z

Figure 1 :

 1 Figure 1: Ground state properties of 4 typical systems attainable with our open package : H 2 O, C 3 , Na +41 and C 60 . We show ionic structure and electronic cloud, sequence of single particle energies (with the HOMO emphasized in blue and the LUMO in red dashes) and the optical response in some specific spatial directions as indicated.

4 density--Figure 2 :

 42 Figure 2: Sketch of the scheme for performing one large time step t 0 -→ t 1 = t 0 +∆t in solving the RTA equations. The numbers in open circles indicate the steps as outlined in the text.

Figure 3 :

 3 Figure 3: Sketch of the scheme for solving the DCMF equations (23). This scheme expands step 4 from the RTA scheme presented in figure 2. The symbol O stands for orthonormalization of the new set of s.p. wave functions and D for the damping operator in the accelerated gradient step, see Eqs. (29), however, with different damping parameters E 0,damp → E (RTA) 0,damp and δ → δ (RTA) .

3 .

 3 Numerical aspects: methods, inputs and outputs 3.1. General numerical schemes 3.1.1. Grid representation and derivatives 3.1.1.1 Grid definition

Figure 4 :

 4 Figure 4: Logarithmic plot of three different convergence criteria for static convergence: the variance of s.p. energy (31) and the two density criteria (32), versus iteration number n. Test case is the H 2 O molecule.

Figure 5 :

 5 Figure 5: Comparison of electronic propagation schemes taking as test case the ground state of Na + 9

'

  ↵ ( r , t ) = 0 Measuring point used in PES mask function M(r)

  l b s b 9 5 j 7 a n u N q K / k 3 r 5 x E o M i f 1 L N 8 3 8 r 0 7 V I t H H s a 7 B p Z o i z a j q W O q S 6 K 6 o m 5 t f q p L k E B G n c I / i g j D T y m m f T a 2 J d e 2 q t 7 a O v + l M x a o 9 S 3 M T v K t b 0 o B L P 8 c 5C + r l Y u m g W D 4 / L F R O 0 l F n s Y N d 7 N M 8 j 1 D B G a q o k b f A I 5 7 w b F w a I + P O u P 9 M N T K p Z h v f l v H w A W I A l X E = < / l a t e x i t >Angular segment for PAD Integration point used in PAD r n < l a t e x i t s h a 1 _ b a s e 6 4 = " w i 2 S z p / S U y k C s / 7 m x I 7 w X 3 J z T o 0 = " > A A A D B 3 i c j V H L S s N A F D 2 N 7 / q q u n Q T L I I L K U k V 7 F L q Q p c q V g u t l C S d a j A v J h N B Q j 7 A P 3 H n T t z 6 A 2 7 V P 9 C / 8 M 4 4 B R + I T k h y 5 t x z z s y d c ZP A T 4 V l v Z S M k d G x 8 Y n J q f L 0 z O z c f G V h 8 T i N M + 6 x l h c H M W + 7 T s o C P 2 I t 4 Y u A t R P O n N A N 2 I l 7 s S P r J 5 e M p 3 4 c H Y m r h J 2 G z l n k D 3 z P E U T 1 K s 1 u n w 3 I q 5 J y U o o i P 9 x t F r m 1 b t c b 6 1 Z R 7 n 4 q m d 3 Q E e f u w O Q 9 j f K o 6 F W q V s 1 S w / w J b A 2 q 0 G M / r j y h i z 5 i e M g Q g i G C I B z A Q U p P B z Y s J M S d I i e O E / J V n a F A m b w Z q R g p H G I v6 H t G s 4 5 m I 5 r L z F S 5 P V o l o J e T 0 8 Q q e W L S c c J y N V P V M 5 U s 2 d + y c 5 U p 9 3 Z F f 1 d n h c Q K n B P 7 l 2 + o / K 9 P 9 i I w Q E P 1 4 F N P i W J k d 5 5 O y d S p y J 2 b n 7 o S l J A Q J 3 G f 6 p y w p 5 z D c z a V J 1 W 9 y 7 N 1 V P 1 V K S U r 5 5 7 W Z n i T u 6 Q L t r 9 f 5 0 9 w X K / Z G 7 X 6 w W Z 1 u 6 m v e h L L W M E a 3 e c W t r G H f b Q o + w a P e M K z c W 3 c G n f G / Y f U K G n P E r 4 M 4 + E d P n W q 0 g = = < / l a t e x i t > r M < l a t e x i t s h a 1 _ b a s e 6 4 = " s W 8 C 3 S 1 f 7 g W a h B g 1 9 y Y y w G z z 0c E = " > A A A C 5 3 i c j V H L S g M x F D 2 O 7 / q q u n Q z W A R X Z a q C L o t u 3 A g K 1 h a s l E x M N Z h 5 k M k I U r p 3 5 0 7 c + g N u 9 U / E P 9 C / 8 C a O o B b R D D N z c u 4 9 J 7 n 3 h q m S m Q m C l y F v e G R 0 b H x i s j Q 1 P T M 7 V 5 5 f O M q S X H P R 4 I l K d C t k m V A y F g 0 j j R K t V A s W h U o0 w 4 s d G 2 9 e C p 3 J J D 4 0 V 6 k 4 i d h Z L L u S M 0 N U p 7 z c d h 6 9 U O W i 7 7 c j Z s 7 D r q 8 7 D n G m e n v 9 T r k S V A O 3 / E F Q K 0 A F x d p P y s 9 o 4 x Q J O H J E E I h h C C s w Z P Q c o 4 Y A K X E n 6 B G n C U k X F + i j R N q c s g R l M G I v 6 H t G u + O C j W l v P T O n 5 n S K o l e T 0 s c K a R L K 0 4 T t a b 6 L 5 8 7 Z s r 9 5 9 5 y n v d s V / c P C K y L W 4 J z Y v 3 S f m f / V 2 V o M u t h y N U i q K X W M r Y 4 X L r n r i r 2 5 / 6 U q Q w 4 p c R a f U l w T 5 k 7 5 2 W f f a T J X u + 0 t c / F X l 2 l Z u + d F b o 4 3 e 0 s a c O 3 n O A f B 0 V q 1 t l 5 d O 9 i o 1 L e L U U 9 g C c t Y p X l u o o 5 d 7 K N B 3 t d 4 w C O e P O n d e L f e 3 U e q N 1 R o F v F t e f f v n P S d + Q = = < / l a t e x i t >

Figure 6 :

 6 Figure 6: Left: Schematic view of a coordinate space grid with an absorbing layer (red shaded zone), a sampling direction rn for accumulating the Photo-Angular Distribution (PAD), see Sec. 3.2.6.3, and a measuring point r M for the computation of a Photo-Electron Spectrum (PES), see Sec. 3.2.6.4. Right: mask function M(r) plotted as a function of r, the distance to the center of the numerical grid.

8 ω

 8 las = 2.6 eV T pulse = 25 fs I = 5 × 10 11 W/cm 2 10

Figure 7 :

 7 Figure 7: Illustration of convergence of the RTA step for the test case Na 8 with jellium background excited by a laser pulse with frequency ω las = 2.6 eV, duration T pulse = 25 fs, and intensity I = 5 × 10 11 W/cm 2 . Left column: as functions of time, z-dipole moment for RTA dynamics compared with mere TDLDA (bottom) and number of DCMF iterations in each RTA step (top). Right column: as functions of the number of DCMF iterations, deviation of local density from goal (bottom) and average variance of s.p. energies (top), both for three times as indicated and highlighted by vertical lines in the upper left panel.

Figure 8 :

 8 Figure 8: Illustration of ionic motion in H 2 O following an instantaneous electronic boost along y direction with initial energy E boost = 27.2 eV, with the molecule lying in the xy plane (top middle inset).Propagation was done with δt = 0.6 as and δt I = 10δt = 6 as. Plotted as functions of time are: left column, ionization (top) and electronic dipole moment in each spatial direction (bottom); middle panel, ionic kinetic energy E kin , total and its three spatial components; right column, ionic velocities (bottom) and positions (top) in x direction.

ω 3 D

 3 las = 5.4 eV, I = 10 14 W/cm 2 , T pulse = 12 fs C

Figure 9 :

 9 Figure 9: The C 3 chain molecule excited by a laser pulse with ω las = 5.4 eV, I = 10 14 W/cm 2 , T pulse = 10 fs, and with polarization along the z direction (= molecule axis), leading to an average ionization of 0.8 electrons. As functions of time: (bottom left) Dz, dipole moment in z direction, (top left) total ionization Nesc and ionizations per state Nesc,α and (right) charge-state probabilities P (+k) in logarithmic scale.

H 2 O I = 3 ×

 23 Photo-Electron Spectrum (arb. units)

Figure 10 :

 10 Figure 10: Photo-Electron Spectrum obtained in H 2 O, in the xy plane (see inset of Fig.8), irradiated by a laser pulse with characteristics as indicated. The laser polarization is along the z direction. The single particle energy spectrum is shown as an inset. The short vertical lines at the bottom indicate the single particle energies shifted by multiples of ω las .

Figure 11 :

 11 Figure 11: Schematic flow diagram of the QDD code.

EFigure 12 :

 12 Figure 12: Time evolution of the one-body entropy (bottom), the total ionization (middle) and the envelop of the dipole moment in z direction (top, vertical log scale) of Na + 11 after various initial boosts delivering excitation energies E * as indicated. Full curves are obtained in RTA while dashed curves are obtained with TDLDA.

ωω

  laser = 11.4 eV I laser = 5.6 × 10 13 W/cm 2 T pulse = 36 fs laser = 10.2 eV I laser = 5.6 × 10 13 W/cm 2 T pulse

Figure 13 :

 13 Figure 13: Comparison of TDLDA (blue) and RTA (red) time evolution of three basic observables for H 2 O: dipole moment perpendicular to molecular plane (bottom), absorbed energy (middle) and total ionization (top). The water molecule was excited by a laser pulse with frequency ω las = 11.4 eV for the off-resonance case (left column) and ω las = 10.2 eV for the on-resonance case (right column). For both frequencies, the laser intensity is I = 5.6 × 10 13 W/cm 2 , and the total pulse length T pulse = 36 fs. The water molecule is in the xy plane and the laser polarization is along z, see inset of figure 8.

  Photo-Angular Distribution

Figure 14 :

 14 Figure 14: Photo-Angular Distributions after initial boost along z applied in Na + 11 , for three different boosts delivering three excitation energies E * as indicated. The dashes show the results obtained in TDLDA, the solid curves those obtained in RTA.

C 1 ) 2 )= 5 )

 125 . Handling of basis sets in RTA C.1. Transformation to natural orbitals Before entering the RTA step, we have a one-body density matrix obtained from a couple of TDLDA steps ρ ≡ ρmf =The absorbing boundary conditions remove electron content from the wavefunctions which thus lose their strict ortho-normality leading to a non-unit norm matrixN αβ = ϕ (mf)We first need to represent the actual one-body density matrix in terms of natural orbitals δ ij . The loss of electron content is taken up by a change of the occupations to w (nat) i . The natural orbitals are defined by the eigenvalue equation ρ|ϕ The task is now to determine the transformation matrix W αi . We insert the expansion (C.5) into the eigenvalue equation (C.4) and obtain a matrix equation ρ αβ W βi = N αβ W βi w

7 )

 7 and use that to build the inverse root of N as(N -1/2 ) αl = V αl n -1/2 l . (C.8)With that, we can map the eigenvalue equation (C.6) into standard form(N -1/2 ) † lα ρ αβ (N -1/2 ) βl U l i = U li w eigenvectors U li , we realize that ρ αβ (N -1/2 ) βl U l i = (N +1/2 ) αl U li w (nat) i = N αβ (N -1/2 ) βl U li w (nat) i . (C.10)and conclude that the searched transformation for the expansion (C.6) is finally composed asW αi = (N -1/2 ) αl U li . (C.11)The new occupation numbers w (nat) iwere already given from solving the eigenvalue equation (C.9)C.2. Mixing of two one-body matricesAt the end of one RTA step, we mix two one-body density matrices as ρ(mix) = (1η)ρ (mf) + η ρ(eq) , (C.12) see step 5 in section 2.3.2.2. The aim is to produce eventually ρ(mix) in the basis of its natural orbitals. The problem is that the constituents ρdifferent sets of s.p. states from which {ϕ (mf) α, α = 1...Ω} is not even ortho-normalized.We combine the two sets of s.p. states to one super-set {φ A , A = 1...2Ω} asφ A = ϕ (mf) α for A = α ∈ {1...Ω} ϕ (eq) α for A = α + Ω , α ∈ {1...Ω} (C.13a) W A = w (mf) α for A = α ∈ {1...Ω} w (eq) α for A = α + Ω , α ∈ {1...Ω} (C.13b)This yields a one-body density in double-spaceR = 2Ω A=1 |φ A W A φ A | (C.14)with a basis {φ A } which is not northo-normalized. That is exactly the task to find natural orbitals which was solved in section C.1. We undertake the same steps (of course, in doubled space) and end up with

Acknowledgments:

We thank Institut Universitaire de France, French ANR contract CIRCE (Contract No. ANR-16-CE30-0012-02), French ANR contract BAMBI (Contract No. ANR-18-CE30-0009-01), Institut de Recherche Technologique Saint Exupéry, Sciences et Technologies pour l'Aéronautique et l'Espace Foundation, the CNRS and the Midi-Pyrénées region (doctoral allocation number 2017 001852) for support during the realization of this work. One of us (ES) thanks Queen's University of Belfast for its hospitality during the realization of part of this work. It was also granted access to the HPC resources of CalMiP (Calcul en Midi-Pyrénées) under the allocation P1238, and of RRZE (Regionales Rechenzentrum Erlangen).

QDD CPC Library link to program files: (to be added by Technical Editor) Developer's repository link: None Code Ocean capsule: (to be added by Technical Editor)

as :

A α (θ, φ) = n∈abso.zone dr r 2 S(re rr n ) n esc,α (r n ) ,

where e r = (sin θ cos φ, sin θ sin φ, cos θ) is the unit vector in the direction of the wanted angles. The smoothing is done by tent functions which comply with the integration rule used in the normalization. The angular segments in Fig. 6 schematically symbolize this smoothing which collects (weighted) information in the vicinity of a ray. An example of PAD is given in section 4.3 in connection with RTA dynamics.

Photo-electron spectra

A Photo-Electron Spectrum (PES) can be deduced from the temporal phase oscillations of the s.p. wave functions ϕ α (r M , t) at measuring points r M close to the absorbing bounds, see Fig.

6 [START_REF] Pohl | [END_REF]. The result for the PES sampled at r M is

Note that the formula does not only give the total PES, but also the PES Y α for emission specifically from state α. In Eq. ( 51b), e 0 = E 0 /|E 0 | is the direction of (linear) polarization of the electrical field E 0 and Ω M the solid angle associated with r M . This form applies for the wave function ϕ α (r M , t) in length gauge as computed in the code. In Eqs. (51c) and (51d),

is the time integrated laser pulse envelop introduced in Eq. (9).

A detailed derivation of Eqs. ( 51) is found in [136]. We summarize here the ideas behind that compact formula. It is deduced under the assumption that the measuring point r M is sufficiently far away from the system (placed around r = 0) such that an outgoing electron wave has direction e k = e M = r M /r M . We expand the wave function at r M into outgoing waves with direction e k and momenta k > 0. These are plane waves e ike k •r M for weak fields and for stronger fields the corresponding electron waves in the time-dependent laser field (Volkov states). The fact that we have only outgoing waves in

A. Iterative correction of total energy

The RTA scheme as summarized in figure 2 The adjustment is done as if an equilibrium Fermi distribution (23d) were to be changed slightly. The Fermi distribution is obtained by maximizing entropy

with the requirement (constraint) to match total particle number and s.p. energy

while the actual status is

Linearizing this for small changes δw α leads to a simply solvable system of linear equations. This scheme was used in our early implementations of RTA [62,107]. However, that tends to overshoot occasionally which calls for costly re-iterations. A more stable scheme is obtained by modifying the maximal entropy principle to a quadratic criterion calling for minimal change of entropy in the form

together, again, with the constraints (A.2b). This amounts to the variational condition

3) Linearizing that yields in straightforward manner the occupation change as

with the Lagrangian parameters obtained by the solution of the linear equations

where δN, δE sp , are the differences between status and target values of N, E sp , as outlined in Eqs. A.2b and A.2c, and s denotes the derivative of s with respect to w (s = ds/dw). Note that this rather robust scheme is applied for the spin-up and spin-down subsets separately.

B. Dimensionless quadrupole moments

In section 3.2.2, we introduced the basic form parameters of a density distribution. They are plagued by the inconvenience to scale with system size. To concentrate on shape, the concept of dimensionless quadrupole moments was used since long in nuclear physics [145,[START_REF] Bohr | Nuclear Structure[END_REF]. It turned out also be useful in cluster physics [128,89] and surely helps to sort the gross structure of molecules. We review here briefly the steps from Cartesian quadrupole moments to dimensionless ones.

The first step is to map Cartesian quadrupole moments into spherical quadrupole moments Q 2µ [START_REF] Edmonds | Angular Momentum in Quantum Mechanics[END_REF] Q 20 = 5 2π