Ramzi Ben Mhenni

Sébastien Bourguignon

Jérôme Idier

A GREEDY SPARSE APPROXIMATION ALGORITHM BASED ON L1-NORM SELECTION RULES

Keywords: Sparse optimization, greedy algorithm, 1 -norm regularization, deconvolution

We propose a new greedy sparse approximation algorithm, called SLS for Single L1 Selection, that addresses a least squares optimization problem under a cardinality constraint. The specificity and increased efficiency of SLS originate from the atom selection step, based on exploiting 1 -norm solutions. At each iteration, the regularization path of a least-squares criterion penalized by the 1 norm of the remaining variables is built. Then, the selected atom is chosen according to a scoring function defined over the solution path. Simulation results on difficult sparse deconvolution problems involving a highly correlated dictionary reveal the efficiency of the method, which outperforms popular greedy algorithms when the solution is sparse.

INTRODUCTION

Sparse approximation is a very active research topic, in which we search for a sparse vector x ∈ R M (i.e., with a large number of zero components), to approximate data y ∈ R N as a linear combination of a small number of atoms, i.e., columns from a dictionary H ∈ R N ×M . Applications range from sparsity-enhancing inverse problems in geophysics [START_REF] Taylor | Deconvolution with the l1 norm[END_REF][START_REF] Mendel | Optimal seismic deconvolution: An estimationbased approach[END_REF] or nondestructive testing [3,4], compression and denoising [START_REF] Eldar | Compressed Sensing: Theory and Applications[END_REF], or subset selection in Statistics [START_REF] Miller | Subset selection in regression[END_REF].

The sparse approximation problem can be formulated by minimizing the least-squares approximation error and imposing an upper bound on the number of non-zero elements of the solution [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF]:

min x 1 2
y -Hx 2 2 subject to (s.t.)

x 0 ≤ K, (1)
where x 0 denotes the 0 "norm": x 0 := Card{x j |x j = 0}.

Finding the best K-sparse solution essentially resorts to combinatorial optimization, which is often not possible for high-dimensional problems-the problem (1) is NP-hard [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF]. Therefore, many works in signal processing and statistics have proposed computationally efficient, suboptimal, methods, that can be classified into two categories [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF].

On the one hand, continuous relaxation methods replace the 0 norm by a continuous surrogate function. Problem (1) is then converted into a continuous optimization problem. In particular, many works have considered the convex case with the 1 norm x 1 := n |x n | [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF]. On the other hand, partial combinatorial exploration strategies generate a sequence of sparse iterates by performing simple transitions on the support (the set of non-zero components) of x. Among this category fall the well-known Matching Pursuit (MP) [START_REF] Mallat | Matching Pursuits with Time-Frequency Dictionaries[END_REF], Orthogonal Matching Pursuit (OMP) [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF] and Orthogonal Least Squares (OLS) [START_REF] Chen | Orthogonal least squares methods and their application to non-linear system identification[END_REF] greedy methods. However, the selection step of a new atom in such methods is very sensitive to interferences between the different atoms, in particular in the case of highly correlated dictionaries [START_REF] Soussen | From Bernoulli Gaussian Deconvolution to Sparse Signal Restoration[END_REF]. More sophisticated strategies try to circumvent this issue (see for example [START_REF] Selesnick | Sparse signal estimation by maximally sparse convex optimization[END_REF]). Among them, the Single Best Replacement (SBR) algorithm [START_REF] Soussen | From Bernoulli Gaussian Deconvolution to Sparse Signal Restoration[END_REF] includes possible removal steps allowing to correct erroneously detected components, CoSAMP [START_REF] Needell | Cosamp: Iterative signal recovery from incomplete and inaccurate samples[END_REF] considers including 2K components at each iteration and selects the K most likely ones, and A OMP [START_REF] Karahanoglu | A* orthogonal matching pursuit: Best-first search for compressed sensing signal recovery[END_REF] performs a partial tree-search exploration step for the selection step.

While the 1 -norm convex formulation exactly solves an approximate problem, the second class of methods returns a local optimum for the exact problem. In the compressed sensing framework [START_REF] Eldar | Compressed Sensing: Theory and Applications[END_REF], conditions on H have been established, for which such approaches are ensured to solve the initial problem [START_REF] Taylor | Deconvolution with the l1 norm[END_REF]. However, in the case of inverse problems where H is ill-conditioned, or if the mutual correlation of its columns is high, such guarantees are lost and, in practice, they often fail in finding the global optimum [START_REF] Bourguignon | Exact sparse approximation problems via mixed-integer programming: Formulations and computational performance[END_REF]. The relative performance of methods is then evaluated under a compromise criterion between the quality of the solution and its computational cost.

