A N Artemyev
email: artemyev@physi.uni-heidelberg.de

A Surzhykov

P Indelicato

G Plunien

Th Stöhlker

Finite basis set approach to two-center Dirac problem in Cassini coordinates

Keywords: numbers: 34, 10, +x, 34, 50, -s

We lay out a new approach to solving the two-center Dirac eigenvalue problem. It is based on application of B-spline basis sets constructed in Cassini coordinates. The approach provides a very promising way for the theoretical description of various atomic processes such as charge transfer, excitation, ionization and electron-positron pair production, that accompany slow collisions of heavy ions. Moreover, yielding a (quasi) complete set of eigensolutions, the finite-basis-set method allows a systematic analysis of the quantum electrodynamics (QED) corrections to the energy levels of heavy quasi-molecules formed in these collisions.

Introduction

Owing to recent advances in accelerator and storage ring techniques, new possibilities arise today to study formation of quasi-molecules in (relatively) slow heavy-ion collisions. In contrast to experimental studies performed during the last century, ion storage and cooler rings do allow nowadays to perform such experiments at arbitrary low beam energies for fully stripped or H-like ions interacting with high-Z neutral targets (Verma et al. 2005, Verma et al. 2006). For the future one may even anticipate slow collisions between high-Z bare and one-electron ions (FAIR 2001, Gumberidze et al. 2005, Stöhlker et al. 2007, Gumberidze et al. 2009). These quasi-molecules are known to be a versatile tool for investigating a number of fundamental problems in heavy-ion atomic physics (Betz et al. 1976, Müller 1976, Soff et al. 1977, Rafelski et al. 1978, Soff, Müller & Greiner 1978, Kirsch et al. 1978, Reinhardt et al. 1979, Soff et al. 1979, Müller et al. 1983, Soff et al. 1987, Müller-Nehler & Soff 1994). In particular, analysis of molecular spectra can provide unique knowledge on the relativistic, many-body and quantum electrodynamics (QED) effects in the non-perturbative domain of high nuclear charges and super strong electromagnetic fields.

Theoretical understanding of the critical-field phenomena in heavy quasi-molecules requires, in general, solution of the two-center time-dependent Dirac equation. For low velocities of colliding ions this equation may still be treated adiabatically and, hence, can be traced back to the static (two-center) eigenproblem:

Ĥ2C Ψ = Ψ , ( 1
)
where Ψ is the electronic wavefunction and the Hamiltonian Ĥ2C is given, for the simplest case of the one-electron-two-nuclei system, by: Ĥ2C = Ĥfree + V ( r, R) .

(2)

Here, Ĥfree is the free electron Dirac Hamiltonian and V describes Coulomb interaction between an electron and two nuclei with R = R 1 -R 2 being the inter-nuclear vector.

A large number of methods have been proposed during the last two decades for solving the eigenproblem (1) and, thus, for performing propagation of the electron wave packet in (slowly varying) nuclear potential. Most of these methods are based on the well-elaborated linear-combination-of-atomic-orbitals (LCAO) approach (Sepp et al. 1986, Toshima & Eichler 1988, Momberger et al. 1993, Rumrich et al. 1993[START_REF] Eichler | Relativistic Atomic Collisions Academic Press[END_REF], Schulze et al. 1998, Gail et al. 2003, Tupitsyn et al. 2010). Within the LCAO approach, the (quasi-) molecular wavefunctions Ψ are constructed from two sets of atomic orbitals centered on two nuclei. Despite the relatively simple treatment of dimer quasi-molecules at large inter-nuclear distances, the LCAO method may encounter computational difficulties if ions are placed close to each other and the orbitals become much more molecular than atomic.

For rather small inter-nuclear distances, the two-center Dirac problem (1)-( 2) can be solved efficiently in prolate spheroidal coordinates as discussed, for example, in Ref. (Kullie & Kolb 2001). Unlike LCAO, the elliptic coordinates, being essentially two-center, are very practical for the accurate description of molecular terms. In contrast, they do not fit very well for the description of the "atomic-type" wavefunctions. Therefore, to the best of our knowledge, the application of the prolate spheroidal coordinates becomes less efficient with the increase of the distance between ions.

Besides the application of the elliptic coordinates, the structure and dynamics of heavy quasi-molecules at small inter-nuclear distances can be relatively well understood within the so-called monopole approximation (Müller-Nehler & Soff 1994, Ackad & Horbatsch 2008). Such an approximation, that employs solely the spherically symmetrical part of the two-center potential, allows one to reduce the problem to the spherically-symmetric one, with a total nuclear charge distributed homogeneously on the spherical shell with the diameter equal to the internuclear distance. This method was found to be very useful for description of strong-field phenomena in close ion collisions but, again, its accuracy is known to be poor when Dirac centers appear to be far from each other.

At the end of our brief review we should mention also the direct numerical integration of the two-center Dirac equation. Towards this line, the various lattice methods have been developed over the last decades (Becker et al. 1986, Strayer et al. 1990, Thiel et al. 1992[START_REF] Momberger | [END_REF], Wells et al. 1996, Ionescu & Belkacem 1999, Busic et al. 2004, S ¸engül et al. 2009). Being potentially very precise these methods are, however, rather time-consuming, which puts restrictions on their broad application.

As seen from the discussion above, many of the methods available nowadays for solving the eigenproblem (1) and performing the time-propagation of electron (or positron) wave packets have restricted ranges of application. Moreover, the proper description of the positive-as well as the negative-continuum solutions of the (twocenter) Dirac equation remains also a rather difficult task. In this contribution, we present an efficient theoretical approach for the treatment of the relativistic two-center problem at arbitrary inter-nuclear distance. Within this approach we employ the finite (discrete) basis solutions, as constructed in Cassini coordinates from the B-spline sets, in order to solve Eq. ( 1) and to find a (quasi) complete set of eigenvalues and eigenfunctions. On the basis of such a set that includes not only bound-but also (positive and negative) continuum-solutions of Eq. ( 1), further analysis can be performed both, for the quantum electrodynamics (QED) corrections to energy levels of heavy quasi-molecules and for the dynamical behavior of these molecules in (slow) ion-ion collisions.

The paper is organized as follows. In the next section we shall provide the detailed description of the method. In particular we will discuss the choice of the coordinates, describe the finite basis set approach and introduce the basis functions. The details of computations are given then in Sec. 3. In order to illustrate the application of the method, detailed calculations are performed for the energy spectra of H + 2 , Th 179+ 2 and U 183+ 2 quasi-molecules. In Sec. 4 we will discuss the results of these calculations and will

0 a -a r 2 r 1 Z θ 2 θ 1 R 1 R 2 Figure 1.
The explanation of the Cassini coordinates (left) and the lines with constant w and δ variables (right). On the right picture the solid lines are with constant value of w and the dashed ones -with constant δ.

compare them with the previously reported data. Finally, short summary and outlook will be given in Sec. 5.

The natural system of units (c = m e = h = 1; the energy unit is the electron rest mass m e c 2 ) is used throughout the paper.

Theoretical background

Cassini coordinates

As mentioned already in the Introduction, the proper choice of the coordinate system may significantly simplify theoretical treatment of the two-center Dirac problem. In our analysis below, we shall solve the Dirac equation in the so-called Cassini coordinates. Since the properties of the Dirac equation in this orthogonal coordinate system have been investigated in Refs. (Schlüter et al. 1983, Wietschorke et al. 1983), here we may restrict ourselves to a rather short account of the basic relations. The Cassini coordinates w, δ, and φ of the electron relative to the nuclei are defined as follows:

w = √ r 1 r 2 a , (3)
δ = θ 1 + θ 2 2 , ( 4
) φ = φ , (5)
where r i is the distance between the electron and the i-th nucleus and a

= | R 1 -R 2 |/2
is the half of internuclear distance (see Fig. 1). Moreover, the angle between the internuclear axis (Z-axis) and the position vector of the electron relative to the i-th nucleus is denoted as θ i while, finally, φ is the azimuthal angle.

The Cassini coordinates introduced in Eqs.

(3)-( 5) are known to obey several important properties. For example, the curves corresponding to constant values of the coordinate w form closed lines around the first (0 < δ < π/2) or the second (π/2 < δ < π) center if w < 1. In contrast, for the values of w greater than one, the (w-) equidistant curves surround both centers as it can be seen from the right panel of Fig. 1. Finally, the line w = 1 is lemniscate.

Another important feature of Cassini coordinates concerns the behavior of the (twocenter) Coulomb potential along w-equidistant lines. For instance, as it has been shown in Ref. (Wietschorke et al. 1983) for the particular case of two point-like nuclei with charges Z 1 = Z 2 , the electrostatic potential as calculated along the curve w = const deviates from its mean value by less than 1.5 % for w < 0.25 and w > 4 and by about 17.2 % in the worst case of w = 1. Such a behavior selects out Cassini coordinates as the most adequate choice of coordinates for the solution of the two-center problem, since the curves w = const almost coincide with the equipotential lines.

After having discussed the basic properties of Cassini coordinates, we shall employ them to re-write the Dirac Hamiltonian (2) in the form:

Ĥ2C = -i D 1/4 aw α 3 ∂ ∂w + α 1 1 w ∂ ∂δ -iα 2 1 ρ ∂ ∂φ + V (w, δ; a) + βm , (6)
where α i and β are the Dirac matrices and V is the two-center Coulomb potential. Explicit form of this potential depends on the particular choice of nuclear charge distribution and will be discussed later in Section 3. In Eq. ( 6) ρ is the distance between the electron and the internuclear axis:

ρ = a √ 2 √ D -1 -w 2 cos 2δ , ( 7
)
and D is the discriminant

D = w 4 + 2w 2 cos 2δ + 1 . ( 8
)
Because of the axial symmetry of the problem, it is convenient to seek the eigenfunctions of the Hamiltonian (6) in the form

Ψ(w, δ, φ) = ψ µ (w, δ) exp(iµφ) , (9)
with µ = ±1/2, ±3/2, ±5/2, ... being the projection of the total angular momentum of electron onto the internuclear axis. For such a choice of the total wavefunction, the Dirac Hamiltonian Ĥ2C reduces to

Ĥ(µ) 2C = -i D 1/4 aw α 3 ∂ ∂w + α 1 1 w ∂ ∂δ + α 2 µ ρ + V (w, δ; a) + βm , (10)
and the problem becomes two-dimensional. The solution of this eigenproblem is, however, not a simple task owing to the fact that the Hamiltonian (10) is complex if the standard representation of Dirac matrices is used. It is often more convenient, therefore, to use an alternative representation in which matrices α 1 and α 3 are purely imaginary while α 2 is real. This can be easily achieved if one defines (new) Dirac matrices as (Schlüter et al. 1983):

α i = γ 0 γ i , (11)
where

γ 0 = σ 3 ⊗ I , (12)
γ 1 = iσ 2 ⊗ σ 2 , (13)
γ 2 = iσ 2 ⊗ σ 1 , (14)
γ 3 = -iσ 1 ⊗ I , (15)
with σ i and I being the Pauli and unity matrices, correspondingly.

Finite basis set method

