René Kalus
email: rene.kalus@vsb.cz

Ivan Janeček
email: janecek@ugn.cas.cz

Florent X Gadéa
email: gadea@irsamc.ups-tlse.fr

Florent Xavier Gadéa

to ionic rare-gas clusters photodissociation, post-ionization fragmentation, and collisions

Keywords: combining Ehrenfest, decoherence, and multiscale approaches applied non-adiabatic dynamics, mean-field approach, diatomics-in-molecules, rare-gas cluster cations

Introduction

Ionic rare-gas clusters have been, during the last decades, a wonderful playground and a veritable laboratory for exploring and questioning non-adiabatic dynamics. Initially, mainly for fundamental purposes since it was shown that the charge was localized in a small subunit of the clusters [1] solvated by almost neutral atoms. Although made of identical atoms, the ionic rare-gas clusters are representatives of complex systems bound by heterogeneous forces, strong chemical ones within the core, intermediate polarization forces between the core and the surrounding atoms, and weak van der Waals forces between the remaining neutral atoms. Many experiments have also been performed for these attractive systems [START_REF]Clusters of Atoms and Molecules[END_REF] addressing the question how this complexity is reflected during their fragmentation after energy is pumped into the system, either by photon absorption, collisions, or sudden ionization. More recently, applications of ionized rare gases and cold rare-gas plasmas of medical interest [START_REF] Laroussi | [END_REF]4], for UV light [5], space craft propulsion [6], astrophysics [7], or even in environmental sciences and geosciences [START_REF] Burnard | The Noble Gases as Geochemical Tracers[END_REF] have led to a renewed interest in the field.

