Jacqueline Bergès

Dominik Domin

Julien Pilmé
email: pilme@lct.jussieu.fr

Benoît Braïda

Chantal

Chantal Houée-Levin

•OH oxidation of methionine in the presence of discrete water molecules: DFT, QTAIM and valence bond analyses

Keywords: Oxidation of Methionine by OH radicals, explicit water molecules, DFT calculation, QTAIM, VB analysis, two center-three electron bonded radical

Electronic supplementary material. This article contains supplementary material.

Introduction

Sulfur containing molecules are major targets of reactive oxygen species (ROS) in biological media. A wealth of experimental results has shown that •OH radicals oxidize quickly and efficiently these biomolecules. The oxidation of these compounds has numerous consequences in vivo. The primary mechanism begins by the formation of an >S-• OH adduct which would dissociate into a sulfur centered cation radical and a hydroxide molecule (Reactions (1) and ( 2)).

(1) >S +  OH  >S- OH (2) >S- OH  >S + + OH - The oxidation mechanism of the methionine (Met) amino acid or residue in peptides by • OH radicals has been largely studied experimentally from the very first steps (reaction ( 1)) [START_REF] Asmus | The Study of Fast Processes and Transient Species by Electron Pulse[END_REF][START_REF] Asmus | Radiation Chemsitry: present status and future trends[END_REF][START_REF] Bobrowski | [END_REF][4][5] to the final compounds [6][7][8]. It is well known that Reaction (1) is dominantly diffusioncontrolled for the amino acid as well as for Met residues in peptides and proteins [9] and that it is followed by several rearrangements triggered by Reaction [START_REF] Asmus | Radiation Chemsitry: present status and future trends[END_REF]. Due to the complexity of the following steps, this wealth of experimental data had to be completed by numerous computational works [5,8,[10][11][12][13][14][15][16][17]. These studies were mainly focused on deciphering the nature of transient species, their stability and their eventual reactivity. Indeed, it has long been suggested [START_REF] Asmus | Radiation Chemsitry: present status and future trends[END_REF] that the sulfur centered radical cation may undergo some kind of cyclization by forming a two-center three electron bond (2c-3e) with a heteroatom that has a lone pair (see [4][5] and references therein for reviews on the matter). Recent Density Functional Theory (DFT) calculations (using the BH&HLYP functional) and Quantum Chemical Topology analysis confirmed the stability of cyclic species containing 5 to 9 atoms in dipeptides [13][14]. The eventual reactivity of the methionine radical cations has also been studied in the view of predicting possible intra-or intermolecular electron transfers in proteins. It has been shown that the redox potential of methionine radicals in peptides varies with the nature of the vicinal residue, the geometry and the size of the pseudo cyclic group [11][12]. This finding provides the basis for understand, for instance, why in some cases long range intramolecular electron transfer occurs from tyrosine to the methionine cation radical. In Met enkephalin, the electron transfer pathway has been computed [18]. Finally, many authors have considered the oxidation reaction in the presence of oxygen [19].

In these simulations, the main difficulty is with the description of the solvent and in the choice of the terminal functional groups of the amino acid in order to provide a more realistic view of protein oxidation and avoiding charge. The solvent was frequently modeled by a continuum (CPCM or IEFPCM [13,19], or by a dielectric constant [20]).This crude description lacks explicit molecular solvation layers that may participate in the oxidation mechanisms. Chu and co-workers added some explicit water molecules in the study of a related reaction, the oxidation of dimethylsulfide (DMS) and Met by H 2 O 2 [21].

Moreover, for only Met embedded in a continuum, several terminal functional groups were tested : the simple zwitterion [22] -CO 2 CH 3 and -NH 3 + end groups, neutral Met with -CO 2 H and -NH 2 terminal functional groups in the gas phase [19], two peptide groups -CO 2 NH 2 and -NH 2 CO 2 [20].

In this work, we reinvestigated the electronic factors controlling both reactions (1) and ( 2 Methionine was methylated to avoid the influence of the terminal groups in the zwitterions (electrostatic interactions between -NH 3 + /-CO 2 -and the water molecules) i.e. the end groups were -CO 2 CH 3 and -NH 2 CH 3. To test the importance of the terminal functional groups, we replaced the amine methyl group by the amide group (mimicking a peptide bond) in some calculations.

Whenever possible we described the cations resulting from OH -elimination with the same computational methodology and we drew the energetic landscape of their various conformations. Since the nature of bonding in the different compounds is of interest, a

Valence Bond (VB) type of analysis can be revealing (for review see [START_REF] Shaik | DMP (ed) Chemical Bond II: 100 Years Old and Getting Stronger[END_REF] and references within). Previous VB calculations have indeed shown that the interaction between a • OH radical and one of the sulfur lone pair electrons of a methionine derivatives, may lead to either a weak electrostatic adduct, or to a stronger two-center three-electron bond, depending on the chemical groups linked to the sulfur [14,23].

Computational methods

In the first step, MD simulation with geometries of the solute molecules kept fixed, were performed in an environment containing a very large number of explicit water molecules, i.e. more than 2000 water molecules in each periodic boundary condition unit cell. The CHARMM CGenFF 3.0.1 force field was used for these solute molecules and the CHARMM modified rigid TIP3P water model [START_REF] Danovich | [END_REF][26][27] was used for the solvent. The molecules were solvated by water in 20 Å in each direction. For the charged species, the force field parameters were obtained for the analogous neutral methionine and the net formal charge was placed on the sulfur atom. Subsequently from these equilibrated MD trajectories randomly selected snapshots where used for extracting nanodroplets containing ~70 solvent water molecules for further geometry reoptimization using the BH&HLYP functional, 6-31+G** basis set and a PCM explicit solvation model to account for the rest of the bulk water. From the latter optimized microsolvated structures, only the few water molecules in the first solvation shell that were directly in interaction with the Met-OH adduct were retained and used in subsequent highest level quantum chemistry calculations. Further geometry optimizations were carried out with BH&HLYP/6-311+g(d,p) on the clusters of the solutes surrounded by explicit water molecules selected in the previous MD steps. Note that it was previously shown [14] that optimizations using the DFT functional BH&HLYP, seemed specially adapted to describing radicals that contain 2c-3e bonds. The clusters were embedded in an implicit solvent (CPCM) that modeled the bulk aqueous environment. All the purely QM (DFT) calculations were carried out using the Gaussian 09 suite of programs [28].

Quantum Theory of Atoms-In-Molecules (QTAIM). The aim of this approach is to answer questions about the chemical bonding in molecules and solids, and predict chemical reactivity trends. The methodology relies on the theory of gradient dynamical systems that enables a partitioning of the 3D molecular space into topological basins via the properties of the gradient vector field of the total electron density [29]. As its use has been widely reviewed, the formalism is only briefly recalled here. The QTAIM atoms can be defined as the union of a nucleus with its associated electron density basin. The basins are delimited by zero-flux surfaces and the integration of the electron density over each basin directly provides its corresponding population. The QTAIM atomic charge, q(A), is calculated by subtracting the electron population of the topological atom A from its atomic number. The atomic spin density is calculated as the difference between the both alpha and beta populations.

Introduced by Bader and Stephens [30], and later recovered by Fradera et al. [31], the QTAIM delocalization index (DI) is a measure of the electron-sharing between two atoms and can be sometimes compared to other bond order indices.

Valence Bond Analysis. In a previous article [23], three of us have shown that, in the specific case of a (2c-3e) bond, the weights of the VB structures could be faithfully estimated from a spin density analysis obtained after a correlated MO-based calculation. The equations derived can be easily generalized in the present case by equations 1-3, where ω i are the estimated weight of structure i, and ρ A the spin density value on atom A. A proper description of a (2c-3e) bond has been shown to be potentially problematic both for DFT [32] and post-HF calculations [33,34]. However, levels of calculation such as BH&HLYP functional, or the Coupled-Cluster CCSD method, can typically provide estimated ω i that have been shown to be in good agreement with those computed using the high level BOVB Valence Bond method.

ω ρ ρ ρ ρ ω ρ ρ ρ ρ ω ρ ρ ρ ρ

Results and discussion on the adducts Met-OH

We focused on two aspects of the complexes: the position of the •OH radical relative to the molecule and the proximity of the sulfur atom to either nitrogen or to one of the oxygen atoms (MetSN or MetSO respectively). We selected four types of initial complexes coming directly from MD snapshots. These complexes are depicted on Scheme 1 where the  OH locations close to atoms S, N or O of Met have been surrounded by circles. The optimized geometries of these complexes, surrounded by a continuum CPCM, are depicted on Figure 1.

In the complexes MetSN-OH (I) and MetSO-OH (IV), the S-OH bond is typically close to 2.3-2.4 Å. In contrast, the complex II (II), •OH radical is close to both S and N atoms. The spin densities are found to be fairly distributed between S and N, with a S-N distance of 2.5 Å. In the complex III (III), the •OH radical is close to nitrogen only when S-N distances are short. Another complex was found with short S-N and S-• OH distances (2.4 Å) (not displayed on Figure 1). This complex was found instable and after the optimization process, the calculation led back to I. S1, column 1). In the alternative methodology, the solvent was described by a continuum modeled by CPCM (Table S1, column 2 & 3).

As observed in Table S1, both methods led to similar bond lengths and spin density distributions. In addition, the adducts MetSN- OH and MetSO- OH had very close electronic energies (E ≈ 2kJ/mol, Table S1). In both cases I and IV, spin densities and atomic charges were compared with three different schemes: Mulliken, QTAIM and NBO analysis. The three different models of atomic charges and spin densities led to similar conclusions. S and  OH always shared the spin densities, i.e. were linked by 2c-3e bond as observed for DMS-OH [23]. No SN or SO bond was formed. The partial charges remain mainly distributed on both the S (δ+) atom and OH (δ-) group. This is consistent with the spin density analysis and also with the delocalization index, which supports this statement since the sulfur interacts mainly with the oxygen atom of OH while the S-N interaction remains negligible. Thus, further calculations were performed using CPCM continuum model and further bonding analyses were preformed with the QTAIM method.

Continuum CPCM model only. Starting from points selected from MD where various positions of • OH around Met have been found, fully optimized calculations were performed using the BHHLYP/6-311+g(d,p) + CPCM level of theory. Table 1 gives the energetic and structural parameters together with QTAIM analysis for the I, II, III and IV complexes. Conversely, the S-N distance was initially short in III, while S and OH were far from each other. The optimized complex III corresponds to a hydrogen transfer that took place from the amino group to the  OH radical formation of a water molecule (see Figure 1), the spin density being mainly localized on the N center ( NH = 0.96 a.u). For the optimized complex IV, OH radical is close to both S and one of the oxygen atoms (2.39 Å) and the spin densities are found distributed between S and OH ( S = 0.32 a.u. and  OH = 0.66 a.u). As observed in Tables 1 and2, III displays a OH charge close to -1.0 e and a typical strong spin density close to 1.0 only located on the nitrogen center. These results are consistent with the formation of a single N-centered radical and a water molecule bonded to the N-center with a hydrogen bond (Figure 1). When the methyl group of the N-terminus was replaced by the amide group NHCOCH 3 the H transfer from -NH and formation of a water molecule was also observed but a conformational change led to MetSN + 2c-3e bond like in II (see Table 1) In the presence of these discrete 4 water molecules, both I and II (Figure 2) have now rather close to relative energies (ΔE less than 8 kJ/mol) (Table 3), contrary to what was obtained when the aqueous environment was only modeled by CPCM (Table 1). The values of QTAIM descriptors (Table 3), are rather similar to those provided in Table 1 for all three complexes I, II and III (MetSN-OH + 4 H 2 O).As for III + 4 H 2 O (Figure 2), the result is also an H transfer from the amino group to  OH, leading to a N-centered radical and H 2 O. Thus, the spin density  NH reaches a strong value of 0.97, fully located on the N-center. Note that the charge Q(S) goes to almost zero (-0.04 e). The  OH radical was either interacting with the amino group or moved thanks to hydrogen bonds with another water molecule. This complex remains the most stable configuration and the relative energy has decreased because of the presence of water molecules (ca. 100 kJ/mol).

MetSN-

Addition
Three snapshots of the optimization process leading to III were extracted. and delocalization indexes δ S-N, δ S-OH of S-N and S-OH bonds are given in a.u..

MetSN+4 H 2 O III

First point

During the H transfer process, the distances d S-N logically increased from 2.52 Å to 3.20 Å whereas simultaneously, the d S-OH decreased from 5.96 Å to 4.52 Å. Consequently, the delocalization index δ S-N strongly decreased from 0.48 to 0.06. In addition, the spin density  S fell to zero and  NH reached a maximal value close to 1 at the end of the H-transfer process.

This confirms that complex III exhibits a single radical center located on the nitrogen.

Actually, only II exhibits a SN 2c-3e bond. There is no 2c-3e SO bond. In the socalled MetSO-OH IV, the most stable complex belongs to family I whatever the conformations of the carboxylate group (see Table S2): In I 1 the sulfur is closer to the ester group while in I 2 , it is closer to the carbonyl group. Both exhibit common features with an electronic energy gap ΔE of less than 18 kJ/mol and a 2c-3e bond between S and  OH with distances of ca. 2.33Å (Table 3). These complexes are similar to the one obtained using only the CPCM model (Tables 1 and2; Figures 1 and2), the four water molecules forming a network of hydrogen bonds surrounding the adduct MetSO-OH.

To conclude, upon addition of four discrete water molecules the most important change is that now II has electronic energies close to that of I (E 7 kJ/mol instead of 80 kJ/mol with CPCM only). The presence of 4 water molecules stabilizes the MetSN •+-OH cluster. Nevertheless, the most stable complex is still III, which undergoes H transfer from -NH to • OH leading to formation of a N-centered radical. Only two complexes of the adducts are likely to contribute to the heterolytic dissociation of the adduct leading to a radical cation MetS + : II that could be considered as a precursor of MetSN + plus OH -anion and I that, like in DMS-OH, can have a good probability to undergo this rupture depending on the chemical environment.

Valence Bond Analysis. The three-electron bond in the sulfur-OH case is characterized, within a VB formalism, by the resonance between VB structures 1 and 2 (see Scheme 2), with comparable weights for both structures. These two VB structures differ from a shift of one electron, and thus by a shift of charge, between the two bonded centers. Besides, the stability of the (2c,3e) bond arises from the resonance energy due to the mixing of the two structures, thus its classification as a "complete" charge-shift bond [START_REF] Shaik | DMP (ed) Chemical Bond II: 100 Years Old and Getting Stronger[END_REF][START_REF] Danovich | [END_REF]. A weak electrostatic adduct would correspond to the case where one of the two structures would be largely dominant in the description, and the other one negligible (~10% weight). In the case of the system studied here the situation is slightly more complicated, though. In conformations of Met where the nitrogen of the methionine is close enough either to the • OH radical or to the sulfur, it could also be involved in a (2c-3e) bond with either of the two.

Therefore, a complete VB description involved the three VB structures displayed in scheme 2,

where the radical could be delocalized on three different centers. Now five different cases could in principle occur: i) a weak electrostatic adduct, characterized by a largely dominant VB structure 1 and low weights for both VB structures 2 and 3; ii) a sulfur-oxygen (2c-3e) bond, with both VB structures 1 and 2 significantly resonance contributors while VB structure 3 has a negligible weight; iii) a nitrogen-oxygen (2c-3e) bond, characterized by a 1⟷3 resonance and a low weight for structure 2; iv) a sulfur-nitrogen (2c-3e) bond, characterized by a 2⟷3 resonance and a low weight for VB structure 1 ; v) when the three VB structures have comparable weights, which would correspond to a radical delocalized between the three atoms.

