Yannick Jeanvoine
email: yannick.jeanvoine@univ-evry.fr

Riccardo Spezia
email: riccardo.spezia@sorbonne-universite.fr

The formation of urea in space. II. MP2 vs PM6 dynamics in determining bimolecular reaction products

In the present work, we have investigated the possibility of forming protonated urea in the gas phase by means of chemical dynamics simulations. Based on previously published highly correlated quantum chemistry calculations (Astron. Astrophys. 610, A26, 2018), we have considered the reaction between the high energy tautomer of protonated hydroxylamine (NH2OH2 + ) and neutral formamide. Simulations were made at MP2 level and using three semi-empirical Hamiltonians which allow better statistics. In particular we have considered the PM6 method and two different dispersion corrections. These more approximated methods show results which are in relatively good agreement with MP2, in particular for the reaction which is potentially responsible for the urea synthesis. Results show that precursor of protonated urea can be formed but this species will evolve with difficulty into the structure of urea in ultra-vacuum conditions. It is likely that the presence of mantle ice would facilitate the overall reaction process.

Introduction

Over the last years, thanks to the development of astrophysical observations, the presence of many complex organic molecules (COMs) has been reported or suggested in space, in particular in the interstellar medium (ISM) and in comets. [START_REF] Menten | Interstellar Molecules[END_REF][START_REF] Petrie | [END_REF]3,4 Interesting cases of such interstellar and cometary COMs are those containing one (or more) C-N bonds, which can be seen as elementary blocks of a peptide bond, like formamide, 5 cyanamide, 6 acetamide, 7,8 cyanoformaldehyde, 9 methylamine, 10,11 aminoacetonitrile, 12 N-methylformamide, 13 urea 14 and benzonitrile. 15 Due to their nature, and thanks to some experimental and theoretical studies, they are considered as prebiotic molecules which may be at the origin of life. 16,17 From a theoretical point of view, one important question is to point out the possible reactions which may be at the basis of their synthesis. Space conditions are extreme: low density and temperature, limited number of species present, etc. In this context, bi-molecular reactions, both neutral and ionic, may be at the origin of the formation of several elementary species. 18,19 In addition to gas phase reactions, different authors have proposed that they occur in the ice mantles of the interstellar grains, inside or at the surface. 20,21,22,23 In this case, the temperature is still very low, but the density is that of ice: more complex reactions are thus possible, and energy can be dissipated in the solid matrix. The aspect of energy dissipation is crucial. In fact, if we schematize two possible bi-molecular reactive events:

A + B → C (I) A + B → C + D (II)
where C is the desired final product, the reaction (I) will have, in the gas phase, an excess energy which must be dissipated to avoid the unimolecular decomposition of C back into reactants. In the ultra-vacuum conditions this cannot occur via collisional deactivation with an environment but only through the spontaneous emission of a photon in the infra-red region.

The ratio between the rate constant of unimolecular decomposition and that of spontaneous emission will determine if C can be stabilized. However, for organic molecules the spontaneous emission of a photon is often much slower than unimolecular decomposition, such that reactions like (I) are unlikely to be at the origin of a COM in ultra-vacuum conditions. 24 They are possible, on the other hand, in the presence of an environment, like an interstellar grain. Reaction type (II) does not suffer the problem of excess energy in the product C because it can be dissipated as relative translational energy of C + D, thus leading to stable molecules. Bimolecular reactions studied as possible synthetic routes in the gas phase should be, thus, of such a kind. One important condition is that the reaction pathway connecting reactants with products should be barrierless. In fact, at very low temperatures typical of ISM (in molecular clouds, like the Sgr(B2) where urea was detected, 14 the typical temperature is in the 10 -20 K range as reported by Ferriere 25 ), thermal reactions are basically not possible. One crucial aspect is thus to find exothermic reactions with an activation energy which is lower than the energy of reactants.

Generally, bimolecular reactions can occur between two neutrals or a neutral and an ion.

Neutral-neutral reactions have often a high barrier, which can be lowered if one (or both) is a radical. For example, Barone and co-workers have recently shown that formamide, glycolaldehyde, acetic acid and formic acid can be formed in gas phase conditions by such radicalic reactions. 26,27,28 Another way to decrease the activation barrier is in the case of ionneutral reactions: in fact, since organic molecules are generally polar, the dipole-charge interaction is favorable, thus lowering the barrier to make the reactants close to each other.

However, it is important to study in detail the potential energy surface (PES) of each single reaction to determine if it is plausible or not under ISM conditions. This was done, for example, to investigate the possible formation of formamide 29,30 and glycine. 31,32 Recently, we have investigated the PES relative to the synthesis of urea in gas phase, by studying the formation and activation energies of a series of radical, neutral and ion-molecule reactions, finding that the reaction (III) NH2OH2 + + HCONH2 → UreaH + + H2O proceeds via a barrierless pathway. 33 Note that in this reaction, similarly to others which may be at the origin of the formation of protonated formamide and glycine via ion-molecule reactions, 31,32,34,35 one reactant is protonated hydroxylamine. The presence of NH2OH in the ISM is still an open question, since it has not been detected yet. Some years ago, Pulliam et al. 36 provided a detailed study on the non-detection of this neutral species suggesting that possibly the protonated and highly reactive one may be formed. In particular, the high energy tautomer (NH2OH2 + ) could be formed and the tautomerization to the more stable NH3OH + form prevented due to a barrier of about 20 kcal/mol, which are difficult to be passed under ISM conditions.

The bi-molecular reaction is a complex phenomenon, where dynamical aspects can take place. At this aim, we have recently studied the ion-molecule reactions which may be at the origin of the formation of formamide 34 and glycine 35 via chemical dynamics simulations. In this approach, the system is prepared in order to bring the reactants together, for example by giving some collisional energy, and the system is let to evolve following Newton's equations of motion in which energy and gradients are calculated on-the-fly. The statistical treatment of the reaction is important to obtain converged quantities as reaction products probabilities and/or cross sections, which can be related to experiments in some cases. This means that high-level calculations are computationally undoable, and the theory level should be lowered either decreasing the basis set, either using semi-empirical Hamiltonians. In particular, we have shown that for formamide synthesis, MP2 with a relatively small basis set (6-31G(d,p)) could be successfully used, 34 while for the synthesis of glycine the PM6-D semi-empirical Hamiltonian was able to provide a good description. 35 The aim of this work is thus to investigate the feasibility of reaction (III) via chemical dynamics simulations and comparing MP2/6-31G(d,p) results with those using PM6 without dispersion correction and adding dispersion in two ways, by the simple D method 37 or by the more recent D3 38 (both developed by Grimme and co-workers). Reaction dynamics were made (i) by performing explicit collisions with some given relative energy between the ion and the neutral and (ii) by generating clusters of the reactants with random orientations at low temperature. The first model considers the reaction as a linear collision between the ion and the neutral: in space, this can occur if the ion has some linear velocity, which can be generated from acceleration of an ion with an initial (even small) velocity by interacting with a magnetic field and/or from cloud-cloud collisions. 39,40 The cluster model aims to consider all possible interaction orientations at given starting distances, as would occur in a relatively dense and/or confined medium. Furthermore, in this way we can increase the interaction time between reactants with respect to collisional simulations, thus giving the possibility of other reaction pathways to occur. Chemical dynamics results were also compared with kinetic branching ratio calculated using the results of the PES studied previously. 33

Methods 2.1 Quantum chemistry calculations

The geometries of different stationary points (minima and transition states, TSs) were obtained at different levels of theory, and notably: MP2 with 6-31G(d,p) basis set, PM6, 41 PM6-D and PM6-D3, where D stands for the D1 method by Grimme 37 and D3 for the new dispersion correction developed again by Grimme. 38 Note that the MP2 calculations were made with the relatively small basis set which corresponds to the higher theory level doable to perform reaction dynamics of urea formation with reasonable statistics (see section 2.2). Furthermore, this same method was successfully used recently to study reaction dynamics responsible for the formation of formamide 34 and static results presented here will be compared with high level calculations made by us recently. 33 Stationary points were verified by calculating frequencies, being all positive for minima and with only one imaginary frequency for TSs. Gaussian 09 was used for MP2 calculations [START_REF] Frisch | Gaussian09, Revision D.01[END_REF] while MOPAC for semi-empirical Hamiltonians. [START_REF] Stewart | MOPAC, version 5.022mn based on MOPAC 5[END_REF]

Chemical dynamics simulations

Bimolecular reactions were performed using as starting structures of the reactant molecules, NH2OH2 + and formamide (HCONH2), the minimum energy structures obtained by the different methods. Simulations are performed using, as previously, MP2/6-31G(d,p), PM6, PM6-D and PM6-D3 methods to calculate energies and forces on-the-fly.

