M Longsworth

M Fivel

The effect of stress on the cross-slip energy in face-centered cubic metals: a study using dislocation dynamics simulations and line tension models

Keywords: activation energy, cross-slip, dislocation dynamics, line tension model, Escaig stress, Schmid stress

 were in good quantitative agreement with the simulation results. Hence, both could be used to calculate the activation enthalpy of screw segments in dislocation dynamics simulations.

Introduction

Cross-slip is a thermally-activated process by which screw dislocations change their glide plane, allowing them to overcome obstacles and populate other slip planes [START_REF] Anderson | Theory of Dislocations[END_REF] . It plays an important role in two competing processes during plastic deformation [START_REF] Jackson | The role of cross-slip in the plastic deformation of crystals[END_REF] . On the one hand, it causes dynamic recovery during stage III by dislocation annihilation [START_REF] Hirsch | The Physics of Metals[END_REF] , and other, it contributes to work hardening by dislocation multiplication [START_REF] Jackson | Dislocation modelling of shear in f.c.c. crystals[END_REF][START_REF] Sudmanns | Dislocation multiplication by cross-slip and glissile reaction in a dislocation based continuum formulation of crystal plasticity[END_REF] .

Since cross-slip is a thermally activated process, its rate is controlled by an Arrhenius-like equation. Many dislocation dynamics (DD) simulations have successfully model plastic deformation in crystals [START_REF] Zbib | On plastic deformation and the dynamics of 3D dislocations[END_REF][START_REF] Weygand | Aspects of boundary-value problem solutions with three-dimensional dislocation dynamics[END_REF][START_REF] Verdier | Mesoscopic scale simulation of dislocation dynamics in fcc metals: Principles and applications[END_REF][START_REF] Déprés | 3D Discrete Dislocation Dynamics Investigations of Fatigue Crack Initiation and Propagation[END_REF][START_REF] Chaussidon | Dislocation dynamics simulations of plasticity in Fe laths at low temperature[END_REF] . Most of them use the rate equation proposed by L. Kubin et al. [START_REF] Kubin | Dislocation Microstructures and Plastic Flow: A 3D Simulation[END_REF] . In their model, the cross-slip probability explicitly depends on the resolved shear stress at the onset of stage-III plastic deformation, which is a macroscopic parameter obtained from the stress-strain curve of the material.

This probability law has demonstrated to be robust enough to allow the formation of persistent slip bands in the low-cycle fatigue regime [START_REF] Déprés | Crack initiation in fatigue: experiments and threedimensional dislocation simulations[END_REF] , or to study the recovery dynamics in copper [START_REF] Devincre | Collinear interactions of dislocations and slip systems[END_REF] . However, it is not sufficient in many situations. As an example, in the case of irradiated materials, the local level of stress is very high because of the large number of obstacles associated with the defect density. In this case, the cross-slip events are too frequent [START_REF] Nogaret | Clear band formation simulated by dislocation dynamics: Role of helical turns and pile-ups[END_REF] . Contrarily, when the stress level is very low, the probability of inducing cross-slip events will be very low or nonexistent.

Based on the work of L. Kubin et al. [START_REF] Kubin | Dislocation Microstructures and Plastic Flow: A 3D Simulation[END_REF] , A. Hussein et al. [START_REF] Hussein | Microstructurally based cross-slip mechanisms and their effects on dislocation microstructure evolution in fcc crystals[END_REF] recently proposed a more physical rate equation that depends only on microstructure parameters, e.g. the energy barrier required to form a constriction point on a screw dislocation, the Debye frequency of the material, the activation volume and the local stress. This allowed them to implement DD simulations considering the different cross-slip mechanisms identified from molecular dynamics (MD) simulations. Quantifying the effect of stress on the activation barrier is thus of great importance to generate realistic DD simulations.

The two most plausible cross-slip models for Face-Centered Cubic (FCC) crystals are the Fleischer (FL) and Friedel-Escaig (FE) mechanisms.

The FL model was developed by R. Fleischer [START_REF] Fleischer | Cross slip of extended dislocations[END_REF] . In his model, the leading partial in a stackingfault dissociates into two partials; one glides on the cross-slip plane and the other remains sessile at the intersection of the two planes (stair-rod dislocation). The trailing partial then reacts with the stair-rod dislocation to form another Shockley partial in the conjugate plane, which completes the cross-slip process.

The FE model was developed by B. Escaig [START_REF] Escaig | Sur le glissement dévié des dislocations dans la structure cubique à faces centrées[END_REF] based on J. Friedel's ideas [START_REF] Friedel | Dislocations: International Series of Monographs on Solid State Physics[END_REF] . Moreover, he used the line tension (LT) model of cross-slip developed by A. Stroh [START_REF] Stroh | Constrictions and Jogs in Extended Dislocations[END_REF] to study the effect of the widening stress (Escaig stress) on the energy barrier. In the FE model, the extended dislocation stops moving in the glide plane due to an obstacle. An already-existing constriction then splits into two halves and then separate in the cross-slip plane. This process proceeds spontaneously due to the stress on the cross-slip plane until the dissociation is completed.

Although there have been some controversy on which cross-slip model prevails [START_REF] Püschl | Models for dislocation cross-slip in close-packed crystal structures: a critical review[END_REF] , recent atomistic simulations confirm that both FL and FE mechanisms can occur. The determination of which process would dominate depends on many factors.

C. Jin et al. [START_REF] Jin | Dislocation cross-slip mechanisms in aluminum[END_REF] used the nudged-elastic-band method to study FCC Al. They concluded that there is a critical dislocation length below which the screw dislocation cross-slips only via FL mechanism, whereas longer ones would cross-slip either via FL or FE mechanisms, depending on the initial and final positions of the dislocation. The critical length they found was of the order of 22b, where b is the magnitude of the Burgers vector.

S. Xu et al. [START_REF] Xu | Shear stress-and line length-dependent screw dislocation cross-slip in FCC Ni[END_REF] reported similar findings in FCC Ni using dynamic concurrent atomistic-continuum and MD simulations. They found that screw dislocations of length 6.47 nm cross-slip via FL mechanism, whereas longer dislocations (of length 12.94 nm) can cross-slip via either the FE or FL mechanisms. They also suggested that the critical shear stress for both FE and FL mechanisms depends on the dislocation length.