From the theoretical point of view, interest has also been nourished by the tractable complexity of the ionic rare-gas clusters. Complexity because of their heterogeneous nature typical for various fields of research from solvation to solid state, tractable because it was shown that the electronic Hamiltonian can be, for these species, very effectively modeled by rather simple diatomics-in-molecules (DIM) models of rather small dimensions [1,[START_REF] Amarouche | [END_REF], 3N and/or 6N (N being the number of atoms in the cluster) if the spin-orbit (SO) interaction is not or is taken into account, respectively. Note, in particular, that the inclusion of the SO interaction, which is unavoidable for the heavier rare gases, is rather simple within the DIM approach [START_REF] Amarouche | [END_REF].

In the 1980s, non-adiabatic dynamics was mainly based on trajectory surface hopping (TSH) approaches [10,11,12,13]. At that time, the surfaces needed to be previously investigated, computed, sampled, and the possible areas for surface hopping located. Only later, on-the-fly TSH calculations become practicable. Among the first on-the-fly approaches, a hemiquantal dynamics method [14] was proposed with the use of the whole DIM basis set (HWD) coupling the mean-field (Ehrenfest) approach [15] with a classical treatment of nuclei. Within this approach, the need for a preliminary potential energy surfaces exploration was avoided. Very interesting results where found for example for the argon trimer photodissociation [16,17,18]. However, the well known problem of the mean-field approach of unphysical mixing of reaction (fragmentation) channels, reflected here in fractionally charged fragments, limited the approach to the direct fragmentation resulting in separate monomers. Later, our team introduced a periodic decoherence mechanism [19] to fix this difficulty and to deal with sequential fragmentation processes, particularly found in the post ionization dynamics of rare-gas clusters. Further on, clusters remaining in excited electronic states for extremely long lifetimes faced the theoretical treatment to another computational limitation. As an approximate solution, we proposed a multiscale treatment [20,21] including radiative and non-radiative transitions.

In the present contribution, our theoretical approaches and main obtained results will be recalled. First, a methodology overview will be presented in Sec.

2, while Sec. 3 will be devoted to a review of main results we obtained for photoabsorbtion, photodissociation, post-ionization fragmentation, and collisions of rare-gas cluster cations. Some new results will be also presented for the luminescence spectrum obtained from the multiscale treatment and for collisions.

Methods and computations

Dynamics

The dynamical method used in our calculations relies on a classical treatment of heavy atomic nuclei and a full quantum treatment of light electrons. As a consequence, it is called the hybrid method throughout this paper. In the literature, the approach used in this work is often referred to as the Ehrenfest (or mean-field ) method. However, we adhere to the hybrid denomination here to emphasize its mixed, quantum-classical nature.

In addition, the inclusion of quantum decoherence (as described below in subsection 2.2) goes well beyond the usual mean-field approaches.

Note also that alternative names like hemiquantal method [14,16,22,23,24,25,26,27] and semiclassical method [19,28,29,20,30,21,31] were also used in our previous studies.

Within the hybrid method, the equations of motion are written using the mean-field formalism [15], which results in a set of coupled classical Hamilton equations for the nuclei,

qi = p i m i , ṗi = ψ| - ∂ Ĥ ∂q i |ψ , (1)
and time dependent Schrödinger equation for the electrons,

ih ∂|ψ ∂t = Ĥ|ψ . (2)
In this paper, small Latin indices are used to label nuclear degrees of freedom (i = 1, ..., 3N where N is the total number of atoms)1 and q i and p i denote nuclear coordinates and momenta, respectively. Further, Ĥ is an electronic Hamiltonian, which is supposed to parametrically depend on the nuclear positions, |ψ is a time dependent electronic wave function, and the angle brackets in the second part of Eq. 1 denote an integration over electronic positions.

For a numerical treatment, the (time-dependent) electronic wave function and the electronic Hamiltonian are expanded against a properly chosen basis set, |Φ α ,2 which may also parametrically depend on nuclear coordinates:

|ψ(t) = α c α (t)|Φ α (q i (t)) , (3)
Ĥ = β,γ H βγ |Φ β Φ γ |. (4)
In our calculations, the basis set wave functions have been assumed orthonormal,

Φ α |Φ β = δ αβ ,
ṗi = - α,β c * α c β ∂H αβ ∂q i . (5)
Similarly, the electronic Schrödinger equation can be rewritten in an expanded form after Eqs. 3 and 4 are inserted in Eq. 2,

ih ċα = β H αβ c β . (6)
A note is probably worth to add here on Eq. 5. Since the Hamiltonian matrix and the expansion coefficients of the electronic wave function are in general complex-valued, the form of Eq. 5 may become rather involved if their real and imaginary parts,

c α = c (re) α + ic (im) α and H αβ = H (re) αβ + iH (im) αβ , are explicitly considered, 4 ṗi = - α,β c (re) α c (re) β + c (im) α c (im) β ∂H (re) αβ ∂q i + c (re) α c (im) β -c (im) α c (re) β ∂H (im) αβ ∂q i . (7)
However, even if the Hamiltonian matrix is complex-valued, its imaginary part may be independent of nuclear coordinates (it holds, e.g., for ionic clusters of rare gases discussed in Subsec. 2.5, see also Ref. [32] for details) and

∂H (im) αβ ∂qi = 0.
In this case, the second term on the right-hand-side of Eq. 7 vanishes and the

equation is considerably simplified, ṗi = - α,β c (re) α c (re) β + c (im) α c (im) β ∂H (re) αβ ∂q i . ( 8
)
As mentioned for example in Ref. [14], the electronic Schrödinger equation, Eq. 6, must be treated with care since rapid time oscillations of the electronic wave function (i.e., expansion coefficients c α ) may lead, if the equation is integrated numerically, to extremely short time steps and huge computational demands. Such short time steps are, however, unnecessary for the much slower motion of nuclei. A procedure was thus proposed [32] to treat these rapid electronic oscillations in a computationally efficient way. More specifically, for a time step short enough for the slow moving nuclei to be considered, within the time step, fixed in space, the electronic equation of motion is integrated between times t and t + ∆ analytically,

|ψ(t + ∆t) ≈ exp - i h Ĥ(t)∆t |ψ(t) , (9)
which finally leads to

c α (t + ∆t) ≈ β,γ Q * βγ (t)Q αγ (t)exp - i h E γ (t)∆t c β (t), (10)
with E γ (t) being the eigenvalues of the electronic Hamiltonian calculated at time t and matrix Q contains corresponding eigenvectors, expanded against basis set |Φ α , as its columns. Though Eqs. 9 and 10 are only approximately valid within the approximation of fixed nuclei, it can be shown (see Ref. 32 for details) that the error introduced by this approximation is negligible if time step ∆t is sufficiently short.

As is clear form the preceding explanation, eigenvalue/eigenvector problem has to be solved repeatedly during a dynamical calculation. To facilitate respective evaluations in systems which exhibit fast dissociation, an improved simulation scheme has been proposed [33]. The scheme consists in periodical checks of the fragmentation state of the simulated system. If well separated fragments are detected (on the base of a simple distance criterion, for example), they are evolved separately and the rank of the Hamiltonian matrix to be diagonalized is reduced.

Quantum decoherence

It is well known that the Mean-Field (MF) approach described above suffers from a serious drawback consisting in unphysical mixture of adiabatic electronic states during the classical evolution of nuclei. For an ionic system (e.g., raregas cluster cations discussed in Subsec. 2.5), it may lead to fragmentary (noninteger) charges observed for MF trajectories upon the system fragments. In a preceding work [19], we proposed to remove this defect from the MF approach by introducing periodic attempts to collapse (or quench) the time-dependent electronic wave function to one of available adiabatic electronic states. The extended dynamical approach was denoted by an MFQ (Mean-Field with Quenching) acronym and was interpreted in terms of quantum decoherence 5 . Note also that the mean-field may be more appropriate for such a model inclusion of quantum decoherence since it completely neglects the decoherence effects, while quantum decoherence is partially, but to an unknown extent, included in the trajectory surface-hopping algorithms by considering changes of adiabatic surfaces in gradient calculations needed for the classical propagation of nuclear degrees of freedom.

The quenching algorithm comprises basically two steps. As a first step, probabilities are calculated for collapses of the current electronic wave function into all available adiabatic states and a particular collapse is proposed using the roulette-wheel algorithm. Secondly, nuclear velocities are adjusted so that the total energy of the system is conserved. The proposed electronic transition can be, in principle, rejected, mainly in the second step of the present algorithm due to insufficient kinetic energy of atomic nuclei needed to compensate upward electronic jumps. If this applies, the system resumes the coherent evolution until next hop attempt [34].

Two alternative approaches for calculating the probability of the electronic wave function collapse have been proposed [19]. Firstly, a numerically cheap calculation employing the adiabatic amplitudes of the current electronic wave function [34] has been considered. Following our preceding studies, we denote this approach by acronym "AMP" (AMPlitudes). More specifically, the normalized probability of collapsing the current electronic state, |ψ , to a particular adiabatic state, |φ α , is calculated from [34] g AMP ψ→φα = a αα ,

where a αα are diagonal elements of the electronic density matrix ( a βα = C β C * α ) and C α denote amplitudes of the current electronic wave function, ψ, with respect to the adiabatic wave functions, |ψ = α C α |φ α .

Secondly, Tully's fewest switches (TFS) algorithm [11] has been proposed to calculate the quenching probabilities 6 . In this algorithm, the probability of a transition from adiabatic electronic state |φ α to adiabatic state |φ β occurring within a time interval ∆t Q (the period of quenching attempts) is given by (see Ref. 11)

g TFS φα→φ β = 2∆t Q a αα 1 h Im(a * βα H βα ) -Re(a * βα Ṙ • d βα ) , (12)
where H βα ≡ φ β | Ĥ|φ α are matrix elements of the electronic Hamiltonian, R represents current nuclear configuration, and d βα ≡ φ β |∇ R φ α are nonadiabatic coupling vectors. For adiabatic states φ α and φ β , the general formula of Eq. 12 simplifies since the first term on its right-hand-side disappears 7 ,

g TFS φα→φ β = - 2∆t Q a αα Re(a * βα Ṙ • d βα ) = 2∆t Q a αα Re(a αβ Ṙ • d * αβ ). ( 13
)
6

The TFS algorithm may be considered more realistic than the AMP approach, the latter being, however, computationally less demanding. Since the AMP method leads to results which usually compare well with available experimental data (see section 3), it has routinely been used in our calculations as a computationally cheaper, but still sufficiently accurate alternative. 7 The following identities are used in the following equation to derive the second equality:

a * βα = a αβ and d βα = -d * αβ .
If a general transition from a general coherent state, ψ = α C α φ α is considered, related quenching probability, g TFS ψ→φ β , is calculated from Eq. 13 by weighting over participating adiabatic states. In our implementation [19], it is done in two steps. Firstly, the initial adiabatic state is chosen randomly with the probability calculated from the amplitudes of the current electronic wave function using Eq. 11, and secondly, Eq. 13 is employed for calculating the TFS transition probabilities from the selected initial state to all possible final states. The probabilities are then used to propose a particular transition. Note, however, that the probabilities given by Eq. 13 do not, in general, sum up to 1 and that the electronic jump can be banned in the TFS approach even at this stage. Then, the coherent evolution is resumed until next quenching attempt.

After the quenching step has been proposed via either the AMP algorithm or the TFS method, nuclear velocities have to adjusted so that the total energy of the system is preserved,

v ≡ Ṙ → v = v + ξx. ( 14
)
Several schemes of the nuclear velocity adjustment have been reported in the literature: e.g., adjusting along the non-adiabatic coupling vector (the real part of the vector in the present work) [11,35],

x = Re[d βα ], (15)
or adjusting along the vector given by the gradient of energy difference between the initial and final states, E β -E α [36,37,38],

x = ∇ R [E β (R) -E α (R)] . (16)
Both methods have been used in our calculations. In line with our previous papers, we denote them by acronyms "C" (non-adiabatic Coupling vector) and "G" (Gradient), respectively.8

In addition to the two literature methods of the post-quenching velocity adjustment, we have also used another, computationally much cheaper method which consists in scaling the velocities by a positive multiplicative factor,

v = ζv. (17)
In this way, the method follows some implementations of a thermostat in constanttemperature molecular dynamics calculations [START_REF] Allen | Computer Simulation of Liquids[END_REF]. We denote this approach by acronym "S" (Scaling). However, a drawback of this additional approach is that it, like the C method, does not conserve the total angular momentum of the nuclei unless the momentum is zero [19] and is thus inappropriate for calculations of collisions with a non-zero impact parameter. For this case, the S procedure has been modified in a way that the nuclear velocities are first transformed to a frame rotating around the center-of-mass of the collision complex, in which the total angular momentum of nuclei is zero, and only then scaled using the S method. Finally, the velocities are transformed back to the initial frame and coherent evolution is continued. We denote this modified S method by an "SR" acronym.

Combining the two methods for the calculation of the collapse probability, TFS and AMP, and the three (four) methods for an after-collapse velocity adjustment, C, G and S (SR), we totally obtain six (eight) MFQ models. In general, they are denoted by acronyms MFQ-XXX/Y, where XXX = TFS or AMP, and Y = C, G or S (SR). 9

In our implementations, periodically evoked quantum decoherence with a quenching period fixed throughout a particular calculation is used. In practice, the period of the decoherence attempts varies, however, with the quenching mechanism employed. For the MFQ-AMP conserve the total angular momentum of nuclei [39]. The G method, on the other hand, does not suffer from such a defect. See also a detailed discussion provided in Refs. 19, 39.

9 For example, MFQ-AMP/S denotes a model within which the collapse probabilities are obtained from adiabatic amplitudes of the current electronic wave function, Eq. 11, and the velocities are afterward adjusted by simple multiplicative scaling, Eq. 17.

approach, it is typically about one attempt per ∆t Q ≈ 10 -100 fs, while for the MFQ-TFS method, it is usually by an order of magnitude shorter (an attempt per ∆t Q ≈ 1 -10 fs) because the rejection rate is much higher (see, e.g., references 19,31). Note also that if a quenching attempt is performed each time step, the MFQ-TFS approach becomes very close to the Tully's Fewest Switches one.

Radiative transitions

The implementation of radiative transitions basically follows in our codes the scheme proposed for including quantum decoherence as described in the preceding subsection. It consists in periodical attempts to change an (adiabatic) electronic state of the simulated system in the downward direction (φ α → φ β , β < α), without adjusting nuclear velocities, however, since the released electronic energy is supposed to be emitted in the form of electromagnetic radiation. The probability for a φ α → φ β radiative transition is calculated from Fermi's golden rule [START_REF] Fermi | Nuclear Physics[END_REF],

g (R) φα→φ β = ∆t R 3πε 0 h4 c 3 (E α -E β ) 3 | µ αβ | 2 , (18)
where E α and E β are adiabatic electronic energies prior and after the transition, respectively, µ αβ is the transition dipole moment between corresponding states,10 and ∆t R is an (arbitrary) time period at which radiation transitions are attempted. In our calculations, ∆t R has been chosen equal to the period of electronic quenching in the MFQ algorithms. 11 In the present work, calculations with the radiative transitions included will be denoted by suffix "-R" added to the acronyms introduced in the preceding subsection (e.g., MFQ-AMP/S-R).

Multiscale method

The non-adiabatic dynamics is primarily intended for the investigation of fast processes with typical times up to t = 1 ns and its extension to much longer times requires rather huge computational demands. As a consequence, if metastable intermediates are formed with life-times by many orders of magnitude longer, dynamical simulations become impracticable, in particular for large systems.

As an example, let us mention the experimental observation of the long-time decay of metastable ionic rare-gas clusters (see, e.g., Ref. [21] and references reported therein).

Mainly for this reason, an approximate extension of the direct non-adiabatic dynamics calculations, which overcomes this limitation, has been recently proposed for radiative transitions [20] as well as for both radiative and non-radiative transitions considered [21]. This approximate methodology consists in extending the short-time non-adiabatic dynamics calculations by a kinetic rate theory treatment. More specifically, each trajectory resulting in long-lived metastables is treated as an ensemble of systems undergoing a long-time electronic relaxation which is represented by a properly proposed set of first-order "chemical" processes. The time evolution of abundances of resulting relaxed electronic states (and decay products deduced from them) is then governed by a set of kinetic equations which can be easily solved, either numerically or analytically, over arbitrary time intervals and at modest computational costs. Since many time scales can be treated in this way, this approach is called the multiscale method.

A general scheme of this multiscale method can be typically split into two phases and several partial steps [20,21]:

1. phase 1 -direct
Γ (R) αβ = 1 3πε 0 h4 c 3 (E α -E β ) 3 | µ αβ | 2 , ( 19
)
where the symbols used on the right-hand side are the same as in Eq. 18. The rate constant for non-radiative electronic transitions can then be obtained from Tully's fewest switches formula as simplified for adiabatic states (see Eq. 13),

Γ (NR) αβ = 2 a αα Re(a αβ Ṙ • d * αβ ). ( 20
)
with, again, the symbols used in Eq. 20 having the same meaning as in Eq. 13 of Subsec. 2.2.

To illustrate how the multiscale method works during the second, long-time phase, a simple model case is considered when only parallel first-order processes are expected to govern the long-time evolution. In this specific case, the kinetic equations used to calculate the time evolution of the occupation numbers of the initial state, n α , and all possible final states, n β , read as follows, 12 ,

dn α dt = - β Γ (R) αβ n α - γ Γ (NR) αγ n α , dn (R) β dt = Γ (R) αβ n α , dn (NR) γ dt = Γ (NR) αγ n α . (21)
Here, the sums on the right-hand-side of the first equation run over all final (adiabatic) states (φ β and φ γ ) allowed within a particular kinetic model for radiative and non-radiative transitions, respectively, and changes in the occupation numbers are treated separately for radiative (R) and non-radiative (NR)

processes. This is because, even though the electronic states reached after radiative and non-radiative transitions may be the same, nuclear kinetic energies and, consequently, post-transition fragmentation pathways may differ for the two cases.

If the initial conditions for Eqs. 21 are considered in a form of n α (t DD ) = 1

and n β (t DD ) = 0 and if the time elapsed from the beginning of phase 2 of calculations is denoted τ (τ ≡ t -t DD ), one gets from Eqs. 21:

n (R) β (τ ) = (1 -e -Γτ ) Γ (R) αβ Γ , n (NR) β (τ ) = (1 -e -Γτ ) Γ (NR) αβ Γ , n α (τ ) = e -Γτ , (22)
where Γ = β Γ (R)

αβ + γ Γ (NR)
αγ . Using solutions given by Eq. 22, occupation numbers are obtained for each trajectory and for each final electronic state φ β resulting, at time t = t DD + τ , from the φ α → φ β transition. Then, possible decay channels are detected for a particular transition and the abundance of fragments is updated using the occupation numbers provided by Eq. 22.

Ionic rare-gas cluster cations

The methodology described in preceding subsections has mostly been used, in our group, to study post-excitation dynamics of a specific class of model systems, namely the singly charged clusters of rare gases (Rg + N , Rg = He -Xe). The main reason for this choice consists mainly in a couple of favorable features of these species. Firstly, there are computationally cheap and still highly accurate interaction models for them which provide both the electronic ground state as well as excited state potential energy surfaces and couplings between them. Secondly, the models are based on (approximately) diabatic basis sets, so, equations of Subsec. 2.1 directly apply. Finally, a plethora of experimental data available on the rare-gas cluster cations in the literature allows for straightforward tests of the proposed methodology.

In this subsection we provide a compact survey of interaction models (model electronic Hamiltonians) used in preceding studies reported on the rare-gas cluster cations from our group. Since most of the material has already been published elsewhere, only a brief account is provided here for reader's convenience with frequent references to relevant preceding works.

Ne -Xe

For rare-gases with last filled shell being ns 2 np 6 (Ne-Kr), the intra-cluster interactions are described in the Rg + N ions within specifically designed diatomicsin-molecules (DIM) models. The DIM approach, originally developed by Ellison [START_REF] Ellison | [END_REF] and later applied to singly ionized clusters of rare gases by Kuntz and Valldorf [1], consists in, firstly, expanding the electronic Hamiltonian into a sum of diatomic and atomic contributions,

Ĥ = N -1 I=1 N J=I+1 ĤIJ -(N -2) n I=1 ĤI , (23)
where N denotes the number of atoms, and, secondly, in designing a specific valence-bond basis set of wave functions for which the corresponding Hamiltonian matrix can be calculated from the electronic energies of atomic and diatomic fragments. If the spin-orbit coupling is neglected, the atomic contributions are constant and their sum can be identified with the zero of the electronic energy, i.e., the atomic part of Eq. 23 is neglected. 13 The diatomic energies are provided from independent sources, either from experiments or (ab initio) calculations. References to studies from which particular diatomic curves have been taken for specific calculations are given in Sec. The original DIM approach as applied to rare-gas cluster cations can be extended in several directions which are particularly important for heavy rare gases. First of all, the spin-orbit coupling can easily be included [START_REF] Amarouche | [END_REF] via a simple semi-empirical atoms-in-molecules scheme [44], secondly, the most important three-body polarization forces (of the induced dipole-induced dipole type) can also be taken into account [START_REF] Amarouche | [END_REF] as well as three-body dispersion forces [45] acting among neutral atoms in the Rg + N complex. If the spin-orbit coupling is considered, the cardinality of the basis set (and the rank of the electronic Hamiltonian

He

For helium, the situation is a bit more involved since the usual DIM model based on a valence bond basis set 14 does not work properly and does even not reproduce the equilibrium linear symmetric structure of the He + 3 cation [46]. It is mainly due to a lack of interaction anisotropy, important in ionic complexes of rare gases, which cannot be properly modeled by spherical atomic orbitals of the valence shell of the helium atom (1s). To include the needed anisotropy, threebody contributions have been included in the He + N DIM Hamiltonian matrix by Knowles et al. [47,48]. In our calculations, we have used the extended DIM model as implemented for He + N and reviewed in detail in Ref. 49.

Results and discussions

Photodissociation

The first step in the photodissociation is the photoabsorption where the system is excited with photons of a well defined energy to some excited electronic state. In addition to the potential energy surfaces given by the DIM method, it is necessary to evaluate the transitions dipole moments which are necessary to calculate the transition intensity. After some tentatives to derive it from diatomic data, the use of the R vector according to the point-charge model [50] has been widely adopted. Following semiclassical ideas and assuming purely vertical transitions, an histogram is built and according to the potential energy difference the transition probability is accumulated from some appropriate sampling of nuclear configurations in the ground state [26,27]. This sampling may employ classical trajectories calculated under microcanonical conditions (fixed energy) [18], but care should be taken with the ergodicity [26,27,51], or under canonical conditions (fixed temperature) [43], often using Monte Carlo procedures, or from importance sampling of quantum wave functions either coming from a normal mode approach [52], a full quantum calculation [53], or from a quantum Monte Carlo evaluation [54]. In Fig. 1, typical theoretical absorption spectra of the ionic rare-gas trimers (argon, krypton, and xenon) calculated for visible (VIS) wavelengths [43,55,56] are compared with experimentally recorded ones [57,58,59]. Clearly the simulation works well for this initial step, giving confidence in the underlying DIM models and other theoretical approaches and assumptions. Clearly, the photon excitation is properly taken into account in the simulation for the VIS peaks for all the three rare gases as well as for the peak experimentally recorded [60] at ultra-violet (UV) wavelengths for argon.

The UV peaks for krypton and xenon were also predicted [55,56].

In typical photodissociation experiments, the system fragments after photon absorption and various fragments are collected, often by Time-Of-Flight (TOF) techniques, ideally all the fragments, however, mainly only the charged fragments since neutral fragments are much more difficult to detect. Therefore, a detailed information is experimentally obtained for the total Kinetic Energy Released (KER) and the Kinetic Energy Distribution (KED) of the fragments.

From the theoretical point of view, since the system starts, after photon absorption, from an electronically excited state, non-adiabatic dynamics is needed as non-adiabatic transitions to lower electronic states are possible and they may lead to internal conversion. How effective non-radiative transitions and internal conversion are is an important problem in photodissociation dynamics. In addition, the SO coupling also plays an important role as was clearly shown by experimental KER results [61] obtained for the heavier rare-gas trimers excited in the VIS region. In Fig. 2, such KER results measured for Xe + 3 are compared with data resulting from our simulations [62]. The agreement is spectacular and evidences that when there is enough energy, Xe + atomic fragments are formed in the excited fine-structure state and not in the ground state. Moreover, the switching is quite abrupt. It is noteworthy that calculations were even able to estimate the internal energy of the initial parent cluster (in the experiment) since our results obtained for hot clusters agree much better with experimental points than the results calculated for the cold ones. It is further clear that the SO coupling and internal conversion compete in the non-adiabatic dynamics of the ionic rare gases clusters.

Analysis of the KED also gives important information on the fragmentation dynamics. For the Ar + 3 photodissociation, a trimodal distribution was observed [59,61,63,64,65] and our simulations clearly also lead to fast and slow fragments as illustrated in Fig. 3. The relative abundance of the two peaks depends on the charge carried by the fragment. For the ions (Ar + ), fast fragments dominate while for the neutrals the reverse holds, slow fragments are more abundant.

However, the peak positions are more or less the same for both, illustrating that charged and neutral fragments experience similar dynamics. In fact, these reverse abundances reflect the initial charge repartition in the excited state. After VIS photon absorption (E phot = 2.35 eV in Fig. 3), for example, the charge is shared by the two outer atoms and the system mainly explodes leaving the central atom slow and the two outer ones fast. Later in the dissociation, non-adiabatic dynamics leads to some charge relocalization to the central atom, thus contributing to the slow peak for the ions. After UV excitations (E phot = 4.3 eV in Fig. 3), although the charge is initially spread over the trimer following a 0.25-0.50-0.25 pattern in the equilibrium centrosymmetric geometry, the situation is more involved than in visible excitations because asymmetric distortions are now favored by the absorption excitation process.

The ratio of the fast and slow peaks areas of the ions was called the symmetric branching ratio [61]. In Fig. 4, experimental [61] and theoretical [28] results are compared for Ar + 3 and Xe + 3 photodissociation at photon energies from about E phot ≈ 2 eV to E phot ≈ 5 eV. The agreement between the experiment and theory is very good, although it is rather difficult to converge the KED in the simulation and a huge number of dissociation events needs to be considered.

Here again, we could estimate the effective temperature of the initial clusters, which seems to be much larger for xenon than for argon.

The average size of ionic fragments emitted from photoexcited parents of larger sizes is illustrated in Fig. 5 for Ar + N (N = 6 -19) for VIS photodissociation in the main absorption peak (E phot = 2.35 eV) [33] and compared to the experimental results [66,67]. These more global data show the competition between evaporation, direct fragmentation of the core, cage effects, internal conversion, and other complex mechanisms emerging as the parent cluster size increases. More specifically, the photon energy is absorbed mainly by the core, therefore dissociative processes mainly start from the core. However, since the core is gradually surrounded by (almost) neutral atoms as the cluster size increases, collisions with them become frequent and resulting cage effects may lead to an efficient transfer of energy or even charge. However, if this is not the case, the surrounding neutral atoms are less bonded to the rest of the cluster than charged atoms and may easily detach via evaporation process. As each evaporation event yields the cluster with less internal energy, some kinetic energy being taken away by the departing atom, the cluster gets cooler and, when repeated, the evaporative cooling takes place.

Two theoretical approaches are also compared to each other, MFQ-AMP/S and MFQ-TFS/C , and lead to quite similar results. For larger sizes, it is however possible to notice some global increase of the non-adiabatic transition efficiency with the MFQ-TFS approach compared to the MFQ-AMP one. The MFQ approach is mandatory here since sequential fragmentation processes are obviously present as illustrated in the inset.

The agreement with the experimental results is very good, given again confidence in the theoretical simulation even for the highly involved nonadiabatic dynamics in large ions. Below N ≈ 11 -10, mainly dimer fragments are formed as well as some trimers. For larger parent sizes, on the other hand, the final fragment sizes increase with the initial parent size. For the smaller parent clusters, evaporative cooling is not sufficient to avoid the core breaking, but the production of ionic dimer fragments (instead of monomer ones) shows that internal conversion is quite effective.

Post-ionization fragmentation

The post-ionization fragmentation of rare-gas clusters has for long been a matter of controversy between theory and experiment, in particular for heavier rare gases (krypton and xenon). While, for these two heavy rare gases, experiment predicted an absolute dominance of monomer ionic fragments resulting from ionized parents decay [START_REF] Steinbach | Electron impact fragmentation of size-selected krypton clusters[END_REF][START_REF] Poterya | Xe clusters[END_REF], in theory [70,[START_REF] Poterya | Xe clusters[END_REF] ionic dimers dominated for small parent sizes or even larger ionic fragments as the parent size increased. For lighter rare-gases (argon and neon), on the other hand, the agreement between theory and experiment was rather satisfactory [39,71].

The above mentioned controversy was the main motivation for us to enter the game and to test our theoretical approaches in investigating this intriguing problem. To verify the reliability of the methodology, we started our calculations with the post-ionization fragmentation of ionic trimers of argon (Ar + 3 ) for which previous calculations [72] had been able to reproduce available experimental data [73] fairly well. Our calculations using the MF dynamical approach [32] and various interaction models predicted between 24% and 28% of monomer fragments, in a good agreement with both the experimental value (about 30%) and predictions resulting from preceding calculations (see the left-most and right-most columns of the Ar + 3 panel of Fig. 6). However, after we applied the MF method to the two heavier rare gases [19], a catastrophic disagreement emerged between our calculations and experiments [19]. More specifically, ionic dimer fragments dominated over monomers in our calculations, in contrast to the experiments where almost exclusively ionic monomers were observed, but in line, at the same time, with other theoretical calculations [70,[START_REF] Poterya | Xe clusters[END_REF] as well reporting the dominance of compound, dimer ionic fragments. In Ref. [19], we proposed this to be a consequence of the fact that quantum decoherence was completely ignored in the MF method used in our group and was only partly included in other calculations utilizing Tully's surface hopping approach [11]. After quantum decoherence was considered [19] in our calculations (MFQ methods of Subsec. 2.2), the agreement with the experiments was considerably improved for krypton and xenon while, at the same time, an excellent correspondence was preserved for argon (see Fig. 6). The agreement became, for the two heavier rare gases, fully qualitative (more ionic monomers were predicted from calculations than dimers) or even quantitative with an only slight underestimation of the experimental values by the theory.

The remaining quantitative difference between experiment and theory seen for krypton and xenon was finally interpreted in terms of the initial electronic excitation of the parent ion prior to its decay. Since nothing firm is known about this excitation from measurements, various models were adopted in calculations.

For example, an equally weighted linear combination of all the adiabatic states available within the DIM method was used in the surface hopping calculations [39,71,70,[START_REF] Poterya | Xe clusters[END_REF]. In our work, on the other hand, the initial electronic state of the fragmenting parent was originally set [32,19] in such a way that the positive charge was localized on a randomly selected atom occupying an electronic state prepared as a random linear combination of six atomic states included, for each atom in the cluster, in the DIM basis, Φ Kpmsz (see Subsec. 2.5.1), with K fixed to the selected atom and m = x, y, z and s z = ±1/2. Since, in addition, we found from the energy considerations that a (four-dimensional) subspace of this atomic state space corresponding to the P 3/2 term cannot contribute to the production of ionic monomer fragments due to lack of energy, only the (twodimensional) orthogonal complement of the P 3/2 space (corresponding to P 1/2 ) was considered in later calculations [21,31].

However, a quantitative agreement could still not be achieved. Mainly because a considerable amount of metastable, long-lived intermediate fragments resulted from our calculations even though quite long simulation times were used (hundreds of picoseconds). For this reason, we developed an approximate simulation scheme [21] allowing us to extend dynamical calculations to sufficiently long times so that all the intermediate metastable species could undergo expected fragmentation. Since basically all time scales became accessible through this scheme, we called it a multiscale method (see also Subsec. 2.4). After the multiscale method was applied to the post-ionization fragmentation of ionized krypton clusters [31], a close-to-quantitative agreement was achieved with available experimental data. For a summarizing view, see Fig. 7 from which both a good performance of the multiscale method as well as shortcomings of direct dynamical calculations if performed at insufficiently long times are clear.

To further verify the assumption that long times are needed to get converged dynamical simulations of the post-ionization fragmentation of cluster cations of heavy rare gases, a series of long-time calculations were performed [74] for the Kr + 7 ion (the largest parent size for which experimental data are available [START_REF] Steinbach | Electron impact fragmentation of size-selected krypton clusters[END_REF]) using an improved simulation protocol with periodical checks of the fragmentation state of the decaying ionic complex and with detected intermediate fragments treated separately (see Ref. [33] for more details). In principle, this extension can be combined with any dynamical approach described in Subsec.

2.2 and, if this is the case, a suffix "(sf)" (an abbreviation of "separate fragments") is added to the acronyms introduced in Subsec. 2.2. As an illustration of such a long-time calculation, the post-ionization fragmentation of the Kr + 7 ion starting from the highest adiabatic electronic state available within the DIM approach is exemplified in Fig. 8. It is clearly seen that not only the experimental data are closely reproduced within such an treatment, but also a detailed information can be obtained on the processes responsible for the disintegration of the ionic parent to the final, elemental fragments (Kr + and Kr + 2 ) through a series of intermediate, metastable species.

Both the approaches used for direct dynamical simulations as well as the multiscale method could easily be extended to include radiative transitions (see Subsec. 2.3 for details). This allows, among others, to collect the information about the luminescent radiation emitted from metastable species along their decay. As an example, we show, in Fig. 9, calculated luminescent spectra of metastable krypton and xenon dimer cations. Noteworthy, such an additional information resulting from calculations is important for further experimental validation of adopted theoretical approaches. Unfortunately, no related experimental data are presently, to our best knowledge, available in literature.

Collisions

To complete this review about the non-adiabatic dynamics processes modeled using the methods developed by our team, let us discuss two kinds of collisions. The first one is a rather usual collision problem where we considered, at various collision energies, collisions of Xe + ions with neutral xenon trimers (Xe 3 ). In this up to now unpublished study, we focus on the production of ionic xenon dimers (Xe + 2 ), having in mind the ionic clusters growth in supersonic expansions. In order to form a dimer ion, the charge needs to hop from the initial charged monomer to a larger cluster while the subsequent detachment of neutral atoms can further stabilize the final ionic dimer. As can be seen in Fig. 10, at low collision energies (represented by E coll = 0.02 eV in this work), the efficiency of the production of Xe + 2 is much higher if one starts with the monomer ion (Xe + ) in its fine-structure ground state (P 3/2 ), while at intermediate (E coll = 0.1 eV) and higher (E coll = 0.5 eV) collision energies, the initial fine-structure excitation plays a rather lesser role. Moreover, the global efficiency strongly decreases as the collision energy grows up, regardless of the initial monomer ion excitation. These new results clearly show that, to increase the ionic cluster size, collisions at low energies should be privileged. Note also that care has been taken here to use a method which preserves total angular momentum. Particularly because the total angular momentum is closely related to the initial impact parameter, and in a method where it is not conserved, data corresponding to different impact parameters would be nonphysically mixed and this would finally lead to wrong results.

The second illustration was primarily motivated, as already mentioned in Sec. 1, by cold plasmas production and their use in various applications, e.g., medical applications for He [START_REF] Laroussi | [END_REF]4], and UV sources [5] and spacecraft propulsion for heavier rare gases [6]. The cold plasmas are usually produced in a corona discharge apparatus, but the need to make the apparatus more practical and easier to build and/or use requires a detailed knowledge of main microscopic processes governing the properties of these cold plasmas. Among others, Momentum Transfer (MT) cross-sections 15 and integral cross-sections were calculated for respectively Non-Reactive Scattering (NRS) and Collision Induced Dissociation (CID) of various rare-gas dimer ions colliding with corresponding rare-gas atoms. The ions were considered both in their electronic ground state as well as in excited states, and also various vibrational and rotational states of the dimers were investigated. Selected results taken from previous studies carried out in our group are shown in Fig. 11. In these studies, the method used for the non-adiabatic dynamics treatment was the MFQ-AMP/SR one (as a computationally less demanding alternative) since basically tens of thousands of trajectories had to be computed for each collision energy to get a reasonable level of convergence. As is clearly seen in Fig. 11, the behavior of the cross-sections with increasing collision energy is rather similar for all the rare gases considered. The MT cross-sections decrease smoothly from low collision energies to about E coll = 2 eV and then the decrease is more pronounced. At collision energies of E coll ≈ 1.5 eV, i.e., slightly below the energy at which the slope of the MT cross-section curves changes, CID processes start, become effective, and rapidly achieve their saturation. In such a way, the fact that the two processes, NRS and CID, combine and make disappear the initial ionic dimer leads to a global smooth decrease of the MT cross-section from the small collision energies (E coll = 0.01 eV) to the high ones (E coll = 100 eV). Noteworthy, the calculated cross-sections were used in a series of realistic calculations of transport properties of the considered ions [75,76,77,78,79] the results of which mostly very well reproduced available experimental data, extend them to experimentally unreachable conditions, or even provide data which were not recorded experimentally but are, a t the same time, important for macroscopic plasma modelings.

Conclusions

Starting from a pioneering "on-the-fly" method [14,16,22,23,24] and challenged by interesting experimental results, we have improved our hybrid dynamical approach (quantum treatment of electrons and classical description of nuclei) [32] in two directions consisting in: a) the inclusion of periodic decoherence [19] to face sequential fragmentation and fix the fractionary charge problem, and b) a multiscale approach [20] to include long-time evolution involving radiative and non-radiative processes. These developments allowed us to treat complex non-adiabatic dynamics occurring in photodissociation, postionization fragmentation, and collisions of ionic rare-gas clusters. In general, the larger the system, the more efficient the non-adiabatic relaxation processes are and, therefore, internal conversion increases with the size of the cluster and more generally of the system. However, the spin-orbit coupling introduces two groups of states in the rare-gas cluster ions, and inter-group or intra-group nonadiabatic transitions act rather differently [62]. Since the inter-group transitions are rather improbable, a concept of pseudo-ground state for the excited state group could be introduced. Some comparisons made with full quantum results (in the context of cold plasmas modeling) [75,80] are very encouraging, although some differences remain related to the classical treatment of the nuclear motion for helium [75] and involved effects for the heavier rare gases [80].

The methodology developed up to now and described in this paper seems to have a strong potential for further extensions and applications. For example, it could be very interesting to test the methods on other systems and problems, like vibrational predissociation, where many experimental results can be found to test the developed theoretical approaches, and also quantum treatments of the smaller systems. Even for rare gases, there is still much to be investigated, e.g., larger ionic clusters, mixed ionic clusters, atoms and molecules solvated in rare-gas envelopes etc. All the methods we have developed to carry on realistic simulations could be involved to treat a very general problem of non-adiabtic dynamics, namely how the excess energy flows and is dissipated or relaxed in complex molecular systems with many electronic states involved and frequent transitions between them.

for Research, Experimental Development and Innovations project (grant no. LM2015070). Partial financial support from the RINGEN programme (grants no. LM2015084 and CZ.02.1.01/0.0/0.0/16013/0001792) is also acknowledged (IJ). 3 cations (Rg = Ar, Kr, and Xe) in the VIS region. Discrete symbols represent experimental data (circles -Ar [57], squares -Kr [58], triangles -Xe [59]), lines correspond to theoretical results [43,55,56] obtained for the DIM+SO model and based on classical Monte Carlo samplings of cluster configurations performed at effective temperatures best reproducing the experiments (full line -Ar, T eff = 150 K, dashed line -Kr, T eff = 200 K, dotted line -Xe, T eff = 250 K). In the inset, experimental [60] and theoretical [28] data are depicted for the UV spectrum of Ar + 3 . Note that in this case T eff = 230 K has been used in theoretical calculations. [66] and lower triangles [67]); stepwise linescalculations [33] performed for the DIM+SO+ID-ID(dmp) interaction model using the MFQ-AMP/S(sf) (solid line) and/or MFQ-TFS/C (dotted line) dynamical approaches. In the inset, time-dependent abundances of (intermediate and final) ionic fragments resulting from Ar + 9 are depicted as obtained for the MFQ-AMP/S(sf) approach. : Luminescence radiation resulting from the metastable decay of the Kr + 2 (solid curves) and Xe + 2 (dashed curves) ions calculated for rotationally-vibrationally ground-state dimers ([j=0,v=0]) using the multiscale approach of Subsec. 2.4 and interaction potentials taken respectively from Refs. 82 and [83]. Input configurations to the long-time phase of the multiscale calculations have been generated, for both electronic states, from squares of respective [j=0,v=0] nuclear wave functions and spectra have been collected at t → +∞. Momentum-transfer (MT) cross-sections (full symbols) and collision-induced dissociation (CID) cross-sections (open symbols) calculated for Rg + 2 (j=0,v=0)/Rg collisions using the MFQ-AMP/SR dynamical method and the DIM+SO interaction model, Rg = He (squares) [75], Ne (circles) [77], Ar (up triangles) [76], Kr (down triangles) [78], and Xe (diamonds) [79]. MT cross-sections

Tables and figures

  and (at least approximately) diabatic, Φ α | ∂Φ β ∂qi = 0. This results, after inserting the expanded forms of Eqs. 3 and 4 into Eq. 1, into 3

3 .

 3 The basis set allowing this consists of 3N valence-bond Slater determinants,|Φ α ≡ |Rg...RgRg + Rg...Rg| (α = [K, p m ], K = 1, • • • , N ,and m = x, y, z), where Rg and/or Rg + denote respectively electronic configurations of an Rg atom and/or ion, the latter corresponding to the K-th nucleus and having an electron removed from the p m -orbital. A particular wave function of this basis set represents thus a state with the positive charge localized in p m -orbital of atom K. The corresponding Hamiltonian matrix, H αβ ≡ Φ α | Ĥ|Φ β , is then obtained from diatomic inputs as described in Ref. 1 and reviewed in detail in Ref. 43.

  matrix) doubles, α = [K, p m , s z ] with s z = ±1/2 being the spin projection of the electron removed from the valence shell. The inclusion of the three-body forces, either polarization forces or dispersion forces, does not, on the other hand, change the dimension of the electronic Hamiltonian matrix and requires only a modification of its diagonal. For more details see Ref. 9 and Ref. 45. Independent inputs needed to this extended DIM approaches consist of a) spinorbit coupling constants of atomic ions ( Rg + ), rare-gas atom polarizabilities, and three-body dispersion coefficients (for specific settings see Sec. 3). Several interaction models result from the diatomic-in-molecules approach and its extensions: a) the original model of Ref. 1 (denoted here as the DIM model), b) a model augmented with the spin-orbit coupling (denoted as the DIM+SO model), c) a model augmented with both the spin-orbit coupling and the induced dipole-induced dipole interaction terms (the DIM+SO+ID-ID model), and d) a most comprehensive model including, in addition to the previous extensions, the three-body dispersion contributions (the DIM+SO+ID-ID+N3 model). Basically, all of the four models have been used in our calculations to assess the role of various interaction contributions.

