Hamza Labiad

Martin Fournier

Laura A Mertens

Alexandre Faure

David Carty

Thierry Stoecklin

Piotr Jankowski

Krzysztof Szalewicz

Sébastien D Le Picard

Ian R Sims

Absolute measurements of state-to-state rotational energy transfer between CO and H 2 at interstellar temperatures

Experimental measurements and theoretical calculations of state-to-state rate coefficients for rotational energy transfer of CO in collision with H2 are reported at the very low temperatures prevailing in dense interstellar clouds (5 -20 K). Detailed agreement between quantum state-selected experiments performed in cold supersonic flows using time-resolved infrared -vacuum-ultraviolet double resonance spectroscopy and close-coupling quantum scattering calculations confirms the validity of the calculations for collisions between the two most abundant molecules in the interstellar medium.

Carbon monoxide (CO) was first identified in the interstellar medium (ISM) in 1970 towards the Orion nebula [START_REF] Wilson | Carbon monoxide in the Orion Nebula[END_REF] and is the most abundant molecule in the dense ISM after molecular hydrogen (H 2 ). Despite its small dipole moment (0.11 D), it is the easiest molecule to detect in galaxies and is used to estimate the masses of interstellar molecular clouds, whose densest regions (so called prestellar cores) are the progenitors of stars. CO is also the main carbon reservoir in these "star-forming" regions, where the temperature is ∼ 10 K. Emissions from CO rotational state j = 1 to j = 0 at 115 GHz and from j = 2 to j = 1 at 231 GHz are used as probes of such cold environments and, more generally, as tracers of molecular hydrogen which is invisible to radio telescopes mainly because it has no dipole moment. The dominant cooling mechanism in molecular clouds is collisional excitation of CO followed by rotational line emission, which enables the gravitational collapse of the gas that ultimately leads to the formation of stars and planetary systems.

CO molecules are not in local thermodynamic equilibrium in interstellar clouds owing to the very low hydrogen density (∼ 10 4 cm -3 ). To understand the observed CO spectra and thus the cooling rates of molecular clouds, radiative transfer calculations of rotational state populations of CO must be made. Essential input to these calculations includes absolute rate coefficients for the rotational energy transfer (RET) of CO in collisions with the dominant collider, H 2 , at the temperatures found in molecular clouds (∼ 5 -20 K). These state-to-state rate coefficients are typically obtained by quantum scattering calculations on an ab-initio potential energy surface (PES) describing the CO-H 2 interaction. At low temper- * ian.sims@univ-rennes1.fr atures, the absolute values of these rate coefficients become especially sensitive to the details of the PES, but direct experimental measurements on an absolute scale at very low temperatures are completely lacking. Indeed, to the authors' knowledge no such measurements involving H 2 exist for any diatomic or polyatomic collision partner.

This letter reports direct, absolute measurements of total and state-to-state rate coefficients for the RET of CO in normal -H 2 (i.e., a para-H 2 to ortho-H 2 ratio of 1:3) at 5.5 K, 10 K and 20 K. The results are compared to state-of-the-art close-coupling quantum scattering calculations.

The CO-H 2 collisional system has been investigated theoretically [START_REF] Green | Rotational excitation of CO by collisions with He, H, and H2 under conditions in interstellar clouds[END_REF][START_REF] Jankowski | Ab initio potential energy surface and infrared spectra of H2-CO and D2-CO van der Waals complexes[END_REF][START_REF] Jankowski | A new ab initio interaction energy surface and high-resolution spectra of the H2-CO van der Waals complex[END_REF][START_REF] Yang | Rotational quenching of CO due to H2 collisions[END_REF][START_REF] Jankowski | Theory untangles the high-resolution infrared spectrum of the ortho H2-CO van der Waals complex[END_REF][START_REF] Jankowski | A comprehensive experimental and theoretical study of H2-CO spectra[END_REF][START_REF] Chefdeville | Experimental and theoretical analysis of low-energy CO + H2 inelastic collisions[END_REF][START_REF] Yang | Quantum dynamics of CO-H2 in full dimensionality[END_REF][START_REF] Yang | Fulldimensional quantum dynamics of CO in collision with H2[END_REF][START_REF] Faure | On the importance of full-dimensionality in low-energy molecular scattering calculations[END_REF][START_REF] Castro | Inelastic cross sections and rate coefficients for collisions between co and h2[END_REF] and experimentally [START_REF] Chefdeville | Experimental and theoretical analysis of low-energy CO + H2 inelastic collisions[END_REF][START_REF] Bréchignac | Rotational relaxation of CO by collisions with H2 molecules: A comparison between theory and experiment[END_REF][START_REF] Stoecklin | Comparative experimental and theoretical study of the rotational excitation of CO by collision with ortho-and para-D2 molecules[END_REF] in the past. The lowest temperature at which absolute rate coefficients have been measured is 77 K [START_REF] Bréchignac | Rotational relaxation of CO by collisions with H2 molecules: A comparison between theory and experiment[END_REF]. Recently, six-dimensional analytic CO-H 2 potential energy surfaces (PESs) were obtained by Yang et al. [START_REF] Yang | Quantum dynamics of CO-H2 in full dimensionality[END_REF] and by Faure et al. [START_REF] Faure | On the importance of full-dimensionality in low-energy molecular scattering calculations[END_REF]. The latter (known as V 15 ) was calculated using coupled cluster singles, doubles, triples and non-iterative quadruples [CCSDT(Q)] theory with very large basis sets, and was used in full-dimensional time-independent quantum scattering calculations. In contrast to Yang et al. [START_REF] Yang | Quantum dynamics of CO-H2 in full dimensionality[END_REF], Faure et al. [START_REF] Faure | On the importance of full-dimensionality in low-energy molecular scattering calculations[END_REF] concluded that the full 6D approach was not needed at low temperature. Using a 4D treatment in which the 6D potential was averaged over ground-state monomer vibrations, an equally good agreement was obtained with the relative state-to-state cross sections measured in the crossed-beam experiments of Chefdeville et al. [START_REF] Chefdeville | Experimental and theoretical analysis of low-energy CO + H2 inelastic collisions[END_REF]. Excellent agreement between theory and experiment was also observed in CO+D 2 crossedbeam experiments [START_REF] Stoecklin | Comparative experimental and theoretical study of the rotational excitation of CO by collision with ortho-and para-D2 molecules[END_REF] (using an earlier V 12 potential [START_REF] Jankowski | Theory untangles the high-resolution infrared spectrum of the ortho H2-CO van der Waals complex[END_REF]) and in the measurement of CO-H 2 virial coefficients [START_REF] Garberoglio | All-dimensional H2-CO potential: Validation with fully quantum second virial coefficients[END_REF].

The current measurements were made using a CRESU (Cinétique de Réaction en Ecoulement Supersonique Uniforme or Reaction Kinetics in Uniform Supersonic Flow) apparatus [START_REF] Rowe | Study of the reactions N + 2 + 2N2 → N + 4 + N2 and O + 2 + 2O2 → O + 4 + O2 from 20 to 160 K by the CRESU technique[END_REF][START_REF] Sims | Ultralow temperature kinetics of neutral-neutral reactions -the technique and results for the reactions CN + O2 down to 13 K and CN + NH3 down to 25 K[END_REF]. The CRESU technique uses the supersonic expansion of a buffer gas through convergent-divergent Laval nozzles to produce uniform supersonic flows of cold gas. Each nozzle works under specific gaspressure conditions to cool the gas to a given temperature. Special double-walled nozzles designed for H 2 buffer gas allow precooling to 77 K by liquid nitrogen to obtain a flow temperature as low as 5.5 K. These have only been used once before, down to 11 K by the Rennes group to study the F + H 2 reaction [START_REF] Tizniti | Measurement of the rate of the F + H2 reaction at very low temperatures[END_REF]. Further details of these special Laval nozzles and their characterization can be found in the Supplemental Material [START_REF]IRVUVDR technique, low temperature uniform supersonic flows of hydrogen gas, and spectroscopic[END_REF]. CO molecules are prepared in selected initial rovibrational states (v = 2, j i = 0, 1 or 4) by absorption of pulsed infrared radiation around 2350 nm, generated by a combined optical parametric oscillator -optical parametric amplifier (LaserVision) pumped by an injectionseeded Nd:YAG laser (Continuum Surelite EX) operating at a repetition rate of 10 Hz. Collisions with H 2 in the cold flow on the microsecond timescale of the experiment cause population transfer to final rovibrational states (v = 2, j f ). These states are probed by pulsed laser-induced fluorescence (LIF) via the A 1 Π(v = 0) ← X 1 Σ + (v = 2) vacuum ultraviolet (VUV) vibronic band at 166 nm. Further details of this so-called infrared -vacuum ultraviolet double resonance (IRVUVDR) technique [START_REF] Carty | Rotational energy transfer in collisions between CO(X 1 Σ + , v = 2, J = 0, 1, 4 and 6) and He at temperatures from 294 to 15[END_REF][START_REF] Mertens | Rotational energy transfer in collisions between CO and Ar at temperatures from 293 to 30 K[END_REF], along with a schematic diagram of the experimental arrangement, may be found in the Supplemental Material. Using this IRVUVDR technique, two types of experiment are possible, kinetic and spectroscopic.