Equation ( 10) displays the reduced Dirac Hamiltonian that describes-in Cassini coordinates-the motion of a single electron in the field of two nuclei. So far, a few approaches have been proposed to find the eigenvalues and eigenfunctions of such a Hamiltonian. In Refs. (Schlüter et al. 1983, Wietschorke et al. 1983), for example, authors made use of the semi-discrete Galerkin method in order to solve the relativistic two-center eigenproblem. An alternative route can be developed based on the finite basis set method. Although both, the Galerkin and the finite-basis-set methods are variational ones, the latter approach seems to be preferable since it provides a (quasi-) complete set of eigenvectors. The knowledge of such a set is substantial for the theoretical analysis of the dynamical properties of heavy quasi-molecules. Moreover, the summation over the complete spectrum is required for the computation of the radiative corrections to the molecular energy levels. In this Section, therefore, we will briefly recall the main ideas of the finite basis set method and will discuss how this method is applied to find (a complete set of) eigenfunctions and eigenvalues of the Dirac Hamiltonian (10).

In order to construct the finite basis sets, we shall start from the principle of least action:

δS = 0 , (16)
from which the Dirac equation can be derived. In this expression, the action S is defined as:

S = ψ µ | Ĥµ 2C -ε|ψ µ , ( 17
)
where ε is the Lagrange multiplier, introduced to ensure the normalization condition, which can be written in Cassini coordinates as

ψ µ |ψ µ = dw dδ aw D 1/4 ψ † µ ψ µ = 1 . (18)
Here, the wave function ψ µ is the standard Dirac four-spinor, whose (four) components ψ i µ , i = 1, ..4 we write as a finite expansion

ψ (i) µ (w, δ) = N α=1 C (i) µ,α B α (w, δ) , (19)
over some basis functions B α (w, δ). The explicit form of these functions is not crucial for the following discussion and will be specified later in Section 2.3. Mathematically, the implication of Eq. ( 19) is equivalent to the projection of the whole space of wave functions ψ µ onto the subspace, which basis is given by a set {B α }.

By inserting now the wave function Eq. ( 19) into the least action principle ( 16) and by evaluating the variation δS with respect to change of expansion coefficients C (i) µ,α , we obtain the system of differential equations:

(δS) j = ∂S ∂( C) j = 0 , (20)
where, for the sake of shortness, we use notation ( C) j ≡ C (i) µ,α . From this system immediately follows the matrix equation to determine the vector

C = {( C) 1 , ...( C) 4N }: A C = εK C , (21)
where A and K are 4N × 4N matrices, constructed from N × N blocks:

A =       (V + mB) 0 W M-∆ 0 ( V + mB) M + ∆ W -W M-∆ (V -mB) 0 M + ∆ -W 0 ( V -mB)       , (22)
K =       B 0 0 0 0 B 0 0 0 0 B 0 0 0 0 B       . (23)
Here, elements of the matrices B, V , M, W , and ∆ are given in terms of the integrals involving the basis functions B i :

B i,k = dwdδ aw D 1/4 B i B k , (24)
V i,k = dwdδ aw D 1/4 B i B k V (w, δ) , (25)
M i,k = dwdδ aw D 1/4 B i B k µ ρ , (26)
W i,k = dwdδ B i ∂B k ∂w , (27)
∆ i,k = dwdδ B i w ∂B k ∂δ . (28)
Equation ( 21) is known as a generalized eigenvalue problem that can be solved by employing the standard techniques from linear algebra. In the present work, for example, we have used the well-established Linear Algebra PACKage Ref. (LAPACK 2008), which provides fast and effective routines for solving systems of linear equations.

Basis functions

By solving the generalized eigenvalue problem one obtains both, the energies µ of the molecular levels and the (sets of) expansion coefficients C (i) µ,α . Thus, any theoretical analysis of the structure and dynamics of "one-electron-two-nucleus" system can be traced back to the system of equations ( 21). As seen from Eqs. ( 22)-( 28), however, in order to define such a system one has to agree first about the explicit form of the basis functions B α . In the present work we will construct these functions from the Bspline sets. During the last decades, the B-spline sets have been already successfully applied for the treatment of the one-center Dirac (Johnson et al. 1988[START_REF] Fischer | [END_REF] as well as the two-center Schrödinger problems (Artemyev et al. 2004). In these works it has been argued that the methods based on B-splines provide a quasi-complete set of solutions of the Dirac (or Schrödinger) equation, which allows perturbative calculations, implying summation over whole (or part) of the spectrum.

The properties of B-splines are described in details in Ref. [START_REF] De Boor | A practical guide to splines[END_REF]. Here we like to recall briefly their main features that are crucial for an accurate solution of the Dirac problem. The B-splines of the k-th order are defined by the recurrence relation

B k i (x) = x -t i t i+k-1 -t i B k-1 i (x) + t i+k -x t i+k -t i+1 B k-1 i+1 (x) , (29)
where the B-splines of the first order read as

B 1 i (x) = 1, t i ≤ x < t i+1 0, otherwise . (30)
In these expressions, the number of arbitrary positioned knots t i is by k larger than the number of splines. Moreover, the first and the last k knots must be equal.

Before we discuss the construction of the basis functions B α , let us first recall the analytical properties of the B-splines. Namely, it follows from Eqs. ( 29)-( 30) that: (i) B-splines of k-th order are piecewise polynomials of (k -1) order, (ii) the i-th B-spline of k-th order is defined on the interval from t i to t i+k , (iii) the splines together with their (k -1) derivatives are continuous, and (iv) the sum of all B-splines of given order at any point is equal to 1. These properties make the B-splines remarkably convenient for the numerical approximation of the basis functions from Eq. ( 19), which can be written as:

B α = f (w, δ) B nw l(α) (w) B n δ m(α) (δ) . (31)
Here B nw l (w) and B n δ l (δ) are two independent sets of B-splines. The function f has been chosen to describe the behavior of the Dirac wavefunction in regions, where the polynomial approximation fails. For an accurate treatment of two-center Dirac problem this function reads:

f (w, δ) = ρ √ D , (32)
where ρ and D are given by Eqs. ( 7) and (8), correspondingly. Finally, for evaluation of Eq. ( 31) we have used equidistant knots for constructing the δ-splines and an exponential grid (t i+1 = λt i , λ > 0) for the w-splines.

After the description of the finite-basis-set approach we are ready now to discuss its numerical implication. In the following section, therefore, we shall present details of calculations leading to an accurate determination of energy spectrum and (a set of) eigenfunctions of heavy quasi-molecules.

Details of calculations

To solve the generalized eigenvalue problem (21) and, hence, to find energy values and functions of the two-center Dirac Hamiltonian, we first need to compute the matrix elements ( 24)-( 28), which involve the basis functions B α . In the present work such a computation has been performed numerically by making use of the Gauss-Legendre quadratures. While the numerical integration can be performed rather easily in almost the entire space, it becomes a complicated task in the vicinity of the middle of the internuclear axis (w = 1 , δ = π/2). The power series expansion of the distance between an electron and internuclear axis ρ converges badly at this point leading to a lack of accuracy. However, noticing that √ D is the generating function of the Gegenbauer polynomials [START_REF] Erdélyi | Higher transcendental functions[END_REF]) and taking into account Eq. ( 7), there is the possibility for accurate evaluation of ρ by making use of expansions:

ρ =        aw 2 √ 2 ∞ n=0 C -1/2 n+2 (-cos 2δ)w 2n w < 1 a √ 2 (w 2 -1)(1 -cos 2δ) + w 2 ∞ n=2 C -1/2 n (-cos 2δ)w -2n w > 1 , ( 33
)
where C λ n (x) are Gegenbauer polynomials defined by

C λ 0 (x) = 1, (34)
C λ 1 (x) = 2λx , (35)
C λ (n+1) (x) = 2(n + λ) x C λ (n) (x) -(n + 2λ -1)C λ (n-1) (x) n + 1 . ( 36
)
Here we would like to point to a misprint in the expansion for ρ as given by Eq. ( 33) in Ref. (Wietschorke et al. 1983).

As seen from Eq. ( 25), apart from the basis functions B α , the knowledge on the two-center potential V (w, δ; a) is also required for the integration of the generalized eigenproblem. Explicit form of this operator, that enters the Dirac Hamiltonian (10), depends on the nuclear charge distribution employed. For example, the motion of the electron in the field of two point-like nuclei with the same charge number Z 1 = Z 2 = Z is defined by the potential:

V (w, δ; a) = -αZ √ 2 aw 2 √ D + 1 + w 2 . ( 37
)
More complicated expression is required to account for the finite nuclear size effects or in the case of hetero-nuclear problem. To derive this, we should start from the general definition of two-center potential:

V = V 1 (| r -R 1 |) + V 2 (| r -R 2 |) , (38)
where V i is the potential, generated by i-th nucleus. As usual in atomic structure calculations, such a single-center potential can obtained as

V i (r) = ρ i (r 1 ) r > r 2 1 dr 1 , (39)
with ρ i being the nuclear charge density distribution. By employing the various models of the charge distribution ρ i one can explore the isotope effects on the structure and dynamics of heavy quasi-molecules. In the present work, for example, we shall use the well-known Fermi distribution:

ρ i (r) = N 1 + exp[(r -c)/a] , ( 40
)
where N is normalization factor and the parameters a and c can be expressed in terms of the nuclear root-mean square radius (see for example (Shabaev 1993)). Moreover, in Eqs. ( 38) and ( 39), the distance between the i-th nucleus and the electron is given in Cassini coordinates by:

| r -R 1,2 | = (z ± a) 2 + ρ 2 (41) z = a 2 w 2 sin 2δ 2ρ , ( 42
)
where the plus sign refers to | r -R 1 | and the minus sign to | r -R 2 |.

Results and discussions

Spectra of heavy ions

Before starting the discussion on the predictions of the finite basis set approach for heavy quasi-molecules, let us first demonstrate the validity of the method if applied to much simpler, one-center Dirac problem. Reduction to such a (purely atomic) problem can be easily achieved if one takes one of the nuclear charges in Eq. ( 38) to be zero.

The second nucleus with charge Z = 0 can be treated as point-like which bound-state energies given by a well-known expression:

ε =   1 + (αZ) 2 (j + 1/2) 2 -(αZ) 2 + n -j -1/2 2    -1/2 , ( 43
)
where n is the principal quantum number and j is the total momentum of electron.

To perform an extensive check of our theoretical approach, detailed calculations of the energy values of hydrogen-like ions were carried out in the Cassini coordinates for various values of inter-center distance a from 10 -4 to 20 a.u. and were compared against the predictions of Eq. ( 43). Special attention in these calculations has been paid to heavy nuclei with charge 92 ≤ Z ≤ 100. For this high-Z regime and for relatively large (a > 2 a.u.) as well as small (a < 5×10 -3 ) inter-center distances, we were able to reproduce the energies of the first 20 bound states with a relative accuracy of about 10 -9 by employing in Eq. ( 31) only 40 w-splines of fifth order and 20 δ-splines of fourth order. In contrast, the number of B-splines had to be increased up to 45 and 28, respectively, in order to achieve a slightly worse accuracy of 10 -7 for the intermediate region with a ≈ 10 -2 a.u., .