In this paper, we propose an algorithm which gathers advantages of the two classes of methods. It essentially consists of a greedy strategy, where the selection rule at each iteration is based on exploiting 1 -norm solutions. The number of iterations is then controlled by the sparsity level K of the searched solution, limiting the computational burden. Moreover, the selection of each new atom, based on solving a convex optimization problem, is expected to be more robust to interferences (high correlations) between the different atoms than standard greedy methods. More precisely, at a given iteration, an 1 -norm-penalized problem is built, where the 1 norm operates on the non-selected variables, and the regularization path is computed by an homotopy continuation algorithm [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Donoho | Fast solution of 1-norm minimization problems when the solution may be sparse[END_REF]. A heuristic selection rule is then proposed, which considers the predominant variable along the regularization path.

Section 2 presents classical greedy forward algorithms and discusses the limitations of their variable selection rules. In Section 3, we introduce our new selection rule, and implementation issues are addressed. The resulting algorithm, called SLS for Single 1 Selection, is then evaluated and compared to other sparse approximation algorithms in Section 4.

LIMITATIONS OF FORWARD GREEDY ALGORITHMS

Forward greedy methods start from an empty set and iteratively construct a sparse solution, by alternating between two steps: a new atom is selected by maximizing a score function, denoted by F, and then the model is updated. Let S denote the index set of the variables already selected (the current support of the solution), and let H S denote the matrix composed of the columns of H indexed by S. Similarly, S indexes the remaining variables. We also denote the size of S (respectively, of S) by k (respectively, k). The general principle of forward selection algorithms is given in Algorithm 1. In the sequel, we suppose that all columns in H have unit norm. For both MP and OMP algorithms, the selected atom is the most correlated to the residual:

F(O)MP(j) = |h T j r|, j ∈ R M , ( 2
)
where hj is the j-th column of H. OMP includes an additional orthogonalization step, which updates the current solution on its support by:

x S = argmin x S ∈R k y -H S x S 2 2 = (H T S H S ) -1 H T S y,
such that the residual is orthogonal to the model H S x S . In this case, the score function ( 2) is computed only for j ∈ S.

For OLS [START_REF] Chen | Orthogonal least squares methods and their application to non-linear system identification[END_REF], the approximation error is minimized among all possible supports including one new component:

 = argmin j∈S min x S ∈R k y -H S∪{j} x S∪{j} 2 2 ,
which amounts to

FOLS(j) = H S∪{j} (H T S∪{j} H S∪{j} ) -1 H T S∪{j} y 2 .
The atom selection stage of MP, OMP and OLS is a one-step procedure, which is a short-term vision of the selection issue. In the case of highly correlated dictionaries, this represents a major limitation of such greedy algorithms. For illustration purposes, let us consider a sparse deconvolution problem, where H is composed of shifted versions of the impulse response of the filter, and adjacent atoms give very close contributions to the model. The toy example in Figure 1 (a) and (b) illustrates a situation where x is composed of two close spikes, giving strongly overlapping echoes in the data y. The score function for the first iteration of both OMP and OLS is F(j) = |h T j y|, and is shown in Figure 1 (c). It is maximal for the index located in the middle of the two true indices, thus selecting a wrong atom -even in the noiseless case.

SINGLE L1 SELECTION (SLS)

In this Section, we build a new forward method, where a more accurate selection rule is proposed, based on a global, convex, formulation of the problem.

Selection rule based on 1 problems

Our selection rule exploits the solutions of 1-norm-based problems. Recall the example in Figure 1. It is clear that a joint approach to the sparse estimation problem would be less sensitive to interferences between the atom contributions, by allowing non-zero weights to the true atom locations. Therefore, we propose to base the selection of a new atom by considering the following 1-norm optimization problem at each iteration of a greedy forward procedure:

min x S ∈R k ,x S ∈R k 1 2 y -H S x S -H S x S 2 2 + λ x S 1, (3)
Similarly to OLS, such a criterion allows the joint re-estimation of the amplitudes of previously selected components (whereas they are fixed in the residual for the selection rule of OMP). It also jointly estimates a sparse vector for the remaining ones, which is not restricted to a single non-zero component as in OLS.

We remark first that the problem in Eq. ( 3) can be recast as an optimization problem in variables x S only: indeed, for a given x S , the solution in x S is explicit. Its expression can be inserted into the least-squares term in (3), so that the problem amounts to the following standard 1-norm-penalized problem:

P λ : min x S ∈R k 1 2 y -H S x S 2 2 + λ x S 1, (4)
where y := Py, H S := PH S and P := I k -H S (H T S H S ) -1 H T S , with I k the k × k identity matrix.

Exploiting the regularization path

We now build a selection rule based on problem P λ in (4). Instead of considering a particular solution for a given value of λ, we propose the following two-step strategy.

Homotopy -We first build the solution path, that is, the set of solutions of P λ for λ ∈ [λ min , λ max ], by the homotopy continuation algorithm, which is acknowledged as a very efficient algorithm for 1-norm optimization problems when the solution is sparse [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Donoho | Fast solution of 1-norm minimization problems when the solution may be sparse[END_REF]. Starting from λ max := H T S y ∞ (such that the solution is identically zero for λ ≥ λ max ), homotopy iteratively identifies the different breakpoints, denoted λ (p) , p = 1, . . . , P , that lead to changes in the support configuration (the set of positive, negative and zero variables). Between two breakpoints, the solution is linear in λ and has an analytical expression. Implementation details can be found, e.g., in [START_REF] Donoho | Fast solution of 1-norm minimization problems when the solution may be sparse[END_REF]. Figure 2 illustrates the iterations performed by the homotopy algorithm.

x * 1 x * 5 x * 2 x * 3 x * 4 λ P λ (0) λ (1) λ (2) λ (3) λ (4)
λ (5) λ (6) x * λ stop after P iterations OMP selection SLS selection Fig. 2. Selection rule computed over the solution path obtained by the homotopy method, as λ decreases from λ (0) = λ max to λ (P ) . Path corresponding to non-zero variables in x * , renamed x * 1 to x * 5 . Although x1 is the first component to be non-zero, x2 is the most prominent variable on [λ (P ) , λ (0) ].

Scoring -Then, we use the score function defined over the set of solutions for λ ∈ [λ min , λ max ] by:

FSLS(j) = P p=0 |x (p) j |, j ∈ S, (5)
where x (p) denotes the solution of P λ at the breakpoint λ = λ (p) . The score function can be updated recursively at each iteration of the homotopy algorithm. This heuristic rule aims at selecting the most prominent component over the regularization path [λ min , λ max ].

In order to control the computing time with a simple rule, we limit the number P of iterations of the homotopy algorithm (which defines λ min ), to P = c(K -k), where K -k is the number of components that still need to be included in the greedy search, and c is a parameter whose tuning will be discussed in Section 4. Such an empirical rule follows the idea that a larger solution path is preferred for the first iterations of the greedy procedure, because more competition may exist between atoms, so the solution may be less stable as a function of λ.

The SLS selection rule is illustrated on a simulated regularization path in Figure 2. Remark that in the limiting case c → 0 (λ → λ max ), it corresponds to selecting the first non-zero component in the regularization path, that is, arg max j∈S h T S j y, corresponds to the OMP selection rule. On the toy example of Section 2, the first iteration of the SLS algorithm allows the correct selection of one true atom -see Figure 1 (e). Then, after two iterations, the true support is correctly estimated, as shows Figure 1 (f). Algorithm 2 summarizes the computation of FSLS via the homotopy algorithm.

At each iteration k of SLS, we call the homotopy algorithm to build the solution path in [λ (P ) , λ (0) ]. From the computational point of view, each step of the homotopy algorithm essentially involves the rank-one update of a linear system to compute d. Our implementation makes use of the blockwise inversion technique to solve the linear systems in an incremental way.

SIMULATION RESULTS

We evaluate the performance of the SLS algorithm, compared to several well-known sparse estimation algorithms: OMP [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF], 1- Algorithm 2: Homotopy algorithm for computing the score function F. Indexation by I returns non-zero variables.

1 Set p = 0, s (0) = 0 and x (0) = 0.

2 Initialize λ (0) = λ max and I = {arg max j∈S h T S j y}. norm regularization or BP (Basis Pursuit), computed here by the homopotopy algorithm [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], OLS [START_REF] Chen | Orthogonal least squares methods and their application to non-linear system identification[END_REF], SBR1 [START_REF] Soussen | From Bernoulli Gaussian Deconvolution to Sparse Signal Restoration[END_REF], Subspace Pursuit2 [START_REF] Needell | Cosamp: Iterative signal recovery from incomplete and inaccurate samples[END_REF], accelerated Iterative Hard Thresholding 3 (IHT) [START_REF] Blumensath | Iterative Thresholding for Sparse Approximations[END_REF] and A OMP 4 [START_REF] Karahanoglu | A* orthogonal matching pursuit: Best-first search for compressed sensing signal recovery[END_REF]. All algorithms are implemented in Matlab and are tuned in order to obtain solutions with the true sparsity level K.

Algorithms are tested on difficult sparse deconvolution problems, with an up-sampled convolution model in order to achieve high-resolution spike locations [START_REF] Carcreff | Resolution enhancement of ultrasonic signals by up-sampled sparse deconvolution[END_REF]. The problem is underdetermined with M = 1, 000 and N = 350. Columns of H are then highly correlated, with mutual coherence max i =j |h T i hj| = 0.81. White Gaussian noise is then added with SNRdB = 10 log Hx realizations of the sparse sequence and of noise. A typical signal is shown in Figure 3. We respectively denote by x and x true the estimated and the true solutions. In order to focus on the capacity of the methods to find the best solution and to detect the correct support, we consider several error measures:

• the quadratic error EQ = y -Hx 2 2 ; • the support error ES = bb true 0, where bj = 1 (respectively, bj = 0) if xj = 0 (respectively, if xj = 0);

• the exact recovery rate Eexact, which is the average number of instances for which the algorithm correctly locates the support of the true sequence.

We first analyze the efficiency of our SLS algorithm as a function of the control parameter c. Figure 4 shows errors EQ and ES by varying c from 0.5 to 5, as a function of the computation time, compared to the previously introduced algorithms, and for sparsity levels K ranging from 10 to 40. For small c, SLS still behaves as a fast and poorly performing algorithm on such difficult problems (recall that for c = 0, SLS identifies with OMP). As could be expected, the quality of the solution improves as c increases, at the expense of the increase of the computation time. In particular, for K ≤ 30, the quality of solutions provided by OLS and SBR is achieved with smaller computation time (c ∼ 1.5), and much better solutions than OLS and SBR are obtained for their respective computation time (c ∼ 3). For less sparse solutions (K ≥ 40), SLS behaves like OLS and SBR on such "quality vs. time diagram", although better solutions can be obtained with higher computation times. A OMP competes with SLS for K = 20, but takes much more time for greater values of K.

From Figure 4, we remark that the value of c = 3 is a good compromise between computation time and solution quality. Thus, we fix it at this value for the next comparisons. Figure 5 presents quantitative results for all algorithms as a function of K. The performance of SLS is clearly the best one in terms of solution quality with respect to the three error criteria. Finally, SLS has a lower computation time than OLS and SBR up to K = 30 and slightly higher for K = 40. Compared to A OMP, SLS has a slightly higher computation time up to K = 20, but A OMP becomes computationally too expensive. Note that other fast algorithms (OMP, BP, IHT and SP) are always much faster than SLS-but always give worse solutions.

CONCLUSIONS

We proposed a novel sparse approximation greedy algorithm, called SLS for Single L1 Selection, whose atom selection rule is based on the 1-norm solutions. This selection rule exploits the regularization path of the least squares criterion penalized by the 1-norm, with a view to improve the robustness of the selection step. We empirically showed that SLS outperforms other standard greedy methods to solve difficult inverse problems where the dictionary columns are highly correlated such as in deconvolution. Moreover, the computing cost of SLS can be made comparable to that of OLS and SBR, exploiting the fact that only partial 1 regularization paths are needed.

Algorithm 1 : 4 5 (

 145 Forward Selection greedy algorithm 1 Initialization: S = ∅ 2 while |S| < K do 3 Variable selection:  = argmax j∈S F(j) ; Support update: S ← S ∪ {}; If needed) update the estimate x S and the residual r = y -H S x S ; 6 end Result: support S and solution x S

Fig. 1 .

 1 Fig. 1. a) Toy example where standard greedy selection fails. (a) Contribution of two atoms and (b) corresponding noise-free data. (c) Scoring function for OMP and OLS at first iteration and (d) OLS solution after two iterations. (e) Scoring function at the first iteration of the proposed SLS algorithm and (f) SLS solution after two iterations. Circles (resp. stars) locate the true (resp. estimated) spikes.

3 while p ≤ P do 4 p ← p + 1 5% 6 d 7 % 8 γ 9 % Update the solution % 10 x 11 % 12 F 13 %

 41678910111213 Compute the change of direction % Compute the step length γ (p) % (p) : smallest decrease in λ to reach a new breakpoint Recursive computation of the score function % Update I % 14 end 15 Return the score F.

2 2=Fig. 3 .

 23 Fig. 3. Typical sparse deconvolution data y.

Fig. 4 .

 4 Fig. 4. Performance of SLS (*) as a function of the computation time, by varying parameter c from 0.5 (left) to 5 (right), compared to different sparse approximation methods. The tuning with c = 3 is highlighted with a red marker.

Fig. 5 .

 5 Fig. 5. Estimation errors and computation time as a function of K for SLS with c = 3 (*) and other sparse approximation methods.

© 2020 IEEE. This is the author's version of an article that has been published in ICASSP 2020. Changes were made to this version by the publisher prior to publication. The final version of record is available at http://dx.doi.org/10.1109/ICASSP40776.2020.9054670 Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

Implementation was taken from the multimedia material linked to[START_REF] Soussen | From Bernoulli Gaussian Deconvolution to Sparse Signal Restoration[END_REF].

B. Sturm's implementation taken at http://media.aau.dk/null_


Qingyong Meng
email: qingyong.meng@nwpu.edu.cn

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

Kejie Shao

Jean-Christophe Loison

Dong H Zhang
email: zhangdh@dicp.ac.cn

Theoretical and Experimental Investigations of Rate Coefficients of O ( 1 D) + CH 4 at Low Temperature

Keywords: RPMD, NN PES, VUV-LIF, Rate coefficients, O ( 1 D) + CH

Rate coefficients of the barrierless O ( 1 D) + CH 4 reaction are determined both theoretically and experimentally at 50 K -296 K. For the calculations, ring polymer molecular dynamics (RPMD) simulations are performed on the basis of a new neural network potential energy surface (PES) in reactant asymptotic part. Only reactant asymptotic part PES is constructed because of its barrierless and exothermic properties.

Experimentally, the reaction rate coefficients are measured using a supersonic flow reactor. Pulsed laser photolysis of O 3 molecules is used as the source of O ( 1 D) atoms, which are detected directly through vacuum ultraviolet laser induced fluorescence at 115 nm. The branching ratio for H atom production is measured by comparing the H atom yields of the O ( 1 D) + CH 4 and O ( 1 D) + H 2 reactions. At T ≥ 75 K, good agreement between theoretical and experimental rate coefficients is found, while at 50 K, the larger difference is discussed in detail.

I. INTRODUCTION

Due to its significance in both atmospheric and combustion chemistry, over the past decade the reaction of methane with O ( 1 D) has attracted great attention both experimentally [1][2][3][4][5][6][START_REF] Sander | Chemical Kinetics and Photochemical Data for Use in Atmospheric Studies[END_REF] and theoretically [START_REF] González | [END_REF][9][10][11], helping to unravel the dynamics and kinetics of this typical polyatomic complex-forming barrierless reaction. These studies [1][2][3][4][START_REF] González | [END_REF][9][10] have determined an approximately temperature-independent rate coefficient of ∼ 1.5 (unless otherwise specified, the rate is given in 10 -10 cm 3 • s -1 ) at the room temperature.

On the basis of a pseudo-triatomic three-dimensional (3D) ground-state ( 1 A ) potential energy surface (PES) [START_REF] González | [END_REF], González et al. [START_REF] González | [END_REF] and Ben Bouchrit et al. [10] reported rate coefficients of the O ( 1 D) + CH 4 → OH + CH 3 reaction from quasi-classical trajectory (QCT) [START_REF] González | [END_REF] and timeindependent quantum mechanical (TIQM) [10] reactive scattering calculations, respectively. They [START_REF] González | [END_REF]10] found (i) that quantum effects do not strongly affect the resulting rate coefficients and (ii) a rather good agreement between computations and experiments was obtained. Despite the existence of agreement at room temperature, it is still necessary to perform full dimensional (12D) quantum dynamics calculations because the CH 3 fragment was treated as a pseudo-triatom in previous calculations [START_REF] González | [END_REF]10]. Furthermore, below room temperature there are neither experimental measurements nor full dimensional quantum simulations available in the literature for this reaction.