In Table 5, the estimated VB structure weights for I and II using the QTAIM spin densities, with or without explicit quantum description of the four water molecules corresponding to the first solvation shell, are displayed. 5. Weights ω i (in %) of the three VB complexes computed at the BH&HLYP level of theory (see scheme 1) for different selected geometries, as inferred from equations 1-3, using the QTAIM spin densities.

From these values, and the previous analysis, the nature of the interactions in the different complexes can readily be inferred. Let's start with geometry I. With a PCM model only as the solvent (first line in Table 5), the VB structure 2 comes out as the major one, with structures 1 and 3 as minor but not negligible. The electronic structures change substantially when the four water molecules are added (second line in Table 5), with a marked increase of the weight of the VB structure 1, from 24% to 45%, and a collapse from 11% to 1% of the weight VB structure 3, which now becomes totally negligible. All in all, almost equal VB weights for 1 and 2 are obtained for geometry I in the presence of four water molecules, which enables us to unambiguously conclude to the existence of a (2c-3e) bond between S and OH in this geometry. This is unprecedented, because complete (2c-3e) bonds in overall neutral compounds in general, and in the S-OH case in particular, are very rare, the difference in charge between the two resonance contributors (1 and 2 here) usually leads to a dominant neutral SOH Valence Bond structure. Quite different is the case of II (third and fourth lines in Table 5). For the latter, the weight of VB structure 1 is zero, which is consistent with the OH being remote from the sulfur atom. The weights of VB structures 2 and 3 do not change much with the introduction of the four water molecules, and come out as rather balanced, indicative of a (2c,3e) bond between the sulfur and nitrogen atoms within the cationic MetSN + fragment, in electrostatic interaction with an OH -.

The oxidation mechanism of sulfur containing molecules has led to hot debates. Two dissociation channels are typically invoked, the "homolytic" pathway with a dissociation to

•OH and neutral fragment, meaning that no oxidation took place, and the "heterolytic" pathway with dissociation to OH -and a cationic counterpart. The "homolytic" channel is theoretically the most preferred one, because the separate neutral products are thermodynamically more stable, even in water environment. However, we have found here the electronic structure of complex I contains a substantial amount of the VB structure 2, corresponding to a MetSN + cation in interaction with a OH -; and that the complex almost isoenergetic to complex I in presence of explicit water molecules, fully corresponds to an interaction between MetSN + and OH -(mixing of VB structures 2 and 3). Therefore, in the case of a diabatic dissociation process, the heterolytic pathway appears as a fully plausible one for MetSN-OH in agreement with the wealth of experimental data (see introduction).

To conclude, all these first results shown that the formation of complexes I, II. and IV represents a key step leading to the heterolytic dissociation of Met- OH in Met + + -OH (radical cations formations).

Results and discussion on the radical cations Met •+

The radical cation either alone or the cluster radical -water molecules were embedded in a continuum. We have first compared both cases to evaluate the importance of the discrete water molecules, following the same methodology as for the adducts.

MetSN + and MetSO + + CPCM. First, we have performed various scan optimization of d SN and d SO i.e. one of these distances was frozen while the other one was gradually raised, the remaining of structural parameters were fully optimized at each step at the BHHLYP/6-311+g(d,p)+ CPCM level of theory.

Conclusions

The goal of this study was focused on the analysis of the electronic factors governing the first steps of the oxidation, formation of OH-adducts and possibilities of heterolytic ruptures. We compared several methods describing the solvent: CPCM only, molecular mechanics description of water molecules (QM/MM) taken from MD snapshots, CPCM with several water molecules coming also from MD snapshots. We found that a small number of water molecules was sufficient to get a more precise description (4 for neutral species with large basis set; 7-11 for charged species with smaller basis set).

The second problem we had to face was the endgroups. We chose to add methyl groups to get rid of the electrostatic interactions of the zwitterion. In the literature, the problem of end groups of Met was solved differently: by letting crude charged groups [22]; more recently by putting peptidic bonds [20].

We performed some calculations for the complex III where the NH 2 CH 3 endgroup was replaced by either -NH 3 + or amide. When the -NH 3 + endgroup was used, the trends were identical to those where -NH 2 CH 3 endgroup was used ,.i.e. H transfer between OH and NH 3 .

In contrast, there is reorganization of Met and formation of a (2c-3 e) SN bond (Table 2) if the amide endgroup is used.

Major new insights were obtained:

1) Several positions of OH radical vs. Met were found concomitant with the different positions of water molecules. The probability of OH radical to attack any atom is well known in pulse radiolysis although the rate constants are different. For instance OH could be close to -NH, far from S etc..

2) When OH was close to N and far from S (complex III) a H transfer from the amino group to  OH took place leading to the formation of water molecule and N-centered radical. This radical was the most stable one. This H transfer could occur directly or through a network of hydrogen bonds linking water molecules. When the amine was charged (-NH 3 + ) or replaced by an amide, the same result was obtained. It was independent of the method of solvent description. The N-centered radical did not undergo folding, excluding the formation of a 2c-3e S-N.

3) The first solvation shell has a great importance for the stabilization energies. With discrete water molecules I and II have very close energies while they differed by ca.

80 kJ/mol using only a CPCM.

4) VB analysis allowed obtaining probability of heterolytic rupture, which was modified by the presence of water molecules. It is notably important but not equal to 100% hence a mixture of cations and OH adducts should be present in solutions. The heterolytic rupture is specially probable with MetS- OH I, like in DMS adducts. In II, the 2c-3e S-N bond is prone to this heterolytic rupture because of stabilization by the presence of the 4 water molecules.