While performing calculations for the energy values and eigenfunctions of the onecenter Dirac Hamiltonian we were faced with the well known problem of spurious states, i.e., non-physical solutions of Eq. ( 21). Being the general problem of the finite-basis-set methods applied to Dirac equation, the spurious states (as well as the ways for their elimination) have been discussed in the literature [START_REF] Drake | [END_REF], Johnson et al. 1988, Shabaev et al. 2004, Igarashi 2006, Igarashi 2007). In Ref. (Shabaev et al. 2004), for example, authors have argued that application of the identical basis function sets for description of the large and small components of relativistic wavefunctions leads to the appearance of the non-physical solutions of the (one-center) Dirac-Coulomb problem. It has been also demonstrated, that one can get rid of these spurious states by using the so-called dual kinetically balanced basis (DKB) set constructed from functions that provide correct relation between the large and small components of the wavefunction in the non-relativistic limit. However, such a DKB approach is most efficient when applied to the spherically symmetric Dirac problem, which can be easily reduced to a one-dimensional radial equation and two-component wave function. In the case of a two-center potential, the construction of the kinetically balanced set is not a simple task. The studies towards the extensions of DKB method to the two-center Dirac problem are currently underway and will be presented elsewhere.

Apart from applying the DKB approach, one can eliminate the spurious, nonphysical solutions of Eq. ( 21) if analyzes their oscillatory behavior. Namely, the spurious eigenfunctions with energies in the range -mc 2 < < +mc 2 differ from the physical bound solutions by significantly faster oscillations. Since all the basis functions (31) are positive-defined, the only possibility for the function ( 19) to oscillate fast, is to have alternating-sign coefficients. Therefore, for bound state solutions one can easily determine the spurious states by evaluating the value

S i µ = α:δ<π/2;C (i) µ,α =0 |C (i) µ,α | C (i) µ,α , (44)
where the summation runs over all non-zero expansion coefficients C (i) µ,α from Eq. ( 19), corresponding to basis functions defined at δ < π/2 ‡. For basis sets constructed from 800 to 2500 functions the (absolute value of) coefficient S i µ was found to be less than 15 if ψ (i) µ is spurious state and greater than 100 for the physical solutions.

Quasi-molecular spectra

So far, we have discussed the application of the finite basis set approach to the onecenter Dirac problem. Calculations carried out for such an "atomic scenario" were used by us for estimating the inherent inaccuracies, associated with the method. Now we are ready to perform the detailed analysis of the energy spectra of (quasi-) molecules.

Up to date, several attempts have been made to solve the two-center Dirac equation and to determine the binding energies of homonuclear dimers. See e.g. (Schulze et al. 1998, Kullie & Kolb 2001, Tupitsyn et al. 2010). In Ref. (Kullie & Kolb 2001), for ‡ The summation over the entire space (0 < δ < π) in Eq. ( 44) in homonuclear case will lead to vanishing of S i µ also for all physical states, which are odd with respect to the change δ → π -δ. molecules. The calculations have been done at the "chemical" internuclear distance of R = 2/Z a.u. and by employing the point-like nuclear charge distributions. Predictions, obtained in Ref. (Kullie & Kolb 2001) for such a model case, are presented in Table 1 and are compared with the solutions of the generalized eigenvalue problem (21) for the "point-like" potential (37). As seen from the table, the ground-state energies of both molecules are well reproduced (up to 7 decimal digits) within the final basis set approach. Apart from the H + 2 and Th 179+ 2 dimers, we have also explored the spectrum of the U 183+ 2 molecule. Perfect agreement between our prediction and data reported recently in Ref. (Tupitsyn et al. 2010) has been again found for the energy 1σg (see right column of Table 1).

While over the past decades the spectra of heavy quasi-molecules have been studied rather intensively for the point-like nuclei, not much is known about the nuclear-size corrections to the molecular levels. In order to estimate these corrections we made use of Eqs. ( 38)-( 40) with the proper values of root-mean square radii from (Angeli 2004, Kozhedub et al. 2008) and evaluated the ground-state energies of the H + 2 , Th 179+ 2 and U 183+ 2 dimers. As can be expected, we found that while for low-Z molecules the effect of finite nuclear size is almost negligible, it becomes sizable in the high-Z domain. For example, as seen from Table 1, the absolute value of the 1σg binding energy of U 183+ 2 quasi-molecule reduces by about 206 eV (7.5 a.u.) if Fermi-type nuclear charge distribution is taken into account.

Up to here, we have discussed the (quasi-) molecular energies as well as the nuclearsize effects for a single inter-nuclear distance R = 2/Z a.u. In future experiments, however, formation of heavy quasi-molecules is to be observed in (slow) collisions of high-Z ions, the distance between which will vary with time: R = R(t). To understand better the R-dependence of the molecular spectra, detailed calculations for the energies of the ground 1σ g and the first excited 1σ u states of the U 183+ 2 dimer have been performed for a wide range of internuclear distances. Results of these calculations are presented in Fig. 2 for both, point-like (dashed line) and extended (solid line) nuclei. Moreover, we plot also the energy of the 1s state, calculated in monopole approximation (dashdotted line) which corresponds to the solution of the Dirac equation with the spherically symmetric part of the two-center potential. As seen from the figure, the 1σ g and 1σ u binding energies tend-for the large inter-nuclear distances-to the energy of the 1s ground state of hydrogen-like uranium U 91+ . Such a well-understood behavior of molecular terms is not reproduced, however, by the monopole approximation for which the energy approaches the value of = mc 2 if R → ∞. As the nuclei come closer to each other, the difference between the energy of even (also known as "gerade") and odd ("ungerade") symmetry states becomes noticeable. Again, this splitting is well described by our two-center Dirac calculations but can not be obtained within the (spherical) monopole model. Finally, at the very small, so-called critical, internuclear distances the energy levels of U 183+ 2 quasi-molecule "dive" into the Dirac negative continuum, leading to an instability of the physical vacuum against electron-positron pair production. By inspecting Fig. 2 and Table 2 one can see that the finite-nuclear size correction to the energy levels of U 183+ 2 varies significantly with the distance between the nuclei. For the ground, 1σ g , state, for example, this correction does not exceed 0.1% of the total energy at the "chemical" inter-nuclear distance R = 2/92 a.u. (see also right column of Table 1) but increases to almost 10% if R ≈ 40 fm. Even larger effect can be observed at this distance for the first excited state for which the absolute value of the binding energy 1σu reduces by almost factor of two if one accounts for the finite nuclear size effects (see Table 2). These corrections to the energy levels of heavy quasimolecule result, of course, in a significant dependence of the critical radius R crit , i.e., of the distance at which these levels "dive" into the negative continuum, on the nuclear charge distribution. For the ground state, for example, our finite-basis-set calculations suggest that the R crit decreases from 38.4 fm for the point-like to 34.7 fm for the finite- (Tupitsyn et al. 2010). Moreover, our data are in accordance with the value R crit = 35 fm obtained in Ref. (Soff, Reinhardt & Betz 1978) within the monopole approximation.

In contrast to the 1σ g ground state for which the finite-size effects just lead to a reduction of the critical radius R crit , they prevent completely the excited 1σ u orbital from the "diving" into the negative continuum before the collision of nuclei. This can be seen from Fig. 2 whose left border indicates the nuclear diameter of Uranium. At this distance, the binding energy 1σu =-0.947109 n.u. is still far from the negative continuum threshold.

Summary and outlook

In summary, we have elaborated a finite-basis-set approach to deal with the two-center Dirac problem. Being constructed in Cassini coordinates, the basis functions allow one to solve the generalized eigenvalue equation ( 21) for all the range of inter-nuclear distances and to find a (quasi) complete set of eigensolutions. In order to illustrate the application of the method and to verify its accuracy, detailed calculations have been performed for the binding energies of H + 2 , Th 179+ 2 and U 183+ 2 dimers. Special attention in these calculations has been paid to the finite nuclear corrections that were found to be of a paramount importance for an accurate description of molecular spectra especially at small distances between the nuclei. For the U 183+ 2 quasi-molecule, for example, these corrections not only reduce (the absolute value of) the ground-state energy by about 10 % but also prevent the first excited level 1σ u from reaching the Dirac negative continuum.

The finite-basis-set method, yielding a complete set of eigenvalues and eigenfunctions of the two-center Dirac Hamiltonian, provide a powerful tool for studying the dynamics of (slow) ion collisions including the electron-positron pair production. Moreover, it can be efficiently used for an accurate evaluation of radiative (QED) corrections to molecular energy levels, that requires the summation over the complete spectrum. Theoretical investigations towards these two lines are currently underway.

Figure 2 .

 2 Figure 2. The energy of 1σ states in U 183+ 2 as a function of internuclear distance.The solid and dashed lines correspond to finite-size and point-like nuclei, respectively. The dash-doted line represents the energy of the 1s state in spherically-symmetrical mono-pole approximation, centered in the middle point of inter-nuclear axis.

Table 1 .

 1 Binding energy of 1σ g quasi-molecular state in atomic units.

	System	H + 2 , R=2 a.u. T h 179+ 2	, R=2/90 a.u. U 183+ 2	, R=2/92 a.u.
	This work, point-like nucleus -1.1026409	-9504.752	-9965.375
	Others, point-like nucleus	-1.1026416 a	-9504.757 a	-9965.307 b
	This work, finite-size nucleus -1.1026401	-9498.448	-9957.796
	a Ref. (Kullie & Kolb 2001)
	b Ref. (Tupitsyn et al. 2010)
	example, accurate data were reported for the values of the 1σ g ground state energy
	of the H + 2 and Th 179+ 2

Table 2 .

 2 1σ g and 1σ u energies in natural units as functions of internuclear distance R. These estimates were found to be in a perfect agreement with the previous calculations by Tupitsyn and co-workers who reported values of 38.43 fm and 34.72 fm, correspondingly

	R[fm]	ε 1σg	ε 1σu	ε 1σg	ε 1σu	ε 1s
		Point nucleus Point nucleus Finite nucleus Finite nucleus Mono-pole
	26.46	Neg. cont.	-0.46415	Neg. cont.	-0.32584	Neg. cont.
	38.5	-0.99842	-0.14956	-0.92831	-0.08908	-0.90882
	52.92	-0.76208	0.04524	-0.72506	0.07514	-0.68872
	264.59	0.00632	0.52711	0.00851	0.52835	0.05200
	529.18	0.24414	0.60211	0.24499	0.60280	0.29291
	2645.89 0.63382	0.65132	0.63418	0.65174	0.73055
	5291.77 0.69193	0.69219	0.69231	0.69258	0.84359
	26458.86 0.73134	0.73134	0.73172	0.73173	0.96572
	52917.72 0.73623	0.73624	0.73662	0.73662	0.98531
	size nuclei.

Acknowledgments