Recently, Shao et al. [11] constructed a global full-dimensional ground-state PES of the O ( 1 D) + CH 4 system by permutationally invariant fitting to roughly 3.4 × 10 5 ab initio energy points with the maximum fitting error of dozen of meV. Since the fitting error of this global PIP PES is rather large, the dynamics simulations for the title reaction may only predict the rate coefficients with rather a large uncertainty. In this work, we re-construct a new, more accurate PES on the basis of the previous ab initio energy points [11]. Based on this new PES, ring-polymer molecular dynamics (RPMD) [12][13][14][15][16] simulations are performed to accurately and efficiently compute the rate coefficients. The RPMD method [12][13][14][15][16] exploits the isomorphism between the statistical properties of the quantum system and those of a fictitious ring polymer consisting of many replications of the original system, which are called beads, connected by harmonic springs. It has been shown that RPMD can be used to describe complex-forming reactions [17][18][19] and that RPMD is capable of providing fairly accurate estimates of the rate coefficients of gas phase polyatomic reactions involving methane [16,20,21].

Experimentally, the title reaction is studied over the range 50 K -296 K using the contin-uous supersonic flow method [22,23] by following the kinetics of O ( 1 D) atom decay directly [24]. Temperature dependent branching ratios are also measured by following the production of H atoms. In this work, the temperature range of 50 K -296 K is chosen because temperature in the Earth's atmosphere [25] is in the range of 100 K -300 K. This paper is organized as follows; in Section II, we will describe the experimental and theoretical methods. Sections III and IV presents the results and discussions, respectively. Finally, Section V concludes with a summary.

