Kevin M Hickson
email: km.hickson@ism.u-bordeaux1.fr

Yury V Suleimanov

An experimental and theoretical investigation of the C( 1 D) + D 2 reaction

In a previous joint experimental and theoretical study of the barrierless chemical reaction C( 1 D) + H 2 at low temperatures (300-50 K) [K.

Introduction

Gas-phase reactions involving neutral carbon atoms are considered to be important in combustion systems, in the chemistry of planetary atmospheres and in the interstellar medium.

Many experimental and theoretical studies [1][2][3][4][5][6][7][8][9][10][11][12][13] of the reactivity of ground state C( 3 P) atoms have been performed, addressing their kinetic and dynamic aspects over a wide range of temperatures 2-6, 8, 9 and collision energies [9][10][11][12][13] and for a variety of different coreagent molecules. The reactivity of the excited 1 D state of atomic carbon is generally less well understood, with the exception of the C( 1 D) + H 2 reaction; a process which is simple enough to allow comparison between precise state-selected experimental measurements [14][15][16] , quasi classical trajectory calculations [17][18][19] and detailed quantum mechanical methods, [20][21][22][23][24][25][26][27][28][29] even for several isotopic forms of the H 2 molecule (H 2 , HD, D 2 ) 18,19,[30][31][32] . Although there are several experimental works following the dynamics of other C( 1 D) atom reactions under a range of conditions 13,[33][34][35][36] , few kinetic investigations are reported in the scientific literature. Of these, the majority is confined to room temperature [37][38][39][40][41][42] . Indeed, the only temperature dependent measurements of the rate constants for C( 1 D) reactions to have been performed to date are those between C( 1 D) + CH 3 OH 7 , C( 1 D) + H 2 43 and the nonreactive quenching of C( 1 D) atoms with N 2 44 . Extending such studies to include other isotopic forms of the same molecule would allow us to look for interesting reactivity differences, thereby providing a rigorous test of the underlying potential energy surfaces (PESs).

Previous room temperature kinetic studies of the C( 1 D) + H 2 , HD and D 2 reactions 42 demonstrated a noticeable isotope dependence of the rates of these processes, with the rate constant decreasing from H 2 to HD to D 2 . Lin and Guo 24 were able to reproduce the experimental results of Sato et al. 42 using a wavepacket based statistical model over the lowest 1 Aʹ PES, attributing the reactivity difference to a kinematic effect related to the increasing reduced mass. The recent joint experimental and theoretical investigation of the C( 1 D) + H 2 reaction by Hickson et al. 43 has cast some doubt on the reliability of these earlier experiments, measuring rate constants which are significantly larger than those measured by Sato et al. 42 . On the theoretical side, novel rate theory based on ring polymer molecular dynamics (RPMD) 45 was implemented that demonstrated very reliable and predictable behavior in multifarious prototype cases as discussed in a recent review 46 . Using the RPMD method, these authors 43 were also able to demonstrate that the new experimental results could be well described if reaction occurred adiabatically over both the lowest 1 Aʹ and first excited 17,20 using a quantum mechanical time dependent wavepacket method, also considering the permutation-inversion symmetry and the Coriolis interactions 22,32 . They found that while the upper 1 Aʹʹ plays an important role in all these reactions, the reactive contribution of this surface is particularly large in the case of C( 1 D) + D 2 . Clearly, an investigation of the temperature dependent kinetics of the reactions of C( 1 D) atoms with other deuterated forms of hydrogen would provide an ideal opportunity to validate the earlier hypotheses of both Defazio et al. 32 and Hickson et al. 43 with regard to the role of the 1 Aʹʹ PES in these systems.

Here, we report the results of an experimental and theoretical investigation of the C( 1 D) + D 2 reaction at low temperature. A supersonic flow reactor employing pulsed laser photolysis for C( 1 D) generation coupled with vacuum ultraviolet laser induced fluorescence detection of the deuterium atom products allows us to follow the reaction kinetics over the 50-296 K range. In parallel, RPMD simulations were carried out for the title reaction in the same temperature range.

Sections 2 and 3 describe respectively the experimental and theoretical methods used in this work and the results which are further summarized in section 4.

Experimental Results

Measurements were performed using a continuous flow reactor, employing axisymmetric Laval nozzles to generate the cold supersonic flows. As this method is well described elsewhere 47, 48 , only features specific to the current work will be presented here. For the present investigation, it was impossible to use any of our Laval nozzles based on N 2 as the carrier gas, as the C( 1 D) + N 2 quenching reaction has been shown to become more efficient as the temperature falls 44 . As a result, three nozzles employing the carrier gas argon were used, allowing temperatures of 50 K, 75 K and 127 K to be attained. Details of the flow characteristics of these three nozzles can be found in Table 1 of Grondin et al. 49 . By removing the Laval nozzle and by reducing the flow velocity, it was also possible to perform room temperature measurements (296 K) with argon as the carrier gas.

CBr 4 was used as the source of C( 1 D) atoms in these experiments. These precursor molecules were carried into the reactor by a small argon flow passing over solid CBr 4 held in a separate vessel at a known fixed pressure and temperature. This allowed us to estimate the concentration of CBr 4 molecules within the supersonic flow to be less than 2 × 10 13 molecule cm -3 . The pulsed multiphoton dissociation of CBr 4 at 266 nm produced both C( 1 D) and C( 3 P) atoms along the entire length of the supersonic flow [4][5][6][7][8]44 . Although the major photolysis product is C( 3 P) atoms 7 , the C( 3 P) + D 2 → CD + D reaction can be neglected at room temperature and below due to its high endothermicity of 120 kJ mol -1 . Due to inefficient gas-phase spin conversion, the D 2 used in the present experiments was characterized by a fixed ortho/para ratio of 2:1 at all temperatures.

This can be contrasted with the expected value at 50 K, which should approach 4:1 assuming thermal equilibrium is attained. The D 2 concentration was held in large excess with respect to C( 1 D) atoms in this work, so that pseudo-first-order conditions could be assumed. As C( 1 D) atoms could not be followed directly in the present experiments, product D( 2 S) atoms were detected instead through vacuum ultraviolet laser induced fluorescence (VUV LIF) at 121.534

nm. Light at this wavelength was generated by frequency doubling the 729.2 nm output of a pulsed dye laser, yielding a UV beam at 364.6 nm, which was then focused into a cell containing 210 Torr of Kr and 540 Torr of Ar, producing the required tunable VUV light by third harmonic generation. A magnesium fluoride lens positioned at the exit of the tripling cell allowed the VUV beam to be collimated and steered into the reactor where it was allowed to cross the supersonic flow. The on-resonance VUV emission from excited D-atoms within the flow were collected by a lithium fluoride lens and focused onto the photocathode of a solar blind photomultiplier tube.

None of the gases used in the experiments (Ar 99.999%, D 2 99.8%, Kr 99.99%) were purified prior to use. Instead, these gases were flowed directly from cylinders into mass-flow controllers which provided precise control over the carrier gas, reagent gas and precursor concentrations within the supersonic flow. All mass-flow controllers were calibrated using the pressure rise at constant volume method for each individual gas used.

Typical product D( 2 S) VUV LIF intensity profiles, recorded as a function of delay time between photolysis and probe lasers at 50 K are displayed in Figure 1.

These curves were readily described by a functional fit of the form

𝐼 $ = 𝐴 𝑒𝑥𝑝 -𝑘 -($) 𝑡 -𝑒𝑥𝑝 -𝑘′𝑡 (1)
where Decay curves were recorded at several D 2 concentrations for any particular experiment, to yield a wide range of values for the first-order rate constant 𝑘′. These 𝑘′ values were then plotted against the corresponding D 2 concentration to yield the second-order rate constant at a given temperature from a weighted linear least squares fit to the data (weighting was performed using the statistical uncertainties generated by the biexponential fitting procedure outlined above).

𝑘′ = k[D 2 ] + k L(C) ,
Examples of such plots obtained at 50 K and at 296 K are shown in Figure 2.

The measured rate constants for the C( 1 D) + D 2 reaction are summarized in Table 1 and displayed as a function of temperature in Figure 3 alongside the present RPMD rate constants and the previous results, both theory and experiment for this process. We also include our previous experimental results for the C( 1 D) + H 2 reaction 43 showed a slight positive temperature dependence 43 . Although these effects are relatively minor, the differences could be due to the use of non-equilibrated mixtures of ortho/para-H 2 ,-D 2 in both studies. If the reactivity of the ortho and para forms of these isotopologues towards C( 1 D) is significantly different, this should induce the largest impact on the measured rate at the lowest temperatures where deviations from equilibrium are greatest. In the case of the C( 1 D) + H 2 reaction, for such an effect to occur would require the j = 0 level of H 2 (para-H 2 ) to react more rapidly than the j = 1 level (ortho-H 2 ) at the lowest temperatures. Consequently, the measured rates would be lower than expected due to an excess of ortho-H 2 in these earlier experiments.

Although an earlier time-dependent wavepacket study of the C( 1 D) + H 2 reaction over uncoupled

1 Aʹ and 1 Aʹʹ PESs 22 concluded that the rotational state specific rate constant for the j = 0 level was approximately 15 % higher than the corresponding j = 1 level at 300 K, the rate constant for the j = 0 level is seen to increase at 100 K whereas the value for the j = 1 level remains constant.

As a result, it may be possible that the relative reactivity trends are inverted at even lower temperatures. In the case of the C( 1 D) + D 2 reaction, a time-dependent wavepacket study indicated that the j = 0 (ortho) and j = 1 (para) levels of D 2 displayed the same reactivity at 300 K 32 . Unfortunately, no further information was provided with regard to the temperature dependence of the state specific rates for this system, making it impossible to determine whether the present measurements deviate from the expected values at equilibrium. Nevertheless, the temperature dependence of the RPMD results (which assume Boltzmann statistics for the rotational levels of D 2 ) is in good agreement with the experimentally determined one, suggesting that the relative reactivity of j = 0 and j = 1 levels of D 2 could be similar below 100 K. We will continue the discussion of the results presented in Figure 3 in Section 4.

Theoretical Results

The present theoretical simulations of the title reaction were performed using the RPMD method. This recently proposed 45 semiclassical method takes into account quantum mechanical effects of nuclear motions via the so-called "classical isomorphism" between the quantum mechanical formalism of statistical mechanics and the classical one for the ring polymer with beads being classical copies of the original system connected with neighbors in the necklace via fictitious temperature-dependent harmonic interaction. The classical isomorphism is exact for various static properties but it was also shown to provide a reliable approximation to real-time Kubo-transformed correlation functions responsible for various dynamical properties 50 . It has been found that RPMD provides very reliable estimates for real-time correlation functions used to calculate thermal rate constants 46,51,52,53 . Using simple gas-phase chemical reactions 43,[54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69][70] , the method was thoroughly compared with rigorous quantum mechanical results (when possible), other dynamics approximations and experiment and its accuracy and reliability has been confirmed as described in a recent review by one of us (YVS). 46 In particular, RPMD exhibits very reliable and predictable behavior for barrierless chemical reactions with deviation from the rigorous quantum dynamics results within the convergence error of the computational procedure (typically below 15 %) 60,61 . As mentioned in the introduction, RPMD demonstrated an excellent agreement with experiment in our previous calculations for the C( 1 D) + H 2 reaction in the temperature range 50-300 K 43 .

In the present RPMD calculations we used the same computational strategy as in our previous study of the C( 1 D) + H 2 reaction. All calculations were performed using the RPMDrate code 71 .

The parameters of the simulations are summarized in Table S1. They are similar to the previous parameters for the C( 1 D) + H 2 system, except that fewer beads were required to converge the RPMD rate constants for the title D-transfer reaction. As previously 43 , we used the multireference configuration-interaction PESs of Bañares et al. 17 for the ground state ( 1 Aʹ) and of Bussery-Honvault et al. 21 for the first excited state ( 1 Aʹʹ). Both PESs were obtained from ab initio data and are of insertion type with deep wells (4.32 and 3.46 eV relative to the reactants for 1 Aʹ and 1 Aʹʹ, respectively) and barrierless minimum energy paths in the entrance channel (perpendicularly constrained approach and bent approach of C towards D 2 around 60° for 1 Aʹ and 1 Aʹʹ, respectively). As for the C( 1 D) + H 2 reaction, we used the 1 Aʹ and 1 Aʹʹ saddle points to initiate the RPMDrate calculations, though we note that the RPMD rate is rigorously independent of our choice of the initial transition state structure 52 .

The results of the RPMDrate simulations are summarized in Figure 4 and 5 and in Table 2.

Figure 4 shows that for both PESs, ring polymer free energy profiles demonstrate small thermodynamically induced barriers near the entrance to the well which diminish as we decrease the temperature. The 1 Aʹʹ surface exhibits a more pronounced decrease. Figure 5 shows that for the 1 Aʹ PES, recrossing dynamics are significantly enhanced as the temperature is decreased.

Though the free energy barrier decreases with the temperature, this does not compensate the decrease due to recrossing dynamics. As a result, the RPMD rate coefficient for the 1 Aʹ PES decreases with decreasing temperature, see Table 2. For the 1 Aʹʹ PES, the temperature dependences of the ring polymer free energy barriers and the transmission coefficient are qualitatively the same, though the decrease in the plateau values of the ring polymer transmission coefficients is less pronounced and, as a result, the RPMD rate constant increases with the temperature. As a result, the rate constants for the 1 Aʹ PES are higher than those for the 1 Aʹʹ one only at 300 K. At temperatures below 75 K, the rate constants for the 1 Aʹʹ PES are twice as large as those for the 1 Aʹ PES.

When compared to the previous theoretical results for the C( 1 D) + H 2 chemical reaction 43 , we note that the above temperature tendencies are similar to the previously observed ones. The final RPMD rate constants are smaller for C( 1 D) + D 2 and, as expected this difference decreases with decreasing temperature (see also Figure 3). However, the kinetic isotope effect depends on the PES: the rate constants are smaller for the C( 1 D) + D 2 reaction compared to the C( 1 D) + H 2 reaction by 20-40 % for the 1 Aʹ PES and by only 13-16% for the 1 Aʹʹ PES. The present results

show that the 1 Aʹʹ PES contribution is more important for the C( 1 D) + D 2 reaction, in line with the previous quantum mechanical time dependent wavepacket results 22,32 . Clearly, the chemical dynamics of the title reaction exhibit high sensitivity to the underlying PESs which therefore require very accurate treatment.

Discussion and Conclusions

Following our previous work on the C( 1 D) + H 2 chemical reaction 43 , in this paper we performed a joint experimental and theoretical study of its deuterated counterpart, C( 1 D) + D 2 , using the same experimental and theoretical techniques. Measurements were made using a continuous supersonic flow apparatus over the 50 -296 K range. Due to the lack of an appropriate method for the direct detection of C( 1 D) atoms, rate constants were derived instead from product D( 2 S) atom formation curves. Theoretical results were obtained using novel rate theory based on ring polymer molecular dynamics (RPMD) which has deserved high confidence in the previous intensive examination of the method using prototype chemical reactions 43 ,

including those with profiles of potential energy surface similar to the title reaction.

Figure 3 shows that the present RPMD and experimental results are in good agreement, with the former being higher than the latter by only 15-25 % but correctly reproducing the very slight temperature dependence observed experimentally. Such deviations are slightly higher than the convergence error of the computational procedure 43,71 or the typical error of RPMD observed for prototype insertion chemical reactions 43,60,61 . We note that inaccuracy in the present RPMD results can be attributed to the sharp probability resonances found for the C( 1 D) + D 2 reaction 32 and which are not taken into account by the RPMD formalism 43 . It can also be attributed to inaccuracies of the underlying potential energy surfaces (PESs) and the adiabatic limit used in the present description of the electronic structure (i.e., the two PESs were treated as being uncoupled), as was discussed in our previous study of C( 1 D) + H 2 43 . However, we note that the previous attempts to include non-adiabatic effects in calculations of the C( 1 D) + H 2 reaction led to a strong reduction in the overall rate constant and an inversion of the temperature dependence compared to the adiabatic limit 43 . Nonadiabatic RPMD studies for the C( 1 D) + X 2 reactions with X = H or D are clearly desirable in the future, when the nonadiabatic RPMD rate theory 72,73 will become as robust and reliable as its original counterpart for single PESs.

Despite the small discrepancy, the present results demonstrate high consistency between theory and experiment. They confirm that both of the two lowest single PESs actively participate in the title reaction. In fact, the contribution from the first excited surface 1 Aʹʹ prevails at temperatures below 300 K, being twice as large as the contribution from the 1 Aʹʹ surface at 75-50 K. We note that the previous quantum mechanical time dependent wavepacket calculations, also obtained using two uncoupled PESs, 32 exhibit similar accuracy though underestimating experimental rates, as can be seen from Figure 3. Overall, the present RPMD results confirm the validity of this method for low-temperature studies of complex-forming reactions. T / K

  where k represents the second order rate constant for the C( 1 D) + D 2 reaction, k L(C) represents the first-order loss of C( 1 D) by other processes such as secondary reactions and diffusion, [D 2 ] is the D 2 concentration, t is time and A is the signal amplitude. The first term in expression (1) represents D-atom losses with a first-order rate constant k L(D) .

  for comparison. Error bars on the present experimental values were derived by combining the statistical uncertainties obtained from second-order fits such as those shown in Figure 2 with an estimated systematic error of 10 % of the nominal rate constant value. This systematic error was considered to potentially arise from calibration errors in various pieces of equipment such as mass-flow controllers and pressure gauges in addition to possible errors in the flow density and velocity calculations. While the measured rate constants for the C( 1 D) + D 2 reaction display only a very slight negative temperature dependence (independent of temperature considering the experimental error bars), it is interesting to note that our earlier measurements of the C( 1 D) + H 2 reaction

Figure 1

 1 Figure 1 Typical D( 2 S) atom formation curves from the C( 1 D) + D 2 → CD + D reaction recorded at 50 K. (Blue solid squares) [D 2 ] = 6.8 × 10 14 molecule cm -3 ; (red solid circles) [D 2 ] = 6.0 × 10 13molecule cm -3 . Solid lines represent biexponential fits to the data of the form given by equation

Figure 2 Figure 3 2 Figure 4

 2324 Figure 2 Pseudo-first order rate constants for the C( 1 D) + D 2 reaction as a function of [D 2 ],

2 Figure 5

 25 Figure 5 Ring polymer transmission coefficients for the C( 1 D) + D 2 reaction at 50, 75, 127 and 300 K over two potential energy surfaces 1 Aʹ (lower panel) and 1 Aʹʹ (upper panel).

  1

Table 1

 1 Measured rate constants for the C( 1 D) + D 2 reaction. Uncertainties on the calculated temperatures represent the statistical (1σ) errors obtained from Pitot tube measurements of the impact pressure; b Number of individual measurements; c Uncertainties on the measured rate constants represent the combined statistical and systematic errors as explained in the text.

	T / K	N b	[D 2 ] / 10 14 molecule cm -3	k C(1D)+D2 / cm 3 molecule -1 s -1
	296	40	0.5 -12.3	(1.9 ± 0.2) c ×⋅10 -10
	127 ± 2 a	40	0.4 -9.9	(2.0 ± 0.2) ×⋅10 -10
	75 ± 2	41	0.4 -5.7	(2.0 ± 0.2) ×⋅10 -10
	50 ± 1	35	0.6 -7.3	(2.1 ± 0.2) ×⋅10 -10
	a

Table 2

 2 Summary of the RPMDrate calculations for the C( 1 D) + D 2 reaction over the 1 Aʹ and 1 Aʹʹ PESs at 50, 75, 127 and 300 K: k QTST -centroid-density quantum transition state theory rate coefficient; κ plateau -ring polymer transmission coefficient; k RPMD -RPMD rate coefficient. The final RPMD rate constants (last column) are corrected by the electronic partition function Q el = 1/5. The parentheses denote powers of ten.

	T(K)	k QTST	κ plateau	k RPMD	k RPMD (corrected by Q el )
		PES 1 Aʹ
	300 1.20(-09) 0.489 5.84(-10)	1.17(-10)
	127 9.43(-10) 0.529 4.99(-10)	9.98(-11)
	75 7.48(-10) 0.594 4.44(-10)	8.88(-11)
	50 6.16(-10) 0.652 4.02(-10)	8.04(-11)
		PES 1 Aʹʹ
	300 2.08(-09) 0.247 5.14(-10)	1.03(-10)
	127 2.51(-09) 0.277 6.95(-10)	1.39(-10)
	75 2.51(-09) 0.323 8.11(-10)	1.62(-10)
	50 2.47(-09) 0.358 8.85(-10)	1.77(-10)

Acknowledgments

KMH acknowledges support from the French INSU/CNRS Programs 'Physique et Chimie du Milieu Interstellaire' (PCMI) and 'Programme National de Planétologie' (PNP). YVS thanks the European Regional Development Fund and the Republic of Cyprus for support through the Research Promotion Foundation (Project Cy-Tera ΝΕΑ ΥΠΟΔΟΜΗ/ΣΤΡΑΤΗ/0308/31).


Loic Brancheriau
email: loic.brancheriau@cirad.fr

Sylvie Sabatier

Marc Jaeger

Nicolas Hemery

Christ Mougani

Philippe De Reffye

Ali Abed Alsater

Serge Palu

Modelling the Guayule plant growth and development with a Functional Structural Plant Model

Keywords: FSPM, Guayule, Natural rubber, GreenLab model, Organ sinks

Introduction

The Guayule (Parthenium argentatum, Asteraceae), is a branched shrub native to northern Mexico showing a growing interest in research and agriculture [START_REF] Ray | Guayule: A source of natural rubber[END_REF]. However, few studies have been done on the structure or functioning of the species and the production itineraries in relation to natural rubber production are still to assess. This study aims to propose a first FSPM of the species using the GreenLab model, calibrated with data issued from two varieties in different environmental conditions. The objective is to understand the interaction between the development of Guayule and natural rubber production.

Materials and Methods

Experimental plots were conducted on two varieties (Fig. A), named CL1 (parent USDA 11591 triploid) and CLA1 (parent USDA AZ 101 tetraploid). The CL1 variety has small stem and leaf size compared to the CLA1 variety. These varieties were planted in the field under various environmental conditions related to density (9091 and 62500 plants/ha), and hydric pressure (100%, 66% and 33% of optimal supply). The measurements carried out to count number of phytomers per module as well as the number of relay axes. We also clustered the total dry weight of leaves and stems per module in order to build the axis organic series [START_REF] Buis | Relations dimensionnelles dans une série organique en croissance chez une plante supérieure[END_REF] and to estimate the content of natural rubber from a NIRS analysis [START_REF] Taurines | Determination of natural rubber and resin content of guayule fresh biomass by near infrared spectroscopy[END_REF]. The field measurements from the series constituted a target to calibrate the GreenLab model functional parameters such as the organ sink functions [START_REF] Kang | Estimating Sink Parameters of Stochastic Functional-Structural Plant Models Using Organic Series-Continuous and Rhythmic Development[END_REF]. We chose such a cohort model, since the structure (here the module rank) can be analysed in the biomass production without detailed explicit geometrical computations.

Results and Discussion

Experimental results. The Guayule shows a sympodial development composed of modules with a terminal inflorescence. Its architecture corresponds to the Leeuwenberg's model (F. Hallé et al. 1978) producing eight to nine modules per year. The number of internodes within the modules is almost stable regardless of variety and environmental conditions. The architecture of the two varieties CL1 and CLA1 is very similar, including the duration of the functioning of the organs (leaves and inflorescences). In the experiment, the biomass production per module (and thus per plant) was shown to be dependent on both density and water supply, according to the variety. Natural rubber content was found higher on young modules, and the water regime affected the natural rubber content of modules according to the variety (Fig. B).

Modelling results. We have modelled the number of internodes per module with a constant binomial law. The number of relay axis per module was modelled with a damped binomial law, simulating the branching rate decreases on young modules. These laws allowed simulating the plant structure (Fig. C). We then retrieved the organ sinks parameters, fitting the target file, and found significant differences in the leaf sink functions, stronger for the variety CLA1. As a result, the simulation of the functional model allowed the module diameters and the total plant height of the two varieties to be recovered correctly (Fig. D.). We still need to validate the relation between the internode growth and natural rubber production from new field analysis.

Acknowledgement

Authors thank the Agropolis Fondation supporting the "modelling guayule growth as an alternative source of natural rubber" project, referenced as PAI-10605-026. Figure: FSPM Guayule model. A: Field view with CL1 & CLA1 varieties. B: Measured values of the natural rubber content (%) in the optimal (R1) and restricted (R3) irrigation scenarios. C: Simulation of the CL1 structure model. D: Functional structural simulation of primary and secondary module growth on both varieties.


Benjamin Bocquillon

Philippe Feyel

Guillaume Sandou

Pedro Rodriguez-Ayerbe

Computation of Neural Networks Lyapunov Functions for Discrete and Continuous Time Systems with Domain Of Attraction Maximization

Keywords: Lyapunov function, Domain of attraction, Optimization, Neural network, Nonlinear system

This contribution deals with a new approach for computing Lyapunov functions represented by neural networks for nonlinear discrete-time systems to prove asymptotic stability. Based on the Lyapunov theory and the notion of domain of attraction, the proposed approach deals with an optimization method for determining a Lyapunov function modeled by a neural network while maximizing the domain of attraction. Several simulation examples are presented to illustrate the potential of the proposed method.

INTRODUCTION

Lyapunov theory, introduced in the late nineteenth century [START_REF] Lyapunov | The general problem of the stability of motion[END_REF], is a classical way to investigate for the stability of an equilibrium point for a dynamical system. The method relies on the search for a function that exhibits three important properties that are sufficient for establishing the Domain Of Attraction (DOA) of a stable equilibrium point : (1) it must be a local positive definite function; (2) it must have continuous partial derivatives, and (3) its time derivative along any state trajectory must be negative semi-definite. Although efficient to prove stability once the so-called Lyapunov function is known, there is no general method for constructing such a function.

The Lyapunov function construction is still an open problem, but several methods, often based on optimization, have emerged in the literature. One can cite [START_REF] Panikhom | Numerical approach to construction of Lyapunov function for nonlinear stability analysis[END_REF], where the best quadratic Lyapunov function is looked for. However, these methods are too conservative in case of industrial complex systems. The work of [START_REF] Argáez | Iterative construction of complete Lyapunov functions[END_REF] proposes a new iterative algorithm that aims to avoid obtaining trivial solutions when constructing completive Lyapunov functions. This algorithm is based on mesh-free numerical approximation and analyses the failure of convergence in certain areas to determine the chain-recurrent set. Once again, the method appears too difficult for being used in an industrial context where flexibility is needed. Finally, the survey [START_REF] Giesl | Review on computational methods for Lyapunov functions[END_REF] has brought different methods and gave a wide overview of the methods that can be used for the Lyapunov function computation. It proposes conservative methods when the system is complex and highly non-linear.

However, to the authors' mind, Artificial Intelligence, Machine Learning and Neural Network bring a great opportunity to design powerful tools to justify and quickly certificate complex industrial systems such as in the aerospace field for instance. One of the first papers using Artificial Intelligence to compute Lyapunov function is [START_REF] Prokhorov | A Lyapunov machine for stability analysis of nonlinear systems[END_REF], where a so called Lyapunov Machine, which is a special-design artificial neural network, is described for Lyapunov function approximation. The author indicates that the proposed algorithm, the Lyapunov Machine, has substantial computational complexity among other issues to be resolved and defers their resolution to future work. The work of [START_REF] Banks | Searching for Lyapunov functions using genetic programming[END_REF] suggests a Genetic Programming for computing Lyapunov functions. However, the computed Lyapunov functions may have locally a conservative behavior. In this study, the use of Neural Networks allows to overcome these limits. Neural Networks are known to be powerful regressors that can approximate any nonlinear function. As a result, they appear as a good candidate for the construction of a Lyapunov function. In the literature, one can find other works using neural network to construct or approximate a Lyapunov function [START_REF] Serpen | Empirical approximation for Lyapunov functions with artificial neural nets[END_REF] and the paper [START_REF] Petridis | Construction of neural network based Lyapunov functions[END_REF] where the authors propose an interesting and promising approach for the construction of Lyapunov functions represented by neural networks.

In [START_REF] Bocquillon | Efficient construction of neural networks Lyapunov functions with domain of attraction maximization[END_REF], the authors propose to use a new constrained optimization scheme such that the weights of this neural network are calculated in a way that is mathematically proven to result in a Lyapunov function while maximizing the DOA.

Although most of works relative to Lyapunov computation deal with continuous time systems, only few ones related to the discrete time case. Numerical methods to compute Lyapunov functions for nonlinear discrete time systems have been presented, for example in [START_REF] Giesl | On the determination of the basin of attraction of discrete dynamical systems[END_REF], where collocation is used to solve numerically a discrete analogue to Zubov's partial differential equation using radial basis functions. For nonlinear systems with a certain structure there are many more approaches in the literature. To name one, in [START_REF] Milani | Piecewise-affine Lyapunov functions for discrete-time linear systems with saturating controls[END_REF] the parameterization of piecewise-affine Lyapunov functions for linear discrete systems with saturating controls is discussed. However, although effective, these approaches either look difficult to implement for real complex systems or are too conservative.

In this paper, we extend the proposition in [START_REF] Bocquillon | Efficient construction of neural networks Lyapunov functions with domain of attraction maximization[END_REF] dedicated to continuous systems to the case of discrete-time ones. We enhance the usability of the algorithm by extending this approach for discrete time systems. Our contribution is to unify the continuous and discrete time cases towards an efficient formulation of the optimization problem. Note that the goal of this study is not to stabilize a given plant, but to analyze its stability, computing an approximation of the potential domain of attraction. Our paper does not deal of the stabilization of the system itself.

The paper is composed of 4 sections: the first part is related to the design of neural network Lyapunov function for the continuous part. The second part, which is the main contribution of this paper, proposes an extension to the discrete time case. The third part is where we present the proposed algorithm to compute a candidate Lyapunov function through optimal neural network weights calculation while maximizing the DOA. And the last one describes examples showing efficiency of the paper.

Continuous Time Case

Theoretical Background

Notations and definitions used in the paper are called in this section. Let R denote the set of real numbers, R + denote the set of positive real numbers, . denote a norm on R n , and X ⊂ R n , be a set containing X = 0. Consider the autonomous system given by (1).

Ẋ = f (X) (1)
where f : X → R n is a locally Lipschitz map from a domain X ⊂ R n into R n and there is at least one equilibrium point X e , that is f(X e ) = 0.

Theorem 2.1 (Lyapunov Theory) [START_REF] Khalil | Nonlinear systems[END_REF]. Let X e = 0 be an equilibrium point for (1). Let V : D → R be a continuously differentiable function,

V (0) = 0 and V (X) > 0 in D -{0} (2) V (X) ≤ 0 in D (3) then, X e = 0 is stable. Moreover, if V (X) < 0 in D -{0} (4)
then X e = 0 is asymptotically stable.

Where D ⊂ X ⊂ R n is called Domain Of Attraction (DOA) and the system will converge to 0 from every initial point X 0 belonging to D.

Neural Network Formalism

The work in this section is based on the approach developed in the paper [START_REF] Bocquillon | Efficient construction of neural networks Lyapunov functions with domain of attraction maximization[END_REF] and is reported here for sake of clarity.

Let us consider the autonomous system in (1), in which we assume without any loss of generality, that the equilibrium point considered for the stability analysis is the point 0 (X=0). Therefore,

f (0) = 0 (5)
Suppose V(X) is a scalar, continuous and differentiable function and its derivative with respect to time is given in (6).

G(X) = dV dt = n ∑ j=1 ∂V ∂x j f j (X) (6) with X = [x 1 ,...,x n ] T .
The Lyapunov function is modeled by a neural network with a single hidden layer whose weights are calculated in such a way that is proven mathematically that the resulting neural network implements indeed a Lyapunov function, showing the asymptotic stability in the neighborhood of 0. We assume that the Lyapunov function V(X) is represented by a neural network where the x i are the inputs, w ji are the weights of the hidden layer, a i the weights of the output layer, h i are the biases of the hidden layer, θ is the bias of the output layer; i=1,...,n and j=1,...,K where K is the number of neurons of the hidden layer and σ is the activation function of the neural network.

Therefore, V(X) can be expressed as:

V (X) = K ∑ i=1 a i σ(ν i ) + θ (7) ν i = n ∑ j=1 w ji x j + h i (8)
From ( 2) and ( 4), sufficient conditions for the point 0 of system (1) to be asymptotic stable in the sense of Lyapunov are:

(a1) V(0) = 0. (a2) V(X) > 0 for all nonzero X in a neighbour- hood of 0. (a3) G(0) = 0.
(a4) G(X) < 0 for all nonzero X in a neighbourhood of 0.

From (a1) and (a2), sufficient conditions for V(X) to have a local minimum at 0 are [START_REF] Petridis | Construction of neural network based Lyapunov functions[END_REF]:

(v1) V(0) = 0. (v2)
∂V ∂x j X=0 = 0 for all j=1,2,...,n.

(v3) H V (the matrix of 2 nd derivatives of V at X=0) is positive definite.

In the same way from (a3) and (a4), sufficient conditions for G(X) to have a local maximum at 0 are [START_REF] Petridis | Construction of neural network based Lyapunov functions[END_REF]:

(d1) G(0) = 0. (d2) ∂G ∂x j X=0
= 0 for all j=1,2,...,n.

(d3) H G (the matrix of 2 nd derivatives of G at X=0) is negative definite.

Then, the second derivative of V(X) and G(X) are computed as functions of the neural network :

V qr = ∂ 2 V ∂x q ∂x r X=0 = K ∑ i=1 a i d 2 σ(ν i ) dν 2 i X=0 ∂v i ∂x r X=0 w qi = K ∑ i=1 a i d 2 σ(ν i ) dν 2 i X=0 w ri w qi (9) G l p = ∂ 2 G ∂x l ∂x p X=0 = n ∑ j=1 K ∑ i=1 a i d 2 σ(ν i ) dν 2 i X=0 w ji w li J j p + + n ∑ j=1 K ∑ i=1 a i d 2 σ(ν i ) dν 2 i X=0
w ji w pi J jl (10) where J qr = ∂ f q ∂x r X=0 q=1,...,n; r=1,...,n; l=1,...,n; p=1,...,n.

Therefore,

H V = [V qr (X = 0)] V qr is given by (9) (11)
H G = [G l p (X = 0)] G l p is given by (10) (12)
Assuming that the Lyapunov function is represented by a neural network, conditions (v1) -( v3) and (d1) -(d3) reduce to (we choose here σ(ν) = tanh(ν)):

(t1) K ∑ i=1 a i σ(h i )+θ = 0. ( 13
) (t2) K ∑ i=1 a i (1 -tanh 2 (h i ))w qi = 0 for q=1,...,n. (14)
(t3) H V as given by ( 11) is positive definite.

(t4) H G as given by ( 12) is negative definite.

We only deal with differentiable activation functions.

Discrete Time Case

Additional Theoretical Background

Consider the autonomous discrete time system given by (15).

X(k + 1) = f (X(k)) (15)
where f : X → R n is a locally Lipschitz map from a domain X ⊂ R n into R n and there is at least one equilibrium point X e , that is f(X e ) = X e .

Theorem 2.2 (Lyapunov Theory Discrete Time) [START_REF] Khalil | Nonlinear systems[END_REF].

Consider the autonomous system given by ( 15).

Let X e = 0 be an equilibrium point for (15). Let V : D → R be a continuous function,

V (0) = 0 and V (X(k)) > 0 in D -{0} (16) V (X(k + 1)) -V (X(k)) ≤ 0 in D (17)
then, X e = 0 is stable. Moreover, if

V (X(k + 1)) -V (X(k)) < 0 in D ( 18
)
then X e = 0 is asymptotically stable.

Adaptation of the Neural Network Formalism

Let us consider the autonomous system in ( 15), in which we assume that the equilibrium point of interest for the stability analysis is the point 0. Therefore,

f (0) = 0 (19)
Suppose V (X(k)) the Lyapunov function for the discrete time case and,

G(X(k)) = V (X(k + 1)) -V (X(k)) (20)
with

X(k) = [x 1 (k), ..., x n (k)] T .
In the discrete time case, sufficient conditions (a1) -(a4) for the point 0 of system to be asymptotically stable in the sense of Lyapunov are still the same than in the continuous part explained above. Conditions (v1) -( v3) and (d1) -(d3) which prove that V (X(k)) to have a local minimum at 0 and G(X(k)) to have a local maximum at 0 are also equivalent. However, the derivative of the function V (X) is not defined in the same way for the continuous and the discrete time case. Therefore, in this extension we will present how to define the new H G matrix. First and second derivatives of G(X(k)) can be calculated from (20) :

G l = ∂G ∂x l = n ∑ h=1 ∂V ∂x h ∂ f h ∂x l - ∂V ∂x l (21) G l p = ∂ 2 G ∂x l ∂x p = n ∑ h=1 ∂ 2 V ∂x h ∂x p ∂ f h ∂x l + n ∑ h=1 ∂V ∂x h ∂ 2 f h ∂x l ∂x p - ∂ 2 V ∂x l ∂x p with l = 1, ..., n and p = 1, ..., n (22)
In view of condition (v2), ( 21) and ( 22) result in, G l (0) = 0 for all l = 1, 2, ..., n

and

G l p (0) = n ∑ h=1 ∂ 2 V ∂x h ∂x p X(k)=0 ∂ f h ∂x l X(k)=0 - ∂ 2 V ∂x l ∂x p X(k)=0 (24)
On the basis of ( 24) we can write,

H G = [G l p (X(k) = 0)] = H V J -H V (25) since H V is symmetric. H G is a matrix whose ele- ments are G l p (0) and J is the Jacobian of f (X(k)) at X(k) = 0.
Let us assume that the Lyapunov function, V (X(k)), is represented by the neural network defined in ( 7) and ( 8). Therefore, the second derivative of G(X(k)) is computed as functions of the neural network weights,

G l p (0) = n ∑ j=1 ( K ∑ i=1 a i ∂ 2 σ(ν i ) ∂ 2 ν i w ji w li )J j p - n ∑ j=1 ( K ∑ i=1 a i ∂ 2 σ(ν i ) ∂ 2 ν i w ji w pi ) (26)
where J qr = ∂ f q ∂x r X(k)=0 q=1,...,n; r=1,...,n; l=1,...,n; p=1,...,n.

Therefore, t4) hold with the new definition of the matrix H G .

