E Leroux

Frédéric Briend
email: briend@cyceron.fr.

GeodesicSlicer: A Slicer toolbox for targeting brain stimulation

Keywords: NIBS, rTMS, 3D Slicer, EEG, target determination, correction factor

NonInvasive Brain Stimulation (NIBS) is a potential therapeutic tool with growing interest, but neuronavigation-guided software and tools available for the target determination are mostly either expensive or closed proprietary applications. To address these limitations, we propose GeodesicSlicer, a customizable, free, and open-source NIBS therapy research toolkit. GeodesicSlicer is implemented as an extension for the widely used 3D Slicer medical image visualization and analysis application platform. GeodesicSlicer uses cortical stimulation target from either functional or anatomical images to provide functionality specifically designed for NIBS therapy research. The provided algorithms are tested and they are accessible through a convenient graphical user interface. Modules have been created for NIBS target determination according to the position of the electrodes in the 10-20 system electroencephalogram and calculating correction factors to adjust the repetitive Transcranial Magnetic Stimulation (rTMS) dose for the treatment. Two illustrative examples are processing with the module. This new open-source software has been developed for NIBS therapy: GeodesicSlicer is an alternative

for laboratories that do not have access to neuronavigation system. The triangulation-based MRI-guided method presented here provides a reproducible and inexpensive way to position the TMS coil that may be used without the use of a neuronavigation system.

Introduction

A key issue in the field of NonInvasive Brain Stimulation (NIBS) is to determine an accurate localization on the scalp to correctly target cortical areas knowing the great anatomical variability of the brain. Since personalized medicine for the treatment of psychosis allows for the consideration of substantial inter-individual variability, recent findings claim that brain stimulation can be guided in a personalized manner (Briend et al. Under Review; Kraus and Gharabaghi 2015; Lahti 2016; I. E. [START_REF] Sommer | Toward personalized treatment of hallucinations[END_REF].

Most clinical applications of the NIBS are based on probabilistic targeting methods which do not account for individual anatomical variability (e.g. for major depressive episodes the so called "5-cm rule" (U. [START_REF] Herwig | Transcranial magnetic stimulation in therapy studies: examination of the reliability of "standard" coil positioning by neuronavigation[END_REF] or the International 10-20 electroencephalogram (EEG) [START_REF] De Witte | Left prefrontal neuronavigated electrode localization in tDCS: 10-20 EEG system versus MRI-guided neuronavigation[END_REF][START_REF] Herwig | Using the international 10-20 EEG system for positioning of transcranial magnetic stimulation[END_REF] or derivative system [START_REF] Beam | An efficient and accurate new method for locating the F3 position for prefrontal TMS applications[END_REF]). This may lead to suboptimal clinical responses when compared to individualized targeting techniques based on structural brain scanning. Research and clinical studies require accuracy and precision not offered by these probabilistic targeting methods [START_REF] Herbsman | Anatomically based targeting of prefrontal cortex for rTMS[END_REF]. For example, a common and easy method for the positioning of the coil in psychiatric therapies uses the standardized T3P3 site according to the International 10-20 system of EEG electrode positioning [START_REF] Jasper | The ten twenty electrode system of the international federation[END_REF]). However, this method is known to be an inaccurate estimation, especially given its variable projections on the individual brain [START_REF] Briend | A new toolbox to compare target localizations for non-invasive brain stimulation: An application of rTMS treatment for auditory hallucinations in schizophrenia[END_REF][START_REF] Herwig | Using the international 10-20 EEG system for positioning of transcranial magnetic stimulation[END_REF]. It is why, there is a need for personalized target method that uses the participant's own anatomical or functional images to guide target placement.

As a personalized target method, the combination of brain imaging and a neuronavigation system in the field of NIBS may improve the efficacy of stimulation treatment (U. [START_REF] Herwig | Transcranial magnetic stimulation in therapy studies: examination of the reliability of "standard" coil positioning by neuronavigation[END_REF][START_REF] Sommer | Toward personalized treatment of hallucinations[END_REF]; I. E. C. [START_REF] Sommer | Can fMRI-guidance improve the efficacy of rTMS treatment for auditory verbal hallucinations[END_REF], however, there are some disadvantages. These include the high cost of these systems, which can exceed $50,000, the complexity to use, the space consumed by the device and the difficulty in using these systems for the study of posterior brain areas located in the blind spot of the neuronavigation system [START_REF] Vaghefi | MRI Guided Brain Stimulation without the Use of a Neuronavigation System[END_REF].

In order to propose an alternative that can combine the accuracy and the simplicity of the two previous methods, we developed an open-source tool "GeodesicSlicer", which facilitates the stimulation site determination, allowing users to manually posit the target of repetitive Transcranial Magnetic Stimulation (rTMS) coil over a cortical target derived from functional or anatomical images. This module creates a 3D mesh morphed to the structural MRI head data of the participant then projected an individualized 10-20 system EEG and the cortical stimulation target on it. Then, the module calculates the geodesic distances between the projected stimulation target and the position of the 3 nearest electrodes in the individualized 10-20 system EEG in order to guide the stimulation. Our technique takes triangulation-based MRI-guided method as Andoh et al. [START_REF] Andoh | A triangulation-based magnetic resonance image-guided method for transcranial magnetic stimulation coil positioning[END_REF]) that devised a method of targeting NIBS using an anatomical scan only. Moreover, it was proposed that rTMS inter-individual variability in its efficacy for treating patient could be attributed to variations in the cortical anatomy in Schizophrenia (Ralph E. [START_REF] Hoffman | Transcranial magnetic stimulation of Wernicke's and right homologous sites to curtail "voices:" a randomized trial[END_REF] or in major depression disorders (MDD) (Trojak et al. 2012). It is why, we implemented in this module, correction factors, according to scalp-to-cortex distance [START_REF] Summers | BrainRuler-a free, open-access tool for calculating scalp to cortex distance[END_REF], to adjust the rTMS dose for the treatment.

We propose GeodesicSlicer as a common easy-to-use software for NIBS site determination, thanks to its implementation in 3D Slicer, a powerful tools for neuroimaging [START_REF] Pieper | The NA-MIC Kit: ITK, VTK, pipelines, grids and 3D slicer as an open platform for the medical image computing community[END_REF].

Methods and materials

Implementation of GeodesicSlicer

Platform

Geodesic Slicer was implemented in 3D Slicer [START_REF] Pieper | The NA-MIC Kit: ITK, VTK, pipelines, grids and 3D slicer as an open platform for the medical image computing community[END_REF], a software which is freely downloadable from the website http://www.slicer.org. 3D Slicer provides an immense amount of functionality to visualize and analyze a wide range of datasets, such as anatomical/functional images, image segmentation results and surface models. Also, it supports import and export data from a wide range of standard data formats. In addition, the 3D Slicer has a widespread use in project research and is more and more downloaded [START_REF] Pinter | SlicerRT: radiation therapy research toolkit for 3D Slicer[END_REF].

Implementation

The GeodesicSlicer module is written in Python. Python is a very popular and easily interpreted language, which allows multiple programming paradigms, including objectoriented, imperative and functional programming styles. Inheriting from 3D Slicer, GeodesicSlicer is available for Windows, Linux, and Mac OS X platforms.

Our implementation of the algorithm in 3D Slicer consists of a graphical user interface front-end to enable interactions of the user with the image and several algorithms back-end. It allows the generations of head surface mesh reconstruction and individualized 10-20 system EEG. Moreover, it can compute the geodesic distances between the target and electrodes landmarks and compute 2 correction factors to adjust rTMS dose for the treatment.

The geodesic distances (i.e. the shortest path between two points in a curved space) to draw the individualized 10-20 system EEG or compute the distances between the target and electrodes landmarks are calculated on a 3D mesh morphed to the structural MRI head data of the participant thanks to the implemented Dijkstra's algorithm [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF], which calculates the shortest path between the vertex of triangle mesh.

Two correction factors to adjust rTMS dose for the treatment for individual subjects are given by the software. First, Stokes et al. [START_REF] Stokes | Distance-adjusted motor threshold for transcranial magnetic stimulation[END_REF], proposed increasing the stimulation intensity by about 3% for each additional millimeter between the coil and the scalp surface. Second, [START_REF] Hoffman | Transcranial magnetic stimulation of Wernicke's and right homologous sites to curtail "voices:" a randomized trial[END_REF]) also take into account of the skin-surface-to-cortical-surface, but their adjustment also reflected the fact that magnetic field strength falls off exponentially relative to distance to the center of the coil (supplementary material of their article). These corrections factors are, for Stokes et al. [START_REF] Stokes | Distance-adjusted motor threshold for transcranial magnetic stimulation[END_REF],

where [AdjMT% = 2,7*(SCDx -SCDm) + rMT] and according to Hoffman et colleagues (Ralph E. [START_REF] Hoffman | Transcranial magnetic stimulation of Wernicke's and right homologous sites to curtail "voices:" a randomized trial[END_REF], where [AdjMT% = 0.90*rMT*e0.036*(SCDx-SCDm)], where AdjMT is the adjusted motor threshold in percent (%), rMT is the unadjusted resting motor threshold in % of stimulator output, SCDx is the scalp-to-cortex distance between the scalp and the cortical stimulation site and SCDm is the scalp-to-cortex distance between the scalp and the primary motor cortex (M1).

Licensing and distribution

GeodesicSlicer (WikiPage) is distributed under a CeCill license. The software may be used not only for research purposes but also in clinical and commercial projects. Note, however, that validation for a particular clinical purpose is an onus of the user. The brain stimulation guidelines and safety procedures are dependent on each neurostimulation therapy used [for example in TMS: [START_REF] Rossi | Safety, ethical considerations, and application guidelines for the use of transcranial magnetic stimulation in clinical practice and research[END_REF]], but not directly to the use of this software.

GeodesicSlicer modules can be downloaded as an extension for 3D Slicer 4.10.0 or higher. All the presented software is open-source and the source code is available on GitHub, which contains detailed guides for user's installation and usage. The authors declare that they have no conflict of interest.

GeodesicSlicer, a MRI-guided method

First, GeodesicSlicer provides realistic and accurate 3D representations of the head scalp. The, the MRI-guided method uses participant's individual MRI to determine the TMS coil position onto the head surface. We describe in the following section the workflow of GeodesicSlicer to position the TMS coil over the participant's head surface from her/his T1weighted anatomical/functional images.

Procedure of GeodesicSlicer

1) Loading of the T1-weighted image into 3D Slicer.

2) 3D representation of the head (the head surface mesh or more accurately, triangle meshes) was individually reconstructed in native space from the T1-weighted wholebrain anatomical image using 3D Slicer software ("editor toolbox" version 4.8).

3) Manual identification of four anatomical landmarks for the essential positioning of the electrodes on the head surface mesh: the nasion, the inion, the left and right tragi (in this order). The Dijkstra's algorithm automatically reconstructed the 10-20 system EEG with T3P3 in the middle of the segment delimited by T3 and P3. For that, the shortest paths between the nasion and inion and the left and right tragi that passed through the center point of the head (electrode Cz). Then, always with the shortest path algorithm, all electrodes are located according to their standardized that represent proportions of the measured distance from the nasion to the inion and from the left to the right tragi [START_REF] Klem | The ten-twenty electrode system of the International Federation[END_REF]. 4) Manual placement of the cortical stimulation target on the T1-weighted image, the projection of it onto the head surface mesh was made by using a classical 3D Euclidean distance, i.e. √ (x2-x1)² + (y2-y1)² + (z2-z1)².