5) After having explored a large zone of geometries of the cation going from folded to extended, a small stability zone was found corresponding to 2c 3e SN bonded radical, consistently with our previous calculations [10,13]. Conversely, our former calculations [13,14] on peptides (MetGly, GlyMet and MetMet) led to several stable cyclic conformations (from 5 to 9 members) with MetSN  + as well asMetSO  + bonds.

The future of these free radical is to give birth to the final products that have been detected [15][16][17]. One should keep in mind that subsequent reactions take place in solution (e.g. disproportionation, hydration, decarboxylation…). For instance a sulfoxide can be created either from an OH adduct or from disproportionation and hydration. The CC double bond necessitates electron transfer from sulfur to adjacent carbon. Decarboxylation could come from electron transfer between N centered radical and carboxylate group.

Compliance with ethical standards

Conflict of interest: The authors declare that they have no competing interests.

  ) for the oxidation of methionine by • OH radicals. We focused on the role of solvation layer by adding discrete water molecules to a CPCM continuum description. In a recent study[23] some of us have simulated using Molecular Dynamics (MD) the positions of discrete water molecules and • OH radicals surrounding the simplest organic sulfur compound: dimethyl sulfide DMS. Here we have used the same simulation methodology for the methionine amino acid. The positions of the water molecules as well as various conformations of Met and the position of •OH radicals were obtained from MD simulations. Like in the previous study of the oxidation of DMS, waters within a radius of 4 Å around methionine were selected from several randomly selected snapshots of an equilibrated MD simulations. The • OH radical could be close to O and N atoms in different conformations of Met. We especially focused on two initial types of conformations: folded with the sulfur atom close to the nitrogen atom (in the MetSN species) or to one of the oxygen atoms (in the MetSO species). Thus, we explored the stabilities of the folded conformations of various Met derived • OH adducts (MetSN-OH and MetSO-OH).

Scheme 1 .

 1 Scheme 1. Schematic complexes of MetSX-OH (x = NH, O) extracted from the MD snapshots. Circles identify the location of the  OH radical close to atoms S, N or O of Met.

Figure 1 .

 1 Figure 1. Conformations of MetSX-OH (x = NH, O) extracted from the MD snapshots and optimized at the BH&HLYP/6-311+g(d,p) + CPCM level of theory. For all these conformations the bonds of 2.4 Å are represented by dotted lines. Circles identify the location of the  OH radical.

  of four discrete water molecules to MetSX- OH (X=N or O). We selected four water molecules directly linked to the atoms S, N, O of Met and to the oxygen of  OH radical from the first solvation shell of Met- OH, following the procedure used for the DMS adduct[23]. Figure2displays I, II, III and IV complexes surrounded by the 4 water molecules and

Figure 2 .

 2 Figure 2. Stable conformations MetSX- OH + 4H 2 O.MetSN-OH (:I; II and III) Color code: gray dotted lines: hydrogen bonds, black: S-N and S-OH 2c 3e bonds. Circles identify the location of the  OH radical.

Figure 3 displaysIVFigure 3 .

 33 Figure 3. Conformations of the initial point (obtained by MD) and an intermediate point of the optimization leading to III. Dotted lines: hydrogen bonds. Circles identify the location of the  OH radical.

Scheme 2 .

 2 Scheme 2. Schematic representation of the three possible VB resonance conformations for the MetSN, MetSO and MetExt adducts with • OH.

Figure 4 .

 4 Figure 4. Conformations of the optimized radical cations MetSX+y H 2 O (y= 11 or 4). Dotted lines: hydrogen bonds. The S-N bond is underlined by bold dotted line.

Figure 5 .

 5 Figure 5. Energy clusters Met + + x H 2 O with x= 7 (left) and x= 11 (right) with respect to both distances d SN (blue dots) and d SO (red dots). Conformations were optimized at the BHHLYP/6-31g(d)+ CPCM level of theory for each distance d SX from a folded conformer (d SX = 2.5Å) to unfolded conformer (d SX = 5.5Å) using a step of 0.25 Å.

Table 1 .

 1 MetSXE 1 is the relative energy with respect to the structure III given in kJ/mol. E 2 is the relative energy with respect to the structure I given in kJ/mol. Selected distances are given in Å. QTAIM Spin densities:  S ;  XH ;  OH and charges: Q S ; Q OH . δ S-N and δ S-OH are the delocalization indexes of the S-X and S-OH bonds are given in a.u.. As observed in Table1, the optimized complex I displays a small S-OH distance (2.34 Å) and spin densities are distributed between S and OH ( S = 0.38 a.u. and  OH = 0.61 a.u.). In the starting point of II, the hydroxyl radical was close to S and to N with S-N and S-OH distances were relatively short. After optimization, the SN bond became shorter. The spin densities are fairly distributed between S and N ( S = 0.56 a.u. and  NH = 0.38 a.u.). With regards to the delocalization index δ S-N , it remains lower than 0.5. These results are in good agreement with a 2c-3e bonding scheme for the S-N bond. Complex II is less stable than I by ≈ 80 kJ/mol.

			MetSN-OH		MetSO-OH
	E +954	-0.76835	-0.73997	-0.80062	-0.76909
	Complexes	I	II	III	IV
	E 1	84.6	158.1	0.0	82.7
	E 2	0.0	73.5		1.9
	d S-X	3.38	2.48	3.20	3.51
	d S-OH	2.34	3.10	4.20	2.39
	d N-H			1.86
	 S ;  XH ;  OH 0.38 0.0 0.61 0.56 0.38 0.0 0.0 0.96 0.0 0.32 0.00 0.66
	Q S ; Q OH	0.19 -0.63	0.30 -0.94	0.30 -0.94	0.14 -0.32
	δ S-X δ S-OH	0.04 0.62	0.48 0.10	0.06 0.02 0.00 0.56

-OH folded complexes (with X=NH for MetSN-OH or X= O for MetSO-OH) computed at the BHandHLYP/6-311+G(d,p) + CPCM level of theory. Electronic energies (E) are given in a.u..

Table 2 .

 2 MetSN-OH folded complexes III have been optimized at the BHandHLYP/6-311+G(d,p) + CPCM level of theory. Electronic energies are given in a.u. Selected distances are given in Å. QTAIM Spin densities  S ;  NH and charges Q S ; Q OH are given in a.u.. δ S-N and δ S-OH are the delocalization indexes of S-N and S-OH bonds given in a.u.. Columns 1 to 2: the end groups are methyls and column 3: the end group of N is the amide.

	 OH	III init	III	III amide
	E + 954	-0.73388	-0.80062	E = -1068.106410
	d S-N	2.52	3.20	2.56
	d S-OH	5.96	4.20	2.86
	d N-OH	3.44	1.86	4.12
	 S ;  NH	0.52 0.41	0.0 0.96	0.36 0.52
	Q S ; Q OH	0.24 -0.98	0.30 -0.94	0.15 -0.73
	δ S-N δ S-OH	0.10 0.00	0.06 0.02	0.44 0.00

Table 3

 3 gathers optimized energetic and structural parameters together with QTAIM analysis.

Table 3 .

 3 MetSX-OH folded complexes (with X = NH or O) + 4H 2 O are optimized at the BHHYP/6-311+G(d,p) + CPCM level of theory. Electronic energies (E) are given in a.u.. E is the relative energy with respect to the structure III given in kJ/mol. Selected distances are given in Å. QTAIM Spin densities  S ,  X ,  NH and charges Q S , Q OH . Delocalization indexes δ S-X, δ S-OH of S-X and S-OH bonds are given in a.u..

			MetSN + 4 H 2 O		MetSO + 4 H 2 O
	Complexes	I	II	III	IV
	E+1260	-0.48570	-0.48276	-0.52360	-0.49264
	E	99.4	107.1	0.0	81.2
	d S-X ; d S-OH	3.18 2.31	2.49 4.14	3.23 3.62	3.39 2.33
	 S ;  X ;  OH	0.46 0.00 0.51	0.52 0.21 0.00	0.00 0.97 0.00	0.48 0.00 0.50
	Q S ; Q OH	0.26 -0.54	0.23 -0.85	-0.04 -0.74	0.25 -0.44
	δ S-X δ S-OH	0.08 0.64	0.50 0.00	0.06 0.04	0.04 0.62

Table 4 .

 4 Electronic (E) of III + 4H 2 O structure along the H transfer has been computed at the

	Intermediate	(H-transfer)
	point

Table 6 .

 6 Table 6 gathers obtained structural parameters together with QTAIM analysis for the most stable structures. Most stable conformations of the scans SN/SO for MetSN  + and MetSO  + structures computed at the BHHYP/6-311+G(d,p) + CPCM level of theory. Electronic energies (E) are given in a.u.. ∆E is the relative energy given in kJ/mol with respect to the structure MetSN  + . Intramolecular distances S-N and S-O (both O of the carboxylate terminal) are given in Å.QTAIM spin densities  S ;  X=N or O and charges Q S ; Q X=N or O for radical cations MetS  + are given in a.u. Like for adducts, the 4 H 2 O mimic the first solvation shell and were selected from their short distances to the polar groups of Met, >COCH3, >CO, -NH and >S.

		MetSN  +	MetSO  +
	E + 878	-0.82585	-0.80412
	E	0.0	57.0
	d SN	2.49	3.45
	d SO	3.34	2.72
	d SO(CH3)	4.80	4.59
	 S ;  X=N or O	0.52 0.40	0.84 0.06
	Q S ; Q X=N or O	0.27 -1.06	0.46 -1.22
	The most stable cation is MetSN  + , with d SN =2.49 Å, a SN2c 3e bond and a relatively small
	distance between the sulfur and the oxygen of the carbonyl group d SO (3.34 Å). The
	equivalent with d SO = 2.7 Å has a much higher energy (more than 50 kJ/mol above that of
	MetSN  + ). The spin densities do not in agreement with a SO 2c-3e picture. Indeed, the
	structure is more consistent with a typical MetS  + cation in agreement with the strong spin
	density mainly localized on the sulfur atom (0.84 a.u.).
	Met + +4 H 2 O + CPCM. The energy of the cluster Met + + 4 H 2 O was calculated at the BH&HLYP/ 6-311+g(d,p) +
	CPCM level. The 4 water molecules were connected together by a H-bonds network. In
	Table 7, are summarized main structural parameters together with QTAIM analysis of
	optimized conformations of the MetSX

 + cations. The optimized conformations are displayed in Figure

4

. The found most stable conformer is still MetSN + , MetSO  + becoming MetS + .

Both conformations of the cations are quasi-isoenergetic (only 10 kJ/mol between them).