II. EXPERIMENTAL AND THEORETICAL METHODS

A. Experimental Method

All experiments are performed using a continuous supersonic flow apparatus which has been described in detail in previous work [19,[22][23][24]26] at specified temperatures of 127 K, 75 K, and 50 K through the use of different argon based Laval nozzles, with experiments at 296 K being performed without a nozzle. As the non-reactive quenching of O ( 1 D) by N 2 is rapid at low temperature [24], it was impossible to use our nitrogen based nozzles for this investigation. The measured and calculated flow characteristics can be found in Table 1 of Grondin et al. [24].

O ( 1 D) atoms are created in the present experiments by the pulsed UV photo-dissociation of O 3 in situ in the supersonic flow at 266 nm with energies around 21 mJ. O 3 itself is produced by the continuous UV irradiation of O 2 as described by Grondin et al. [24]. The O ( 1 D) atoms are followed directly through resonant VUV-LIF using the 3s 1 D → 2p 1 D transition at 115.215 nm. Tunable radiation at this wavelength is produced by frequency doubling 691.29 nm light before focusing the doubled output into a cell containing 100 Torr of Xe and 230 Torr of Ar for third harmonic generation. For the branching ratio measurements leading to H and OH products, tunable light at and around the H ( 2 S ) Lyman α transition at 121.567 nm is generated by a similar procedure to the one described above for the detection of O ( 1 D) atoms, using a fundamental wavelength of 729.4 nm and 210 Torr of Kr and 540 Torr of Ar in the tripling cell. The VUV light generated in this way is collimated by a magnesium fluoride lens at the cell exit, with the tripling cell situated perpendicular to both the observation axis and the supersonic flow. Fluorescence from cold atoms within the supersonic flow is collected perpendicular to the supersonic flow axis and the tripling cell by a solar blind photomultiplier tube (PMT) which is isolated from the flow reactor by a lithium fluoride (LiF) window. The evacuated region between the PMT and the LiF window contained an LiF lens to focus the atomic fluorescence emission onto the PMT photocathode. The PMT output is fed into a boxcar integrator coupled to a computer for signal acquisition. O ( 1 D) and H ( 2 S ) VUV-LIF signals are recorded as a function of time between photolysis and probe lasers with the timing being controlled by a delay generator. 30 laser shots are averaged at each time point, with each kinetic decay consisting of at least 40 intervals. To establish the baseline signal, several time points are recorded with the probe laser firing prior to the photolysis laser.