5) The fifth step was the computation of the geodesic distances between this projected target on the surface and the electrodes of the 10-20 system EEG corresponding of the participant's head. These three distances (in cm) were then used to triangulate and to position the TMS coil manually over the participant's head.

6) The last step is needed to adjust the rTMS dose for treatment. After placing cortical landmark in M1 according to the Yousry's method [START_REF] Yousry | Localization of the motor hand area to a knob on the precentral gyrus. A new landmark[END_REF], the brain area to determine the motor threshold in rTMS, and giving rMT of stimulator output, the software gives two AdjMT of stimulator output.

The duration of this workflow lasted about 10 minutes per participant.

Illustrative Examples Using GeodesicSlicer

Use cases are presented to demonstrate the capabilities of GeodesicSlicer extension for addressing clinically relevant rTMS site determination. We chose brain imaging data of two patients from previous study of our team [START_REF] Dollfus | High-Frequency Neuronavigated rTMS in Auditory Verbal Hallucinations: A Pilot Double-Blind Controlled Study in Patients With Schizophrenia[END_REF]), but the rTMS stimulation proposed here is just theoretical. The subjects have previously written informed consent and these studies were approved by a local ethical committee. All coordinates are given in the MRI native space.

Accuracy of the measure

Eight controls (35.77 ± 5.29 years; 2 women) were included to assess the validity of the GeodesicSlicer method. The placement of the nasion, inion, and the two tragi determines the position of the electrodes in the 10-20 system EEG. We measured the distances from the nasion to the inion and from the left tragus and the right tragus in the controls with a measuring tape and compared them to the same distances calculated by GeodesicSlicer. Bland-Altman plots assess retest reliability of two measures and were used to test the stability across these distances [START_REF] Bland | Measuring agreement in method comparison studies[END_REF].

Results

Illustrative Examples Using GeodesicSlicer Case 1: Determination of the projected stimulation target in one patient with schizophrenia with auditory verbal hallucinations

The rTMS can be used as treatment for auditory verbal hallucinations (AVHs), notably, in the case of refractory to treatments (R. E. [START_REF] Hoffman | Transcranial magnetic stimulation of left temporoparietal cortex in three patients reporting hallucinated "voices[END_REF]). In the case described below, we will consider that patient receive a treatment by rTMS applied over a precise anatomical site in the left temporal region, that significant effects in AVH reduction.

The patient is a 35-year-old man, diagnosed with schizophrenia (based on the DSM-V, Diagnostic and Statistical Manual of Mental Disorder 5 th edition) and who suffers of constant AVHs. He was recruited from the University Hospital (Caen, France).

rTMS site determination

The patient underwent a structural MRI on a 3T scanner (Intera Achieva 3T, Philips Medical System, the Netherlands) with a three-dimensional (3D), high-resolution T1-weighted structural volume (T1 TFE sequence, 256 x 256 matrix size with 180 contiguous slices, field of view (FOV) = 256 mm, 1 mm isotropic resolution, antero-posterior slice orientation, repetition time = 6.914 ms, echo time = 3.16 ms, flip angle = 6, inversion time = 940 ms).

Using this system, the cortical stimulation target was localized at the crossing between the projection of the ascending branch of the left lateral sulcus and the left superior temporal sulcus. To do that, using a sagittal section of the structural MRI to visualize the upper sylvian and temporal sulcus, we propose to consider the intersection between the orthogonal projection of the verticalization of the Sylvius fissure and the upper left temporal sulcus (see Supplementary Data Video 1 in [START_REF] Dollfus | High-Frequency Neuronavigated rTMS in Auditory Verbal Hallucinations: A Pilot Double-Blind Controlled Study in Patients With Schizophrenia[END_REF]).

GeodesicSlicer results

After the different steps from GeodesicSlicer described above , the results are described below (see Fig. 1). The 10-20 system EEG electrodes were generated after determining the nasion (x = -0.87, y = 111.20, z =-22.57), the inion (x = 8.30, y = -101.80, z =-37.17), the left pre-auricular (x = -76.31, y = 7.92, z =-42.99) and the right pre-auricular (x = 4.30, y = 22.20, z =-47.49) on the head surface mesh. After placing the cortical stimulation target on the patient's T1weighted anatomical image (x = -40.66, y = -23.26, z = 3.57, corresponding to the superior temporal sulcus), the stimulation target was projected on his head surface mesh. We then obtained the following 3 nearest electrodes around projected stimulation target and their geodesic distance (in cm) with the projected stimulation target: Electrode 1: T3 at 3.65, electrode 2: T5 at 4.47 and electrode 3: C3 at 8.39. These three distances were then used to triangulate and to position the TMS coil manually over the patient's head.

After placing M1 on the patient's T1 anatomical image (x = -30.92, y = -11.77, z =52.58) and choosing the stimulation intensity of the resting motor threshold (rMT) by default (100%), we found the SCDx = 3.34 cm, and SCDm = 3.09 cm and the following two AdjMT (in %

Case 2: Determination of the projected stimulation target in one patient with major depressive disorder

The rTMS can also be used as treatment for MDD that represent one of the most common psychiatric diseases with a prevalence in the general population general of 10-15%. A large number of depressed patients are resistant to drug treatment and, actually, the rTMS is proving to be the greatest therapeutic efficacy [START_REF] Mcgirr | Effectiveness and acceptability of accelerated repetitive transcranial magnetic stimulation (rTMS) for treatment-resistant major depressive disorder: an open label trial[END_REF]). In the case described below, we will consider that patient could, for example, receive a treatment by rTMS in stimulating a key region involved in the MDD: the left dorsolateral prefrontal (DLPFC) [START_REF] Mcgirr | Effectiveness and acceptability of accelerated repetitive transcranial magnetic stimulation (rTMS) for treatment-resistant major depressive disorder: an open label trial[END_REF].

The patient is a 63-year-old woman, diagnosed with MDD (based on the DSM-V, Diagnostic and Statistical Manual of Mental Disorder 5 th edition). She was recruited from the University Hospital (Caen, France) and realized cerebral MRI on a 3T scanner in the same way that in the Case 1.

rTMS site determination

Using the structural MRI, after determining the position of the upper frontal sulcus and the lower frontal sulcus, the target can be defined as being equidistant from the upper and lower frontal sulcus in the coronal plane which crosses the anterior extremity of the temporal pole.

GeodesicSlicer results

The results from GeodesicSlicer are described below (see Fig. 2). The 10-20 system EEG electrodes were generated after determining the nasion (x = -3.65, y = 97.45, z = -5.77), the inion (x = -3.14, y = -89.86, z = -16.90), the left pre-auricular (x = -78.02, y = 2.93, z = -14.58) and the right pre-auricular (x = 71.20, y = 2.88, z = -14.87) on the head surface of mesh.

After placing the cortical stimulation target on the patient's T1-weighted anatomical image (x = -37.62, y = 40.13, z = 55.85, corresponding to the patient's DLPFC), the stimulation target was projected on his head surface (x = -54.22, y = 49.50, z = 66.96). Then, we obtained the following 3 nearest electrodes around projected stimulation target and their geodesic distance (in cm) with the projected stimulation target: Electrode 1: F3 at 2.30, electrode 2: F7 at 3.73 and electrode 3: C3 at 4.03. With this these three distances, after triangulation, it could be possible to position the TMS coil manually over the patient's head.

After placing M1 area on the patient's T1 anatomical image (x = -32.49, y = -3.51, z = 71.66) and have choose the stimulation intensity of the rMT by default (100%), we found SCDx = 22.06 cm, and SCDm = 25.37 cm and the following two AdjMT% (in % stimulator):

Accuracy of the measure

The limits of agreement of Bland-Altman plots of the nasion-inion distances were from -0.82 to 1.96 cm, with a mean difference of 0.57 ± 0.70 cm (range, -0.64 to 1.59 cm), and from -0.52 to 1.82 cm for the tragus-tragus with a mean difference of 0.65 ± 0.59 cm (range, -0.16 to 1.53 cm), showing that the reliability between the geodesic distances was consistent between those calculated by GeodesicSlicer and the manual measures.

Discussion

Personalized and guided stimulation [START_REF] Fox | Identification of reproducible individualized targets for treatment of depression with TMS based on intrinsic connectivity[END_REF][START_REF] Sommer | Toward personalized treatment of hallucinations[END_REF] seems to be an efficient way to improve the discrepant efficacy results previously reported in NIBS studies (Briend et al. Under Review). In this context, we developed GeodesicSlicer, a novel MRI-guided method using individual brain imagery to position the TMS coil reliably on the participant's head. GeodesicSlicer aims to become a complete and easy-to-use toolkit for stimulation researchers by providing accurate target according to 10-20 system EEG electrode positioning and a correction factor to adjust the rTMS dose for the treatment.

In comparison with neuronavigation system, the current method may have some advantages because it is relatively inexpensive and does not require any additional experimental setting. In addition, the MRI-guided method could be particularly useful for therapeutic protocols, because the result of the software allows easy targeting of the same stimulation site across multiple sessions and in multicenter trials, as the three distances can be used offline to position the TMS coil [START_REF] Andoh | A triangulation-based magnetic resonance image-guided method for transcranial magnetic stimulation coil positioning[END_REF]). Moreover, it only needs 3D Slicer, a software freely available online and is relatively user-friendly.