Elucidating discussions with V. M. Shabaev, I. I. Tupitsyn, Yu. S. Kozhedub, C. Kozhuharov, and V. A. Yerokhin are gratefully acknowledged. The authors are grateful to I. I. Tupitsyn for providing the results of his calculations prior to publication. We are also grateful to V. A. Yerokhin for useful remarks. A.A. and A.S. acknowledge support from the Helmholtz Gemeinschaft and GSI (Nachwuchsgruppe VH-NG-421). The support by the ExtreMe Matter Institute EMMI in the framework of the Helmholtz Alliance HA216/EMMI is also gratefully acknowledged. Laboratoire Kastler Brossel is Unité Mixte de Recherche du CNRS, de l' ENS et de l' UPMC No. 8552.


Zhenyu Tian

Tao Yuan

Rene Fournet

Pierre-Alexandre Glaude

Baptiste Sirjean

Dr Frédérique Battin-Leclerc
email: frederique.battin-leclerc@ensic.inpl-nancy.fr

Kuiwen Zhang

Fei Qi

Prof Fei Qi

An experimental and kinetic investigation of premixed furan/oxygen/argon flames

ABSTRACT: The detailed chemical structures of three low-pressure (35 Torr) premixed laminar furan/oxygen/argon flames with equivalence ratios of 1.4, 1.8 and 2.2 have been investigated by using tunable synchrotron vacuum ultraviolet (VUV) photoionization and molecular-beam mass spectrometry. About 40 combustion species including hydrocarbons and oxygenated intermediates have been identified by measurements of photoionization efficiency spectra. Mole fraction profiles of the flame species including reactants, intermediates and products have been determined by scanning burner position with some selected photon energies near ionization thresholds. Flame temperatures have been measured by a Pt-6%Rh/Pt-30%Rh thermocouple. A new mechanism involving 206 species and 1368 reactions has been proposed whose predictions are in reasonable agreement with measured species profiles for the three investigated flames. Rate-of-production and sensitivity analyses have been performed to track the key reaction paths governing furan consumption for different equivalence ratios. Both experimental and modeling results indicate that few aromatics could be formed in these flames. Furthermore, the current model has been validated against previous pyrolysis results of the literature obtained behind shock waves and the agreement is reasonable as well. KEYWORDS: Premixed laminar flame; Kinetic modeling; Furan; Tunable synchrotron VUV photoionization; Molecular-beam mass spectrometry.

Introduction