Table 7 .

 7 Most stable conformations of the scans SN/SO for MetSX + +4 H 2 O calculated at the

		MetSN + +4 H 2 O	MetSO + +4 H 2 O
	E + 1184	-0.53960	-0.52437
	E	0.0	39.9
	d SN	2.50	4.81
	d SO	5.09	2.39
	d SO(CH3)	5.25	4.49
	 S ;  X=N or O	0.50 0.42	0.82 0.12
	Q S ; Q X=N or O	0.23 -1.08	0.45 -1.21

BHHYP/6-311+G(d,p) + CPCM level of theory. Electronic energies (E) are given in a.u.. ∆E is the relative energy with respect to the structure MetSN + + 4 H 2 O given in kJ/mol. Intramolecular distances S-N, S-O (both O of the carboxylate terminal) are given in Å.

QTAIM spin densities  S ;  X and charges Q S, Q X for radical cations MetS  + are given in a.u..

These curves are very similar whatever the number of water molecules taken into account. The energy minimum for both curves is for a 2c-3e SN bond of 2.5 Å. Further data have been obtained with a crude analysis of the landscape of the potential electronic energy surface (PES), carried out for the clusters Met + + 11 H 2 O displayed in Figure 6. Three regions on the PES have been identified (see Figure 6): a small zone (purple to dark blue) corresponds to the most stable conformation. The SN bond is short (ca. 2.4 Å) and the SO one is long (ca. 5 Å). Interestingly the "point 1", which corresponds to the energetic minimum with only 4 water molecules, is situated in the purple zone.

The second one (green-blue) is divided in several parts, the SN bond of ca. 2.5 Å and the SO bond of ca. 4 to ca. 5.5 Å ; the SN and SO bonds around 4-5 Å; the SO bond smaller (ca. 3 Å)and SN around 4.5 Å. The energies have been found around 60 kJ/mol higher than that of the minimum. The "points 2 and 3" fall into these zones of high energy. The third one (greenyellow) is constituted of a large plateau with energies higher than 100 kJ/mol. This plateau is too high to allow transformation between SN and SO.

To conclude, the stable cations display S-N bonds rather than S-O bonds, these latter being much higher in energy. It is worth noting that a very high plateau separates the two different compounds of radical cations which suggests that no reasonable transformation from SO to SN bond is kinetically expected.


Riku Otono
email: otono.riku.ot6@is.naist.jp

Adélaïde Genay
email: adelaide.genay@inria.fr

Monica Perusquía-Hernández

Naoya Isoyama
email: isoyama@is.naist.jp

Hideaki Uchiyama
email: hideaki.uchiyama@is.naist.jp

Martin Hachet
email: martin.hachet@inria.fr

Anatole Lécuyer
email: anatole.lecuyer@inria.fr

Kiyoshi Kiyokawa

Studying "Avatar Transitions" in Augmented Reality: Influence on Sense of Embodiment and Physiological Activity

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

ABSTRACT

The impact of avatar transitions, in which an avatar's appearance is changed, has been little studied. We implemented an avatar transition system that is decomposed into two parts: an active one and a passive one. In the active transition, the avatar's appearance is transformed according to the user's physical action. In the passive transition, users automatically experience the same transition without physical action. The influence of transitions on the user's physical performance and the sense of embodiment was explored. Our results show that active transitions increase the sense of agency and magnitude of surface electromyography (sEMG) compared to passive transitions, and both transitions increase the sEMG magnitude following the virtual Reality embodiment.

Index Terms: Human-centered computing-Visualization-Visualization techniques-Treemaps; Human-centered computing-Visualization-Visualization design and evaluation methods

INTRODUCTION

Transition scenes are often used in video games to show the transformation of an avatar's appearance. However, few studies have investigated such transformations in Virtual Reality (VR) or Augmented Reality (AR) where a user embodies the avatar. Therefore, the effect of avatar transitions on the user's sense of embodiment (SoE) is not well known. Only one VR paper has studied transitions between being in one's body and being in one's avatar [START_REF] Jung | In limbo: The effect of gradual visual transition between real and virtual on virtual body ownership illusion and presence[END_REF]. They found that such transitions can have positive effects on body ownership illusions. We propose to continue exploring such transitions, but this time to change the appearance of one's avatar while embodying it. We call these transitions "avatar transitions."

To investigate the effects of avatar transitions, we implemented an AR embodiment system letting users transition from a normal avatar to a muscular avatar representing them in a third-person perspective (3PP), (see Fig. 1). We compared two types of avatar transitions by measuring the user's SoE and surface electromyography (sEMG) magnitude using a self-report questionnaire and sEMG sensors. We measured sEMG to see if the change in representation caused an increase in a physiological activity.

The two transition types we tested were: 1) an active transition where the user had to perform an action to trigger a visual transition scene, and 2) a passive transition where this scene was automatically triggered. Our results showed that 1) The active transitions had a positive impact on the sense of agency and the physiological activity when embodied in a muscular avatar, and 2) Both transitions positively influenced the physiological activity after embodiment.

METHODOLOGY

In the active transition, the users were asked to firmly squeeze a pair of dumbbells and then to relax while still holding it repeatedly. While doing this the user's avatar appeared to become more muscular every time they squeezed the dumbbells, (referred to as "visual feedback"). The avatar's muscularity morphed to increase by 20% when the sEMG values exceeded the users' average values, which were measured beforehand. After the users flexed their muscles up to five times, the avatar's body will reach its maximum muscle mass.

In the passive transition, the users could see their avatar's body becoming more muscular without any action on their part: their avatar gained 20% of muscle mass every five seconds and reached its maximum muscle mass after 25 seconds.

Users experience these transitions wearing a head-mounted display (VIVE Pro Eye), in which the AR view is rendered in Unity 2019.3.11f running on Windows 10. A 360 camera (Ricoh Theta V) was used to capture real-time 360 images.

EVALUATION

27 volunteers (11 females) were recruited (age: M=25.33, SD=3.54). We conducted a within-subjects design (with/without action × with/without visual feedback) and measured the difference between the SoE scores and sEMG from the following conditions.

• Condition VA: visual feedback and action. Users observed their avatar's appearance progressively become muscular through their actions. • Condition VN: visual feedback and no action. Users see their avatar's appearance automatically become progressively muscular. • Condition BA: no visual feedback and action. Users try to make their avatar more muscular, but a black screen was displayed until the end of the transformation. • Condition BN: no visual feedback and no action. Users see a black screen without any action until a muscular avatar appears at the end of the transformation.

The order of the conditions was counterbalanced with a Latinsquares design. We evaluated the SoE six questions on a 7point Likert scale related to its sub-components (agency, ownership, and self-location) [START_REF] Galvan Debarba | Characterizing first and third person viewpoints and their alternation for embodied interaction in virtual reality[END_REF]. In each condition, participants answered these questions twice: while embodied in the normal body avatar, and when embodied in the muscular body avatar. We computed scores for the SoE and then obtained the difference between the normal and muscular avatars by averaging the score difference.

To measure sEMG, EMG electrodes were attached on the extensor carpi ulnaris (ECU) muscle and on the Biceps. The average RMS was calculated as sEMG values per epoch. The sEMG values were measured in each condition (during embodiment in a normal, muscular avatar, and after embodiment). The muscular sEMG values ratio (EV R muscular ) was calculated by dividing muscular sEMG by normal sEMG. The post-hoc sEMG values ratio (EV R post hoc ) was calculated by dividing the "after" sEMG by the "normal" sEMG.

RESULTS

We performed a three-way ANOVA (gender: female vs. male; scene: visual feedback vs. black screen; and action: with action vs. without action) followed by a post-hoc analysis with the Holm-Bonferroni correction for all results. In terms of sense of agency, we observed the interaction effect of the scene and the action (F(1, 22) = 5.75, p < .05. See Fig. 2). The post-hoc analysis showed a significant difference between VA and BA (F(1, 23) = 12.34, p < .01) and BA and BN (F(1, 23) = 7.16, p < .05). This result implies that the active transition improved the sense of agency. We could not find significant differences between the active or passive transitions in terms of ownership and self-location.

The three-way ANOVA also revealed the interaction effect of the scene and the action on the left bicep EVR (F(1, 22) = 5.77, p < .05. See Fig. 3). Post-hoc analysis showed a significant difference between VA and BA (F(1, 23) = 8.48, p < .01) and VA and VN (F(1, 23) = 9.81, p < .01). Regarding the right bicep EVR, the interaction Therefore, these findings show that active transitions amplify the physiological activity during embodiment in a muscular avatar.

For both the left and right bicep's post-hoc EVR, a significant main effect of the scene was revealed (left bicep: F(1, 22) = 4.52, p < .05, right bicep: F(1, 22) = 6.64, p < .05). These results show that the visual feedback (instead of a seeing black screen) improves the physiological activity after embodiment in a muscular avatar even when physical action is not performed while transformation.

CONCLUSIONS AND FUTURE WORK

In this research, the impact of avatar transitions in AR was investigated. The implementation includes visual feedback of the changes in the avatar's appearance and triggering this change based on the user's physical activity. Using this system, we found that active transitions positively impacted the sense of agency and the physiological activity during embodiment. Furthermore, both transitions improved the physiological activity after embodiment. Our findings imply that the active transitions help users feel that they are in control of their avatars and feel empowered by the avatar's appearance. In addition, both transitions amplify the Proteus effect of a transformed avatar. This phenomenon could be applied to AR applications such as sports training using a muscular avatar. In the future, we will study avatar transitions using a different appearance avatar.

Figure 1 :

 1 Figure 1: We explored the effects of graphical transitions to embody an avatar. From the upper left to the upper right: female normal body avatar, avatar transition, and muscular body avatar. From the bottom left to the bottom right: a male version of upper avatars.

Figure 2 :

 2 Figure 2: From left to right, the results of the difference of agency, ownership, and self-location between the normal and muscular avatars.

Figure 3 :

 3 Figure 3: The upper plots depict the EV R muscular , showing the active transitions improve the muscular activity during embodiment. The lower plots depict the EV R post hoc , indicating the visual feedback improves the muscular activity after embodiment. effect of the scene and the action (F(1, 22) = 4.70, p < .05). Post-hoc analysis found a significant difference between VA and BA (F(1, 23) = 7.29, p < .05) and VA and VN (F(1, 23) = 5.62, p < .05).Therefore, these findings show that active transitions amplify the physiological activity during embodiment in a muscular avatar.For both the left and right bicep's post-hoc EVR, a significant main effect of the scene was revealed (left bicep: F(1, 22) = 4.52, p < .05, right bicep: F(1, 22) = 6.64, p < .05). These results show that the visual feedback (instead of a seeing black screen) improves the physiological activity after embodiment in a muscular avatar even when physical action is not performed while transformation.

ACKNOWLEDGMENTS

This research was supported in part by a Grant in aid for KAKENHI 22K12109 and the Start Up Fund for female researchers at NAIST.


