% this file is called up by thesis.tex % content in this file will be fed into the main document %: ----------------------- name of chapter ------------------------- %\newcommand{\boldm}[1] {\mathversion{bold}#1\mathversion{normal}} \ifpdf \graphicspath{{4/figures/PNG/}{4/figures/PDF/}{4/figures/}} \else \graphicspath{{4/figures/EPS/}{4/figures/}} \fi \chapter{Dynamical Simulation of Collision-Induced Dissociation} \label{chap:collision} This \textbf{fourth chapter of this thesis} merges two independent studies relating to the dynamical simulation of the collision-induced dissociation of (H$_2$O)$_n$UH$^+$ clusters and pyrene dimer cation Py$_2^+$. The two studies in this chapter share the same methodology to generate the collision trajectories for the collision of argon with (H$_2$O)$_n$UH$^+$ and Py$_2^+$. The collision process of the two studies involves dynamical simulations carried out at a QM/MM level where argon is treated as a polarisable MM particle and the lowest-energy targets (H$_2$O)$_n$UH$^+$ and Py$_2^+$ are treated as SCC-DFTB level. The dynamical simulations performed in these two studies allow to visualise the collision trajectories from which it is possible to analyse in details a number of properties. The theoretical results are compared with the CID experimental results conducted on the same systems, \textit{i.e.} (H$_2$O)$_n$UH$^+$ and Py$_2^+$, by S. Zamith and J.-M. l’Hermite, which facilitates their interpretation and complete the CID experiments. \section{Experimental Methods} \label{exp_cid} The stability of cluster can be investigated from dissociation experiments. Clusters can be dissociated in electric field, magnetic field, high pressure environment, or by heating (such as absorption of photons) or colliding with energetic particles and so on. For instance, the sodium cluster ions and lithium cluster cation were dissociated with a pulsed UV laser source.\cite{Brechignac1989, Brechignac1994} Gaseous hydrated trivalent metal ions were dissociated using blackbody infrared radiative dissociation (BIRD).\cite{Wong2004, Bush2008} The collision between cluster and high or low energetic particles at different pressure also have been investigated. Collisions between the high energetic projectile ions (such as 3 keV Ar$^+$, 22.5 keV He$^{2+}$) and neutral targets were investigated by Gatchell and Holm.\cite{Holm2010, Gatchell2014, Gatchell2017} Collisions between clusters and projectile have been also explored at low collision energy, which allows for the derivation of dissociation energies and the thermal evaporation and stability of clusters. \cite{Boering1992, Wells2005, Zamith2019thermal} By colliding a molecule, or a molecular aggregate, with a non-reactive rare gas atom (neon, argon) or a small molecule such as H$_2$O or N$_2$, it is possible to monitor the parent ions and collision products by use, for instance, of tandem mass spectrometry (MS/MS).\cite{Ma1997, Chowdhury2009} The resulting mass spectra provide a wealth of information about the structure of the parent and product ions from which one can infer, for instance, dissociation mechanisms \cite{Nelson1994, Molina2015} or bond and hydration enthalpies \cite{Carl2007}. The overall process of collisional activation followed by dissociation/fragmentation is commonly referred to as the collision‐induced dissociation (CID) that is also named collisionally activated dissociation (CAD). CID is a mass spectrometry technique to induce dissociation/fragmentation of selected ions in the gas phase, which is one of standard methods for the determination of dissociation/fragmentation pathways. \cite{Sleno2004ion, Wells2005} The CID technique consists of accelerating a given ion into a collision gas thereby the ion getting energy and inducing fragmentation. The produced ionic fragments are then mass analyzed, yielding essentially a mass spectrum.\cite{Cody1982} The CID technique has been applied in different context. Higher-energy C-trap dissociation is a CID technique specific to the orbitrap mass spectrometer in which dissociation/fragmentation occurs outside the trap \cite{Olsen2007higher, Hart2011} Sustained off-resonance irradiation collision-induced dissociation (SORI-CID) is a CID technique used in Fourier transform ion cyclotron resonance mass spectrometry which involves accelerating the ions in cyclotron motion, in a circle inside of an ion trap, in the presence of a collision gas.\cite{Gauthier1991, Laskin2005} %Another technique, multiple-collision induced dissociation (MCID), was used for the study of dissociation energies of singly charged silver cluster, Ag$_n^+$ (2 $\leqslant$ n $\leqslant$ 25) by Kr{\"u}ckeberg and collaborators. \cite{Kruckeberg1999} CID and the dissociated/fragmented ions produced by CID are used for several purposes: First, partial or complete structural determination can be achieved. Second, CID can simply achieve more sensitive and specific detection. By detecting a unique dissociated/fragmented ion, the precursor ion can be detected in the presence of other ions of the same m/z value, mass-to-charge ratio, which reduces the background and increases the limit of detection. CID has been applied to a variety of systems, in particular hydrated atomic ions \cite{Mcquinn2009, Carl2013, Hofstetter2013, Coates2017, Coates2018} and molecular ions \cite{Graul1989, Wei1991, Goebbert2006, Haag2009}. In the second case, it has been used to understand the impact of high-energy radiations on living cells and DNA or RNA \cite{Liu2006, Nguyen2011, Shuck2014}, as well as low-energy collisions on molecules of biological interest \cite{Castrovilli2017, Bera2018}. Theoretical and experimental studies devoted to fragmentation of hydrated molecular aggregates are scarce, \cite{Li1992, Bobbert2002, Liu2006, Bakker2008, Markush2016, Castrovilli2017} although CID has been applied to water clusters containing an atomic ion \cite{Carl2013, Hofstetter2013, Coates2018} and on charged water clusters \cite{Dawson1982, Bakker2008, Mcquinn2009, Zamith2012}. This is a real lack as understanding hydration of molecules and biomolecules is of paramount importance to get insights into their structure, stability, dynamics and reactivity in aqueous medium. In that respect, CID investigations could play an important role in understanding those properties in a environment free from long-range solvent effects but also for different hydration degrees or protonation states. This can be evidenced by the experimental study of Liu \textit{et al.} on the fragmentation of the singly-charged adenosine 5'-monophosphate (AMP$^-$) which shows two different fragmentation channel depending on the solvation state of AMP$^-$.\cite{Liu2006} However, to the best of our knowledge, no modelling was performed to complement these experiments except for a few static calculations.\cite{Carl2013, Hofstetter2013, Coates2018}. %colliding with a high or low energetic particles at high or low pressure,\cite{Kambara1977, Kruckeberg1999, Spasov2000, Holm2010, Gatchell2014, Gatchell2017, Zamith2019thermal} and so on has been explored in experiment. %magnetic field and electric field are usually used in the experiment setup for the dissociation of clusters. %Or write heating, colliding with particles, high pressure, electric field, and magnetic field can induce the dissociation of clusters. %I can not find literature which describe only electric field or magnetic field can lead to dissociation of clusters. You know some ?} Threshold collision-induced dissociation (TCID) method has also been used, for instance to study the fragmentation patterns and to measure the dissociation energies of clusters.\cite{Spasov2000, Armentrout2008} Zamith \textit{et al.} did a CID study of the mass-selected protonated uracil water clusters with water molecules and noble gases, respectively.\cite{Braud2019} In addition, they also reported the TCID study of pyrene cluster cations. \cite{Zamith2020threshold} For these two projects, the single collision event is the predominant process. In this chapter, MD simulations based on a quantum chemical formalism are able to model such complex dissociation mechanism to provide an atomic-scale description for these collisions to explain and complete these experiments. %Mass spectrometry has been the workhorse for studies of collisions involving PAHs, fullerenes and their clusters. %measure the evaporation rates of mass selected pyrene clusters cations as a function of their initial temperature. The dissociation energy of the clusters are extracted from the experimental data using a statistical decay model, namely the Phase Space theory (PST). %\subsection{Principle of TCID} \subsection{Principle of TCID} \label{principleTCID} In usual TCID setups, experiments are done in ion guides, allowing to perform collisions with large mass atoms such as xenon without losing ions by deflection due to the collision. In order to unambiguously determine dissociation energies, one has to take care of a number of experimental parameters. First, the number of collisions should be as low as possible in order to insure single collision conditions. This can be achieved by performing experiments at various pressures and extrapolating results to zero pressure. Second, one has to consider the possible so-called kinetic shifts that can alter the dissociation energy measurement. Indeed, at threshold collision energy, the system under study might not dissociate during the timescale of the experiment. The apparent threshold has, therefore, to be corrected. This is usually done by extrapolating the experimental values using Rice-Ramsperger-Kassel-Marcus (RRKM) dissociation rates.\cite{Klippenstein1992, Baer1996} Third, the initial thermal energy distribution has to be taken into account. Finally, TCID experimental results are usually fitted assuming a given form for the CID cross section, which can be expressed as \cite{Armentrout2008}: \begin{eqnarray} \label{CIDcross} \sigma(E_{col})=(\sigma_{_0} n/E_{col})\sum_{i}{g_{i}} \int_{E_0-E_i}^{E_{col}}[1-e^{-k(\varepsilon + E_i)\tau}] \times (E_{col}-\varepsilon)^{n-1}d\varepsilon \end{eqnarray} where $\sigma_{_0}$ is the collision cross section, $n$ is the energy dependence of the reaction cross section, and $E_{col}$ is the collision energy. The populations $g_i$ of rovibrational states with energies $E_i$ are used to carry out the thermal averaging. The dissociation rate $k$ is usually calculated using RRKM type theories, and $\tau$ is the typical experimental time between the collision and detection. For comparison with experimental curves, eq \ref{CIDcross} is further convolved with the kinetic energy distributions of both the ion and neutral reactants. If one needs to incorporate sequential fragmentation and/or competitive channels, these can also be included.\cite{Rodgers1998, Armentrout2007} In this method, ion guides are not used. Therefore, it needs to simulate the full ion trajectories in order to ensure that ion losses are correctly taken into account. Collisions are, thus, described with a microscopic model rather than with the average curve given by eq \ref{CIDcross}. This approach allows to quite naturally include sequential dissociation and to potentially test energy transfer models. One advantage of the setup resides in the fact that the systems under study are thermalized at low temperature prior to collisions. This implies that averaging over thermal energies of the parent ion plays a minor role, thus leading to minor uncertainties. %For instance, for the largest size under study, n = 6, the average internal energy is 13 meV at 25 K. %\subsection{Experimental setup} \subsection{Experimental Setup} \label{EXPsetup} The experimental setup of Laboratoire Collisions Agr{\'e}gats R{\'e}activi{\'e} (LCAR) by Zamith {\textit et al.} for the collision of protonated uracil water clusters or pyrene dimer cation with noble gas atoms is set up as follows: \figuremacrob{experiment-setup}{Schematic view of the experimental setup. (a) Cluster gas aggregation source. (b) Thermalization chamber. (c) First Wiley–McLaren acceleration stage. (d) Massfilter. (e) Energy focusing. (f) Deceleration. (g) Collision cell. (h) Second Wiley–McLaren acceleration stage. (i) Reflectron. (j) Micro-channel plate detector.} %%%%%%%%%% %from uracil exp paper: %Charged clusters are produced in a gas aggregation source.20 A flow of helium carrier gas and water vapor enters the source in a cell where uracil in solid form (Alpha Aesar, 99\%) is heated at 378 K to get a vapor pressure of about 6.6 $\times$ 10$^{-2}$ mbar.21 A miniature electron gun located at the cell exit allows us to ionize the clusters in the source. The typical kinetic energy of the electrons is 100 eV. The clusters formed in the source at the temperature of liquid nitrogen are carried out of the source by the helium gas flow. The produced clusters then enter the thermalization chamber where tens of thousands collisions between the helium gas and the clusters cool them down to the helium gas temperature. A thermalization temperature of 25 K is obtained by cooling the thermalization chamber with a closed cycle helium cryostat. This low temperature reduces the amount of internal energy prior to collision compared to experiments usually performed at room temperature. For instance, using the harmonic frequencies !i calculated for the neutral uracil,22 its internal energy Eint is estimated by: %%%%%%%%%%% Clusters are produced in a gas aggregation source \cite{Braud2017}(a) and then thermalized (b) at 25 K through thousands of collisions with helium. The experimental setup can be used in two modes. In the first mode, only the first Wiley-McLaren acceleration stage (c) is used to work together with the reflectron (i). Clusters are detected using dual micro-channel plates (MCPs) (j) biased at-10 kV. This allows to perform regular Time of Flight Mass Spectrometry (TOFMS) and to optimize the cluster production. In this mode, the mass filter (d), the electrodes for energy focusing (e) and deceleration (f), and the second Wiley-McLaren acceleration stage (h), are grounded. In the second mode, all the electrodes were used to mass-select the clusters. In order to perform collisions between the mass-selected clusters and the rare gas atoms, precisely delayed high voltage pulses were applied to electrodes (c)-(f). Pulsed high voltages applied to the first Wiley–McLaren electrodes (c) accelerate all the clusters, giving them an average kinetic energy of 622 eV. The applied voltages and the spacing between the electrodes of the Wiley-McLaren are chosen such that, 26 cm downstream, there is a linear relation (to first order) between the position of clusters and their kinetic energy. Using a pulsed high voltage, an electric field is created in this region (e) that compensates this linear kinetic energy dispersion, and all clusters finally have the same kinetic energy within a few electron volts. The time at which this pulsed high voltage is applied determines which cluster size is correctly energy focused. After this kinetic energy focusing, ions are decelerated by a potential barrier (f). At the end of the potential barrier, the potential is shut down in a field free zone and the mass-selected clusters then fly freely through the collision cell (g) up to the second Wile-McLaren acceleration stage (h). Clusters are then mass-analyzed using the reflectron (i) and the MCP detector (j). High voltage is applied on the mass filter (d) when the mass of interest enters the cylinder and shut down before it comes out. This allows us to eliminate part of the neighboring masses. In the experiments of pyrene clusters, the kinetic energy of the clusters in the laboratory frame is varied between 5 eV and 200 eV. For the experiments of protonated uracil water clusters, the kinetic energy of the clusters in the laboratory frame is 100 eV. Kinetic energies of the ions can be easily deduced from experimental parameters. Indeed, since the distances in the apparatus are well-known, measuring, for instance, the time the ions take to travel from the end of the slowing down stage to the second acceleration stage gives the speed of the ions. More precise kinetic energy calibration is obtained by recording the signal of the ions as a function of delays and/or voltages. These curves are then reproduced by simulations to obtain the kinetic energy distribution of the ions. \cite{Chirot2006new} The simplified experiment setup is shown in Figure \ref{exp-setup}. Clusters are produced in a gas aggregation source and thermalized at a temperature of 25 K. Clusters are then mass-selected with a chosen kinetic energy, which collide with argon atoms in a collision cell. The collision products are then analysed by TOFMS. %Further details about the experiment can be found elsewhere \cite{zamith_thermal_2019, braud_gas_2017, chirot_new_2006, zamith_CID_2020}. \figuremacrob{exp-setup}{Schematic of the simplified experimental setup.} \section{Computational Details} \label{Comput_meth} \subsection{SCC-DFTB Potential} \label{DFTBpotential} For the work in this chapter, the SCC-DFTB in combination with the mio-set for the Slater-Koster tables of integrals is applied. \cite{Elstner1998, Porezag1995, Seifert1996, Frenzel2004, Elstner2014, Spiegelman2020} The SCC-DFTB potential for protonated uracil water clusters is shown in section \ref{sec:structure-methods}. DFTB is an efficient tool to perform MD simulations, in particular addressing the evaporation/dissociation processes in various chemical systems. \cite{Simon2017, Korchagina2017, Rapacioli2018, Simon2018} The dynamics simulations of the collision process were performed with a QM/MM scheme \cite{Warshel1976} whose details can be found in the original paper\cite{Cui2001, Iftner2014}. The projectile argon is treated as a polarisable MM particle interacting with the target (protonated uracil water cluster or pyrene dimer cation Py$_2^+$), the latter being treated at the DFTB level. The oscillation problem often appears for dissociated or close to dissociation systems. For the collision trajectories described below, a Fermi distribution (Fermi temperature 2000~K) was applied to avoid oscillation issues during the self-consistent procedure.\cite{Kukk2015} The Fermi distribution allows to recover the continuity in energy and gradients in the case of level crossing. \cite{Kukk2015} It should be mentioned that, in order to keep a low computational cost, no correction has been used to improve the DFTB charge resonance description. However, this charge delocalization issue has been specifically addressed in the case of PAH cation dissociation and it was shown to have a minor effect on the final computed mass spectra \cite{Simon2017}. I also mention that the collision energy is in principle high enough to have electronic excitation in the system, which is taken into account at a crude level by the use of a Fermi temperature. Finally, nuclear quantum effects are not taken into account. Although this may affect the results at very low collision energies, the effect is expected to be small for the experimental collision energies of 7.2 eV and 17.5 eV. Although all these limits should be kept in mind, I would like to emphasize that, recently, the dissociation of PAH molecules has been simulated and a good agreement with experimental results was obtained despite similar crude approximations, namely neglect of non-adiabatic and nuclear quantum effects, improper treatment of charge delocalization and use of a Fermi temperature \cite{Simon2017, Simon2018, Rapacioli2018atomic}. \subsection{Collision Trajectories} \label{makingtrajectories} The preparation for the collisional trajectories for the collision of protonated uracil water clusters or Py$_2^+$ with Ar is the same. The schematic example for the collision of Py$_2^+$ with Ar is shown in Figure \ref{howinputs}. Starting from the optimized Py$_2^+$ geometry \cite{Dontot2019}, a preliminary thermalization run of 200 fs at 25 K (maintained by a Nos\'e-Hoover chain thermostat \cite{Nose1984, Hoover1985}) is performed. Then, the argon atom projectile is introduced in the simulation with a velocity determined to reproduce a given collision energy. The target dimer Py$_2^+$ was positioned at the origin of the simulation referential and randomly rotated to allow all possible impact points on the dimer. The argon atom is initially positioned at x=10, y=$b$ and z=0~\AA{}, with $b$ being the impact parameter. At each center of mass collision energy $E_{col}$, a series of 300 collision trajectories were conducted (the center of mass of the aggregate was kept at position (0, 0, 0)) for each of the 13 $b$ values which are evenly distributed (interval being 0.5 \AA) between 0 and R+0.5 \AA. R refers to the radius of Py$_2^+$. 600 collision trajectories were performed per isomer of protonated uracil water clusters. Trajectory calculations have been performed with a time step of 0.5~fs and a total duration of 15 ~ps and 3~ps for the collision of argon with protonated uracil water clusters and Py$_2^+$, respectively. For the collision of Py$_2^+$ with argon, we have checked that for high collision energies such as 20 and 25 eV, a time step of 0.1 fs does not change significantly our numerical results, which will be shown in section \ref{sec:MDanalysis}. It should be noted that the quaternion was used for rotation process in the generation of initial inputs. %In some cases it may be advisable to use the rigid-body approximation. This approximation allows us to go from $3n-6$ degrees of freedom to $6N-6$, where $n$ and $N$ are the number of atoms and the number of molecules (3 degrees of translation and 3 degrees of rotation) in a system, respectively. The complex quaternion formalism ($\textbf{\textit{q}}$ = $q_0, q_1, q_2, q_3$) was used o describe the orientation of a solid body with respect to Euler angles (($\theta, \phi, \psi$)) formalism. %because the divergence problems appear if we are satisfied with the traditional Euler angles ($\theta, \phi, \psi$) formalism. The quaternions involve an additional degree of freedom (similar to a homothety), which can be offset by using a normalization constraint on the vector \textbf{\textit{q}}: \begin{align} q_0^2+ q_1^2+ q_2^2+ q_3^2=1 \label{vectorq} \end{align} \figuremacrob{howinputs}{Schematic of the generation of the initial inputs.} \\ \subsection{Trajectory Analysis} \label{trajecanylysis} During the results collection, the final snapshot was extracted for each trajectory. For collision of between Py$_2^+$ and argon, when the center of mass of the two pyrene monomer is more than 10 \AA, a dissociation is defined. For the dissociaton definition of protonated uracil water clusters, it is a little more complicated than Py$_2^+$. A fragment is defined as a group of atoms in which the distance of any pair of adjacent atoms is less than 5.0 \AA . The number of hydrogen, nitrogen and oxygen atoms in one fragment is denoted by $k$, $l$ and $m$, respectively. For instance, a fragment characterised by $l=0$ and $k=2m+1$ is a pure water cluster containing the excess proton. Identifying such a fragment at the end of the trajectory means that a neutral uracil fragment exists, otherwise the excess proton is located on a uracil containing fragment. In practice, at each time step, the fragments are identified on the basis of their $k$, $l$, $m$ values, allowing to record their time-dependent evolution. The mass spectrum is built retaining only the fragments containing the excess proton, as only charged fragments are detected in the experiment. The opacity $P(b,E_{col})$, {\textit i.e.} the dissociation probability as a function of impact parameter at a given collision energy is computed by averaging the results over the simulations corresponding to these conditions. The cross sections are then derived from the following formula: \begin{eqnarray} \label{integ} \sigma_{frag}(E_{col}) &=& \int_0^{b_{max}} 2\pi P(b,E_{col})bdb \\ & \simeq& \sum_{i=0}^{b_{max}} \frac{P(b_i,E_{col})+P(b_{i+1},E_{col})}{2}\pi(b_{i+1}^2-b_{i}^2) \nonumber \end{eqnarray} Mean values are computed using the same approach, followed by a division by $\pi b_{max}^2$. When mean values are restricted to trajectories leading to dissociation (noted $-d$) or not (noted $-ud$), additional normalisation by the total number of dissociated or undissociated trajectories is also necessary. \section{\label{sec:collisionwUH}Dynamical Simulation of Collision-Induced Dissociation of Protonated Uracil Water Clusters} \subsection{Introduction} Motivated by the recent CID experiments conducted by Braud \textit{et al.} consisting in (H$_2$O)$_{1-15}$UH$^+$ clusters colliding with an impacting atom or molecule M (M = H$_2$O, D$_2$O, neon, and argon) at a constant center of mass collision energy of 7.2~eV,\cite{Braud2019} the dynamical simulations of the collision between the protonated uracil water clusters (H$_2$O)$_{n=1-7, 11, 12}$UH$^+$ and an argon atom were performed. %First, it appears important to understand the interaction of DNA or RNA basis with water seeing %their relevance in living organisms. They can also be subject to radiation damages which is still a medical challenge and thus %needs to be further investigated. In that context, a number of studies have been conducted on molecules deriving from uracil %or on uracil with only a few water molecules. \cite{Rasmussen2010,Imhoff2007,Abdoul2000,Champeaux2010,Delaunay2014, %Bacchus2009,Kossoski2015} \cite{Maclot2011, Domaracka2012, Markush2016, Castrovilli2017}. Second, a very recent The low collision energy (7.2~eV) only leads to intermolecular bond breaking, without any electronic excitation, rather than intramolecular bond breaking. The branching ratios for different charged fragments were determined in experiment, which allows to deduce the fragmentation cross section for all (H$_2$O)$_{n=1-15}$UH$^+$ species and the location of the excess proton after collision: on a uracil containing cluster or on a pure water cluster. This allows to determine the proportion of neutral uracil loss (corresponding to cases where the excess proton is located on pure water clusters) as a function of the number $n$ of water molecules. A sharp increase of neutral uracil loss was observed for $n$ = 5-6 (2.8\% and 25.0\% for n = 4 and 7, respectively). % Those experiment were complemented by theoretical calculations that aim at characterizing the lowest-energy isomers of (H$_2$O)$_{n}$UH$^+$ ($n$ = 1-7, 11, 12) clusters (see section \ref{structureUH}), which %They show that (i) For $n$ = 1-2, the uracil is protonated; (ii) For $n$ = 3-4, the excess proton is still on the uracil but has a tendency to be displaced towards adjacent water molecules; (iii) When $n$ is larger than 4, the excess proton is transferred to the water molecules. shows that the location of the proton after collision recorded in the CID experiment is determined by its position in the lowest-energy parent isomer. In other words, a shattering mechanism occurs after collision. Despite these findings, static calculations can not provide a full picture for the fragmentation process and some issues are still not properly understood: (i) What is the main path of the fragmentation mechanisms? (ii) What are the fragments after collision? (iii) How does the proportion of fragments change according to time? (iv) Is the proportion of neutral uracil molecules loss only determined by the nature of the lowest-energy isomers? %(v) What is the impact of nuclear quantum effects (NQE) for such process that occurs at very low temperature? To answer these questions, this simulations present a complete MD study of the fragmentation process for (H$_2$O)$_{n=1-7, 11, 12}$UH$^+$ aggregates colliding with an argon atom. Section~\ref{resul_disc} discusses the statistical convergence of collision trajectories, theoretical time-dependent proportion of fragments, proportion of neutral uracil loss, total fragmentation cross sections and mass spectra of fragments bearing the excess proton. These data are compared to available experimental results in order to discuss in details dissociation mechanism as a function of $n$. The main outcomes are summarized in section~\ref{Concl}. \subsection{Results and Discussion} \label{resul_disc} %\subsubsection{The proper number of impacting orientation} \subsubsection{Statistical Convergence} \label{convergence} In order to ensure that the statistical convergence is reached in the collision trajectories, initial conditions have to reproduce all possible collision orientation with good statistics. The procedure to generate a set of collision trajectories is described in section \ref{makingtrajectories}. As a visual proof, pictures a, b and c in Figure ~\ref{3b-sphere} represent 200, 400 and 600 random argon orientations with impact parameter being 0 for cluster (H$_2$O)$_3$UH$^+$, respectively. In these pictures, (H$_2$O)$_3$UH$^+$ is fixed and all initial positions for argon are orientated which leads to distribution maps of the initial positions of argon with respect to fixed (H$_2$O)$_{3}$UH$^+$. It is worth pointing out that in the collision trajectories, argon is fixed and uracil is rotated. Picture d in Figure~\ref{3b-sphere} presents 200 random argon orientations with impact parameter being 0.0 and 6.0, respectively. The similar pictures for cluster (H$_2$O)$_{12}$UH$^+$ are displayed in Figure \ref{12f-sphere}. These Figures demonstrate that the more collision trajectories are performed, the more colliding opportunities of argon at all possible orientations are obtained. \figuremacro{3b-sphere}{Schematic representation of random argon orientations for the collision with the second lowest-energy isomer of cluster (H$_2$O)$_3$UH$^+$. 200 (a), 400 (b) and 600 (c) random argon orientations are generated with impact parameter being 0. ~200 orientations are generated with impact parameter being 0 and 6 (d), respectively.} \\ \figuremacro{12f-sphere}{Schematic representation of random argon orientations for the collision with the second lowest-energy isomer of cluster (H$_2$O)$_{12}$UH$^+$. 200 (a), 400 (b) and 600 (c) random argon orientations are generated with impact parameter being 0. ~200 orientations are generated with impact parameter value being 0 and 7 (d), respectively.} \\ In addition, to confirm that statistical convergence is reached for the properties discussed in sections \ref{time}, \ref{small}, \ref{large}, and \ref{mass-spectra}, Tables~\ref{tab:converge-1w-5w} and \ref{tab:converge-6w-12w} present the $P_{NUL}$ involved from the following formula \ref{PNUL} \begin{eqnarray} \label{PNUL} %&&P_{NUL}(E_{col})= \int_0^{b_{max}} 2\pi \frac{N_{NUL}(b,E_{col})}{N_{frag}(b,E_{col})}bdb \nonumber \\ & P_{NUL}(E_{col})&=& \int_0^{b_{max}} N_{NUL}(b,E_{col}) 2\pi bdb / \int_0^{b_{max}} N_{frag}(b,E_{col})2\pi bdb \nonumber \\ & \simeq& \frac{\sum\limits_{i=0}^{b_{max}}\frac{1}{2}(N_{NUL}(b_i,E_{col})+N_{NUL}(b_{i+1},E_{col}))\pi(b_{i+1}^2-b_{i}^2)} {\sum\limits_{i=0}^{b_{max}}\frac{1}{2}(N_{frag}(b_i,E_{col})+N_{frag}(b_{i+1},E_{col}))\pi(b_{i+1}^2-b_{i}^2)} \nonumber \\ \end{eqnarray} and $\sigma_{frag}$ of two isomers (the first lowest energy isomer and the one whose $P_{NUL}$ fits best to the experiment results (in bold)) of each cluster (H$_2$O)$_{n=1-7, 11, 12}$UH$^+$ obtained from 200, 400, and 600 random argon orientations per impact parameter value. Whatever the considered isomer, the three $P_{NUL}$ and $\sigma_{frag}$ values from 200, 400, and 600 random argon orientations are very close. Indeed, the largest difference is observed for isomer 7a which has $P_{NUL}$ values of 29.5 and 31.3~\% for 200 and 600 random orientations, respectively. This demonstrate that even for 200 initial random orientations, simulation are close to statistical convergence. In the present study, all results discussed in the main text were obtained with 600 initial random argon orientations per impact parameter value which ensures statistical convergence of the results independently of cluster size. %\begin{eqnarray} %\label{PNUL} % P_{NUL}(E_{col}) &=& \int_0^{b_{max}} 2\pi \frac{N_{NUL}}{N_{frag}}(b,E_{col})bdb \\ &% % \simeq& % \frac{\sum\limits_{i=0}^{b_{max}}(N_{NUL}(b_i,E_{col})+N_{NUL}(b_{i+1},E_{col}))\frac{(b_{i+1}^2-b_{i}^2)}{b_{max}^2}} % {\sum\limits_{i=0}^{b_{max}}(N_{frag}(b_i,E_{col})+N_{frag}(b_{i+1},E_{col}))\frac{(b_{i+1}^2-b_{i}^2)}{b_{max}^2}} \nonumber %\end{eqnarray} \begin{table} \begin{center} \caption{The proportions of $P_{NUL}$ and $\sigma_{frag}$ of first lowest-energy isomer and the isomer whose $P_{NUL}$ fits the experiment (in bold) of (H$_2$O)$_{1-5}$UH$^+$ with simulations of 200, 400, and 600 as initial conditions.} \label{tab:converge-1w-5w} \begin{tabular}{|c|c|c|c|} \hline \textbf{Cluster} & \textbf{Simu} & \textbf{\boldm{$P_{NUL}$} (\%)} & \textbf{\boldm{$\sigma_{frag}$ (Å$^2$)}}\\ \hline \rowcolor{lightgray} \bf{1a} & 200 & 0.1 & 28.4 \\ \rowcolor{lightgray} \bf{1a} & 400 & 0.1 & 28.3 \\ \rowcolor{lightgray} \bf{1a} & 600 & 0.2 & 28.9 \\ 1b & 200 & 0.2 & 26.3 \\ 1b & 400 & 0.1 & 25.7 \\ 1b & 600 & 0.1 & 25.9 \\ \rowcolor{lightgray} 2a & 200 & 0.0 & 35.9 \\ \rowcolor{lightgray} 2a & 400 & 0.0 & 36.5 \\ \rowcolor{lightgray} 2a & 600 & 0.0 & 36.3 \\ \bf{2b} & 200 & 0.0 & 34.7 \\ \bf{2b} & 400 & 0.1 & 34.8 \\ \bf{2b} & 600 & 0.1 & 34.9 \\ \rowcolor{lightgray} 3a & 200 & 5.4 & 37.4 \\ \rowcolor{lightgray} 3a & 400 & 5.2 & 36.2 \\ \rowcolor{lightgray} 3a & 600 & 5.7 & 36.3 \\ \bf{3b} & 200 & 0.0 & 41.2 \\ \bf{3b} & 400 & 0.0 & 41.5 \\ \bf{3b} & 600 & 0.0 & 41.9 \\ \rowcolor{lightgray} 4a & 200 & 26.9 & 40.1 \\ \rowcolor{lightgray} 4a & 400 & 28.2 & 40.3 \\ \rowcolor{lightgray} 4a & 600 & 29.4 & 40.1 \\ \bf{4b} & 200 & 2.7 & 45.3 \\ \bf{4b} & 400 & 2.6 & 45.6 \\ \bf{4b} & 600 & 2.6 & 45.2 \\ \rowcolor{lightgray} 5a & 200 & 37.2 & 45.7 \\ \rowcolor{lightgray} 5a & 400 & 37.8 & 46.1 \\ \rowcolor{lightgray} 5a & 600 & 38.2 & 46.6 \\ \bf{5d} & 200 & 0.1 & 47.3 \\ \bf{5d} & 400 & 0.1 & 47.3 \\ \bf{5d} & 600 & 0.1 & 47.5 \\ \hline \end{tabular} \end{center} \end{table} \begin{table} \begin{center} \caption{The proportions of $P_{NUL}$ and $\sigma_{frag}$ of first lowest-energy isomer and the isomer whose $P_{NUL}$ fits the experiment (in bold) of (H$_2$O)$_{6, 7, 11, 12}$UH$^+$ with simulations of 200, 400, and 600 as initial conditions.} \label{tab:converge-6w-12w} \begin{tabular}{|c|c|c|c|} \hline \textbf{Cluster} & \textbf{Simu} & \textbf{\boldm{$P_{NUL}$} (\%)} & \textbf{\boldm{$\sigma_{frag}$ (Å$^2$)}}\\ \hline \rowcolor{lightgray} 6a & 200 & 38.0 & 46.6 \\ \rowcolor{lightgray} 6a & 400 & 38.0 & 45.6 \\ \rowcolor{lightgray} 6a & 600 & 39.3 & 45.8 \\ \bf{6f} & 200 & 18.9 & 54.2 \\ \bf{6f} & 400 & 19.0 & 55.2 \\ \bf{6f} & 600 & 18.5 & 55.0 \\ \rowcolor{lightgray} 7a & 200 & 29.5 & 54.8 \\ \rowcolor{lightgray} 7a & 400 & 31.3 & 53.4 \\ \rowcolor{lightgray} 7a & 600 & 31.3 & 53.4 \\ \bf{7d} & 200 & 22.6 & 55.3 \\ \bf{7d} & 400 & 22.9 & 54.3 \\ \bf{7d} & 600 & 23.0 & 54.0 \\ \rowcolor{lightgray} 11a & 200 & 26.7 & 53.8 \\ \rowcolor{lightgray} 11a & 400 & 28.2 & 53.5 \\ \rowcolor{lightgray} 11a & 600 & 28.3 & 52.9 \\ \bf{11d} & 200 & 14.5 & 55.2 \\ \bf{11d} & 400 & 15.4 & 56.1 \\ \bf{11d} & 600 & 15.6 & 56.5 \\ \rowcolor{lightgray} 12a & 200 & 8.0 & 59.2 \\ \rowcolor{lightgray} 12a & 400 & 7.5 & 60.5 \\ \rowcolor{lightgray} 12a & 600 & 7.6 & 60.2 \\ \bf{12c} & 200 & 10.4 & 55.3 \\ \bf{12c} & 400 & 10.8 & 55.8 \\ \bf{12c} & 600 & 10.8 & 55.4 \\ \hline \end{tabular} \end{center} \end{table} \subsection{Time-Dependent Proportion of Fragments} \label{time} The time-dependent proportion of each fragment was extracted from collision trajectories. To illustrate the change in behavior resulting from the difference in cluster size, Figures~\ref{proporEachFrag-1a2a}-\ref{proporEachFrag-7d12c-zoom} display the time-dependent proportion of fragments obtained from the dissociation of the low-lying energy isomers of (H$_2$O)$_{1-7, 11, 12}$UH$^+$. I will discuss the time-dependent proportion of one small cluster (7a) and one big cluster (12a) in detail as an example. For the sake of clarity, only the fragments displaying significant proportion, higher than 0.035 and 0.015 for 7a and 12a, respectively are considered in Figure \ref{proporEachFrag-7a12a-zoom}. This corresponds to the eight and ten most prominent fragments for 7a and 12a, respectively. From Figure~\ref{proporEachFrag-7a12a-zoom}, it is clear that for both aggregates, the proportion of H$_2$O has the sharpest increase after collision and then stay almost constant as a function of time. For 7a, $\sim$3~ps after collision, the proportion of almost all fragments does not change any more. Before that, the proportion of the (H$_2$O)$_6$UH$^+$ fragment increases first and then decreases, which indicates a sequential dissociation of water molecules. For 12a, proportion of (H$_2$O)$_{11}$UH$^+$ and (H$_2$O)$_{10}$UH$^+$ fragments displays a sharp increase quickly after collision which is then followed by a fast decrease, and finally it keeps a minute decrease up to the end of the simulations. The decrease of proportion of (H$_2$O)$_{10}$UH$^+$ and (H$_2$O)$_{11}$UH$^+$, and the increase of proportion of (H$_2$O)$_{6}$UH$^+$, (H$_2$O)$_{7}$UH$^+$ and (H$_2$O)$_{8}$UH$^+$ indicate sequential dissociation after collision is occurring. It is worth noting that, in contrast to 7a, the proportions of the main fragments of 12a do not tend to be a constant at the end of the simulations. This implies that, for this large aggregate, structural rearrangements are more likely to occur prior to complete dissociation. Proportions of the main fragments of clusters 7d and 12c shown in Figure \ref{proporEachFrag-7d12c-zoom} display similar behavior as for 7a and 12a. As a first conclusion, Figure~\ref{proporEachFrag-7a12a-zoom} suggests that clusters with 7 water molecules experience a direct dissociation mechanism as was hypothesised by Braud \textit{et al.}.\cite{Braud2019} A similar conclusion can be drawn for smaller cluster sizes as supported by Figures \ref{proporEachFrag-1a2a}-\ref{proporEachFrag-5a6a-zoom}. In contrast, cluster with 11 (see Figure \ref{proporEachFrag-11a-zoom}) and 12 water molecules shows a behavior compatible with a certain amount of statistical dissociation, namely a long-time evolution that allows structural rearrangements. These important observations can now be refined by looking at more detailed properties. \figuremacrob{proporEachFrag-1a2a}{Time-dependent proportions of the main fragments obtained from the dissociation of the lowest-energy isomers of (H$_2$O)$_1$UH$^+$ (left) and (H$_2$O)$_{2}$UH$^+$ (right).} \\ \figuremacro{proporEachFrag-3a4a-zoom}{Time-dependent proportions of the main fragments obtained from the dissociation of the lowest-energy isomers of (H$_2$O)$_3$UH$^+$ (left) and (H$_2$O)$_{4}$UH$^+$ (right). Bottom panels correspond to a zoom over the lower proportions.} \\ \figuremacro{proporEachFrag-5a6a-zoom}{Time-dependent proportions of the main fragments obtained from the dissociation of the lowest-energy isomers of (H$_2$O)$_5$UH$^+$ (left) and (H$_2$O)$_{6}$UH$^+$ (right). Bottom panels correspond to a zoom over the lower proportions.} \\ \figuremacro{proporEachFrag-11a-zoom}{Time-dependent proportions of the main fragments obtained from the dissociation of the lowest-energy isomer of (H$_2$O)$_{11}$UH$^+$. Right panel corresponds to a zoom over the lower proportions.} \\ \figuremacro{proporEachFrag-7a12a-zoom}{Time-dependent proportions of the main fragments obtained from the dissociation of the lowest-energy isomers of (H$_2$O)$_7$UH$^+$ (left) and (H$_2$O)$_{12}$UH$^+$ (right). Bottom panels correspond to a zoom over the lower proportions.} \\ \figuremacro{proporEachFrag-7d12c-zoom}{Time-dependent proportions of the main fragments obtained from the dissociation of the the third lowest-energy isomer of (H$_2$O)$_7$UH$^+$ (left) and the third lowest-energy isomer (H$_2$O)$_{12}$UH$^+$ (right). Bottom panels correspond to a zoom over the lower proportions.} \\ \subsection{Proportion of Neutral Uracil Loss and Total Fragmentation Cross Sections for Small Clusters} \label{small} In order to get more insights in the fragmentation, molecular dynamics trajectories were analysed in terms of proportion of neutral uracil loss ($P_{NUL}$) defined in section \ref{convergence} and total fragmentation cross sections ($\sigma_{frag}$) defined in section~\ref{trajecanylysis}. These two properties are also accessible from experiments. Another property extracted from the MD simulations, but not accessible from experiment, is the proportion of protonated uracil ($P_{PU}$) which is equal to the ratio of the number of simulations leading to a protonated uracil molecule over the number of simulations leading to a fragment containing the uracil and the excess proton. In order to correlate the outcome of the collision and the structure of the aggregate undergoing the collision, all considered low-energy isomers are characterized by there relative energies ($E_{rel.}$) and the location of the excess proton (LEP). For the latter, three distinct configurations were considered: The excess proton is bounded to the uracil molecule (noted U-H); The excess proton is bounded to a water molecule that is adjacent to an oxygen atom of the uracil molecule (noted W-H-U); The excess proton is bounded to a water molecule that is separated by at least one other water molecule from the uracil molecule (noted W-H). All these data are gathered in Table~\ref{tab:full} and we first discuss the behavior of the small species (H$_2$O)$_{1-7}$UH$^+$. \begin{table*} \begin{center} \caption{Relative energy $E_{rel.}$ (in kcal.mol$^{-1}$) at the MP2/Def2TZVP level, LEP, $P_{PU}$ (in \%), $P_{NUL}$ (in \%), $\sigma_{frag}$ (in \AA$^2$) of the considered low-energy isomers of (H$_2$O)$_{1-7, 11, 12}$UH$^+$ clusters. Isomers which $P_{NUL}$ fit best to the experimental value are indicated in bold. $P_{{NUL}_{exp}}$ and $\sigma_{{frag}_{exp}}$ are the experimental values for $P_{NUL}$ and $\sigma_{frag}$, respectively. For (H$_2$O)$_{12}$UH$^+$, experimental values were obtained for collision with Ne, whereas all other theoretical and experimental data are for collision with Ar.}\label{tab:full} \begin{tabular}{|c|c|c|c|c|c|l|c|} \hline \textbf{Isomers} & $E_{rel.}$ & LEP & $P_{PU}$ & $P_{NUL}$ & $P_{{NUL}_{exp}}$ & $\sigma_{frag}$ & $\sigma_{{frag}_{exp}}$ \\ \hline \rowcolor{lightgray} \bf{1a} & 0.0 & U-H & 100 & 0.2 & & 28.9 & \\ \rowcolor{lightgray} 1b & 0.7 & U-H & 100 & 0.1 &\multirow{-2}*{0.9} & 25.9 & \multirow{-2}*{12.3} \\ 2a & 0.0 & U-H & 100 & 0.0 & & 36.3 & \\ \bf{2b} & 0.2 &U-H & 100 & 0.1 &\multirow{-2}*{0.4} &34.9 & \multirow{-2}*{22.8}\\ \rowcolor{lightgray} 3a & 0.0 & U-H & 100 & 5.7 & & 36.3 & \\ \rowcolor{lightgray} \bf{3b} & 0.3 & U-H & 100 & 0.0 &\multirow{-2}*{1.7} & 41.9 & \multirow{-2}*{31.2} \\ 4a & 0.0 & W-H-U & 98.0 & 29.4 & & 40.1 & \\ \bf{4b} & 0.9 &U-H & 99.7 & 2.6 &\multirow{-2}*{2.8} &45.2 & \multirow{-2}*{43.4}\\ \rowcolor{lightgray} 5a & 0.0 &W-H & 78.5 & 46.6 && 38.2 & \\ \rowcolor{lightgray} 5b & 0.3 &W-H-U & 89.0 & 28.5 && 38.7 & \\ \rowcolor{lightgray} 5c & 2.0 &W-H-U & 87.8 & 27.1 && 44.6 & \\ \rowcolor{lightgray} \bf{5d} & 2.4 &U-H & 100 & 0.1 &\multirow{-4}*{7.5}& 47.5 & \multirow{-4}*{48.0}\\ 6a & 0.0 &W-H & 44.1 & 39.3 & & 45.8 & \\ 6b & 0.2 &W-H & 43.5 & 33.8 & & 58.6 & \\ 6c & 0.3 &W-H & 46.4 & 36.6 & & 46.1 & \\ 6d & 0.9 &W-H & 64.6 & 34.7 & & 42.6 & \\ 6e & 2.5 &W-H & 45.9 & 34.9 & & 50.5 & \\ \bf{6f} & 2.7 &W-H-U & 76.2 & 18.5 &\multirow{-6}*{18.0}& 55.0 & \multirow{-6}*{54.3} \\ \rowcolor{lightgray} 7a & 0.0 &W-H & 28.2 & 31.3 && 53.4 & \\ \rowcolor{lightgray} 7b & 0.3 & W-H-U & 52.4 & 21.4 && 51.7 &\\ \rowcolor{lightgray} 7c & 0.3 & W-H & 41.3 & 31.1 && 49.5 &\\ \rowcolor{lightgray} \bf{7d} & 0.8 & W-H-U & 40.9 & 23.0 &\multirow{-4}*{25.0}& 54.0 & \multirow{-4}*{59.7} \\ 11a & 0.0 &W-H & 4.6 & 28.3 & & 52.9 & \\ 11b & 1.4 &W-H & 3.2 & 28.5 & & 54.7 & \\ 11c & 1.5 &W-H & 4.2 & 22.8 & & 55.2 & \\ \bf{11d}& 1.9 &W-H & 6.8 & 15.6 &\multirow{-4}*{11.8} & 56.5 & \multirow{-4}*{63.8} \\ 11e & 1.9 &W-H & 5.4 & 22.7 & & 52.6 & \\ 11f & 2.3 &W-H & 7.9 & 24.3 & & 52.0 & \\ % 11g & 2.5 &W-H & 3.7 & 29.8 &\multirow{-6}*{11.8}& 55.0 & \multirow{-6}*{63.8} \\ \rowcolor{lightgray} 12a & 0.0 & W-H & 6.7 & 7.6 && 60.2 & \\ \rowcolor{lightgray} 12b & 0.6 &W-H & 34.0 & 22.4 && 52.2 & \\ \rowcolor{lightgray} \bf{12c} & 0.7 &W-H & 48.7 & 10.8 && 55.4 & \\ \rowcolor{lightgray} 12d & 1.3 &W-H-U & 5.4 & 9.7 && 54.3 & \\ \rowcolor{lightgray} 12e & 1.8 &W-H-U & 67.5 & 6.0 && 54.2 & \\ \rowcolor{lightgray} 12f & 2.4 &W-H-U & 55.0 & 17.1 &\multirow{-6}*{12.2}& 54.1 & \multirow{-6}*{77.0} \\ \hline \end{tabular} \end{center} \end{table*} Various information can be inferred from these properties. Firstly, one observes a general increase of $\sigma_{frag}$ as a function of cluster size with values ranging from 25.9~\AA$^2$~for isomer 1b to 60.2~\AA$^2$~for isomer 12a. Interestingly, only slight variations of $\sigma_{frag}$ are observed for different isomers of the same aggregate. In contrast, $P_{NUL}$ is much more sensitive to the nature of the considered isomers, in particular when these isomers display different LEP values. For instance, $P_{NUL}$ is 46.6 ~\% for 5a (W-H) while it is only 0.1~\% for 5d (U-H). More interestingly, there seems to exist a strong correlation between $P_{NUL}$ and LEP. Indeed, $P_{NUL}$ values below 1.0~\% are characterized by an excess proton initially bounded to uracil (U-H type). This suggests that when uracil is protonated, probability for deprotonation after collision is very low and thus $P_{NUL}$ is close to 0\%. $P_{NUL}$ values between 9.7 and 29.4~\% are obtained from W-H-U configurations while larger $P_{NUL}$ values, above 31.1~\%, arise from W-H configurations in clusters (H$_2$O)$_{5-7}$UH$^+$. This demonstrates that, from the excess proton point of view, the outcome of the collision is highly sensitive to the nature of the isomer undergoing the collision as hypothesised by Braud \textit{ al.} \cite{Braud2019} This important finding can be of help to determine which isomer, or set of isomers, is likely to undergo collision by comparing experimental and theoretical $P_{NUL}$ as this is not necessarily the lowest-energy isomer as discussed below. For (H$_2$O)$_{1-2}$UH$^+$, the theoretical and experimental $P_{NUL}$ values, close to zero, are in good agreement regardless of the considered isomer. For (H$_2$O)$_3$UH$^+$, the experimental $P_{NUL}$ is 1.7~\% which is well reproduced by both isomers 3a and 3b although 3b is the one closer to the experimental value, 0.0~\% against 5.7~\% for 3a. This was expected as they are very close in energy, only 0.3 kcal.mol$^{-1}$ difference, and in structure, as displayed in Figure~\ref{fig-1a-3b}, both being of U-H type structure. Consequently, in the experiment, each one of them could be at the origin of the experimental signal. (H$_2$O)$_4$UH$^+$ behaves differently. The two low-energy isomers, 4a and 4b, display very different $P_{NUL}$ values, 29.4 and 2.6~\%, respectively. The experimental value is 2.8~\% which suggests that 4b, although slightly higher in energy by 0.9 kcal.mol$^{-1}$, is the isomer prevailing during the collision process. The difference in behavior can be explained by the U-H configuration of 4b, in which the excess proton is bounded to the uracil, whereas in 4a, it is bounded to a water molecule adjacent to uracil (see Figure~\ref{fig-4a-5d}). The case of (H$_2$O)$_5$UH$^+$ is more complex as this is the first species displaying the three types of LEP configuration among its four lowest-energy isomers as can be seen on Figure~\ref{fig-4a-5d}. This implies very different $P_{NUL}$ values: 46.6~\% for 5a, 28.5 and 27.1~\% for 5b and 5c, respectively, while it is only 0.1~\% for 5d. The experimental $P_{NUL}$ value for (H$_2$O)$_5$UH$^+$ is still relatively low, 7.5~\%, which suggests that a U-H type structure prevails during the collision process. Although 5d is 2.4~kcal.mol$^{-1}$ higher in energy than 5a, this isomer is thus expected to undergo the collision. \figuremacro{fig-1a-3b}{Selected low-energy configurations of (H$_2$O)$_{1-3}$UH$^+$. Relative energies at the MP2/Def2TZVP level are in kcal.mol$^{-1}$.} \figuremacro{fig-4a-5d}{Selected low-energy configurations of (H$_2$O)$_{4-5}$UH$^+$. Relative energies at the MP2/Def2TZVP level are in kcal.mol$^{-1}$.} (H$_2$O)$_6$UH$^+$ and (H$_2$O)$_7$UH$^+$ are the first two aggregates for which no low-energy isomer belongs to the U-H type structure. As a consequence, in contrast to smaller species, the theoretical $P_{NUL}$ values are all higher than 15~\%. This is in line with the experimental values which display a net increase at $n$ = 6. Isomers 6a, 6b, 6c, 6d, and 6e (see Figure~\ref{fig-6a-6f}) are all W-H type structures which leads to $P_{NUL}$ values almost twice higher than the experimental one. Consequently, as for (H$_2$O)$_5$UH$^+$, one can assume that the isomer of (H$_2$O)$_6$UH$^+$ undergoing the collision is more likely to be a W-H-U type structure although it is higher in relative energy. Isomer 6f can be such a candidate as it displays of $P_{NUL}$ value of 18.5\% which is in agreement with the experimental value, 18.0\%. Due to its increasing size, (H$_2$O)$_6$UH$^+$ displays W-H configurations with the excess proton at various distances from the recombining oxygen. Indeed, in 6a, 6c and 6d this distance is 1.774, 1.745 and 1.804~\AA, while in 6b, 6e and 6f, it is shorter: 1.660, 1.614, and 1.494~\AA, respectively. However, no net correlation is observed between this distance and the value of $P_{NUL}$: 39.3, 33.8, 36.6, 34.7, 34.9 and 18.5\% for 6a, 6b, 6c, 6d, 6e and 6f, respectively. In particular, the behavior of 6e is striking. It has almost the same relative energy as 6f and they are structurally similar (see Figure~\ref{fig-6a-6f}) but display different $P_{NUL}$ values. This suggests that, for $n$ larger than 5, the ability of the water molecule network to stabilise the excess proton, \textit{i.e} to promote or prevent its diffusion toward the uracil molecule, starts to be competitive with the configuration type of the isomer. In 6e, the excess proton is in a configuration close to the Zundel ion which may explain its high $P_{NUL}$ value as compared to 6f. For (H$_2$O)$_7$UH$^+$, a W-U-H type configuration is also expected to fit best the experimental result. And indeed 7d, a W-H-U type structure, which is only 0.8~kcal.mol$^{-1}$ above the lowest-energy isomer (see Figure~\ref{fig-7a-7d}), has a $P_{NUL}$ value of 22.9~\% as compared to 25.0~\% experimentally. Isomers 7a and 7c have a W-H configuration and their $P_{NUL}$ values (31.3 and 31.1~\%, respectively) are higher than the ones of 7b and 7d which have a W-H-U configuration. Finally, it is worth noting that even when the excess proton is initially bounded to a water molecule, \textit{i.e.} when a W-H type structure is considered, the maximum $P_{NUL}$ that has been obtained is only 46.6~\%. This demonstrates that for small aggregates such as (H$_2$O)$_{5-7}$UH$^+$ ((H$_2$O)$_{1-4}$UH$^+$ do not display low-energy W-H type structures), dissociation mainly lead to protonated uracil containing fragments. This is in line with the experimental results. Analysis of $P_{PU}$ values also show that uracil is protonated in a significant amount of these protonated uracil containing fragments. $P_{PU}$ has a clear tendency to decrease with cluster size, but can be quite high even for W-H type structures, for instance 5a, 6d and 7c in Table~\ref{tab:full}. This demonstrates that upon collision, the excess proton is likely to transfer to uracil on a rather short time scale. \figuremacro{fig-6a-6f}{Selected low-energy configurations of (H$_2$O)$_{6}$UH$^+$. Relative energies at the MP2/Def2TZVP level are in kcal.mol$^{-1}$.} \figuremacro{fig-7a-7d}{Selected low-energy configurations of (H$_2$O)$_{7}$UH$^+$. Relative energies at the MP2/Def2TZVP level are in kcal.mol$^{-1}$.} The clusters discussed above are characterized by complex potential energy surfaces characterized by several low-energy isomers, with relative energies that can be lower than 1~kcal.mol$^{-1}$, and which get more complex as the number of water molecules increases. Consequently, the exact energetic ordering between the low-energy isomers can not be precisely known as this is below chemical accuracy and we thus can not claim here to have found the lowest-energy structure of each aggregate, or the isomer undergoing the collision. Nevertheless, what we show is that $P_{NUL}$ is mainly determined by the initial position of the proton in the isomer undergoing the collision. Consequently, for the collision energy and the range of cluster size we have considered, the structure of the aggregate undergoing the collision plays a key role in determining the dissociation process and collision outcomes much more than energetics. This is consistent with the analysis of the time-dependent proportion of fragments which suggests a direct dissociation mechanism. This is further highlighted on Figure~\ref{neutralUloss-Ne-Ar}, which presents the experimental $P_{NUL}$ for collision with Ar and Ne, respectively as a function of $n$ and the corresponding theoretical values obtained from the lowest-energy isomers as well as from the isomers for which $P_{NUL}$ matches best to the experimental data. As can be seen, a very good agreement can be obtained with the experimental data by considering a specific set of isomers. Interestingly, if a similar plot is drawn for $\sigma_{frag}$ considering the same isomers (see Figure~\ref{cross-section-Ne-Ar}), a good agreement with the experimental data and much better than $\sigma_{geo}$ (calculated from formula \ref{cross-section-geo})is also obtained with the two sets of isomers which confirms the weaker dependence upon isomer of $\sigma_{frag}$. \figuremacro{neutralUloss-Ne-Ar}{Theoretical (green and blue lines) and experimental (red line) $P_{NUL}$ values for the (H$_2$O)$_{1-7, 11, 12}$UH$^+$ clusters. Theory 1 (green line) is obtained from the isomers which $P_{NUL}$ matches best to the experimental data while Theory 2 (blue line) is obtained from lowest-energy isomers.} \\ \figuremacro{cross-section-Ne-Ar}{Theoretical (green and blue lines) and experimental (red line) $\sigma_{frag}$ values for the (H$_2$O)$_{1-7, 11, 12}$UH$^+$ clusters. Theory 1 (green line) is obtained from the isomers which $P_{NUL}$ matches best to the experimental data while Theory 2 (blue line) is obtained from lowest-energy isomers.} \\ \subsection{Behaviour at Larger Sizes, the Cases of (H$_2$O)$_{11, 12}$UH$^+$} \label{large} In the experiments conducted by Braud \textit{et al.},\cite{Braud2019} $P_{NUL}$ starts to decrease at $n$=8. This decrease is not consistent with the above argument of a direct dissociation mechanism and larger species more likely characterized by W-H and W-H-U type structures. This apparent discrepancy motivated us to extend the present study to a larger cluster, namely (H$_2$O)$_{11, 12}$UH$^+$. For (H$_2$O)$_{12}$UH$^+$, the only available experimental data is for collisions with Ne instead of Ar, although for the same center of mass collision energy. As shown in Figure~\ref{neutralUloss-Ne-Ar}, experimental $P_{NUL}$ values for Ne or Ar, although not equal, display similar trend. In the following, I thus discuss the experimental data of (H$_2$O)$_{12}$UH$^+$ colliding with Ne. For cluster (H$_2$O)$_{1-7, 11}$UH$^+$, keep discussing the experimental data from colliding with Ar. The behaviour for (H$_2$O)$_{11}$UH$^+$ and (H$_2$O)$_{12}$UH$^+$ is rather different when looking at detailed properties. Indeed, for (H$_2$O)$_{11}$UH$^+$, $P_{NUL}$ values for isomers 11a, 11b, 11c, 11e, and 11f are very similar as they range from 22.7 to 29.8~\%. For 11d, $P_{NUL}$ equal 15.6~\% which fits best to the experiment. These $P_{NUL}$ values are lower than those of (H$_2$O)$_{6}$UH$^+$, as observed experimentally, and in the same range as (H$_2$O)$_{7}$UH$^+$. All isomers display a W-H type configuration as seen in Figure~\ref{fig-11a-f}. $P_{PU}$ is very small for all (H$_2$O)$_{11}$UH$^+$ isomers which shows that on the time scale of the simulations, protonation of uracil hardly occurs. For (H$_2$O)$_{12}$UH$^+$, 12c isomer, which has a W-H type configuration (see Figure~\ref{fig-12a-f}), has a $P_{NUL}$ value which fits best to the experiment, 10.8~\% against 12.2~\%, while isomer 12a, also a W-H type configuration (see Figure~\ref{fig-12a-f}), has a $P_{NUL}$ value equal to 7.6~\%. Overall, $P_{NUL}$ values calculated for (H$_2$O)$_{12}$UH$^+$ isomers are lower than those of (H$_2$O)$_{6}$UH$^+$ and (H$_2$O)$_{7}$UH$^+$, which is in line with the experiment. The main difference with the (H$_2$O)$_{1-7}$UH$^+$ aggregates is that no clear relation exist between the $P_{NUL}$ value and the initial localisation of the excess proton. Indeed, 12a, 12b and 12c are all W-H type configurations but with $P_{NUL}$ values ranging from 7.6 to 22.4~\%. The same is observed for 12d, 12e and 12f although they are all W-H-U type configurations. Similarly, no difference in behaviour is obtained between W-H and W-H-U type configurations. This can be explained by assuming that the dissociation mechanism in (H$_2$O)$_{12}$UH$^+$ involves some amount of structural rearrangement that softens the impact of the isomer undergoing the collision. Indeed, as (H$_2$O)$_{12}$UH$^+$ has more degrees of freedom, it can more easily accommodate the kinetic energy transferred by the colliding atom prior to dissociation which thus takes place on a longer time scale. This excess of internal energy thus fosters structural rearrangements, in particular proton transfers toward the uracil, explaining the smaller $P_{NUL}$ value for (H$_2$O)$_{12}$UH$^+$. This is in full agreement with the conclusions obtained in section~\ref{time} from Figures~\ref{proporEachFrag-11a-zoom},\ref{proporEachFrag-7a12a-zoom} and \ref{proporEachFrag-7d12c-zoom}. To further support this conclusion, we conducted 200 MD simulations in the micro-canonical ensemble in which the whole kinetic energy of Ar was randomly distributed in all the vibrational modes of isomer 12c by drawing initial velocities in a 1185~K Boltzmann distribution. Among them, 166 simulations display dissociation with one or two water molecules dissociating from the main cluster. No neutral uracil loss is observed. To conclude, although the present simulations are too short to assert that (H$_2$O)$_{12}$UH$^+$ undergoes a statistical dissociation mechanism, they clearly show that a direct mechanism is not sufficient to account for the theoretical and experimental results. Consequently, structural rearrangements are very likely to occur prior to dissociation and the experimental results for $P_{NUL}$ and $\sigma_{frag}$ values can not result from a single (H$_2$O)$_{12}$UH$^+$ isomer. In contrast, similarities in both $P_{NUL}$ and $P_{PU}$ values for all considered (H$_2$O)$_{11}$UH$^+$ isomers, as well as $P_{NUL}$ values close to (H$_2$O)$_{7}$UH$^+$ ones, do not evidence structural rearrangements in this species although they could be present. \figuremacro{fig-11a-f}{Selected low-energy configurations of (H$_2$O)$_{11}$UH$^+$. Relative energies at the MP2/Def2TZVP level are in kcal.mol$^{-1}$.} \\ \figuremacro{fig-12a-f}{Selected low-energy configurations of (H$_2$O)$_{12}$UH$^+$. Relative energies at the MP2/Def2TZVP level are in kcal.mol$^{-1}$.} \\ \subsection{Mass Spectra of Fragments with Excess Proton} \label{mass-spectra} In this section, in order to analyse collision products in more details, the branching ratios of the different fragments containing the excess proton were extracted from the collision simulations of clusters (H$_2$O)$_{1-7, 11, 12}$UH$^+$ and compared with the experimental ones shaped as mass spectra.\cite{Braud2019} For each cluster size, only simulations corresponding to the isomer which $P_{NUL}$ value fits best to the experiment were considered (1a, 2b, 3b, 4b, 5d, 6f, 7d, 11d, 12c). The results are presented in Figures~\ref{MS-BR-1w-4w-Ne-Ar-branch}, \ref{MS-BR-5w-11w-Ne-Ar-branch} and \ref{MS-BR-12w-Ne-branch}. For cluster (H$_2$O)$_{12}$UH$^+$, there is no experimental data for collision with argon. For (H$_2$O)$_{12}$UH$^+$, experimental results were obtained for collision with neon. From the experimental results for argon and neon for 1a, 2b, 3b, 4b, 5d, 6f, 7d, and 11d, it shows the branch ratios for collision with argon and neon are close and have the same trend. So it should be reasonable to compare the simulated branch ratios of 12c with the ones of experimental data from the collision of (H$_2$O)$_{12}$UH$^+$ with argon. \figuremacro{MS-BR-1w-4w-Ne-Ar-branch}{Simulated mass spectra (positive area) of the charged fragments after 15~ps simulation time (fragments (H$_2$O)$_n$H$^+$ in red and (H$_2$O)$_n$UH$^+$ in blue for argon; (H$_2$O)$_n$H$^+$ in pink and (H$_2$O)$_n$UH$^+$ in green for neon) from isomers (a) 1a, (b) 2b, (c) 3b, (d) 4b. The counterparts in experiment are plotted (negative area).} \\ \figuremacro{MS-BR-5w-11w-Ne-Ar-branch}{Simulated mass spectra (positive area) of the charged fragments after 15~ps simulation time (fragments (H$_2$O)$_n$H$^+$ in red and (H$_2$O)$_n$UH$^+$ in blue) from isomers (e) 5d, (f) 6f, (g) 7d, and (h) 11d. The counterparts in experiment are plotted (negative area).} \\ \figuremacro{MS-BR-12w-Ne-branch}{Simulated mass spectra (positive area) of the charged fragments after 15~ps simulation time (fragments (H$_2$O)$_n$H$^+$ in red and (H$_2$O)$_n$UH$^+$ in blue) from isomers 12c. The counterparts in experiment obtained for collision with neon are plotted in negative area (H$_2$O)$_n$H$^+$ in pink and (H$_2$O)$_n$UH$^+$ in green).} \\ Overall, the experimental and theoretical spectra present the same general trends: (i) The mass spectra present a broad distribution of sizes, without prominence of a particular peak; (ii) All the spectra are dominated by the heaviest protonated uracil containing fragment (loss of a single water molecule) with the exception of the simulated mass spectrum for 2b; (iii) Fragments containing protonated uracil prevail over pure protonated water fragments, as already observed from the $P_{NUL}$ values provided in Table~\ref{tab:full}; (iv) Pure protonated water fragments only appear for the largest cluster sizes. Indeed, although very minor contributions are observed in both the simulated and experimental spectra for parent clusters with $n$=3-5, significant contributions of these species only appear when the parent cluster contains at least 6 water molecules. A more detailed discussion of the simulated and experimental mass spectra will be made as follows. For (H$_2$O)$_{2}$UH$^+$, fragments (H$_2$O)UH$^+$ and UH$^+$ are observed in both experiment and theory although their relative ratio is different. For (H$_2$O)$_{3}$UH$^+$, the simulated and experimental spectra agree quite well with a dominant peak for (H$_2$O)$_{2}$UH$^+$. For (H$_2$O)$_{4}$UH$^+$, in both experimental and theoretical spectra, the peak intensity of the fragments containing protonated uracil increases with the number of water molecules. For (H$_2$O)$_{5}$UH$^+$ and (H$_2$O)$_{6}$UH$^+$, this is also the case except for the UH$+$ fragment which is overestimated when compared to the experimental result in Figure \ref{MS-BR-5w-11w-Ne-Ar-branch} (e). For (H$_2$O)$_{7, 11, 12}$UH$^+$, the intensities for the heaviest fragments are overestimated. From Figures \ref{MS-BR-1w-4w-Ne-Ar-branch}, \ref{MS-BR-5w-11w-Ne-Ar-branch} and \ref{MS-BR-12w-Ne-branch}, it is clear that the smaller the cluster (except Figure \ref{MS-BR-1w-4w-Ne-Ar-branch} (b)) is, the better the agreement between the simulated and experimental branching ratios is. This trend indicates that for small clusters, \textit{i.e.} for $n=1-6$, short simulation time is enough to capture the full dissociation pattern, in other words, the dissociation mechanism is direct with no noticeable contribution of long term evolution. However, for larger clusters, starting at $n$=7, owing to the larger number of degrees of freedom, short simulation time does not capture the full dissociation pattern, \textit{i.e.} long term statistical dissociation is more likely to play a noticeable role. This is fully in line with the conclusions obtained in section~\ref{large} for (H$_2$O)$_{12}$UH$^+$ and refine the interpretation given in section~\ref{small} for (H$_2$O)$_{7}$UH$^+$. This also shows that although the data presented in section~\ref{large} for (H$_2$O)$_{11}$UH$^+$ do no evidence the contribution of structural re-arrangements on the short time scale, they are very likely to occur as in (H$_2$O)$_{12}$UH$^+$. One has to keep in mind that modeling the complete duration of the experiment (up to $\mu$s) is out of reach with MD/SCC-DFTB simulations. In this work, the simulation time was 15~ps, for all cluster sizes. Large fragments such as (H$_2$O)$_{6-12}$UH$^+$ may lose more water molecules if long enough simulation time were available, as suggested from the time dependent evolution of selected trajectories in section~\ref{time}. To certify this, the total energy of (H$_2$O)$_6$UH$^+$ fragments at SCC-DFTB level is calculated originating from the dissociation of (H$_2$O)$_7$UH$^+$ (7d) from all the 1421 trajectories producing fragment (H$_2$O)$_6$UH$^+$ over the total 600 $\times$ 15 trajectories. Then the energies of the lowest-energy isomer of (H$_2$O)$_5$H$^+$ and H$_2$O at SCC-DFTB level are subtracted. The deduced relative energies $\Delta E$ are reported in Table~\ref{tab:fragenergy} for five cases. When $\Delta E$ is greater than zero, it is possible for the (H$_2$O)$_6$UH$^+$ fragment to lose a water molecule. The percentage of $\Delta E$ being positive in all the trajectories leading to fragment (H$_2$O)$_6$UH$^+$ is 53.0~\%, which indicates that many (H$_2$O)$_6$UH$^+$ fragments have still the potential to lose one more water molecule after the end of the simulation. \begin{table} \begin{center} \caption{Energies of different (H$_2$O)$_6$UH$^+$ fragments selected from the dissociation of 7d at SCC-DFTB level, and the lowest energies (H$_2$O)$_5$UH$^+$ and (H$_2$O) at SCC-DFTB level. The relative energy $\Delta E$ = $E_{(H_2O)_6UH^+}$ -($E_{(H_2O)_5UH^+}$ +$ E_{H_2O}$). All energies here are given in eV.} \label{tab:fragenergy} \begin{tabular}{|c|c|c|c|} \hline \textbf{\boldm{$E_{(H_2O)_6UH^+}$}} & \textbf{\boldm{$E_{(H_2O)_5UH^+}$}} & \textbf{\boldm{$E_{H_2O}$}} & \textbf{\boldm{$\Delta E$}}\\ \hline -44.310 & -40.312 & -4.057 & 1.605 \\ \hline -44.322 & -40.312 & -4.057 & 1.279\\ \hline -44.307 & -40.312 & -4.057 & 1.687 \\ \hline -44.344 & -40.312 & -4.057 & 0.680 \\ \hline -44.373 & -40.312 & -4.057 & -0.109\\ \hline \end{tabular} \end{center} \end{table} \subsection{Conclusions about CID of (H$_2$O)$_{n}$UH$^+$} \label{Concl} Collision-induced dissociation of protonated uracil water clusters (H$_2$O)$_{1-7, 11, 12}$UH$^+$ at constant center of mass collision energy has been investigated by molecular dynamics simulations using the SCC-DFTB method. The very good agreement between the simulated and measured $P_{NUL}$ and $\sigma_{frag}$ as well as branching ratios indicate that the essence of the dissociation induced by collisions is well captured by the simulations. The $P_{NUL}$ values from the different isomers of the (H$_2$O)$_{1-7}$UH$^+$ cluster show that the localization of the excess proton after dissociation is strongly determined by the initial configuration of the isomer undergoing the collision. This suggests that (H$_2$O)$_{1-7}$UH$^+$ aggregates primarily engage a direct dissociation path after collision that takes place on a very short time scale, \textit{i.e.} lower than 15~ps. More strikingly, in most cases, the proposed lowest-energy isomer does not lead to the best fit to the experiment. However, the relative energy between the lowest-energy isomers and the isomers best fitting to the experiment is less than 1.0 kcal.mol$^{-1}$ for (H$_2$O)$_{1-4, 7}$UH$^+$ clusters and less than 2.7 kcal.mol$^{-1}$ for (H$_2$O)$_{5,6}$UH$^+$ clusters. This is in line with the strong sensitivity of the collision outcome with the nature of the isomer undergoing the collision. This even suggests that the LEP can help in determining the main characteristic of the isomer involved in the collision. For (H$_2$O)$_{11, 12}$UH$^+$, these conclusions do not apply any more which shows that significant structural rearrangements occur after collision. This is confirmed by the time-dependent proportion of fragments which continue to vary even at 15~ps for (H$_2$O)$_{11, 12}$UH$^+$ whereas it is almost flat for (H$_2$O)$_{1-7}$UH$^+$. Analysis of the fragment branching ratios helps in clarifying these points. Indeed, for the smallest clusters, (H$_2$O)$_{1-5}$UH$^+$, the short simulation time well reproduces the corresponding experimental results which is in line with a direct mechanism. In contrast, for (H$_2$O)$_{6-7}$UH$^+$, although $P_{NUL}$ is well reproduced by the simulations, the experimental and theoretical branching ratios differ which show that more time is needed to properly describe the dissociation. For (H$_2$O)$_{11, 12}$UH$^+$, neither theoretical nor experimental branching ratios and $P_{NUL}$ are in agreement which is a strong indication that a significant contribution of structural rearrangements occur; This suggests that a contribution of a statistical mechanism is more likely to occur for larger species such as (H$_2$O)$_{11, 12}$UH$^+$. This work demonstrates that explicit molecular dynamics simulations performed at a quantum chemical level can provide a wealth of information about collision-induced mechanism in molecular clusters, in particular, hydrated molecular species. Such simulations thus represent a key tool to complement CID experiments and hope the present study will motivate similar computational studies on future CID experiments of hydrated molecular aggregates. %In a near future, we think it would be of great interest to pursue this study by looking at the influence of collision energy, both lower or higher, on the dissociation mechanism as a function of the cluster size. Furthermore, inclusion of nuclear quantum effects in the simulations could also help to increase the accuracy of the model and improve the comparison with the experiments. \section{Dynamical Simulation of Collision-Induced Dissociation for Pyrene Dimer Cation} \subsection{Introduction} PAH clusters have been investigated in several scientific fields. In combustion science, the role of PAH clusters in combustion processes is still under debate, in particular they might or not be the intermediate systems in the growth of soot particles. \cite{Chung2011, Saggese2015, Eaves2015, Mao2017, Wang2018} In atmospheric and environmental science, PAHs are known as the pollutants, which is harmful to human health. For instance, the carcinogenic PAHs associated with particulate matter in air pollution has showed clear evidence of genotoxic effects, such as DNA adducts, chromosome aberrations. \cite{Kyrtopoulos2001, Farmer2003} In new energy resources field, for the understanding of the properties of organic crystal or the design of new organic solar cell devices, PAH stacks are investigated as the prototypes.\cite{Aumaitre2019} In astrophysics, PAHs species are believed to be ubiquitous and abundant in the interstellar medium because of their compact and stable structure. \cite{Tielens2008} The PAH clusters are important contributors to the diffuse interstellar bands and UV-visible absorption bands. PAH clusters have been proposed to be the origin of a series of infrared emission bands, which are ubiquitous in the Universe. \cite{Leger1984,Allamandola1985} The broadening of these bands in regions protected from the star's UV flux suggests the following scenario: PAHs are trapped in clusters in UV-protected regions and photo-evaporated by star's UV photons in the so-called photodissociated region. \cite{Rapacioli2005, Berne2008} For all these topics, it is necessary to make a better understanding of the fundamental properties of PAH clusters. The crucial quantities are the stability, molecular growth processes, dissociation energies and their evolution with PAH charge, species, cluster size. The investigation of PAH clusters has been performed in experiment. Many studies focused on the investigation of structural properties of these clusters at the most stable geometrical configurations \cite{Eschenbach1998, Schmidt2006, Goulart2017, Wang2018, Lei2019}. Their energetic properties such as ionisation potentials have been recorded \cite{Joblin2017} as well as their spectral properties \cite{Roser2015,Lemmens2019}. clusters may evaporate, breaking the PAH units themselves or leading to chemical reactivity between the different units, which shows the role of PAH clusters in the growth of PAHs themselves. If free flying PAHs are possibly from the evaporation of larger clusters, then this calls for more experimental data on their dissociation properties. The evolution of PAH clusters has been explored from experiments following the evaporation after absorption of UV photons, collision with low or high energetic particles or in a high-pressure environment. \cite{Schmidt2006, Holm2010, Gatchell2015, Joblin2017, Gatchell2017, Zamith2019thermal} The range of collision energies considered experimentally is quite large, ranging from eV to high energy collision at a few keV. Low energy collision experiments allow for the derivation of dissociation energies \cite{Zamith2019thermal} whereas the oligomerization of PAHs within the cluster induced by high energy collisions \cite{Delaunay2015} or photoabsorption \cite{Zhen2018} suggests the possible role of clusters in the interstellar PAHs growth process \cite{Chen2018}. The quantitative data from experiments of PAH clusters are still rather limited, which motivates the modeling studies of them. In the calculation of PAH clusters, the size of the systems limits the use of {\it ab initio} wave function methods to the investigation of properties of the smallest clusters, namely dimers \cite{Piacenza2005, Birer2015}, whereas larger clusters can be addressed either at the DFT level or with more semi-empirical schemes \cite{Zhao2008truhlar, Rapacioli2009corr, Mao2017, Bowal2019}. Many of these studies, focused on structural properties, evidence a stacking growth process in agreement with experimental results. In addition, IR properties were also reported at the DFT level \cite{Ricca2013}. Most of the theoretical studies involve neutral clusters, mostly due to the fact that treating charge resonance process in ions is a challenging task for DFT based methods \cite{Grafenstein2009}. The singly charged PAH clusters are more stable than their neutral counterparts due to charge resonance stabilization.\cite{Rapacioli2009} Cationic PAH clusters are expected to be abundant in the photo-dissociation regions because the ionization energy of the PAH cluster is lower than that of the isolated PAH, which leads to the efficient formation of cationic PAH clusters. In addition, the ionized PAH clusters are easier to control, so it is more important to study them. It should be mentioned the recent studies computing ionisation potentials \cite{Joblin2017} as well as structural \cite{Dontot2019} and spectral (electronic \cite{Dontot2016} and vibrational \cite{Dontot2020}) properties of cations, performed with an original model combining DFTB \cite{Porezag1995, Seifert1996, Elstner1998, Spiegelman2020} with a configuration interaction scheme\cite{Rapacioli2011}. With respect to these studies, very few is known about the dynamical aspects of PAH clusters carrying internal energy. High energy collisions of PAH clusters with energetic ions have been simulated by Gatchel {\textit et al.} \cite{Gatchell2016, Gatchell2016knockout} at the semi-empirical and DFTB levels. Recently experiments at lower collision energies were performed by Zamith \textit{et al.} from LCAR \cite{Zamith2020threshold} (the principle of this experiment and the experimental setup were shown in sections \ref{principleTCID} and \ref{EXPsetup}), which were analysed by treating statistically the dissociation after collision energy deposition. Namely, the dissociation rate of pyrene clusters has been computed using phase space theory (PST)\cite{Zamith2019thermal}. A fair agreement with experimental results was obtained concerning the collision energy dependence of the dissociation cross section. However, the employed model failed at reproducing in details the shape of the peaks in the time-of-flight (TOF) spectra. In this section, it is aimed at extending the description of such low energy collision processes (less than several tens of eV) combining a dynamical simulations to describe the fast processes in addition to the statistical theory to address dissociation at longer timescales. With this approach, (i) good agreement between simulated and experimental mass spectra will be shown, thus validating the model, (ii) dissociation cross sections as a function of the collision energy is derived, (iii) the kinetic energy partition between dissociative and non-dissociative modes will be discussed and (iv) the energy transfer efficiency between intra and intermolecular modes will also be discussed. % \ref{sec:compapp}, \subsection{Calculation of Energies} In the analysis, we will discuss the kinetic energy contributors, applying the following decomposition of the total kinetic energy $E^k_{tot}$ of the dimer: \begin{align} \label{Eparti} E^k_{tot}&=E^k_{Ar} + E^k_{td} + E^k_{Py^1} + E^k_{Py^2} + E^k_{Re} \nonumber \\ E^k_{tot} &=\frac{1}{2}\sum_{i=1}^{52}m_i(\Vec{v}_i)^2 \nonumber \\ E^k_{Ar} &= \frac{1}{2}m_{Ar}\Vec{v}_{Ar}^2 \nonumber \\ E^k_{td} & =\frac{1}{2}m_{Py_2}\Vec{v}_t^2(Py_2) \\ E^k_{Re} &= \frac{1}{2} \frac{m_{Py^1}m_{Py^2}}{m_{Py^1}+m_{Py^2}}(\Vec{v}_t(Py^2)-\Vec{v}_t(Py^1))^2 \nonumber \\ E^k_{Py^n} &= \frac{1}{2}\sum_{i=1}^{{26}}m_i^n (\Vec{v}_i^n-\Vec{v}_t(Py^n))^2 \nonumber %\label{Eparti} \end{align} In these equations and in the following, $Py_2$ refers to the pyrene dimer (possibly dissociated) whereas $Py^1$ and $Py^2$ refer to the first and second monomers, respectively. $E^k_{tot}$ can be also calculated from the masses $m_i^n$ and velocities $\Vec{v}_i$ of its atoms. $E^k_{Ar}$ refers to the kinetic energy of the argon (with mass $m_{Ar}$ and velocity $\Vec{v}_{Ar}$). $E^k_{td}$ is the translation kinetic energy of the dimer (with mass $m_{Py_2}$ and velocity $\Vec{v}_t(Py_2)$). $E^k_{Re}$ is the relative kinetic energy of the two pyrene monomers, computed from their masses of $m_{Py^1}=m_{Py^2}$ and monomer global translation velocities $\Vec{v}_t(Py^{n=1,2})$. $E^k_{Py^{n}}$ is the rovibrational kinetic energy of the monomer $n$ computed from the masses and velocities of its atoms ($m_i^n$ and $\Vec{v}_i^n$, respectively). The intramolecular vibrational kinetic energy ($E^k_{intra^{n}}$) of monomer $n$ obtained after removing the contributions associated to the monomer translation and rotation modes is calculated as follows: \begin{align} E^k_{intra^n}&= \frac{1}{2}\sum_{i=1}^{{26}}m_i^n (\Vec{v}_i^n-\Vec{v}_t(Py^n)-\Vec{v}_{ir}^n)^2 \label{Eintra} \end{align} where $\Vec{v}_{ir}^n$ is the velocity of atom $i$ associated to the monomer global rotation. In addition, the dimer intermolecular kinetic energy ($E^k_{inter}$) is calculates as follows: \begin{align} E^k_{inter}&=E^k_{tot}-E^k_{Ar}-E^k_{td}-E^k_r-E^k_{intra^1}-E^k_{intra^2} \label{Einter} \end{align} where $E^k_r$ refers to the rotation kinetic energy of the dimer. $\Vec{v}_{ir}^n$ and $E^k_r$ are calculated using the following formulas. \begin{align} \Vec{L}(Py^n) &= \sum_{i=1}^{26}m_i^n (\Vec{r}_i^n-\Vec{r}_{_{CM}}(Py^n)) \times (\Vec{v}_i^n-\Vec{v}_t(Py^{n})) \nonumber \\ \Vec{L}(Py_2) &= \sum_{i=1}^{52} m_i (\Vec{r}_i-\Vec{r}_{_{CM}}(Py_2)) \times (\Vec{v}_i-\Vec{v}_t(Py_{2})) \nonumber \\ I &= mr^2 \nonumber \\ \Vec{\omega} &= [I]^{-1} \times \Vec{L} \nonumber \\ \Vec{v}_{ir}^n &= \Vec{\omega}(Py^n) \times (\Vec{r}_i^n - \Vec{r}_{_{CM}}(Py^n)) \\ E^k_{tot^n} &=E^k_{r^n} + E^k_{td^n} + E^k_{intra^n} \nonumber \\ E^k_{tot^n} &=\frac{1}{2}\sum_{i=1}^{26}m_i^n(\Vec{v}_i^n)^2 \nonumber \\ E^k_{td^n} &=\frac{1}{2}m_{Py^n}\Vec{v}_t^2(Py^n) \nonumber \\ E^k_{r} &= \frac{1}{2} \Vec{\omega}(Py_2) \times [I](Py_2) \times\Vec{\omega}(Py_2) = \frac{1}{2} \Vec{L}(Py_2) \times [I]^{-1}(Py_2) \times \Vec{L}(Py_2) \nonumber \\ \label{Erotation} \end{align} $E^k_{tot^n}$ is the total kinetic energy of monomer $n$. $E^k_{td^n}$ is the translation kinetic energy of pyrene monomer $n$. $\Vec{L}(Py^n)$ is the angular momentum of pyrene monomer $n$. $\Vec{L}(Py_2)$ is the angular momentum of the dimer. $[I]$ refers to the moment of inertia tensor. $[I]^{-1}$ is the inverse of $[I]$. $\Vec{\omega}$ is the angular velocity. $\Vec{r}_{i}^n$ and $\Vec{r}_{_{CM}}(Py^n)$ denote the coordinates of atom $i$ and center of mass of dimer of monomer $n$, respectively. $\Vec{r}_{i}$ and $\Vec{r}_{_{CM}}(Py_2)$ and denote the coordinates of atom $i$ and center of mass of dimer, respectively. From the endpoint of the simulation, we can also compute the total energy transferred towards internal rovibrational modes of the pyrene dimer as: \begin{equation} \Delta E_{int}^{Py_2} = E^{k,0}_{Ar} - E^k_{Ar} - E^k_{td} \end{equation} where $E^{k,0}_{Ar}$ is the initial argon kinetic energy whereas $E^k_{Ar}$ and $E^k_{td}$ correspond to kinetic energies at the end of the MD simulation. In the case of dissociated dimers at the end of the simulations, we can also deduce the energy deposited in the rovibrational modes of the monomers as: \begin{equation} \Delta E_{int}^{Py^1+Py^2} = E^{k,0}_{Ar} - E^k_{Ar} - E^k_{td} - E^k_{Re} \end{equation} \subsection{Simulation of the Experimental TOFMS} The experimental TOFMS are reproduced by simulating the ion trajectories through the experimental setup in the presence of the electric fields. These are calculated by solving numerically the Laplace equation. Equations of motion are integrated using the fourth order Runge-Kutta method with adaptive step size. The occurrence of collision or dissociation is decided at each time step of the ion trajectory based on the collision and dissociation probabilities. In the work of Zamith {\textit et al.}, \cite{Zamith2020threshold} the energy transfer was treated upon collision by using the Line of Center model (LOC) \cite{Levine1987}. In the LOC model, the transferred energy is the kinetic energy along the line of centers. Evaporation rates were then estimated using PST, in which only statistical dissociation to be possible after energy deposition in the cluster by collision was conside. Although this approach, which will be referred to as PST in the following, has been proved to be able to satisfactorily reproduce CID cross section experiments,\cite{Zamith2020threshold} it fails to reproduce in details the shape and position of the fragment peaks in the TOFMS, as will be shown in section \ref{sec:MS}. In order to better reproduce the position and peak shapes, the MD and PST methods were combined. The outputs of the MD simulations were used to treat the collisions in the ion trajectories. At each time step the probability for a collision is evaluated. The principle of MD+PST is displayed in Figure \ref{MDPST}. \figuremacrob{MDPST}{Principle of MD+PST.} One MD trajectory (with proper weighting of the $b$ values) was randomly picked from all outputs of MD simulations at a given collision energy. Then two cases have to be considered. First, if the dissociation occurred during the picked MD calculation (short time dissociation), then we use the MD final velocities of the fragments to further calculate the ion trajectories. On the other hand, if the pyrene dimer is still intact at the end of the picked MD calculation, then we update the dimer velocity and use the collision energy transfer $\Delta E_{int}^{Py_2}$ deduced from the MD calculation to increase the internal energy of the cluster. The dissociation rate resulting from this new internal energy is then evaluated using PST. In the latter case, if dissociation occurs (long time dissociation), the relative velocities of the fragment are evaluated using the PST outcome. The whole process of MD+PST is performed many times for different trajectories to ensure the reliability of the final obtained data. For each time, the TOFMS of Py$^+$ or Py$_2^+$ is updated. %Then we can obtain the update TOFMS histogram of the monomer cation Py$^+$ from the short and long time dissociation. Here we emphasize that, due to the short time scale of the MD calculations (3 ps), only direct dissociation can be captured by the MD simulations. Therefore, one has to evaluate the probability of dissociation at longer time scales after the energy deposition by collision. This is done here by considering that at longer time scales, dissociation occurs statistically and is treated by using PST. \subsection{\label{sec:results}Results and Discussion} In the following, we will discuss the dissociation at short, experimental and infinite timescales. The first two ones correspond to dissociation occurring during the MD simulation only or with the MD+PST model. The dissociation at infinite time accounts for all MD trajectories where the amount of energy transferred to the internal dimer rovibrational modes $\Delta E_{int}^{Py_2}$ is larger than the dissociation energy of 1.08~eV (value from references \cite{Dontot2019, Zamith2020threshold}). It can be regarded as the dissociation occurring after an infinite time neglecting any cooling processes, such as thermal collisions or photon emissions. \subsubsection{\label{sec:MS}TOFMS Comparison} An example of TOFMS is given in Figure~\ref{expTOF}. Figure~\ref{expTOF}(a) is centered around the intact parent mass (Py$_2^+$) whereas in (b) is displayed the region around the fragment peak (Py$^+$). In Figure~\ref{expTOF}(a) are displayed three curves corresponding to the experimental one and the results of the two simulations (PST and MD+PST) for the parent ion. One can see that the peak shape and position are properly reproduced using the simulations; therefore, the essential of the ion propagation is captured by the simulations. Although some of the detected parent ions have undergone a collision without dissociation, no difference is seen in the parent peak since the collision rate is kept very small. \figuremacrob{expTOF}{Normalized time of flight mass spectra of the parent pyrene dimer cation (a), and the pyrene fragment Py$^+$ (b) resulting from the collision of Py$_2^+$ with argon at a center of mass collision energy of 17.5~eV. The black line is for the experimental result whereas red and green curves are the MD+PST and PST model results. The blue curve is the PST subcontribution of the MD+PST model.} %\figuremacro{2}{} In Figure~\ref{expTOF}(b), the experimental result is compared to the PST and MD+PST simulations. Clearly, the PST based simulation fails to reproduce both the position and the shape of the peak. On the other hand, a much better agreement is found when using the output of the MD+PST simulations. This agreement is a good indication that this scheme captures the essence of the pyrene dimer cation dissociation induced by argon collisions at this collision energy. Actually, in this scheme, the largest contribution to the TOFMS results from dimers dissociating on short timescales, $i.e.$ during the MD simulation. The remaining contribution, $i.e.$ resulting from dimers dissociating at longer timescales and computed from the second step PST calculation, is minor and represented in blue in Figure \ref{expTOF}(b). \subsubsection{\label{sec:MDanalysis}Molecular Dynamics Analysis} \textbf{Description of selected trajectories} A first qualitative description of the collision processes can be obtained from the analysis of some arbitrarily selected MD trajectories. Figure \ref{collisions} (top and bottom) reports some snapshots extracted from two trajectories with the same collision energy (17.5~eV) and impact parameter (3.5~\AA). Only the top one leads to the Py$_2^+$ dissociation. During the results collection, we extract the final snapshot for each trajectory and consider that the dimer is dissociated if the distance between the two monomers molecular mass centers is larger than 10~\AA. Figures \ref{collisions}-1/1* represent the system after its preliminary thermalization, when the argon atom introduced in the simulation with its initial velocity. Figures \ref{collisions}-2/2* and 3/3* represent the beginning and end of the collision. From these points, the two trajectories show different behaviors. For the top trajectory in Figure \ref{collisions}, snapshot 5 corresponds to the step where the two pyrene monomers start to go away from each other. After this, the intermolecular distance continues to increase further in snapshot 6. For the bottom trajectory in Figure \ref{collisions}, Figures \ref{collisions}-5* and 6* correspond to the middle and ending snapshots of the simulation, respectively. The snapshots 4*,~5* and 6* show the process of energy redistribution within the clusters. In particular, the soft modes associated to global deformation of the molecular planes appear to be excited. \figuremacrob{collisions}{Snapshots for two different molecular dynamics trajectories. Top and bottom: trajectories with impact parameter of 3.5~\AA{} and a collision energy of 17.5~eV, leading to dissociation and non-dissociation (top and bottom, respectively).} From these two particular cases, it can be seen that the evolution of the trajectory either toward a dissociation or a redistribution of the transferred energy strongly depends on the process of energy transfer during the collision. In the top trajectory in Figure \ref{collisions}, the argon atom is pushing the two monomers far away from each other, $i.e.$ the transferred energy is mostly localised in an intermolecular dissociative mode. On the opposite, in the bottom trajectory in Figure \ref{collisions}, the collision mostly involves an intramolecular soft vibrational mode. The transferred energy is then redistributed over all the other modes. The statistical distribution of this energy is then hardly favorable to the dissociation due to the large number of intramolecular modes (72 per pyrene) with respect to the 6 intermolecular modes, only 3 of them (1 breathing and 2 parallel displaced modes) being dissociative modes. The amount of transferred energy is also a major ingredient for the fate of the cluster. Depending on the details of the collision such as impact parameter or cluster orientation, very different amounts of energy can be transferred. This is illustrated in Figure~\ref{distriPerc-Etf-175eV-d-bin03} where the distribution of transferred energy $\Delta E_{int}^{Py_2}$ restricted to trajectories that would dissociate after infinite time, is plotted for simulations at the experimental collision energy of 17.5~eV. This distribution could hardly be guessed without a dynamical description of the collision at the atomic level. Indeed, a simpler model such as the LOC model (used in the pure PST approach) would lead to a constant distribution between the binding energy and the maximum collision energy as shown in Figure \ref{distriPerc-Etf-175eV-d-bin03}. In the distribution resulting from MD simulations, lower transferred energies are favored with respect to the distribution extracted from the LOC model. All these effects are intrinsically taken into account in the MD simulations on the opposite to the pure PST model, explaining the better agreement of the MD+PST scheme with the experimental results. % 1/17.5=x*(17.5-1.08) \figuremacro{distriPerc-Etf-175eV-d-bin03}{Distribution of transferred energy in rovibrational modes $\Delta E_{int}^{Py_2}$ for trajectories leading to dissociation at the end of MD (center of mass collision energy of 17.5~eV). The dashed line shows the distribution of transferred energy used in the LOC model.} Finally, we note that the pyrene monomers remained intact (no fragmentation) up to collision energies of 25~eV. The snapshots of a fragmentation trajectory at collision energy of 27.5~eV are shown in Figure \ref{fragmentation}. It can be seen that the pyrene molecule impacted by the argon undergoes an opening of an aromatic cycle and the loss of two hydrogen atoms, leaving as a H$_2$ molecule. As the study of monomer's fragmentation is beyond the scope of the present paper, we will focus on trajectories with collision energies below this fragmentation threshold energy in the following. \figuremacrob{fragmentation}{Snapshots for molecular dynamics trajectory with impact parameter of 0.5~\AA{} and a collision energy of 27.5 eV leading to intramolecular fragmentation.} \textbf{Dissociation cross section} The opacity curves are presented in Figure \ref{opacitycurves} for various collision energies. At low impact parameters, the dissociation is very efficient even at low collision energy. At the lowest collision energy of 2.5~eV, the opacity curve presents a smooth decrease from 2 to 5~\AA, whereas for collision energies larger than 10~eV, all curves are very similar. These high energy curves show high dissociation probability below 3.5~\AA, reach 50\% at about 4.5~\AA{} and drop to zero for values larger than 5.5~\AA. These values can be compared to the van der Waals radius of argon (1.88~\AA) plus half of (i) the distance between the two molecules centers of masses (3.04~\AA), (ii) the smallest (6.82~\AA) or (iii) largest pyrene axes (9.25~\AA) leading to distances of 3.40, 5.29 and 6.50~\AA, respectively. Below 3.40~\AA, all trajectories involve a frontal impact of the argon on the dimer carbonaceous system and almost all of them lead to dissociation. Unexpectedly, the opacity curve drops to zero at 5.5~\AA{} which is lower than the largest computed value of 6.5~\AA. Interestingly, taking the largest distance between carbon atoms in pyrene (7.0~\AA) instead of that between hydrogen atoms (9.25~\AA) leads to a value of 5.4~\AA{} which is in line with the opacity curves. This suggests that the dissociation is efficient only if the carbonaceous skeleton area is impacted, the impact in the region of external hydrogen atoms resulting mostly in an intramolecular C-H mode excitation at the expense of dissociative modes. As a conclusion, it seems that for energies larger than 10~eV, the opacity curves are similar as they are driven by simple geometric rules, in other words, if the dimer receives a direct impact of the argon on the carbonaceous skeleton area, it will dissociate. Interestingly, this seems to be in agreement with previous works \cite{Chen2014, Gatchell2016knockout} which also pointed out the efficient nuclear stopping power of carbon atoms in a very different context (higher energy collisions leading to knock-out process). \figuremacro{opacitycurves}{Opacity curves as a function of the impact parameter $b$ for several selected center of mass collision energies.} The blue curve in Figure \ref{cross-section} shows the MD dissociation cross sections of pyrene dimers obtained from the opacity curves following eq. \ref{integ}. It presents a steep increase for energies bellow 7.5~eV before remaining almost constant around 65~\AA$^2$ for collision energies greater than 10-15~eV. This is a direct consequence from the already discussed similarity of opacity curves for the high collision energies. The purple curve corresponds to dissociation at infinite timescales. Figure \ref{cross-section} also reports the cross sections computed from the MD+PST model. It can be seen that, for low collision energies, the MD and MD+PST cross sections are very close, indicating that most of the dissociations occur on the short timescales. On the opposite, at high collision energies, a non-negligible fraction of the dimers, which are not dissociated at the end of the MD simulation, carry enough energy to evaporate on the experimental timescales. At the experimental center of mass collision energy of 17.5~eV, the MD+PST cross section (about 70~\AA$^2$) is slightly above the pure MD dissociation ratio, which indicates that the dissociation at long timescales represents a small fraction of the dissociated pyrene dimers as already seen from the TOF spectra analysis (see Figure \ref{expTOF}). We have also plotted in Figure~\ref{cross-section} the model cross section $\sigma_\infty$ that successfully reproduced the threshold collision-induced dissociation experimental results \cite{Zamith2020threshold}. This model cross section is obtained by considering that the collision energy transfer is given by the LOC model and the expression for the cross section is given by: \begin{equation} \sigma_{LOC}(E_{col}) = \sigma_0 (E_{col} - D)/(E_{col}). \end{equation} where $D=1.08$~eV is the dissociation energy \cite{Dontot2019, Zamith2020threshold} and $\sigma_0=63$~\AA$^2$ is a scaling factor usually thought as the geometrical cross section. This model cross section is usually further convolved with dissociation rates, collision energy distributions and internal energy distributions in order to be compared with experimental results. However, since here for the theoretical calculations there is no collision energy distribution, this curve could in principle be directly compared with the purple one in Figure~\ref{cross-section}, namely the cross section for infinite time. One can see that the MD, MD+PST results and the model cross section have similar collision energy dependence. The magnitude of the two cross sections is rather different at high collision energy, with about 60~\AA$^2$ and 74~\AA$^2$ for the model and infinite timescale cross sections respectively. Nevertheless, this difference is probably within the error bars of the experimental cross section measurement. The dissociation cross sections for MD timescales with time step being 0.1 fs at collision energy of 20 and 25 ~eV ($\sigma_\mathrm{MD}$(0.1)) in Figure~\ref{cross-section} are close to the ones of time step being 0.5 fs ($\sigma_\mathrm{MD}$), which indicates a time step of 0.1 fs used in the MD simulation does not change significantly the corresponding dissociation cross section. \figuremacro{cross-section}{Dissociation cross sections of Py$_2^+$ after collision with argon as a function of center of mass collision energy for the short (MD), experimental (MD+PST) and infinite timescales. Cross sections resulting from the LOC model are also plotted. $\sigma_\mathrm{MD}$ (0.1) denotes the dissociation cross section for short (MD) timescale with a time step of 0.1 fs.} \textbf{Energy partition} The mean value obtained for the transferred energy after removing the translation kinetic energy of the dimer, namely $\Delta E_{int}^{Py_2}$, is plotted in Figure \ref{transferredE-Ar-300} as a function of the collision energy. Although this quantity evolves almost linearly with the collision energy, the curves are different when one considers only the trajectories leading to dissociation or non-dissociation. For trajectories where the dimer does not dissociate, $\Delta E_{int-ud}^{Py_2}$ remains small for all collision energies below 20~eV and shows a very slight increase for collision energies larger than 20~eV. For trajectories leading to dissociation, $E_{int-d}^{Py_2}$ grows almost linearly, but above 10-15~eV most of the absorbed energy is actually used to heat the individual monomers (the green curve) whereas the energy given in the dissociative mode (difference between the blue and green curves) remains almost constant. We note that, despite the trends of the mean energy values derived from all simulations or restricted to the undissociated cases are interesting, their absolute values have small meaning as they depend on the arbitrarily chosen $b_{max}$ value, {\it i.e.} increasing $b_{max}$ would result in more undissociated trajectories with less and less energy transfer. On the opposite, absolute values of mean energies for the dissociation trajectories are relevant, as increasing the $b_{max}$ value would not result in new dissociation trajectories. For MD simulations with time step being 0.1 fs at collision energy of 20 and 25 ~eV, the corresponding energies in Figure~\ref{transferredE-Ar-300} are close to the ones of time step being 0.5 fs, which indicates a time step of 0.1 fs used in the MD simulation does not change significantly the corresponding deposition of the total transferred energy. \figuremacro{transferredE-Ar-300}{At the end of the MD collision simulations with a time step of 0.1 and 0.5 fs, the total transferred energy $\Delta E_{int}^{Py_2}$ to the rovibrational modes or restricted to the sole dissociated ($\Delta E_{int-d}^{Py_2}$) or undissociated ($\Delta E_{int-ud}^{Py_2}$) pyrene dimers as a function of collision energy. The transeferred energy to the monomers rovibrational modes for the dissociated dimers $\Delta E_{int-d}^{Py^1+Py^2}$ is also plotted.} It is also interesting to focus on the kinetic energy partition, in particular because its decomposition in sub-contributions (dissociative vs non-dissociative modes) is easier (see eqs. \ref{Eparti}) than that of the potential (and consequently total) energy. For each simulated collision energy, the values for the kinetic energy sub-contributions (eqs. \ref{Eparti}) are averaged over all the trajectories and reported in Table \ref{tab:table1} and Figure \ref{Epartition-Ar-300-SP}. In addition, the averaged kinetic energy sub-contributions are also calculated over the dissociated and undissociated trajectories separately for each simulated collision energy. Then sum the contributions of dissociated and undissociated trajectories (black curves in Figure \ref{Epartition-Ar-300-SP}) calculated by the following eqs \ref{separately}, \begin{align} \label{separately} E_{1}^k &=E_{td-d}^k*P+E_{td-ud}^k*(1-P) \nonumber \\ E_{2}^k &=E_{Re-d}^k*P+E_{Re-ud}^k*(1-P) \\ E_{3}^k &=(E_{Py^1}^k + E_{Py^2}^k)_{-d}*P+(E_{Py^1}^k + E_{Py^2}^k)_{-ud}*(1-P) \nonumber \end{align} where $P$ refers to the dissociation probability at a given center of mass collision energy. The results of sum the contributions of dissociated and undissociated trajectories are the same with the ones over all the trajectories (see Figure \ref{Epartition-Ar-300-SP}). This ensures our calculations for the mean kinetic energy sub-contributions are right. \begingroup \setlength{\tabcolsep}{6pt} \renewcommand{\arraystretch}{1.1} \begin{table} \begin{center} \caption{The kinetic energy partition after the collision of pyrene dimer with argon at different collision energies $E_{col}$. All energies are in eV.} \label{tab:table1} \begin{tabular}{|c|c|c|c|c|c|} \hline \boldm{$E_{col}$} & \boldm{$E^k_{td}$} & \boldm{$E^k_{Ar}$} & \boldm{$E^k_{Py^1}$} & \boldm{$E^k_{Py^2}$} & \boldm{$E^k_{Re}$} \\ [0.08cm] \hline 2.5 & 0.17 & 1.67 & 0.25 & 0.25 & 0.19\\ \hline 5.0 & 0.30 & 3.52 & 0.40 & 0.41 & 0.41\\ \hline 7.5 & 0.42 & 5.38 & 0.55 & 0.56 & 0.64\\ \hline 10.0 & 0.53 & 7.32 & 0.71 & 0.70 & 0.80\\ \hline 12.5 & 0.65 & 9.20 & 0.89 & 0.88 & 0.95\\ \hline 15.0 & 0.74 & 11.20 & 1.03 & 1.03 & 1.06\\ \hline % He & 17.5 & 0.15 & 15.40 & 0.56 & 0.55 & 0.26\\ % Ne & 17.5 & 0.52 & 11.95 & 1.07 & 1.08 & 0.84\\ 17.5 & 0.82 & 13.16 & 1.16 & 1.22 & 1.18\\ \hline 20.0 & 0.89 & 15.12 & 1.37 & 1.32 & 1.30\\ \hline 22.5 & 0.96 & 17.09 & 1.52 & 1.51 & 1.36\\ \hline 25.0 & 1.01 & 19.18 & 1.61 & 1.69 & 1.45\\ \hline % 27.5 & 1.09 & 21.16 & 1.78 & 1.83 & 1.49\\ % 30.0 & 1.16 & 23.01 & 2.09 & 1.98 & 1.54\\ \end{tabular} \end{center} \end{table} \endgroup \figuremacro{Epartition-Ar-300-SP}{Mean kinetic energy partition at the end of the MD simulations.} The comparison of these kinetic energy sub-contributions between the time step being 0.1 and 0.5 fs used in the MD simulations is shown in Table \ref{tab:table2} and Figure \ref{Epartition-Ar-300-Tstep-01}, which indicate a time step of 0.1 fs almost didn't affect the results of kinetic energy sub-contributions. \begingroup \setlength{\tabcolsep}{6pt} % Default value: 6pt \renewcommand{\arraystretch}{1.1} % Default value: 1 \begin{table} \begin{center} \caption{The kinetic energy partition and cross section at the end of MD simulations with time step being 0.1 and 0.5 at different collision energies of 20 and 25. All energies are in eV. Time step ($Tstep$) is in fs. Cross section $\sigma_{_{MD}}$ is in ~\AA.} \label{tab:table2} \begin{tabular}{|c|c|c|c|c|c|c|c|} \hline \boldm{$E_{col}$} & \boldm{$Tstep$} & \boldm{$E^k_{td}$} & \boldm{$E^k_{Ar}$} & \boldm{$E^k_{Py^1}$} & \boldm{$E^k_{Py^2}$} & \boldm{$E^k_{Re}$} & \boldm{$\sigma_{_{MD}}$} \\ [0.08cm] \hline 20.0 & 0.1 & 0.89 & 15.18 & 1.34 & 1.35 & 1.28 & 64.11\\ \hline 20.0 & 0.5 & 0.89 & 15.12 & 1.37 & 1.32 & 1.30 & 64.45\\ \hline 25.0 & 0.1 & 1.03 & 19.04 & 1.71 & 1.67 & 1.44 & 62.86\\ \hline 25.0 & 0.5 & 1.01 & 19.18 & 1.61 & 1.69 & 1.45 & 64.77\\ \hline % 30.0 & 1.16 & 23.01 & 2.09 & 1.98 & 1.54\\ \end{tabular} \end{center} \end{table} \endgroup \figuremacro{Epartition-Ar-300-Tstep-01}{Mean kinetic energy partition at the end of the MD simulations with time step being 0.5 fs at the center of mass collision energy from 2.5 to 25 eV. The mean kinetic energy partition with time step being 0.1 fs at center of mass collision energies of 20 and 25 eV are plotted with filled round circles.} In Figure~\ref{prot-Ar-300} are reported the ratios of the pyrene dimer translational kinetic energy $E^k_{td}$, relative kinetic energy $E^k_{Re}$ and monomers rovibrational kinetic energies $E^k_{Py^1}+E^k_{Py^2}$ over the total pyrene dimer kinetic energy $E^k_{tot}-E^k_{Ar}$. It clearly appears that, whereas the contribution of the dimer translation kinetic energy ($E^k_{td}$) remains almost constant (very slight decrease from about 20\% to 18\% of the dimer kinetic energy), this is not the case for the other two contributions. For collision energies below 7.5~eV, the proportion of the kinetic energy associated to the center of mass relative velocities increases whereas the opposite is observed for the monomers rovibrational kinetic energy. These two trends are reversed above 7.5~eV. \figuremacro{prot-Ar-300}{Kinetic energy proportion after collision of Py$_2^+$ with argon as a function of collision energy.} Again, it is convenient to analyse separately the kinetic energy partition for trajectories leading to dissociation or not as done in Figure \ref{Epartition-Ar-300-d-ud} for $E^k_{td}$, $E^k_{Re}$ and $E^k_{Py^1}+E^k_{Py^2}$. For both dissociated and undissociated trajectories, the total energy in the system computed from the initial energy at 25 K (0.32 eV) plus the transferred energy is twice the final kinetic energy computed from velocities shown in Figure \ref{Epartition-Ar-300-d-ud} (black curves). This is exactly what one should expect from the Virial theorem. %which proves that the calculations of kinetic energy are as good as the dynamics. In the absence of dissociation, the transferred energy is either small or redistributed over all the vibrational modes of the dimer, leading to small values for $E^k_{Re-ud}$ (mean value always below 0.04~eV). The monomers rovibrational kinetic energies remain constant with an increase for collision energies above 20~eV, indicating that the slight increase of transferred energy results in a heating of the monomers, as already inferred from Figure \ref{transferredE-Ar-300}. Once a dimer dissociates, the two pyrene molecules relative kinetic energy $E^k_{Re-d}$ can not be transferred back to the intramolecular modes and its mean value is never negligible with respect to the monomers rovibrational kinetic energies $(E^k_{Py^1}+E^k_{Py^2})_{-d}$. However, although the slope of $(E^k_{Py^1}+E^k_{Py^2})_{-d}$ remains constant with collision energies, that of $E^k_{Re-d}$ decreases clearly. This is in line with the analysis of Figure \ref{transferredE-Ar-300}, which shows the amount of energy transferred to the dissociative modes remains constant for high collision energies whereas the monomers are getting more internal energy. \figuremacro{Epartition-Ar-300-d-ud}{Kinetic energy partition for dissociated (-d) and undissociated (-ud) trajectories at the end of the MD simulation as a function of collision energy.} Finally, some characteristic timescales are computed, which are presented in Figure \ref{figuretimescale}. They correspond to the timescales for the argon with its initial velocity to travel across some characteristic distances, namely, a C-H (1.10~\AA) or a C-C bond (1.40~\AA) and the largest molecular axis (9.25~\AA). These timescales can be compared with those of the pyrene dimer vibrational modes as an efficient energy transfer would be favored by similar orders of magnitudes. The intermolecular dimer modes possibly mixed with very soft folding modes are lying within the 70-120 $cm^{-1}$ spectral range \cite{Dontot2020} with corresponding half-periods of 130-240 fs. These timescales are of the same order of magnitude as the time for the argon to travel across the largest pyrene axis for collision energies below 10~eV. Typical frequencies for intramolecular non-soft modes are lying from 500 $cm^{-1}$ to 3000 $cm^{-1}$ (C-H stretching modes), leading to half-periods of 5-33 fs. For all the simulated collision energies, the characteristic times required for the argon to travel across typical C-H or C-C bond distances belong to the same order of magnitude as some of the intramolecular hard modes. Therefore, it appears from this qualitative description that the collision energy transfer toward the intermolecular modes is easier at collision energies lower than 10 eV whereas the transfer toward intramolecular modes is efficient for all the simulated collision energies. This is actually in line with the fact that the part of the absorbed collision energy taken by the non-soft intramocular modes is increasing with the collision energy at the expense of that taken by the intermolecular and soft intramolecular modes, which is in agreement with the previous energy analysis (Figures~\ref{transferredE-Ar-300}, \ref{prot-Ar-300} and \ref{Epartition-Ar-300-d-ud}). \figuremacro{figuretimescale}{Timescales, as a function of center of mass collision energy, for argon to travel across some typical distances: a carbon-carbon bond (green), a carbon-hydrogen bond (purple) or the largest axis of the pyrene molecule (blue).} \textbf{Efficiency of energy transfer within the dimer} \figuremacro{T-time-zoom_abcdef}{Instantaneous kinetic temperatures as a function of time for intra and intermolecular modes of the pyrene dimer at a collision energy of 22.5~eV. Impact parameters $b$ are (a) 2, (b) 3, (c) 0, (d) 2.5, (e) 2, and (f) 2~\AA{}. In cases (a) and (b) dissociation takes place whereas in the other cases the dimer remains undissociated at the end of the MD simulation. In (c) to (f) the lower panel is a vertical zoom of the corresponding intramolecular parts in upper panel.} \figuremacrob{E-time-abcdef}{Instantaneous kinetic energies as a function of time for intra and intermolecular modes of the pyrene dimer at a collision energy of 22.5 eV. Impact parameters $b$ are (a) 2, (b) 3, (c) 0, (d) 2.5, (e) 2, and (f) 2 \AA{}. In cases (a) and (b), dissociation takes place whereas in the other cases the dimer remains undissociated at the end of the MD simulation.} In this section, we address how the energy is shared inside the dimer after the collision. In particular, we look at the efficiency of energy transfer between the intramolecular modes of each unit and the intermolecular modes. The amount of deposited energy as well as its partition between the intramolecular modes of each molecule and the intermolecular modes is strongly dependent on the collision condition: the impact parameter, the orientation of the dimer, whether a head on collision occurs with one of the dimer atoms (and its nature, carbon or hydrogen). This results in very different evolutions of the subsequent energy flows for which precise values concerning timescales can hardly be derived. Nevertheless, the analysis of the trajectories allows to identify some characteristic behaviors. In order to estimate the thermalization process efficiency, the instantaneous intra and intermolecular kinetic temperatures are evaluated using the following formula: \begin{align} \label{kineticT} T^k=2 \frac{}{nk_b} \end{align} %$$ T^k=2 \frac{}{nk_b} $$ where $k_b$ refers to the Boltzmann constant. $n$ is the number of involved modes and $E^k$ is the kinetic energy for the intra or intermolecular modes (see eqs.~\ref{Eparti}). $T^k$ is plotted in Figure \ref{T-time-zoom_abcdef} for some selected trajectories obtained for collision energies of 22.5 eV and various impact parameters. The evolution of the corresponding energies ($E^k_{intra^1}$, $E^k_{intra^2}$ and $E^k_{inter}$) are presented in Figure \ref{E-time-abcdef}. In Figures \ref{T-time-zoom_abcdef} and \ref{E-time-abcdef}, simulations (a) and (b) correspond to trajectories for which dissociation occurred, whereas the dimer remained intact in the other simulations. In the simulation (a), a larger amount of energy is deposited in the first monomer with respect to the second one. The dissociation occurs before an efficient energy transfer takes place between the two monomers, leading to one hot monomer and one cold monomer at the end of the simulation. The situation is slightly different in the dissociation trajectory (b): there is a much smaller difference between the energies received during the collission by each of the monomers. One can observe that the equilibration of the two monomers intramolecular energies can take place before dissociation, leaving the two monomers with similar energies/kinetic temperatures. In the other four pictures (c, d, e, and f), corresponding to undissociated trajectories, one can see that the thermalization between the two monomers intramolecular modes occurs with timescales from 0.2 to 1.5 ps shown in Figure \ref{T-time-zoom_abcdef}. On the other side, the energy equilibration between intra and intermolecular modes takes more time. Indeed, the thermalization is almost achieved in simulations (c) and (d) at 1.5 and 2.5 ps, respectively, but would take more than the simulated duration 3 ps for trajectories (e) and (f) displayed in Figure \ref{T-time-zoom_abcdef}. As a conclusion of these trajectories analyses, it seems that the thermalization between intramolecular modes of the two monomers is relatively efficient (on the order of ps). On the other hand, the thermalization with the intermolecular modes is less efficient and sometimes is not observed during the simulated time of 3 ps. The direct dissociation of the dimer is a fast process (on the order of a few tenths of ps) which may prevent the thermalization taking place, leading to monomer temperatures reflecting the initial energy deposition. \subsection{Conclusions about CID of Py$_2^+$} A QM/MM dynamics study of the collision of Py$_2^+$ with argon at various collision energies were carried out. Argon was treated as a polarisable MM particle and Py$_2^+$ was treated using the SCC-DFTB method. In the dynamical simulations, a time step of 0.5 fs is proper even for high collision energies 25 eV. The TOF mass spectra of parent Py$_2^+$ and dissociation product Py$^+$ were simulated by the PST using the MD outputs at a centre of mass collision energy of 17.5~eV. With respect to TOF mass spectra extracted from pure PST simulations, considering non-statistical dissociation processes that take place before the energy redistribution from MD simulations improves the match between experimental and theoretical TOF spectra. The agreement between the measured and simulated mass spectra peak shapes and positions shows that the essence of the collision-induced dissociation is captured by the simulation. It appears that the TOF spectra mostly result from dimers dissociating on short timescales (during the MD simulation) and the remaining minor contribution is from dimers dissociating at longer timescales (the second step, during PST calculation). This indicates that Py$_2^+$ primarily engages a direct dissociation path after collision. The extraction of snapshots from the MD simulations allows to visualize the collision processes. It shows that the evolution of the trajectories either toward a dissociation or a redistribution of the transferred energy strongly depends on the initial collision conditions. Intramolecular fragmentation of the monomers occurred only for collision energies above 25 eV. The dissociation cross sections show a steep increase for collision energies below 7.5~eV and remain almost constant for collision energies greater than 10~eV. The dissociation cross section of Py$_2^+$ increases when dissociation occurring on longer timescale is included. As such, the dissociation cross section computed from the MD+PST model at the centre of mass collision energy of 17.5~eV is slightly higher than the value derived from pure MD simulations. The analysis of the partition of the final kinetic energy as a function of the collision energy shows how the absorbed energy is shared between the dissociative modes and the heating of individual monomers. It shows that above 7.5~eV, increasing the collision energy mostly results in an increase of the intramolecular energy. The qualitative analysis of the different timescales involved in the collision further supports the kinetic energy partition analysis. Finally, the analysis of energy transfer efficiency within the dimer suggests that direct dissociation is too fast to allow significant thermalization of the system. On the other hand, when there is no dissociation, thermalization can occur with a faster equilibration between the intramolecular modes of the two units than with the intermolecular modes. The present results can be compared with experimental and theoretical works discussing the direct and indirect fragmentation of PAH and PAH clusters submitted to higher energy collisions \cite{Chen2014, Gatchell2016knockout}. These authors showed that the nuclear stopping power dominates over the electronic one below 1 keV, giving a justification to our approach based on classical MD and PST. They also showed that the direct non-statistical PAH fragmentation (knock-out) is an efficient process above 20 eV. This is in line with the fact that monomer fragmentation was only observed in our MD simulations above 25 eV. Our work shows that, for PAH clusters, a regime exists below this collision energy where the dimer dissociation is governed by non-statistical processes. In this study, the collision process, dissociation path, energy partition and distribution, and the efficiency of energy transfer were deeply explored for the Py$_2^+$ system, which can provide valuable reference for the CID study of larger PAH cation clusters. %Beyond the specific study of Py$_2^+$ dissociation, the methodology paves the way for future analysis of CID experiments like dissociation of molecular clusters such as water-Uracil clusters. %: ----------------------- paths to graphics ------------------------ % change according to folder and file names \ifpdf \graphicspath{{X/figures/PNG/}{X/figures/PDF/}{X/figures/}} \else \graphicspath{{X/figures/EPS/}{X/figures/}} \fi