The gases used in the experiments (Ar 99.999%, O 2 99.999%, Xe 99.998%, CH 4 99.9995%, H 2 99.9995%) are flowed directly from cylinders into calibrated mass flow controllers, allowing us to calculate the reagent and precursor concentrations from their flow ratios using the calculated supersonic flow densities.

B. Theoretical Method

In this work, the rate coefficients of the O ( 1 D) + CH 4 reaction are computed based on a local PES in the reactant asymptotic part because of its barrierless and exothermic properties [START_REF] González | [END_REF]9,11].

In order to construct this local PES, the fundamental invariant neural network (FI-NN) fitting method [27] is used to fit a 81-15-40-1 NN function with 1911 parameters to a total of 36690 ab initio points, which were computed at the MRCI+Q/aug-cc-pVTZ level [11]. In the NN function, the activation function of x √ 1+x 2 is used. Among these ab initio points, 90% are used as the training set, while the other 10% are used as the validating set. The training and validating deviations of the NN PES are 0.208 and 0.239 meV, respectively. Since the RPMD method has been well reviewed elsewhere [12][13][14][15][16]28], only brief notes are given here. Consider an N-atomic (i.e. 3N-dimensional) system with momenta p = {p i } N i=1 , positions r = {r i } N i=1 , and potential energy

V(r 1 , • • • , r N ). The RPMD Hamiltonian is proved [12, 13] in the following (3N × n)-dimensional form H n (p, r) = n j=1        N i=1        |p ( j) i | 2 2m i + 1 2 m i ω 2 n r ( j) i -r ( j-1) i 2        + V(r ( j) 1 , • • • , r ( j) N )        , (1)
where 1), extensive classical molecular-dynamics (MD) simulations are performed to compute the RPMD rate coefficients [14][15][16], which can be written as the Bennett-Chandler factorization [START_REF] Bennett | Molecular dynamics and transition state theory: The simulation of infrequent events[END_REF],