Jérémie Vaubaillon
email: vaubaill@imcce.fr

Charlotte Loir

Clara Ciocan

Mathuran Kandeepan

Maxime Millet

Adrien Cassagne

Lionel Lacassagne
email: lionel.lacassagne@lip6.fr

Pedro Da Fonseca

Fabian Zander

David Buttsworth

S Loehle

J Toth

S Gray

A Moingeon

N Rambaux

A 2022 τ -Herculid meteor cluster from an airborne experiment: automated detection, characterization, and consequences for meteoroids

Keywords: meteoroids -Comets, general -method, data analysis

published or not. The documents may come

Introduction

Meteor clusters are the occurrence of many (typically more than three) meteors detected in a restricted portion of the sky within a few seconds (typically less than 5s; for a review see Koten et al. 2017, and references). The existence of meteor clusters has been suspected for decades, but only a handful (namely six) of observations have been reported [START_REF] Koten | [END_REF]. In particular, the Leonid storms occurring between 1998 and 2002 raised the question as to whether or not such observations happened simply by chance, given the high number of meteors recorded during each event. Evidence of their genuine existence has been reported that disfavors their observation simply being the result of statistical fluctuations [START_REF] Watanabe | [END_REF]Tóth & Klačka 2004). Given the ever growing number of meteor cameras surveying the sky around the globe every night, such events are expected to be more frequently reported. However, even with more than a thousand meteor detection cameras running today (see [START_REF] Koten | Meteoroids: Sources of Meteors on Earth and Beyond[END_REF], for a review of all the networks), meteor cluster observations are still very rare events.

The exact origin of meteor clusters is poorly known. A probable process is thermal stressing of very fragile comet dust (Watanabe et al. 2003), a hypothesis recently confirmed by Čapek et al. (2022) for the case of the 2016 September ϵ-Perseid (SPE) cluster. The meteoroid disruption drives the level of meteor showers, which itself depends on the structure of the comet. Jenniskens et al. (2008) reported a lack of fluffy meteoroid in an old Leonid trail, which these authors suggested is possibly explained by meteoroid self-disruption in the interplanetary space (although, other hypotheses might explain this observation). In order to explain the present quasi-steady-state of the level of sporadic meteors and of the amount of zodiacal dust, models must take the meteoroid life expectancy into account as well as the replenishment mechanism (Wiegert et al. 2009;[START_REF] Levasseur-Regourd | [END_REF]. Such mechanisms include the gravitational perturbation of long-period comets, the structure and population of the Oort cloud, the role of giant planets (especially Jupiter) in removing or accreting small bodies in the inner Solar System, and so on. Therefore, the frequency of meteoroid self-fragmentation in space has implications for our current understanding of the Solar System.

Since [START_REF] Koten | [END_REF], only one meteor cluster observation has been reported (U.Hawaii 2021), although an extensive search was recently performed among the Geminids (Koten et al. 2021). One open question refers to the frequency of spontaneous meteoroid breakup in interplanetary space, which would lead to meteor clusters and how these breakups would influence the lifetime expectancy of meteoroids. Here, we report another unambiguous detection of a meteor cluster of 34 fragments, detected within 7.5 s during the 2022 τ-Herculids outburst caused by the 1995 trail ejected from Jupiter family comet 73P/Schwassmann-Wachmann 3 (hereafter 73P). The detection was realized using a novel computer-vision application that was able to detect 100 % of the meteors that a human eye can see in the video. The results enable a discussion of the origin, frequency, and implications of such events.

Observations

Campaign and instrument

Ye & Vaubaillon (2022) predicted that the 2022 τ-Herculids meteor shower would be visible on 31 May 2022, and indeed it was successfully observed during an airborne observation campaign led by the University of Southern Queensland (F.Z., D.B.) and supported by Rocket Technologies International (S.G.). On board the aircraft, a suite of low-light scientific cameras were installed at several windows, and the sky was monitored continuously. A detailed description of the whole campaign is beyond the scope of this paper, but will be published in a dedicated paper. In addition to the imaging systems, spectroscopic systems were mounted in parallel. Windows on both sides of the aircraft were equipped with cameras, and the flight path was chosen according to the predictions (Ye & Vaubaillon 2022).

In this paper, we present results from data collected by a Basler acA1920-155um camera equipped with a Basler 6mm f/1.4 lens. The gain was set at maximum value (36) and 20 images per second were taken during the 4 hours of the flight. In order to compensate for the airplane (Phenom 300) roll motion, a G6-Max camera stabilizer was used. The camera was controlled with a RaspBerry-4 mini-computer, running the "RMS" acquisition and meteor detection software [START_REF] Vida | International Meteor Conference Egmond[END_REF](Vida et al. , 2021)). In addition, an AMOS-Spec-HR camera (Comenius Univ.) was mounted at another plane window. The hardware was a DMK 33UX252 (resolution of 2048 × 1536 px and set to 14fps) equipped with a 6 mm, F/1.4 lens, providing a FOV of 60 × 45 deg.

Detection of the meteor cluster

With the described settings, we detected 165 τ-Herculids meteors and five sporadic meteors. At the time of the shower outburst maximum (around 05:00 UT), we detected about one meteor per minute. When the level of the shower was decreasing, starting at 06:48:56 UT, the Basler camera detected 34 meteors, all coming from the τ-Herculids radiant. The AMOS camera, being slightly more sensitive, allowed the detection of 38 meteors within 11.3 s. Figure 1 shows a composite image of the meteor cluster (as detected by the Basler camera). This cluster observation was not reported by any of the ground-based meteor networks.

Meteor cluster characterization

The whole meteor cluster characterization was performed with the data from the Basler camera and is detailed in Table 1. Further characterizations of the whole shower are ongoing. The total time duration of the event is 11.3 s, but the Basler camera detected it for 7.5 s only. The maximum angular distance between all the meteors is ∼ 50 deg. The average airplane position during the cluster was lat=34.20 deg, lon=-101.88 deg, alt=14201 m and the camera was pointing towards its left hand side. The entry velocity was computed using the algorithm developed by Neslušan et al. (1998). Given the low entry velocity of the τ-Herculids (12.2 km.s -1 ), it is reasonable to assume an average meteor altitude of 90 km. Individual meteor azimuth and elevation (above the horizon) were measured. The relative apparent angular distance between each fragment is in the range [0.25; 50.2] deg, the measured elevation is within [41.2; 70.6] deg, and the physical distance between two fragments is within [0.4; 90.5] km. Adding the total duration of the event, the maximum possible physical distance between all the fragments is D m = [227; 244] km.

Following the methodology of [START_REF] Koten | [END_REF], we find that, assuming a Poissonian distribution of meteors (this assumption is discussed in Sect. 5.2), the probability of such clustering by chance is ∼ 5.5 × 10 -22 at best. As a result, we consider that the chance observation of this number of meteors during such a short time period is highly improbable and conclude that the disintegration of a parent τ-Herculid meteoroid took place in the interplanetary space.

Assuming a zero ejection velocity for all the fragments, the maximum time between the parent meteoroid disintegration in interplanetary space and the Earth atmosphere entry strongly depends on the considered meteoroid size. We converted the apparent magnitude into an absolute magnitude (assuming a meteor altitude of 90 km), then converting this latter into an equivalent photometric mass (using Hughes 1995) and radius. The latter ranges from 7.5 mm to 22.4 mm (assuming a density of 2500 kg.m -3 ). The total mass of the initial meteoroid is estimated to be 1.16 kg. The maximum age of the cluster is computed using the smallest size, as this is the most sensitive to the solar radiation pressure, and is found to lie in the range [13.3; 13.8] days.

Figure 2 shows the absolute magnitude distribution of the cluster. The population index is r = 2.01, corresponding to a differential size distribution of s = 3.09. This feature is discussed in Sect. 5.1.

Computer-vision detection