Collisional simulations are made as follows. Minimum energy structures as obtained with the respective methods were considered as reference structures. The initial coordinates and momenta of reactant molecules were obtained separately from vibrational normal mode sampling. Normal modes quantum numbers were obtained from a Boltzmann sampling at T = 15 K. Rotational energy is then added giving RT/2 to each principal axis of inertia at the same T = 15 K. Initial coordinates and velocities are then transformed in Cartesian coordinates. Then the collisional system is set. First the two reactants were placed at a given distance (here 7 Å) at which they do not interact. Then, a relative kinetic energy between the two collision partners is set in the center-of-mass framework. The two are placed initially on a straight line, with velocity vectors in opposite orientation in order to make them colliding. Then the velocity vector of one partner is shifted by a given value with respect to the center of mass of the other molecule, this defining the so-called impact parameter, b. Details and explanations of initial conditions for bi-molecular collisions can be found elsewhere. 44,[START_REF] Hase | [END_REF] Different values of collision energy were considered, and notably: 5, 10 and 15 kcal/mol. The impact parameter, b, was randomly sampled in the interval [0bmax]. Two sets of trajectories were considered, with bmax = 0.5 Å (to maximize the impact between the reactants) and bmax = 2.0 Å (to check the effect of increasing the centrifugal barrier).

Once the initial conditions are set as described previously, trajectories were propagated by solving Newton's equation of motion with velocity Verlet algorithm 46 with a time step of 0.1 fs for PM6, PM6-D and PM6-D3 and 0.2 fs for MP2, which are a good compromise between accuracy and computing time. Trajectories were stopped when the two final products were at a distance larger than 7.5 Å (resulting trajectories are in the 1.5-2.0 ps range as average) and, in any case, they are run for a maximum simulation time of 5 ps for PM6, PM6-D and PM6-D3 and 2 ps for MP2. In the case of PM6, PM6-D and PM6-D3 we run 6,000 trajectories for each value of collision energy, bmax value and method, while for MP2 we run about 100 trajectories for each set.

In addition, we have performed simulations in which the two reactant molecules were placed at a given distance with random orientation, forming complexes which will have some potential energy which can be used for reactivity. Two values of initial distance were used: 3 and 3.2 Å. These values were set in order to have an initial complex which is not too repulsive (closest distances will result in too big repulsive overlap between the reactants) neither too weakly interacting (which will result in no reactivity). As previously, the reactants were prepared with a vibrational and rotational temperature of 15 K. The same conditions to stop a trajectory were applied and the same number of trajectories were made for the different sets of simulations.

All chemical dynamics simulations were performed by using the VENUS software 47,48 coupled with Gaussian09 for MP2 trajectories and with MOPAC for PM6, PM6-D and PM6-D3 calculations.

Results

Key intermediates and transition states

In a previous work, 33 we have shown that the NH2OH2 + + HCONH2 reaction in the gas phase can produce urea (but also other products) via a barrierless mechanism. The key steps responsible for the first different reaction steps are shown in Figure 1. In particular, the system can proceed once the complex I3A is formed towards the formation of protonated urea (via structures TS3A, I4A, TS4A, I5A and TS5A) or to proton transfer products (via structures TS2A, I2A, TS1A and I1A). More details are in Figure 2 of Ref. 33. One particular aspect is that the TSs are slightly negative at MP2/aug-cc-pVTZ level of theory. We have thus compared the energetics of these key structures in the initial part of the PES as resulting from the methods used in chemical dynamics simulations with the highly correlated calculations made previously. Results are summarized in Table 1. 1.