Previous studies have calculated geodesic distances between scalp landmarks using surface mesh representations of the head in order to guide NIBS [START_REF] Andoh | A triangulation-based magnetic resonance image-guided method for transcranial magnetic stimulation coil positioning[END_REF]). Others have worked on the head surface mesh with vectors linking key anatomical landmarks drawn on the mesh and used it to calculate the precise distances on the scalp corresponding to these vectors [START_REF] Vaghefi | MRI Guided Brain Stimulation without the Use of a Neuronavigation System[END_REF]. Just one study has used a semi-automatic approach to generate the 10-20 system EEG correlates external skull locations [START_REF] Xiao | Semiautomatic 10/20 Identification Method for MRI-Free Probe Placement in Transcranial Brain Mapping Techniques[END_REF]), but no study has combined triangulation system based on individualized 10-20 system EEG that morphed to the head surface mesh reconstruction of participant.

Others teams have developed a feasible low-cost solution to track coil positions during rTMS procedures, but the setup and run of the clinical experiment are time-consuming processes [START_REF] Dayan | 3D-printed head models for navigated non-invasive brain stimulation[END_REF][START_REF] Rodseth | A Novel Low-Cost Approach for Navigated Transcranial Magnetic Stimulation[END_REF][START_REF] Washabaugh | A low-cost system for coil tracking during transcranial magnetic stimulation[END_REF]. Although, our method that adjusts the 10-20 system EEG for the participant's skull size, conversely to other [START_REF] Beam | An efficient and accurate new method for locating the F3 position for prefrontal TMS applications[END_REF], have the benefit of taking into account differences in cortical anatomy or skull sizes between each individual in a quick and low-cost way. This triangulation-based MRI-guided method is an alternative to more complicated and costly stereotaxic targeting paradigms. In addition, this software package could be useful for laboratories that do not have access to neuronavigation system.

Limitations and Future Works

Our module has some limitations caused by the method itself. First, unlike probabilistic targeting techniques, Geodesic Slicer requires brain imaging. Indeed, MRI scanning is expensive and not always available in many institutes [START_REF] Xiao | Semiautomatic 10/20 Identification Method for MRI-Free Probe Placement in Transcranial Brain Mapping Techniques[END_REF], but this personalized method using MRI with the participant's own anatomical or functional images to guide is very accurate [START_REF] Briend | A new toolbox to compare target localizations for non-invasive brain stimulation: An application of rTMS treatment for auditory hallucinations in schizophrenia[END_REF][START_REF] Kraus | Projecting Navigated TMS Sites on the Gyral Anatomy Decreases Inter-subject Variability of Cortical Motor Maps[END_REF][START_REF] Lahti | Making Progress Toward Individualized Medicine in the Treatment of Psychosis[END_REF][START_REF] Sommer | Toward personalized treatment of hallucinations[END_REF].

Second this method faces problems in time cost, mainly due to the manual positioning of the international 10-20 system EEG and the measurement procedure [START_REF] Xiao | Semiautomatic 10/20 Identification Method for MRI-Free Probe Placement in Transcranial Brain Mapping Techniques[END_REF].

However by its facilitating approaches for rTMS target localization, the International 10-20 system of EEG represents the gold standard in clinical uses (Uwe [START_REF] Herwig | Using the international 10-20 EEG system for positioning of transcranial magnetic stimulation[END_REF].

Third, as Andoh and colleagues already mentioned [START_REF] Andoh | A triangulation-based magnetic resonance image-guided method for transcranial magnetic stimulation coil positioning[END_REF]), the current MRI-guided method is not a neuronavigation system and therefore cannot provide online monitoring such as a real-time control of the coil angle. However, it is noteworthy that the influence of the coil angle remains debated [START_REF] Niyazov | Functional magnetic resonance imaging and transcranial magnetic stimulation: Effects of motor imagery, movement and coil orientation[END_REF]. Despite these limitations, there are several areas for future work. In particular, some steps of the MRI-guided method workflow under 3D Slicer can be automated as the landmarks placement including the nasion, the inion and the left and the right tragi. Moreover, the path precision calculated by the Dijkstra's algorithm produces "jagged" lines depends on the length of the triangle edges determined during the mesh creation. Further development of GeodesicSlicer should include these issues to resolve this possible variability, by adding a smoothing procedure that fitted the Dijkstra's algorithm path [START_REF] Vaghefi | MRI Guided Brain Stimulation without the Use of a Neuronavigation System[END_REF].

Moreover, we did not assess the effectiveness of this costless method against neuronavigation on the one hand and on the other hand about the effectiveness of correction factors to adjust rTMS dose. The goal of this paper was to test the measurement accuracy of this software, but further studies with larger samples will be necessary to determine the clinical benefit of this method. Although, one paper that use triangulation-based MRI-guided method [START_REF] Andoh | A triangulation-based magnetic resonance image-guided method for transcranial magnetic stimulation coil positioning[END_REF]) has already show the accuracy of it.

Conclusions

In summary, Geodesic Slicer is an alternative to the time-consuming process of neuronavigation system. The triangulation-based MRI-guided method presented herein provides a reproducible and inexpensive way to position the TMS coil that may be used in case of unavailability of online neuronavigation, for instance, in a clinical setting. This MRI-guided method can use cortical landmarks from all MRI scans. Representation of patient's head surface mesh generate by GeodesicSlicer with the nasion and the left pre-auricular in turquoise, notably placed to generate the 10-20 system EEG in red. The stimulation target is then projected onto the head surface mesh (in blue) and localized near these three nearest electrodes with their geodesic distances: Here, the electrodes F3, F7 and C3. In brain stimulation session, these three distances could be used to triangulate and to position the TMS manually over the patient's head. 2. Views sagittal, coronal and axial of the cortical stimulation target (in red). 3. Views sagittal, coronal and axial of the M1 area (in green).

  stimulator): according to Strokes (Stokes et al. 2007): 107.22 and according to Hoffman (Ralph E. Hoffman et al. 2013): 98.75.

  according to Strokes (Stokes et al. 2007): 90.73 and according to Hoffman (Ralph E. Hoffman et al. 2013): 79.88.

Figure with captions Fig. 1 .

 captions1 Figure with captions

Fig. 2 .

 2 Fig. 2. Results of GeodesicSlicer in a patient with major depressive disorder. 1.Representation of patient's head surface mesh generate by GeodesicSlicer with the nasion and the left pre-auricular in turquoise, notably placed to generate the 10-20 system EEG in red. The stimulation target is then projected onto the head surface mesh (in blue) and localized near these three nearest electrodes with their geodesic distances: Here, the electrodes F3, F7 and C3. In brain stimulation session, these three distances could be used to triangulate and to position the TMS manually over the patient's head. 2. Views sagittal, coronal and axial of the cortical stimulation target (in red). 3. Views sagittal, coronal and axial of the M1 area (in green).

Acknowledgement:

The authors would like to thank Drs A. Lasso and K. Yoshimi as well as A. Nourry for their valuable help in 3D Slicer and VTK library, and William P. Armstrong for the English rereading.

Information Sharing Statement: GeodesicSlicer has been put into a toolbox and can be download as an extension of 3D Slicer or directly from Github.

The available implementation was uncoupled from the ethics protected image data used in the case studies. One brain imaging data example can be downloaded here. However, clinical data might be obtained upon request by contacting the corresponding author.

Conflict of interest:

The authors have no conflict of interest to declare.

Funding source: This work was supported by the French Health Ministry (Programme Hospitalier de Recherche Clinique), the Fondation Fondamentale, the Association Perceneige, the Region Normandie and the University Caen Normandie.


Dr Luigi Liquori
email: luigi.liquori@inria.fr

This document examines the use of IoT technology in contact tracing and introduces the concept of Asynchronous Contact Tracing (ACT). ACT identifies contacts with IoT connected objects that have been contaminated by the SARS-CoV-2 virus and works in synergy with solutions designed for manual and digital contact tracing to identify and alert people who may have been infected by the virus.

This shifts the paradigm from synchronously tracing the contacts of the people infected by Covid-19 to asynchronously tracing of contacts of materials (such as infected surfaces, waste-water, airconditioning filters, etc.) that are hosting the SARS-CoV-2 virus.

. The level of contamination depends on the nature of the surface the concentration of the virus, the ambient temperature, the season of the year, the level of humidity, and exposure to sun light. The period of contamination can span from a few hours to several days.

The ACT process uses existing, ready-to-market IoT-based technology and well-established wireless network techniques. The process is not dependent on achieving a certain number of tests, or of people adopting it, in order for the results to be useful. Moreover, it does not require the transmission of any personal information by the user, thus respecting both EU GDPR and public sensibility to personal privacy.

 or the Law of Parsimony (Latin: Lex Parsimoniae), that states that entities and theories useful to solve a problem should not be multiplied unless necessary. On the contrary, simpler entities and theories are preferable to more complex ones because they are easier to test and more likely to be true.