p u m p -p
In kinetic experiments, an initially prepared j i state is probed and the LIF intensity is measured as a function of time, as shown in Fig. 1. The rate coefficient for total RET from j i to all j f states, k tot , is extracted from the resulting exponential decay, taking into account transfer back into the probed state as described in previous work [START_REF] Carty | Rotational energy transfer in collisions between CO(X 1 Σ + , v = 2, J = 0, 1, 4 and 6) and He at temperatures from 294 to 15[END_REF][START_REF] Mertens | Rotational energy transfer in collisions between CO and Ar at temperatures from 293 to 30 K[END_REF][START_REF] James | A combined experimental and theoretical study of rotational energy transfer in collisions between NO(X 1 Π 1/2 , v = 3, J) and He, Ar and N2 at temperatures down to 7 K[END_REF].

In spectroscopic experiments, an initial j i state is prepared and the wavelength of the VUV laser is scanned to acquire a spectrum that probes all possible j f states. As shown in Fig. 1, two spectra are obtained, one at 15-50 ns after the IR pulse, when only a small amount of RET has occurred, and one at 2-3 µs when thermal equilibrium of rotational state populations has been established. Voigt lineshapes are fitted to each rotational line in both spectra to extract line intensities, which are then used to obtain state-to-state rate coefficients according to the formula:

k ji→j f = I δt j f f eq j f I eq j f 1 δt[H 2 ] ( 1
)
where I δt j f is the line intensity corresponding to the final rotational state j f at short time delay δt, f eq j f is the Boltzmann factor (fractional population at rotational equilibrium) for j f , I eq j f is the line intensity for j f at equilibrium, and [H 2 ] is the H 2 concentration. The validity of this approach has been discussed in previous work [START_REF] Carty | Rotational energy transfer in collisions between CO(X 1 Σ + , v = 2, J = 0, 1, 4 and 6) and He at temperatures from 294 to 15[END_REF][START_REF] Mertens | Rotational energy transfer in collisions between CO and Ar at temperatures from 293 to 30 K[END_REF][START_REF] James | A combined experimental and theoretical study of rotational energy transfer in collisions between NO(X 1 Π 1/2 , v = 3, J) and He, Ar and N2 at temperatures down to 7 K[END_REF], but briefly, it relies on the short time delay being small enough that transfer via multiple levels may be neglected, while the long time delay is sufficient to ensure full rotational relaxation, but still short enough to avoid significant vibrational relaxation. The latter is assured as vibrational relaxation in this system is many orders of magnitude slower than rotational relaxation [START_REF] Yang | Quantum dynamics of CO-H2 in full dimensionality[END_REF]. CO-CO self-relaxation can be neglected as [CO]/[H 2 ] is maintained at ≲0.05.

Full-dimensional scattering calculations were performed with the DIDIMAT code [START_REF] Faure | On the importance of full-dimensionality in low-energy molecular scattering calculations[END_REF] using the 6D V 15 PES for the v = 0, 1 and 2 states of CO. These were compared to close-coupling calculations performed using the MOLSCAT [START_REF] Hutson | MOLSCAT: MOLecular SCATtering[END_REF] scattering code. In the latter, the rigidrotor approximation was applied and the 6D V 15 PES was averaged over the v = 2, j = 0 vibrational wave function of CO and the v = 0, j = 0 vibrational wave function of H 2 to give a 4D PES, ⟨V 15 ⟩ 20 . In both types of calculations, the (time-independent) close-coupling approach was used to solve a set of coupled radial equations derived from the Schrödinger equation. In the rigid-rotor equations, the rotational constants for CO and H 2 were fixed at their experimental values B 2 (CO)=1.8875 cm -1 and B 0 (H 2 ) = 59.322 cm -1 [START_REF] Huber | Molecular Spectra and Molecular Structure. IV. Constants of diatomic molecules[END_REF]. Convergence of calculations was carefully checked with respect to the size of the basis sets and the propagation parameters. The rotational basis set of CO thus included levels j = 0 -20 while that of H 2 included j H2 = 0, 2 for para-H 2 and j H2 = 1, 3 for ortho-H 2 . Cross sections σ ji→j f (E c ) were computed as function of the collision energy E c from the appropriate T-matrices, as described in Stoecklin et al. [START_REF] Stoecklin | Comparative experimental and theoretical study of the rotational excitation of CO by collision with ortho-and para-D2 molecules[END_REF]. Fig. 2 shows cross sections as a function of collision energy for j = 0 → 1 RET of CO upon collision with para-H 2 (v = 0, j = 0). The cross sections resulting from full-dimensional scattering calculations for CO(v = 0, 1 and 2) are shown as solid lines, while the data points result from the approximate rigid-rotor (4D) close-coupling calculations for CO(v = 2). They agree extremely well with the full-dimensional scattering calculations. It can also be seen that while the positions and magnitudes of resonance peaks differ for different vibrational states, the background cross-sections are similar in magnitude. Therefore, while it is relevant to perform calculations for the vibrational state populated in the experiments (especially at the lowest temperature), we note that the shift and decrease of resonance peaks will affect the rate coefficient values by a few percent only. In the following, only the rigid-rotor 4D calculations are presented and discussed. The rotational rate coefficients were computed by averaging the corresponding cross sections over FIG. 3. Measured state-to-state rate coefficients (squares) and theoretical state-to-state rate coefficients (circles), calculated using the 4D ⟨V15⟩20 PES for ji = 0 → j f and ji = 1 → j f at 5. Table I and Fig. 3 compare measured state-to-state rate coefficients and theoretical values (calculated using the ⟨V 15 ⟩ 20 PES) for j i = 0 → j f and j i = 1 → j f at 5.5 K, 10 K and 20 K. Table I also displays values for j i = 4 → j f calculated at 5.5 K, 10 K and 20 K and measured at 20 K. At 5.5 K and at 10 K, j i = 4 was not accessible in the rovibrational state preparation as rotational states above j = 2 are not populated appreciably at such low temperatures. The agreement between experiment and theory is excellent and the results are all self-consistent within the principle of detailed balance [START_REF]i.e. k1→0 = k0→1 (g1/g0 ) exp (-∆E0-1/kT ) where k0→1 and k1→0 are the forward and backward rate coefficients between levels 0 and 1, g0 and g1 are the degeneracies associated with those levels and ∆E0-1 is the difference in energy. B. H. Mahan, Microscopic reversibility and detailed balance. an analysis[END_REF]. Table I also contains measured total RET rate coefficients, k tot , and the sum of the state-to-state rate coefficients, Σk; these agree within experimental error, with the marginal exception of j i = 0 at 10 K.