Table 1. Relative energies of key intermediates and transition states characterizing the initial steps (see

Overall, the different methods provide qualitatively similar results, with differences which we will comment on the following. First, we should note that MP2/6-31G(d,p) calculations provide a description which is very similar to the one with the very large basis set and CCSD(T) results. Energies for the formation of the intermediates (I1A, I2A, I3A, I4A and I5A) are similar or very similar, while differences in the TSs (in particular in the TS3A) are small and preserve the barrierless character of the reaction. Semi-empirical Hamiltonians provide a formation energy of the intermediates which is underestimated with respect to MP2 calculations, as well as for TSs. Only TS5A has different energetic. However, even if the energies of TSs are now less negative, the reaction is still barrierless and the relative stability order of the products is conserved. The addition of dispersion increases the quality of results slightly, but differences between D and D3 are very small. We should note that in some cases (I1A, I4A, TS1A, TS2A and TS5A) the geometries obtained by semi-empirical methods are very different from MP2 ones, such that the comparison is somehow fuzzy. Thus, we have calculated the relative energies also using MP2 geometries: results do not change too much except for the TS5A structure, providing a better agreement.

Overall, the differences between PM6, PM6-D and PM6-D3 are small. In some cases, for example for TS3A, PM6-D provides an energy which is in better agreement with reference calculations than PM6 and PM6-D3. In this complex, the dispersion interaction is clearly important to correctly describe the N-N and N-O interactions, and thus PM6-D provides better results than PM6. PM6-D3 also has dispersion, which is in principle calculated in a more rigorous way, but results are slightly worse than PM6-D. In any case the difference is small (about 1 kcal/mol) and largely inside the uncertainty of the methods.

We will now see how the dynamics made with MP2/6-31G(d,p), PM6, PM6-D and PM6-D3 will differ, eventually.

Reaction products in collisional dynamics

We will now report results obtained by ion molecule collisions at different relative energies. First we will discuss MP2/6-31G(d,p) simulations. By analyzing the final products, we observed that the following reactions occur: Also, a fraction of trajectories does not give any reaction products (i.e. the final structures are the same as the reactants) and in other cases other products (which correspond to further fragmentation of first products) are obtained. The abundance of the main reaction products is reported in Figure 2 while details for all reaction products are given in Table 2. Results as obtained with both bmax = 0.5 and 2.0 Å are reported as a function of collision energy for the different methods used.

NH2OH2 + + HCONH2 → NH3OH + + HCONH2 (a) → NH3COH + + NH2OH (b) → NH2CHOH + + NH2OH (c) → CON2H5 + + H2O (d) → (NH2OH:HCONH2)H + (e)
The most abundant product is the formation of NH2CHOH + and NH2OHreaction (c). Alternatively, protonated formamide can be obtained as N-protonated tautomerreaction (b) which is less abundant (this structure is less stable than the other tautomer of about 12 kcal/mol, see Table 3). Finally, the NH3OH + + HCONH2 reaction product is also observed with less abundance. While energetics of reactions (a) and (c) are very similar (with a slightly preference for (a) at MP2/6-31G(d,p) level of theory), reaction (c) is much simpler to carry out, since it needs only a single proton transfer, while reaction (a) needs an exchange of the protonation site of hydroxylamine from the oxygen to the nitrogen atom.

Reaction (d), which can be responsible to the formation of urea (loss of H2O), is observed even if in low abundance. However, the CON2H5 + ion obtained from simulations, does not have the structure of the intermediates identified in the barrierless minimum energy path connecting the reactants with the products (neither does the one of urea). The structure obtained from MP2 dynamics is shown in Figure 3, where we report also some snapshots (and key distances) extracted from one prototypical trajectory showing the simple reaction event which is responsible for the formation of reaction product (d). Based on our previous study on the PES, this structure cannot evolve to urea without a considerable amount of energy, which would not be likely in cold conditions. Finally, the limited statistic, allowed by the number of trajectories doable in a reasonable computational time at MP2 level of theory, provides a statistical uncertainty corresponding to reaction (d) which is relatively high (on the same order of magnitude of the same reaction product percentage, see Table 2), such that one could question if these observations have a statistical ground. This is why it is important to confirm this result by employing a theoretical method which is much faster thus allowing a much larger number of trajectories with smaller uncertainties on the product abundances. The effect on results of bmax parameter is negligible, while the dependence with the collision energy of the abundance of different reaction pathways is instructive. Reaction (a) probability decreases with collision energy, as typical of barrierless reactions: increasing the collision energy corresponds to decreasing the time the reactants are in the vicinity and the reactive event can occur. In this reaction two proton transfers should occur, and thus lowering the time available for the two events will correspond in decreasing the reaction probability. On the other hand, reactions (b) and (c) correspond to a simple proton transfer from the NH2OH2 + ion to the neutral formamide, forming two tautomers of neutral formamide. Reaction (b) corresponds to the formation of the high energy tautomer and thus its relative abundance increases with collision energy. Finally, reaction (d) is almost independent on the collision energy. As noticed previously, a better evaluation of the relative abundance corresponding to this reaction will be made using trajectories done at semiempirical level which allow to reduce the uncertainties associated to the simulation results. We should also point out that at the end of simulations only a negligible fraction of surviving ion-molecule clustersproduct (e)is obtained, which is an indication that the simulation time is enough to dynamically describe at least the first steps of bi-molecular reactions. The statistical limit using kinetics theory is reported in detail in section 3.4.

From the data reported in same Table 2, we can compare in detail the performances of semiempirical Hamiltonians and how the dispersion affects the results. Concerning the reaction (d), which can be related to urea synthesis, we observe that all the methods are able to provide this reaction with similar ratio. Now, given the larger number of trajectories done, the associated uncertainties are much smaller, providing a statistical ground to the results showing the formation of the ion CON2H5 + . Furthermore, the structure of CON2H5 + ion obtained from semi-empirical simulations is mostly the same as the one of MP2 simulations, and only few trajectories (less than 0.2 %) result in another isomer (the one that is obtained by random cluster initial conditions, see section 3.3). Note that the simulation time-length even in semi-empirical Hamiltonian simulations is surely not enough to observe further isomerization of this structure. However, as remarked previously, the isomer of ion CON2H5 + obtained in collisional simulations cannot evolve in cold conditions to urea: the aim of simulations is to understand if the key steps of a bimolecular reaction can lead (at least) to well identified structures able to evolve into urea (or other interesting species). Reaction (c) is also the most abundant one, with percentages similar to MP2 results. On the other hand, reaction (b) is underestimated while reaction (a) is overestimated. The main difference is that when using semi-empirical Hamiltonians, the number of trajectories forming more than two products (which here are classified as "other reactions") is higher. Thus, the ion NH3COH + which corresponds to reaction (b) can further fragment.

It is interesting to note that even if reaction (a) is more exothermic at MP2 level than with semi-empirical Hamiltonians (see Table 3), this is not reflected by a lower abundance as simulation products, and conversely for reaction (b). Finally, the behavior of reaction products abundances as a function of collision energy is very similar to what obtained by MP2 trajectories (in Figure 1 we report the case of PM6-D simulations).

Reaction products of random complexes

Collisional simulations model the fact that the two reactants come close each other with a given translational energy. This can result in short interaction time with the consequence that some reaction pathways are rarely followed. In space, the conditions which can bring two species together can be of various origin and the density is not uniform, such that some dense conditions can be at the origin of formation of clusters with different orientations. Thus, to have a more exhaustive picture of the bimolecular reaction dynamics, we also consider initial conditions in which the reactants are randomly placed close together, forming what we call "random complexes". Two initial distances were considered (3.0 and 3.2 Å) in these simulations, and we obtained reactions similar to collisional simulations ones, but with some interesting differences. The percentages of the different reaction products are reported in Table 4.

+ + H2O (II) 3 ± 2 (4 ± 2) 3.7 ± 0.2 (3.6 ± 0.2) 3.6 ± 0.2 (3.4 ± 0.2) 3.8 ± 0.2 (3.6 ± 0.2) Product (e) (NH2OH:HCONH2)H + 1 ± 1 (1 ± 1)
0.5 ± 0.1 (0.7 ± 0.1) 0.6 ± 0.1 (0.9 ± 0.1) 0.5 ± 0.1 (0.7 ± 0.1) Other reactions (III) 10 ± 3 (5 ± 2) 20.0 ± 0.5 (13.9 ± 0.4) 21.2 ± 0.5 (15.4 ± 0.4) 20.0 ± 0.5 (14.1 ± 0.4) (I) Formally they are the same products as when any reaction occurs, but here we have a proton exchange between the reactants. (II) The ion CON2H5 + was observed in two different forms, see text for details. (III) "Other reactions" correspond to formation of three or more products.

As in collisional simulations we obtain reaction products (a)-(e). Other reaction products (corresponding to the formation of three or more products) were found with higher abundances and in about 5% of trajectories the same reactants are formally obtained, but with proton exchanges. This last reaction can occur since in cluster simulations the reactants have, by construction, enough time to stay close to each other and thus two proton exchanges are possible, while in collisions, because of the scattering event, this is not observed. The increasing abundance of reaction products composed by more than two molecules can be explained by the fact that some random initial orientations have high internal energy making first product molecules unstable which can thus break more easily. Reactions (a)-(e) are the same as in collisional simulations, but with some differences. Reactions (a), (b) and (e) have almost the same abundance, when using MP2, while in the case of semi-empirical Hamiltonians reaction (b) is now more probable. Reaction (c) is still the most abundant one but the percentages are now much smaller than in collisional simulations. This is due to the fact that there are additional pathways as discussed previously. Finally, reaction (d) is slightly more abundant. It is interesting to note that in this case the dominant structure of the ion CON2H5 + is different from the one obtained in collisional simulations; one geometry is shown in Figure 4A. Interestingly this corresponds to the structure I4A on the PES which is the first key intermediate product where H2O is formed and that can evolve into protonated urea. In Figure 4B we show a prototypical trajectory showing how this intermediate can be formed together with the evolution of two distances able to describe the progress of the reaction (Figure 4C). However, in this structure the water molecule is no longer interacting with the ion. Thus, the system could not evolve on the exact PES which leads to urea with a barrierless pathway but on a slightly different one. This PES was also studied in our previous work (see Figure 5 of Ref. 33) and, as we have remarked, the absence of water molecule increases the energy of transition states, such that the reaction is no more barrierless. In other terms, due to the absence of any medium which helps the water molecule to stay close to the ion, once the intermediate is formed it cannot isomerize to the urea structure.

Concerning the comparison between MP2 and semi-empirical Hamiltonians, also in this case the differences are relatively small. As for collisional simulations, the semi-empirical Hamiltonians have the tendency of overestimating the products composed by more than two species. This is at the origin of the decreasing of the relative abundance of reaction (c). It is important to notice that, as for collisions, all the semi-empirical Hamiltonians show that products of reaction (d) can be formed and that now it is possible to form also a structure similar to I4A as in MP2 dynamics.

Reaction Kinetics

Chemical dynamics have shown that reaction (III) can be at the origin of the formation of precursors of protonated urea, even if with relatively low probability. It is now interesting to compare the dynamical results with the kinetic description. At this end, we used the information in terms of minima and transition states on the PES for reaction (III) at CCSD(T)/aug-cc-pVTZ//MP2/aug-cc-pVTZ level of theory, schematically shown in Figure 1 (more details are in Figure 2 of Ref. 33). Based on molecular properties it is possible to determine, under the statistical limit, the reaction rate constants and branching ratios.

The first step is the ion-molecule capture, whose rate constant can be evaluated by the usual Su-Chesnavich approach. 49 This is based on the well-known temperature independent Langevin rate coefficient: (1) where q is the charge,  the scalar polarizability of the neutral and M the reduced mass of the reactants. In the case of NH2OH2 + + HCONH2 reaction, we have kL = 2.89 • 10 -9 cm 3 s -1 . To recover the temperature dependence, the parameter x is defined as (2) where D is the dipole moment, T the temperature and kB the Boltzmann constant. For reaction (III) at 15 K, x = 8.6, while at 300 K (as typically reported in the KIDA database 50 ), x = 1.7. The temperature dependence of the rate coefficient, kD, is then:

(3) (4)
In NH2OH2 + + HCONH2 bi-molecular reaction the I3A complex is formed first, and then the reaction can proceed towards two pathways: (i) I4A, I5A etc. … eventually up to urea. As reported in details in Ref. 33, this pathway can also branch towards the formation of HNCO + H2O + NH4 + ; (ii) I2A and then the proton transfer products NH3OH + + HCONH2, NH2OH + HCONH3 + and NH2OH + HCOHNH2 + , corresponding to reaction pathways (a), (b) and (c), respectively. For simplicity, we will refer to "urea pathway" for pathway (i) and "proton transfer pathway" for pathway (ii).

Based on energies and frequencies of minima and transition states, it is possible to calculate the micro-canonical rate constant using the well-known Rice-Ramsperger-Kassel-Marcus (RRKM) theory: 51 (5) where  is the reaction path degeneracy, h the Planck constant, N † (E-E0) the sum of states at the transition state, (E) the density of states at the minimum, E the internal energy and E0 the barrier height. The quantum sum and density of states were calculated using the direct count algorithm [START_REF] Beyer | [END_REF] with the Zhu and Hase code. 53 Two rate constants correspond to the first steps connecting the I3A with the two reaction pathways: k34, corresponding to passing the first transition state, TS3A, for the "urea pathway" and k32, corresponding to passing the transition state TS2A, for the "proton transfer pathway". Rate constant values are reported in Table 5 for different energies. Note that the I3A complex is stabilized by about 50 kcal/mol with respect to reactants, such that once formed in a bi-molecular process it will have at least this amount of energy. This corresponds to E = 0 in Table 5. As we can notice, k32 >> k34, such that we can consider that under statistical conditions the pathway leading to urea will not be followed. On the other hand, chemical dynamics simulations, which do not presuppose any statistical redistribution of the energy in the formed complex, show a substantially different picture. This is an indication that non-statistical effects can be important in the appearance of less statistically favored products. In this same Table 5 we also report the rate constants for the subsequent reactions. We should note that the k34 rate constant is two or three orders of magnitude smaller than the ones of the subsequent processes, k45 and k56, while it is larger than the back-reaction, k43. Note also that the back reactions are generally smaller than the corresponding forward ones, such that one can consider that once this pathway is followed, it can proceed towards the products.

From the Laplace transform of rate constants it is possible to estimate the branching ratio between urea and proton transfer pathways as a function of temperature, UP(T):

where k34(E) and k32(E) are the RRKM micro-canonical rate constants and I3A(E) is the density of states of the I3A intermediate. In the same way, it is possible to calculate the canonical rate constant from the micro-canonical one: (7) and in the present case k34(T=15K) = 6.10 • 10 +3 s -1 and k32(T=15K) = 1.30 • 10 +12 s -1 . Note that these rate constants are calculated with quantum density of states on the CCSD(T)/augcc-PVTZ//MP2/aug-cc-PVTZ PES. Direct dynamics simulations are, however, purely Newtonian and done on the MP2/6-31G(d,p) on-the-fly PES. Using classical density of states on the MP2/6-31G(d,p) PES, the classical rate constants result to be k34(T=15K) = 1.28 • 10 -10 s -1 and k32(T=15K) = 1.78 • 10 +12 s -1 . The classical value for k34 is substantially different from the quantum one because the classical barrier between I3A and TS3A is higher (46 vs 41 kcal/mol) and for low energies, the sum of states at the TS is considerably smaller. It is well known that classical sum and density of states are a valid approximation only at high energies. [START_REF] Baer | Unimolecular Reaction Dynamics : Theory and Experiments[END_REF] From the information on the different rate constants reported previously, it is possible to evaluate the branching ratio of the different products and the temperature dependent rate constants, kD(T). Results for low temperature (15K) and ambient temperature (300K) are reported in Table 6. Ambient temperature results are shown since normally in the KIDA database the 300K results are reported (other temperatures can be extrapolated further) and they can be eventually used in future room temperature experiments. Using the classical density of states, corresponding to the direct dynamics performed here, the formation of urea is even less probable. Overall, statistical results show that at 15 K only the NH2OH + HCOHNH2 + products should be formed. This is not what was observed in direct dynamics, suggesting that non-statistical effects can play an important role, as observed already in different contexts. 54 The reason for this difference can reside in the fact that the statistical limit approach assumes that the intramolecular vibrational-energy redistribution (IVR) is much faster than subsequent reactions. While more quantitative studies are surely needed, the comparison between direct dynamics and statistical modeling suggests that the reaction time scale is faster than IVR in the present case.

Conclusions

In this work, we have studied the reaction dynamics between protonated hydroxylamine (NH2OH2 + ) and formamide (HCONH2), which was shown in a previous static study to be potentially responsible to the formation of urea in gas phase. 33 In particular, this was motivated by the observation of urea in the ISM 14 and reactants and initial conditions were chosen in order to mimic reasonable space conditions.

We have run a set of trajectories at MP2/6-31G(d,p) level of theory and results compared with simulations done with PM6, PM6-D and PM6-D3 semi-empirical Hamiltonians, being computationally much less expensive. Simulations show that these methods can provide results in agreement with MP2 ones (for some reaction products the agreement is more quantitative than for others). In particular, we did not notice relevant differences between PM6, PM6-D and PM6-D3. We suggest in future studies (on this or related systems) to include dispersion. Between D and D3 method we noticed that D3 is computationally slower than D (due to the fact that in D3 also three body interactions are calculated), such that if there are not particular reasons to use D3, the simple D correction can be sufficient and the resulting computational cost is basically the same as simple PM6. It is important to notice that the use of semi-empirical Hamiltonians allows the running of many more trajectories (6000 vs 100 for each set) and thus significantly decreases the uncertainty associated to reaction products.

Concerning the dynamical process which would possibly lead to protonated urea under ISM conditions, we have obtained that in collisional simulations (also with relatively small collisional energy) an ion with the mass of protonated urea is formed, but with a structure which will difficultly evolve to the desired target. On the other hand, the simulations made with random clusters show that it is possible to form a product which was identified as a first step toward the protonated urea synthesis. However, two problems arise: (i) the random cluster simulations show that to form the intermediate which can lead to urea, the reactants should be forced to have the two N atoms relatively close to each other at the beginning, which would be unlikely if the gas phase is not relatively dense; (ii) the intermediate formed (I4A) is now de-hydrated, and it was shown that a barrierless process needs a water molecule to stay in the vicinity of the I4A species. These two conditions ("compression" of the reactants and the presence of a water molecule in the vicinity of the first reaction product) suggest that the reaction leading to urea can be favored in a dense environment, like ice and/or the vacuum/ice intermediates, similarly to what is proposed for other species. 18,18,19,21 By comparing simulation results with kinetics calculations, we observed that some relevant pathways are accessible only in a dynamical regime, showing that non-statistical processes are important.

Concluding, it would be interesting in the future to investigate, by chemical dynamics, how this reaction will proceed on ice mantle models. At this purpose, the good performances of semi-empirical Hamiltonians are a promising result, since simulations of the reaction with a surface will be computationally heavy such that methods like MP2 cannot be used, but PM6-D, for example, can. An important aspect, missing in the present dynamical model, concerns quantum nuclear effects, like tunneling and zero-point energy conservations. To this aim recent improvements, like centroid or ring polymer molecular dynamics, 55,56 could be considered for future studies. Given that these methods often require increasing computational effort, the possibility of using relatively fast methods, like PM6-D, to calculate energies and forces, paves the way for using them for these kinds of reactions.

Figure 1 .

 1 Figure 1. Schematic potential energy surface corresponding to the initial steps of NH2OH2 + + HCONH2 reaction. Energy levels are schematic, values as obtained with the different levels of theory are reported inTable 1.

Figure 2 .

 2 Figure 2. Evolution of reaction products % as a function of the collision energy (Erel) as obtained from MP2/6-31G(d,p) (solid lines) and PM6-D (dashed lines) simulations. Results for the two bmax values employed are shown.

Figure 3 .

 3 Figure 3. Panel A: structure of the ion CON2H5 + formed as product of reaction (d) I in collisional simulations. Panel B: snapshots taken from a prototypical trajectory providing reaction (d) I . Panel C: atomic distances showing how the reaction proceeds in time.

Figure 4 .

 4 Figure 4. Panel A: structure of the ion CON2H5 + formed as product of reaction (d) II in simulations in which random complexes at 3.2 A are set as initial conditions. Panel B: snapshots taken from a prototypical trajectory providing reaction (d) II in this same kind of simulations. Panel C: atomic distances showing how the reaction proceeds in time.

Table 2 .

 2 Percentage of different reaction products obtained from collisional simulations at different bmax values as a function of collision energy (Erel) as obtained with MP2, PM6, PM6-

	D and PM6-D3 dynamics.
	bmax (Å) Erel (kcal/mol) Reaction	MP2	PM6	PM6-D	PM6-D3
	0.5	5	(a)	8 ± 3	14.4 ± 0.5 17.7 ± 0.5	14.0 ± 0.5
			(b)	12 ± 3	1.3 ± 0.1	1.9 ± 0.2	1.4 ± 0.1
			(c)	76 ± 4	79.0 ± 0.5 74.0 ± 0.5	79.2 ± 0.5
			(d)	1 ± 1	1.1 ± 0.1	1.3 ± 0.1	0.8 ± 0.1
			(e)	-	0.5 ± 0.1	0.7 ± 0.1	0.5 ± 0.5
			No reaction 3 ± 2	1.6 ± 0.1	1.8 ± 0.1	1.9 ± 0.1
			other	-	2.1 ± 0.2	2.6 ± 0.2	2.2 ± 0.2
		10	(a)	5 ± 2	10.5 ± 0.4 12.7 ± 0.4	11.1 ± 0.4
			(b)	10 ± 3	1.7 ± 0.1	2.4 ± 0.2	1.6 ± 0.1
			(c)	76 ± 4	77.6 ± 0.5 74.3 ± 0.5	78.3 ± 0.5
			(d)	1 ± 1	2.7 ± 0.2	2.8 ± 0.2	2.6 ± 0.2
			(e)	1 ± 1	-	0.4 ± 0.1	0.4 ± 0.1
			No reaction 5 ± 2	5.6 ± 0.3	4.8 ± 0.3	4.5 ± 0.3
			Other	2 ± 1	1.9 ± 0.1	2.6 ± 0.2	1.5 ± 0.1
		15	(a)	-	7.8 ± 0.3	9.6 ± 0.4	9.0 ± 0.4
			(b)	17 ± 3	3.4 ± 0.2	4.6 ± 0.3	3.5 ± 0.2
			(c)	72 ± 3	72.6 ± 0.5 69.3 ± 0.5	71.7 ± 0.5
			(d)	1 ± 1	3.9 ± 0.3	4.4 ± 0.3	3.8 ± 0.3
			(e)	-	-	0.3 ± 0.1	-
			No reaction 9 ± 3	10.3 ± 0.4 8.8 ± 0.4	9.8 ± 0.4
			Other	1 ± 1	2.0 ± 0.2	3.0 ± 0.2	2.2 ± 0.2
	2.0	5	(a)	6 ± 2	13.7 ± 0.4 16.5 ± 0.4	14.1 ± 0.4
			(b)	12 ± 3	1.3 ± 0.1	2.1 ± 0.2	1.2 ± 0.1
			(c)	76 ± 3	80.1 ± 0.5 75.0 ± 0.5	79.8 ± 0.5
			(d)	2 ± 1	0.9 ± 0.1	1.2 ± 0.1	0.7 ± 0.1
			(e)	1 ± 1	0.5 ± 0.1	0.7 ± 0.1	0.4 ± 0.1
			No reaction 3 ± 2	2.1 ± 0.2	1.8 ± 0.2	2.0 ± 0.2
			Other	-	1.4 ± 0.1	2.7 ± 0.2	1.8 ± 0.2
		10	(a)	7 ± 2	9.8 ± 0.4	11.8 ± 0.4	10.5 ± 0.4
			(b)	15 ± 3	1.6 ± 0.2	2.4 ± 0.2	1.9 ± 0.2
			(c)	67 ± 5	79.8 ± 0.5 76.3 ± 0.5	79.6 ± 0.5
			(d)	2 ± 1	2.3 ± 0.2	2.6 ± 0.2	2.1 ± 0.2
			(e)	1 ± 1	-	0.3 ± 0.1	-
			No reaction 6 ± 2	4.8 ± 0.3	4.3 ± 0.3	4.5 ± 0.3
			Other	2 ± 1	1.7 ± 0.2	2.3 ± 0.2	1.4 ± 0.1
		15	(a)	6 ± 2	7.5 ± 0.3	8.5 ± 0.3	7.8 ± 0.3
			(b)	16 ± 3	3.4 ± 0.2	4.3 ± 0.3	3.4 ± 0.2
			(c)	66 ± 5	73.9 ± 0.5 71.8 ± 0.5	74.4 ± 0.5
			(d)	2 ± 1	3.7 ± 0.2	4.2 ± 0.3	3.5 ± 0.2
			(e)	-	-	-	-
			No reaction 8.74 ± 3 9.7 ± 0.4	8.6 ± 0.4	9.2 ± 0.4
			Other	1 ± 1	1.8 ± 0.2	2.6 ± 0.2	1.7 ± 0.2

Table 3 .

 3 Reaction energies considering as reactants NH2OH2 + + HCONH2. Values are in kcal/mol.

	Reaction	MP2/6-31G(d,p)	PM6	PM6-D	PM6-D3
	(a) NH3OH + + HCONH2	-27.36	-23.18	-23.09	-23.17
	(b) NH3COH + + NH2OH	-15.32	-22.16	-22.06	-22.15
	(c) NH2CHOH + + NH2OH	-27.04	-23.16	-23.11	-23.17
	(d) I CON2H5 + + H2O	-37.40	-30.02	-29.97	-30.23
	(d) II CON2H5 + + H2O	-33.16	-32.59	-33.22	-32.69
	(d) III CON2H5 + + H2O	-95.66	-70.80	-71.00	-70.88
	I the geometry of CON2H5 + is the optimized geometry of the isomer resulting from MP2 collisional simulations
	(see Figure 3.A).
	II the geometry of CON2H5 + is the optimized geometry of the isomer resulting from some MP2 cluster reactivity
	and which coincides with geometry I4A without the water molecule (see Figure 4.A).
	III the geometry of CON2H5 + is the optimized geometry of protonated urea.

Table 4 .

 4 Percentages

	Reaction	MP2	PM6	PM6-D	PM6-D3
	No reaction	6 ± 2	5.3 ± 0.3	4.1 ± 0.3	5.2 ± 0.3
		(4 ± 2)	(5.4 ± 0.3)	(4.0 ± 0.3)	(5.2 ± 0.3)
	NH2OH2 + + HCONH2 (I)	15 ± 4	12.2 ± 0.4	11.2 ± 0.4	12.2 ± 0.4
		(12 ± 3)	(8.7 ± 0.4)	(8.3 ± 0.4)	(8.8 ± 0.4)
	Product (a)	8 ± 3	13.1 ± 0.4	15.7 ± 0.4	13.0 ± 0.4
	NH3OH + + HCONH2	(6 ± 2)	(16.6 ± 0.4)	(20.6 ± 0.5)	(17.4 ± 0.4)
	Product (b)	12 ± 3	8.7 ± 0.4	8.5 ± 0.4	8.7 ± 0.4
	NH3COH + + NH2OH	(10 ± 3)	(10.5 ± 0.4)	(9.6 ± 0.4)	(10.3 ± 0.4)
	Product (c)	45 ± 5	36.5 ± 0.6	35.1 ± 0.6	36.6 ± 0.6
	HCOHNH2 + + NH2OH	(58 ± 5)	(40.6 ± 0.6)	(37.8 ± 0.6)	(39.8 ± 0.6)
	Product (d)
	CON2H5

of different reaction products as obtained from simulations starting from random complexes, as obtained with initial distance of 3.0 and 3.2 Å (these lasts are reported in parenthesis).MP2, PM6, PM6-D and PM6-D3 simulation results are compared.

Table 5 .

 5 Rate constants for the elementary uni-molecular reaction events once the I3A complex is formed. Values are reported for four different internal energy values. k(E) are in s -1 and E in kcal/mol.

		E = 0	E = 5	E = 10	E = 15
	k32	1.29 • 10 +12	1.33 • 10 +12	1.37 • 10 +12	1.40 • 10 +12
	k34	6.48 • 10 +3	8.40 • 10 +4	5.32 • 10 +5	2.22 • 10 +6
	k43	3.76 • 10 +2	5.13 • 10 +3	3.38 • 10 +4	1.46 • 10 +5
	k45	6.78 • 10 +5	3.27 • 10 +6	1.13 • 10 +7	3.13 • 10 +7
	k54	3.43 • 10 +5	2.03 • 10 +6	8.20 • 10 +6	2.61 • 10 +7
	k56	2.01 • 10 +6	1.19 • 10 +7	4.96 • 10 +7	1.62 • 10 +8

Table 6 .

 6 Branching

	ratio (in %) and temperature dependent rate constant, kD(T) in cm 3 s -1 for
	the different reaction products of NH2OH2 + + HCONH2 reaction (see section 3.2 for the
	definition of the different reactions).
		Statistical branching ratio (in %) kD (cm 3 s -1 )
	Product	T = 15K	T = 300K	T = 15K	T = 300K
	(a)	9 • 10 -17	11	1 • 10 -26	4 • 10 -10
	(b)	4 • 10 -174	1 • 10 -7	5 • 10 -184	5 • 10 -18
	(c)	100	89	1 • 10 -8	3 • 10 -9
	(d)	5 • 10 -7	8 • 10 -7	6 • 10 -17	3 • 10 -17
	HNCO + H2O + NH4 +	4 • 10 -303	1 • 10 -21	-	4 • 10 -32

Acknowledgments

We thank ANR DynBioReact (Grant No. ANR-14-CE06-0029-01) and CNRS program INFINITI (project ASTROCOL) for financial support. This article is dedicated to the memory of Dr. János G. ÁNGYÁN, who directed the Ph. D. Thesis of one of us (YJ). We hope we will pass on the values and commitment we learned from him to those we work with and hopefully we will be able to do it with János's humor and joy.


Corentin Tallec
email: corentin.tallec@u-psud.fr

Yann Ollivier

UNBIASED ONLINE RECURRENT OPTIMIZATION

) that bypasses the need for model sparsity and makes implementation easy in current deep learning frameworks, even for complex models. Like NoBackTrack, UORO provides unbiased gradient estimates; unbiasedness is the core hypothesis in stochastic gradient descent theory, without which convergence to a local optimum is not guaranteed. On the contrary, truncated BPTT does not provide this property, leading to possible divergence. On synthetic tasks where truncated BPTT is shown to diverge, UORO converges. For instance, when a parameter has a positive short-term but negative long-term influence, truncated BPTT diverges unless the truncation span is very significantly longer than the intrinsic temporal range of the interactions, while UORO performs well thanks to the unbiasedness of its gradients.

)), the current standard for training recurrent architectures, is well suited to many short training sequences. Treating long sequences with BPTT requires either storing all past inputs in memory and waiting for a long time between each learning step, or arbitrarily splitting the input sequence into smaller sequences, and applying BPTT to each of those short sequences, at the cost of losing long term dependencies.