The real-time open-source software detection chain named Fast Meteor Detection Toolbox (FMDT) was applied in order to detect the meteor in the imaging data 1 . FMDT is derived from software designed to detect meteors on board the ISS or a Cubesat (Millet et al. 2022a,b;[START_REF] Petreto | IEEE International Conference on Design and Architectures for Signal and Image Processing (DASIP)[END_REF]. FMDT is foreseen to be applied to airborne camera systems; for example in atmospheric balloons or aircraft. It is robust to camera movements thanks to a motion-compensation algorithm.

Figure 3 presents the whole FMDT detection chain. For each pair of images, an intensity hysteresis threshold, a connected component labeling, and an analysis algorithm [START_REF] Lemaitre | ACM Workshop on Programming Models for SIMD/Vector Processing (PPoPP)[END_REF][START_REF] Lacassagne | IEEE International Conference on Image Analysis and Processing (ICIP)[END_REF] are applied to get a list of connected components (CCs) with their bounding boxes and surface S . Moreover, it also provides the first raw moments to compute the centroid (x G , y G ) = (S x /S , S y /S ) of each blob of pixels. A morphological threshold is then set on the surface S to reject small and big CCs. A k-nearest neighbor matching process is then applied to extract pairs of CCs from images I t+0 and I t+1 with t the image number in the video sequence. These matches This motion estimation is used to classify the CCs into two classes, namely nonmoving stars and moving meteors, according to the following criterion: |e k -ēt | > σ t where e k is the compensation error of the CC number k, ēt the average error of compensation of all CCs of image I t , and σ t the standard deviation of the error. A second motion estimation is carried out with only nonmoving star CCs in order to obtain a more accurate motion estimation and a more robust classification. Finally, piece-wise tracking is carried out by extending the (t + 0, t + 1) matching with (t + 1, t + 2) matching to reduce the number of false-positive detections. For the present video data, the geometric mean error e t for the whole sequence is 0.91 pixels for the first estimation and 0.18 for the second one. The apparent speed varies from 3 up to 10 pixels/frame.

For the considered video sequence, FMDT was able to detect and track 100% of the meteors in the video that are visible to the naked eye, with only four false positives. The proposed solution was compared with a manual detection (where an expert watched and labeled the entire video). This manual detection constitutes the "ground truth" and was first able to detect 28 meteors, with meteors 4, 8, 16, 18, 20, and 29 being missed. The ground truth was then enhanced thanks to the automatic detection chain. This demonstrates the need for an automated system for meteor detection. Figure 4 shows the overlap between meteors detected automatically and those of the ground truth. We define the tracking rate T r as the ratio of the cumulative duration of the automatically detected meteors and of the cumulative duration of the ground-truth meteors:

T r =         34 m=1 D auto-detect m         /         34 m=1 D ground truth m        
, where D m is the duration of the considered meteor m. In the observed video sequence, T r = 80.4%. A video of the sequence with meteor tracking is available online2 . For most of the meteors, the automatic detection is very close to reality. Moreover, the minimum time required for a manual detection is close to the full time of the video sequence while the proposed application is real-time and compatible with the cubesat power consumption constraint. Moreover, FMDT is able to leverage multi-core processor architectures through a task graph description and the use of the AFF3CT multi-threaded runtime library (Cassagne et al. 2019[START_REF] Cassagne | Proceedings of the International Symposium on Topics in Coding (ISTC) (IEEE)[END_REF][START_REF] Cassagne | [END_REF]. AFF3CT was designed for digital communication systems but is well adapted to real-time image processing.

Discussion

Meteoroid properties

The measured cluster differential size distribution index s = 3.1 is slightly lower than expected from a collision cascade (3.5; see e.g., O'Brien & Greenberg 2005). We examine how the derived s value compares to other measurements. Reanalyzing the 2016 SPE cluster, of which the parent body is unknown, Čapek et al. ( 2022) find a shallow s = 1.85. Similarly, for decameter-size fragments ejected by 73P, Reach et al. (2009) found two relatively low size distributions of s = 1.84 and 2.56 for the smallest and the largest fragments, respectively. However, given the suspected rocket effect involved in such an event, the physical process is presumably different from what is at play in a meteor cluster. Measurements for comet 67P/Churyumov-Gerasimenko provide a wide range of values, depending on the comet heliocentric distance and meteoroid size range (see Güttler et al. 2019, and references for a review). For sizes comparable to meteoroids responsible for visual meteors, extreme low values of s = 1.8 were derived when the comet was at high heliocentric distance (Merouane et al. 2016). Higher values of s > 3.5 were found after perihelion [START_REF] Fulle | [END_REF]Moreno et al. 2017;Fulle et al. 2010). Last but not least, an extreme high value of s = 6.4 was found by Vida et al. (2021) for a meteoroid fragmenting in the atmosphere. Interestingly, extreme values of s are derived for two drastically different meteoroid environments. The lowest s values are found when the meteoroid breaks up in interplanetary space, and very high values are found when this takes place in the Earth's atmosphere. Whether this difference reflects the way meteoroids interact with gaseous environment is unclear, and investigating this matter would require additional work.

The value reported here is high compared to [START_REF] Koten | [END_REF], but is still at the low end of all reported values. Güttler et al. (2019) provide a review of the literature for 67P, and recall that meteoroids do not fragment in the coma. It is worth pointing out that the size distribution changes as a function of heliocentric distance. If this is true for all comets, and as the meteoroids of a given trail are ejected at different heliocentric distances, a meteoroid trail is therefore composed of a wide variety of meteoroid subtrail each described with a unique size distribution. In addition, within a trail, the meteoroids are mixed because of the relatively wide range of sizes and ejection velocity vectors. The portion of meteoroid subtrail sampled by the Earth during a meteor shower is therefore a mixture of all these size distributions. The size distribution s has a tremendous influence on the level of a shower (number of meteor per unit of time, Vaubaillon et al. 2005a,b). Future work is needed to quantify the influence of modeling a variable size distribution for the prediction of the meteor showers, and how this might reconcile past post-predictions with observations (e.g., the 2006 Leonids).

Meteoroid cluster frequency

Meteoroids are known to fragment in the atmosphere with very high probability (90%; see Subasinghe et al. 2016). The reason for a cluster observation is (see above) an interplanetary fragmentation event. However, this is very rarely observed: only seven clusters have been reported in the past ∼ 40 years of meteor observations. With the new data from the τ-Herculids campaign, we find that the probability of such a cluster observation is P ∼ 5.5 × 10 -22 (see Sect. 3). [START_REF] Sampson | [END_REF] points out that the assumption of a Poisson distribution might not be appropriate, and generally underestimates P. However, if in the case presented here an extreme error of a factor 10 9 is assumed, this still leads to P ∼ 10 -13 , showing the extreme rarity of the phenomenon.

Computing a meteor cluster observation frequency would require the consideration of the limiting magnitude as a function of time, the software detection efficiency, and the camera running efficiency, among others, for each meteor-detection camera. As such a thorough study is out of the scope of this paper, we attempt to provide an order-of-magnitude estimate. In the past 20 years, ∼ 2 × 10 6 meteors were observed by the IMO during a total effective observation time of ∼ 8 × 10 6 hours (Molau 2021). The EDMOND database currently3 counts ∼ 4.6 × 10 6 meteors gathered between 2000 and 2016 [START_REF] Kornoš | Proceedings of the International Meteor Conference[END_REF]. The SonotaCo and GMN networks have recorded totals of respectively ∼ 3.5 × 10 5 and ∼ 2.2 × 10 5 meteors over the past 14 years (SonotaCo et al. 2021;Vida et al. 2021). During this time, only a handful of clusters were reported [START_REF] Koten | [END_REF]). From our experience, meteor-detection software from video data (RMS, UFOCapture and FreeTure) is able to detect more than one meteor in a given frame. This is enough to conclude that the occur- rence of a meteor cluster happens with a frequency of less than one in a million meteors.

All reported clusters happened during a meteor shower (Watanabe et al. 2003;Koten et al. 2017, and this work). Their orbits cover all possible cometary orbits: τ-Herculids for Jupiter family comet(JFC)-type, Leonids for Halley type (HT), and September Perseids for long-period (LP) type. The time between the disruption in the interplanetary space and the entry in the atmosphere was estimated to be only a few days. This represents less than 0.3% of the orbital period of a JFC meteoroid. A short heliocentric distance presumably increases the chances of selffragmentation of meteoroids, given the higher micro-meteoroid space density, higher temperature and thermal stresses, and generally higher influence of radiation on their rotation state.

Origin of meteoroid self-fragmentation

The cluster presented here became visible nearly 2 h after the expected maximum τ-Herculids shower outburst caused by the trail ejected from comet 73P in 1995. No encounter with any other trail was expected at this time. The age of the parent meteoroids cannot be pinpointed, but given the lifetime expectancy of Jupiter family streams, this is likely to be a few hundred years at most [START_REF] Vaubaillon | Meteoroids: Sources of Meteors on Earth and Beyond[END_REF]. Out of the currently eight meteor cluster detections (including this study), only [START_REF] Piers | [END_REF] was not related to a known meteor shower. The extreme fragility of some cometary meteoroids (Hornung et al. 2016) might explain this feature. The often-quoted physical processes responsible for meteoroid fragmentation in interplanetary space are thermal stresses, collision, rotational outburst, and outgassing of volatile material. The processes involved in the natural release of meteoroids from an active asteroid were described by [START_REF] Jewitt | Asteroids IV (Asteroids IV[END_REF]. Čapek et al. (2022) found that thermal stress was most probably responsible for the 2016 SPE meteor cluster.

Future application of the developed algorithm

In addition to a cluster detection, we present a first application of the new processing chain for meteor detection named FMDT. This toolbox is derived from the CubeSat project Meteorix dedicated to the detection of meteors and space debris from space [START_REF] Rambaux | ESA NEO and Debris Detection Conference-Exploiting Synergies -ESA/ESOC[END_REF](Rambaux et al. , 2021)). This detection chain allow the realtime identification of meteors and enable autonomous selection of scientific data to be sent back to Earth from on board a Cube-Sat. The full chain also contains an optical flow algorithm for accurate motion estimation. The agreement in detections between the "RMS" software [START_REF] Vida | International Meteor Conference Egmond[END_REF](Vida et al. , 2021) ) and the new approach proposed by our team (Millet et al. 2022a,b;[START_REF] Petreto | IEEE International Conference on Design and Architectures for Signal and Image Processing (DASIP)[END_REF] allows us to test and validate the approach implemented and to increase the Technology Readiness Level to 5. Such a tool might be used for future detection of meteors from orbiting spacecraft (using e.g., the SPOSH camera; Bouquet et al. 2014;Oberst et al. 2011), or more generally from mobile observation platforms [START_REF] Vaubaillon | [END_REF].

Conclusion

We describe and fully characterize the eighth meteor cluster ever reported. Based on our analysis of the observation data, we conclude that the probability of a cluster meteor observation is less than one in a million observed meteors. The derived differential size distribution index s = 3.1 is relatively shallow. This index varies with heliocentric distance for regular comet outgassing. Future meteor-shower-prediction models might take this phenomenon into account for better accuracy.

We developed an open-source computer-vision-based toolbox, namely the Fast Meteor Detection Toolbox (FMDT) to detect and track meteors. In spite of the acquisition camera instability caused by the aircraft, it was able to detect 100 % of the meteors that are detectable in the video with the naked eye, even those of high magnitude.

Fig. 1 .

 1 Fig. 1. Composite closeup view of the detection of the 34 τ-Herculids meteor cluster from an airborne observation campaign. The center of the field of view points toward the little dipper (Ursa Minor) constellation.Table 1. τ-Herculid meteor cluster characterization: beginning time in UT, duration (D) in seconds, and apparent magnitude (M). Both groundtruth and automatic detection results (see Sec. 4) are reported. Ground truth Automatic detection Ground truth Automatic detection # Beginning D Beginning D M # Beginning D Beginning D M 1 06:48:55.959 0.30 06:48:55.959 0.30 -0.19 19 06:48:59.009 0.20 06:48:59.059 0.15 0.62 2 06:48:56.359 0.50 06:48:56.409 0.35 -1.50 20 06:48:59.209 0.10 06:48:59.209 0.10 1.79 3 06:48:56.359 0.80 06:48:56.359 0.80 -2.06 21 06:48:59.359 0.25 06:48:59.409 0.20 -0.21 4 06:48:56.909 0.10 06:48:56.909 0.10 0.87 22 06:48:59.559 0.35 06:48:59.559 0.30 -0.56 5 06:48:57.209 0.10 06:48:57.209 0.10 1.19 23 06:48:59.759 0.35 06:48:59.759 0.10 0.73 6 06:48:57.259 0.20 06:48:57.309 0.10 0.92 24 06:48:59.809 0.15 06:48:59.809 0.10 1.33 7 06:48:57.509 0.45 06:48:57.509 0.40 -1.07 25 06:48:59.809 0.50 06:48:59.809 0.50 -0.65 8 06:48:57.559 0.15 06:48:57.599 0.15 0.86 26 06:48:59.859 0.25 06:48:59.859 0.20 0.41 9 06:48:57.559 0.45 06:48:57.599 0.25 -1.40 27 06:49:00.009 0.35 06:49:00.009 0.30 -1.72 10 06:48:57.609 0.15 06:48:57.659 0.10 0.27 28 06:49:00.559 0.15 06:49:00.559 0.15 0.70 11 06:48:57.709 0.45 06:48:57.809 0.25 -0.88 29 06:49:00.609 0.30 06:49:00.709 0.10 -1.09 12 06:48:57.809 0.20 06:48:57.809 0.15 0.92 30 06:49:00.709 0.25 06:49:00.759 0.10 0.26 13 06:48:57.859 0.50 06:48:57.859 0.50 -0.68 06:49:00.809 0.15 14 06:48:58.159 0.15 06:48:58.159 0.15 1.44 31 06:49:00.809 0.35 06:49:00.859 0.25 -0.31 15 06:48:58.559 0.55 06:48:58.659 0.45 -1.32 32 06:49:02.009 0.30 06:49:02.009 0.30 0.27 16 06:48:58.609 0.20 06:48:58.659 0.10 -0.31 33 06:49:02.059 0.25 06:49:02.059 0.20 -0.82 17 06:48:58.709 0.30 06:48:58.709 0.25 -0.01 34 06:49:03.309 0.15 06:49:03.309 0.15 0.56 18 06:48:58.859 0.15 06:48:58.859 0.15 0.85

Fig. 2 .

 2 Fig. 2. Cumulative absolute magnitude distribution of the cluster fragments.

Fig. 4 .

 4 Fig. 4. Overlapping duration of the ground truth and the automatic detection for the meteor cluster. The automatic detection bars are placed relative to the beginning time of each meteor.

  Computer-vision detection detailed chain. Plain gray boxes correspond to input data, plain white boxes are the processing, and the italicized brown texts are the processing outputs.

	image (t + 0) image (t + 1)	intensity hysteresis thresholding intensity hysteresis thresholding blob of pixels	connected components (CCs) labeling & analysis connected components (CCs) labeling & analysis CCs = list of regions, surface S & centroid (x G , y G )	morphological thresholding morphological thresholding with S ∈ [S min , S max ] sub-list of CCs	k-nearest neighboor matching list of (t + 0,t + 1) matchings	global motion estimation translation & rotation	residual motion computation CCs speed (v x , v y ) & stats ēt , σ t	motion thresholding initial list of meteors / stars	piece-wise tracking final list of meteors
	0 0.1 0.2 0.3 0.4 0.5 0.6 Fig. 3. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 Duration (s) 0.8 Ground truth (ref.) 0.7 Automatic detection.
					Meteor id (#)

FMDT repository: https://github.com/alsoc/fmdt

Meteor cluster sequence with highlighted detection: https://lip6.fr/adrien.cassagne/data/tauh/tracks.mp4

https://www.meteornews.net/edmond/edmond/edmond-database/, accessed on 28 July 2022

Acknowledgements. MoMet is supported by "Programme National de Planetologie" and IMCCE / Observatoire de Paris / PSL. Airborne campaign organized and funded through RTI and supported by University of Southern Queensland. Maxime Millet PhD grant is funded by Region Ile-de-France. Dr. Fabian Zander is funded by the Australian Research Council through the DECRA number DE200101674. Dr. Juraj Toth was supported by ESA contract No. 4000128930 /19/NL/SC, the Slovak Research and Development Agency grant APVV-16-0148, the Slovak Grant Agency for Science grant VEGA 1/0218/22.


Angela Parise

Aurelio Alvarez-Ibarra

Xiaojing Wu

Xiaodong Zhao

Julien Pilmé

Aurélien De La Lande
email: auréliendelalande@u-psud.fr

Quantum Chemical Topology of the Electron Localization Function in the Field of Attosecond Electron Dynamics

We report original analyses of attosecond electron dynamics of molecules subject to collisions by high energy charged particles based on Real-Time Time-Dependent-Density-Functional-Theory simulations coupled to Topological Analyses of the Electron Localization Function (TA-TD-ELF). We investigate irradiation of water and guanine. TA-TD-ELF enables qualitative and

quantitative characterizations of bond breaking and formation, of charge migration within topological basins, or of electron attachment to the colliding particle. Whereas the Lewis-VSEPR structure of gas phase water is blown out within a few attoseconds after collision, that of guanine is far more robust and reconstitutes rapidly after impact even though the molecule remains electronically excited. This difference is accounted by the presence of the electron bath surrounding the impact point which enables energy relaxation within the molecule. Our approach should stimulate future studies to unravel the early steps following irradiation of various types of systems (isolated molecules, biomolecules, nanoclusters, solids…) and is also readily applicable to irradiation by photons of various energies. TOC GRAPHICS KEYWORDS quantum chemical topology, attosecond electron dynamics, radiation damages, DNA base, Real-Time TDDFT.

High energy charged particles (electrons, positrons, protons, alphas, hadrons…) interact with matter in turn inducing electronic excitations or ionization of the constituting molecules [START_REF] Mozumder | Interaction of Fast Charged Particles with Matter[END_REF] . A fine understanding of irradiation is essential in various scientific domains. In astrochemistry ionizing radiations (IoR) deeply impact the complex chemistry occurring in gases or on grains [START_REF] Dartnell | Ionizing Radiation and Life[END_REF] . In biology IoR cause deleterious damages to genomes, proteins or signaling pathways that normally regulate cell activity, having harmful consequences like accelerated aging, tumors and cancers [START_REF] Kryston | Mutation Research/Fundamental and Molecular Mechanisms of Mutagenesis. Role of oxidative stress and DNA damage in human carcinogenesis[END_REF] .

Let us also mention that protection of astronauts from IoR is the major obstacle for planning human interplanetary travels. On the opposite, IoR have beneficial effects in proton-or hadrontherapies or for sterilization purposes in industry. Interaction of IoR with matter and their subsequent ultrafast consequences continues to challenge both experimental and computational approaches. Massive efforts push the development of attosecond sources, which in turn have started to stimulate new developments in computational physical chemistry [START_REF] Lepine | Attosecond molecular dynamics: fact or fiction?[END_REF] . In this Letter we are interested in deciphering the attosecond response of molecules' electron clouds submitted to ionizing radiations. We emphasize that collisions with charged particles are very different from absorptions of X-ray or UV/visible photons or from interactions of molecules with strong laser fields. Collision with keV-MeV alpha produces superposition of highly excited electronic states.

The energy deposited locally into the electronic cloud is typically of several tens of eV. Although useful, the analysis of electron dynamics based on the occupations of molecular orbitals (MO) is cumbersome for these processes because of the very large number of involved MOs. The Electron Localization Function (ELF) [START_REF] Becke | A simple measure of electron localization in atomic and molecular systems[END_REF][START_REF] Savin | Electron Localization in Solid-State Structures of the Elements: the Diamond Structure[END_REF] and its extension to the time domain [START_REF] Burnus | Time-dependent electron localization function[END_REF] are alternative appealing tools. Two-or three dimensional plots of TD-ELF have been reported to describe alterations of molecular density fields when molecular systems interact with ultra-short laser pulses or collide with slow ions [START_REF] Castro | Time-dependent electron localisation function: A tool to visualise and analyse ultrafast processes[END_REF][START_REF] Penka | Ionization and harmonic generation in CO and H 2 CO and their cations with ultrashort intense laser pulses with timedependent density-functional theory[END_REF][START_REF] Bruner | Attosecond Charge Migration with TDDFT: Accurate Dynamics from a Well-Defined Initial State[END_REF] . TD-ELF actually incorporates a valuable wealth of information on the time evolution of the chemical structure that is potentially accessible by quantum topological analyses [START_REF] Castro | Time-dependent electron localisation function: A tool to visualise and analyse ultrafast processes[END_REF][START_REF] Pilmé | Topological analyses of time-dependent electronic structures: application to electron-transfers in methionine enkephalin[END_REF] . We extend the analysis of ELF within the framework of gradient field theory which was initially developed for stationary electronic structures [START_REF] Piquemal | What can be learnt on biologically relevant systems from the topological analysis of the electron localization function?[END_REF][START_REF] Silvi | Classification of chemical bonds based on topological analysis of electron localization functions[END_REF] to nonstationary electronic structures and apply it for the first time to irradiation of molecules by high energy alpha particles. Our approach enables, for instance, to qualitatively and quantitatively characterize formation/breaking of bonds between nuclei, electron flows (charge migration) among topological basins, or the attachment of electron density to the colliding particle. To illustrate the strength of the approach we successively consider collision of water or of a DNA base (guanine) by alpha particles and analyze the evolution of the ELF topological basins in combination with other descriptors of the electron dynamics.

Our electron dynamics simulations are based on Real Time TD-DFT [START_REF] Runge | Density-Functional Theory for Time-Dependent Systems[END_REF] (Time-Dependent Density Functional Theory) which is a powerful theoretical approach to address attosecond dynamics of matter submitted to irradiation by laser field or charged particles [START_REF] Zhang | Theoretical investigation of He 2+ -Ar collisions in the energy range of 4-300 keV/amu[END_REF][START_REF] Hong | H + -H 2 O collisions studied by timedependent density-functional theory combined with the molecular dynamics method[END_REF][START_REF] Baxter | Time-dependent spin-density-functional-theory description of He + -He collisions[END_REF][START_REF] Nagano | Application of the time-dependent local density approximation to collision between a highly charged ion and an atom[END_REF][START_REF] Dinh | Dynamics of clusters and molecules in contact with an environment[END_REF][START_REF] Wang | Microscopic studies of atom-water collisions[END_REF][START_REF] Gao | Theoretical study on collision dynamics of H+ + CH4 at low energies[END_REF][START_REF] Wang | Collision dynamics of proton with formaldehyde: Fragmentation and ionization[END_REF][START_REF] Covington | Time-dependent density-functional-theory investigation of the collisions of protons and $\ensuremath{\alpha}$ particles with uracil and adenine[END_REF] . We focus on inelastic collisions that are known to be largely dominant for kinetic energies of few tens of keV and higher [START_REF] Mozumder | Interaction of Fast Charged Particles with Matter[END_REF] . The colliding particle is described as a point charge moving in space at constant velocity. The TD Coulomb potential created by the particle is included in the external Kohn-Sham potential. We use an implementation of RT-TDDFT in the deMon2k software [START_REF] Köster | deMon2k Version[END_REF] (see Reference [START_REF] Wu | Simulating Electron Dynamics in Polarizable Environments[END_REF] for details). All calculations have been done with the GGA (Generalized-Gradient-Approximation) CAP-PBE exchange-correlation potential [START_REF] Carmona-Espíndola | Generalized gradient approximation exchange energy functional with correct asymptotic behavior of the corresponding potential[END_REF][START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] . CAP stands for Correct Asymptotic Potential and exhibits a correct asymptotic behavior [START_REF] Carmona-Espíndola | Generalized gradient approximation exchange energy functional with correct asymptotic behavior of the corresponding potential[END_REF] . The capabilities of various XC functionals to fully capture the physics involved ionization by high energy particles

are not yet fully understood. Global hybrids or range-separated hybrids would provide better description of Rydberg or charge separated states than GGA which constitutes a limitation of the present protocol. More technical details on the simulation protocol can be found in SI.

The quantum chemical topology of ELF relies on the analysis of the gradient field of a scalar function by applying the theory of dynamical systems [START_REF] Silvi | Electron Density Analysis A2 -Reedijk[END_REF] . The latter provides a partitioning of the molecular space into non-overlapping volumes (basins). The ELF topology displays atomic and non-atomic maxima that can be related to chemical concepts issued from the Lewis theory. The core basins, labeled C(X) are located at nuclear centers, whereas the valence basins are located in the remaining space. These valence basins are characterized by a synaptic order, i.e. the number of core basins with which they share a common boundary [START_REF] Silvi | Classification of chemical bonds based on topological analysis of electron localization functions[END_REF] . Monosynaptic basins, labeled V(X), usually associated to the lone-pair regions of X atoms, while disynaptic basins, labeled V(X, Y), correspond to two-center covalent X-Y bonds. Basins with a synaptic order of zero (termed as asynaptic) are sometimes found in the molecular space. They are labeled as V(A) in this work.

The integration of the electron density over the basin volume provides a basin population or intrinsic moments [START_REF] Pilmé | Advancing beyond charge analysis using the electronic localization function: Chemically intuitive distribution of electrostatic moments[END_REF] . Topological analyses have been performed with a modified version of the TOPCHEM package [START_REF] Kozlowski | New insights in quantum chemical topology studies using numerical grid-based analyses[END_REF] .

Figure 1 depicts the total ionization cross section for water as a function of the kinetic energy of the colliding particle. Each point results from an average over one hundred trajectories with different initial positions of the alpha particle and impact parameter (see SI for details). RT-TDDFT results are in excellent agreement with the available experimental data [START_REF] Rudd | Ionization cross sections for 10-300-keV/u and electron-capture cross sections for 5-50-keV/u 3 He 2+ ions in gases[END_REF] . Notably the maximum ionization probability around 0.1 MeV is recovered. RT-TDDFT results are also in good agreement with perturbative approaches such as the first Born Approximation Continuum Distorted Wave (FBA-CDW) or the Continuum Distorted Wave-Eikonal Initial State (CDW-EIS) methods [START_REF] Martinez | The continuum distorted wave-eikonal initial state model for single electron capture in ion-atom collisions[END_REF][START_REF] Galassi | The continuum distorted wave eikonal initial state model for transfer ionization in H + , He 2+ + He collisions[END_REF][START_REF] Rudd | Ionization cross sections for 10-300-keV/u and electron-capture cross sections for 5-50-keV/u 3 He 2+ ions in gases[END_REF][START_REF] Quinto | Single ionization and capture cross sections from biological molecules by bare projectile impact[END_REF] . At higher energies the slope of the RT-TDDFT curve follows that of standard models. The apparent capability of our RT-TDDFT implementation based on the CAP-PBE functional to simulate ionization by charged particles is in line with a few previous analogous simulations of collisions of water with fast protons [START_REF] Hong | H + -H 2 O collisions studied by timedependent density-functional theory combined with the molecular dynamics method[END_REF] . These are encouraging results

for future applications to large molecular systems. seem to be the most severe, which corresponds to collision energies for which the amount of energy transferred to the molecule reaches a maximum (Figure S2). Note that the core basin C(O) remains undisturbed whatever the impact energy. On contrary to the stationary topology, the ELF valence regions seems to be scattered within the molecular space because many asynaptic valence basins (V(A)), not clearly associated to atoms, are observed (see Figure 3).

These asynaptic basins appear as a typical signature of the intensity of collisions, although they are associated to small topological patches containing little electron densities. As a second example, we report the collision of a guanine molecule impacted by an alpha particle travelling along a line perpendicular to the aromatic cycles. The kinetic energy of the particle is chosen to be 0.1 or 1MeV, the former being close to the maximum ionization probability as predicted by CDW-EIS models [START_REF] Quinto | Single ionization and capture cross sections from biological molecules by bare projectile impact[END_REF] . We also consider two impact points for guanine, either close to the central C=C bond (case i) or to the N-H group of the 6-member aromatic cycle (case ii). By convention the collision time corresponds to snapshot a. We focus first on case i) with a 1MeV particle and on the very moment of collision (snapshot a-b). We also Collision induces partial depopulation of all initially doubly occupied MOs, except those corresponding to core electrons. The latter point is corroborated by the ELF analysis on core basins (Figure S12). Meanwhile a wide manifold of initially unoccupied MOs acquire electronic populations, typically around 0.05 e -. We observe continuous fluctuations of occupation numbers in the rest of the simulation, indicating that the molecule remains in a highly excited state.

Actually only a moderate fraction of the energy initially deposited by the collision is dissipated by ionization (Figure S3). The electronic kinetic energy integrated over atoms brings further insights into the relaxation mechanism. As expected, it is the atoms closest to the impact point that are the first ones to have their electronic kinetic energy increased, followed by fast (ca 100 as) energy damping and ultimately persistent oscillations on longer times. The other atoms respond with a delay depending on the distance from the impact point, the closer the atom to the impact the faster the local increase of electronic kinetic energy. The aforementioned features are signatures of intense dynamics taking place within the electron cloud even several hundreds of attoseconds after collision. This is also apparent in the movies showing the fluctuating TD-ELF isosurfaces. Somewhat surprisingly the Lewis structure of the molecule appears to be overall robust (see Figure 4 f ). Simulations corresponding to case ii) (impact at the level of 6-members ring NH group) with 1MeV particles lead to qualitatively similar conclusions. The main difference is that we observe initial disruption of the Lewis-VSEPR structure around the nitrogen atom, before relaxation. With 0.1 MeV alpha particles the energy delivered to the system is much higher (by a factor of almost two, Figure S3) and the electron cloud is clearly more affected. For instance the number of electron loss 2fs after collision, goes from ca 0.3 e-(1MeV) to 1.2 e-(0.1

MeV). Here again the ELF topological structure remains qualitatively robust (Figure S13). Future studies will have to thoughtfully investigate the influence of environment on the very early steps following interaction of matter by IoR for example to characterize charge migrations or other ultrafast phenomena [START_REF] Alexander | Ultrafast correlation-driven electron dynamics[END_REF] . Recent implementations of RT-TDDFT within hybrid QM/MM schemes [START_REF] Wu | Simulating Electron Dynamics in Polarizable Environments[END_REF][START_REF] Morzan | Electron dynamics in complex environments with real-time time dependent density functional theory in a QM-MM framework[END_REF][START_REF] Chen | QM/MM Study of Photoinduced Reduction of a Tetrahedral Ag 20 + Cluster by a Ag Atom[END_REF] should find natural applications in that context. The herein presented analyzing tools should find a prominent place in such future investigations. We finally note that the herein developed approaches are readily application to other types of irradiations for example by photons.

ASSOCIATED CONTENT

Computational details on RT-TDDFT simulations. Movies showing ELF isosurfaces for selected trajectories of water and guanine molecules (mpeg format). The following files are available free of charge.

Figure 1 :

 1 Figure 1: ionization cross-section for gas phase water. Ionization probability is calculated 1fs after collision. RT-TDDFT curve normalized to the experimental value for 0.1MeV. Blue line: RT-TDDFT; yellow: CDW-EIS model 34 ; green: FBA-CDW model, orange triangles: experiments from Rudd et al 31 .On the stationary electron density of isolated water obtained by a SCF (Self-Consistent-Field)

Figure 2 :

 2 Figure 2: ELF isosurfaces (0.88) on the stationary electron densities of H 2 O isolated (Top) and in presence of He nucleus ( particle, in yellow, Bottom). ELF topological basin populations are given in electron numbers.Interestingly for energies of 1 and 100 MeV no ELF basins encompassing the helium nucleus are

Figure 3 :

 3 Figure 3: TD-ELF for water molecule impacted by alpha particles of various kinetic energies. For each topological basin we indicate, in order, the volume (Å 3 ) and the electronic population (in green). Some asynaptic basins are not labeled because they hold very little electron density and are thus not meaningful. Values of ELF isosurfaces are given for each picture.

  encourage the reader to visualize the movies in SI. The ELF basins around the region of impact are drained by the +2 charge of the alpha particle. The V(C, C) basin is split into two monosynaptic V(C) basins, one located well above the aromatic plane in the direction of the helium atom, the other situated in opposite direction, i.e. below the plane. Together these two basins hold 3.5 electrons, which is close to the electronic population of the V(C, C) basins before collision (3.23, snapshot a). The intrinsic dipole moments within the monosynaptic V(C) basins are large, indicating strong polarization of the electron cloud caused by the collision. We also note the splitting of the unique V(C, O) basin associated on snapshot a) to the C=O bond into two basins comprising each, in average, 1.21 electrons. This is characteristic of a non-conjugated double C=O bond suggesting a weakening of electron delocalization among the aromatic cycle of guanine. As seen from a close examination of the central C-C region (see movies and Figure4), the original bonding between central C-C atoms progressively reforms within a few hundreds of attoseconds after collision. The V(N) basins are also insightful. Those of nitrogens pertaining to the C=N-C groups belong to the  system as reflected by ELF attractors lying within the aromatic plane. These basins are strongly affected by collision as seen from the large fluctuations of their electronic populations after disruption of the central C=C bond (FigureS11).