Fig. 4 shows a plot of thermally-averaged cross sections for RET from j i = 0 and 1 to all j f states at all temperatures measured. The cross sections were obtained by dividing the measured and theoretical Σk values, for a given j i and temperature, by the mean thermal velocity at that temperature. The same was done using measured k tot values from the kinetic experiments. The cross sections tend to decrease with decreasing temperature. As the temperature, and thus the relative velocity between the CO and H 2 , decreases, the probability of a collision occurring increases as the long-range van der Waals attraction brings the molecules together. This promotes an increase in cross section as the temperature is lowered, so a second mechanism that results in a decrease must dominate. To understand this, consider the stateto-state rate coefficients again. Positive ∆j collisions are endothermic, whereas negative ∆j collisions are exothermic. For j i = 0, RET is purely endothermic. It can be seen in the top row of Fig. 3 that the state-to-state rate coefficients for ∆j = +1 increase with increasing temperature as the proportion of collisions that are energetic enough to cause RET increases. The same is true for ∆j = +2 RET. Only at 20 K is the proportion of collisions energetic enough to cause ∆j = +3 and ∆j = +4 RET. For j i = 1, the rate coefficients for ∆j = -1 RET are almost constant with temperature due to this process being exothermic. For the endothermic, positive ∆j RET from j i = 1, the trend is similar to that for j i = 0, but the evidence for a ∆j = +2 propensity is weak owing to a suppression of ∆j > 1 RET as a result of the greater spacing of rotational energy levels at higher j. The net result is that there is an increase in k (which equals k tot ), and hence cross section, as the temperature increases that more than cancels out the opposite temperature dependence due to the long-range van der Waals attraction.

At 20 K, it can also be seen that the rate coefficient for j i = 0 → j f = 2 is larger than that for j i = 0 → j f = 1. Such a ∆j propensity is due to quantum mechanical effects in scattering involving "almost homonuclear" molecules and is related to the anisotropy of the potential. It can be understood in terms of the partial-wave expansion of the potential V where for a homonuclear diatomic molecule the term l = 1 is missing, so the matrix element of V between the initial j i function and the final j i + 1 function is zero. This selection rule, ∆j = even, becomes only a propensity rule for an "almost homonuclear" molecule like CO. This effect was first investigated by McCurdy and Miller [START_REF] Mccurdy | Interference effects in rotational state distributions -propensity and inverse propensity[END_REF] in the context of classical S-matrix theory.

In conclusion, absolute state-to-state rate coefficients and thermally-averaged cross-sections for rotational energy transfer of CO in collision with normal -H 2 have been measured directly using infrared vacuum ultraviolet double resonance spectroscopy in a CRESU apparatus at interstellar temperatures. The CO states measured were v = 2, j i = 0, 1 and 4 at 20 K and, v = 2, j i = 0 and 1 at 10 K and at 5.5 K. The excellent and detailed agreement on an absolute basis between direct experimental measurements and theoretical calculations using both the 6D V 15 PES and the approximate, 4D ⟨V 15 ⟩ 20 PES provides rare benchmarking of theory. In particular, we believe that the extensive set of collisional data computed by Yang et al. [START_REF] Yang | Rotational quenching of CO due to H2 collisions[END_REF] and based on the previous V 04 PES of Jankowski and Szalewicz [START_REF] Jankowski | A new ab initio interaction energy surface and high-resolution spectra of the H2-CO van der Waals complex[END_REF] can be employed reliably. Indeed, as shown by Faure et al. [START_REF] Faure | On the importance of full-dimensionality in low-energy molecular scattering calculations[END_REF], the difference between V 15 and V 04 has only a moderate impact on the low-energy integral cross sections, and even less so on the rate coefficients (see also the discussion in Chefdeville et al. [START_REF] Chefdeville | Experimental and theoretical analysis of low-energy CO + H2 inelastic collisions[END_REF]). The calculations can now be used with complete confidence in the interpretation of spectra of astronomical CO in star-forming molecular clouds and in the understanding of molecular-cloud cooling that leads to star formation by gravitational collapse.

PFIG. 1 .FIG. 2 .

 12 FIG. 1. Representative IRVUVDR data from the kinetic and spectroscopic experiments at T = 5.5 K and [H2] = 3.2×10 16 cm -3 .Raw data are shown in black. The short delay VUV LIF spectrum of CO(v = 2), with a fit shown in blue, was acquired at a time delay between the IR pump-laser pulse and the VUV probe-laser pulse of 50 ns after initial preparation of the ji = 1 rotational state. The intensities of neighboring rotational lines are proportional to the state-to-state rate coefficients for RET into those levels. The spectrum at long time delay, with a fit shown in red, was acquired once thermal equilibrium of rotational states was reached, 2 µs after the pump-laser pulse. The spectra have been fitted using Voigt line profiles. The decaying signal is the intensity at the peak of the R(1) spectral line as a function of time delay, fit (in red) to a single-exponential decay.

  FIG.3. Measured state-to-state rate coefficients (squares) and theoretical state-to-state rate coefficients (circles), calculated using the 4D ⟨V15⟩20 PES for ji = 0 → j f and ji = 1 → j f at 5.5 K, 10 K and 20 K. Error bars correspond to 2σ.

  Maxwell-Boltzmann velocity distributions at 5.5 K, 10 K and 20 K.

TABLE I

 I

			T = 20 K
			ji
	j f	0	1	4
	0	-	2.4 ± 0.6 (3.2)	1.0 ± 0.8 (0.6)
	1	8.3 ± 0.7 (7.4)	-	2.2 ± 1.2 (1.3)
	2 11.8 ± 0.9 (10.3) 6.5 ± 1.5 (6.2)	8.5 ± 1.1 (8.4)
	3	1.0 ± 0.6 (0.9)	4.0 ± 1.7 (4.4)	7.9 ± 1.1 (7.5)
	4	0.6 ± 0.5 (0.03) 0.5 ± 1.3 (0.3)	-
	5	(0.03)	(0.08)	1.7 ± 0.7 (2.3)
	k 21.7 ± 1.4 (18.7) 13.4 ± 2.7 (14.2) 21.3 ± 2.2 (20.1)
	ktot	16.1 ± 4.0	13.4 ± 2.6	20.1 ± 4.1
			T = 10 K
	j f	0	1	4
	0	-	3.7 ± 0.4 (3.4)	(0.5)
	1	7.6 ± 2.4 (6.0)	-	(1.1)
	2	4.0 ± 1.9 (4.9)	3.8 ± 0.6 (3.9)	(8.3)
	3	(0.02)	1.1 ± 0.7 (1.2)	(8.1)
	k 11.6 ± 3.1 (11.0) 8.6 ± 1.0 (8.5)	(18.0)
	ktot	6.6 ± 0.9	7.9 ± 0.9	-
			T = 5.5 K
	j f	0	1	4
	0	-	4.8 ± 0.6 (3.6)	(0.5)
	1	5.4 ± 0.5 (4.0)	-	(1.0)
	2	1.9 ± 0.5 (1.3)	1.1 ± 0.3 (1.7)	(7.8)
	3	(0.01)	0.08 ± 0.08 (0.1)	(8.5)
	k 7.3 ± 0.7 (5.3)	6.0 ± 0.7 (5.4)	(17.8)
	ktot	5.8 ± 1.1	6.0 ± 0.7	-

. Measured state-to-state rate coefficients and theoretical values (in parentheses) for RET in CO -H2 collisions, calculated using the 4D ⟨V15⟩20 PES, for ji = 0, 1, 4 → j f at 5.5 K, 10 K and 20 K. Units are in 10 -11 cm 3 s -1 and error limits correspond to 2σ statistical errors. The sum of the state-to-state rate coefficients, Σk, and the equivalent rate coefficients, ktot, from kinetic experiments are also given.

This work was supported by the Agence Nationale de la Recherche (ANR-HYDRIDES), contract ANR-12-BS05-0011-01, by the Région de Bretagne and the French National Programme "Physique et Chimie du Milieu Interstellaire" (PCMI) of CNRS/INSU with INC/INP cofunded by CEA and CNES. DC would like to thank the Université de Rennes 1 for support via the Visiting Professor programme. This work was also supported in part by the National Science Foundation under Grant No. CHE-1900551. PJ acknowledges support from Polish National Science Centre grant no. 2017/25/B/ST4/01300.


This document serves as a reference for the formal specification of the HILECOP model-to-text transformation. HILECOP is a methodology for the design and synthesis of safety-critical digital systems [START_REF] Andreu | A Distributed Architecture for Activating the Peripheral Nervous System[END_REF]. Many elements used in the following definitions are introduced in the following thesis [START_REF] Iampietro | Vérification formelle d'une méthodologie pour la conception et la production de systèmes numériques critiques[END_REF].