This paper introduces Unbiased Online Recurrent Optimization (UORO), an online and memoryless learning algorithm for recurrent architectures: UORO processes and learns from data samples sequentially, one sample at a time. Contrary to BPTT, UORO does not maintain a history of previous inputs and activations. Moreover, UORO is scalable: processing data samples with UORO comes at a similar computational and memory cost as just running the recurrent model on those data.

Like most neural network training algorithms, UORO relies on stochastic gradient optimization. The theory of stochastic gradient crucially relies on the unbiasedness of gradient estimates to provide convergence to a local optimum. To this end, in the footsteps of NoBackTrack (NBT) [START_REF] Ollivier | Training recurrent networks online without backtracking[END_REF], UORO provides provably unbiased gradient estimates, in a scalable, streaming fashion.

Unlike NBT, though, UORO can be easily implemented in a black-box fashion on top of an existing recurrent model in current machine learning software, without delving into the structure and code of the model.

The framework for recurrent optimization and UORO is introduced in Section 2. The final algorithm is reasonably simple (Alg. 1), but its derivation (Section 3) is more complex. In Section 6, UORO is shown to provide convergence on a set of synthetic experiments where truncated BPTT fails to display reliable convergence. An implementation of UORO is provided as supplementary material.

A widespread approach to online learning of recurrent neural networks is Truncated Backpropagation Through Time (truncated BPTT) [START_REF] Jaeger | Tutorial on training recurrent neural networks[END_REF], which mimics Backpropagation Through Time, but zeroes gradient flows after a fixed number of timesteps. This truncation makes gradient estimates biased; consequently, truncated BPTT does not provide any convergence guarantee. Learning is biased towards short-time dependencies. 1 . Storage of some past inputs and states is required.