The interest for studying the gas-phase reactions of biofuels, mainly alcohols, ethers and biodiesels, has increased in recent years because they offer the long-term promise of fuel-source regenerability and reduced climatic impact [1]. In engine tests, Avantium Company demonstrated that furan-based biofuels, or "furanics" were potential biofuels and had significant advantages over first-generation biofuels [2,[START_REF] Avantium | Avantium steps ahead with its biofuels program-engine test demonstrates potential of "furanics[END_REF]. For example, 2,5-dimethylfuran (DMF) has an energy density 40% greater than ethanol, which is comparable to that of gasoline. Furan and its derivatives have also been identified among the emissions produced from the biomass burning and the combustion of fossil fuels [START_REF] Andreae | [END_REF][5][6]. Investigations with pyrolysis of leaf and wood indicate that furan is one of the major degradation products of these compounds [7]. In addition, furan has also been proved to be a significant component of tobacco smoke [8] and is selected as a model oxygenated reburn fuel for NO reduction [9]. Therefore, it is necessary to understand the combustion chemistry of furan.

Previous studies relevant to the gas phase reactions of furan were mainly focused on experiments and modeling of its thermal decomposition [10][11][12][13][14][15][16][17][18]. In 1972, Cullis and Norris [10] studied the pyrolysis of furan within the temperature range of 1173 -1323 K at 1 atm. The main products observed were methane (CH 4 ), acetylene (C 2 H 2 ), ethylene (C 2 H 4 ), and benzene (C 6 H 6 ).

With a heated flow reactor and an on-line mass spectrometry (MS), Grela et al. [11] investigated the low-pressure (10 -3 Torr, 760 Torr = 1 atm) pyrolysis of furan over the temperature range of 1050 -1270 K and found the presence of C 3 H 4 (allene, aC 3 H 4 /propyne, pC 3 H 4 ) and CO. They also proposed a high-pressure Arrhenius expression k ∞ = 10 15.6 exp(-73.5 kcal•mol -1 /RT) s -1 as the rate constant for furan decomposition. Later, Bruinsma et al. [13] studied the decomposition of furan within the temperature range of 960 -1085 K at 1 Torr. CO and C 3 H 4 were identified and a very different Arrhenius equation, k = 10 12.9 exp(-65.7 kcal•mol -1 /RT) s -1 , was reported. Three groups [12,14,15] carried out the thermal decomposition of furan behind shock waves. Lifshitz et al. [12] studied the decomposition of furan over the temperature range 1050 -1460 K at gas densities of 3 10 -5 mol/cm 3 (2.6 -3.6 atm) and proposed that the two main pathways of furan decomposition were due to CO + pC 3 H 4 and C 2 H 2 + ketene (CH 2 CO). By gas chromatography/mass spectrometry (GC/MS) analysis, they identified CO, pC 3 H 4 , C 2 H 2 and indirectly, CH 2 CO, as the major products and CH 4 , aC 3 H 4 , C 2 H 4 , diacetylene (C 4 H 2 ), vinylacetylene (C 4 H 4 ), 1,3-butadiene (C 4 H 6 ) and C 6 H 6 as minor products [12]. Based on their measurements, they obtained the rate constant k = 10 15.43±0.45 exp[-(78.3±2.0) kcal•mol -1 /RT] s -1 for the overall decomposition of furan in the temperature range 1060 -1260 K [12]. Organ and Mackie [14] extended the temperature range to 1100 -1700 K (20 atm) and obtained the same conclusion as that of Lifshitz et al. [12] by detecting CO, pC 3 H 4 , C 2 H 2 and CH 2 CO directly. The global rate of furan decomposition was given as k = 10 15.3±0.3 exp[-(77.9±1.9) kcal•mol -1 /RT] s -1 and the initiation was postulated to be the homolytic C-O bond scission leading to a biradical [14]. With laser-schlieren (LS) and time-of-flight (TOF) MS detection techniques, Fulle et al. [15] confirmed the two dominant dissociation channels in the temperature range 1300 -3000 K (600 Torr) and k = 10 15.3 exp[-78.7 kcal•mol -1 /RT] s -1 was given. Later, Winkler et al. [16] investigated the continuous flow pyrolysis of furan at 1173 K (1 atm) and up to 22 products were identified by GC-MS. They suggested that aromatics were generated from C 2 -and C 3 -species and found no tendency for the oxygen containing fragments of furan to form larger molecules. Recently, Hore et al. [17] studied the pyrolysis of furan using a combination of infrared laser powered homogeneous pyrolysis, chemical and physical trapping of radicals and precursors specifically designed to generate selected radicals in the temperature range 950 -1000 K (~7 Torr) and confirmed the key role of 1,2-H shifts in the theoretical studies (see discussion below). More recently, Vasiliou et al. [18] performed furan decomposition with a high-temperature nozzle as a tubular reactor. They also demonstrated the importance of (CO + pC 3 H 4 ) and (C 2 H 2 + CH 2 CO) routes; however, furan could also decompose to produce propargyl radical at a temperature higher than 1550 ± 100 K [18].

Despite the numerous experimental investigations of furan thermal decomposition, there were only three theoretical calculations concerning these reactions. Liu et al. [19,20] explored the thermal decomposition of furan with density functional (B3LYP) techniques for geometries and QCISD(T) for energies. They reported the energetics of furan decomposition channels, but gave no information about rate constants. In 2000, Sendt et al. [21] obtained the thermochemistry and rate parameters of several key reactions related to furan at CASSCF, CASPT2 and G2-(MP2) levels. They proposed a detailed kinetic model and validated it against furan pyrolysis data measured by Organ and Mackie [14]; they concluded that the formation of cyclic intermediates and the formation of decomposition products (CO + pC 3 H 4 ) initiated by 1,2-H transfer were the dominant pathways for furan pyrolysis. However, compared to such studies on the thermal decomposition of this molecule, no literature exists on its combustion.

In the present work, low-pressure premixed furan/O 2 /Ar flames with three equivalence ratios ( ) of 1.4 (flame A), 1.8 (flame B) and 2.2 (flame C) have been investigated with the molecular-beam mass spectrometry (MBMS) and tunable VUV synchrotron photoionization techniques [22][23][24].

Many flame species have been detected and identified by comparing the measured ionization energies (IE) with literature ones; their mole fraction profiles have been determined as well.

Moreover, a detailed kinetic model has been developed based on the experimental measurements and reaction pathways were proposed with particular attention to the formation of the major intermediates.

Experimental Section

The experiments were performed at National Synchrotron Radiation Laboratory, Hefei, China.

The instrument and data analysis procedures have been reported previously [25,26]. Briefly, the apparatus is composed of a flame chamber which contains a movable flat-flame burner (McKenna Burner) with 6 cm diameter, a differentially pumped chamber with a molecular-beam sampling system, and a photoionization chamber installed with a home-made reflectron time-of-flight mass spectrometer (RTOF-MS). Combustion species are sampled by a quartz cone with a 40° included angle and a ~500 μm orifice at the tip. The resulting molecular beam passes into a differentially pumped ionization region through a nickel skimmer, and then crosses with the tunable VUV synchrotron light in the photoionization chamber. The photoions are collected and analyzed by the RTOF-MS with an approximate mass resolution of 1400 [22]. The integrated ion intensities for a specific mass are normalized by the photon flux, and then plotted as a function of the photon energy or distance from the burner surface, which then yields photoionization efficiency (PIE) or mole fraction profiles, respectively. Considering the cooling effect of molecular beam [27], the experimental error for determining IE is within 0.05 eV for stable species and slightly higher for some radicals because of the weak signal-to-noise ratio. The uncertainties of the mole fractions are ±5-10% for the major species, ±25% for intermediates with known photoionization cross sections and a factor of two for those with estimated cross sections.

Three flames from slightly rich (Ф = 1.4, C/O = 0.538) to nearly sooting conditions (Ф = 2.2, C/O = 0.786) were investigated. Flow rates of the gas reagents are controlled separately by MKS mass flow controllers, while the flow rates of liquid fuels into the vaporizer are controlled by a syringe pump (ISCO 1000D, USA). Furan was purchased from Sinopharm Chemical Reagent Limited Co., Shanghai, China with the purity of ≥ 99.5%. No further purification was carried out for this study. The operating conditions are listed in Table 1. For comparison, the pressure was kept at 35 Torr; the flow rate of Ar has a constant value of 1.000 standard liter per minute (SLM) and the cold-gas velocity is 33.49 cm/s for all flames. The temperature profiles have been measured by using a Pt-6%Rh/Pt-30%Rh thermocouple, 0.100 mm in diameter and 15 mm upstream from the sampling cone, with an Y 2 O 3 -BeO anti-catalytic coating [27]. Radiation heat losses have been considered and calibrated for the temperature profiles [START_REF] Fristrom | Flame Structure and Processes[END_REF], the detailed procedures of which were described in [START_REF] Tian | [END_REF].

Based on the reproducibility of the measurement, the accuracy of the temperature measurements was estimated to be ± 100 K, within which the modeling results have a 3% uncertainty for the major species and about 10% for intermediates [START_REF] Tian | [END_REF].

Kinetic model

To simulate the species profiles measured in the three furan flames, a new kinetic mechanism has been developed. Table 2 lists the reactions related to furan and derived species. This mechanism was coupled with our recent mechanism which has been developed for the oxidation of toluene [30].

The complete mechanism (206 species and 1368 reactions) is available as supplementary material and only the primary and secondary reaction pathways are described here. Most hermochemical data have been estimated by the software THERGAS based on group additive method [31]. Some intermediate species derived from furan and those involving a missing group in the THERGAS software were estimated by ab initio calculation. The ab initio calculations were performed at the CBS-QB3 level [32] implemented in Gaussian 03 [START_REF] Frisch | Gaussian03, revision B05[END_REF]. Enthalpies of formation at 298 K have been obtained from isodesmic reactions, while thermochemical properties as a function of temperature have been calculated using the CHEMRATE software [START_REF] Mokrushin | Chemrate v.1.5.2[END_REF] according to statistical mechanical principles. Hindered rotors have been taken into account and relaxed scans have been performed at the B3LYP/6-31G+(2d,p) level to determine the barrier of rotation. Table 3 presents the names, the formulae and the enthalpies of formation (Δ f H°) of furan and its derived species which were not presented in our previous work [30].

Primary mechanism

The primary mechanism contains the reactions of furan and of the radicals directly derived from it. The unimolecular initiations of furan have been considered and Fig. 1 presents the simplified potential energy surface. According to Sendt et al.'s study [21], the most important unimolecular pathways for furan decomposition are the reactions involving 1,2-H migrations, which then lead to the formation of two different cyclic carbenes. These species can decompose to give the major products as CO + pC 3 H 4 and C 2 H 2 + CH 2 CO. For the first route, Sendt et al. [21] showed that two rate-determining steps must be taken into account. The first one corresponds to the initial 1,2-H transfer from furan to form formyl allene (FA, , CH 2 CCHCHO, 1,2-butadienal),

(reaction 1 (R1) in Table 2), while the second step is related to the decomposition of FA by a concerted reaction (TS2) in which a 1,4-H transfer and a simultaneous C-C bond fission occurs to form CO and pC 3 H 4 (R52). It must be noted that the decomposition of FA to propargyl and CHO radical (R51) has been also considered in our mechanism in accordance with the work of Sendt et al. [21].

The detailed mechanism for the second unimolecular reaction of furan involves three elementary steps. The first one corresponds to the formation of a cyclic carbene (1,2-H shift) as discussed previously. The second step is related to the ring opening of the cyclic carbene to produce a diradical while the third step involved two successive bond cleavages of this latter species to form C 2 H 2 and CH 2 CO. It is possible, as mentioned by Sendt et al. [21], to calculate the overall rate constant by considering the last transition state TS3 which corresponds to the highest energy barrier involved in the overall reaction (Fig. 1).

All the calculations have been performed with the CBS-QB3 method [32] and the activation energies found are in a good agreement with those obtained by Sendt et al. [21] at the G2B2 and CASS-SCF level. Frequency analysis makes it possible to point out one imaginary frequency for each transition state (TS). High-pressure limit unimolecular rate constants were calculated using canonical transition state theory. The kinetic parameters were obtained with a modified Arrhenius form (Table 2).

An interesting feature, resulting from previous calculations, is the relatively low barriers involved in the unimolecular initiations. Consequently, these reactions become very important and their rate constants are sensitive parameters. Thus, we performed RRKM calculations in order to take into account pressure effect under flame conditions. Time dependent master equation has been solved with the CHEMRATE program [START_REF] Mokrushin | Chemrate v.1.5.2[END_REF] for the two pathways involved in Fig. 1. Calculations have been performed at 4.7 kPa and T = 700-2000 K. A value of 260 cm -1 was used in the master equation analysis, with Ar as third body [START_REF] Joshi | [END_REF]. Lennard-Jones parameters for furan and other isomers (FA) were estimated as σ = 5.18 Å and = 357 K in accordance with Joshi and Wang [START_REF] Joshi | [END_REF]. Rate constants, k(E) have been estimated at 25 cm -1 (0.07 kcal mol -1 ) increments. The lowpressure rates (4.7 kPa) were used in the simulation of flame results, while the high pressure ones in pyrolysis conditions.

The kinetics of the unimolecular decompositions to give H atoms and furyl radicals (R3 and R4) have been deduced from those of the reverse reactions, k = 1.0×10 14 s -1 according to Allara and Shaw [36]. For the bimolecular initiations with oxygen molecules (R5 and R6), A-factors have been obtained using the correlation proposed by Ingham et al. [37] and the activation energies have been taken equal to the enthalpies of reaction. The rate constants of the additions of H and OH to furan (R7-R9) have been deduced from value proposed for the similar reactions in case of 1,3-butadiene [38] and our previous work on the oxidation of alkenes [39], respectively. The rate for the reaction of furan with CH 3 forming 2-methylfuran ( ) and H (R10) has been taken from the reverse reaction proposed by Lifshitz et al. [40]. The kinetics for the abstractions of vinylic H atoms (R11-R18) were deduced from the correlations proposed by Heyberger et al. [39].

Two furyl radicals, furyl-2 ( ) and furyl-3 ( ) are formed by H-abstractions and two The decomposition reactions of furyl radicals have been also investigated at the CBS-QB3 level of theory [32]. Even if the bond dissociation energy of C-H is very high (about 119 kcal mol -1 according to our calculations), the decomposition of these radicals could play a significant kinetic role under flame conditions where the high temperature and concentration of radicals favor Habstractions from furan. The rate constants for these reactions were estimated as proposed by Allara et al. [36].

Secondary mechanism

The secondary mechanism includes the reactions of the primary products which are not considered in the mechanism of the oxidation of toluene [30], namely FA, methylfuran, but-3-ynal, 3-butenal, 2-methyl-but-3-ynal, 2-methyl-3-butenal, furan23H, furan25H and tetrahydrofuran ( , THF). As discussed previously, FA can decompose to pC 3 H 4 + CO (R52), the rate constant of which has been determined using quantum calculations with CBS-QB3 method [32]. FA is also a source of C 3 H 3 and CHO radicals (R51), the kinetics of which is from Sendt et al. [21] with the A-factor divided by 2.2 to take into account the pressure effect (for R1, the k ∞ /k 4.7kPa = 2.2 at 1000 K). In the case of methylfuran, the decomposition by breaking a C-H or C-O bond (R53-R60), bimolecular initiations with oxygen molecules (R61) and H-abstractions (R62-R71) are considered. The rate constants of these reactions were taken from Lifshitz et al. [40], except those of the bimolecular initiation with O 2 [46], and those reactions with O or C 2 H 3 [START_REF] Pitz | Proceedings of the 2nd Joint Meeting of the U.S. Sections of the Combustion Institute[END_REF], with OH [START_REF] Knispel | [END_REF] and with HO 2 [START_REF] Mehl | Fall Technical Meeting of the Western States Section of the Combustion Institute[END_REF] which were obtained by analogy with toluene. For but-3-ynal, 3-butenal, 2-methylbut-3-ynal and 2-methyl3-butenal, the H-abstractions by H, OH and CH 3 radicals followed by decompositions (R79-R90) were taken into account, the kinetics of which were deduced from the similar reactions of acetaldehyde [45,[START_REF] Warnatz | Rate coefficients in the C/H/O system[END_REF][START_REF] Cavanagh | [END_REF].

For dihydrofuran and THF, the rates for their bimolecular reactions with oxygen (R91-R93, R118 and R119) were estimated from the correlation proposed by Ingham et al. [37] and those of metatheses reactions (R94-R108, R120-R129) were deduced from the correlations proposed by Heyberger et al. [39]. H-additions to dihydrofuran would give THF-2-yl (R109) and THF-3-yl (R110 and R111), respectively, and OH-additions would produce the same products HCHO + C 2 H 4 + CHO (R112 and R113), the kinetics of which were from the correlations proposed by Heyberger et al. [39].

By β-scission, THF-2-yl and THF-3-yl could decompose to C 2 H 4 + CH 2 CHO (R114) and HCHO + C 3 H 5 (R115),respectively, with the rates proposed by Sirjean et al. [52]. Terminations of the two tetrahydrofuryl radicals would form THF (R116 and R117), the kinetic data of which were estimated from the work of Allara and Shaw [36].

Results and Discussion

Experimental results

In this study, about 40 combustion intermediates, including hydrocarbons and oxygenated species, were detected and identified in the flames characterized by Ф = 1.4, 1.8 and 2.2. These species are listed in Table 4 along with the literature and measured IEs, peak concentrations and their locations in the flames. The detailed experimental data of mole fractions are available in supplemental material. To illustrate the effect of the equivalence ratio on the flame structure, Figures H 2 O, CO and CO 2 have similar tendencies for the position with equilibrium value as that of Ar. In the post-flame zone, the mole fractions of Ar and CO 2 steadily decrease with the equivalence ratio increasing, while those of H 2 and CO show reverse trends. As evidenced from Fig. 5, the post flame H 2 O levels were not sensitive to equivalence ratio. However, equivalence ratio has a significant impact on the X CO /X CO2 ratio. In the post-flame zone, the X CO /X CO2 ratio is about 1.2 for flame A.

This ratio increases to 2.8 in flame B and reaches 8.3 in flame C. The modeling of the chemical kinetics was performed using the CHEMKIN software [53], with the experimental temperature profiles introduced as input parameters for simulation. To account for the perturbations induced by the quartz probe and the thermocouple [54] and fit the experimental results, the temperature profile was shifted 1.5, 1.7 and 2.2 mm away from the burner surface for the flames with Ф = 1.4, 1.8 and 2.2, respectively. The perturbation effect is known in molecular-beam mass spectroscopic investigations, and a shift equaling 4-5 times the orifice diameter has been reported for premixed low-pressure flames (4.0-5.0 kPa) [START_REF] Tian | [END_REF]54,55].

The comparisons between the experimental measurements and predictions for the three flames are presented in Figs. 5678, in which lines stand for simulated results. There is generally a reasonable agreement between measurements and predictions. For the major species shown in Fig. 5

Rate-of-production and sensitivity analysis for furan consumption in flames

To illustrate the major reaction channels of furan, flux analysis has been performed for two equivalence ratios 1.4 (flame A, at 3.7 mm from the burner, T = 1470 K, 88.9% furan conversion) and 2.2 (flame C, at 7.9 mm from the burner, T = 1609 K, 85.0% furan conversion). As presented in

Aromatic formation

Among the products of furan pyrolysis analyzed by Winkler et al. [16], a lot of aromatic species such as naphthalene, phenanthrene, pyrene and acenaphthylene were identified. However, even under the richest condition (flame C), only benzene, toluene, phenylacetylene (C 8 H 6 ) and styrene (C 8 H 8 )were identified in the present work. This could be due to the relatively low temperature of the studied flames, as well as to the existence of O atom in furan, hindering the formation of much larger aromatics. In the following section, we will focus on the formation pathways of these species, based on the flux analysis taken at the same positions as mentioned in 4.

Validation of the present mechanism using pyrolysis results

To extend the validity of the present mechanism, previous pyrolysis results from Organ and Mackie (20.0 atm with initial furan concentration of 2.0%, 1.0% and 0.2% dilute in Ar) [14] and Fulle et al. (0.26 atm with furan concentration of 2.0% dilute in Ne) [15] were modeled. Considering that the pressures used in pyrolysis experiments are higher than in the flames, high pressure rate constants for reactions R1, R2, R51 and R52 that were calculated using quantum calculations with CBS-QB3 method [32] have been used in simulations.

Figure 11 compares the predicted furan and product profiles with the experimental results (shock tube, residence time 300 μs) reported by Organ and Mackie [14], in which symbols are measurements and solid, dash dot and dashed lines are respective predictions by the current model for 2.0%, 1.0% and 0.2% initial furan concentration. Also, the predictions for 2.0% initial furan concentration by the mechanism of Sendt et al. [21] ). The model of Sendt et al. [21] shows worse prediction of C 3 H 4 than ours in the temperature range of 1450 -1525 K, but better prediction when the temperature is higher than 1525 K.