Recently, W. Kuykendall et al. [START_REF] Kuykendall | Stress effects on the energy barrier and mechanisms of cross-slip in FCC nickel[END_REF] used atomistic simulations to examine the effect of stress on the energy barrier of cross-slip via FE and FL mechanisms in FCC Ni. They concluded that the FE mechanism prevails when the Escaig stress on the glide plane is dominant, and that increasing the Schmid and Escaig stress on the cross-slip plane promotes the FL mechanism.

The material of study in this work is FCC Cu due to its broad literature. The mechanism of cross-slip in Cu has been fully examined using atomistic simulations and it was reported that dislocations in this material cross-slip via FE mechanism [START_REF] Rasmussen | Atomistic Determination of Cross-Slip Pathway and Energetics[END_REF][START_REF] Vegge | Atomistic simulations of crossslip of jogged screw dislocations in copper[END_REF] . The activation energy of Cu has been calculated using both elasticity theory and atomistic simulations. However, the vast majority of energy barrier calculations employ atomistic simulations [START_REF] Vegge | Atomistic simulations of crossslip of jogged screw dislocations in copper[END_REF][START_REF] Rao | Atomistic simulation of cross-slip processes in model fcc structures[END_REF][START_REF] Rao | Atomistic simulations of cross-slip nucleation at screw dislocation intersections in face-centered cubic nickel[END_REF][START_REF] Rao | On the Escaig obstacle hypothesis for cross-slip in face-centered-cubic materials[END_REF][START_REF] Rao | Screw dislocation cross slip at cross-slip plane jogs and screw dipole annihilation in FCC Cu and Ni investigated via atomistic simulations[END_REF][START_REF] Rasmussen | Simulations of the atomic structure, energetics, and cross slip of screw dislocations in copper[END_REF] .

Alternatively, B. Ramírez et al. [START_REF] Ramírez | Ab initio continuum model for the influence of local stress on cross-slip of screw dislocations in fcc metals[END_REF] used a network-based formulation of dislocation dynamics to study the effect of a general stress on the energy barrier of cross-slip via FE mechanism. They found that the stress-free activation energy in Cu is 1.9 eV when the core is represented by two partials, and that it converges to 1.43 eV when the core is distributed by 20 Volterra partial fractional dislocations, the later value being in better agreement with the experimental results. They also found that the cross-slip energy reduces to 0.62 eV in presence of a Lomer-Cottrel junction.

In recent years, K. Kang et al. [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] extended Escaig's LT model to include the effect of gliding stress (Schmid stress) on the cross-slip plane. Contrary to Escaig's claim that the Schmid stress had a negligible effect on the energy barrier, they found that it had a comparable effect to the Escaig stress. They calculated the energy barrier in FCC Ni using atomistic methods at different stress conditions and found good qualitative agreement with their LT model.

More recently, A. Malka-Markovitz et. al. [START_REF] Malka-Markovitz | Cross-slip in face centred cubic metals: a general full stressfield dependent activation energy line-tension model[END_REF] solved Strohs's LT model of cross-slip [START_REF] Stroh | Constrictions and Jogs in Extended Dislocations[END_REF] exactly by linearizing the interaction force between the partials. They obtained a general expression for the activation enthalpy as function of the elastic constants and local stress.

The DD simulations of fatigue mostly rely on the energy barrier, which is often calculated using atomistic simulations. However, the interatomic potential quantitatively affects the simulation results [START_REF] Rao | Atomistic simulation of cross-slip processes in model fcc structures[END_REF] . Hence, they may not be reliable. On the other hand, in DD simulations, the core-width used to compute the elastic energy can be determined using the theory developed by G. Schöck [START_REF] Schöck | Atomic dislocation core parameters[END_REF] , as suggested by B. Ramírez et al. [START_REF] Ramírez | Ab initio continuum model for the influence of local stress on cross-slip of screw dislocations in fcc metals[END_REF] . Moreover, DD simulations naturally enable the possibility to study the effect of complex structures on the energy barrier. In this way, they offer an attractive alternative to large-scale atomic simulations.

The contribution of this work is twofold. First, the effect of stress on the cross-slip activation enthalpy is carefully calculated using DD simulations, which is an interesting alternative to the common atomistic approach found in the literature. Second, the results are compared with those obtained from two recent LT models of cross-slip in order to verify their consistency. The objective is to determine whether the aforementioned LT models could be reliably used in DD simulations to calculate the activation enthalpy of screw segments.

The rest of this paper is organized as follows. Section 2 introduces the DD simulation methodology, K. Kang's et al. [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] LT model of cross-slip and the general expression for the activation enthalpy proposed by A. Malka-Markovitz et. al. [START_REF] Malka-Markovitz | Cross-slip in face centred cubic metals: a general full stressfield dependent activation energy line-tension model[END_REF] . Section 3 describes the effect of stress on the activation enthalpy obtained using DD simulations and the two LT models. Finally, Section 4 presents the summary and conclusions of this work.

Methodology

DD simulations 2.1.1. Force computation

Dislocation dynamics is a modeling approach to study crystal plasticity, in which dislocation lines are discretized in linear segments and their motion is determined [START_REF] Sills | Fundamentals of Dislocation Dynamics Simulations[END_REF] . A segment is defined by two nodes on its ends and a tangent line ξ that connects them. The Burgers vector b associated to it indicates the displacement's magnitude and direction as the segment moves.