Online, exact gradient computation methods have long been known (Real Time Recurrent Learning (RTRL) [START_REF] Williams | A learning algorithm for continually running fully recurrent neural networks[END_REF]; [START_REF] Barak | Gradient calculations for dynamic recurrent neural networks: A survey[END_REF]), but their computational cost discards them for reasonably-sized networks.

NoBackTrack (NBT) [START_REF] Ollivier | Training recurrent networks online without backtracking[END_REF] also provides unbiased gradient estimates for recurrent neural networks. However, contrary to UORO, NBT cannot be applied in a blackbox fashion, making it extremely tedious to implement for complex architectures.

Other previous attempts to introduce generic online learning algorithms with a reasonable computational cost all result in biased gradient estimates. Echo State Networks (ESNs) [START_REF] Jaeger | Tutorial on training recurrent neural networks[END_REF]; [START_REF] Jaeger | Optimization and Applications of Echo State Networks with Leaky-Integrator Neurons[END_REF] simply set to 0 the gradients of recurrent parameters. Others, e.g., [START_REF] Maass | Real-time computing without stable states: A new framework for neural computation based on perturbations[END_REF][START_REF] Jochen | Backpropagation-decorrelation: online recurrent learning with O(N) complexity[END_REF], introduce approaches resembling ESNs, but keep a partial estimate of the recurrent gradients. The original Long Short Term Memory algorithm [START_REF] Hochreiter | Long short-term memory[END_REF] (LSTM now refers to a particular architecture) cuts gradient flows going out of gating units to make gradient computation tractable. Decoupled Neural Interfaces [START_REF] Jaderberg | Decoupled neural interfaces using synthetic gradients[END_REF] bootstrap truncated gradient estimates using synthetic gradients generated by feedforward neural networks. The algorithm in [START_REF] Movellan | A Monte Carlo EM approach for partially observable diffusion processes: Theory and applications to neural networks[END_REF] provides zeroth-order estimates of recurrent gradients via diffusion networks; it could arguably be turned online by running randomized alternative trajectories. Generally these approaches lack a strong theoretical backing, except arguably ESNs.