Preliminary definitions

The HILECOP model-to-text transformation (HM2T) generates a H-VHDL design [2, p.65] out of an SITPN model [2, p.37]. It also generates a structure that relates the elements of the SITPN model to the elements of the H-VHDL design. This structure is called a SITPN-to-H-VHDL binder. Its formal definition is as follows:

Definition 1 (SITPN-to-H-VHDL design binder). Given a sitpn ∈ SIT P N and a H-VHDL design d ∈ design, a SITPN-to-H-VHDL binder γ ∈ W M(sitpn, d) is a tuple <P Map, T Map, CMap, AF Map> where:

• P Map ∈ P → {id p | comp(id p , place, g, i, o) ∈ d.beh}

• T Map ∈ T → {id t | comp(id t , transition, g, i, o) ∈ d.beh}

• CMap ∈ C → {id c | (in, id c , bool) ∈ d.ports}
• AF Map ∈ A ∪ F → {id af | (out, id af , bool) ∈ d.ports} Notation 1. For a given binder γ and an element of an SITPN structure e ∈ P ∪ T ∪ C ∪ A ∪ F , we write γ(e) where e is looked up in the appropriate function. For instance, for a given f ∈ F , γ(f ) is a shorthand notation for AF Map(f ) where γ = < . . . , AF Map>.

While presenting the specification of the HM2T, we refer to two pre-defined designs, namely: the place and transition designs. The details of the implementation of the place and transition designs in H-VHDL are given in Appendices A and B. The HILECOP design store has been specifically defined for the elaboration and simulation the H-VHDL designs generated by the HM2T. It holds the definitions of the place an transition designs.

3 Formal specification

In the following section, we present the formal specification of the HM2T. The specification is a relation between the inputs of the HM2T, namely a SITPN model and a bounding function, and the possible outputs which can be either a pair composed of a H-VHDL design and a SITPN-to-H-VHDL binder or the error value. The relation is written

HM2T spec ⊆ SIT P N × (P → N) × ((design × W M(sitpn, d)) ∪ {err}).

The number of input ports of the design is equal to the number of conditions of the SITPN model

: |{id | (in, id, τ ) ∈ d.ports}| = |C|.
3. 12. For all place of the SITPN model with no output transition, the input port map and output port map of the corresponding PDI includes the following associations:

∀p ∈ P, g p , i p , o p , output(p) = ∅ ⇒ comp(γ(p), place, g p , i p , o p ) ∈ d.beh ⇒ {(oaw(0) ⇒ 0), (oat(0) ⇒ basic), (otf(0) ⇒ false)} ⊆ i p ∧ {(oav ⇒ open), (pah ⇒ open), (rtt ⇒ open)} ⊆ o p .
13. For all place of the SITPN model with no action, the marked output port is left unconnected in the output port map of the corresponding PDI:

∀p ∈ P, g p , i p , o p , acts(p) = ∅ ⇒ comp(γ(p), place, g p , i p , o p ) ∈ d.beh ⇒ (marked ⇒ open) ∈ o p
where acts(p) = {a | A(p, a) = true}, the set of actions associated with place p.

14. For all transition of the input SITPN model, there exists a corresponding TDI identified through γ in the behavior of the output design: ∀t ∈ T, ∃g t , i t , o t s.t. comp(γ(t), transition, g t , i t , o t ) ∈ d.beh.

15. For all transition of the input SITPN model and its corresponding TDI, the generic map of the TDI holds the following associations:

∀t ∈ T, g t , i t , o t , comp(γ(t), transition, g t , i t , o t ) ∈ d.beh ⇒ g t = {(tt ⇒          NOT_TEMPORAL if t / ∈ dom(I s ) TEMPORAL_A_A if I s (t) = [a, a] TEMPORAL_A_B if I s (t) = [a, b] TEMPORAL_A_INF if I s (t) = [a, ∞] ), (mtc ⇒      1 if t / ∈ dom(I s ) b if I s (t) = [a, b] a if I s (t) = [a, ∞] ), (ian ⇒ 1 if input(t) = ∅ |input(t)| otherwise ), (cn ⇒ 1 if conds(t) = ∅ |conds(t)| otherwise )}.
where input(t) = {p | pre(p, t) = ⌊(ω, a)⌋}, the set of input places of transition t, and conds(t

) = {c | C(t, c) = 1 ∨ C(t, c) = -1}
, the set of conditions associated with transition t.

16. For all transition of the input SITPN model and its corresponding TDI, the input port map of the TDI holds the following associations:

∀t ∈ T, g t , i t , o t , comp(γ(t), transition, g t , i t , o t ) ∈ d.beh ⇒ {(A ⇒ 0 if t / ∈ dom(I s ) l(I s (t)) otherwise ), (B ⇒ 0 if t / ∈ dom(I s ) ∨ u(I s (t)) = ∞ u(I s (t)) otherwise )} ⊆ i t .
17. For all transition of the input SITPN model with no input place, the input port map and output port map of the corresponding TDI holds the following associations:

∀t ∈ T, g t , i t , o t , input(t) = ∅ ⇒ comp(γ(t), transition, g t , i t , o t ) ∈ d.beh ⇒ (∃id s s.t. (id s , bool) ∈ d.sigs ∧ (rt(0) ⇒ id s ) ∈ i t ∧ (fired ⇒ id s ) ∈ o t ) ∧ {(iav(0) ⇒ true), (pah(0) ⇒ true)} ⊆ i t .
18. For all transition of the input SITPN model with no condition, the input port map of the corresponding TDI holds the following association:

∀t ∈ T, g t , i t , o t , conds(t) = ∅ ⇒ comp(γ(t), transition, g t , i t , o t ) ∈ d.beh ⇒ (ic(0) ⇒ true) ∈ i t .
19. For all post arc of the input SITPN model, the TDI and PDI corresponding to the source transition and target place of the arc are connected as follows:

∀t ∈ T, p ∈ P, g t , i t , o t , g p , i p , o p , ω ∈ N * , post(t, p) = ⌊ω⌋ ⇒ comp(γ(t), transition, g t , i t , o t ) ∈ d.beh ⇒ comp(γ(p), place, g p , i p , o p ) ∈ d.beh ⇒ ∃i ∈ [0, |input(p)| -1] s.t. (iaw(i) ⇒ ω) ∈ i p ∧ ∃id s s.t. (id s , bool) ∈ d.sigs ∧ (fired ⇒ id s ) ∈ o t ∧ (itf(i) ⇒ id s ) ∈ i p .
20. For all pre arc of the input SITPN model, the PDI and TDI corresponding to the source place and target transition of the arc are connected as follows:

∀t ∈ T, p ∈ P, g t , i t , o t , g p , i p , o p , ω ∈ N * , a ∈ {basic, test, inhib}, pre(p, t) = ⌊(ω, a)⌋ ⇒ comp(γ(t), transition, g t , i t , o t ) ∈ d.beh ⇒ comp(γ(p), place, g p , i p , o p ) ∈ d.beh ⇒ ∃i ∈ [0, |output(p)| -1] s.t. {(oaw(i) ⇒ ω), (oat(i) ⇒ a)} ⊆ i p ∧∃j ∈ [0, |input(t)| -1], id av , id rt , id f rd , id pah s.t. {(id av , bool), (id rt , bool), (id f rd , bool), (id pah , bool)} ⊆ d.sigs ∧ (oav(i) ⇒ id av ) ∈ o p ∧ (iav(j) ⇒ id av ) ∈ i t ∧ (rtt(i) ⇒ id rt ) ∈ o p ∧ (rt(j) ⇒ id rt ) ∈ i t ∧ (otf(i) ⇒ id f rd ) ∈ i p ∧ (fired ⇒ id f rd ) ∈ o t ∧ (pah(i) ⇒ id f rd ) ∈ o p ∧ (a = test ∨ a = inhib ∨ confl(p) = ∅ ⇒ (pah(j) ⇒ true) ∈ i t ) ∧ (a = basic ∧ confl(p) = ∅ or err ⇒ (pah(j) ⇒ id pah ) ∈ i t ).
where confl ∈ P → 2 T ∪ {err} takes a place p as input and yields either an error or an ordered set of transitions computed as follows:

(a) If all conflicts between the output transitions of p are solved by mutual exclusion, or if the set of conflicting transitions of p is a singleton, then confl returns an empty set.

(b) Otherwise, the function tries to establish a total ordering over the set of conflicting transitions of p w.r.t the firing priority relation:

• If no such ordering can be established (in that case, the firing priority relation is ill-formed, and the input SITPN is not well-defined), confl returns the err value. • Otherwise, the function returns the set in a decreasing priority order.

For all place of the input SITPN model for which conflicts in its output transitions are not solved by mutual exclusion, the port indices of the corresponding PDI reflect the priority order established between

the conflicting output transitions:

∀p ∈ P, t, t ′ ∈ confl(p), g p , i p , o p , g t , i t , o t , g t ′ , i t ′ , o t ′ , t ≻ t ′ ⇒ comp(γ(p), place, g p , i p , o p ) ∈ d.beh ⇒ comp(γ(t), transition, g t , i t , o t ) ∈ d.beh ⇒ comp(γ(t ′ ), transition, g t ′ , i t ′ , o t ′ ) ∈ d.beh ⇒ (∀i, j ∈ N, id f rd , id f rd ′ , id s , id s ′ , name t , name t ′ , id out , (fired ⇒ id f rd ) ∈ o t ⇒ (fired ⇒ id f rd ′ ) ∈ o t ′ ⇒ {(otf(i) ⇒ id f rd ), (otf(j) ⇒ id f rd ′ )} ⊆ i p ⇒ (name t ⇒ id s ) ∈ i t ⇒ (name t ′ ⇒ id s ′ ) ∈ i t ′ ⇒ {(id out (i) ⇒ id s ), (id out (j) ⇒ id s ′ )} ⊆ o p ⇒ i < j)
22. There exists an actions process that assigns a value to the output port representing the activation status of the actions (referred to as action ports) of the input SITPN model: (d) Otherwise, the value of the action port is the result of the Boolean sum between the marked output port of all PDIs representing the places associated with the corresponding action:

∃ss ra ,
∀a ∈ A,

For all association between a transition and a condition, the input port reflecting the Boolean value of the given condition (referred to as a condition port) is connected as follows to the input port map of the TDI corresponding to the associated transition:

∀t ∈ T, g t , i t , o t , c ∈ C, comp(γ(t), transition, g t , i t , o t ) ∈ d.beh ⇒ (C(t, c) = 1 ⇒ ∃i ∈ [0, |conds(t)| -1] s.t. (ic(i) ⇒ γ(c)) ∈ i t ) (C(t, c) = -1 ⇒ ∃i ∈ [0, |conds(t)| -1] s.t. (ic(i) ⇒ not(γ(c))) ∈ i t ).
To conclude the definition of the HM2T spec relation, we have HM2T spec (sitpn, b, err) if sitpn is not a well-defined SITPN model [START_REF] Iampietro | Vérification formelle d'une méthodologie pour la conception et la production de systèmes numériques critiques[END_REF].

A The place design in abstract H-VHDL syntax (in, input_arcs_weights, array (nat(0, 255), 0, sub(input_arcs_number, 1))), 13 (in, output_arcs_types, array (nat(0, 2), 0, sub(output_arcs_number, 1))), 14 (in, output_arcs_weights, array (nat(0, 2), 0, sub(output_arcs_number, 1))), 15 (in, input_transitions_fired, array (bool, 0, sub(input_arcs_number, 1))), 16 (in, output_transitions_fired, array (bool, 0, sub(output_arcs_number, 1))), 17 (out, output_arcs_valid, array (bool, 0, sub(output_arcs_number, 1))), 18 (out, priority_authorizations, array (bool, 0, sub(output_arcs_number, 1))),

C Generic constant and signal names reference

Formal

  specification of the HILECOP model-to-text transformation

Definition 2 (

 2 HILECOP model-to-text transformation specification). For all SITPN model sitpn ∈ SIT P N, bounding function b ∈ P → N, H-VHDL design d ∈ design, and SITPN-to-H-VHDL binder γ ∈ W M(sitpn, d), we have HM2T spec sitpn, b, (d, γ) if: 1. Design d has an empty generic clause: d.gens = ∅.

  ss a s.t. ps(actions, ∅, rst{ss ra }else{falling{ss a }}) ∈ d.beh and (a) The length of the ss ra and ss a sequences is equal to the number of actions of the input SITPN model: |ss ra | ; = |ss a | ; = |A| where |ss| ; = |ss 1 | ; + |ss 2 | ; if ss = ss 1 ; ss 2 1 otherwise (b) During the initialization, all action ports are assigned to false by the actions process: ∀a ∈ A, γ(a) ⇐ false ∈ ss ra (c) An action port corresponding to an action associated with no place is assigned to false during the falling edge phase: ∀a ∈ A s.t. pls(a) = ∅, γ(a) ⇐ false ∈ ss a where pls(a) = {p | A(p, a) = true}, the set of places associated with action a.

  initial_marking, nat(0, maximal_marking)), 12

B

  19(out, reinit_transitions_time, array (bool, 0, sub(output_arcs_number, 1))), 20 (out, marked, bool)}, 21 22 --Internal signal declarations 23 sigs={(s_input_tokens_sum, nat(0, maximal_marking)), 24 (s_marking, nat(0, maximal_marking)), 25 (s_output_tokens_sum, nat(0, maximal_marking))}, , nat(0, maximal_marking))}, 31 v_tokens_sum := 0; 32 for (i, 0, sub(input_arcs_number, 1)){ 33 if (input_transitions_fired(i)) { 34 v_tokens_sum := add(v_tokens_sum, input_arcs_weights(i)) 35 } else { null } , 0, sub(output_arcs_number, 1)) { 43 if (and(output_transitions_fired(i), eq(output_arcs_types(i), 0))) { 44 v_tokens_sum := add(v_tokens_sum, output_arcs_weights(i)) 45 } else { null } 46 }; 47 s_output_tokens_sum ⇐ v_tokens_sum) 8 The transition design in abstract H-VHDL syntax input_conditions, array(bool, 0, sub(conditions_number, 1)), 13 (in, time_A_value, nat(0, maximal_time_counter)), 14 (in, time_B_value, nat(0, maximal_time_counter)), 15 (in, input_arcs_valid, array(bool, 0, sub(input_arcs_number, 1))), 16 (in, reinit_time, array(bool, 0, sub(input_arcs_number, 1))), 17 (in, priority_authorizations, array(boolean, 0, sub(input_arcs_number, 1))), 18 (out, fired, bool)}, , {(v_internal_condition, bool)}, 33 v_internal_condition := true; 34 for (i, 0, sub(conditions_number, 1)) { 35 v_internal_condition := and(v_internal_condition, input_conditions(i)) , {(v_internal_enabled, bool)}, 40 v_internal_enabled := true; 41 for (i, 0, sub(input_arcs_number, 1)) { 42 v_internal_enabled := and(v_internal_enabled, input_arcs_valid(i)) (reinit_time_counter_evaluation, {(v_internal_reinit_time_counter, bool)}, 47 v_internal_reinit_time_counter := false; 11 and ge(s_time_counter, sub(time_A_value, 1)) 77 and (s_time_counter < time_B_value)) 78 or (eq(transition_type, 2) 79 and eq(s_time_counter, sub(time_A_value, 1))) 80 or (eq(transition_type, 3) 81 and ge(s_time_counter, sub(time_A_value, 1))))) 82 or (s_reinit_time_counter 83 and neq(transition_type, 0) 84 and eq(time_A_value, 1)))) 85 || 86 ps (priority_authorization_evaluation, {(v_priority_combination, bool)}, 87 v_priority_combination := true; 88 for (i, 0, sub(input_arcs_number, 1)) { 89 v_priority_combination := and(v_priority_combination, priority_authorizations(i)) s_firable ⇐ false } 95 else { falling { s_firable ⇐ s_firing_condition } }) 96 || 12 ps (fired_evaluation, ∅, fired ⇐ and(s_firable, s_priority_combination)) Listing 2: The transition design in H-VHDL abstract syntax.

  The number of output ports of the design is equal to the number of actions and functions of the SITPN

	7. All the fields of the SITPN-to-H-VHDL binder are bijective functions: P Map(γ) is bijective, T Map(γ) is bijective,. . . 8. For all place of the input SITPN model, there exists a corresponding place design instance (PDI) iden-tified through γ in the behavior of the output design: ∀p ∈ P, ∃g p , i p , o p s.t. comp(γ(p), place, g p , i p , o p ) ∈ d.beh. 9. For all place of the input SITPN model and its corresponding PDI, the generic map of the PDI holds the following associations: ∀p ∈ P, g p , i p , o p , comp(γ(p), place, g p , i p , o p ) ∈ d.beh ⇒ g p = {(mm ⇒ b(p)), (ian ⇒ 1 if input(p) = ∅ |input(p)| otherwise ), (oan ⇒ 1 if output(p) = ∅ |output(p)| otherwise )} where input(p) = {t | post(t, p) = ⌊ω⌋}, the set of input transitions of place p, and output(p) = {t | pre(p, t) = ⌊(ω, a)⌋}, the set of output transitions of place p. 10. For all place of the input SITPN model and its corresponding PDI, there is an association between the im input port and the initial marking of the place in the input port map of the PDI: includes the following associations: ∀p ∈ P, g p , i p , o p , input(p) = ∅ ⇒ comp(γ(p), place, g p , i p , o p ) ∈ d.beh ⇒ ∀p ∈ P, g 11. For all place of the SITPN model with no input transition, the input port map of the corresponding PDI {(iaw(0) ⇒ 0), (itf(0) ⇒ false)} ⊆ i p .