14

  For helium cluster cations, He + N , a basis set is used in the form of |Φ k = |He...HeHe + He...He| (with He and He + denoting respectively 1s 2 and/or 1s 1 configurations of the helium atom and/or ion, and He + occurs in the K-th position, K = 1, • • • , N ). Each particular wave function represents thus a state of the cluster with the positive charge localized in 1s-orbital of atom K. As a consequence, totally N basis set wave functions are used for He + N and the spin-orbit coupling is neglected since non-zero angular momentum states are not considered.

Figure 1 :

 1 Figure 1: Photoabsorption spectra of Rg + 3 cations (Rg = Ar, Kr, and Xe) in the VIS region. Discrete symbols represent experimental data (circles -Ar [57], squares -Kr [58], triangles -

34 Figure 2 : 3 ion (Xe + 3 →

 34233 photon energy[eV]

2 )Figure 3 :Figure 4 :

 234 photon energy[eV]

  photon energy[eV]

Figure 5 :

 5 Figure 5: Average size of ionic fragments for Ar +N parents photoexcited by E phot = 2.35 eV: symbols -experimental data (upper triangles[66] and lower triangles[67]); stepwise lines -

2 Figure 9

 29 Figure9: Luminescence radiation resulting from the metastable decay of the Kr + 2 (solid curves) and Xe + 2 (dashed curves) ions calculated for rotationally-vibrationally ground-state dimers ([j=0,v=0]) using the multiscale approach of Subsec. 2.4 and interaction potentials

2 Figure 10 :

 210 photon energy[eV]

Figure 11 :

 11 Figure 11:

Whereas needen, capital Latin indices will be used to label particular atoms.

Throughout this paper, small Greek indices always label electronic basis set wave functions and/or electronic states.

For a more general derivation see reference 32.

Since the electronic Hamiltonian matrix is hermitian, its real part and imaginary part must be symmetric and/or antisymmetric, respectively.

The quenches of the time-dependent electronic wave function towards an adiabatic electronic state can be understood as instantaneous collapses of the wave function upon a classical measurement performed on quantum electrons by atomic nuclei. As a consequence, they model the process of dissipation of quantum information from the electronic subsystem to the classical environment.

Note that the C method leads to a slightly unphysical behavior since it does not fully

Both parameters are calculated at the current nuclear configuration reached on the classical nuclear trajectory

.11 Radiative transitions have always been attempted immediately after a successful quenching event, the adiabaticity of the starting state, φα, is thus secured.

Note that the occupation numbers are calculated for each particular trajectory entering phase 2 of the multiscale treatment and the results obtained for all the trajectories are then used to update fragment abundances as described below.

In this case, the zero of the energy of an Rg + N cluster corresponds to a fully dissociate state, Rg + + (N -1)Rg.

The MT cross-section is defined as follows: σ MT = 4π σ d (1 -cos φ) dΩ, where σ d is the usual differential cross-section, φ denotes the deflection angle, and Ω is solid angle.

Acknowledgments

Traveling support from the Barrande programme of French-Czech collaboration (grant. no. 40688TH and 8J18FR031) is warmly acknowledged. The calcu-

 [32]and MFQ [19] calculations using the DIM+SO model at time t = 100 ps after the ionization event (when the dynamical calculations are well converged [19]) and compared with available experimental data reported in references [81] (Ar + 3 ), [START_REF] Steinbach | Electron impact fragmentation of size-selected krypton clusters[END_REF] (Kr + 3 ), and [START_REF] Poterya | Xe clusters[END_REF] (Xe + 3 ). For comparison, results of independent Tully's Surface Hopping (TSH) calculations are also shown [71,70,[START_REF] Poterya | Xe clusters[END_REF]. fragments [31] resulting from the post-ionization fragmentation of parent Kr + N ions as obtained from the multiscale method used to extend dynamical data calculated, up to t = 200 ps after the ionization event, via the MFQ-AMP/S approach. For comparison, the MFQ-AMP/S dynamical data recorded at t = 200 ps (dashed and dotted lines are used for monomer and dimer ionic fragments, respectively) and experimental data [START_REF] Steinbach | Electron impact fragmentation of size-selected krypton clusters[END_REF] (full symbols) are also shown.

In the theory, the DIM+SO+ID-ID(dmp) model has been used for small parents (N ≤ 5) and the DIM+SO+ID-ID(dmp)+N3(DDD) model have been applied to larger parents (N ≥ 6). resulting from the post-ionization fragmentation of the Kr + 7 ion initially excited to the highest electronic level (L21) as obtained [74] from the MFQ-TFS/S(sf) dynamical method and the DIM+SO+ID-ID(dmp)+N3(DDD) interaction model. Horizontal lines correspond to respective experimental values of monomer and dimer abundances as reported in reference [START_REF] Steinbach | Electron impact fragmentation of size-selected krypton clusters[END_REF]. Sequential fragmentation of electronically highly excited Kr + 7 ion upon electronimpact ionization of its neutral precursor (Kr 7 ).