Figure 12 shows the flux analysis result for 2.0% initial furan concentration at 1527 K, corresponding to 64.8% of furan conversion, at which both propyne and benzene have peak mole fractions experimentally. As noted in the result of Sendt et al. [21], their model supported a mechanism involving two parallel furan decomposition routes to the major products CO + pC Figure 13 presents the comparison of the predicted furan and product profiles with the experimental results (2% furan/98% Ne, 1533 K, 198 Torr) reported by Fulle et al. [15], in which symbols are measured data and solid and dash lines are respective predictions by the current model and the mechanism of Fulle et al. [15]. Although CO and propyne are underpredicted slightly at the end of the observation period, a globally correct prediction can be observed for the current model. It should be noted that no other products were detected within the period of ~800 μs and only the two unimolecular dissociation reactions were involved in the model of Fulle et al. [15]. They also determined the ratio [C 2 H 2 ]/[CO] and gave an expression [C 2 H 2 ]/[CO] = (5.5 10 -9 ) T 2.5 by fitting their own results and that of Lishitz et al. [12]. In this work, the predicted ratios of [C 2 H 2 ]/[CO] are 0.45 at 400 μs and 0.44 at 800 μs, which are very close to respective 0.42 and 0.40 that are calculated using the fitted expression of Fulle et al. [15].

To analyze the difference of the rate constants calculated in this work for furan pyrolysis, it is necessary to compare them with previous measurements and calculations, as displayed in Fig. 14.

The unimolecular reactions for furan pyrolysis are commonly accepted as dissociation to CO + pC 3 H 4 and C 2 H 2 + CH 2 CO [12,14,15,18,21], while CO + pC 3 H 4 is dominant at temperatures lower than 1270 K [11,13]. At temperatures higher than 1550 K, furan can generate C 3 H 3 [18]. As seen in Fig. 14, only the total reaction, including routes of CO + pC 3 H 4 and C 2 H 2 + CH 2 CO are considered.

In the current model and mechanism of Sendt et al. [21], the reaction furan = FA is used since it is the rate-determining step to form CO + pC 3 H 4 . Generally, the rate constant of this work is in good agreement with previous ones predicted by Lishitz et al. [12], Bruinsma et al. [13], Organ and

Mackie [14] and Sendt et al. [21]. The low-temperature rate constants given by Grela et al. [11] are higher than the current model and measurements of Organ and Mackie [14], which was also observed by Sendt et al. [21]. The model of Fulle et al. [15] showed good agreement with the present mechanism at temperatures higher than 1300 K, but displays a slightly lower rate constant than the current ones at lower temperatures. and reasonable agreement has been obtained between experimental measurements and simulations.

Conclusions

The current mechanism has also been successfully validated against previous data obtained for furan pyrolysis. However, further studies still remain to improve this first attempt of mechanism for the combustion of furan.

In flames, the width of the reaction zone is increased and the peak positions of most intermediates are shifted downstream with the equivalence ratio increasing; furan is consumed mainly by metatheses leading to furyl radicals, OH addition to the double bond of furan giving 2-propenal and CH=CH-CHO radicals by β-scission, isomerization forming formyl allene via 1,2hydrogen migration; benzene is formed mainly by isomerization of fulvene which comes from the combination of the propargyl radicals. A sensitivity analysis indicates that the reaction furan = CH 2 CCHCHO has a strong promoting effect, while the decomposition of CH Note: the rate constants are given in the form k = A T n exp(-E a /RT) where A has units of cm, mol and s, T has units of K, and E a has units of kcal/mol. a Rate constant calculated in this work. b Rate constant taken equal to that of the recombination of H atoms with alkyl radicals as proposed by Allara and Shaw [36]. c Rate constant calculated as proposed by Ingham et al. [37], but E a corrected with enthalpies of the bimolecular reaction. d Rate constant estimated by analogy with the values proposed by Gueniche et al. [38] for the similar reaction of 1,3butadiene. e Rate constant estimated by using the correlations proposed by Heyberger et al. [39] in the case of alkenes. f Rate constant taken equal to the values proposed by Lifshitz et al. [40] for methylfuran. g Rate constant taken equal to that of the isomerization of iC 4 H 5 and nC 4 H 5 radicals as proposed by Leung and Lindstedt [41]. h Rate constant estimated by analogy with values proposed by Mebel et al. [42] for reactions of O 2 + C 2 H 3 . i Rate constant estimated by analogy with values proposed by Dagaut et al. [43] for reactions of O + C 2 H 3 . j Rate constant taken equal to the values proposed by Gueniche et al. [38] for the similar reaction of lC 4 H 7 -1 and lC 4 H 7 Y. k Rate constant estimated by using the correlations proposed by Touchard et al. [44] in the case of low temperature oxidation of 1-pentene. l Rate constant estimated by analogy with values proposed by Baulch et al. [45] for the similar reaction of O + C 2 H 5 . m Rate constant calculated using the modified collision theory at 1200 K using software KINGAS [56]. n Rate constant estimated by analogy with values proposed by Baulch et al. [45] for reaction of C 2 H 2 + CH 3 = sC 3 H 5 . o Rate constant taken equal to the values proposed by Sendt et al. [21] for formyl allene, with the A-factor divided by 2.2 for flame conditions. p Rate constant estimated by analogy with values proposed by Oehlschlaeger et al. [46] for reaction of toluene + O 2 . q Rate constant estimated by analogy with values proposed by Pitz et al. [START_REF] Pitz | Proceedings of the 2nd Joint Meeting of the U.S. Sections of the Combustion Institute[END_REF] for reaction of toluene + O. r Rate constant estimated by analogy with values proposed by Knispel et al. [START_REF] Knispel | [END_REF] for reaction of toluene + OH. s Rate constant estimated by analogy with values proposed by Mehl et al. [START_REF] Mehl | Fall Technical Meeting of the Western States Section of the Combustion Institute[END_REF] for reaction of toluene + HO 2 . t Rate constant estimated by using the correlations proposed by Brezinsky et al. [57] in the case of benzyl radical. u Rate constant estimated by analogy with values proposed by Bounaceur et al. [58] for reaction of C 6 H 5 CH 2 O. v Rate constant estimated by analogy with values proposed by Warnatz [50] for the similar reaction of CH 3 CHO + H. w Rate constant estimated by analogy with values proposed by Cavanagh et al. [START_REF] Cavanagh | [END_REF] for the similar reaction of CH 3 CHO + OH.

x Rate constant estimated by analogy with values proposed by Baulch et al. [45] for the similar reaction of CH 3 CHO + CH 3 . y Rate constant estimated by analogy with values proposed by Sirjean et al. [52] for the ring opening reaction of C 5 H 9 #. (*) This energy has been taken from Sendt et al. [21].