use, as opposed to (Synchronous) Contact Tracing protocols [START_REF]ROBERT protocol, v1.1, Candidate proposal for PEPP-PT (Pan European Privacy-Preserving Proximity Tracing)[END_REF][START_REF]DP3T (Decentralized Privacy-Preserving Proximity Tracing)[END_REF][START_REF]PACT (Private Automated Contact Tracing[END_REF][START_REF]Privacy-Preserving Contact Tracing[END_REF] which tend to be employed when society is put on an urgent, war footing in reaction to an acute problem 1 [START_REF] Cédric | Discours prononcé à l' Assemblé nationale, 27 Mai[END_REF].

The ACT process is not only applicable to the current pandemic wave. The parameters can be adapted to any other virus in a future pandemic.

ACT will be socially and economically acceptable to people who should view it as a National Service (i.e. offered by, the Public Health Authorities) and not be viewed as an obligatory requirement.

ACT will promote individual testing only in the unfortunate event of the user receiving official notifications that he/she may be potentially at risk. It can be applied to all the contexts where people share the same physical space, such as a supermarket, schools, restaurants, hotels, gyms, offices, working plants, hospitals, hospices, etc. It can also be applied to an object that is encountering people while it is in movement, such as a public transportation network.

ACT tests materials that have been in contact with humans and uses wireless and IoT technologies to notify a particular contamination event to Public Health Authorities that, in turn, will inform users who were potentially close to those infected areas.

The process is intrinsically asynchronous because it does not require people to be in the same place and at the same time, and, more importantly, it does not require any information exchange between humans, since the virus has been detected by Group Testing on materials and not on humans.

For many communities, this type of tracing will support an elaborate form of selective lockdown, i.e. the surgical closure of specific areas following a forecast announcing a new spike of infection. It is without doubt a process that will naturally benefit many social and industrial organizations, cities, tourism, education, commerce, and travel etc.

Introduction Lockdown

The main weapon against Covid-19 currently available is social distancing and it requires full/partial/selective lockdown of cluster areas in the event of a high peak in the spread of the virus. This creates severe damage to an economy and to the personal life of its citizens. While we are all waiting for a vaccine, the other relevant tool in the fight against the virus is testing. Unfortunately, widespread testing of large populations in a very short time remains unpracticable. A single testing would require some 1ML tests per week and about 60 weeks for a nation such as France. This is clearly an unworkable solution to the problems raised by the current pandemic and is common to most countries across the world.

It is also clear that many people are unwilling to be Covid-19 tested for social reasons, such as job restriction, economic consequences, violation of private life, or even fear of quarantine. It is now well understood that without widespread testing of the population, the only weapon against Covid-19 is lockdown and subsequent severe economic and social disruption.

Contact Tracing

Contact Tracing (CT) has been actively used in Europe since the 16th century to contain epidemic disease. The principles remain the same today whether carried out by phone, mail or personal contact. The aim is to identify the origin of the infection and to where, or to whom, it has been transferred. If receipt of this information is followed immediately by isolation, treatment and aggressive decontamination, it can lead to containment and the gradual elimination of the disease itself. This year has seen an explosive demand for information about Covid-19. It rapidly focused on the possible use of new technology, and particularly on the capabilities of the mobile phone to automate the process of 'track and trace', producing more accurate and timely information flows -to the advantage of public health, governments and the patient and their personal contacts. The potential of 5G and AI seemed to beckon to a future where silent, rapid transmission of data would protect us all from an unseen virus that observes no political or social boundaries.

There was an immediate counter reaction -this degree of protection might also be intrusive. It threatened the security of our private lives, our right to move freely and our right to confidentiality, particularly in the medical field. There have been many public statements asserting that current practice of CT observes Human Rights as expressed by the EU, but the resulting focus on the individual subject of CT has led to much confusion as to the purpose and direction of the resulting data collection.

Who Needs to Know? And What do they Need to Know?

These are questions that initially were considered to have obvious answers. Design of ICT systems for CT would follow the traditional, manual methodology but should make it faster, more accurate, more useful. Solutions were found to protect the individuals' personal rights and freedoms, but little attention was paid to the 'back-end' -the ultimate destination of the data and the use to which it would be put.

Put simply, it has always been true that security and transparency share a trade-off. One threatens the other. If the potential patient, the user of the 'track and trace' enabled phone is to be fully protected then the program must not store their name or contact details in any way that could lead to their identification. This debate has distracted attention from an important issue. Track and trace exercises in the past have been used principally to support government and public health initiatives not just to inform the individual as to their potential risk. Our modern focus on the rights of the individual threatens the management of track and trace data, which could seriously minimize its usefulness to central authorities which we hold responsible for containing the infection.

Testing for the SARS-CoV-2 virus

Testing is essential to the fight against the SARS-CoV-2 virus and it must be accurate. At the present moment, there are different products based on two methodologies, namely virologic tests (e.g. RT-PCR) and antibodies tests. Unfortunately the accuracy of the tests is not absolute, and this lead to test repetitions and especially to interpret inconsistent results, e.g. VIRAL=POSITIVE and ANTIBODY=NEGATIVE.

There is a further problem caused by the number of commercial products available with varying standards of accuracy as well as the high cost of some testing kits. e.g. see [START_REF]Guider votre choix parmi les 253 tests Covid-19 du marché[END_REF].

Today, as we write, there is a plethora of test solutions (ex: TUL, TDR, TROD, PCR ...). This "abundance" is neither helpful nor economically sustainable. Some of the tests require access to specific processing equipment which may be expensive or not generally available in local hospitals or laboratories. The last frontier is the Rapid Antibody Test (< 30 minutes) intended for the use as an aid in identifying individuals with an adaptive immune response to SARS-CoV-2, indicating prior infection. This rapid lateral flow test is intended for the professional use in laboratory and near patient-testing environment and qualitatively detects IgM and IgG specific to SARS-CoV-2 in serum, plasma, and whole blood.

Dorfman's Group Testing Process

During the 2nd World War Professor Dorfman invented a group testing method to screen for syphilis amongst American soldiers. The testing process we aim to use is inspired by the principle of Dorfman's Group Testing described in his seminal paper, "The detection of defective members of large population" [START_REF] Dorfman | The detection of defective members of large population[END_REF]. The rationale of Group Testing is simple. A sample of blood is taken from, say, five people. These samples are mixed together and one test for the virus is made on the combined fluids. If this test result is negative then all five individual samples are considered to be negative, thus saving four test kits. If the result is positive then the original five samples must be tested individually.

The testing methodology applied in the ACT process is a subset of the Dorfman Group Testing, namely, we are interested only in the presence or absence of the SARS-CoV-2 virus in order to make the ACT process work; the reason for this is trivial because we cannot identify the 'infectors' of waste-water, or air-filtering or other materials (the waste-water or the air-filtering act like an hashfunction, i.e. it would be computationally impossible to come back to the infectors. Note that infector here is used to denote either humans or things hosting the virus.

Use of this sub-kind of group testing (sometimes called pooling testing), in combination with modern digital IoT technology, can provide a new and effective 'forecast' for the introduction of 'selective lockdown'; this means that whenever we test an infected location, the ACT technique will communicate to the population a precise geographical area which can be defined immediately as safe or un-safe.

Asynchronous Contact Tracing

Today, Synchronous Contact Tracing (SCT) [START_REF]ROBERT protocol, v1.1, Candidate proposal for PEPP-PT (Pan European Privacy-Preserving Proximity Tracing)[END_REF][START_REF]DP3T (Decentralized Privacy-Preserving Proximity Tracing)[END_REF][START_REF]PACT (Private Automated Contact Tracing[END_REF][START_REF]Privacy-Preserving Contact Tracing[END_REF] protocols are all being designed to exploit the following digital information; namely, physical and spatial proximity, whilst preserving personal anonymity, in order to track the chain of transmission of the virus.

The challenge of Asynchronous Contact Tracing (ACT) is to remove the time and space limitation and to trace only surfaces that have been contaminated by the virus, rather than just the persons potentially infected. The SARS-CoV-2 virus is known to stay viable for a variable time on a hard surface, depending from the nature of the surface, the virus concentration, the temperature, the humidity conditions, the exposition to sun light. This time can vary from few hours to several days, some examples are 2-3 hours on paper, 4 hour on copper, 3-4 days on plastic and steel, 7 days on face masks, and even more in specific climate conditions [START_REF] Awh | Stability of SARS-CoV-2 in different environmental conditions[END_REF][START_REF] Van Doremalen | Aerosol and Surface Stability of SARS-CoV-2 as Compared with SARS-CoV-1[END_REF][START_REF] Kampf | Persistence of coronaviruses on inanimate surfaces and their inactivation with biocidal agents[END_REF][START_REF] Harvey | Longitudinal monitoring of SARS-CoV-2 RNA on high-touch surfaces in a community setting[END_REF].

Typical hard surfaces at risk include items such as a plastic drink bottle, a can of beer, a milk carton, a water tap or electric towel in a public toilet, a park bench, a metro train's doors, seats and buttons, ATM machines etc. In reality this includes everything around us.

ACT is a process that exploits the following digital information:

• The possibility of testing material (waste-water / air-conditioning filters / hard surfaces) using standard, reliable, testing techniques and producing results in reasonable time [START_REF]Spare links about how extracting group testing from waste-water and air-conditioning: About Waste Water[END_REF]. By results here we mean the percentage of SARS-CoV-2 virus on that material without adding information about who infected it. The accent is therefore on the infected object and not on the person. This is an interesting application of Dorfman's Group Testing techniques. Note that we are not interested to go further in the investigation of who (as a person) really is infected: here the focus is on the fact that a given material or location is infected and, as such, it need to be closed to the public and people whose visited those locations or touched those material should be informed and suggested to go fast for a RT-PCR test.

• The result of testing materials are first interpreted by the National Public Health authorities into a precise forecast, which is geo-localized and then widely communicated to the population using WiFi2 and well-established Web techniques. Access to the forecast may be achieved by a mobile application or by querying a Web-accessible data-base.

Note that the user does not need to communicate his position nor his identity or pseudonyms but has just to listen to the forecast or query the Web-accessible data-base.

This is why we consider ACT to be an intrinsically 'asynchronous'.

To resume, the main objectives of ACT are:

• to trace the IoT that "hosts" the SARS-CoV-2 virus instead of the people that got infected by the SARS-CoV-2 virus;

• to inform National (and Inter-National) Health Authorities so that they can establish common

Policies and local, selective, surgical lockdown in a timely and effective manner, and avoid the severe economic, social and psychological disadvantages of full generic lockdown;

• to forecast for the public information on the clinical status of areas which may be infected and new clusters where the virus is spreading.

Two Simple Scenarios for Asynchronous Contact Tracing Scenario I

Let us imagine that Alice is infected without knowing and goes shopping at the supermarket at 10:15am; she touches some items or sneezes or coughs on them or just breathes close to some items. She then leaves the supermarket. The corridors that Alice walked, its caddy and the goods that Alice touched are surely infected by SARS-CoV-2 virus. Now let us imagine that Bob, who is not infected at all, visit the same supermarket corridors, (unfortunately take the same caddy) and touches those same items at 11:00am, unfortunately and involuntarily infecting him by SARS-CoV-2 virus.

The following considerations are now in play:

• Synchronous Contact Tracing Protocols [START_REF]ROBERT protocol, v1.1, Candidate proposal for PEPP-PT (Pan European Privacy-Preserving Proximity Tracing)[END_REF][START_REF]DP3T (Decentralized Privacy-Preserving Proximity Tracing)[END_REF][START_REF]PACT (Private Automated Contact Tracing[END_REF][START_REF]Privacy-Preserving Contact Tracing[END_REF] cannot allow the exchange between Alice and Bob of pseudonyms, because their respective mobile phones and contact tracing apps were not present and active at the same place (let's say less than 1 meter) and at the same time (let's say for at least 15 minutes).

• If Alice is later tested Covid-19 positive, then there are no means by which Bob will ever be informed and yet this is a likely to be a common means of transmission of the virus.

• It is not necessary to test Alice to inform Bob that he was potentially infected: it would suffice to test the items that Alice touched and the corridors that Alice walked, because Bob did probably the same.

• It is not even necessary for Alice to exist as a physical person, because the virus could transfer on to a surface as a result of coming into contact with other contaminated surfaces, such as articles being sent from an infected region.

• If Alice were infected in the same conditions as Bob but two days before in another supermarket, and if Alice had been informed about the potential risk of infection, then she probably would have done a unitary test (RT-PCR) and, if positive, she would not then touch the items in the supermarket because she would have stayed at home, and probably not infecting Bob.

Scenario II

Alice buys chocolate from the supermarket in Paris. The outer packaging is contaminated and she may receive a dose of the virus. That same day, Bob buys the same chocolate in Lyon, from the same supermarket chain as does Marie in Antibes and Mario in Turin the day after. The supermarket chain has adopted the ACT technology. It follows that all four supermarket corridors selling chocolate are measured positive with the SARS-CoV-2 virus during the cleaning that takes place every two hours in those supermarkets. The consequences of those four measures are as follows;

• Isolate the fours corridor, stop sales of chocolate items in the four supermarkets;

• Verify and warn factory, country of origin and logistics.

In conclusion, ACT can be a practical virus-tracing tool, following the supply chain related to chocolate in that supermarket chain.

ACT in action

ACT performs tests on all materials. If the test is positive, then the following actions could be performed:

• Immediate publishing of the test result on a Public Web-accessible Database with consequent immediate publishing of the isolation/closure of the potentially infected location on a Public Web-accessible Database.

The tasks performed by the smart application associated with ACT are to:

• listen and record all the Peripheral Service and memorize their identifier;

• match all of the identifiers sent by the National Control Service, each of one associated with a precise forecast, with the one he/she has recorded, and, in case of non-empty intersection, notify immediately a potential contamination to the owner of the smartphone. At that point he/she may decide (or not) to go to a laboratory for an individual test;

• ask its National Control Service to be informed about the forecast by other foreign Inter-National Control Services.

The tasks performed by the PC Application associated to ACT are to:

• access the forecast produced by the National Control Service;

• with a simple GUI, also to access precise forecast relative to the local area, denoted via suitable coordinates eg. latitude, longitude, and radius.

Functional Units

Asynchronous Contact Tracing process can be explained as follows by introducing the following SmartM2M/oneM2M functional units as described in the following Figure describing the ACT method architecture [START_REF]SmartM2M; Asynchronous Contact Tracing System; Fighting pandemic disease with Internet of Things[END_REF]. • A Peripheral Service, placed in a specific location, such as a supermarket corridor, a public toilet, a metro station, a fitness room etc. has the capability to transmit, using WiFi technology, the following information such as its unique identifier (also known as the MAC identifier BSSID).

• The Smartphone Application, called IoTAntiCovid , is the main digital tool available to the user to monitor the level of contamination. It collects the identifier from the one or many Peripheral Services and periodically compares it the identifier-related forecast published by the National Control Services, with an associated time and forecast, and about potential contaminated areas. It also receives the local forecast when approaching one or many Peripheral Services and will periodically warn the user according to the Public Health Policies received by means of the one or many National Control Service. The IoTAntiCovid application also offers the opportunity to request and receive forecast about other Peripheral Services (not necessary the one he/she was close) information always without communicating its user's location. The periodicity of the request will be configurable by the IoTAntiCovid application.

• The PC Application, called IoTAntiCovid -PC, is an user connected device (tablet, personal computer, smart TV, etc.) that communicate with the National Control Services to query and receive back information about the status of the situation in specific geographical areas. By delegation it can also query other Inter-National Control Services. It basically supports the user about the discover of its own risk according to the past visited areas and/or to plan his movement and behaviours in advance. The purpose of IoTAntiCovid -PC is similar to the Smartphone Application, namely be informed about the forecast of the virus, but its UX/UI is rather different.

• The Detection Service informs the one or many Local Services about the detection of the contamination in one of the proximity areas associated to one or many Peripheral Service. The chemical process and the kind of test used to check the material for contamination (e.g. air-filter, waste-water filter, dirty cleaning tools and water, etc.) are not strictu sensu pertinent to ACT Technical Specification. Each Detection Service can be associated with one or many Peripheral Services and it shall communicate with the one or many Local Services the detection of a contamination. It will exchange its monitoring and configuration information with one or many Local Services, including means to exchange the information necessary to the management of the specific detection tool.

• The Local Service receives the information from one or many Detection Services belonging to it and propagates this information to the National Control Service. It receives the necessary indications, in order to behave according to the National Authorities policies for the one or many Peripheral Services configuration (e.g. for the broadcast frequency of the identifier or for the emission of the forecast information). It will also transmit to the National Control Service the identifiers of the one or many Peripheral Services involved in a positive detection, together with their location and the approximate range of the detection. Finally it will exchange monitoring and configuration information with the National Control Service.

• The National Control Service will receive the information from the one or many Local Services related to contamination and will provide to it information according to the Public Health Authorities policies on the one or many Peripheral Services configuration (e.g. for the broadcast frequency of the identifier or for the emission of the forecast information). The National Control Service also represents a point of interaction with the users. It provides, on request, the identifier of the one or several Peripheral Services announcing contamination to the many Smartphone Application IoTAntiCovid and additional information about indications and suggestions according to the Public Health policies (e.g. the suggestion of avoiding certain areas or to perform a human test verification if certain conditions of exposure are met, according to the detection of the virus in the areas covered by one or many Peripheral Services). The National Control Service will exchange its monitoring and configuration information with the one or many Local Service and with the many Smartphone Application IoTAntiCovid . It will also coordinate and communicate the information about the areas of detection of the virus with one or many Peripheral Services in other Nations or Regions. Finally, the National Control Service will be in contact with others National Control Services in order to exchanges forecast, common Policies so informing users of different Nations, Regions or Municipalities, etc., when they are traveling from home. This enables ACT to be a genuine, international forecast tool.

Conclusion

To resume the main contribution of ACT is to use a sub-kind of group testing, in combination with modern digital IoT technology, to provide a new and effective 'forecast' for the introduction of 'selective lockdown'; this means that whenever we test an infected location, the ACT technique will communicate to the population a precise geographical area which can be defined immediately as safe or un-safe. We are interested only in the presence or absence of the SARS-CoV-2 virus in order to make the ACT process work; the reason for this is trivial because we cannot identify the 'infectors' of waste-water, or air-filtering or other materials. Note that infector here is used to denote either humans or things hosting the virus.

Figure 1 :

 1 Figure 1: ACT Method Architecture

"C'est un projet qui répond à une crise historique sans laquelle il n'existerait pas et au-delà de laquelle il n'existera pas : l'épidémie de Covid-19".[START_REF] Cédric | Discours prononcé à l' Assemblé nationale, 27 Mai[END_REF].

It is a work in progress.

Acknokledgements

The authors acknowledge Dr. Patrick Guillemin (ETSI) for hit technical suggestions and friendly support.


Steve Ndengue

Yohann Scribano

Fabien Gatti

Steve Ndengué

Richard Dawes

MCTDH calculations of 5-atoms inelastic collision State-to-state inelastic rotational cross sections in five-atom systems with the MCTDH method

come

I. INTRODUCTION

Our understanding of energy transfer between colliding atoms and molecules is crucial for several fields like chemical processes in Earth's atmosphere, combustion chemistry processes for industrial and aeronautic applications, and also in outside Earth's environnement, such as in the InterStellar Medium (ISM), stellar and planetary (and exo-planetary) atmospheres. For example, determination of molecular abundances in proto-stellar clouds requires one to solve radiative transfer equations which themselves require knowledge of state-to-state rotational excitation (and de-excitation) rate constants. The most accurate technique to determine accurate rotational inelastic cross sections (and thus rate constants) is certainly the time independent quantum method based on Close Coupling (CC) which solves the nuclear Schrödinger equation. The CC-methodology is a powerful and intuitive method to solve for the rovibrational states and obtain the collisional cross-sections using potential energy function typically obtained by means of accurate ab initio methods. This approach was intensively used in the past for the simulation of rotationally inelastic collisions of nonreactive diatomic and polyatomic molecules with a) dawesr@mst.edu atoms and molecular colliders such as H 2 which is the most abundant molecule in space. The other well-known technique is the Quasi-Classical trajectory method which solves the classical Hamilton's equations of motion, and is usually processed by quantizing the classical action of trajectories using binning histogram methods. This approach is numerically very efficient since computational time decreases as the collision energy increases. However, the QCT approach is best suited for the high temperature regime where quantum effects are expected to be less important, and whereas for low temperature, its validity could be questionable. More recently, the mixed quantum-classical trajectory (MQCT) 1 and statistical adiabatic channel model (SACM) 2,3 have been proposed as alternative approaches to the CC method and they have both shown very promising results in certain cases, but less accurate behavior in other cases, with their domain of validity not yet fully known and needing to be investigated. In this work, we want to address another alternative approach which can have some numerical advantages compared to CC approaches and especially when we consider high-dimensional molecular systems for which the density of internal rovibrationnal state is large leading to unfeasable CC simulations. Here we extend a recently reported MCTDH investigation 4 of Ar colliding with a rigid triatomic rotor (H 2 O) by studying the collision between diatomic rigid rotors with triatomic rigid rotors using the MCTDH methodology. This five-dimensional system corresponds to the description of rigid rotors inelastic rotational excitation, an example of which is H 2 O colliding with H 2 , a system of interest for ISM and astrochemistry. In the past, the spectroscopy of the weakly bounded complex H 2 O-H 2 was intensively studied [5][6][7][8][9] but also the collisional dynamics between those two molecules with various methods. Those methods used Time Independent CC formalism 7,[10][11][12][13][14][15] , as well as quasi-classical approach 16 which make this system a very good five-dimensional benchmark. The MCTDH method has emerged in recent years as one of the methods of choice for the study of the quantum molecular dynamics (spectroscopy or collisional dynamics) for system in high dimensionality. We are particularly focused here on the low-energy collision region in order to most stringently test the method. This is because of the resonance structures mentioned before, and also because of the known difficulty for time-dependent approaches to describe low-energy collisions, since the wavepacket in that case moves slowly, and the propagation then takes a very long time (which may lead to a propagation of errors). The paper is organized as follows: in Sec.II we will first give a brief introduction to the MCTDH method and the Hamiltonian (the kinetic energy operator and the potential energy surfaces) used for the dynamics, in Sec.II D we will describe how the cross-sections are obtained and discuss the convergence of the calculations. In Sec.III we will present our results and their comparison to CC simulations. We will finally conclude in Sec.IV and give some perspectives of this method for more complex molecular systems.

II. COMPUTATIONAL METHODOLOGY

A. MCTDH method

The MultiConfiguration Time Dependent Hartree (MCTDH) [17][18][START_REF]Multidimensional Quantum Dynamics: MCTDH Theory and Applications[END_REF][START_REF] Beck | [END_REF] is a time-dependent method in which each degree of freedom is associated with a small number of orbitals or single particle functions (SPFs) which, through their time dependence, allow an efficient description of the molecular dynamic process. The total MCTDH wave function is expanded in Hartree products, that is, products of single-particle functions

Ψ(Q 1 , . . . , Q f , t)= n1 j1=1 • • • n f j f =1 A j1•••j f (t) f κ=1 ϕ (κ) jκ (Q κ , t) = J A J Φ J , (1)
where f is the number of degree of freedom of the system, Q 1 , . . . , Q f are the nuclear coordinates, A j1•••j f denotes the MCTDH expansion coefficients, and ϕ (κ) jκ (Q κ , t) are the n κ SPFs associated with each degree of freedom κ. The subsequent equation of motion for the coefficients and single particle functions, obtained by applying the Dirac-Frenkel variational principle to the wave function ansatz, are a set of coupled nonlinear differential equations which conserve the norm and the total energy (for time-independent Hamiltonian).

The MCTDH equations of motion are solved at every time step and require to build the mean f ields and the Hamiltonian matrix. As the Hamiltonian matrix and mean f ields require evaluation of integrals of the type

ϕ (1) j1 • • • ϕ (f ) j f | Ĥ|ϕ (1) l1 • • • ϕ (f ) l f . (2)
Because of the considerable time those calculations take in the process, the MCTDH approach relies in transforming the Hamiltonian operator Ĥ to a sum-of-products of single particle operators

Ĥ = s r=1 c r f κ=1 ĥ(κ) r , (3)
which then significantly simplifies the evaluation of the Hamiltonian matrix and the mean f ields ϕ

j1 • • • ϕ (f ) j f | Ĥ|ϕ (1) l1 • • • ϕ (f ) l f = s r=1 c r f κ=1 ϕ (κ) jκ | ĥ(κ) r |ϕ (1)
(4) In the following two subsections we will describe the expression of the Kinetic Energy Operator and the Potential Energy Operator as a sum-of-products of single particle operators.

B. MCTDH Kinetic Energy Operator

We described quite extensively in the previous paper 21 on this system the Hamiltonian (Kinetic Energy Operator (KEO) and Potential Energy Surface (PES)). The Kinetic Energy Operator (KEO) expressed in Jacobi coordinates is used in this work, following the subsystem KEO derivation presented by Brocks et al. 22 and generalized by Gatti and Iung 23 . Here, just as in our previous work 21 on the same system we work in the E 2 frame obtained by rotation of the first two Euler angles of the SF frame. Althrough this section we express the KEO in atomic units.

2 T = - 1 µ R ∂ 2 ∂R 2 + 2 T1 + 2 T2 + 1 µR 2 ˆ J † ˆ J + ( ˆ L 1 + ˆ L 2 ) 2 -2( ˆ L 1 + ˆ L 2 ) ˆ J E2 , (5)
where µ R is the reduced mass of the H

T1 = A 2 L2 1,+ + L2 1,-+ L1,+ L1,-+ L1,-L1,+ BF1 - C 2 L2 1,+ + L2 1,--L1,+ L1,--L1,-L1,+ BF1 +B L2 z BF 1 , (6)
where the rotational constants 6 have values A = 27.8572 cm -1 , B = 14.5145 cm -1 and C = 9.2799 cm -1 . The rigid rotor kinetic energy of the H 2 fragment is written simply as T2 = B H2 ˆ L 2 2 where the rotational constant B H2 = 59.2434 cm -1 . The final form of the KEO obtained after some analytical development and which is implemented in the MCTDH code is then

2 T = - 1 µ R ∂ 2 ∂R 2 + 2 T1 + 2 T2 + 1 µR 2 J(J + 1) + ˆ L 2 1 + ˆ L 2 2 -2 L2 1,z -2 L2 2,z E2 + 1 µR 2 L1,+ L2,-+ L1,-L2,+ -2 L1,z L2,z E2 + 1 µR 2 C + (J, K)( L1,+ + L2,+ ) E2 + 1 µR 2 C -(J, K)( L1,-+ L2,-) E2 , (7)
with

C ± (J, K) = (J(J + 1) -K(K ± 1). (8)
In the previous equations, ˆ L X (X = 1, 2) is the total angular momentum of fragment 1 or 2 and LX,± are their corresponding creating and lowering operators expressed as LX,± = LX,x ± LX,y . In terms of angles, we can write for their expressions in the E 2 frame

L2 1 = - 1 sinβ ∂ ∂β sinβ ∂ ∂β + k α sin 2 β , (9)
L2 2 = - 1 sinθ ∂ ∂θ sinθ ∂ ∂θ + k φ sin 2 θ (10)
and

L1± = ± ∂ ∂β -k α cotβ, (11)
L2± = ± ∂ ∂θ -k φ cotθ, (12)
with an additional shift k → k ± 1 in the equations (11) and (12). In the BF 1 frame (for H 2 O) we have similar expressions with k α changed into k γ and the corresponding additional shift

k γ → k γ ± 1 L2 1 = - 1 sinβ ∂ ∂β sinβ ∂ ∂β + k γ sin 2 β ( 13
)
and

L1± = ± ∂ ∂β -k γ cotβ. (14)

C. MCTDH Potential Energy Operator

We use as in our previous work 21 the same PES of Valiron et al 25 and transform it to the appropriate set of coordinates 21 , but now for scattering calculations. Using the multipolar description of Valiron PES and the transformation suggested by Avoird and Nesbitt 5 , we can have an analytically exact representation of the PES in the set of coordinates used for the dynamics and thus be able to extend the PES to any distance without doing an additional refitting with MCTDH. This analytical representaion is particularly useful as the range of the PES which is significantly extended along the R coordinate in order to describe accurately the low energy region of the cross section. The surface from Ref. 25 after transformation from its original set of coordinates to the E 2 frame can be expressed as

V (R, β 1 , γ 1 , α 1 , θ 2 , φ 2 ) = r β ,rγ rα,r θ Ṽr β ,rγ rα,r θ (R)frβ,rγ rα,r θ (ω 1 , ω 2 ) (15) where frβ,rγ rα,r θ (ω 1 , ω 2 ) = D (r β ) rα,rγ (β 1 , γ 1 , α 1 ) C r θ ,-rα (θ 2 , φ 2 ).
( 16)

D (r β ) rα,rγ (β 1 , γ 1 , α 1 ) is the Wigner D-matrix and C r θ ,-rα (θ 2 , φ 2 )
is the Racah normalized spherical harmonics. As we work in the momentum representation for the degrees of freedom γ, α and φ we replace γ 1 , α 1 and φ 1 by their momentum representation k γ , k α and k φ where we dropped the 1 and 2 indices to simplify the notation.

D. Inelastic cross-sections calculations

The H 2 O+H 2 scattering calculations were designed in a similar way as our previous scattering studies 4,9,26 where we used the MCTDH program: we will repeat those steps here.

(i) First a wavepacket defining the initial state has to be constructed: it is here (since we are in the rigid rotor approximation) the product of the initial rotational state of H 2 O, the initial rotational state of H 2 and a Gaussian function along the dissociative coordinate R, starting far away from the interaction region. The Gaussian's parameters (with a negative impulsion) are selected in order to cover the energy range of interest. The initial wavefunction is thus expressed as

Ψ i (R, β, k γ , k α , θ, k φ )= ψ i (β, k γ , k α , θ, k φ )χ(R) (17) = τ 1,i (β, k γ , k α )τ 2,i (θ, k φ )χ(R)
where ψ i is a product of the initial rotation states of H 2 O (τ 1,i ) and H 2 (τ 2,i ). The initial rotational state of H 2 O are linear combination of Wigner-D matrices usually written in the simplified notation j KaKc where j is the molecular rotational angular momentum and K a , K c the projection on the BF-axis along the prolate and oblate limits. If, for example, we write the Wigner D-matrices as |j, k, m , the rotational state for j = 0 is |0 00 = |0, 0, 0 with eigenvalue 0. The rotational states for j = 1 are 6 :

|1 10 = 1 √ 2 [|1, 1, m + |1, -1, m ], |1 11 = |1, 0, m , |1 01 = 1 √ 2 [|1, 1, m -|1, -1, m ], (18)
with respective eigenvalues A + B, A + C and B + C (clearly with a 2j + 1 = 3 degeneracy) with A, B and C being the rotational constants of H 2 O given above and m = -1, 0, 1 the magnetic quantum number projections. For j ≥ 2, we used the symbolic arithmetic program Sympy 27 and obtained the normalized eigenfunctions and related eigenvalues that helped to generate analytically the initial but also the final rotational functions. While it is obviously not compulsory to obtain the initial and final functions analytically, it nevertheless helps to verify that these results match the known analytical results for this rotational Hamiltonian. The initial state of H 2 is an associated Legendre polynomial that we write simply as l m where l is the molecular angular momentum of the H 2 fragment and m the projection along the BF-axis. The Gaussian χ writes

χ(R) = 1 √ 2πw exp - R -R 0 2w 2 e ip0(r-R0) (19)
with R 0 , p 0 and w being respectively its center in coordinate and momentum space and its width. For this work, we tested various parameters and two sets of energy distributions (for which a state-to-state cross section will be presented later) are shown in Figure 2.

When the initial rotational states of H 2 O or H 2 have a non-zero angular quantum number (j = 0 in j KaKc or l = 0 in l m ), then calculations for each magnetic quantum number projection have to be done also in order to obtain the state-to-state cross-sections from a specific initial rotational state. The lowest para and ortho rotational levels of H 2 O and H 2 are presented in Table I. The 2j +1 (and 2l+1) repetition of the calculations is avoided

(x) = η|x -x c | b Θ(x -x c )
, where Θ is the Heaviside step function, x c the starting point of the CAP was selected such that the length of the CAP is 20 bohrs, η its strength is 4.76×10 -9 and b its order is 4.

(iii) We then have to select a primitive and SPF bases for the calculation. The primitive basis, its range and the number of the SPFs used for the wavepacket propagation are summarized in Table II. A primitive basis composed of Fast Fourier Transform (FFT) functions for the intermolecular distance R was coupled with a Wigner-DVR 24 and a two dimensional Legendre-K DVR 28 to describe respectively the orientation of the H 2 O and the H 2 fragments in the E 2 frame.

(iv) As we are interested in the low-energy region of the cross-sections, we need to propagate the wavepacket for a very long time. In order to control the propagation and the convergence of the calculation, we allowed the wave packet to propagate for a maximum of 200 ps (with the wavefunction printed every 20 fs) or when 99.0% of the wavepacket has been absorbed by the CAP. For higher values of J the propagations completes very quickly: for example, for J = 35 the calculation stops after about 10 ps while for J = 0 a complete 200 ps of propagation is needed and still by that time only about 98.5% of the wavepacket has been absorbed by the CAP. Once the propagations are completed, two procedures can be used to obtain the various state-to-state transition probabilities: it is either the flux [START_REF] Beck | [END_REF]29,30 approach or the Tannor and Weeks 31 method. We tested both methods but we will mainly present results from the Tannor and Weeks method as both results were similar. Unless specified otherwise, in the following equations and text E = E tot that is the total energy in the collisional process. The kinetic or collision energy will be labelled E kin . Using the Tannor and Weeks method, the state-to-state reaction probabilities are expressed as 31,32

P J f ←i (E) = |S J f i (E)| 2 (20) = 1 4π 2 |∆ i (E)| 2 |∆ j (E)| 2 T 0 e iEt C f i (t)dt 2 ,
where S f i is the S-matrix element, C f i is the autocorrelation function obtained in term of the propagated function Ψ i and the final state function Ψ f ,

C f i (t) = Ψ f |e -iHt |Ψ i . ( 21
)
where ∆ i and ∆ f are the energy distribution 33,34 of the initial and final wavepackets, respectively. The functions Ψ i and Ψ f are expressed as in Eq. ( 17) where we can rewrite the labelled transition i ← j as

j 1 ζ 1 m 1 , j 2 l 2 ← j 1 ζ 1 m 1 , j 2 l 2 . ( 22
)
Here j 1 (j 1 ') and j 2 (j 2 ') are respectively the orbital angular momentum quantum numbers of the H 2 O and H 2 molecule, ζ 1 (ζ 1 ) and m 1 (m 1 ) are respectively an index along the BF-axis which identifies specific rotational states and the Space-Fixed projection of H 2 O running from -j 1 to j 1 . l 2 (l 2 ) is the Space-Fixed projection of H 2 running from -j 2 to j 2 . In the literature, ζ 1 (ζ 1 ) is often replaced by K a K c where K a and K c are respectively the projection of the molecular orbital angular momentum along the Body-Fixed (Molecular) axis in the prolate and oblate limits.

As mentioned in a previous publication 4 , while the flux and Tannor and Weeks methods lead to similar probabilities, the flux method usually produces a small, but non-negligible flux in an energetically forbidden region (at energies lower than the total energy threshold). This numerical issue can easily be corrected by zeroing the probabilities in the energetically forbidden region when summing the cross-sections (Eqn. ( 23)). Conversely, the Tannor and Weeks method requires a broad energy distribution of the final wavefunctions (wavefunctions on which the propagated wavefunction is projected) to make the correlation function disappear more quickly. Following the calculation of the state-to-state transition probabilities, a weighted sum (by the 2J BF + 1 factor) of these probabilities then produces the inelastic scattering crosssection from the relation

σ j 1 ζ 1 j 2 j1ζ1j2 (E) = π 2 2µ R (2j 1 + 1)(2j 2 + 1)E kin × j1 m1=-j1 j2 l2=-j2 j 1 m 1 =-j 1 j 2 l 2 =-j 2 Jmax J BF =0 × (2J BF + 1)P J BF j 1 ζ 1 m 1 ,j 2 l 2 ←j1ζ1m1,j2l2 (E), (23)
where

E kin = E -E int = E -H2 -H2O with E int , H2
and H2O the internal initial energy of the whole system and the initial rotational energies of the H 2 and H 2 O fragments respectively.

III. RESULTS AND DISCUSSIONS

Before presenting and discussing the results, it is worth pointing out that the purpose of this work was not simply to demonstrate precise agreement between the MCTDH results and those obtained with the Time-Independent Close-Coupling method (if any could be found), but rather to highlight the computational efficiency with which MCTDH makes it possible to obtain accurately the cross section with reasonably limited resources. However, as the following results will show, we were able to obtain overall a good agreement with Close Coupling results. For all of the calculations reported in this work we used a Primitive and Single Particle Functions basis as reported in Table II. Between 1,000 to 6,000 Single Particle Functions were used to describe the dynamics of the system compared with more than 1,000,000 Primitive Basis functions that would have been needed overall for a standard wave packet calculation.

A. Convergence of the calculations

The convergence of the MCTDH calculations depends on various parameters: the primitive basis, the SPF basis, the range of the propagation, the intensity of the CAP and the duration of the propagation. For this work we used a primitive basis with angular parameters similar to the one used by Scribano et al 15 : the number of DVR points for the radial degree of freedom was selected to be sufficiently dense for the range selected.

Convergence: radial coordinate range

To determine the range of integration, we performed 2 scattering calculations: one going from 2 to 52 bohrs (with 768 radial DVR points) and the other from 2 to 72 bohrs (with 960 radial DVR points), that is with about the same density of points along the radial coordinate for the 2 calculations, but also the same number of SPF (1,000 SPF for these tests). Those 2 calculations were compared with the available close-coupling calculations (considered exact) as depicted in Figure 1. As expected, extending the grid and displacing the position of the CAP (which however keeps the same length) improves the quality of the results. Also, the contraction of the grid induces a very small blue shift in the position of the resonance structures of the cross-section. We computed the transition probabilities for each individual J from 0 to 35 (the maximum J value considered for the time independent close coupling calculations) in parallel using )

E total (cm -1 ) MCTDH -R max = 72 MCTDH -R max = 52 MOLSCAT -CC
32 processors for each J-calculation. These calculations with 1,000 SPF lasted from less than a day for the higher J-values to about 7 days for the lower J-values which require a longer propagation.

Convergence: the energy distribution

An important aspect of the calculation is the choice of the initial wavefunction. As stated before, the width and initial momentum of the Gaussian are connected to the energy range that will be covered by the scattering calculation. We perform 2 separate scattering calculations to see the influence of the final cross-sections on the choice of the energy range covered by the scattering. For the first calculation we selected a Gaussian centered at 52 bohr with a width of 1.2 a.u. and an initial momentum of -1.6 a.u. (where the negative sign emphasizes the fact that the initial wave packet goes towards the interaction region) and for the second calculation, the Guassian is centered at 52 bohr but with a width of 1.0 a.u. and an initial momentum of -2.1 a.u. The energy distribution of the 2 calculations at J = 0 are displayed in Figure 2 and a comparison between the cross-section for the 0 00 → 1 11 transition are compared with the time independent close coupling calculation on Figure 3. The comparison shows that the calculation with the shorter energy distribution range describes better the low energy than one with the larger range. Therefore, for calculations spanning a large collisional energy range, it would be necessary to perform different sets of calculations with two or more energy distribution ranges in order to cover accurately a wide range of collision energies. An advantage of the MCTDH method is that the approach, just like the Close Coupling method, is variational. Therefore, we know that once a suitable choice of a primitive basis is made, we can improve the accuracy of the calculation by further increasing the size of the SPF basis. A characteristic of the MCTDH method is that it is designed to reproduce very well the most intense fea- tures of a process with a limited number of functions: however, the less intense features can be recovered easily by increasing the size of the Single Particle Functions. In the limit N SP F → N P BF , we can recover the exact numerically exact calculations, where N SP F and N P BF are the total number of Single Particle and Primitive Basis Functions respectively. We were able to verify this feature by looking specifically at the 0 00 → 2 11 transition with an increasing SPF basis, and comparing the results with the exact Time Independent results. One aspect of that transition is that it is about 50 times less intense than the 0 00 → 1 11 transition which is the most intense one, and also retrieved from the same set of calculations. The comparison is shown in Figure 4. The parameters of the 4 basis are presented in Table III for the 3 combined modes and are such that the larger combination of SPF basis is the one mentioned in the Table, but for tion a higher J, the combination of numbers is likely to be smaller. For the final results reported in this work, we had to decide on convergence criteria. As we had to do several propagations (for each J), and since the SPF basis size decreases as J increases, we chosed an SPF basis such that for each (combined) mode, the weight of the last SPF (arranged in decreasing order of weight) should ) E total (cm -1 ) mctdh molscat be less than 10 -5 for the calculations. This uniform criteria allowed to varying basis which will however maintain a consistent accuracy of the calculation.

B. Rotational excitation cross sections

The calculated inelastic scattering cross-sections for the transitions from the ground state to the first excited states, i.e. 0 00 → 1 11 , 0 00 → 2 02 and 0 00 → 2 11 are displayed on Figures 5, 6 and 7 respectively. The figures show a very good agreement between the MCTDH calculations and the CC calculations at these low energies. In particular, the position of the resonances are well reproduced in each calculations and an even better agreement will certainly be met with a larger SPF basis. The position and the intensity/width of the peaks are connected in the calculations to the range and the duration of the propagation respectively. A too-narrow range of the propagation shifts the peaks to higher energies because of confinement effects and a shorter duration of the propagation would make the resonances too broad or even nonexistent. In fact, we actually see in the course of the propagation that some resonances start to appear only after a specific amount of propagation time has elapsed. Similarly, just like with a rovibrational calculation, increasing the size of the (SPF) basis improves the agreement on the peaks position with the time independent calculations which are considered here to be numerically exact.

Similarly to the calculations involving the para-states, we perfomed calculations for the ortho-states. Here however, as we mentiond before, we used the same basis sets (made of both para and ortho rotational states ) E total (cm -1 ) mctdh molscat of H 2 O) to run the propagation and obtain the crosssections. The excitattion cross-sections for 2 ortho transitions 1 01 → 1 10 and 1 01 → 2 12 are displayed respectively in Figures 8 and9. Here again, we observe a very good agreement between the MCTDH and the time independent results, both with respect to the position and the intensity of the peaks. We notice however some small issues at the start of the cross-sections where some peaks are displayed with a smaller intensity than the time independent results, an issue that is likely linked to the ) E total (cm -1 ) mctdh molscat duration of the propagation. Also towards the end of the spectrum, it appears that the MCTDH calculations have a higher amplitude than the time independent an issue that would be corrected if a wider energy distribution appropriate to describe higher energy values were to be selected. We performed additional calculations to determine the state-to-state cross-sections from the first excited parastate. This allowed to obtain a de-excitation cross-section 1 11 → 0 00 in addition to an excitation cross-section 1 11 → 2 02 . These cross-sections are displayed in Figures 10 and 12. Figure 11 is similar to Figure 10 with however a logarithmic scale used for the y -axis (crosssections intensity). We observe as in the previous cases that the excitation cross-section reproduces quite well the position of the peaks/resonances of the spectrum; this is also the case for the de-excitation process until we get to very low kinetic energy values. However, while the intensities are usually quite well reproduced for the excitation processes, the logarithmic scale cross-sections (Figure 11) highlights how difficult it is to accomplish it for the de-excitation process. This difference in accuracy between the 2 types of calculations (excitation and deexcitation) using the same SPF basis can be understood from the cross-sections equation (Eq. ( 23)) but also from the calculation of state-to-state transition probabilities in the flux formalism. Within this formalism, the transition probabilities are obtained as a ratio of the flux going towards a specific channel to the total flux. When the calculation is exact, the total flux is well approximated by the energy distribution of the propagation. While for the excitation processes the starting point of the state-to-state cross-sections is often in an energy region where the absolute value of the energy distribution is quite significant, for de-excitation processes, this often corresponds to an energy region where the energy distribution is very small, 2 order of magnitude or more lower than the peak of the energy distribution. As the ) E total (cm -1 ) mctdh molscat flux to that channel at low kinetic energy is also very small, it is quite difficult to describe with a small contracted basis designed for the whole energy range, the small features that mainly influence the low kinetic energy domain. Additionally, the 1/E c term amplifies the disagreement between the MCTDH and time independent calculations at low collision energies where its role is alleviated as energy increases as we can see from the picture. However, we realized that by increasing the SPF FIG. 13. Comparison of MCTDH calculations-with the smaller (old) and larger (new) basis-with time independent close coupling results for the 111 → 000 rotational transition of H2O. Note the logscale applied to the y -axis. ) E total (cm -1 ) mctdh-new mctdh-old molscat basis size we were able to significantly improve the comparison with the time independent results. For example, while for most of the calculations displayed in this work we used a SPF basis with parameters described in Table II, for the 1 11 → 0 00 we used up to 35 SPFs for the H 2 O angular basis and up to 25 SPFs for the FFT on the R coordinate on some J propagation that we identified as the most problematic. Despite this significant increase, we noted some (modest) change on the low energy region of the de-excitation cross-section for the 1 11 → 0 00 transition while only a modest to nonexistant one on the excitation cross-section for the 1 11 → 2 02 transition as we can see on Figures 13 and 14 respectively. So essentially we only observe a problem in the de-excitation cross-section of the 1 11 → 0 00 transition within 10 -20 cm -1 above the treshold energy, i.e. here, below ≈ 55 cm -1 of total energy E.