BACKGROUND

UORO is a learning algorithm for recurrent computational graphs. Formally, the aim is to optimize θ, a parameter controlling the evolution of a dynamical system

s t+1 = F state (x t+1 , s t , θ) (1) o t+1 = F out (x t+1 , s t , θ) (2)
in order to minimize a total loss L := ∂θ by unfolding the network through time, and backpropagating through the unfolded network, each timestep corresponding to a layer. BPTT thus requires maintaining the full unfolded network, or, equivalently, the history of past inputs and activations. 2 Truncated BPTT only unfolds the network for a fixed number of timesteps, reducing computational cost in online settings [START_REF] Jaeger | Tutorial on training recurrent neural networks[END_REF]. This comes at the cost of biased gradients, and can prevent convergence of the gradient descent even for large truncations, as clearly exemplified in Fig. 2a.

UNBIASED ONLINE RECURRENT OPTIMIZATION

Unbiased Online Recurrent Optimization is built on top of a forward computation of the gradients, rather than backpropagation. Forward gradient computation for neural networks (RTRL) is described in [START_REF] Williams | A learning algorithm for continually running fully recurrent neural networks[END_REF] and we review it in Section 3.1. The derivation of UORO follows in Section 3.2. Implementation details are given in Section 3.3. UORO's derivation is strongly connected to [START_REF] Ollivier | Training recurrent networks online without backtracking[END_REF] but differs in one critical aspect: the sparsity hypothesis made in the latter is relieved, resulting in reduced implementation complexity without any model restriction. The proof of UORO's convergence to a local optimum can be found in [START_REF] Massé | Autour de l'Usage des gradients en apprentissage statistique[END_REF].

FORWARD COMPUTATION OF THE GRADIENT

Forward computation of the gradient for a recurrent model (RTRL) is directly obtained by applying the chain rule to both the loss function and the state equation (1), as follows.

Direct differentiation and application of the chain rule to t+1 yields

∂ t+1 ∂θ = ∂ t+1 ∂o (o t+1 , o * t+1 ) • ∂F out ∂s (x t+1 , s t , θ) ∂s t ∂θ + ∂F out ∂θ (x t+1 , s t , θ) . (4)
Here, the term ∂s t /∂θ represents the effect on the state at time t of a change of parameter during the whole past trajectory. This term can be computed inductively from time t to t + 1. Intuitively, looking at the update equation ( 1), there are two contributions to ∂s t+1 /∂θ:

• The direct effect of a change of θ on the computation of s t+1 , given s t .

• The past effect of θ on s t via the whole past trajectory.

With this in mind, differentiating (1) with respect to θ yields

∂s t+1 ∂θ = ∂F state ∂θ (x t+1 , s t , θ) + ∂F state ∂s (x t+1 , s t , θ) ∂s t ∂θ . (5)
This gives a way to compute the derivative of the instantaneous loss without storing past history: at each time step, update ∂s t /∂θ from ∂s t-1 /∂θ, then use this quantity to directly compute ∂ t+1 /∂θ. This is how RTRL [START_REF] Williams | A learning algorithm for continually running fully recurrent neural networks[END_REF] proceeds.

A huge disadvantage of RTRL is that ∂s t /∂θ is of size dim(state) × dim(params). For instance, with a fully connected standard recurrent network with n units, ∂s t /∂θ scales as n 3 . This makes RTRL impractical for reasonably sized networks.

UORO modifies RTRL by only maintaining a scalable, rank-one, provably unbiased approximation of ∂s t /∂θ, to reduce the memory and computational cost. This approximation takes the form st ⊗ θt , where st is a column vector of the same dimension as s t , θt is a row vector of the same dimension as θ , and ⊗ denotes the outer product. The resulting quantity is thus a matrix of the same size as ∂s t /∂θ. The memory cost of storing st and θt scales as dim(state) + dim(params). Thus UORO is as memory costly as simply running the network itself (which indeed requires to store the current state and parameters). The following section details how st and θt are built to provide unbiasedness.

RANK-ONE TRICK: FROM RTRL TO UORO

Given an unbiased estimation of ∂s t /∂θ, namely, a stochastic matrix Gt such that E Gt = ∂s t /∂θ, unbiased estimates of ∂ t+1 /∂θ and ∂s t+1 /∂θ can be derived by plugging Gt in ( 4) and ( 5). Unbiasedness is preserved thanks to linearity of the mean, because both (4) and ( 5) are affine in ∂s t /∂θ.

Thus, assuming the existence of a rank-one unbiased approximation Gt = st ⊗ θt at time t, we can plug it in (5) to obtain an unbiased approximation Ĝt+1 at time t

+ 1 Ĝt+1 = ∂F state ∂θ (x t+1 , s t , θ) + ∂F state ∂s (x t+1 , s t , θ) st ⊗ θt . (6)
However, in general this is no longer rank-one.

To transform Ĝt+1 into Gt+1 , a rank-one unbiased approximation, the following rank-one trick, introduced in Ollivier et al. ( 2015) is used: Proposition 1. Let A be a real matrix that decomposes as

A = k i=1 v i ⊗ w i . ( 7
)
Let ν be a vector of k independent random signs, and ρ a vector of k positive numbers. Consider the rank-one matrix

à := k i=1 ρ i ν i v i ⊗ k i=1 ν i w i ρ i (8)
Then à is an unbiased rank-one approximation of A:

E ν Ã = A.
The rank-one trick can be applied for any ρ. The choice of ρ influences the variance of the approximation; choosing

ρ i = w i / v i (9)
minimizes the variance of the approximation, E A -Ã 2 2 [START_REF] Ollivier | Training recurrent networks online without backtracking[END_REF].

The UORO update is obtained by applying the rank-one trick twice to (6). First, ∂Fstate ∂θ (x t+1 , s t , θ) is reduced to a rank one matrix, without variance minimization. 3 Namely, let ν be a vector of independant random signs; then,

∂F state ∂θ (x t+1 , s t , θ) = E ν ν ⊗ ν ∂F state ∂θ (x t+1 , s t , θ) . ( 10
)
This results in a rank-two, unbiased estimate of ∂s t+1 /∂θ by substituting (10) into ( 6)

∂F state ∂s (x t+1 , s t , θ) st ⊗ θt + ν ⊗ ν ∂F state ∂θ (x t+1 , s t , θ) . (11)
Applying Prop. 1 again to this rank-two estimate, with variance minimization, yields UORO's estimate Gt+1