H V = [V qr (X(k) = 0)] V qr is given by (9) (27) H G = [G l p (X(k) = 0)] G l p is given by (26) (28) Point X(k) = 0 is asymptotically stable if condi- tions (t1) -(

The Proposed Algorithm

Optimization Scheme

First, for an appropriate Lyapunov function to be determined, values of the weights of the neural network should be calculated such that the conditions (t1)-( t4) are satisfied. To this end, a cost function, Q , should be selected so that positivity and negativity respectively of H V and H G are constrained. The symmetric matrix H G is negative definite if all its eigenvalues are negative.

Denote λ v i , i=1,...,n v the set of the n v eigenvalues of H V and λ g i , i=1,...,n g the set of the n g eigenvalues of H G .

DOA Maximization problem

Conditions (a2) and (a4) prove the asymptotic stability only "in a neighbourhood of 0". Consider that a Lyapunov function V(X) is given, by definition of D in ( 2) and (3), the system will converge to 0 from every initial point X 0 belonging to D. Thus, to maximize the DOA we search for an approximation D as large as possible, where the system is stable. Therefore, D = D would be the best possible case.

Consider Z a set of points obtained from a hypercube whose faces are gridded in order to cover a sufficiently large enough domain ⊂ X. Note that, Z has to be chosen large enough such as contains D. We denote P = max (ratio v , ratio dv ) where ratio v are the number of points X ⊂ Z where evaluated V(X)< 0 and ratio dv are the number of points X ⊂ Z where evaluated V (X)> 0. In order to maximize D, we have to minimize P.

Note that, there are other more intelligent methods to determine the gridding but so far, we used the one presented previously. For example, in the future, a possible approach is to use the gradient of the current Lyapunov function candidate during the optimization. An idea could be to analyze the first run of the optimization algorithm and the larger the gradient is, the thinner the gridding has to be in the next runs.

Constrained Implementation

We now formalize the problem as a constrained one to avoid the use of barrier function which would lead to a suboptimal problem [START_REF] Petridis | Construction of neural network based Lyapunov functions[END_REF], whose solution needs to be a posteriori verified. The scheme proposed here is flexible so that more complex problems such as exponential stability, robust stability or Input-to-state stability (ISS), will efficiently be tackled in future works.

The problem can be expressed in the general form of an optimization problem in which the cost function Q needs to be minimized.

To this purpose, we set down:

H V = H V × -1 Denote λ v as the eigenvalues of H V . - λ v = max(real(λ v )) - λ g = max(real(λ g )) - λ = max ( - λ v , - λ g )
According to ( 7) and ( 8), we denote α as the decision variables where:

α = [w ji , a i , h i , θ] (29)
Then, the cost function to be minimized has the following form:

min α Q If - λ > 0 Q = - λ Else Q = - 1 P + 1
which is a similar formulation of the cost function that can be find in [START_REF] Feyel | Robust control optimization with metaheuristics[END_REF] and has proven its efficiency. According to the definition of the problem

Q , a neural network candidate is a Lyapunov function if Q < 0. Assuming D ⊂ Z, best case, Q = -1, refers
to D = D. The P+1 avoids singularities when P = 0. The benefit of this formulation is its great flexibility: easy adaptation for a multitude of complex problems, extension to other type of stability and no additional parameter to tune for the penalty function. Besides, we have the guarantee that in the domain D the eigenvalues of H V and H G have respectively real positive and negative values. Finally, no parameter is needed.

Simulation Results

In this section, we apply our approach to one continuous time system and to two different discrete time systems to validate the extension. For clarity purposes, we have chosen some two dimensional examples so that the results can be easily plotted.

The entire test was performed on a machine equipped with an Intel Core i5 -8400H (2.5 GHz) processor and 16 GB RAM.

Parameter Settings

The parameters settings used in the tests are as follows:

• We consider 1 hidden layer and the number of neurons of this hidden layer is arbitrarily set to K = 12.

• In both continuous and discrete time cases, we set x as a rectangle of 21 × 21 centered at 0. Therefore, V (X) and V (X) are evaluated in 441 points in the range of each system.

• The number of searched variables α (29) is: K × (2n) + 1 = 49, and each of them has its search space interval arbitrarily defined by [-4; 4].

• The optimization method to calculate the weights of the neural network to compute a Lyapunov function is the Genetic Algorithm from the Global Optimization Toolbox in Matlab, used, for instance, in [START_REF] Krishna | Optimal design of roof-truss using ga in matlab. i-Manager[END_REF]. Other parameters that are not mentioned are the default value in GA, like the mutation rate for instance.

Continuous Time System

Example 1

Let us consider the following system:

ẋ1 = -tan x 1 + x 2 2 ẋ2 = -x 2 + x 1
The ranges for x 1 and x 2 are x 1 ∈ [-1, 1] and x 2 ∈ [-1, 1]. The stability of the origin is considered and the Figures 1 and2 show the result. The proposed method is compared to the method from [START_REF] Banks | Searching for Lyapunov functions using genetic programming[END_REF], based on genetic programming techniques. From Figures 1 and2, we can easily check that our method provides a larger D than the one found in the paper, which shows the efficiency of our approach. Therefore, the origin of the system is asymptotically stable.

Discrete Time Systems Example 2

Let us consider the following system:

     x + = 1 2 x + x 2 -y 2 y + = - 1 2 y + x 2
The ranges for x and y are x ∈ [-1, 1] and y ∈ [-1, 1]. The stability of the origin is considered and the Figures 3 and4 show the result. The Figure 4 shows the function V (X(k + 1)) -V (X(k)) of the system. We can easily check that the underlying function expressed by the network input-output relation is a Lyapunov function for this system. Therefore, the origin of the system is asymptotically stable.

Example 3

Let us consider the following system:

x + = -0.125y -0.125(1 -x 2 -y 2 )x y + = 0.125x -0.125(1 -x 2 -y 2 )y
The ranges for x and y are x ∈ [-0.75, 0.75] and y ∈ [-0.75, 0.75]. The stability of x e = [0; 0] is considered and the Figures 5 and6 show the result. The Figure 6 shows the function V (X(k + 1)) -V (X(k)) of the system. We can easily check that the underlying function expressed by the network input-output relation is a Lyapunov function for this system. Therefore, the equilibrium point x e = [0; 0] of the system is asymptotically stable.

Performance Measurement

Since the tested optimization algorithm is stochastic, a statistical analysis of the results is required. Thus, the performance measurement rule is as follows.

The cost function Q is defined in III.B. Each test of each system is subjected to 10 successive runs: we note the minimum value of Q obtained (Q min ), the mean value of Q (Q mean ), its standard deviation (Q std ) and finally, the average calculation time (t cpu ) taken to perform these 10 runs. The results are presented in Table 1. According to the definition of problems Q , a neural network candidate is a Lyapunov function if Q < 0 in Table 1. In these cases, we have the best domain if Q = -1.Therefore, we see that our extension of our approach to discrete time cases has very good results. The optimization algorithm finds a Lyapunov function 26 times out of 30 runs with these 3 examples. The best runs lead to the figures presented in section 5.

                 Q min = min i=1,...,nruns Q i Q mean = 1 n nruns ∑ i=1 Q i Q std = 1 n n runs ∑ i=1 (Q i -Q mean ) 2 (30)
In this paper, we investigate a new approach for the construction of a Lyapunov function modelled by a Neural network with the optimization of the domain of attraction for discrete time systems. We propose to use a constrained method such that the weights of the neural network are calculated in a way that is proven mathematically that the result function is a Lyapunov function while maximizing the DOA. Future work deals with proving more complex stability properties, such as the exponential, the robust stability and the Input-to-state stability (ISS). Besides, future work deals with developing this algorithm for real complex systems that we can find in industrial frameworks.

Figure 1 :

 1 Figure 1: The constructed Lyapunov function for this sys-

Figure 2 :

 2 Figure 2: The time derivative of the constructed Lyapunov function.

Figure 3 :

 3 Figure 3: The constructed Lyapunov function for this system.

Figure 4 :

 4 Figure 4: The difference function of the constructed Lyapunov function.

Figure 5 :

 5 Figure 5: The constructed Lyapunov function for this system.

Figure 6 :

 6 Figure 6: The difference function of the constructed Lyapunov function.

Table 1 :

 1 Algorithm Performance Measurement. Q min Q mean Q std t cpu/run(mn)

	Example 1	-1	-0.95 0.16	68.3
	Example 2	-1	-0.85 0.23	154.4
	Example 3	-1	-0.82 0.32	149.7


Dianailys Nuñez-Reyes

Kevin M Hickson
email: kevin.hickson@u-bordeaux.fr

The Reactivity of C( 1 D) with Oxygen Bearing Molecules NO and O 2 at Low Temperature

Keywords: temperature, rate constants, Laval nozzle reactor. Highlights

Gas-phase reactions; low

Introduction

The gas-phase reactions of atomic carbon are important processes in a variety of different environments and over a wide temperature range. In this respect, numerous experimental studies of the kinetics [1][2][3][4][5][6][7][8] and dynamics [9][10][11][12][13] of ground state atomic carbon reactions have been performed, providing detailed information regarding the mechanistic aspects of both radical-molecule and radical-radical reactions involving C( 3 P). The reactivity of excited state C( 1 D) atoms is less well understood, except for certain prototypical reactions such as C( 1 D) + H2 ® CH + H; a process which is simple enough to allow rigorous comparison between the most elaborate theoretical calculations and experiment [14][15][16][17][18][19][20][21]. In addition, a few detailed dynamics investigations of the reactivity of C( 1 D) atoms with hydrocarbons have been performed [22][23][24]. To address this issue, we have initiated a series of measurements to investigate both reactive and non-reactive processes involving C( 1 D) atoms, with a particular focus on their interactions with closed shell molecules [7,18,19,25,26]. Even less is known about the reactivity of C( 1 D) towards open shell species such as NO or towards biradicals such as O2, although the kinetics of the equivalent ground state carbon atom reactions have already been studied [3].

The ground state reagents C( 3 P) and NO(X 2 P) correlate adiabatically with CN(X 2 S + )/CN(A 2 P) + O( 3 P), CO(X 1 S + ) + N( 2 D) and CO(X 1 S + ) + N( 4 S) as exothermically allowed products. Indeed, several of these channels have been identified experimentally [27] and the reaction is characterized by a fast rate constant that increases to low temperature [3,[27][28][29][30][31]. For the C( 1 D) + NO(X 2 P) reaction, although an additional channel yielding CN(X 2 S + ) + O( 1 D) products is expected to be open, ground state CO(X 1 S + ) + N( 4 S) product formation is spin-forbidden despite the highly exothermic nature of this process (-567 kJ mol -1 ); a fact which might influence the overall reactivity. Interestingly, earlier kinetic studies of the C( 1 D) + NO reaction at room temperature [32] indicate that it could be significantly slower than the reaction between ground state reagents.

In the case of the C( 3 P) + O2(X 3 Sg -) reaction, the spin-allowed exothermic products are CO(X 1 S + ) + O( 3 P) and CO(X 1 S + ) + O( 1 D), although experimental evidence suggests that the favored exit channel is not the one leading to ground state products [31]. The equivalent reaction with C( 1 D) could lead to the formation of CO(X 1 S + /a 3 P) + O( 3 P) as possible products because the CO(X 1 S + ) + O( 1 D) channel is spin-forbidden. Previous studies of the C( 3 P)/C( 1 D) + O2 reactions also indicate that the C( 1 D) reaction occurs less efficiently than the C( 3 P) one [32] at room temperature, although no information exists regarding the relative reactivity of these processes over a range of temperatures.

Here we report the results of an experimental investigation of the C( 1 D) + NO and C( 1 D) + O2 reactions in the gas-phase. Rate constants for these processes were measured over the 50 -296 K temperature range using a Laval nozzle reactor coupled with pulsed laser photolysis to produce C( 1 D) atoms, employing a chemical tracer method to follow the reaction kinetics.

Section 2 presents the method used to perform this study. Section 3 describes the outcome of these experiments and discusses the validity of the results through modeling studies of the processes occurring within the reactor. Our conclusions are presented in Section 4.

Experimental Section

The experiments reported here were performed using a supersonic flow reactor operating in continuous mode [33,34], with the low temperature flows being generated by axisymmetric Laval nozzles. Although a wide range of nozzles is available in our laboratory employing a variety of carrier gases (Ar, N2 or N2/SF6 mixtures) to access temperatures in the range 50-240 K, only argon based nozzles were used during this study. The fast removal of excited state carbon atoms by molecular colliders excluded the use of any of the N2 based nozzles [25]. The measured and calculated flow properties of the Laval nozzles used here are listed in Table 1. a The errors on the Match number, density, temperature and mean flow velocity (1σ) are calculated from separate measurements of the impact pressure using a Pitot tube as a function of distance from the Laval nozzle and the stagnation pressure within the reservoir.

Measurements were also performed at room temperature, by removing the Laval nozzle and by reducing the flow velocity to eliminate pressure and temperature gradients within the reactor.

The pulsed laser photolysis method was used to generate C( 1 D) in these experiments, through the 266 nm multiphoton dissociation of CBr4 as the atom source using a frequency quadrupled Nd:YAG laser, with the laser beam aligned along the supersonic flow axis. This process produces mostly ground state C( 3 P) atoms with relative yields of C( 1 D) atoms in the 10-15 % range as already observed in a previous study [7]. A small argon flow directed into a vessel containing solid CBr4 maintained at a known pressure and temperature allowed CBr4 vapour to be carried into the Laval nozzle reservoir. CBr4 concentrations did not exceed 1.7 × 10 13 molecule cm -3 within the supersonic flow. C( 1 D) atoms could not be observed directly in the present experiments because no strong transitions originating from the 1 D state were accessible using our current detection method. Instead, small quantities of H2 or CH4 were added to the supersonic flow, allowing us to follow the H-atom products of the C( 1 D) + H2 ® CH + H reaction or through one of the major channels of the C( 1 D) + CH4 reaction (leading to H2CC/C2H2 + H + H as products) by vacuum ultraviolet laser induced fluorescence (VUV LIF) using the 1s 2 S ® 2p 2 P 0 Lyman-α transition at 121.567 nm. Both processes produce atomic hydrogen with yields close to unity [18,26], while C( 3 P) atoms themselves are unreactive with both H2 and CH4 at room temperature and below. A more detailed explanation of the methodology used to extract rate constants for the C( 1 D) + NO/O2 reactions is provided in the results and discussion section. Tunable light around 121.6 nm was produced through third harmonic generation by focusing a UV laser beam at 364.7 nm produced by a Nd:YAG pumped dye laser into a cell containing 210 Torr of krypton and 540 Torr of argon for phase matching purposes [7]. The resulting VUV radiation was collimated by a MgF2 lens that also served as the cell window. The cell was attached to the reactor at the level of the detection axis via a 75 cm sidearm containing baffles that effectively prevented most of the divergent UV light from entering the reactor. Moreover, the sidearm was constantly flushed with N2 or Ar during the experiments to prevent the attenuation of VUV radiation by residual reactive gases in this region. The VUV LIF emission itself was collected by a solar blind photomultiplier tube coupled with a boxcar integrator for signal acquisition and processing.

Gases were flowed directly from cylinders through calibrated digital mass flow controllers into the reactor. The stated purities are (Linde Ar 99.999%, O2 99.999%, Kr 99.99% CH4 99.9995%, Air Liquide N2 99.999%, H2 99.9999% and NO 99.9%).

Results and discussion

For any series of measurements, the H2 or CH4 concentrations were set to a fixed value (in large excess with respect to C( 1 D) atoms) so that in the absence of competing secondary reactions, the temporal evolution of the H-atom VUV LIF signal could be described using the formula

𝐼 # = 𝐴{𝑒𝑥𝑝(-𝑘′ / 𝑡) -𝑒𝑥𝑝(-𝑘′ 2 𝑡)} (1)
where A is the H-atom signal amplitude and 𝑘′ 2 and 𝑘′ / are the pseudo-first-order rate constants for H-atom formation and loss respectively.

𝑘′ / represents the first-order loss rate of atomic hydrogen and has contributions from diffusion 𝑘 4566,# as well as secondary reactions 𝑘 #89 [𝑋].

𝑘′ 2 governs the rise time of the H-atom signal and comprises the following terms:

= 𝑘 >( ? @)8A [𝑌] + 𝑘 >( ? @)8D [𝑍] + 𝑘 >( ? @)8>FG H [𝐶𝐵𝑟 L ] + 𝑘 >( ? @)8MG [𝐴𝑟] + 𝑘 4566,>( ? @)
𝑘 >( ? @)8A [𝑌] is the pseudo-first-order rate for the reaction of C( 1 D) with Y where Y is H2 or CH4, 𝑘 >( ? @)8>FG H [𝐶𝐵𝑟 L ] is the pseudo-first-order rate for the reaction of C( 1 D) with CBr4, 𝑘 >( ? @)8MG [𝐴𝑟] is the pseudo-first-order rate for the removal of C( 1 D) by Ar and 𝑘 4566,>( ? @) is the first-order diffusional loss rate of C( 1 D). All these terms are constant for any single series of measurements. 𝑘 >( ? @)8D [𝑍] is the pseudo-first-order rate for the reaction of C( 1 D) with Z where Z is NO or O2. NO or O2 are systematically varied so that 𝑘′ 2 varies as a function of [𝑍]

for any single series of measurements. Both NO and O2 were used in large excess with respect to C( 1 D) atoms so that the pseudo-first-order approximation was valid.

The measured H-atom VUV LIF signal recorded at 50 K as a function of delay time between lasers is displayed in Figure 1 in the presence and absence of NO. profiles such as those shown in Figure 1. The second-order rate constants were derived from weighted linear least squares fits to the data (blue and red solid lines).

The derived second-order rate constants are listed in Table 2 and displayed as a function of temperature in Figure 3. The quoted uncertainties represent the combined statistical and systematic errors. The statistical errors were derived from weighted fits such as those shown in T / K room temperature results of Braun et al. [35] for both reactions is quite poor. While the C( 1 D) + O2 reaction displays only small variations as a function of temperature in the 50 -296 K range, the C( 1 D) + NO reaction increases rapidly below 100 K, with the rate constant reaching a value of (9.7 ± 1.1) × 10 -11 cm 3 s -1 at 50 K; a factor of two increase with respect to the measured value at 296 K.

As the kinetics of the C( 1 D) + NO and C( 1 D) + O2 reactions were only followed indirectly in the present experiments, it was necessary to examine the potential effects of secondary reactions on the measured rate constants. Moreover, as C( 3 P) atoms are also generated in the reactor as a byproduct of CBr4 photolysis at 266 nm, it was important to evaluate the influence of the C( 3 P) + NO and C( 3 P) + O2 reaction products on the chemistry. In particular, any secondary H-atom sources could interfere with our rate constant measurements of the primary reactions. It is also possible that some C( 3 P) atoms might be produced through electronic quenching of C( 1 D) by O2 or NO, although no experimental or theoretical evidence exists for these non-reactive processes.

In the case of the C( 1 D)/C( 3 P) + O2 system, previous work [36] indicates that O( 1 D) is likely to be a major product for at least one of these processes and should therefore lead to H-atom production through the O( 1 D) + H2 ® OH + H reaction. To check for the possible influence of secondary H-atom production in our studies of the C( 1 D) + O2 reaction, we employed CH4 as the H-atom source in additional experiments performed at 127 K and at 296 K. While the Hatom yield of the C( 1 D) + CH4 reaction is close to unity over a wide temperature range [26], the corresponding O( 1 D) + CH4 reaction leads to only a 20 % yield of secondary H-atom products (with CH3O/CH2OH as the coproduct) [37], in stark contrast to the O( 1 D) + H2 reaction with a 100 % H-atom yield [38]. Consequently, if H-atom formation through secondary reactions of O( 1 D) atoms is important, we would expect to observe a noticeable difference between those experiments employing CH4 and those experiments employing H2 as the H-atom source. The lower panel of Figure 2 shows the results of two such experiments conducted at 127 K, with both methods essentially leading to the same value for the second-order rate constant.

In the case of the C( 1 D)/C( 3 P) + NO system, we also had to account for the reaction of ground state C( 3 P) atoms with NO which is known to lead to several different bimolecular products. ). The resulting pseudo-first-order rates were plotted against the corresponding NO concentration as shown in Figure S1. As expected, the y-axis intercept value increases with increasing H2, due to the increased rate of H-atom formation. Although the derived second-order rate constant is seen to increase slightly with H2 as shown by the slopes of the three datasets in Figure S1, the values are essentially identical when the experimental error bars are considered, as shown in Figure 3.

Previous work by

As a secondary check of the overall chemistry, numerical simulations of the reactions potentially involved in the C + O2 and C + NO systems were also performed using the differential integrator FACSIMILE. The reactions used in the simulations are listed in Tables S1 andS2 alongside their corresponding references when available [18,26,27,31,[37][38][39][40][41][42][43][44][45][46][47]. The branching ratios and/or temperature dependences for several of these processes are unknown, so it was necessary to estimate these quantities for a meaningful comparison with the experimental results. The aim of the simulations was to generate H-atom formation profiles over a similar range of reagent concentrations to those used in the experiments. These H-atom profiles were then fitted using expression (1), yielding a range of values for 𝑘′ 2 at a given temperature so that simulated second-order rate constants kSIM(T) for the C( 1 D) + O2 and C( 1 D) + NO reactions could be extracted from the slopes of simulated second-order plots. A comparison between the simulated and experimental H-atom profiles is shown in Figure S2.

The simulated second-order rate constants were then compared with the assumed input value used in the simulation, to check for any significant deviations brought about by secondary chemistry. Given the large uncertainties in the branching ratios and rate constants for several important reactions used in the simulations, we estimate an uncertainty of at least 20 % on the resulting simulated rate constants. This value was determined through a sensitivity analysis by varying the rate constants of most important secondary processes. For the C + O2 system, simulations were performed with both H2 and CH4 as the primary H-atom sources in a similar manner to the experiments. The excess reagent concentrations were obtained directly from their partial pressures and the initial C( 1 D) and C( 3 P) concentrations were estimated from the gasphase concentration of CBr4 and its absorption cross-section along with the photolysis laser fluence. For the C( 1 D) + O2 reaction, the simulated second-order rate constants at 50, 75 and 127 K are within 4 % of the input value when using H2 as H-atom source. At room temperature, the simulated rate constant is 21 % lower than the input value, indicating that secondary Hatom formation could lead to an error in the measured rate at this temperature. When the same simulations are performed at 296 K and 127 K using CH4 as the H-atom source, the deviations are reduced to less than 13 % at 296 K (the value at 127 K is within 4 % of the input value).

The larger discrepancy in simulations employing H2 arises mostly from the reaction of CH (produced by the C( 1 D) + H2 reaction) with O2, leading to atomic hydrogen as a major product (channels 7a and 7c in Table S1). O( 1 D) atoms produced by reaction 3a were also thought to be an important source of H-atoms in the simulations through the O( 1 D) + H2 reaction (reaction 6). However, when the branching ratio of this channel is set to zero, the discrepancy between the input and output rate constants for reaction 1 increases slightly at all temperatures, clearly indicating the influence of non-linear effects in this complex system of coupled reactions. When CH4 is used as the H-atom source, the only secondary process leading to H-atom production is the O( 1 D) + CH4 reaction (reaction 11b) with a much smaller yield than the O( 1 D) + H2 one [37]. Nevertheless, as our room temperature measurements of the C( 1 D) + O2 reaction with H2 and CH4 result in very similar second-order rate constants, this might indicate that the real branching ratios for the H-atom production channels of the CH + O2 reaction (7a and 7c) are smaller than the ones used in the simulations.

For the C( 1 D) + NO reaction, where only H2 was used as the H-atom source, the discrepancy between the simulated second-order rate constants and the input value increases as the temperature decreases, from 5 % at 296 K to 25 % at 50 K. The reactions that are responsible for most of this effect are those leading to O( 1 D) formation, as these atoms react rapidly with H2 to form atomic hydrogen (reaction 6). Indeed, if reaction 6 is switched off in the simulations, the discrepancy between the input and output second order rate constants decreases to a maximum of 10 % at 50 K. In particular, the reaction of N( 2 D) atoms with NO (reaction 17a) is the major source of O( 1 D) atoms in these simulations, with N( 2 D) atoms being mostly produced by the reaction of ground state carbon atoms with NO (reaction 13b). As earlier measurements by Bergeat et al. [27] only yielded the sum of N-atom production channels (reactions 13a + 13b) for this reaction, we set the rate constants for each of these channels to overestimated in these simulations. Moreover, the branching ratios for the individual channels of the N( 2 D) + NO reaction have never been measured [45].

We also tested the effect of varying the H2 concentration for the C( 1 D) + NO simulations at 127 K, to validate the measurements performed at three different H2 concentrations. With [H2] = 0.7 and 1.7 × 10 14 cm -3 , the simulated second-order rate constants are identical, with a value of 4.6 × 10 -11 cm 3 s -1 (using the input value of 5.1× 10 -11 cm 3 s -1 ). With [H2] = 3.5 × 10 14 cm -3 , the simulated second-order rate constant rises to 5.2 × 10 -11 cm 3 s -1 . The main experiments at 127 K were all performed using [H2] = 1.7 × 10 14 cm -3 .

It is interesting to compare the present results for the C( 1 D) + NO and C( 1 D) + O2 reactions with previous work on the equivalent ground state carbon atom reactions. Chastaing et al. [3] and Geppert et al. [31] investigated the kinetics of the C( 3 P) + O2 and C( 3 P) + NO reactions over an extended temperature range (15-295 K). They found that the C( 3 P) + O2 reaction becomes faster as the temperature falls with rate constants ranging from (4.09 ± 0.16) × 10 -11 cm 3 s -1 at 295 K to (11.2 ± 0.7) × 10 -11 cm 3 s -1 at 15 K, in contrast with our measured rate constants for the C( 1 D) + O2 reaction which do not display a marked temperature dependence.

Chastaing et al. [2] hypothesized that the increase in reaction rate as the temperature decreases could originate from the changing population distribution over spin-orbit levels of C( 3 P) which will be at thermal equilibrium in Laval nozzle type flows due to collisions with the carrier gas.

If reaction occurs adiabatically, only from specific spin-orbit states, this could lead to large changes in the overall reaction rate as the temperature is varied. In the case of the C( 1 D) + O2 reaction, the fine structure states of the C( 1 D2) atom are degenerate so that the population distribution is unaffected by temperature. In this instance, a comparable effect will not occur for the excited state carbon atom reaction with O2.

For the C( 3 P) + NO reaction, Chastaing et al. [3] showed that the rate constants increased slightly with decreasing temperature from (1.2 ± 0.09) × 10 -10 cm 3 s -1 at 295 K to (2.28 ± 0.17) × 10 -10 cm 3 s -1 at 15 K in a similar manner to our results for the corresponding C( 1 D) atom reaction which displays a twofold increase in reactivity over the 50 K -296 K range (although the C( 3 P) atom reaction is 3 to 4 times faster than the C( 1 D) atom one). In the case of the C( 3 P) + NO reaction, the situation is particularly complex with C( 3 PJ) and NO( 2 PW) reagents correlating with 36 different states, 18 of which correlate with exothermic products O( 3 P) and CN(X 2 S + ) [48]. Consequently, the changing spin-orbit populations of both C( 3 PJ) and NO( 2 PW) could result in a wide range of temperature dependences, as certain adiabatic product channels become more or less favored as the temperature changes. In the case of the C( 1 D2) + NO( 2 PW) reaction, only the population distribution of the spin-orbit levels of NO changes as a function of temperature. As a result, we might expect to observe a negative temperature dependence for this process if reaction between C( 1 D2) and ground electronic state NO( 2 P1/2) is one of the preferred pathways.

Interestingly, the reactions of C( 1 D) with the radical species O2 and NO display rate constants that are smaller than those of the reactions of C( 1 D) atoms with closed shell molecules in general [32]. This is likely to be related to the large number of electronic states correlating adiabatically with both sets of reagents, leading to small electronic degeneracy factors (that is the fraction of the thermal population on reactant electronic states that lead to product formation) and a correspondingly lower reactivity.

Conclusions

We have measured rate constants for the reactions of excited state carbon atoms, C( 1 D), with oxygen bearing molecules O2 and NO over the 50 -296 K temperature range using a supersonic flow reactor. Pulsed laser photolysis was employed to generate C( 1 D) atoms and a chemical tracer method, based on the addition of either H2 or CH4 to the flow, was used to follow the progress of the reaction. H-atoms produced by the tracer reaction were detected by vacuum ultraviolet laser induced fluorescence. The rate constants for the C( 1 D) + O2 reaction vary only slightly over the experimental temperature range. In contrast, the rate constants for the C( 1 D) + NO reaction increase as the temperature falls. As these reactions could not be followed directly, various test measurements were performed to validate the measured rate constants. In addition, numerical simulations were performed to evaluate the impact of secondary reactions on the rate constant measurements. Taken together, the test measurements and simulations indicate that secondary reactions have only a minor influence on the derived rate constant values.

Figure 1 H

 1 Figure 1 H-atom formation curves recorded at 50 K for the C( 1 D) + NO reaction using the H2

Figure 2

 2 Figure 2Pseudo-first-order rate constants as a function of coreagent concentration. Top panel:

Figure 2 .Table 2

 22 Figure 2. The systematic errors are estimated at the level of 20 % of the nominal rate constant

Figure 3

 3 Figure 3 Rate constants for C( 1 D) reactions as a function of temperature. C( 1 D) + NO: (solid

  Bergeat et al. [27] has established that this reaction leads to 60 % production of atomic nitrogen (N( 4 S) + N( 2 D)) with the remaining 40 % leading to O( 3 P) formation at room temperature. Of these channels, the formation of N( 2 D) will lead to secondary H-atom formation through the N( 2 D) + H2 ® NH + H reaction. Moreover, the C( 1 D) + NO reaction itself could lead to O( 1 D) formation (in addition to the N( 2 D) channel), as CN(X 2 S + ) + O( 1 D) products are also adiabatically and energetically accessible for these reagents. Consequently, the reactions of O( 1 D) and N( 2 D) with H2 could both lead to secondary H-atom formation. To test for the possible influence of these reactions, we ran additional test experiments at 127 K to determine the rate constant for the C( 1 D) + NO reaction. In these experiments, H-atom formation profiles were recorded for a range of NO concentrations using three different fixed H2 concentrations ([H2] = 0.7, 1.7 and 3.5 × 10 14 cm -3

Table 1

 1 Characteristics of the supersonic flows

	Mach number	2.0 ± 0.03 a	3.0 ± 0.1	3.9 ± 0.1
	Carrier gas	Ar	Ar	Ar
	Density (´10 16 cm -3 )	12.6 ± 0.3	14.7 ± 0.6	25.9 ± 0.9
	Impact pressure (Torr)	10.5 ± 0.2	15.3 ± 0.5	29.6 ± 1.0
	Stagnation pressure (Torr) 13.9	34.9	113
	Temperature (K)	127 ± 2	75 ± 2	50 ± 1
	Mean flow velocity (ms -1 ) 419 ± 3	479 ± 3	505 ± 1
	Chamber pressure (Torr)	1.5	1.2	1.4

  Uncertainties on the measured rate constants represent the combined statistical and systematic errors as explained in the text. d Rate constants determined with H2 as the chemical tracer. e Rate constants determined with CH4 as the chemical tracer.

								(2.8 ± 0.6)
							9-9.8	(2.7 ± 0.6)
	75±2	37	0.8	0-8.0	(6.7 ± 1.4) 36	0.9 d	0.4-7.2	(2.7 ± 0.5)
	50±1	36	1.2	0-8.9	(9.7 ± 2.0) 41	1.2 d	0.6-9.7	(3.4 ± 0.7)

a Uncertainties on the calculated temperatures represent the statistical (1s) errors obtained from Pitot tube measurements of the impact pressure. b Number of individual measurements;

c

Acknowledgements K.M.H. and D.N.R acknowledge support from the French program "Physique et Chimie du Milieu Interstellaire" (PCMI) of the CNRS/INSU with the INC/INP co-funded by the CEA and CNES.