IV. CONCLUSION

We performed the first benchmark 5D calculations for the inelastic scattering of an asymmetric top molecule (H 2 O) with a linear molecule (H 2 ), all in the rigid rotor approximation, using the MultiConfiguration Time Dependent Hartree procedure and without any additional approximation on a recently constructed PES. We find as expected a very good agreement between the time independent close coupling results and our calculations which are done with a relatively modest amount of resources.

Time Independent Close Coupling calculations are frequently regarded as the benchmark for these type of atom/molecule or molecule/molecule collisional pro- cesses. However, they quickly become untractable when either the density of rotational states of the partners or the dimensionality of the collisional system increases. It then becomes essential to look for alternate means of computing collisional cross-sections. We anticipate that the MCTDH procedure (and its variants ML-MCTDH, ...) could be one of those alternate methods to describe fully quantum mechanically the collisional dynamics in these challenging cases. As the MCTDH inelastic scattering calculations usually generate cross-sections for all transitions starting from an initial state, they will also be very relevant for comparison with experimental results where quite often only selected transitions could be measured.

With respect to astrophysical applications where one is often interested in the state-to-state rates, instead of the specific state-to-state cross-sections, one can infer the state-to-state rate of a de-excitation process (where the cross-sections calculations appear to be challenging at low kinetic energy using the MCTDH method) from the reverse excitation process, which is known as detailed balance.