Gt+1 = ρ 0 ∂F state ∂s (x t+1 , s t , θ) st + ρ 1 ν ⊗ θt ρ 0 + ν ρ 1 ∂F state ∂θ (x t+1 , s t , θ) (12)
which satisfies that E ν Gt+1 is equal to (6). (By elementary algebra, some random signs that should appear in (12) cancel out.) Here

ρ 0 = θt ∂Fstate ∂s (x t+1 , s t , θ) st , ρ 1 = ν ∂Fstate ∂θ (x t+1 , s t , θ) ν (13)
minimizes variance of the second reduction.

The unbiased estimation ( 12) is rank-one and can be rewritten as Gt+1 = st+1 ⊗ θt+1 with the update

st+1 ← ρ 0 ∂F state ∂s (x t+1 , s t , θ) st + ρ 1 ν (14) θt+1 ← θt ρ 0 + ν ρ 1 ∂F state ∂θ (x t+1 , s t , θ). (15)
Initially, ∂s 0 /∂θ = 0, thus s0 = 0, θ0 = 0 yield an unbiased estimate at time 0. Using this initial estimate and the update rules ( 14)-( 15), an estimate of ∂s t /∂θ is obtained at all subsequent times, allowing for online estimation of ∂ t /∂θ. Thanks to the construction above, by induction all these estimates are unbiased. 4We are left to demonstrate that these update rules are scalably implementable.

IMPLEMENTATION

Implementing UORO requires maintaining the rank-one approximation and the corresponding gradient loss estimate. UORO's estimate of the loss gradient ∂ t+1 /∂ θ at time t + 1 is expressed by plugging into (4) the rank-one approximation ∂s t /∂θ ≈ st ⊗ θt , which results in

∂ t+1 ∂o (o t+1 , o * t+1 ) ∂F out ∂s (x t+1 , s t , θ) • st θt + ∂ t+1 ∂o (o t+1 , o * t+1 ) ∂F out ∂θ (x t+1 ,
∂F state ∂s (x t+1 , s t , θ) • st = lim ε→0 F state (x t+1 , s t + ε st , θ) -F state (x t+1 , s t , θ) ε (17)
computable through two applications of F state . This operation is referred to as tangent forward propagation Simard et al. (1991) and can also often be computed algebraically.

This allows for complete implementation of one step of UORO (Alg. 1). The cost of UORO (including running the model itself) is three applications of F state , one application of F out , one backpropagation through F out and F state , and a few elementwise operations on vectors and scalar products.

The resulting algorithm is detailed in Alg. 1. F.forward(v) denotes pointwise application of F at point v, F.backprop(v, δo) backpropagation of row vector δo through F at point v, and F.forwarddiff(v, δv) tangent forward propagation of column vector δv through F at point v. Notably, F.backprop(v, δo) has the same dimension as v , e.g. F out .backprop((x t+1 , s t , θ), δo t+1 ) has three components, of the same dimensions as x t+1 , s t and θ .

The proposed update rule for stochastic gradient descent (3) can be directly adapted to other optimizers, e.g. Adaptative Momentum (Adam) [START_REF] Diederik | Adam: A method for stochastic optimization[END_REF] or Adaptative Gradient [START_REF] Duchi | Adaptive subgradient methods for online learning and stochastic optimization[END_REF]. Vanilla stochastic gradient descent (SGD) and Adam are used hereafter. In Alg. 1, such optimizers are denoted by SGDOpt and the corresponding parameter update given current parameter θ, gradient estimate g t and learning rate η t is denoted SGDOpt.update(g t , η t , θ).

MEMORY-T UORO AND RANK-k UORO

The unbiased gradient estimates of UORO injects noise via ν, thus requiring smaller learning rates.

To reduce noise, UORO can be used on top of truncated BPTT so that recent gradients are computed exactly.

Formally, this just requires applying Algorithm 1 to a new transition function F T which is just T consecutive steps of the original model F . Then the backpropagation operation in Algorithm 1 becomes a backpropagation over the last T steps, as in truncated BPTT. The loss of one step of F T is the sum of the losses of the last T steps of F , namely t+T t+1 := t+T k=t+1 k . Likewise, the forward tangent propagation is performed through F T . This way, we obtain an unbiased gradient estimate in which the gradients from the last T steps are computed exactly and incur no noise. The resulting algorithm is referred to as memory-T UORO. Its scaling in T is similar to T -truncated BPTT, both in

, θ), (0, st , 0)) (_, _, δθ g ) ← F state .backprop((x t+1 , s t , θ), ν ) /* compute normalizers */ ρ 0 ← θt st+1 + ε + ε , ρ 1 ← δθ g ν + ε + ε with ε = 10 -7 /* reduce */ st+1 ← ρ 0 st+1 + ρ 1 ν, θt+1 ← θt ρ 0 + δθ g ρ 1 /* update θ */ SGDOpt.update(g t+1 , η t+1 , θ)
terms of memory and computation. In the experiments below, memory-T UORO reduced variance early on, but did not significantly impact later performance.

The noise in UORO can also be reduced by using higher-rank gradient estimates (rank-r instead of rank-1), which amounts to maintaining r distinct values of s and θ in Algorithm 1 and averaging the resulting values of g. We did not exploit this possibility in the experiments below, although r = 2 visibly reduced variance in preliminary tests.

UORO'S VARIANCE IS STABLE AS TIME GOES BY

Gradient-based sequential learning on an unbounded data stream requires that the variance of the gradient estimate does not explode through time. UORO is specifically built to provide an unbiased estimate whose variance does not explode over time.

A precise statement regarding UORO's convergence and boundedness of the variance of gradients is provided in [START_REF] Massé | Autour de l'Usage des gradients en apprentissage statistique[END_REF]. Informally, when the largest eigenvalue of the differential transition operator ∂F state /∂s is uniformly bounded by a constant δ < 1 (which characterizes stable dynamical systems), the normalizing factors in ( 14) and ( 15) enforce that the influence of previous ν's decrease exponentially with time.

We hereby provide an experimental validation of the boundedness of UORO's variance in Fig. 1a. To monitor the variance of UORO's estimate over time, a 64-unit GRU recurrent network is trained on the first 10 7 characters of the full works of Shakespeare using UORO. The network is then rerun 100 times on the 10000 first characters of the text, and gradients estimates at each time steps are computed, but not applied. The gradient relative variance, that is

E g t -E [g t ] 2 E [g t ] 2 , ( 18
)
is computed, where the average is taken with respect to runs. This quantity appears to be stationary over time (Fig. 1a).

UORO'S VARIANCE INCREASES WITH THE NUMBER OF HIDDEN UNITS

As the number of hidden units in the recurrent network increases, the rank one approximation that is used to provide an unbiased gradient estimate becomes coarser. Consequently, the relative variance, as defined in (18), should increase as the number of hidden units increases.

This increase is experimentally verified in Fig. 1b. Untrained GRU networks with various number of units are run for 10 timesteps, 100 times for each size, and the UORO gradient estimate after these 10 timesteps is computed (but not applied). The relative variance of these gradients over the 100 runs is evaluated, for each network size. As shown in the figure, the relative variance increases with the number of units. Note the horizontal log scale.

The increase of the variance of the estimate with network size underlines the need for smaller learning rates when training large networks with UORO, compared to truncated backpropagation. This can imply slower learning for the kind of dependencies that truncated backpropagation can learn. The need for lower learning rates with larger networks is exemplified in Fig. 1c. GRU networks of various hidden sizes are trained with UORO on a simple copy task, as presented in [START_REF] Hochreiter | Long short-term memory[END_REF], with a lag of T = 5. The networks are all trained with the same decreasing learning rate, η t = 10 -4 1+3•10 -3 t . For all network sizes except the largest, the error decreases slowly but steadily. For the largest network, the variance is too large compared to the learning rate, and the error jumps sharply midway through.

EXPERIMENTS ILLUSTRATING TRUNCATION BIAS

The set of experiments below aims at displaying specific cases where the biases from truncated BPTT are likely to prevent convergence of learning. On this test set, UORO's unbiasedness provides steady convergence, highlighting the importance of unbiased estimates for general recurrent learning.

Influence balancing. The first test case exemplifies learning of a scalar parameter θ which has a positive influence in the short term, but a negative one in the long run. Short-sightedness of truncated algorithms results in abrupt failure, with the parameter exploding in the wrong direction, even with truncation lengths exceeding the temporal dependency range by a factor of 10 or so.

1 × 10 -12 1 × 10 -10 1 × 10 -8 1 × 10 -6 1 × 10 -4 1 × 10 -2 1 × 10 0 1 × 10 2 1 × 10 4 1 ×
with A a square matrix of size n with A i,i = 1/2, A i,i+1 = 1/2, and 0 elsewhere; θ ∈ R is a scalar parameter. The second term has p positive-θ entries and np negative-θ entries. Intuitively, the effect of θ on a unit diffuses to shallower units over time (Fig. 3a). Unit i only feels the effect of θ from unit i + n after n time steps, so the intrinsic time scale of the system is ≈ n. The loss considered is a target on the shallowest unit s 1 ,