Furthermore theses basins are

  occasionally split into sub-basins. The V(N) basins belonging to the  system (CNCCH 3 and CNHC) are situated above and below the aromatic plane. They are similarly impacted by collision since they transiently disappear and later reform. Altogether, these topological patterns indicate that both  and  systems surrounding the impact point are contributing to the redistribution of the delivered energy through the entire molecule (around 67eV, FigureS3). To analyze more in depth the relaxation mechanism we computed the MO occupation numbers and the electronic kinetic energy integrated over atoms (FiguresS4 and S9).

Figure 4 :

 4 Figure 4: Topological analyses of TD-ELF after collision of guanine by a 1MeV  particle travelling along the z direction. a-e:for each TD-ELF topological basins we indicate, in order, the volume (Å 3 ), the electronic population and the intrinsic dipole moment (D, components in brackets). f: superposition of ELF isosurfaces (=0.8) of snapshots a-d. Color code: green: carbons, white: hydrogens, blue: nitrogen, red: oxygen, yellow: helium ( particle) and orange: off nuclei ELF attractors. See also SI for movies.

ACKNOWLEDGMENT

We thank Federica Agostini, David Lauvergnat and Mehran Mostafavi for several stimulating discussions and assistance in the preparation of the manuscript.

AUTHOR INFORMATION

The authors declare no competing financial interests.