model: |{id | (out, id, τ ) ∈ d.ports}| = |A ∪ F |. 4. Design d is elaborable in the context of the HILECOP design store and given an empty dimensioning function: ∃∆ ∈ ElDesign, σ e ∈ Σ s.t. D H , ∅ ⊢ d elab --→ ∆, σ e . 5. All design instantiation statement in the design's behavior either creates a PDI or a TDI: ∀id c , id e , g, i, o s.t. comp(id c , id e , g, i, o) ∈ d.beh, id e = place ∨ id e = transition. 6. The actions and functions processes are the only two processes in the design's behavior: ∀id p , vars, ss s.t. ps(id p , vars, ss) ∈ d.beh, id p = actions ∨ id p = functions. p , i p , o p , comp(γ(p), place, g p , i p , o p ) ∈ d.beh ⇒ (im ⇒ M 0 (p)) ∈ i p .

  There exists an functions process that assigns a value to the output port representing the execution status of the functions (referred to as function ports) of the input SITPN model: ∃ss rf , ss f s.t. ps(functions, ∅, rst{ss rf }else{rising{ss f }}) ∈ d.beh and (a) The length of the ss rf and ss f sequences is equal to the number of functions of the input SITPN model: |ss rf | ; = |ss f | ; = |F | (b) During the initialization, all function ports are assigned to false by the functions process: ∀f ∈ F , γ(f ) ⇐ false ∈ ss rf (c) An function port corresponding to an function associated with no transition is assigned to false during the rising edge phase: ∀f ∈ F s.t. trs(f ) = ∅, γ(f ) ⇐ false ∈ ss f where trs(f ) = {t | F(t, f ) = true}, the set of transitions associated with function f . (d) Otherwise, the value of the function port is the result of the Boolean sum between the fired output port of all TDIs representing the transitions associated with the considered function:

pls(a) = ∅ ⇒ ∃e or s.t. γ(a) ⇐ e or ∈ ss a ∧ is_bsum(e or , |pls(a)|) ∧∀p ∈ pls(a), g p , i p , o p , comp(γ(p), place, g p , i p , o p ) ∈ d.beh ⇒ ∃id m s.t. (id m , bool) ∈ d.sigs ∧ id m ∈ e or ∧ (marked ⇒ id m ) ∈ o p where is_bsum is defined as follows: e ∈ {id, b} is_bsum(e, 1) is_bsum(e 1 , n) is_bsum(e 2 , m) is_bsum(or(e 1 , e 2 ), n + m) 23. ∀f ∈ F , trs(f ) = ∅ ⇒ ∃e or s.t. γ(f ) ⇐ e or ∈ ss f ∧ is_bsum(e or , |trs(f )|) ∧∀t ∈ trs(f ), g t , i t , o t , comp(γ(t), transition, g t , i t , o t ) ∈ d.beh ⇒ ∃id m s.t. (id m , bool) ∈ d.sigs ∧ id m ∈ e or ∧ (fired ⇒ id m ) ∈ o t

Table 1 :

 1 Constants and signals reference for the H-VHDL transition and place designs. In the Category column, T (resp. P) indicates a generic constant, input port, output port or internal signal defined in the transition (resp. place) design.

	Generic constants and signals reference
	Full name	Alias	Category	Type
	input_arcs_number	ian generic constant (T)	nat(0, NATMAX)
	transition_type	tt	generic constant (T)	{not_temp, temp_a_b,
				temp_a_a, temp_a_inf}
	conditions_number	cn	generic constant (T)	nat(0, NATMAX)
	maximal_time_counter	mtc generic constant (T)	nat(0, NATMAX)
	time_A_value	A	input port (T)	nat(0, mtc)
	time_B_value	B	input port (T)	nat(0, mtc)
	input_conditions	ic	input port (T)	array(bool, 0, cn -1)
	reinit_time	rt	input port (T)	array(bool, 0, ian -1)
	input_arcs_valid	iav	input port (T)	array(bool, 0, ian -1)
	priority_authorizations	pah	input port (T)	array(bool, 0, ian -1)
	fired	f	output port (T)	bool
	s_condition_combination	scc	internal signal (T)	bool
	s_reinit_time_counter	srtc	internal signal (T)	bool
	s_priority_combination	spc	internal signal (T)	bool
	s_firable	sfa	internal signal (T)	bool
	s_enabled	se	internal signal (T)	bool
	s_time_counter	stc	internal signal (T)	nat(0, mtc)
	s_firing_condition	sfc	internal signal (T)	bool
	input_arcs_number	ian generic constant (P)	nat(0, NATMAX)
	output_arcs_number	oan generic constant (P)	nat(0, NATMAX)
	maximal_marking	mm	generic constant (P)	nat(0, NATMAX)
	initial_marking	im	input port (P)	nat(0, mm)
	output_arcs_types	oat	input port (P)	array({basic, test, inhib}, 0, oan -1)
	output_arcs_weights	oaw	input port (P)	array(nat(0, 255), 0, oan -1)
	output_transitions_fired otf	input port (P)	array(bool, 0, oan -1)
	input_arcs_weights	iaw	input port (P)	array(nat(0, 255), 0, ian -1)
	input_transitions_fired	itf	input port (P)	array(bool, 0, ian -1)
	output_arcs_valid	oav	output port (P)	array(bool, 0, oan -1)
	reinit_transitions_time	rtt	output port (P)	array(bool, 0, oan -1)
	priority_authorizations	pah	output port (P)	array(bool, 0, oan -1)
	marked	m	output port (P)	bool
	s_marking	sm	internal signal (P)	nat(0, mm)
	s_output_token_sum	sots	internal signal (P)	nat(0, mm)
	s_input_token_sum	sits	internal signal (P)	nat(0, mm)


Raúl Hernández Sánchez
email: raul.hernandez-sanchez@univ-amu.fr

Bruno Denet

Paul Clavin

Critical Dynamics of Direct Initiation of Spherical Detonations

Introduction

Gaseous detonations can be described as supersonic combustion waves whose structure consists of an inert shock wave followed by a reactive layer [START_REF] Clavin | Combustion Waves and Fronts in Flows: Flames, Shocks, Detonations, Ablation Fronts and Explosion of Stars[END_REF]. The classical Chapman-Jouguet (CJ) theory predicts accurately the overall gas state behind a steady detonation front and its unique self-propagating regime considering the detonation as a discontinuity. However, the unsteady dynamics of the detonation wave are strongly sensitive to small modifications of its inner structure [START_REF] Clavin | Analytical study of the direct initiation of gaseous detonations for small heat release[END_REF]. Hence, anticipate whether a detonation will be initiated based exclusively on the constitutive properties of the reactive mixture represents a rather difficult question. A simplified model describing the critical dynamics of direct initiation of spherical detonations in the small heat release limit is investigated here numerically. This model provides fundamental physics insights into the initiation of detonations phenomenon that may contribute to the definition of more sophisticated design specifications in the context of explosion safety.