t = 1 2 (s 1 t -1) 2 . ( 20
)
Learning is performed online with vanilla SGD, using gradient estimates either from UORO or Ttruncated BPTT with various T . Learning rates are of the form η t = η 1+ √ t for suitable values of η.

As shown in Fig. 2a, UORO solves the problem while T -truncated BPTT fails to converge for any learning rate, even for truncations T largely above n. Failure is caused by ill balancing of time dependencies: the influence of θ on the loss is estimated with the wrong sign due to truncation. For n = 23 units, with 13 minus signs, truncated BPTT requires a truncation T ≥ 200 to converge.

Next-character prediction. The next experiment is character-level synthetic text prediction: the goal is to train a recurrent model to predict the t + 1-th character of a text given the first t online, with a single pass on the data sequence.

A single layer of 64 units, either GRU or LSTM, is used to output a probability vector for the next character. The cross entropy criterion is used to compute the loss.At each time t we plot the cumulated loss per character on the first t characters, 1 t t s=1 s . (Losses for individual characters are quite noisy, as not all characters in the sequence are equally difficult to predict.) This would be the compression rate in bits per character if the models were used as online compression algorithms on the first t characters. In addition, in Table 1 we report a "recent" loss on the last 100, 000 characters, which is more representative of the model at the end of learning.

Optimization was performed using Adam with the default setting β 1 = 0.9 and β 2 = 0.999, and a decreasing learning rate η t = γ 1+α √ t , with t the number of characters processed. As convergence of UORO requires smaller learning rates than truncated BPTT, this favors UORO. Indeed UORO can fail to converge with non-decreasing learning rates, due to its stochastic nature.

DISTANT BRACKETS DATASET (s, k, a). The distant brackets dataset is generated by repeatedly outputting a left bracket, generating s random characters from an alphabet of size a, outputting a right bracket, generating k random characters from the same alphabet, repeating the same first s characters between brackets and finally outputting a line break. A sample is shown in Fig. 3b.

UORO is compared to 4-truncated BPTT. Truncation is deliberately shorter than the inherent time range of the data, to illustrate how bias can penalize learning if the inherent time range is unknown a priori. The results are given in Fig. 2b (with learning rates using α = 0.015 and γ = 10 -3 ). UORO beats 4-truncated BPTT in the long run, and succeeds in reaching near optimal behaviour both with GRUs and LSTMs. Truncated BPTT remains stuck near a memoryless optimum with LSTMs; with GRUs it keeps learning, but at a slow rate. Still, truncated BPTT displays faster early convergence. Plots for a few setups are given in Fig. 4. The learning rates used α = 0.03 and γ = 10 -3 . Numerical results at the end of training are given in Table 1. For reference, the true entropy rate is 0.14 bpc, while the entropy rate of a model that does not understand that the numbers of a's and b's coincide is double, 0.28 bpc.

Here, in every setup, UORO reliably converges and reaches near optimal performance. Increasing UORO's range does not significantly improve results: providing an unbiased estimate is enough to provide reliable convergence in this case. Meanwhile, truncated BPTT performs inconsistently. Notably, with GRUs, it either converges to a poor local optimum corresponding to no understanding of the temporal structure, or exhibits gradient reascent in the long run. Remarkably, with LSTMs rather than GRUs, 16-truncated BPTT reliably reaches optimal behavior on this problem even with biased gradient estimates.

CONCLUSION

We introduced UORO, an algorithm for training recurrent neural networks in a streaming, memoryless fashion. UORO is easy to implement, and requires as little computation time as truncated BPTT, at the cost of noise injection. Importantly, contrary to most other approaches, UORO scalably provides unbiasedness of gradient estimates. Unbiasedness is of paramount importance in the current theory of stochastic gradient descent. Furthermore, UORO is experimentally shown to benefit from its unbiasedness, converging even in cases where truncated BPTT fails to reliably achieve good results or diverges pathologically.

  Figure 1: (a) The relative variance of UORO gradient estimates does not significantly increase with time. Note the logarithmic scale on the time axis. (b) The relative variance of UORO gradient estimates significantly increases with network size. Note the logarithmic scale on number of units. (c) Variance of larger networks affects learning on a small range copy task.

  Figure 2: (a)Results for influence balancing with 23 units and 13 minus; note the vertical log scale. (b)Learning curves on distant brackets (1, 5, 10).

  [a]eecbe[a] [j]fbfjd[j] [c]bgddc[c] [d]gjhai[d] [e]iaghb[e] [h]bigaj[h](b) Distant brackets (1, 5, 10).

  a n b n (1, 32).

Figure 3 :

 3 Figure 3: Datasets.

Figure 4 :

 4 Figure 4: Learning curves on a n b n (1,32)

a

  n b n (k, l) DATASET. The a n b n (k, l) dataset tests memory and counting Gers & Schmidhuber(2001); it is generated by repeatedly picking a random number n between k and l, outputting a string of n a's, a line break, n b's, and a line break (see Fig.3c). The difficulty lies in matching the number of a's and b's.

  For instance, a standard recurrent neural network, with hidden state s t (preactivation values) and output o t at time t, is described with the update equationsF state (x t+1 , s t , θ) := W x x t+1 + W s tanh(s t ) + b and F out (x t+1 , s t , θ) := W o tanh(F state (x t+1 , s t , θ)) + b o ; here the parameter is θ = (W x , W s , b, W o , b o ), and a typical loss might be s (o s , o * s ) := (o so * s ) 2 .Optimization by gradient descent is standard for neural networks. In the spirit of stochastic gradient descent, we can optimize the total loss L =

		t (o t , o * t ) one term at a time and update the
	0≤t≤T
	parameter online at each time step via
	θ ← θ -η t	∂ t ∂θ	(3)
	where η t is a scalar learning rate at time t. (Other gradient-based optimizers can also be used, once
	∂ t ∂θ is known.) The focus is then to compute, or approximate, ∂ t ∂θ .
	BPTT computes ∂ t

0≤t≤T t (o t , o * t ), where o * t is a target output at time t.

  s t , θ). (16) Backpropagating ∂ t+1 /∂o t+1 once through F out returns (∂ t+1 /∂o t+1 • ∂F out /∂x t+1 , ∂ t+1 /∂o t+1 • ∂F out /∂s t , ∂ t+1 /∂o t+1 • ∂F out /∂θ), thus providing all necessary terms to compute (16).

Updating s and θ requires applying (

14

)-(

15

) at each step. Backpropagating the vector of random signs ν once through F state returns _, _, ν ∂F state (x t+1 , s t , θ)/∂θ , providing for (15).

Updating s via (14) requires computing (∂F state /∂s t ) • st . This is computable numerically through

  Algorithm 1 -One step of UORO (from time t to t + 1)Inputs:x t+1 , o * t+1 , s t and θ: input, target, previous recurrent state, and parameters st column vector of size state, θt row vector of size params such that E st ⊗ θt = ∂s t /∂θ -SGDOpt and η t+1 : stochastic optimizer and its learning rate Outputs:-t+1 , s t+1 and θ: loss, new recurrent state, and updated parameters -st+1 and θt+1 such that E st+1 ⊗ θt+1 = ∂s t+1 /∂θ -gt+1 such that E gt+1 = ∂ t+1 /∂θ /* compute next state and loss */ s t+1 ← F state .forward(x t+1 , s t , θ), o t+1 ← F out .forward(x t+1 , s t , θ) Draw ν, column vector of random signs ±1 of size state st+1 ← F state .forwarddiff((x t+1 , s t

	t+1 ← (o t+1 , o * t+1 )
	/* compute gradient estimate */
	(_, δs, δθ) ← F out .backprop (x t+1 , s t , θ), gt+1 ← (δs • st ) θt + δθ	∂ t+1 ∂o t+1
	/* prepare for reduction */

Table 1 :

 1 Averaged loss on the 10 5 last iterations on a n b n (1, 32).

		Truncation	LSTM GRU
		No memory (default) 0.147 0.155
	UORO	Memory-2	0.149 0.174
		Memory-16	0.154 0.149
		1	0.178 0.231
	Truncated BPTT	2	0.149 0.285
		16	0.144 0.207

Arguably, truncated BPTT might still learn some dependencies beyond its truncation range, by a mechanism similar to Echo State Networks[START_REF] Jaeger | Tutorial on training recurrent neural networks[END_REF]. However, truncated BPTT's gradient estimate has a marked bias towards short-term rather than long-term dependencies, as shown in the first experiment of Section 6.

2 Storage of past activations can be reduced, e.g.[START_REF] Gruslys | Memory-efficient backpropagation through time[END_REF]. However, storage of all past inputs is necessary.

Variance minimization is not used at this step, since computing w i v i for every i is not scalable.

In practice, since θ changes during learning, unbiasedness only holds exactly in the limit of small learning rates. This is not specific to UORO as it also affects RTRL.