The current MCTDH implementation of the Wigner-DVR does not yet enable even/odd symmetry specified calculations. Hence it was not possible to selectively perform wavepacket propagation calculations for the para and ortho cases which have a specific symmetry with respect to the H 2 O fragment axis. We thus did not take advantage of the primitive basis size reductions that would occur from the symmetry consideration. Nevertheless, the calculations for this system are still tractable with respect to both the time and memory requirements.

The usage of the MCTDH method for this type of ap-plications is still at a development stage. However, we are working on defining a scheme for which these type of computations are likely to become standard and straightforward for people interested in such applications. Imposing the symmetry for the Wigner-DVR in the MCTDH package will be part of this process and will certainly improve the results.

FIG. 1 .

 1 FIG. 1. Comparison of MCTDH calculations for Rmax = 52 and Rmax = 72 with time independent close coupling results for the 000 → 111 rotational transition of H2O.

FIG. 2 . 2 FIG. 3 .

 223 FIG.2. Energy distribution used for the MCTDH calculations at J = 0. edist1 is generated with a Gaussian of width 1.2 a.u. and an initial momentum of -1.6 a.u. and edist2 with a Gaussian of width 1.0 a.u. and momentum of -2.1 a.u.

FIG. 5 .

 5 FIG. 5. Comparison of MCTDH calculations with time independent close coupling results for the 000 → 111 rotational transition of H2O.