ω n = (β n ) -1 , β n = β/n, β = 1/(k B T ), and r (0) i = r (n) i . Based on equation (
k (n) RPMD (T ) = k QT S T (s 0 , T ) • κ (n) (t p , T, s ‡ ) • exp -β∆W (n) .
Here k QT S T (s 0 , T ) is the surface area of a sphere with large enough radius times the thermally averaged flux of reactants entering this sphere, κ (n) is the transmission coefficient at the peak position of the free-energy curve, and ∆W (n) is the free-energy barrier. s is the scaled reaction coordinate, where s 0 = 0.0 and s = 1.0 represent reactant and product asymptotic groups, respectively. It is worth noting that there exists ∆W (n) for this barrierless reaction. We shall return to this point later.

As given above, only intermediate-forming part of the reaction coordinate, s ∈ [-0.05, 0.26],

is involved in computing ∆W (n) . This part of s is divided into 107 windows, called umbrellasampling geometries. Within each window, extensive RPMD simulations are performed to obtain the mean force ∂W/∂s and then ∆W (n) is determined by numerical integration. To do this, the RPMD system is first equilibrated for 20 ps, followed by a production run of 100 ps. Next, the transmission coefficient κ (n) is computed by first equilibrating for 20 ps, and then sampling the initial configurations once every 2 ps to serve as the initial positions for 100 child trajectories with t p = 7 ps. Here a total of 2.0 × 10 5 trajectories are evolved for κ (n) . All of the RPMD simulations are carried out using the modified version of the RPMDrate package [30].

III. RESULTS

A. Experimental Measurements

As there exist two main product channels for the O ( A discussion of these additional losses can be found in Grondin et al. [24]. To determine the second-order rate coefficients at a specified temperature, the excess CH 4 concentrations are varied and the resulting pseudo-first-order rate coefficients are plotted against the corresponding CH 4 concentration as shown in Figure 2. The gradient of the straight line fit to these data (weighted by the calculated uncertainties of fits to the individual decay profiles) yields the final rate coefficient.

The y-axis intercept value of these plots is large, corresponding to the relaxation of O ( 1 D) atoms mostly through collisions with the carrier gas Ar. We will discuss the temperature-dependent rate coefficients after showing the theoretical calculations. In Table I, the number of individual measurements and the range of CH 4 concentrations (in 10 14 cm -3 ) used in measurements are also given.

B. Theoretical Calculations

In order to accurately compute the overall rate coefficients, let us first consider the extended Hinshelwood type mechanism of the O ( 1 D)

+ CH 4 reaction, A+B k a ---- k -a C k r
-→ P, where A and B are reactants, C is intermediate, and P is product. Applying the steady state approximation to the intermediate complex C, it is easy to obtain the overall rate coefficient in the form k = k r k a k -a +k r = k a k -a /k r +1 , where k a can be computed by RPMD simulations based on a local PES of reactant asymptotic part and the chemical activation ratio k -a /k r can be estimated by the free-energy difference of the reactants and products, ∆ r W. Because of the large value of |∆ r W| for the title reaction (larger than 1.5 eV), k -a /k r must be very close to zero. Indeed, the RPMD simulations on the global PIP PES [11] predict that k -a /k r ∼ 0 and hence k k a , which implies that the Hinshelwood type mechanism is a good approximation. Therefore to compute k a , we construct the local PES using the NN fitting method and the previous MRCI+Q energy database [11]. As shown in Figure 3, this local NN PES shows a very low fitting error (less than 1.0 meV) over the energy range of 0.0 -1.5 eV, which implies that it is highly accurate.

Let us turn to the RPMD simulations. To determine the number of beads n in computing ∆W (n) , RPMD simulations with n = 1, n = 4, and n = 8 are performed at 75 K and 100 K. At 75 K the free-energy barriers from n = 1, n = 4, and n = 8 are 1.8, 1.5, and 1.6 meV, respectively. At 100 K the free-energy barriers from n = 1, n = 4, and n = 8 are 2.6, 2.1, and 2.5 meV, respectively. The errors of exp(-β∆W (n) (T )) associated with n are then smaller than 6%. These free energy results clearly show that one can use n = 4 to compute ∆W (n) . Similarly, n = 1 is found to be enough to compute κ (n) . The convergence of small number of beads should be caused by the small zero point energy (ZPE) differnece (3.16 kcal/mol) [11] between the reactants and the complex. Moreover, the small number of beads for both ∆W (n) and κ (n) implies that quantum effects do not strongly affect the rate coefficients, in a similar manner to previous 3D TIQM [10] simulations.

Next, after inspecting the convergence, the temperature-dependent ∆W (n) are computed and resulting free-energy curves are given in the left panel of Figure 4. Although O ( 1 D) + CH 4 is potential-energy barrierless, it can be found that the RPMD free energy curves do indeed display barriers. This is because the decrease in enthalpy from reactants to the peak of free-energy is somewhat compensated by a decrease in entropy, which clearly indicates the competition between decrease of potential energy and damping from unreactive coordinates of the system. Curves of transmission coefficients versus evolution time t are given in the right panel of Figure 4, where one can find the κ (n) curves oscillate with a frequency of approximately 10 fs. These oscillations may be associated with stretching modes (3130 and 3196 cm -1 [11]) of the C-H bonds in methane, whose periods are roughly 10 fs.

C. Rate Coefficients

The present VUV-LIF and RPMD rate coefficients, together with previous computational [START_REF] González | [END_REF][9][10][11] and experimental [1][2][3][4][5][6][START_REF] Sander | Chemical Kinetics and Photochemical Data for Use in Atmospheric Studies[END_REF] results are given in Table I. For clarity, in Figure 5 we show the rate coefficients as a function of temperature.

At 296 K the present VUV-LIF result of 1.65 ± 0.17 is close to previous experimental results [1][2][3][4] as well as the IUPAC recommended value of 1.50 ± 0.36 [5][6][START_REF] Sander | Chemical Kinetics and Photochemical Data for Use in Atmospheric Studies[END_REF]. The relative differences are smaller than 10%, which implies that the present VUV-LIF measurements are accurate. It should be noted at this point that the rate coefficient values measured in this study, along side those of Blitz et al. [2] were the only ones to be obtained using direct detection of O ( 1 D) atoms. All other previous studies used a range of indirect detection methods to follow the progress of the O ( 1 D) + CH 4 reaction. Now, let us compare the RPMD rate coefficients with the previous QCT and [START_REF] González | [END_REF] and TIQM [10] results. As shown in Table I, the previous QCT results [START_REF] González | [END_REF] at 300 K and 200 K are given to compare with the RPMD ones at 296 K and 210 K. Although there is temperature difference (smaller than 10 K), good agreement between RPMD and QCT [START_REF] González | [END_REF] can be found. The relative differences are smaller than 20%. Next, turning to the TIQM results, since TIQM calculations were performed for the OH + CH 3 channel only, they are first scaled by a branching ratio of 0.78, which was previously determined [10]. After scaling the TIQM rate coefficients, reasonably good agreement between RPMD and TIQM [10] results is found. The relative differences are smaller than 28%, except for the one at 50 K (roughly 37%). These differences may arise from the following three points. First, the reduced dimensional (3D) PES [10] was constructed by the UMP2/6-311G(2df,2pd) energy points, while the present full dimensional (12D) PES is constructed at the higher MRCI+Q/aug-cc-pVTZ level. Second, potential energy functions and fitting errors of these two PESs may influence the final rate results. Finally, the 3D TIQM calculations [10] may lose dynamics information, such as contributions of the umbrella vibration of the CH 3 fragment and symmetric/asymmetric stretching vibrations of CH 4 , even though previous work has shown these contributions to be small [10].

Finally, let us compare the present VUV-LIF and RPMD rate coefficients. At 75 K, 127 K, and 296 K, good agreement between RPMD and VUV-LIF results is found with relative differences being smaller than roughly 20%. At 50 K, however, the relative difference between RPMD and VUV-LIF results is 42%. This larger relative difference could arise from (i) contributions of excited state PESs of the O ( 1 D) + CH 4 system, (ii) an absence of long range (LR) interactions between O ( 1 D) and CH 4 , and (iii) errors in ab initio calculations and construction of the PES. Now, let us discuss these points in detail.

IV. DISCUSSIONS

It has been found that the contribution of excited states, in particular the 1 A state of O ( 1 D), may play a non-negligible role in the reaction. Through their CASPT2 and CASSCF calculations, Hernando et al. [9] expected a significant contribution of the excited 2 1 A PES to the O ( 1 D) + CH 4 reaction. Recently, two of us (K.M.H. and J.-C.L.) found opposite temperature dependence for reaction over the 1 A and 1 A PESs of the C ( 1 D) + H 2 system [19]. The RPMD rate for the 1 1 A ( 1 A ) PES was found to increase with temperature, while that for 2 1 A ( 1 A ) decreases with temperature. The net effect was seen to be a very slight increase in the total rate coefficient [19]. On the other hand, the rather small RPMD κ (n) values might cause the rather small rate coefficient at 50 K. As shown in the right panel of Figure 4, at T > 75 K the κ (n) values in the plateau region decrease with increasing temperature. However, the situation is changed at 50 K, where the κ (n) (50 K) value is smaller than expected. Adding a perturbated LR attractive potential, we might expect that the value of κ (n) (50 K) increases. Consequently, the computed rate coefficient at 50 K increases.

Now let us turn to the LR potential itself. Adding the LR interaction model of α

1 R -3 -α 2 R -4
into the PES of the D + + H 2 system, Lara et al. [32] reported good agreement between theoretical and experimental temperature-independent rate coefficients at 10 K -90 K [32]. In their LR potential model, the dominant contributions involve the charge-quadrupole interaction R -3 and the charge-induced dipole R -4 . Obviously, this LR potential is too attractive for the present radicalmolecule system. We might expect a ∝ R -6 LR model and hence the rate coefficients more slowly approach to zero as T → 0 than those shown in Figure 5. However, the application of the LR potential model for the O ( 1 D) + CH 4 system is still an open question.

Finally, the errors in constructing the PES may be another reason. In the present work, a set of MRCI+Q/aug-cc-pVTZ energy points [11] are used to construct the present NN PES. Although the systematic and random errors are small as shown in Figure 3, errors may still be introduced through the initial MRCI+Q calculations and constructions. At temperatures as low as 50 K, since the free-energy barriers are smaller than 1.5 meV, these PES errors may influence the resulting free-energy barrier and hence the rate coefficients. However, it will be difficult to further improve the present RPMD results because the present PES error of roughly 1 meV is already very small.

V. CONCLUSIONS

In summary, we determine rate coefficients for the barrierless O ( 1 D) + CH 4 reaction at low temperatures (50 K -296 K) both theoretically and experimentally. On the basis of the present local and accurate NN PES, extensive RPMD calculations are performed, predicting similar rate coefficients and a similar temperature-dependence to the experimental VUV-LIF rate coefficients.

Moreover, the fact that a small number of beads can approach convergence clearly implies that quantum effects do not strongly affect the rate coefficient values, a result which was also found by previous TIQM calculations [10]. Experimentally, the rate coefficients at 50 K, 75 K, 127 K, and 296 K are measured using a supersonic flow reactor, where the pulsed laser photolysis of O 3 molecules is used as the source of O ( 1 D) atoms. During the measurements, O ( 1 D) is detected directly by VUV-LIF at 115 nm. A comparison of the present VUV-LIF rate coefficients with the RPMD ones, leads to rather large difference between the experimental and theoretical rate coefficients at 50 K. This may be due to several reasons, including possible contributions of excited state surfaces and the LR interaction. Furthermore, errors in the ab initio calculations and/or construction of the PES may also influence the resulting rate coefficients. Therefore, further work is required to improve the calculated rate coefficients at T < 75 K. However, since the atmosphere temperature in the Earth is in the range of 100 K -300 K [25], the present calculations and experiments provide very useful references for atmospheric chemistry. g The QCT values of González et al. [START_REF] González | [END_REF] at 200 K and 300 K are shown. RPMD and VUV-LIF rate coefficients, respectively. The other symbols with error bars repre-sent previous experimental results [1][2][3][4][5][6][START_REF] Sander | Chemical Kinetics and Photochemical Data for Use in Atmospheric Studies[END_REF]. Finally, the maroon "×" and the magenta "A" symbols represent previous QCT [START_REF] González | [END_REF] and TIQM [10] results, respectively.

  1 D) + CH 4 reaction, leading to OH + CH 3 and H + CH 2 OH/CH 3 O products, it is important to quantify the contribution of these pathways. By comparing the amplitude of the H atom signal from the O ( 1 D) + CH 4 reaction versus the one obtained from a reference process, the O ( 1 D) + H 2 reaction recorded at the same time, (assuming that this reaction produces atomic hydrogen with a 100% yield) it was possible to determine the sum of the H atom branching channels at a given temperature. In the left panel of Figure 1 we show the measured H atom fluorescence intensity as a function of time for both the target O ( 1 D) + CH 4 and reference O ( 1 D) + H 2 reactions at 50 K. At least seven pairs of decays similar to the ones shown are recorded at each temperature, to minimize potential experimental errors. At 50 K and 296 K, mean values of 0.20 ± 0.01 and 0.21 ± 0.01 respectively are obtained with these statistical errors cited at the level of one standard deviation. This result is in excellent agreement with the previous crossed molecular beam experimental determination of the OH + CH 3 productchannel branching ratio of 0.77 at higher energy[31]. Moreover, on the basis of the PIP PES, QCT simulations[11] also predicted a branching ratio of 0.8 for this channel.Example traces of the O ( 1 D) VUV-LIF signal intensity as a function of time recorded at 50 K with the carrier gas Ar alone, and with added CH 4 are shown in the right panel of Figure1. It can be seen that O ( 1 D) is more efficiently removed when CH 4 is present in the supersonic flow, although collisions with Ar also lead to its rapid loss through quenching to the ground state O ( 3 P). Nevertheless, [Ar] is fixed for any series of measurements, so that this process represents a constant pseudo-first order loss for O ( 1 D). Methane is added as the excess reagent, so that simple exponential fits to these temporal profiles allowed us to determine the total pseudo-first-order rate coefficients for O ( 1 D) loss from the time constant. Other potential loss processes for O ( 1 D) atoms including quenching by O 2 and O 3 , and loss through diffusion are all negligibly small when compared to the quenching induced by the carrier gas Ar and by the reactive losses with CH 4 .

The 2 1 A

 1 ( 1 A ) PES of the O ( 1 D) + CH 4 reaction could play the same role as in the case of C ( 1 D) + H 2 . However, this is still an open question since the excited-state PESs of the O ( 1 D) + CH 4 system are not currently available.

• s - 1 )

 1 of the O ( 1 D) + CH 4 reaction, together with the RPMD ones. The first column gives the temperature in Kelvin. The second column gives number of individual measurements. The third column gives the range of concentrations (in 10 14 cm -3) of CH 4 used in the measurements. The rightmost two columns give the relative differences (in percentage) between (i) the present VUV-LIF and RPMD results, and (ii) the present RPMD and previous 3D TIQM[10] results. The other columns present the rate coefficients. For clarity, the rate coefficients as a function of temperature are shown in Figure5.