Two different modes of initiation of a self-propagating detonation are commonly identified: direct initiation and deflagration to detonation transition (DDT). The direct initiation mode corresponds to a fast transition from a decaying blast wave into a self-sustained detonation, while in the DDT a slower transition from a deflagration to a detonation takes place. The fast mode, addressed here, is thus characterised by the absence of a predetonation deflagration.

Numerous efforts have been made to define the detonability (i.e., the potential to detonate under specific conditions) of a particular mixture. In particular, the minimum deposited energy required for direct initiation has been widely investigated. A first criterion was proposed by Zeldovich et al. [START_REF] Zeldovich | An experimental investigation of spherical detonation of gases[END_REF] suggesting the time taken by the decaying blast wave to reach the CJ velocity should be larger than the reaction time. This criterion introduces the concept of a critical radius, of the order of the detonation thickness, at which the blast wave should propagate faster than the CJ velocity to initiate a detonation. However, experimental tests showed later that the critical radius is actually three orders of magnitude larger than the detonation thickness [START_REF] Lee | Dynamic parameters of gaseous detonations[END_REF]. [START_REF] He | On the direct initiation of gaseous detonations by an energy source[END_REF] [START_REF] He | On the direct initiation of gaseous detonations by an energy source[END_REF] analysed the effect of the curvature in the quasisteady evolution of cylindrical and spherical detonations. Their study showed that flow divergence is responsible for the existence of a critical radius larger than the detonation thickness below which there is no possible solution for a self-propagating detonation. The critical radius predicted by this analysis was subsequently confirmed by Direct Numerical Simulations (DNS) [START_REF] He | Theoretical determination of the critical conditions for the direct initiation of detonations in hydrogen-oxygen mixtures[END_REF]. Nevertheless, this model fails to explain the non-steady behaviour also observed in DNS [START_REF] He | Theoretical determination of the critical conditions for the direct initiation of detonations in hydrogen-oxygen mixtures[END_REF][START_REF] Eckett | The role of unsteadiness in direct initiation of gaseous detonations[END_REF][START_REF] Ng | Direct initiation of detonation with a multi-step reaction scheme[END_REF]. The work of Liñán et al. [START_REF] Liñán | Initiation of reactive blast waves by external energy sources[END_REF] builds upon the quasi-steady approximation to study the external flow structure when transitioning from a strong non-reactive blast wave to a steady spherical CJ detonation. Recently, in a series of works by Clavin et al. [START_REF] Clavin | Analytical study of the direct initiation of gaseous detonations for small heat release[END_REF][START_REF] Clavin | Decay of plane detonation waves to the self-propagating Chapman-Jouguet regime[END_REF][START_REF] Clavin | Asymptotic analysis of the critical dynamics of spherical gaseous detonations[END_REF][START_REF] Tofaili | One-dimensional dynamics of gaseous detonations revisited[END_REF] the unsteady dynamics of detonations in the asymptotic limit of small heat release have been revisited. The asymptotic analysis including unsteadiness, curvature, and the burnt gas flow reduces the problem to a single Burgers' like equation with a reactive and a geometrical term. Numerical integration of this equation shows that this simplified model is able to reproduce to a certain extent the unsteady behaviour observed in DNS. Particularly, a deceleration of the front well below the CJ velocity and a later re-acceleration towards the self-sustained regime can be observed in the critical dynamics. Numerical results of a parametric study are presented here with special focus on the relative position of the sonic point as key driver for the critical slowdown prior to successful ignition.

Detonation model

The reactive Euler's equations in spherical geometry governing the dynamics of the inner detonation structure can be written as a set of hyperbolic equations

1 γp ∂ ∂t + (u ± a) ∂ ∂r p ± 1 a ∂ ∂t + (u ± a) ∂ ∂r u = q m c p T ω(T, Y ) t r -2 u r . (1)
∂ ∂t + u ∂ ∂r ln T - γ -1 γ ln p = q m c p T ω(T, Y ) t r , ∂ ∂t + u ∂ ∂r Y = ω(T, Y ) t r , (2)
where p, u, T , and Y are respectively the pressure, the radial velocity of the flow in the laboratory reference frame, the temperature, and the progress variable (Y = 0 in the initial mixture and Y = 1 in the burned gas), and a, γ, q m , t r , and ω are respectively the sound speed, the ratio of specific heat γ ≡ c p /c v , the chemical heat release per unit mass of mixture, the reaction time at the Neumann state of the planar CJ detonation and the non-dimensional heat-release rate. These equations relate the propagation of the disturbances of pressure p and radial velocity u to the rate of heat release ω/t r and the divergence of the flow 2u/r. Assuming an ideal gas model with p = c v (γ -1)ρT , a 2 = γp/ρ and γ = const., when the chemical kinetics ω(T, Y ) are known, the four equations in (1) and (2) form a closed set for p, u, T and Y .

Due to the large activation energy E that characterizes the combustion chemistry, the lead shock can be considered as a discontinuity in the study of the detonation structure. The boundary conditions in the initial mixture are thus given by the Rankine-Hugoniot relations with respect to the unburnt gas

p N p u = 1 + 2γ γ + 1 (M 2 -1), ρ u ρ N = 2 + (γ -1)M 2 (γ + 1)M 2 , u N a u = 1 - ρ u ρ N M, (3)
where the subscript u denotes the unburnt fresh mixture at rest and the subscript N denotes the Neumann state just behind the shock. The Rankine-Hugoniot relations depend exclusively on the propagation velocity of the lead shock D through its Mach number M ≡ D/a u .

The study of the direct initiation process begins with the blast wave assumed to be produced by an intense point explosion (i.e., a sudden release of a large amount of energy concentrated at a point).

The blast wave classical problem accepts a self-similar solution in the strong shock limit in which the velocity distribution evolves roughly in a linear fashion from the respective value behind the shock to the origin [START_REF] Clavin | Combustion Waves and Fronts in Flows: Flames, Shocks, Detonations, Ablation Fronts and Explosion of Stars[END_REF]. The initial flow field is then described by a linear profile from the Neumann state at the lead shock position to zero at the origin.

As in the detonation model of Clavin & Williams (2002) [START_REF] Clavin | Dynamics of planar gaseous detonations near Chapman-Jouguet conditions for small heat release[END_REF], the direct initiation is studied here through an asymptotic analysis in the limit of small heat release 2 ≡ q m /(c p T u ) 1 coupled with the Newtonian approximation (γ -1)/ 1. This distinguished limit introduces two main simplifications: the variation of the sound speed across the detonation structure, being of order 2 , can be neglected to first order; and, two time scales can be identified. The propagation velocities (in the lead shock reference frame) of the acoustic downstream running wave and the entropy wave are of the order of the sound speed which is greater by a factor 1/ than the propagation velocity of the upstream running acoustic wave. Therefore, to first order the effects of the downstream running wave and entropy wave are considered instantaneous, and the dynamics of the inner structure is governed by the larger time scale t r / . Hence, the inner structure of the unsteady and curved detonation must be studied in terms of the non-dimensional coordinate ξ attached to the moving front of the lead shock and the reduced time scale of order unity τ

ξ ≡ r -r f (t) a u t r , and τ ≡ t t r . (4)
In addition, it is convenient to introduce the dimensionless quantities of order unity in the small heat release limit as in [START_REF] Clavin | Decay of plane detonation waves to the self-propagating Chapman-Jouguet regime[END_REF] µ(ξ, τ

) ≡ 1 a u -(D o CJ -u) a u , π(ξ, τ ) ≡ 1 1 γ ln p p u , y(τ ) ≡ b D -D o CJ a u , (5)
where µ(ξ, τ ), π(ξ, τ ) and y(τ ) represent respectively the normalized flow velocity, the pressure and the propagation velocity of the lead shock. A reduced form of the activation energy b ≡ 2 (γ -1)E/(k B T u ) is also defined.