The force exerted by the local stress at each node is calculated using the Peach-Köhler formula. Similarly to other discrete DD codes [START_REF]The ParaDiS program[END_REF] , NuMoDis [START_REF]Numerical Modelling of Dislocations (NuMoDis[END_REF] uses the formalism developed by A. Arsenlis et al. [START_REF] Arsenlis | Enabling strain hardening simulations with dislocation dynamics[END_REF] to compute the nodal forces. This force affects the nodal velocity through a given mobility law. Thus, the dislocation dynamics can be determined by updating the position of all nodes after a simulation step.

In order to study cross-slip in FCC crystals, the dissociated character of dislocations must be considered. The SF produces a force that tends to attract the partials towards each other, compensating the elastic repulsion between them. This results in an equilibrium separation distance in the SF ribbon.

Perfect dislocations dissociate into partials if it decreases its elastic energy. However, discrete DD codes are based on nodal forces and do not consider the systems energetics to determine its evolution. In consequence, E. Martinez et al. [START_REF] Martínez | Atomistically informed dislocation dynamics in fcc crystals[END_REF] proposed to account for the intrinsic SF energy γ i by introducing a SF force (f sf ). They derived that f sf should be of magnitude γ i , and perpendicular to both the tangent line and SF plane n. The equation satisfying these conditions is given by:

f sf = γ i n × ξ (1)
In NuMoDis [START_REF]Numerical Modelling of Dislocations (NuMoDis[END_REF] , the SF ribbon is defined by an ordered list of nodes conforming the dislocation partials and the plane that contains them. The SF forces acting on a partial, in order to be physical, must point towards the opposite partial. Thus, given a slip plane orientation, the node order defining the circuit must be such that the SF force given by Eq. 1 satisfy this condition.

Velocity computation

In FCC metals, the relationship between forces and velocities is linear due to the small intrinsic lattice resistance [START_REF] Bulatov | Computer Simulations of Dislocations[END_REF] . The total force at a node f is given by f = f P K + f sf , where f P K is the Peach-Köhler force. The projections of f P K and f sf on the glide plane act perpendicular to the line segment. Therefore, the mobility law is of the form:

f = B • v (2)
where v is the nodal velocity and B the viscosity matrix:

B( ξ) = B(I -ξ ⊗ ξ) ( 3
)
where B is the viscosity coefficient. Thus, the nodal velocity is obtained by inverting the viscosity matrix from Eq. 2:

v = B -1 • f (4)
At each simulation step, the force at every node due to the influence of all other segments and SF is computed; the nodal velocities are thus calculated and their new positions are updated.

Elastic energy

W. Cai et al. [START_REF] Cai | A non-singular continuum theory of dislocations[END_REF] showed that the non-singular interaction energy E ns between a segment with ends x 1 and x 2 and Burgers vector b with another segment with ends x 3 and x 4 and Burgers vector b is given by:

E ns = - µ 8π x 4 x 3 x 2 x 1 ∂ k ∂ k R a b i b j dx i dx j - µ 4π(1 -ν) x 4 x 3 x 2 x 1 ∂ i ∂ j R a b i b j dx k dx k + µ 4π(1 -ν) x 4 x 3 x 2 x 1 ∂ k ∂ k R a b i b i dx j dx j - µν 4π(1 -ν) x 4 x 3 x 2 x 1 ∂ k ∂ k R a b i b j dx j dx i (5)
where µ is the shear modulus, ν the Poisson's ratio,

R a = √ R • R + a 2 , R = x -x
, and a the core-width parameter.

The analytical expressions resulting from the solution of Eq. 5 for parallel and non-parallel segments were provided by W. Cai et al. [START_REF] Cai | A non-singular continuum theory of dislocations[END_REF] . These formulas were used to calculate the interaction and self-energies of the segments in NuMoDis [START_REF]Numerical Modelling of Dislocations (NuMoDis[END_REF] .

The core-width parameter a affects the elastic energy both qualitatively and quantitatively. G. Schöck [START_REF] Schöck | Atomic dislocation core parameters[END_REF] developed a theory based on the Peierls framework [START_REF] Schoeck | The Peierls model: Progress and limitations[END_REF] to obtain a physical value of a. He showed that the total energy of a single straight dislocation is given by [START_REF] Schöck | Atomic dislocation core parameters[END_REF] :

E P = E L (φ) ln r 2w(φ) + 1 ( 6
)
where φ is the angle between the Burgers vector and the line direction, r is the outer cut-off radius, w(φ) is the average width of its edge and screw components, and E L (φ) is the pre-logarithmic elastic line energy factor:

E L (φ) = µb 2 4π cos 2 (φ) + sin 2 (φ) 1 -ν (7)
G. Schöck [START_REF] Schöck | Atomic dislocation core parameters[END_REF] argued that w(φ) can be approximated as the width of a single dislocation in an isotropic medium passing through two equilibrium positions i.e. with SF energy γ = γ u [START_REF] Schöck | Atomic dislocation core parameters[END_REF] :

w(φ) b = cos 2 (φ) + sin 2 (φ)/(1 -ν) b/c + 2π 2 γ u /(µb) ( 8
)
where c is the interplanar glide plane spacing, γ u the unstable SF energy and b is the magnitude of the Burgers vector.

Using W. Cai's et al. [START_REF] Cai | A non-singular continuum theory of dislocations[END_REF] model, the line energy for screw E S and edge dislocation E E are found to be [START_REF] Schöck | Atomic dislocation core parameters[END_REF] :

E S = µb 2 4π ln r √ e a S (9
)

E E = µb 2 4π(1 -ν) ln r a E (10)
where a S and a E are the core spreading parameters of the screw and edge components, respectively, and e is the base of the natural logarithm.

The continuum description in Eqs. 9 and 10 must be consistent with the Peierl's atomistic treatment in eq. 6. Imposing this condition, a S and a E are required to be [START_REF] Schöck | Atomic dislocation core parameters[END_REF] :

a S = 2w(0) √ e and a E = 2w(π/2) e (11)
Following the work of B. Ramírez et al. [START_REF] Ramírez | Ab initio continuum model for the influence of local stress on cross-slip of screw dislocations in fcc metals[END_REF] , the core-width parameter can be considered as an averaged value of a S and a E :

a ≈ 2w(0) e = 2 e b b/c + 2π 2 γ u /(µb) (12)

Effective stacking-fault energy

Dislocations in FCC crystals dissociate into partials. The dislocation moves due to the action of an external stress. The partial at the front respect to the direction of motion is called the "leading" partial, and the other "trailing" partial. Let L = { P 1 , P 2 , ..., P N } and L = { Q 1 , Q 2 , ..., Q N } represent the set of nodes conforming the leading and trailing partials, respectively (see Fig. 1).

Figure 1: The solid dots denote the nodal positions; the shape of the leading and trailing partials are described by L and L, respectively.

The effective SF energy is given by the SF energy W SF plus the work done by the applied forces W app , which can be calculated as follows:

W SF = γ i N j=1 y j+1 + y j 2 ∆x j - ỹj+1 + ỹj 2 ∆x j (13)
W app = N j=1 F y j ∆x j + F ỹj ∆x j (14)
where:

y j = P j • û, ỹj = Q j • û, ∆x j = ( P j+1 -P j ) • b, ∆x j = ( Q j+1 -Q j ) • b for j = 1, 2, ..., N (15)
û = n × b is an unitary vector contained in the slip plane of normal n, perpendicular to the normalized Burgers vector b. In absence of Schmid stress, the partials are indistinguishable from each other due to the spatial symmetry with respect to the x axis. For that reason, the sign of û does not affect the energy. However, the leading and trailing partials are univocally defined in the presence of a Schmid stress. In that case, û must point towards the bowing-out direction because the energy barrier does not depend on the sign of the Schmid stress [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF][START_REF] Malka-Markovitz | Cross-slip in face centred cubic metals: a general full stressfield dependent activation energy line-tension model[END_REF][START_REF] Kuykendall | Stress effects on the energy barrier and mechanisms of cross-slip in FCC nickel[END_REF] .

The proportionality factors F and F are given by:

F = - bσ E 2 √ 3 - bσ S 2 and F = bσ E 2 √ 3 - bσ S 2 (16)
where the subscripts E and S refer to the Escaig and Schmid stress, respectively. All the nodes were mobile except those two at the extremes, which simulate the constriction points at the equilibrium configurations. Therefore, the partials always satisfied y 1 = ỹ1 = y N = ỹN = 0. Figure 2: An equilibrium configuration during cross-slip via FE mechanism. The process is completed when the constriction separation tends to infinity. In this limit, the dislocations is completely contained in the conjugate plane with separation between the partials l cs for all x (figure adapted from K. Kang et al. [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] ).

Cross-slip simulations

Cross-slip via FE mechanism was assumed in all simulations. The DD code was used to find the equilibrium configurations at different cross-slip stages. The total energy of a configuration consisted of its total elastic energy plus the effective stacking-fault energy, which were computed as described in sections 2.1.3 and 2.1.4. The applied stresses were small enough such that the Frank-Read source remained deactivated in all simulations.

According to the FE mechanism, the dislocation must stop moving before cross-slipping. Therefore, the Schmid stress on the glide plane σ g S was always assumed to be zero. The dislocation width in the glide plane was affected by the Escaig stress σ g E . Soon after stop moving, an alreadyexisting constriction splits into two halves and the dislocation expands in the cross-slip plane until equilibrium is reestablished (see Fig. 2). This occurs in the limit when the distance between the constriction points in the cross-slip plane 2x d tends to infinity.

In DD codes, an infinite dislocation does not really exist. Thus, the initial configuration consisted of a dislocation with length of 1000 Å lying in the (111) plane. The dislocation was pinned at its ends. The Burgers vector of the partials pointed away from each other as shown in Fig. 3(a). This was a long-enough dislocation such that most part of it attained equilibrium i.e the distance l g between the partials.

An intermediate cross-slip configuration consisted of a semi-infinite dislocation on the glide plane with a constricted cross-slipped segment on the conjugate plane (111). The Burgers vector of the partials pointed towards each other on the cross-slip plane. The total configuration length was fixed to 1000 Å, in agreement with the initial state. The intermediate cross-slip configuration was thus characterized by the distance between its constriction points.

In the limiting case when the constrictions points are infinitely apart, the dislocation must be completely contained in the cross-slip plane (see 3(b)). Unlike in the glide plane, the dislocation was allowed to move by bowing-out on the conjugate plane. The sign of the Schmid stress on the cross-slip plane σ cs S was affected only the bowing-out direction, and not the total energy. As in the glide plane, the width of dislocation in the conjugate plane was affected by the Escaig stress σ cs E . Figure 4: The shape of the partials in the cross-slip plane are described by the y 1 (x) and y 2 (x), and satisfy the boundary conditions y 1 (±x d ) = y 2 (±x d ) = 0 (figure adapted from K. Kang et al. [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] ).

In order to study the effect of σ g E , σ cs S and σ cs E , the coordinate system was rotated such that the [111] direction matched with the ẑ axis and the Burgers vector [110] with the x axis. Thus, the stress on the glide and cross-slip planes as function of the tensor components in the Cartesian coordinate system were given by:

σ g S = σ xz (17)
σ g E = σ yz (18)

σ cs S = 2 √ 2σ xy -σ xz 3 ( 19
)
σ cs E = 7σ zy + 2 √ 2(σ yy -σ zz ) 9 (20)
The energy of an intermediate cross-slip configuration was measured with respect to the energy of the initial state i.e. the energy of the dislocation resting at equilibrium in the glide plane under stress σ g E . Thus, the energy barrier was defined as the maximum energy among all intermediate cross-slip configurations as function of the distance between the constriction points.

Line tension model

In K. Kang's et al. [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] LT model of cross-slip, the FE mechanism is assumed. Thus, the dislocation does not move before cross-slipping i.e. the Schmid stress on the glide plane is zero (σ g S = 0). However, the Escaig stress on the glide plane σ g E is not necessarily zero. The equilibrium separation between the partials in the glide plane l g depends on σ g E :

l g = A F g (21)
where for a screw dislocation dissociated into two 30 • partials [START_REF] Cai | Intrinsic Mobility of a Dissociated Dislocation in Silicon[END_REF] :

A = 1 4 - 1 12(1 -ν) µb 2 2π ( 22
)
F g = γ i -bσ g E /(2 √ 3) (23)
The convention used in this framework is that positive σ g E increases the distance between the partials and negative decreases it.

Analogously, if the dislocation does not move in the cross-slip plane, its equilibrium distance l cs is given by:

l cs = A F cs (24)
where

F cs = γ i -bσ cs E /(2 √ 3) (25)
The shape of the partials is described by continuous functions y 1 (x) and y 2 (x) (see Fig. 4). The constriction points are located in the x axis at ±x d and thus y 1 (±x d ) = y 2 (±x d ) = 0 must be satisfied.

The energy contribution from the cross-slipped segment (measured with respect to the energy of an infinite straight screw dislocation resting at equilibrium in the glide plane) is given by [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] :

E cs = W 1 + W 2 + W 3 ( 26
)
where W 1 is the change in elastic repulsion between the partials, W 2 is the change of the line energy and W 3 the change in SF energy and the work done by the stress. Specifically,

W 1 = A x d -x d ln l g y 1 -y 2 dx (27)
W 2 = T x d -x d 1 + (y 1 ) 2 + 1 + (y 2 ) 2 dx -4T x d (28)
W 3 = x d -x d (F cs 1 y 1 + F cs 2 y 2 ) dx -2F g l g x d (29)
where

F cs 1 = F cs -bσ cs S /2 (30)
F cs 2 = -F cs -bσ cs S /2 (31)
The energy functional in Eq. 26 was minimized using the conjugate gradient method as suggested by Kang et al. [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] .

The line energy per unit length was given by [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] :

T = αµb 2 /2 (32)
The variable α was treated as a fitting parameter to match the DD simulations of cross-slip (see section 3).

Formally, the energy barrier is given by:

E b (σ g E , σ cs S , σ cs E ) = max x d E(x d ; σ g E , σ cs S , σ cs E ) (33)
where

E = min y 1 ,y 2 E cs (y 1 , y 2 ; x d ; σ g E , σ cs S , σ cs E ) + E g (σ g E ) (34)
E g being the energy of an isolated constriction on a perfect screw dislocation:

E g = 2l g √ AT 1 0 [-ln t + (t -1)] 1/2 dt (35)
The activation energy can be obtained for a given stress condition by solving Eq. 34 numerically.

Another useful formulation of the LT model of cross-slip was proposed by A. Malka-Markovitz et. al. [START_REF] Malka-Markovitz | Cross-slip in face centred cubic metals: a general full stressfield dependent activation energy line-tension model[END_REF] , who solved Strohs's LT model [START_REF] Stroh | Constrictions and Jogs in Extended Dislocations[END_REF] exactly by linearizing the interaction force between the partials. They obtained the following general expression for the activation enthalpy as function of the elastic constants and local-stress components:

∆H(σ) ∆E 0 = β g 2 + β cs 2 tanh(l C ) - β cs 2 2α L 1.55 l C E * + α 3 Ls 1.55 δ β cs l 3 C ( 36
)
where ∆E 0 = 1.9 eV is the energy barrier of copper obtained in this work using DD simulations and α L ≈ 0.6 is a scale factor.

β g and β cs are functions of the Escaig stress: Property Value a 0 [START_REF] Davey | Precision Measurements of the Lattice Constants of Twelve Common Metals[END_REF] 3.61 Å b [START_REF] Kittel | Introduction to Solid State Physics[END_REF] a 0 / √ 2 γ i [START_REF] Bonneville | Cross-slipping process and the stress-orientation dependence in pure copper[END_REF] 42 mJ/m 2 γ u [START_REF] Ramírez | Ab initio continuum model for the influence of local stress on cross-slip of screw dislocations in fcc metals[END_REF] 182 mJ/m 2 ν [START_REF] Ledbetter | Elastic Properties of Metals and Alloys. II. Copper[END_REF] 0.324 µ [START_REF] Schmauder | Micromechanics and Nanosimulation of Metals and Composites: Advanced Methods and Theoretical Concepts[END_REF] 54.6 GPa B [START_REF] Philibert | Dislocations et déformation plastique[END_REF] 1.5 ×10 -5 Pa-s Table 1: Parameters used for FCC Cu; a 0 is the lattice parameter, b the magnitude of the Burgers vector, γ i the intrinsic SF energy, γ u the unstable SF energy, ν the Poisson's ratio, µ the shear modulus and B the viscosity coefficient at room temperature.

β g = 1 1 + √ 3b 6γ i σ g E and β cs = 1 1 + √ 3b 6γ i σ cs E ( 37
)
whereas δ is a quadratic function of the Schmid stress on the cross-slip plane:

δ = 1 6 bσ cs S γ i 2 (38)
The E * function is a dimensionless off-set in the interaction energy between the glide and crossslip planes [START_REF] Malka-Markovitz | Cross-slip in face-centered cubic metals: a general Escaig stress-dependent activation energy line tension model[END_REF] :

E * = ln β cs β g (39)
The critical length l C in Eq. 36 can be found by solving the following equation numerically [START_REF] Malka-Markovitz | Cross-slip in face-centered cubic metals: a general Escaig stress-dependent activation energy line tension model[END_REF] :

1.55 cosh 2 (l C ) - 3δ β cs l 2 C = 2E * (40)
In contrast with the convention used in K. Kang's et al. [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] LT model of cross-slip, positive σ g E decreases the distance between the partials and negative increases it.

In this work, the sign convention of K. Kang's et al. [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] was used in the plots.

The elastic constants of Cu used in both the DD simulations and LT model are shown in Table 1.

Results and discussion

The calibration of the line tension model

The relative energy of the intermediate cross-slip configurations at zero stress is shown in Fig. 5. The energy uncertainty of the order of ±0.1 eV in all the simulations results can be attributed to the line discretization used in the DD simulations, as well as to the finite core-width size of the nonsingular theory of dislocations. The parameter α in K. Kang's et al. [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] LT model was calibrated such that the equilibration energy matched with the DD simulations at zero stress. The value was found to be α = 0.22. Together with the core-width parameter in the elastic energy, the choice of α was a key parameter for the good agreement between K. Kang's et al. [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] LT model and DD simulations. The effect of compressive Escaig stress was to decrease the SF width. The stressed equilibrium separation obtained from the DD simulations was in great agreement with the LT model, as shown in Fig. 6.

The energy predicted by the K. Kang's et al. [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] LT model and that obtained with the simulations does not match for short cross-slipped segments as seen in Fig. 5. This is expected because the LT model approximates the energy of a differential segment as that of a straight dislocation [START_REF] Dupuy | A study of dislocation junctions in FCC metals by an orientation dependent line tension model[END_REF] . In fact, the LT model breaks down in the limiting case were 2x d = 0 because the elastic repulsion diverges (see Eq. 27). However, this is not relevant since the objective of this work is to quantify the stress effect on the energy barrier, which occurs for longer dislocation lines where the LT model is more physical. In the unstressed case, the energy monotonically augments as the constriction separation increases, but quickly reaches a maximum value of 1.9 eV. This means that once the cross-slip segment reaches a critical length (around 60 Å), no further energy is required to be invested and the cross-slip process continues spontaneously. Therefore, the unstressed energy barrier of copper is equal to 1.9 eV, as reported by B. Ramírez et al. [START_REF] Ramírez | Ab initio continuum model for the influence of local stress on cross-slip of screw dislocations in fcc metals[END_REF] .

The equilibration energy encountered in the unstressed cross-slip is roughly equal to twice the energy required to form an isolated constriction. It can be argued that at the beginning of cross-slip, when the two constrictions are close to each other, the elastic repulsion between them adds up to the energy required to form them in isolation. However, as the distance between them increases, their interaction further decreases. Thus, if they are sufficiently far apart, the total energy with respect to the relaxed dislocation is roughly equal to that of the two independent constrictions. The planar constriction energy at zero stress is 0.7 eV, according to the DD simulations (see Fig. 7).

A. Stroh et al. [START_REF] Stroh | Constrictions and Jogs in Extended Dislocations[END_REF] found the energy of a constriction using elasticity theory. Assuming a shear modulus of G = 45 GPa and an equilibration distance between the partials at zero stress of d 0 = 18 Å, one finds that in his model the energy of a single constriction is roughly 2 eV, which is almost three times higher than that found in this work (0.7 eV).

W. Püschl [START_REF] Püschl | The Energy of Constrictions in Extended Dislocations[END_REF] estimated a constriction energy of 1.1 eV for FCC Cu using elasticity theory; G. Saada [START_REF] Saada | Cross-slip and work hardening of f.c.c. crystals[END_REF] generalized A. Strohs's LT model for arbitrary constricting distance. He obtained different constriction energy values in the range of 0.33 -0.7 eV, depending on the choice r 0 /b, where r 0 is a cut-off radius. Although these values are in better agreement with the experimental energy barrier, the (arbitrary) choice of r 0 was not clearly determined.

T. Rasmussen et al. [START_REF] Rasmussen | Simulations of the atomic structure, energetics, and cross slip of screw dislocations in copper[END_REF] found a constriction energy of 1.6 eV using atomistic simulations, in reasonable agreement with A. Stroh et al. [START_REF] Stroh | Constrictions and Jogs in Extended Dislocations[END_REF] , but too high according to the experimental energy barrier.

More recently, B. Ramírez et al. [START_REF] Ramírez | Ab initio continuum model for the influence of local stress on cross-slip of screw dislocations in fcc metals[END_REF] found a constriction energy of 0.7 eV using dislocation dynamics, in excellent agreement with this work.

The effect of pure compressing stress on the constriction energy was found to be negligible. The slight decrease on the constriction energy observed in Fig. 7 is not conclusive. Nevertheless, one can assert that the average constriction energy did not vary much in the given range of stress.

The small effect of stress on the constriction energy can be explained by using the model developed by W. Püschl [START_REF] Püschl | The Energy of Constrictions in Extended Dislocations[END_REF] . It predicts a marginal decrease in constriction energy (of less than 0.2 eV) when the splitting width between the partials decreases from 8b ≈ 20.4 Å to 6b ≈ 15.3 Å. These are roughly the equilibration distances simulated by NuMoDis [START_REF]Numerical Modelling of Dislocations (NuMoDis[END_REF] at zero (17.4 Å) and -250 MPa (12.1 Å), respectively. Thus, according to this model, one should not expect the constriction energy to vary much by applying a pure compressing stress for the given range of values.

The asymptotic convergence of the cross-slip energy barrier at zero stress observed in Fig. 5 has been thoroughly reported before. M. S. Duesbery et al [START_REF] Duesbery | The mechanics and energetics of cross-slip[END_REF] calculated the constriction pair energy as function of the stackingfault energy using a linear-elastic framework, and found that for constriction separation larger than ∼ 50b, the total energy becomes that of two independent constrictions. In particular, for FCC Cu, they found that the asymptotic value of the cross-slipped configuration was 3.7 eV. They also identified the asymmetric nature of the two constrictions in the cross-slip plane, as they asserted that one had lower energy being screw-like, and the other higher energy being edge-like. These are different from the Stroh-type constriction, which are edge-like on the one side, and screw-like on the other. The effect of stress on the energy barrier: (A) σ g S = σ cs E = σ cs S = 0, (B) σ g S = σ g E = σ cs S = 0 and (C) σ g S = σ g E = σ cs E = 0. The green circles represent the DD simulation results, the blue diamonds correspond to the numerical solution of K. Kang's et al. [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] LT model of cross-slip (see Eq. 34) and the red crosses were obtained from the general expression for the activation enthalpy proposed by A. Malka-Markovitz et. al. [START_REF] Malka-Markovitz | Cross-slip in face centred cubic metals: a general full stressfield dependent activation energy line-tension model[END_REF] (see Eq. 36).

Later on, T. Rasmussen et al. [START_REF] Rasmussen | Simulations of the atomic structure, energetics, and cross slip of screw dislocations in copper[END_REF] corroborated several results of M. S. Duesbery et al. [START_REF] Duesbery | The mechanics and energetics of cross-slip[END_REF] using atomistic simulations, as they obtained that indeed the total energy converged to the value of two independent constrictions for separations larger than ∼ 50b. They measured a screw-like constriction energy of -3.8 eV, and an edge-like constriction energy of 1.1 eV, which added up to the total energy barrier in FCC Cu (2.7 eV). Although this value is lower than that obtained by M. S. Duesbery et al [START_REF] Duesbery | The mechanics and energetics of cross-slip[END_REF] , it is roughly the double of the experimental value obtained by J. Bonneville et al. [START_REF] Bonneville | A study of cross-slip activation parameters in pure copper[END_REF] . S. I. Rao et al. [START_REF] Rao | Atomistic simulation of cross-slip processes in model fcc structures[END_REF] obtained an energy barrier for cross-slip at zero stress in the range of 1.07 -1.28 eV using molecular statics. This energy is significantly smaller than that predicted by T. Rasmussen et al. [START_REF] Rasmussen | Simulations of the atomic structure, energetics, and cross slip of screw dislocations in copper[END_REF] . They argued that it could be due to the difference in interatomic potentials used, as well as the Green's function technique used to relax the boundary forces in the simulations developed by S. I. Rao et al.

B. Ramírez et al. [START_REF] Ramírez | Ab initio continuum model for the influence of local stress on cross-slip of screw dislocations in fcc metals[END_REF] found an unstressed energy barrier of 1.9 eV using DD simulations for dislocations dissociated into two partials. They obtained that the energy plateau was reached for a constriction separation of about 60 Å. Both results are in excellent agreement with the simulations in this work. Furthermore, notice that the critical length is about the half of the value found by previous methods.

The total energy of a solid crystal is all potential energy (assuming of course that the crystal does not move). The potential energy of an elastic system is the difference between the internal energy Π int and the energy due to external forces Π ext [START_REF] Kim | Introduction to Nonlinear Finite Element Analysis[END_REF] . The Π int consists of the stress-strain energy due to the stress field introduced by the dislocations (and their interactions) respect to the initial geometry i.e. the perfect crystal. The Π ext consists of the work done by the applied forces, which is obtained by multiplying the displacement generated by the dislocations as they move with the applied forces. Therefore, the energy barrier can be decreased in presence of an external stress.

In many DD fatigue simulations [START_REF] Zbib | On plastic deformation and the dynamics of 3D dislocations[END_REF][START_REF] Weygand | Aspects of boundary-value problem solutions with three-dimensional dislocation dynamics[END_REF][START_REF] Verdier | Mesoscopic scale simulation of dislocation dynamics in fcc metals: Principles and applications[END_REF][START_REF] Déprés | 3D Discrete Dislocation Dynamics Investigations of Fatigue Crack Initiation and Propagation[END_REF][START_REF] Chaussidon | Dislocation dynamics simulations of plasticity in Fe laths at low temperature[END_REF] , following the pioneering work of L. Kubin et al. [START_REF] Kubin | Dislocation Microstructures and Plastic Flow: A 3D Simulation[END_REF] , the driving force for cross-slip is the Schmid stress. However, the screw dislocation does not really move because an obstacle is assumed to stop its motion in the glide plane. Thus, the effect of σ g S is to compress the dislocation against the obstacle, reducing the equilibration distance between the partials [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] . In this way, the σ g S stress plays a similar role to σ g E , if the dislocations is assumed to stop completely by being compressed against an obstacle. For this reason, σ g S is assumed to be zero in all simulations, following the model developed by B. Escaig [START_REF] Escaig | Sur le glissement dévié des dislocations dans la structure cubique à faces centrées[END_REF] .

Uncoupled stress

In order to study the independent influence of the stress components on the energy barrier, three cross-slip possibilities were considered: (A) pure compression on the glide plane, (B) pure expansion on the conjugate plane and (C) pure bowing-out on the conjugate plane. The stress was increased up to 250 MPa, because the equilibrium configuration breaks down around 260 MPa [START_REF] Davis | Loading estimates of lead, copper, cadmium, and zinc in urban runoff from specific sources[END_REF] .

The activation enthalpy obtained from DD simulations, the numerical solution of K. Kang's et al. [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] LT model of cross-slip and the general expression for the activation enthalpy proposed by A. Malka-Markovitz et. al. [START_REF] Malka-Markovitz | Cross-slip in face centred cubic metals: a general full stressfield dependent activation energy line-tension model[END_REF] are in good quantitative agreement with each other, as shown in Fig. 8.

The largest influence on the energy barrier was obtained by pure compression on the glide plane, which reduced it up to 45% (see Fig. 8 (A)). The minimum energy barrier of 1 eV was obtained by applying a compressive stress of 250 MPa on the glide plane. This is in contrast with the lower energy value of 0.95 eV obtained by B. Ramírez et al. [START_REF] Ramírez | Ab initio continuum model for the influence of local stress on cross-slip of screw dislocations in fcc metals[END_REF] under a pure σ cs E of 220 MPa. On the other hand, K. Kang et al. [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] also concluded that the compressive stress applied on the glide was the most effective way of reducing the energy barrier in FCC nickel. The pure widening stress on the cross-slip reduced it up to 20% (see Fig. 8 (B)) and the pure bowing-out stress on the cross-slip produced a decrease in the energy barrier of about 15% (see Fig. 8 (C)).

The effect of stress was not only to decrease the energy barrier but also the activation length of cross-slip. The maximum of energy at 200 MPa occurred for a constriction separation of about 25 Å. This value is comparable to the activation length of 12b ≈ 30.6 Å at 220 MPa obtained by B. Ramírez et al. [START_REF] Ramírez | Ab initio continuum model for the influence of local stress on cross-slip of screw dislocations in fcc metals[END_REF] .

Coupled stress

In order to study the combined effect of stress on the energy barrier, three cross-slip possibilities were considered: (A) compression on the glide plane and expansion on the conjugate plane, (B) compression on the glide plane and bowing-out on the conjugate plane and (C) expansion and bowing-out on the conjugate plane.

As in the previous case, the activation enthalpy obtained from DD simulations and the two formulations of the LT model of cross-slip analyzed in this work [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF][START_REF] Malka-Markovitz | Cross-slip in face centred cubic metals: a general full stressfield dependent activation energy line-tension model[END_REF] are in good agreement (see Fig. 9).

The effect of applying a coupling compressive stress on the glide plane lead to the largest decrease in the energy barrier (see Figs. 9 (A) and(B)). These results resemble the effect applying an uncoupled compressive stress on the glide plane (see Fig. 8 (A)), which also lead to an energy barrier reduction of roughly 50%.

The effect of applying coupled stresses on the conjugate plane only (see Fig. 9 (C)) was comparable of applying an uncoupled Schmid stress on the conjugate plane (see Fig. 8 (C)). In both cases, the energy barrier decreased up to 20%.

These observations are validated by the atomistic simulations of G. Esteban-Manzanares et. al. [START_REF] Esteban-Manzanares | Influence of the stress state on the cross-slip free energy barrier in Al: An atomistic investigation[END_REF] , who studied the effect of stress in aluminum using molecular dynamics and obtained the same qualitative results.

Summary and Conclusions

In this work, the unstressed constriction energy was equal to 0.7 eV, as reported by B. Ramírez et al. [START_REF] Ramírez | Ab initio continuum model for the influence of local stress on cross-slip of screw dislocations in fcc metals[END_REF] . The effect of pure compressive stress in the range of -50 to -250 MPa on the constriction energy was negligible, in agreement with W. Püschl [START_REF] Püschl | The Energy of Constrictions in Extended Dislocations[END_REF] .

The unstressed energy barrier was found to be 1.9 eV, as obtained by B. Ramírez et al. [START_REF] Ramírez | Ab initio continuum model for the influence of local stress on cross-slip of screw dislocations in fcc metals[END_REF] . However, it does not lie within the experimental range of 1.15±0.37 eV measured by J. Bonneville et al. [START_REF] Bonneville | A study of cross-slip activation parameters in pure copper[END_REF] . On the other hand, the energy barrier could be reduced by applying an external stress, leading to an energy barrier in better agreement with the experimental results. The most effective way of reducing it, was by pure compression on the glide plane. Furthermore, it was found using DD simulations, that both the Schmid and Escaig stress on the conjugate plane have a comparable effect in reducing the energy barrier, in qualitative agreement with the atomistic simulations performed by K. Kang et al. [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] in FCC nickel. The energy barrier of FCC Cu at zero stress has been obtained using the line LT model and atomistic simulations before. Some of them reported an energy barrier in better agreement with experiments [START_REF] Rao | Atomistic simulation of cross-slip processes in model fcc structures[END_REF][START_REF] Saada | Cross-slip and work hardening of f.c.c. crystals[END_REF] . However, the LT model suffers from arbitrariness in the cut-off radius determination [START_REF] Saada | Cross-slip and work hardening of f.c.c. crystals[END_REF] and the atomistic simulations are highly sensitive to the interatomic potential used [START_REF] Rao | Atomistic simulation of cross-slip processes in model fcc structures[END_REF] . In DD simulations, these problems do not exist and the core-width used to compute the elastic energy can be calculated using the theory developed by G. Schöck [START_REF] Schöck | Atomic dislocation core parameters[END_REF] as suggested by B. Ramírez et al. [START_REF] Ramírez | Ab initio continuum model for the influence of local stress on cross-slip of screw dislocations in fcc metals[END_REF] . Moreover, the effect of non-homogeneous stress on the energy barrier, like those generated by pile-ups and forest dislocations, can be more naturally studied using DD simulations, in comparison with other simulation techniques.

The activation enthalpy obtained from DD simulations was in excellent agreement with both the numerical solution of K. Kang's et al. [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] LT model of cross-slip and the general expression for the activation enthalpy proposed by A. Malka-Markovitz et. al. [START_REF] Malka-Markovitz | Cross-slip in face centred cubic metals: a general full stressfield dependent activation energy line-tension model[END_REF] . Since the later is an analytical model, it can be easily implemented in a DD code to calculate the activation enthalpy of screw segments.

Acknowledging a limitation of this work, a perfect screw dislocation in FCC crystals can be dissociated in a set of fractional partial dislocations (FPD) pairs, instead of two Shockley partials with Burgers vectors b 1 and b 2 . Based on the work of S. Banerjee et al. [START_REF] Banerjee | Non-singular descriptions of dislocation cores: a hybrid ab initio continuum approach[END_REF] , B. R. Ramírez et al. [START_REF] Ramírez | Ab initio continuum model for the influence of local stress on cross-slip of screw dislocations in fcc metals[END_REF] generalized the dissociation of perfect screw dislocations into n pairs of fractional partial dislocations (FPD) with Burgers vectors b 1 /n and b 2 /n. By using 20 FPDs, they obtained an energy barrier reduction from 1.9 to 1.43 eV at zero stress, which lies within the experimental error of the energy barrier measured by J. Bonneville et al. [START_REF] Bonneville | A study of cross-slip activation parameters in pure copper[END_REF] .

Regarding the LT model, the constant line energy per unit length T makes it impossible to distinguish between edge-like and screw-like constrictions. Thus, a more physical description of cross-slip can be made by letting the total energy depend on the angle between the Burgers vector and the local tangent vector [START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF][START_REF] Malka-Markovitz | Cross-slip in face-centered cubic metals: a general Escaig stress-dependent activation energy line tension model[END_REF][START_REF] Dupuy | A study of dislocation junctions in FCC metals by an orientation dependent line tension model[END_REF] . In spite of it, the energy barrier calculated using the LT model was found to be in good quantitative agreement with the DD simulations, at least for the case of dissociation into two partials. Hence, the orientation-dependent LT model might be in better agreement with DD simulations that describe the dislocation core using more than a pair of FPDs.

Figure 3 :a 0 6 [

 36 Figure 3: Simulation setup. (a) The reference configuration consists of a dissociated right-handed screw (RHS) dislocation at equilibrium on the glide plane (111) with partials' Burgers vector -→ b g 1 = a 0 6 [121] and -→ b g 2 = a 0 6 [211] pointing away from each other. (b) The final configuration consists of a dissociated RHS dislocation at equilibrium in the cross-slip plane (111) with partials' Burgers vector -→ b cs 1 = a 0 6 [121] and -→ b cs 2 = a 0 6 [211] pointing towards each other.

1 E

 1 [eV] K. Kang et al. NuMoDis

Figure 5 :

 5 Figure 5: The energy barrier as function of the constriction separation at zero stress.

Figure 6 :

 6 Figure 6: The effect of stress on the equilibrium separation between the partials. The segmented line corresponds to equation 21.

Figure 7 :

 7 Figure 7: The effect of compressive stress on the constriction energy according to DD simulations.

  et al. A. Malka-Markovitz et al. NuMoDis

Figure 8 :

 8 Figure 8: The effect of stress on the energy barrier:(A) σ g S = σ cs E = σ cs S = 0, (B) σ g S = σ g E = σ cs S = 0 and (C) σ g S = σ g E = σ cs E = 0.The green circles represent the DD simulation results, the blue diamonds correspond to the numerical solution of K. Kang's et al.[START_REF] Kang | Stress dependence of cross slip energy barrier for face-centered cubic nickel[END_REF] LT model of cross-slip (see Eq. 34) and the red crosses were obtained from the general expression for the activation enthalpy proposed by A. Malka-Markovitz et. al.[START_REF] Malka-Markovitz | Cross-slip in face centred cubic metals: a general full stressfield dependent activation energy line-tension model[END_REF] (see Eq. 36).

  et al. A. Malka-Markovitz et al. NuMoDis

Figure 9 :

 9 Figure 9: The effect of stress on the energy barrier: (A) σ g S = σ cs S = 0, (B) σ g S = σ cs E = 0 and (C) σ g S = σ g E = 0.

Preprint submitted to the Journal of Mechanics and Physics of Solids December 9, 2020

Acknowledgments

The authors acknowledge use of the 3D dislocation dynamics code NuMoDis jointly developed at CEA (L. Dupuy) and Univ. Paris XII (M. Bletry). This work has received funding from the European Research Council under the European Union's ERC grant agreement number 725142.

Conflict of Interest

Margarita Longsworth: Conceptualization, Methodology, Software, Writing -original draft. Marc Fivel: Conceptualization, Supervision, Writing -review & editing.


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