Figure 1 :

 1 Figure 1: VUV-LIF intensity as a function of time at 50 K for H ( 2 S ) atom formation (left panel) and decay of O ( 1 D) atoms (right panel). In the left panel, the blue open circles represent the Hatom signal from the O ( 1 D) + CH 4 reaction with [CH 4 ] = 5.7 × 10 14 cm -3 , while red solid circles represent the H-atom signal from the O ( 1 D) + H 2 reference reaction with [H 2 ] = 5.7 × 10 14 cm -3 . In the right panel, the blue open circles represent the O-atom signal with no added methane, where O ( 1 D) decay is almost entirely due to to quenching with the carrier gas Ar. The red filled circles represent the O-atom signal with [CH 4 ] = 8.4 × 10 14 cm -3 . In both sub-figures, the blue (or red) solid lines are exponential fits to the corresponding data points.

Figure 2 :

 2 Figure 2: Pseudo-first-order rate constants for the O ( 1 D) + CH 4 reaction as a function of the methane concentration. A weighted linear least squares fit yields the second-order rate coefficient. Red open circles represent VUV-LIF measurements of the rate coefficient at 296 K, while blue filled circles represent those recorded at 50 K. The error bars reflect the statistical uncertainties at the level of a single standard deviation obtained by fitting to temporal profiles such as those shown in Figure 1.