FIG. 6 .FIG. 7 .

 67 FIG.6. Comparison of MCTDH calculations with time independent close coupling results for the 000 → 202 rotational transition of H2O.

FIG. 8 .FIG. 9 .

 89 FIG. 8. Comparison of MCTDH calculations with time independent close coupling results for the 101 → 110 rotational transition of H2O.

FIG. 10 . 2 )E

 102 FIG.10. Comparison of MCTDH calculations with time independent close coupling results for the 111 → 000 rotational transition of H2O.

FIG. 11 .FIG. 12 .

 1112 FIG. 11. Comparison of MCTDH calculations with time independent close coupling results for the 111 → 000 rotational transition of H2O. Note the logscale applied to the y -axis (cross-section intensity).

FIG. 14 .

 14 FIG.14. Comparison of MCTDH calculations-with the smaller (old) and larger (new) basis-with time independent close coupling results for the 111 → 202 rotational transition of H2O.

TABLE I .

 I Rotational energy levels (cm -1 ) or ortho and para H2O and H2. The rotational constants used to generate those levels are A=27.8572 cm -1 , B=14.5145 cm -1 and C=9.2799 cm -1 for H2O and BH 2 = 59.2434cm -1 for H2. Ortho levels of H2 are displayed in italic.

		H2O				H2
		ortho		para	ortho & para
	101	23.7995	000	0.00000	0	0.0000
	110	42.4024	111	37.1583	1	118.4868
	212	79.5133	202	70.1329	2	355.4604
	221	135.3220	211	95.2454	3	710.9208
	312	173.5976	220	136.5876	4	1184.8680