Neglecting terms smaller than 2 , the first order of the intrinsic dynamics of the detonation structure in the small heat release limit is governed by a single hyperbolic equation including reactive and geometric terms

∂µ ∂τ + µ - y b ∂µ ∂ξ = 1 2 ω(ξ, y) - 1 + µ rf (τ ) , (6)
where the reduced front radius rf (τ ) ≡ r f /(t r a u ) has been introduced. Similarly, the Neumann state condition (3) on the leading front is reduced to

ξ = 0 : µ = 1 + 2y(τ )/b. (7)
The chemical kinetics for a propagation velocity close enough to the CJ velocity y(τ ) = O(1) are modelled as in [START_REF] He | On the direct initiation of gaseous detonations by an energy source[END_REF]. The unsteady distribution of the rate of heat release ω(ξ, τ ) is expressed in terms of the front propagation velocity and the steady state distribution. Accordingly, the distribution of reaction rate is governed by a scaling law depending on the reaction zone length ξ b . Adopting the distribution of the steady planar CJ wave ω o CJ (ξ) presented in [START_REF] Clavin | Decay of plane detonation waves to the self-propagating Chapman-Jouguet regime[END_REF], the heat release distribution takes the form

ω(ξ, y) = e y(τ ) ω o CJ (ξe y(τ ) ), ξ < ξ b = -e -y : ω = 0. (8)
However, this model shares the limitation of the single step Arrhenius rate law. As pointed out in [START_REF] Ng | Direct initiation of detonation with a multi-step reaction scheme[END_REF], if the induction to reaction zone length ratio is fixed, the system will react until completion. It is required then to set a lower bound in the propagation velocity y below which the reaction remains frozen

y ≤ y c : ω(ξ, y) = 0. (9)
This lower bound corresponds to the well known cross-over temperature that leads to the chemicalkinetics quenching. A value of y c = -b/2 is used here considering typical values of cross-over temperature, see [START_REF] Clavin | Asymptotic analysis of the critical dynamics of spherical gaseous detonations[END_REF].

Numerical integration

In order to analyse the presented simplified model, the hyperbolic equation ( 6) has been integrated numerically. The numerical integration begins with a linear velocity profile as a rough approximation to the blast wave solution. Following a splitting strategy, the three terms of the hyperbolic equation are integrated successively. The convective component is evaluated using the High-Resolution Central Scheme of Kurganov & Tadmor (KT) [START_REF] Kurganov | New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection-Diffusion Equations[END_REF]. Next, both the curvature and the reactive terms are integrated by an explicit Euler method. At the end of the integration step, a new value of the propagation velocity y is computed from the solution obtained through the boundary condition [START_REF] Eckett | The role of unsteadiness in direct initiation of gaseous detonations[END_REF].

The KT central scheme admits a semi-discrete formulation maintaining its high-resolution independent of O(1/∆t) [START_REF] Kurganov | New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection-Diffusion Equations[END_REF]. Thanks to this feature, a sufficiently small time step can be chosen so that the truncation and splitting error do not degenerate the solution.

Several numerical integrations have been performed for a range of initial front radii rfi , and two reduced activation energies b illustrating a stable and a weakly unstable regime. For a fixed initial propagation velocity y i , each initial radius rfi represents an initial deposited energy. A greater initial energy would create in a larger blast wave propagating at equal velocity, thus these two parameters are directly linked. In the chemical kinetics model investigated here, the transition from marginally stable detonations to weakly unstable detonations occurs for a value of the reduced activation energy of b c ≈ 1.27 [START_REF] Clavin | Decay of plane detonation waves to the self-propagating Chapman-Jouguet regime[END_REF].

Results and discussion

The different trajectories 'detonation velocity y vs radius rf ' obtained are shown in Fig. 1. A stable detonation regime with b = 1 is presented on the left, while the evolution of a marginally unstable detonation b = 2 is shown on the right. The results of these stable and marginally unstable detonations are similar, except for the nonlinear oscillation superimposed on the trajectories. In both cases, successful initiation and detonation failure are observed revealing the dependence of the outcome on the initial conditions. The characteristic behaviour of the non-steady critical dynamics is also observed in both situations. After slowing down well below the CJ velocity, the detonation front re-accelerates back to the self-sustained regime. The evolution of the detonation initiation near criticality is closely related to the evolution of the sonic point ξ s , where the flow velocity relative to the front equals the sound speed. In terms of the problem variables, the sonic point is located where the relation µ(ξ s , τ ) = y(τ )/b is satisfied. At this point, the gas flow is isolated from the damping induced by the rarefaction wave. The sonic point trajectory is shown in Fig. 2. For clarity, the stable regime has been chosen to show the sonic point ξ s and the reaction wave exit ξ b positions relative to the lead shock. The difference between figures lies in the initial front radius rfi (i.e., the initial energy). On the left, the smaller initial radius leads to a detonation failure while on the right the detonation succeeds to initiate because of the larger initial energy. The corresponding trajectory 'detonation velocity y vs front radius rf ' is included in the background so that the front behaviour can be analysed in relation to the sonic point location. As the lead shock decays due to both the effect of rarefaction wave and flow divergence, the sonic point ξ s approaches the end of the reaction ξ b . Once the sonic point is close enough to the exit of the reaction wave, the time delay of the detonation inner structure response to the rarefaction wave becomes relevant. Initially, the decay rate of the front is reduced and a mechanism of slowdown can be identified. As long as the chemical-kinetics quenching does not terminate the combustion, the sonic condition will be attained, i.e. the sonic point will catch the end of the reaction layer. Note that the minimum propagation velocity y(τ ) for successful initiation corresponds to the sonic condition. As soon as the sonic point reaches the reaction wave, it isolates the detonation inner structure from the rarefaction wave. Then, the detonation now unaffected by the damping effect of the rarefaction wave, starts to accelerate towards the self-sustained CJ regime.

Conclusions

The asymptotic limit of small heat release reduces the mathematical description of the inner structure of a detonation to a single hyperbolic equation. This distinguished limit leverages the two-time-scale nature of the problem, which is amplified although it also characterises to a lesser extent the real detonations close to the CJ regime. The solution obtained by the numerical integration of the hyperbolic equation reveals the well known existence of a critical energy below which the direct initiation of a detonation fails. Furthermore, the usual critical dynamics, characterised by a decay well bellow the CJ velocity, followed by a re-acceleration towards the self-sustained regime, is also captured by this simple model. This behaviour is explained by the evolution of the sonic point where the velocity relative to the lead shock equals the sound speed. Initially, the propagation velocity descends below the CJ velocity while the sonic point approaches the exit of the reaction wave. As the sonic point gets closer to the end of the reaction zone, the slowdown mechanism of the wave decay is strengthened. The slowdown mechanism arises from an increase in the time delay of the response of the detonation inner structure to the trailing rarefaction wave. A minimum in the propagation velocity is observed when the sonic point reaches the inner structure of the detonation. The deceleration of the decay is followed by a re-acceleration back to the CJ regime in which the reaction zone is isolated from the rarefaction wave through the sonic point. Hence, the sonic point remains within the detonation inner structure during the re-acceleration phase and in the self-sustained regime. The initiation process will fail if the front velocity decays bellow the chemical-kinetics quenching point before the sonic point catches the detonation structure. At small front radius, the curvature is larger and therefore the flow divergence damping is stronger, so the outcome depends ultimately on the lead front radial position at which the velocity of the overdriven detonation crosses the CJ velocity.

Figure 1 :

 1 Figure 1: Trajectories 'detonation velocity y vs front radius rf ' for different initial radius. Left: Stable detonation b = 1. Right: Weakly unstable detonation b = 2.

Figure 2 :

 2 Figure 2: Sonic point ξ s and end of reaction ξ b trajectories in thick line with its corresponding 'detonation velocity y vs front radius rf ' trajectory in the background for a stable detonation b = 1. Left: Failure of detonation initiation for rfi = 0.7. Right: Successful initiation of the detonation for rfi = 1.0.

th ICDERS -June 19 -24, 2022 -Napoli, Italy

Acknowledgements

We thank Prof. Amable Liñán for stimulating discussions. This work was supported by the French National Research Agency under agreement ANR-18-CE05-0030.