Figure 3 :

 3 Figure 3: The NN Fitting errors of energy points as a function of the corresponding ab initio potential energies. These energy points are used to construct the present local NN PES for the O ( 1 D) + CH 4 asymptotic part. All energies are relative to the reactant asymptote.

Figure 4 :

 4 Figure 4: Comparison of the RPMD free-energy curves (left panel) and transmission coefficient curves (right panel) for O ( 1 D) + CH 4 at 50 K (black), 75 K (red), 100 K (green), 127 K (blue), 175 K (yellow), 210 K (brown), and 296 K (violet). The free energy in the left panel is drawn as a function of scaled reaction coordinate s, while the transmission coefficient in the right panel is drawn as a function of evolution time (in fs) of child trajectory.

Figure 5 :

 5 Figure 5: Experimental and theoretical rate coefficients (in cm 3 • s -1 ) of the O ( 1 D) + CH 4 reaction as a function of the temperature (in Kelvin). The black "+" and red " * " symbols represent present

  FIG. 1: RPMD, PIP PES, O + CH4

TABLE I :

 I Summary of the present VUV-LIF rate coefficients (given in cm 3

  Uncertainties on the temperatures represent the statistical (1σ) errors obtained from Pitot tube measurements of the impact pressure. Uncertainties on the measured rate coefficients represent the combined statistical (1σ) and estimated systematic (10%) errors.

	-21.8%	-19.1%	-22.5%	-28.1%	-36.8%
	-	5.6%	-	-20.3%	-42.1%
	f	f	f	f	f
	10 -10	10 -10	10 -10	10 -10	10 -10
	2.06 ×	1.88 ×	1.73 ×	1.53 ×	1.25 ×
		10 -10	10 -10	10 -10	10 -11
		1.52 ×	1.34 ×	1.10 ×	7.99 ×
		10 -10		10 -10	10 -10
		0.15) ×		0.16) ×	0.15) ×
		11 (1.44 ±	-	6 (1.38 ±	8 (1.38 ±
		0 ∼	-	0 ∼	0 ∼
		127 ± 2 34	100 -	75 ± 2 42	50 ± 1 41

a b Present work. c IUPAC recommended value [5-7]. d Chemiluminescence measurements of Vranckx et al. [4]. e PLP-PLIF value by Dillon et al. [3].

f TIQM values by Ben Bouchrit et al.

[10]

.

Acknowledgments K.M.H. and J.-C.L. thank the CNRS programs PCMI and PNP for their support. Q.M., K.S., and D.H.Z. thank Ministry of Science and Technology of P. R. China (Grant No. 2013CB834601), National Natural Science Foundation of China (Grant Nos. 21433009 and 21503214), and Chinese Academy of Sciences (Grant No. XDB17010200). Q.M. also acknowledges the start-up fund from Northwestern Polytechnical University. Q.M. wishes to thank Dr. B. Fu for discussing on the dynamics of the O ( 1 D) + CH 4 system. We thank Prof. Dr. P. Honvault for providing their 3D