[ [ [ [ [ [ [ (flame C, at 7.9 mm from the burner, T = 1609 K, 85.0% furan conversion).

[ [ [

Fig. 12 Rate-of-production analysis for the consumption of furan at a temperature of 1527 K and a 64.8% conversion of furan in the furan pyrolysis.

[

C 4 H 5 O

 45 radicals, dihydrofuryl-2 ( , C 4 H 5 O-2) and dihydrofuryl-3 ( , C 4 H 5 O-3) are produced by H-additions. Isomerization reactions between the two furyl (R19) and the two C 4 H 5 O radicals (R28) are possible. The rate constants of the isomerizations are deduced from the similar reaction of n-C 4 H 5 and i-C 4 H 5 radicals for furyl radicals [41], and from that of CH 2 =CH-CH 2 -CH 2 and CH 2 =CH-CH-CH 3 radicals for C 4 H 5 O radicals [38].

Figure 2

 2 Figure 2 shows the major routes of decomposition of furyl radicals. All the activation energies

4 to 8

 8 present the temperature and concentration profiles. In Fig.4, the temperature profiles are shown. Due to size of the thermocouple and its perturbation on the flame structure, the temperatures close to the burner (1.0 mm above) were extrapolated to the burner surface. As shown in Fig.4, the flame temperatures steadily decrease and their peak locations shift downstream with the equivalence ratio increasing. In flame A, the maximum temperature is 1888 K located at 4.4 mm from the burner, while this value is decreased to 1772 K with a location of 6.2 mm in flame B and reaches 1699 K at 7.1 mm in flame C.Mole fraction profiles of the major species including furan, O 2 , Ar, H 2 , H 2 O, CO and CO 2 are displayed as symbols in Fig.5. The increase of the equivalence ratio in the furan flames causes a noticeable width increase of the reaction zone. In flame A, mole fraction of Ar remains constant after the position of 5.0 mm from the burner surface. This distance where Ar keeps constant shifts downstream to 7.0 mm in flame B and 9.0 mm at flame C. Among the three flames, both furan and O 2 consumed in further distances with the equivalence ratio increasing. The other major species, H 2 ,

Figures 6 to 8

 8 Figures 6 to 8 display the mole fraction profiles of combustion intermediates in the furan

  photoionization cross section for some intermediates can be another issue of the disagreement. As

Fig. 9 , 2 .

 92 Fig. 9, the major routes are quite different between the two flames. Metatheses by H, O and OH to

  photoionization and molecular-beam mass spectrometry. About 40 combustion intermediates and

Color figure in electronic versions only [Figure 1 ]Fig. 1

 11 Fig. 1 Simplified potential energy surface (PES) for the unimolecular initiation of furan. The energies (in kcal mol -1 ) are relative to furan and calculated at the CBS-QB3 level of theory at 0 K.

Fig. 2

 2 Fig. 2 Unimolecular decompositions of furyl radicals considered in this study. Activation energies (in

Fig. 3 [Figure 4 ]Fig. 4

 344 Fig. 3 Unimolecular decompositions of dihydrofuryl radicals considered in this study. Activation

Fig. 5

 5 Fig. 5 Experimental (symbols) and predicted (lines) concentration profiles of Ar, O 2 , furan, H 2 , H 2 O,

Fig. 6

 6 Fig. 6 Experimental (symbols) and predicted (lines) concentration profiles of (a) CH 3 (methyl

Fig. 7 6 ( 1 , 3 -

 7613 Fig. 7 Experimental (symbols) and predicted (lines) concentration profiles of (a) C 4 H 2 (diacetylene),

Fig. 8

 8 Fig. 8 Experimental (symbols) and predicted (lines) concentration profiles of (a) HCHO

Fig. 9

 9 Fig. 9 Rate-of-production analysis for the consumption of furan at two equivalence ratios, (a) Φ =

Fig. 10

 10 Fig. 10 Sensitivity analysis for the conversion of furan at two equivalence ratios, (a) Φ = 1.4 (flame

Fig. 11

 11 Fig. 11 Temperature dependence of designated species in the pyrolysis of furan (dilute in Ar, 20

Fig. 13

 13 Fig. 13 Species mole fractions versus residence time during the pyrolysis of furan (2% furan/98%

  that C 2 H 2 , the major precursor of C 4 H 2 , is more abundantly formed as Фincreases.Figure8shows the mole fraction profiles of the oxygenated intermediates. As Ф increases from 1.4 to 2.2, the peak concentration of 2-methylfuran is approximately doubled, while those of CH 2 CO

	The mole fraction profiles for C 4 -C 7 species presented in Fig. 7 also exhibit similar shift trends for the peak positions. With Ф increasing from 1.4 to 2.2, the maximum mole fractions of 1,2,3-butatriene/vinylacetylene (C 4 H 4 ), 1,3-butadiene (C 4 H 6 ), 1,3-cyclopentadiene (C 5 H 6 ), 1,3,5-(C 6 H 5 OH) by a factor of 10. Comparatively, the level of formaldehyde (HCHO) and ethanol/acetaldehyde (C 2 H 4 O) are not sensitive to equivalence ratio in the flame regime. 4.2 Kinetic modeling results hexatriyne (C 6 H and methylketene/2-propenol (C 3 H 4 O) are multiplied by one and a half, and that of phenol 4.2.1 Simulation of the flame structure

2, the peak concentrations of methane (CH 4 ), acetylene (C 2 H 2 ), propargyl radical (C 3 H 3 ) and allene (aC 3 H 4 ) are approximately doubled, while that of ethylene (C 2 H 4 ) is multiplied by one and a half, and that of propyne (pC 3 H 4 ) by a factor of 5. For allyl radical (C 3 H 5 ), while, the difference in the maximum mole fraction observed in flame A and B is acceptable, the agreement strongly deteriorates in flame C. 2 ) and benzene (C 6 H 6 ) are approximately doubled, while that of diacetylene (C 4 H 2 ) is multiplied by a factor of 5, and those of cyclopentadienyl radical (C 5 H 5 ) and toluene (C 7 H 8 ) are tripled. Thus, the equivalence ratio has the largest effect on the formation of C 4 H 2 . This can result from the fact

  , the current model captures well the experimental values within the error limits. With Ф increasing from 1.4 to 2.2, the model also gives evidence of the increase of the reaction width, which is consistent with the measurements. The present mechanism overpredicts the formation of CO in flames A and B near the burner surface. The flux analysis indicates that CO is mainly formed by the reaction of O 2 with formyl radical (CHO) which is produced by addition of OH radical to the C=C double bond of furan. Figures 6 and 7 summarize the results for some selected C 1 -C 7 intermediates. Overall, the current model reproduces fairly well the mole fraction profiles of C 2 H 2 , C 3 H 3 , pC 3 H 4 , C 3 H 5 , C 4 H 2 , C 4 H 4 , C 4 H 6 , C 5 H 6 , C 6 H 2 , C 6 H 6 and C 7 H 8 . The predicted peak concentration of CH 3 agrees well with the measurement in flame A. The present mechanism tends to overpredict the contributions fromCH 3 in flames B and C. However, the agreement is still within the expected uncertainties (a factor of two for free radicals). In the post-flame zone, CH 3 has a too high predicted concentration, which is produced mainly from the reactions of H atom with pC 3 H 4 , methanol (CH 3 OH) and CH 2 CO based on the rate-of-production (ROP) analysis. Compared to the reasonable capture of the CH 4 mole fraction shape for flames A and B, the current model exhibits a longer tail for flame C, which could be due to the reactions of abundant CH 3 with HCHO and CHO radical. With Ф increasing, the maximum mole fraction of C 2 H 4 becomes larger and the peak position shifts further. While the simulated mole fraction profiles of C 2 H 4 have the same tendencies as the experimental measurements, the peak values are underpredicted in flames A and B. The ROP analysis at the peak positions indicates that C 2 H 4 is mainly formed by the reactions of CH 2 CO+CH 2 and C 2 H 3 CHO+H in all three flames. Note that the formation of ethane (C 2 H 6 ) is negligible in all three flames. As revealed by the ROP analysis, C 2 H 6 is mainly formed through the combination of CH 3 radicals and most of C 2 H 6 converts to C 2 H 5 radical and finally to C 2 H 4 . For aC 3 H 4 , the current model gives fairly reasonable predictions in flames A and B, but a poor one in flame C. According to the ROP analysis result, most of aC 3 H 4 comes from C 3 H 5 . As seen Fig. 6 (h), the maximum mole fraction of C 3 H 5 is somewhat overpredicted in flames B and C, which could then lead to the overprediction of aC 3 H 4 .It should be noted that in flame C, the maximum mole fraction of C 4 H 2 is about 2.5 times overpredicted. By checking the ROP result, the reaction C 2 H 2 + C 2 H accounts for more than 80% of C According to the modeling result, the ratios of their maximum mole fraction X vinylacetylene /X 1,2,3-butatriene are 8.8, 8.2 and 4.9 in flames A, B and C, respectively, which demonstrates that vinylacetylene is the dominant species. The ROP result shows that about 70% vinylacetylene comes from the reaction of C 3 H 3 with triplet CH 2 . As seen in Fig.7(d), the ratio of the experimental and predicted peak mole fraction for C 5 H 5 is within 0.5 ~ 2.0.Considering the experimental and modeling uncertainties, mole fraction profiles of both species are determined accurately enough by the current model. In flame A, the metatheses of 1,3cyclopentadiene with H and decomposition of phenoxy radical accounts for 45% and 33% C 5 H 5 formation, respectively. However, decomposition of phenoxy radical is the major source (81%) of C 5 H 5 in flame C. The predicted peak values of the two C 6 H 6 isomers, benzene and fulvene (vC 6 H 6 ) reveals that benzene is dominant for mass 78, which is in agreement with experimental observations. An increase after the peak value is observed for the predicted C 7 H 8 profile in flame C. This phenomenon is caused by the combination of C 6

The too large formation of CO results from the too fast rate of the reaction OH + furan. The rate constant of this reaction needs to be more accurately determined by further experimental or theoretical studies. On the other hand, measurement uncertainties cannot be ruled out especially since the reaction zone is narrower in flames A and B than that in flame C. 4 H 2 formation. As C 2 H 2 is satisfactorily predicted, this disagreement could result from an experimental error. For mass 52 represented in Fig.

7

(b), the PIE spectra indicate that it involves the contribution of both 1,2,3-butatriene and vinylacetylene.

  2.2.For all the investigated flames, the reaction fulvene + H = C 6 H 6 + H is an important formation channel of C 6 H 6 (36.1% for flame A and 24.2% for flame C), where fulvene is mostly formed by the combination of C 3 H 3 radicals and by reaction between C 3 H 3 and C 3 H 5 radicals. Besides this reaction, the combination of C 3 H 3 radicals is another significant source of C 6 H

	combination reaction C 6 H 5 + CH 3 (56.5% for flame A and 91.8% for flame C), in which C 6 H 5
	radicals are dominantly generated from the combination of C 3 H 3 radicals. In contrast, the reaction
	toluene = benzyl + H is also important for toluene production in flame A (39.6%), but it is a
	consumption route in flame C. In flame C, C 8 H 6 is formed mainly by the reactions of C 2 H 2 and C 4 H 2
	with C 6 H 5 radicals, and C 8 H 8 is produced mainly by the decomposition of 2-phenylethyl radicals and
	by the combination reaction C 6 H 5 + C 2 H 3 .

6 (23.4% for flame A and 71.9% for flame C). The ipso-addition of H atom to C 6 H 5 OH (27.8%) and the reaction C 6 H 5 + H (+M) = C 6 H 6 (+M) (6.2%) also contribute to the C 6 H 6 formation in flame A. However, these two reactions have no contribution to C 6 H 6 production in flame C. The flux analysis indicates that the formation pathway of toluene is not sensitive to the equivalence ratio. It is dominated by the

  are given for comparison. In general, the present model gives acceptable agreement with experimental furan profile. Particularly, in the temperature range 1450 -1550 K, the current predictions are slightly closer to the measurements than those obtained with the model of Sendt et al.[21]. Good agreement is also obtained for C 2 H 2 , C 3 H 4 and C 6 H 6 . When the temperature is higher than 1500 K, CO is somewhat underpredicted by the present mechanism under all inlet conditions. This underprediction is also observed by the model of Sendt et al.[21]. While the model of Sendt et al.[21] gave underprediction of the minor product CH 4 in the temperature range > 1500 K, the current model also tends to underpredict CH 4 under all inlet conditions. As indicated from the ROP analysis, CH 4 is mainly formed by the reaction of CH 3 radicals with C 3 H 4 . Since C 3 H 4 is reasonably well predicted, the underprediction of CH 4 can be due to the consumption of CH 3 radicals through the reaction furan + CH 3 = methylfuran + H. C 2 H 4 is fairly predicted for 0.2% inlet condition with the present mechanism. However, it is underpredicted for 2% and 1% inlet conditions, which also can result from a too large consumption of CH 3 radicals, which are a key source of C 2 H 4 (CH 3 ↔ C 2 H 4 ↔ C 2 H 6

  3 H 4 and C 2 H 2 + CH 2 CO, which is exactly the same as observed in this work. As seen in Fig.12, furan is mainly converted to FA, which is also the same as that observed in flame C. Under pyrolysis conditions, the main products from FA are pC 3 H 4 + CO, while most FA gives C 3 H 3 + HCO in flame C. Besides to FA, decomposition to C 2 H 2 + CH 2 CO is another important consumption pathway of furan pyrolysis (25.2%). However, this channel is a minor pathway in flames even under the richest condition. In addition, 7.8% of furan disappears by reactions with H and CH 3 radicals, forming methylfuran, furyl-2, furyl-3, C 4 H 5 O-2 and C 4 H 5 O-3 radicals. Furyl-2 can decompose to C 2 H 2 + HCCO and CH≡C-CH=C=O + H via CH=CH-CH=C=O radical. By β-scission to CH≡C-CH=CH-O, furyl-3 can give C 3 H 3 + CO and CH≡C-CH=C=O + H. In contrast, the pathway to CH≡C-CH=C=O + H from CH=CH-CH=C=O is not observed from the flame studies. Moreover, C 4 H 5 O-2 and C 4 H 5 O-3 finally produce C 2 H 2 + CH 2 CHO and C 3 H 5 + CO, respectively, while the formations of these two radicals are much more involved in the consumption of furan in the flames.

Table 1

 1 2 CCHCHO to CO and pC 3 H 4 in flame with Ф = 1.4, and CH 2 CHCHCHO to CO and C 3 H 5 radical in flame with Ф = 2.2 have strong inhibiting effects on furan decomposition. Under pyrolysis conditions, the major decomposition products of furan are CO + pC 3 H 4 and C 2 H 2 + CH 2 CO. Moreover, the present rate constants of furan dissociation at high pressures based on ab initio calculations are in good agreement with previous experimental measurements and calculations. Flame characteristics Notes: a Ф refers to the equivalence ratio; b D M is the mass flow rate of the gas mixture.

	Flame Ф a	Flowrate (SLM) Furan O 2 Ar	C/O	b (g/s/cm 2 ) D M
	A	1.4 0.328 1.053 1.000 0.538 2.526E-3
	B	1.8 0.395 0.986 1.000 0.667 2.590E-3
	C	2.2 0.453 0.928 1.000 0.786 2.646E-3

Table 2

 2 Reactions of furan and derived species

	Reaction		A	n	E a	Footnote	No.
			PRIMARY MECHANISM
			Reactions of furan
	Unimolecular initiations by breaking of a C-O or C-H bond
	Furan = FA	(P = 4.7 kPa)	6.0 × 10 58	-13.388 94.0		a	(1)
		(P → ∞)	2.3 × 10 12	0.416	70.891	a	(1')
	Furan = C 2 H 2 + CH 2 CO	(P = 4.7 kPa) (P → ∞)	3.8 × 10 73 1.8 × 10 14	-17.412 113.3 0.534 86.591	a a	(2) (2')
	Furyl-2 + H = furan		1.0 × 10 14	0.0	0.0		b	(3)
	Furyl-3 + H = furan		1.0 × 10 14	0.0	0.0		b	(4)
	Bimolecular initiations
	Furan + O 2 = furyl-2 + HO 2 Furan + O 2 = furyl-3 + HO 2		2.0 × 10 13 2.0 × 10 13	0.0 0.0	70.68 70.88		c c	(5) (6)
	Additions and decompositions
	Furan + H = C 4 H 5 O-2 Furan + H = C 4 H 5 O-3 Furan + OH = C 2 H 3 CHO + CHO Methylfuran + H = furan + CH 3	2.6 × 10 13 2.6 × 10 13 2.7 × 10 12 5.0 × 10 14	0.0 0.0 0.0 0.0	3.2 1.56 -1.04 8.0		d d e f	(7) (8) (9) (10)
	Metatheses
	Furan + H = furyl-2 + H 2 Furan + H = furyl-3 + H 2 Furan + O = furyl-2 + OH		8.2 × 10 5 8.2 × 10 5 1.2 × 10 11	2.5 2.5 0.7	12.28 12.28 8.96		e e e	(11) (12) (13)
	Furan + O = furyl-3 + OH		1.2 × 10 11	0.7	8.96		e	(14)
	Furan + OH = furyl-2 + H 2 O Furan + OH = furyl-3 + H 2 O Furan + CH 3 = furyl-2 + CH 4 Furan + CH 3 = furyl-3 + CH 4	2.2 × 10 6 2.2 × 10 6 1.4 × 10 0 1.4 × 10 0	2.0 2.0 3.5 3.5	2.78 2.78 12.9 12.9		e e e e	(15) (16) (17) (18)
		Reactions of furyl radicals (cyclo-C 4 H 3 O)
	Isomerization
	Furyl-3 = furyl-2		1.5 × 10 13	0.0	67.8		g	(19)
	Decompositions by β-scission
	Furyl-2 = CHCHCHCO		1.8 × 10 13	0.341	34.511	a	(20)
	Furyl-3 = CHCCHCHO		1.1 × 10 13	0.306	31.215	a	(21)
	Oxidations
	Furyl-2 + O 2 = CH 2 CHCO + CO 2 Furyl-2 + O 2 => CHCHCHO + CO + O Furyl-3 + O 2 = CH 2 CHCO + CO 2 Furyl-3 + O 2 => C 2 H 2 + CHO + CO + O	4.5 × 10 16 3.3 × 10 11 4.5 × 10 16 3.3 × 10 11	-1.39 -0.29 -1.39 -0.29	1.0 10.0 1.0 10.0		h h h h	(22) (23) (24) (25)
	Combinations
	Furyl-2 + O = CHCHCHO + CO	3.0 × 10 13	0.0	0.0		i	(26)
	Furyl-3 + O = CHCHCHO + CO	3.0 × 10 13	0.0	0.0		i	(27)
			Reactions of cyclo-C 4 H 5 O
	Isomerization
	C 4 H 5 O-2 = C 4 H 5 O -3		3.3 × 10 9	1.0	39.1		j	(28)
	Decompositions by β-scission
	C 4 H 5 O-2 = C 2 H 2 +CH 2 CHO C 4 H 5 O-3 = CH 2 CHCHCHO		8.1 × 10 12 1.1 × 10 13	0.246 0.085	34.925 22.555	a a	(29) (30)
	Reactions with O 2 C 4 H 5 O-2 + O 2 = furan + HO 2 C 4 H 5 O-3 + O 2 = furan + HO 2	2.6 × 10 11 1.6 × 10 12	0.0 0.0	2.5 15.2		k k	(31) (32)
	Combinations
	C 4 H 5 O-2 + O = C 2 H 3 CHO + CHO C 4 H 5 O-3 + O = C 2 H 2 + HCHO + CHO C 4 H 5 O-2 + H =furan23H C 4 H 5 O-3 + H = furan23H C 4 H 5 O-3 + H = furan25H	5.5 × 10 13 5.5 × 10 13 1.0 × 10 14 1.0 × 10 14 1.0 × 10 14	0.0 0.0 0.0 0.0 0.0	0.0 0.0 0.0 0.0 0.0		l l b b b	(33) (34) (35) (36) (37)
			Reactions of C 4 H 3 O
	Decompositions by β-scission
	CHCHCHCO = C 2 H 2 + CHCO CHCHCHCO = CHCCHCO + H	5.3 × 10 15 8.7 × 10 11	-0.234 0.706	34.234 38.680	a a	(38) (39)

Table 3

 3 Thermochemical data of furan and its derived species.

		Species	Structure	(kcal/mol) (cal/mol/K)	(T) (cal/mol/K)
								298	400	500	600	800 1000 1500 2000
	Furan	*				-8.23	65.2	15.5 21.1 25.6 29.2 34.3 37.7	42.7	45.4
	Furyl-2	*			59.0	66.5	15.1 20.1 23.9 27.0 31.3 34.2	38.2	40.3
	Furyl-3	*			59.2	66.3	15.0 19.9 23.8 26.9 31.2 34.1	38.1	40.3
	Methylfuran	*		-18.9	73.8	21.2 27.8 33.3 37.9 44.8 49.6	56.6	60.5
	FurylCH 2 *			14.4	72.6	21.3 28.0 33.3 37.4 43.3 47.3	53.2	56.4
	C 4 H 5 O-2			18.1	68.2	21.4 26.1 30.2 33.7 39.1 42.7	43.7	44.7
	C 4 H 5 O-3			8.45	68.1	21.2 25.8 29.9 33.4 38.9 42.5	43.6	44.7
	Furan25H			-14.18	65.1	15.9 23.0 28.9 33.6 40.7 45.7	57.3	61.4
	Furan23H			-21.76	67.0	21.7 27.0 31.6 35.7 42.0 46.2	48.2	50.0
	THF-2yl			-4.09	71.9	24.9 31.1 36.2 40.6 47.2 51.3	53.2	55.0
	THF-3yl			-1.59	71.9	24.9 31.1 36.3 40.6 47.2 51.3	53.2	55.0
	THF					-44.08	69.3	25.2 32.0 37.7 42.6 50.0 54.8	57.7	60.2
	FurylCH 2 O			0.25	82.4	26.4 33.4 38.8 43.0 48.8 53.0	67.3	72.7
	CHCHCHCO *		66.4	74.5	19.7 23.5 26.4 28.7 32.2 34.7	38.4	40.5
	CHCCHCHO		67.0	82.8	19.0 22.3 25.1 27.4 30.9 33.3	37.3	40.3
	CH 2 CHCHCHO		23.9	83.1	20.7 25.1 28.8 31.9 36.7 40.0	45.8	50.1
	CHCHCHCHOH		49.3	84.7	22.0 26.3 29.7 32.6 37.1 40.3	45.8	49.9
	CHCCHCO *		47.1	70.9	18.9 21.9 24.2 25.9 28.6 30.5	33.5	35.1
	CHCHCHO			42.4	72.4	13.9 17.5 20.6 23.4 28.3 32.3	36.7	39.5
	CHCCH 2 CHO		18.7	76.3	20.9 24.4 27.4 30.2 34.5 37.3	38.8	40.7
	CH 2 CHCH 2 CHO		-18.9	80.4	22.0 26.7 30.9 34.5 40.4 44.3	46.3	49.0
	CH 2 CCHCHO (FA) *		15.6	73.4	19.7 24.0 27.6 30.6 35.2 38.5	43.1	45.5
	C 5 H 6 O				14.8	83.3	25.8 31.0 35.5 39.4 45.6 49.5	51.4	53.9
	C 5 H 8 O				-24.5	86.9	26.9 33.3 38.9 43.8 51.5 56.5	58.9	62.2

Table 4

 4 Flame species measured with their respective IEs, maximum mole fractions (X max ), and positions (P, mm) in furan/oxygen/argon flames

	Note: a Refers to [59]; b errors for stable species are ±0.05eV, for radicals are ±0.10eV; c the values are the total
	mole fraction of a specific mass, for example, the mole fraction of mass 44 includes both ethenol and
	acetaldehyde; d refers to [60].

ACKNOWLEDGMENTS This work performed by CNRS-Nancy University was funded by the European Commission through the "Clean ICE" Advanced Research Grant of the European Research Council. FQ is grateful for the funding supports from Chinese Academy of Sciences, Natural Science Foundation of China (50925623), and Ministry of Science and Technology of China (2007DFA61310 This work was granted access to the HPC resources of CINES under the allocation 2010-[c2010086030] made by GENCI (Grand Equipement National de Calcul Intensif).

Note: Thermochemical data have been calculated by software THERGAS [31], except for species marked with * for which they have been obtained from theoretical calculation (see text).

List of Supplemental Material

A "supplemental material" is included in this paper. It consists of the following sections: Section 1: Experimental mole fraction data of furan flames Table S1 is the experimental mole fraction data measured in the three furan flames.

Section 2: The detailed reaction mechanism of this work

The reaction mechanism developed in this study, including 206 species and 1368 reactions, is displayed in this section.

Section 3: Isodesmic reactions

This part gives the isodesmic reactions from which the enthalpies of formation of the species listed in Table 3 are calculated. Section 4: Geometries/Frequencies/Moments of inertia In this section, the geometries, frequencies and moments of inertia of 34 species involved in figures 1, 2 and 3 are listed.