TABLE II .

 II Parameters of the primitive basis used for the rovibrational calculations of H2O-H2. FFT stands for the Fast Fourier Transform. Wigner stands for the Wigner DVR. KLeg is the K-Legendre DVR. The units for distance and angle are bohrs and radians respectively.due to the symmetry of H 2 O (and H 2 ). If we write an initial rotational eigenstate of H 2 O as |j, ζ, m (H 2 as |l, m ), then the state-to-state probabilities from the ini-Next we have to use a suitable Complex Absorbing Potential (CAP). The intensity of the CAP was determined using the program plcap included in the MCTDH package. As before the form of the CAP is W

	Coordinate Primitive Number of Range	Size of
		Basis	Points		SPF basis
	R	FFT	960	2.0 -72.0 10-20
	β	Wigner	6	0 -5	10-25
	kγ	K	11	-5 -5
	kα	K	11	-5 -5
	θ	KLeg	3	0 -2	10-12
	k φ	K	5	-2 -2

tial states |j, ζ, m and |j, ζ, -m (|l, m and |l,

-m

) are identical and instead of running for a rotational state j KaKc of H 2 O and l of H 2 (2j + 1) × (2l + 1) calculations, we can simply run (j + 1) × (l + 1) calculations to obtain the cross-sections.

(ii)

TABLE III .

 III FIG. 4. Comparison between MCTDH calculations with 4 different set of basis and Exact (MOLSCAT) time independent close coupling calculations. Each MCTDH calculations is shifted by 0.2 Å2 downward with respect to the previous one and the most accurate calculation is closer to the MOLSCAT result. The parameters of the basis are given in TableIII. Parameters of the SPF basis used in Figure4.

		0.2 0.4				molscat Basis 1 Basis 2 Basis 3 Basis 4
	Cross section (Å 2 )	-0.4 -0.2 0
		-0.6
		-0.8	100	120	140	160	180	200
					E total (cm -1 )
				Basis 1 Basis 2 Basis 3 Basis 4
			nR	10	12	15	20
			n β,kγ ,kα	10	15	20	25
			nθ ,k φ	10	10	12	12

ACKNOWLEDGMENTS This research was supported by the U.S. department of Energy Office of Science, Office of Basic Energy Sciences [Award No. DE-SC0010616] to R.D. S.A.N thanks Hans-Dieter Meyer for various discussions related to the MCTDH implementation of this calculation. Y.S thanks Alexandre Faure for providing us the close-coupling data.