Liangzhao Huang
email: huang.liangzhao@outlook.com

Philippe Maugis
email: philippe.maugis@im2np.fr

Atomistic simulation of the collective carbon motion in body-centered tetragonal iron: A new insight into the martensite ageing

Keywords: Fe-C, Martensite, Thermodynamics, Kinetics, Elasticity, Monte Carlo

Thermodynamic and kinetic properties of carbon (C) atoms in body-centered tetragonal iron (Fe) are investigated by using Monte Carlo (MC) simulations. Pairwise interactions between carbon atoms are obtained by combining the linear elasticity theory and the state-of-the-art results from density functional theory (DFT). This energy database is applied to MC simulations to predict the equilibrium carbon configuration and the collective kinetic motion of carbon atoms in an as-quenched martensite. From the metropolis MC simulation, we obtain a novel equilibrium phase of Fe 6 C 2 structure. However, according to the kinetic MC simulations, it is difficult to reach this equilibrium phase during the martensite ageing because, at room temperature or below, the carbon diffusivity is so slow that it will take an unrealistically long time for the system to achieve the equilibrium; and at higher temperature, even though the kinetics are accelerated, the carbon concentration of the predicted equilibrium phase is so high that other metastable carbides can be formed before such equilibrium is reached. Moreover, the effects of the temperature, the applied stress, and the initial state on the ageing kinetics are highlighted. The evolution of the carbon cluster concentration and the time scale of the ageing kinetics are in good agreement with some existing experimental results.

Introduction

Martensitic steels have a wide range of mechanical applications due to their high strength. It corresponds to a carbonsupersaturated iron with a centered tetragonal (bct) structure. However, it is a metastable phase because it decomposes into carbon-free and carbon-rich nanodomains during ageing, thereby changing its mechanical properties. To fully understand the mechanisms of this decomposition, the thermodynamic and kinetic properties of this process need to be thoroughly investigated.

An in-situ X-ray diffraction study evidenced the set up of tetragonality as the martensite forms [START_REF] Epp | In situ X-Ray Diffraction Analysis of Carbon Partitioning During Quenching of Low Carbon Steel, Metall[END_REF]. It was also reported in Ref. [START_REF] Kurdjumov | Nature of axial ratio anomalies of the martensite lattice and mechanism of diffusionless γ → α transformation[END_REF] that the tetragonality of the quenched martensite increased during its ageing at room temperature. In addition, measurements of the lattice parameters of Fe-C steels with different carbon contents at room temperature have revealed a sharp transition from low-carbon cubic ferrite to high-carbon tetragonal martensite [START_REF] Nagakura | Crystallographic study of the tempering of martensitic carbon steel by electron microscopy and diffraction[END_REF][START_REF] Cheng | Lattice parameters of iron-carbon and iron-nitrogen martensites and austenites[END_REF][START_REF] Fan | Lattice-parameter variation with carbon content of martensite. II. Long-wavelength theory of the cubic-to-tetragonal transition[END_REF][START_REF] Lu | The effect of carbon content on the c/a ratio of as-quenched martensite in Fe-C alloys[END_REF]. As for the carbon distribution, multiple studies showed that the martensite ageing is accompanied by the formation of carbon-rich clusters [START_REF] Choo | Mössbauer measurements on the aging of iron-carbon martensite[END_REF][START_REF] Kusunoki | Modulated structure of iron-carbon martensite studied by electron microscopy and diffraction[END_REF][START_REF] Sherman | The Aging and Tempering of Iron-Nickel-Carbon Martensites[END_REF]. Choo and Kaplow [START_REF] Choo | Mössbauer measurements on the aging of iron-carbon martensite[END_REF] showed that the carbon-rich region in Fe-8.2C (at%) formed at room temperature after 3 months was probably an ordered one-similar to Fe 4 N. Furthermore, a new arrangement replaced the room-temperature aged structure after tempering at 82 • C, indicating the tentative formation of -carbide. Kusunoki and Nagakura [START_REF] Kusunoki | Modulated structure of iron-carbon martensite studied by electron microscopy and diffraction[END_REF] used high-resolution electron microscopic to observe the redistribution of carbon atoms in the martensite tempered from 273 to 363 K. They showed that a modulated structure was formed with small carbon clusters randomly distributed in the 102 planes. Such carbon-rich planar regions were periodically spaced with intervals of 10-20 Å. Taylor et al. [START_REF] Taylor | Spinodal decomposition during aging of Fe-Ni-C martensites[END_REF] indicated that the martensite decomposition followed the spinodal mechanism. They also showed that the carbon concentration in the carbon-rich domain of the Fe-Ni-C martensite aged at room temperature reached up to 11 at.%. In addition, a structural modulation along 203 directions was developed. Recent studies based on atom probe tomography (APT) technique confirmed that the carbon-rich phases contained 10-16 at.% C in Fe-C or Fe-Ni-C martensites tempered at low temperatures [START_REF] Zhu | Carbide characterization in lowtemperature tempered steels[END_REF][START_REF] Taylor | Spinodal decomposition during aging of Fe-Ni-C martensites[END_REF][START_REF] Nagakura | Crystallographic study of the tempering of martensitic carbon steel by electron microscopy and diffraction[END_REF][START_REF] Badinier | The Mechanisms of Transformation and Mechanical Behavior of Ferrous Martensite[END_REF][START_REF] Maugis | Carbon diffusivity and kinetics of spinodal decomposition of martensite in a model Fe-Ni-C alloy[END_REF].

On the theoretical side, the tetragonality of martensites can be understood from the assumption that carbon atoms preferentially occupy one of the three octahedral sublattices. This can be interpreted as a result of an ordering at low temperature (the so-called α phase). A simple model to describe this ordering phenomenon was proposed by Zener [START_REF] Zener | Theory of Strain Interaction of Solute Atoms[END_REF], who stated that the ordered distribution of carbon atoms is energetically more favorable than the disordered one. This idea has then been developed into a more elaborate mean-field theory of carbon ordering by Khachaturyan [START_REF] Khachaturyan | On the theory of the ordering of carbon atoms in a martensite crytal[END_REF][START_REF] Kurdjumov | Nature of axial ratio anomalies of the martensite lattice and mechanism of diffusionless γ → α transformation[END_REF]. This model has been widely used in the thermodynamic modeling of martensitic transformation [START_REF] Taylor | Aging of ferrous martensites[END_REF][START_REF] Chirkov | Tetragonality and the distribution of carbon atoms in the Fe-C martensite: Molecular-dynamics simulation[END_REF][START_REF] Maugis | Ferrite, martensite and supercritical iron: A coherent elastochemical theory of stress-induced carbon ordering in steel[END_REF][START_REF] Maugis | A Temperature-Stress Phase Diagram of Carbon-Supersaturated bcc-Iron, Exhibiting "Beyond-Zener" Ordering[END_REF] and the kinetic studies in Fe-C alloys [START_REF] Maugis | Stress-controlled carbon diffusion channeling in bct-iron: A mean-field theory[END_REF][START_REF] Huang | Effect of substitutional Ni atoms on the Snoek relaxation in ferrite and martensite Fe-C alloys: An atomisitic investigation[END_REF]. There exist many discussions on the stability and the equilibrium crystal structure of the Zener order-ing [START_REF] Sato | Octahedral defects in a b.c.c. lattice examined by lattice theory[END_REF][START_REF] Sinclair | Molecular dynamics study of the ordering of carbon in highly supersaturated α-Fe[END_REF][START_REF] Ruban | Self-trapping of carbon atoms in α-Fe during the martensitic transformation: A qualitative picture from ab initio calculations[END_REF][START_REF] Kandaskalov | A first-principle study of the structural, elastic, lattice dynamical and thermodynamic properties of α"-Fe16C2 and α"-Fe16N2 phases[END_REF][START_REF] Kandaskalov | Thermodynamic stabilities in the Fe-Fe16C2 system: Influence of carbon-carbon interactions studied by DFT[END_REF][START_REF] Maugis | Stability of Zener order in martensite: an atomistic evidence[END_REF]. The α -Fe 16 C 2 phase is one of the most known structures, similar to the α -Fe 16 N 2 phase [START_REF] Jack | The iron-nitrogen system: the preparation and the crystal structures of nitrogen-austenite( γ ) and nitrogen-martensite( α' )[END_REF]. Its stoichiometry corresponds approximately to the experimentally obtained carbon content in the carbon-rich phase. Using the lattice theory based on empirical elastic interaction between interstitial atoms, Sato [START_REF] Sato | Octahedral defects in a b.c.c. lattice examined by lattice theory[END_REF] predicted α to be the most stable ground-state structure in bcc Fe-C. This was confirmed by Sinclair et al. [START_REF] Sinclair | Molecular dynamics study of the ordering of carbon in highly supersaturated α-Fe[END_REF] using molecular dynamics simulations with empirical potentials. However, the calculations of Ruban [START_REF] Ruban | Self-trapping of carbon atoms in α-Fe during the martensitic transformation: A qualitative picture from ab initio calculations[END_REF] based on density functional theory (DFT) predicted an altered Fe 16 C 2 -structure formed by carbon atoms of two different sublattices. This contradicted the idea of Zener ordering. Recently, Maugis et al. [START_REF] Maugis | Stability of Zener order in martensite: an atomistic evidence[END_REF] have revisited the ground states of Fe-C configurations using DFT. It was shown that Zener-ordered configurations are energetically most stable. Nevertheless, the above theoretical studies focused on a homogeneous distribution of carbon atoms in space, while this is not the case in an aged martensite because carbon clusters are formed. Yan et al. [START_REF] Yan | Configurational thermodynamics of C in body-centered cubic/tetragonal Fe: A combined computational study[END_REF] performed metropolis Monte Carlo simulations to study the equilibrium state of the carbon distribution based on EAM and DFT inputs. At low temperature, they obtained a two-phase configuration with ferrite and an ordered compound (α ) different from α . Carbon atoms in α occupy two different sublattices. A Zener-type ordered structure (α ) with a single sublattice only appeared when a tetragonal strain was applied. In another study on the martensite ageing using the atomic fragment theory (AFT) [START_REF] Lavrskyi | Modélisation en fonctionnelle de la densité atomique des transformations de phases dans le système fe-c à basse température[END_REF], Fe 4 C or Fe 2 C superstructures were obtained in Zener-ordered systems depending on the magnitude of the tetragonality imposed on the matrix lattice. However, in this study, the time evolution was expressed in Fourier space and the ageing kinetics in physical time were not presented. In fact, the carbon configurations in the martensite are open to debate because none of these theoretically predicted configurations have been synthesized or observed by experiments. Moreover, martensite ageing is a time-dependent phenomenon. It is therefore unclear whether the carbon-rich phases of martensite obtained in the experiments reach equilibrium or not. Such a numerical study on the ageing kinetics is still in need.

In this article, our goal is to investigate the thermodynamic and kinetic properties of carbon related to the martensite ageing based on Monte Carlo simulations. We will combine the linear elasticity theory (ET) and the existing results from first-principles calculations based on density functional theory (DFT) to generate the energy database of pairwise interactions between carbon atoms. We will perform some parametric studies to highlight the effect of the temperature, the applied stress, and the initial state on the ageing kinetics. In addition, we will rely on our simulation results to re-investigate some experimental results from previous work.

Model

Carbon-induced strain field

An interstitial carbon solute can be on the octahedral or tetrahedral sites of the iron superstructure (cf. Fig. 1). Within the elastic theory, this carbon solute can be described by the elastic dipole P [START_REF] Varvenne | Point defect modeling in materials: Coupling ab initio and elasticity approaches[END_REF][START_REF] Varvenne | Elastic dipoles of point defects from atomistic simulations[END_REF][START_REF] Clouet | Elastic modeling of point-defects and their interaction[END_REF]. It is a second-rank tensor, defined as the first moment of the point-force distribution induced by the carbon atom. The corresponding displacement field is thus given by

u i (r -r 0 ) = - n,m,p ∂G i j ∂r k (r -r 0 -r nmp )P jk , (1)
where G i j is the elastic Green's function, r 0 is the solute position, and r nmp = nR 1 + mR 2 + pR 3 is the positions of the periodic images, with n, m, p ∈ Z and R 1 , R 2 , and R 3 the periodicity vectors. It can be obtained by computationally-intensive atomistic calculation. In a particular case, for example within an infinite medium, the elastic Green's function, G ∞ i j , can be efficiently calculated using the Barnett's numerical scheme [START_REF] Barnett | The precise evaluation of derivatives of the anisotropic elastic Green's functions[END_REF]. In order to avoid the complex atomistic calculation, we rely on G ∞ i j to compute the elastic field in a finite volume. For a volume V without applied stress, the strain field induced by an interstitial solute can be written as [START_REF] Balluffi | Introduction to Elasticity Theory for Crystal Defects[END_REF] = ∞ + cor .

(

) 2
In this equation, ∞ corresponds to the displacement field in a volume V surrounded by an infinite medium. Thus, it is given by

∞ i j (r -r 0 ) = 1 2       ∂u ∞ i ∂r j (r -r 0 ) + ∂u ∞ j ∂r i (r -r 0 )       , (3)
= - n,m,p ∂ 2 G ∞ ik ∂r l ∂r j (r -r 0 -r nmp )P kl , (4)
cor is a term to correct ∞ by applying a virtual stress to make the volume V traction-free. We assume the virtual stress to be homogeneous in space; thereby the corresponding strain field is also homogeneous: cor = cor . One can demonstrate that the homogeneous strain field of the traction-free volume containing a solute is given by [START_REF] Balluffi | Introduction to Elasticity Theory for Crystal Defects[END_REF]

i j = 1 V S (V) u i (r -r 0 ) n j dS = S i jkl P kl V , (5)
where S is the elastic compliance tensor, S (V) is the external surface of the volume V, and n is the normal vector to the surface. Therefore, we have cor

i j = i j -∞ i j (6) = 1 V         S i jkl + S (V) n,m,p ∂G ∞ ik ∂r l (r -r 0 -r nmp ) n j dS         P kl . (7)
Note that i j is independent of the carbon position in the volume. Moreover, we check that ∞ i j is nearly independent of the carbon position if the volume V is large enough. For the sake of simplicity, we assume in Eq. ( 7) that the carbon atom is located at the origin. In case there are multiple interstitial carbon atoms in the volume, the total strain field is the linear superposition of the individual strain field induced by each carbon.

Carbon-carbon interactions

Following the elastic theory (ET) presented in the previous section, the interaction between two carbon atoms [P (1) at r 1 and P (2) at r 2 ] corresponds to the interaction between one atom and the strain field generated by the other, which is given by V ET C-C = -P (1) i j i j (r 1 -r 2 ) = P (1) i j P (2) kl g i jkl (r 1 -r 2 ), [START_REF] Kusunoki | Modulated structure of iron-carbon martensite studied by electron microscopy and diffraction[END_REF] where g i jkl (r 1 -r 2 ) can be deduced from Eqs. (2), ( 4) and [START_REF] Choo | Mössbauer measurements on the aging of iron-carbon martensite[END_REF]:

g i jkl (r 1 -r 2 ) = n,m,p       ∂ 2 G ∞ ik ∂r l ∂r j (r 1 -r 2 -r nmp ) - 1 V S (V) ∂G ∞ ik ∂r l (r 1 -r 2 -r nmp ) n j dS - S i jkl V . (9)
Note that this approach captures the long-range elastic interactions between two carbons, while it becomes incorrect when the carbons are close (within a few lattice parameters) because the nonlinear elasticity, lattice discreteness and chemical bonding are no longer negligible in the latter case. In general, these short-range interactions can be obtained from ab initio calculation (DFT) or atomistic simulations based on EAM potentials. Therefore, when the distance between two carbon atoms is shorter than a cut-off radius r c , we use total interactions from a recent ab initio database (V DFT C-C ) [START_REF] Kandaskalov | Carbon diffusion in bccand bct-Fe: Influence of short-range C-C pair interactions studied from first-principles calculations[END_REF], otherwise from Eq. ( 8):

V C-C =        V DFT C-C , |r 1 -r 2 | ≤ r c V ET C-C , |r 1 -r 2 | > r c . ( 10
)
Note that we consider not only the interaction between two carbon atoms on the octahedral sites [i.e. V C-C (OO)], but also the interaction between a carbon atom on the tetrahedral site and the other one on the octahedral site [i.e. V C-C (TO)]. For V C-C (OO), the cut-off radius r c (OO) is set to 1.23 a 0 , and for V C-C (TO), the cut-off radius r c (TO) is set to 1.60 a 0 . We will justify these values in Section 3.1. As we assume that in the simulation box there is no more than one carbon atom on the tetrahedral site, the interaction between two carbon atoms on the tetrahedral site will not be taken into account.

System enthalpy

From the linear elasticity theory, the elastic energy of a crystal lattice of volume V containing N carbon atoms is

E = 1 2 V C : + 1 2 N i=1 N j=1 V C i -C j , ( 11
)
where C is the elastic tensor and is the tensor of the homogeneous strain resulting from all the carbon interstitials. The first term is the lattice strain energy, and the second term is an Ising-type Hamiltonian summing over all pair interaction energies between two carbon interstitials (denoted by C i and C j ). Under an applied homogeneous stress σ ext , the system enthalpy is derived from the energy by a Legendre transformation:

H = 1 2 V C : + 1 2 N i=1 N j=1 V C i -C j -V σ ext : . (12)
Note that the reference state for E and H is the undistorted iron lattice. The first and third terms of H are related to the homogeneous elastic field, and they describe the lattice response to the homogeneous strain/stress. However, only the second term of H depends on the local elastic field of the system, and it describes the lattice response to the spatial distribution of the carbon atoms. This term is central to our study since it is related to the carbon clustering and the equilibrium cluster structure.

Elasto-diffusion

The carbon migration enthalpy depends on the local chemical environment and the applied stress. The pairwise interaction model presented in the previous section is applied to describe the variation of the carbon migration barriers with the local environments. These pairwise interactions are included in the frequency of a carbon jump from an octahedral site o through a transition site t (o/t):

ω o/t = ν 0 exp        - H m o/t + (E TO C-C -E OO C-C ) -(P T i j -P O i j ) S i jkl σ ext kl k B T        , (13)
where H m o/t is the stress-free carbon migration enthalpy from o through t, k B is the Boltzmann constant, T is temperature, ν 0 is the attempt frequency that is assumed to be strain-and composition-independent, P T and P O are dipole tensors of carbon on tetragonal and octahedral sites (cf. Fig. 1), σ ext is an externally applied homogeneous stress, and

E TO C-C = k z k V k C-C (TO), E OO C-C = k z k V k C-C (OO), (14)
where z k is the number of carbon atoms at the k-th nearestneighbor shell of the considered carbon atom, and

V k C-C (TO) [resp. V k C-C (OO)]
is the pairwise interaction between a carbon atom on an octahedral site and the other one on a tetrahedral (resp. octahedral) site at the k-th nearest-neighbor distance.

Results

Pairwise interaction energies and cut-off distance

Using the elastic dipoles and elastic constants from Eyméoud et al. [START_REF] Eyméoud | Impact of Ni alloying on Fe-C martensite ageing: an atomistic investigation[END_REF], we obtain the C-C pair interactions based on the ET model presented in Section 3.1. These interactions are plotted in Fig. 2 as a function of the C-C distance. In the same figure, we also plot the C-C pair interactions obtained from two other approaches: DFT calculation by Kandaskalov et al. [START_REF] Kandaskalov | Carbon diffusion in bccand bct-Fe: Influence of short-range C-C pair interactions studied from first-principles calculations[END_REF]; atomistic simulations from Zhang et al. [START_REF] Zhang | Mechanism of collective interstitial ordering in Fe-C alloys[END_REF] based on Lau's EAM potentials [START_REF] Lau | Many-body potential for point defect clusters in Fe-C alloys[END_REF].

The pairwise interactions between two carbon atoms on the octahedral site given by both DFT and EAM indicate a strong C-C repulsion (positive value) within one bcc lattice parameter a 0 [cf. Fig. 2 (top)]]. As expected, within 1.23 a 0 , C-C interactions given by ET are very different from those by DFT and EAM. Beyond 1.23 a 0 , both repulsive and attractive interactions exist, and the interactions are much weaker. DFT and EAM give different signs of some pair interactions. Moreover, DFT gives more repulsive interactions than EAM. This difference is due to an effect of the periodic images because the size of the supercell used in DFT is in general small. This effect has also been highlighted in Ref. [START_REF] Yan | Configurational thermodynamics of C in body-centered cubic/tetragonal Fe: A combined computational study[END_REF]. However, the signs of the pair interactions given by ET are in good agreement with those from EAM potentials. Moreover, beyond 1.23 a 0 , the deviation between the results by ET and EAM is satisfactory. Therefore, the pair interactions from ET beyond 1.23 a 0 is of good reliability. This corresponds to the chosen value of r c (OO).

The pairwise interactions between a carbon atom on the octahedral site and the other one on the tetrahedral site are also plotted in Fig. 2 (bottom). However, to the authors' knowledge, there is no existing database from EAM potentials. We compare the interactions from ET to those given by first principles [START_REF] Kandaskalov | Carbon diffusion in bccand bct-Fe: Influence of short-range C-C pair interactions studied from first-principles calculations[END_REF]. Even though the ET approach reproduces relatively well the global variation tendency of pair interactions with the C-C distance, the quantitative agreement with the ab-initio database is relatively poor at small distance (< 1.60 a 0 ). However, beyond 1.60 a 0 , the agreement between the two approaches is satisfactory. Therefore, r c (TO) is set to 1.60 a 0 .

Equilibrium carbon configuration

We carry out on-lattice Metropolis Monte Carlo simulations under the NPT-ensemble to compute the equilibrium carbon configuration using the C-C pair interactions given by Eq. [START_REF] Taylor | Spinodal decomposition during aging of Fe-Ni-C martensites[END_REF]. The system Hamiltonian is given by Eq. [START_REF] Badinier | The Mechanisms of Transformation and Mechanical Behavior of Ferrous Martensite[END_REF].

Here, we consider the bcc Fe-2C (at.%) system at room temperature (300 K). In Fig. 3, we plot the evolution of the selected short-range order (SRO) parameters (b), as well as the system enthalpy variation (c), during the relaxation. We also present the equilibrium carbon configuration obtained at the end of the relaxation (Fig. 3-d). Atomic configurations are visualized using OVITO [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITO-the open visualization tool[END_REF].

The system enthalpy decreases with increasing simulation steps, and at the end it oscillates around a relatively low level (∼ -6.8 meV/Fe), meaning that the system is very close to the equilibrium state. For a chosen carbon atom, the SRO parameter corresponding to the i-th nearest neighbor (NN) shell is defined as:

SRO(i) = N i -N random N max -N random , (15)
where N i is the number of the occupied i-th NN sites, N random is the theoretical number of the occupied i-th NN sites by assuming that the carbon atoms are randomly distributed, and N max is the possibly maximum number of the occupied i-th NN sites (i.e., the coordination number). The result presented in Fig. 3b is the average SRO parameter of all carbon atoms. The six highest SRO parameters at the end of the simulation are the configurations '4A', '10A', '7A', '3A', '7B', and '11A', presented in Fig. 3-a. These configurations correspond to pairs of carbon atoms on the same sublattice. This is in line with the final carbon distribution, where all the carbon atoms occupy the same sublattice (z).

The final carbon distribution is a plate-form cluster. This cluster consists of a plane of normal vector (103) and another plane of normal vector (103). The 'zigzag' pattern of this cluster may be an artificial result of the periodic boundary. The local carbon content (C/Fe) in this cluster is 25-30 at.%. This concentration indicates a Fe 4 C or Fe 3 C stoichiometry. According to the final SRO parameters, we deduce that the crystal structure of this cluster has a Fe 6 C 2 structure as showed in Fig. 4. This structure is composed of the '4A' and '10A' configurations, corresponding to the two highest SRO parameters. It is different from the known structures like the α -Fe 16 C 2 phase. We will explain this difference relying on the C-C interaction database in Section 4.1.

Kinetic evolution of martensite ageing

We rely on the on-lattice kinetic Monte Carlo approach to simulate the evolution of martensite ageing. The residence-time algorithm is applied. Periodic boundary conditions are applied to the bcc 25 × 25 × 25 a 3 0 -supercell, with inserted carbon atoms on the octahedral sites. The number of carbon atoms is adjusted by the carbon atomic fraction. The carbon jump frequency is obtained from the pair interaction model, as presented in Section 2.4.

In this section, we consider the bcc Fe-2C (at.%) system at room temperature (300 K), starting from an initial state where all carbon atoms are on z-sublattices and they are homogeneously distributed in the simulation box. Following the idea of Bain's transformation [START_REF] Bain | The nature of martensite, trans[END_REF][START_REF] Roitburd | The nature of martensitic transformations[END_REF], this initial state represents a simulation of the bct structure just after quenching the austenitic structure from a high temperature, during which all the carbon atoms 'migrate' to one preferential sublattice (e.g., the z-sublattice) through an elastic diffusionless mechanism due to the lattice distortion along z-direction. Note that this is not the unique way to represent the as-quenched state, and we will investigate, in Section 3.6 the effect of the initial state on the evolution of martensite ageing.

In Fig. 5, we plot the evolution of the selected SRO parameters, the system enthalpy, the sublattice occupancy, and the carbon distribution in Fe-2C at 300 K. The selected SRO param-Figure 2: Top: Interactions between two carbon atoms of the same variant on the octahedral site (left) and those between two carbon atoms of different variants on the octahedral site (right). Bottom: Interactions between a carbon atom on the octahedral site and the other one on the tetrahedral site. These interaction energies are obtained from different approaches: DFT calculation [START_REF] Kandaskalov | Carbon diffusion in bccand bct-Fe: Influence of short-range C-C pair interactions studied from first-principles calculations[END_REF], atomistic simulations using EAM potentials [START_REF] Zhang | Mechanism of collective interstitial ordering in Fe-C alloys[END_REF], and elastic theory (ET) based on force dipoles. Negative energies correspond to attractive interactions. The nomenclature of the interaction shell is the same as the one in Ref. [START_REF] Kandaskalov | Carbon diffusion in bccand bct-Fe: Influence of short-range C-C pair interactions studied from first-principles calculations[END_REF]. eters correspond to the probability of the most attractive OOconfigurations presented in Fig. 5-a. These configurations consist of two carbon atoms on the same sublattice. The selected SRO fast increases with time. This also leads to a rapid decrease of the pair-interaction part of the system enthalpy, which mainly contributes to the reduction of the total enthalpy (cf. Fig. 5-b). This is because, in the initial state, the carbon atoms are ran-domly placed in the bcc lattice, and thus there are many repulsive C-C pair configurations. The number of these repulsive configurations decreases rapidly after several hundred jumps per carbon atom, resulting in the sharp decrease of the system enthalpy. It is worthy noting that at this early stage of the evolution (within 1 hour), a disordering is observed: the z-sublattice occupancy decreases to 0.85 and the occupancy of the other two sublattices increases up to 0.10 (cf. Fig. 5-c). However, this disordering does not last for long before the carbon atoms are reordered. After 10 hours, almost all the carbon atoms occupy the z-sublattice. This corresponds well to the variation of the tetragonality along the z-direction (cf. Fig. 5-d). The tetragonality is initially high because all the carbon atoms are set to be at the z-sublattice. Then, it decreases during the rapid disordering and increases up to the initial value as the carbon atoms are reordered. Moreover, from the evolution of the carbon distribution (Fig. 5-e), carbon clustering occurs after 100 hours. In the final state, we observe that several carbon-enriched clusters are separated by carbon-free space. The variations of the SRO parameters and the system enthalpy are not stabilized in the final state of the simulation. In addition, the enthalpy of the final system (-3.4 meV/Fe) is remarkably higher than the equilibrium value (-6.8 meV/Fe from Fig. 3). Therefore, the system does not reach equilibrium after a physical time of 2300 hours (about 3 months). This also indicates that the time required to reach equilibrium at room temperature can be very long (more than 3 months), because the carbon diffusion is slow at relatively low temperature and, in addition, carbon atoms can be kinetically trapped by the martensite variant which is reflected by a very small kinetic correlation factor (smaller than 0.01).

In order to provide a more detailed investigation of the carbon clusters, we show in Fig. 6-a the distribution of the average carbon concentration on different selected planes. The concentration distributions on the (001), (010) and (100) planes are obtained by averaging the carbon concentration along the [001],

[010] and [100] directions, respectively. On the (001) plane, the carbon distribution is relatively homogeneous, while on the (010) and (100) planes the concentration field is very heterogeneous. This means that the cluster separation follows the direction z (i.e., [001]). We are also interested in how high the local carbon concentration can be in the clusters. We select two representative carbon clusters. Their positions are highlighted in the figure (i.e., the two control volumes containing Clusters 1 and 2 indicated in Fig. 6-b). The carbon concentrations in these two volumes increase with time and they are 8-10 at.% at the end of simulation (cf. Fig. 6-c). This evolution supports the spinodal nature of this decomposition. Moreover, by comparing the carbon distribution at 10 and 2300 hours, we observe that the spacing between carbon clusters increases with time. It is also worthy noting that there is not a single orientation for the habit planes of the carbon clusters, which consist of the facet family {102}. For example, the orientation of the habit plane of Cluster 1 is (102) and the one of Cluster 2 is (0 12).

Effect of applied stress

Since residual stresses exist in martensitic microstructure after quench, it is important to investigate the effect of the applied stress on the evolution of the martensite ageing. Two representative types of stress are discussed below.

In Fig. 7, we compare the results of the SRO evolution and the carbon distribution obtained with and without applying an uniaxial tension. This corresponds to a tension of 800 MPa along direction z. The value of this tension is within the range of the experimental ones [START_REF] Macherauch | Residual Stresses After Quenching[END_REF][START_REF] Weise | Martensitic Transformation and Residual Stresses after Thermomechanical Treatment of Heat Treatable Steel 42CrMo4 (SAE 4140)[END_REF]. Such an applied tension slows down the SRO evolution: it takes about 5 times longer for the stressed system to reach the same SRO state as the unstressed system. Moreover, by comparing the carbon distribution at the same physical time, the carbon decomposition is more remarkable in the unstressed lattice. The slower kinetic evolution should be due to the reduced carbon diffusivity under applied stress. We compute the tracer diffusivity of carbon in these cases by measuring the mean squared displacement of carbon atoms by Monte Carlo simulations. The average diffusivity (D) in the unstressed system is 1.2(±0.4) × 10 -25 m 2 /s and the one in the stressed system is 2.9(±0.7) × 10 -26 m 2 /s. Therefore, the carbon diffusion in the stressed system is about 4 times slower than that in the unstressed system. This corresponds well with the deceleration rate of the SRO evolution.

In Fig. 8-a, we investigate the early stage of the carbon evolution in Fe-2C submitted to a pure shear stress such that σ 22 = -σ 33 = 200 MPa. The applied compression along z direction (σ 33 = -200 MPa) disfavors carbon atoms on the zsublattice, while the applied tension along y direction favors carbon atoms on the y-sublattice. The resulting variation of the sublattice occupancy is rapid (cf. Fig. 8-b). Within 100 s, the carbon atoms, which are all initially on the z-sublattice, are located mainly on the y-sublattice. For a longer simulation time, both short-and long-range orderings are expected. They are similar to those presented in Fig. 5-e, except that the decomposition of the carbon clusters is expected to be along the y direction in this case.

In brief, the variation of the spatial distribution and the sublattice occupancy of the carbon atoms with the applied stresses suggests that one can modify the orientation of the variant by applying external stresses.

Effect of temperature

In this section, we study the temperature effect on the kinetics of the carbon evolution in Fe-2C. In Fig. 9 (top), we compare the SRO evolution at different temperatures. Since the carbon kinetics increases with temperature, martensite ageing is faster at higher temperatures. For instance, at 400 K, the SRO rapidly increases: the SRO corresponding to the configuration 10A is up to 0.160 after a physical time of 30 minutes. However, at 250 K, the SRO variation is much slower: the SRO for 10A configuration is about 0.125 after 100 years. In Fig. 9 (bottom), we compare also the final states obtained at different temperatures but at the same Monte Carlo step. At 400 K, a notable carbon decomposition is obtained after a physical time of 30 minutes, while at 250 K, only small carbon clusters are obtained after a physical time of 150 years. The remarkable difference among the carbon clusters obtained in the final state of the simulation at different temperatures is the SRO value for the 4A configuration. It corresponds to the second nearest distance between two carbon atoms of the same sublattice. Therefore, the higher this value, the closer to the equilibrium Fe 6 C 2 phase, the more carbon-rich the clusters are. The substantial increase of the 4A-SRO value can be seen as a sign of the long-range ordering.

In Fig. 10, we provide a quantitative study on the concentration of the carbon clusters obtained in the final state. The planar distribution of the carbon concentration shows that the distance between the carbon-enriched area at 400 K is significantly larger than that at 300 K (Fig. 6). We choose to consider two representative clusters at the positions respectively highlighted by the two control volumes (i.e., Clusters 3 and 4 indicated in Fig. 10-b). The carbon concentration in these two zones increases up to 14 at.% in the final state (cf. Fig. 10-c). Moreover, the habit planes of Clusters 3 and 4 are also in the facet family {102}, which is the same as Clusters 1 and 2 in the case of 300 K.

Note that we have performed another metropolis Monte Carlo simulation of Fe-2C at 400 K (not presented). The obtained equilibrium phase was almost the same as the one from 300 K (Fig. 3-d). This supports the fact that increasing the temperature from 300 K to 400 K accelerates the martensite ageing but hardly changes the equilibrium state of the carbon cluster.

Effect of the initial state

In this section, we investigate the effect of the initial state on the evolution of the carbon distribution. In Fig. 11, we compare the decomposition kinetics at 400 K obtained from two different initial states: one starts from a state with carbon atoms randomly distributed on the z sublattice (cf. Fig. 11-a); the other starts from a state with carbon atoms randomly distributed on the three sublattices such that the occupancy of each sublattice is equal (cf. Fig. 11-b). The initial state of the former case corresponds to the one we used in the simulations presented in previous sections. In this case, the evolution of the carbon distribution at 400 K is similar to that obtained at 300 K (cf. Fig. 5): a rapid partial disordering is followed by short-range ordering and then long-range decomposition. The occupancy of z sublattice is almost always 1 except during the short disorder state at the very beginning of the evolution. However, starting from the initial state with equal occupancy of each sublattice, clusters of carbon atoms on different sublattices coexist at the early stage (≤ 100 s). We choose to call the cluster of carbon atoms on i sublattice (i = x, y, z) the i-variant. Before 100 s, x-and y-variants are bigger than the z-variant. In addition, the size of z-variant decreases with time and disappears after several hundred seconds. Then, the size of x-variant increases with time while the one of y-variant decreases. This evolution is also reflected on the variation of the sublattice occupancy with time [cf. Fig. 12 (left)]. The variation of the occupancy is not stabilized at the end of simulation (∼ 2000 s), but we can predict from the obtained evolution that the occupancy of the x sublattice continues to increase while the one of the y sublattice continues to decrease. Therefore, only one variant will be left after a longer time. However, the coexistence of multiple variants may not be stable in the presence of applied stress because a tension along a chosen direction favors the corresponding variant. In order to verify this point, we performed another simulation starting also from an initial state with equal sublattice occupancy. In addition, an uniaxial tension of 100 MPa along z direction is applied to the simulation box. In this case, the occupancy of z sublattice increases rapidly [cf. Fig. 12 (right)] and only z-variant is found after 100 s.

Discussion

Comparison with other simulation studies

Zener-type carbon ordering has been obtained from our simulations: the carbon cluster is of single variant. This is in line with the C-C pair interactions [cf. Fig. 2 (top)]. The six most attractive (negative) interactions between two carbon atoms are for the 10A (-0.09 eV), 7A (-0.07 eV), 11A (-0.06 eV), 9A/9B (-0.05 eV), 12A/12B (-0.04 eV), and 8A (-0.01 eV) configurations. The first three of them (10A, 7A, 11A) correspond to pairs of carbon atoms of the same variant (aa-type), while the other three (9A/9B, 12A/12B, 8A) are for carbon atoms of different variants (ac-type). In addition, the aa-type configurations are more attractive (negative) than the ac-type configurations. Moreover, the total coordination number of the three aa-type configurations is close to the one of the three ac-type configurations. According to theses results, carbon atoms prefer to be at only one octahedral sublattice. However, according to the C-C interaction database used in Refs. [START_REF] Ruban | Self-trapping of carbon atoms in α-Fe during the martensitic transformation: A qualitative picture from ab initio calculations[END_REF] and [START_REF] Yan | Configurational thermodynamics of C in body-centered cubic/tetragonal Fe: A combined computational study[END_REF], acconfigurations are more attractive than aa-configurations. This explains why carbon atoms in the ordered structures obtained in these studies (e.g., the so-called α phase) occupy two different sublattices. Therefore, the Zener-type ordering tendency and the ordered structure depend much on the chosen C-C interaction database. A possible way to justify the validity of the database is to compare the theoretically obtained carboninduced tetragonality with the experimental one. The tetragonality of a fully ordered martensite can be defined as

t = c a 0 = 1 + γ u C , ( 16
)
where c is the lattice parameter along the tetragonal direction in the fully ordered martensite, a 0 is the carbon-free lattice param-eter, and γ describes the linear dependence of t of the fully ordered phase on the atomic fraction of carbon atoms, u C (C/Fe). In Ref. [START_REF] Yan | Configurational thermodynamics of C in body-centered cubic/tetragonal Fe: A combined computational study[END_REF], from atomistic simulations based on EAM potentials of Veiga et al. [START_REF] Veiga | Comments on "Atomistic modeling of an Fe system with a small concentration of C[END_REF], γ was found to be 0.38. However, from the DFT-based database [START_REF] Eyméoud | Impact of Ni alloying on Fe-C martensite ageing: an atomistic investigation[END_REF] chosen in our work, γ = 0.93, which is close to the experimental values: γ = 0.85 from Cheng et al. [START_REF] Cheng | Lattice parameters of iron-carbon and iron-nitrogen martensites and austenites[END_REF] and γ = 1.03 from Rayne and Chandrasekhar [START_REF] Rayne | Elastic constants of iron from 4.2 to 300[END_REF]. To this extent, the DFT-based C-C pair interactions chosen in the present study is of higher reliability than those in Ref. [START_REF] Yan | Configurational thermodynamics of C in body-centered cubic/tetragonal Fe: A combined computational study[END_REF].

Note that the crystal structure of the carbon cluster obtained in the present work is different from the well known α phase, which is composed of 11A configurations. This difference is due to the choice of the energy database because the simulated equilibrium phase is sensitive to the C-C interactions. According to the pair interactions chosen in this work, the 10A configuration (-0.09 eV) is more stable (negative) than the 11A configuration (-0.06 eV) even though both of them are negative. A recent DFT calculation from Kandaskalov et al. [START_REF] Kandaskalov | Carbon diffusion in bccand bct-Fe: Influence of short-range C-C pair interactions studied from first-principles calculations[END_REF] also showed that the energy of 10A configuration (0.015 eV) is lower than that of 11A (0.030 eV). Moreover, another recent C-C interaction database from Zhang et al. [START_REF] Zhang | Mechanism of collective interstitial ordering in Fe-C alloys[END_REF] based on EAM potential and adjusted by DFT calculation showed that 10A and 11A configurations have almost the same energy (∼ -0.04 eV). These results support the existence of the Fe 6 C 2 -phase (composed of the 10A configurations) obtained in the present study.

Revisiting the experimental results

Our kinetic simulations of the martensite ageing allows for a new interpretation of some previously published experimental results.

Our simulations show that carbon-rich clusters are formed in Fe-2C (at.%), with a local carbon concentration up to 10 at.% after a physical time of 3 months at 300 K. The predicted carbon concentration and the ageing rate correspond well to the experimental results of Kusunoki and Nagakura [START_REF] Kusunoki | Modulated structure of iron-carbon martensite studied by electron microscopy and diffraction[END_REF], Taylor et al. [START_REF] Taylor | Spinodal decomposition during aging of Fe-Ni-C martensites[END_REF], Zhu et al. [START_REF] Zhu | Carbide characterization in lowtemperature tempered steels[END_REF], and Maugis et al. [START_REF] Maugis | Carbon diffusivity and kinetics of spinodal decomposition of martensite in a model Fe-Ni-C alloy[END_REF]. The obtained cluster was believed to be of the α -Fe 16 C 2 phase because its stoichiometry corresponds rather well to the measured carbon concentration. However, the diffraction experiments conducted by Van Genderen et al. [START_REF] Van Genderen | Early Stages of Decomposition in Iron-Carbon and Iron-Nitrogen Martensites: Diffraction Analysis Using Synchrotron Radiation[END_REF] showed no evidence of diffraction peaks due to the α structure. According to our simulation results, the obtained carbon cluster should be a transition phase of the Fe 6 C 2 structure, as presented in Fig. 4. The equilibrium cannot be reached within a realistic experimental time at room temperature due to the low carbon diffusivity and the remarkable kinetic trapping of carbon atoms by the carbon cluster. This also explains why the carbon concentrations obtained in similar room-temperature ageing experiments vary from 10 to 15 at.% but do not exactly corresponds to the Fe 16 C 2 stoichiometry.

Increasing the temperature substantially accelerates the ageing kinetics. An evident carbon decomposition can be observed experimentally in Fe-15Ni-1C (at.%) within one hour at higher temperatures (95-150 • C) [START_REF] Taylor | Spinodal decomposition during aging of Fe-Ni-C martensites[END_REF]. Moreover, within the same ageing time, the higher the temperature, the higher the carbon concentration of the resulting carbon cluster [START_REF] Taylor | Spinodal decomposition during aging of Fe-Ni-C martensites[END_REF]. According to our simulation results, this does not mainly result from a thermodynamic reason. We have checked that the simulated equilibrium structures obtained at 300 and 400 K are nearly identical by metropolis Monte Carlo simulations. However, we showed that the carbon concentration of the cluster obtained at 400 K (123 • C) at a physical time of 30 minutes is higher than that obtained at room temperature (300 K) at a physical time of 3 months. This is because the former cluster is closer to the equilibrium structure (Fe 6 C 2 ) than the latter due to a higher carbon diffusivity. It is worth noting that the increase of the carbon concentration of the cluster obtained at higher temperature can be a precursor of transition carbide formation (e.g., the -carbide). This point requires further investigations.

Taylor et al. [START_REF] Taylor | Spinodal decomposition during aging of Fe-Ni-C martensites[END_REF] have also shown that the ageing kinetics in Fe-Ni-C alloys is slowed down by a higher carbon concentration. This can be understood from two aspects. On one hand, the increase of the carbon concentration leads to an higher carbon-induced tension. According to the simulation results presented in Section 3.4, an uniaxial tension slows down the martensite ageing due to a lower carbon diffusivity. On the other hand, the higher the carbon concentration, the more significant the kinetic correlation resulting from the "labyrinth" effect due to carbon-carbon interactions [START_REF] Kandaskalov | Carbon diffusion in bccand bct-Fe: Influence of short-range C-C pair interactions studied from first-principles calculations[END_REF].

Conclusion

In this article, we have performed Monte Carlo (MC) simulations to investigate the thermodynamic and kinetic properties of carbon atoms in body-centered tetragonal iron using an energy database obtained by combining the linear elasticity theory (ET) and the state-of-the-art results from density functional theory (DFT). The most relevant results are summarized as below.

• According to our database of the C-C interactions, the interaction between two carbon atoms on the same octahedral sublattice is, in general, more attractive than that between two carbon atoms on different sublattices. This suggests a Zener-type ordering of the clusters during the martensite ageing, which has also been confirmed by metropolis MC simulations.

• The equilibrium carbon configuration obtained from metropolis MC simulations is of the Fe 6 C 2 structure, which is different from the well known α -Fe 16 C 2 phase. This is because the so-called '10A' configuration (comprising the Fe 6 C 2 structure) is more attractive than the '11A' configuration (comprising the α phase).

• The kinetic MC simulations show that the time to reach equilibrium state during the martensite ageing at room temperature can be so long that the experimentally ob-served carbon clusters are in transition states. Moreover, in practice, the theoretically predicted equilibrium phase may not be reached because the carbon concentration is so high in this phase that carbide with structures other than bcc may be formed due to severe lattice distortion. However, in the present study, this phase is of bcc structure by construction. This can lead to the difference with those calculated from other crystal structures (e.g., hcp, fcc) and with those observed by experiments with a much longer time or at a higher temperature.

• By applying an uniaxial tension, the ageing kinetics is significantly slowed down due to a smaller carbon diffusivity. These points correspond well to the experimental observations. Moreover, this uniaxial tension favors the formation of carbon variants of the same direction as the applied tension.

• Compared to the room-temperature ageing, increasing the temperature accelerates the carbon decomposition and leads to the formation of carbon clusters of higher concentration within a much shorter time.

Figure 1 :

 1 Figure 1: Body-centered orthorhombic cell with iron atoms (grey), octahedral interstitial sites (large colored spheres, O x,y,z ), and tetrahedral interstitial sites (small colored spheres, T x,y,z ). Red (R), green (G), blue (B) colors refer to site types according to R = x, G = y, B = z.

Figure 3 :

 3 Figure 3: a: Key OO-configurations. The violet sphere indicates the first (static) C atom and blue spheres represent the possible positions of the second C atom to form different C-C configurations. b: Evolution of the short-range order (SRO) parameters. c: Evolution of the system enthalpy. d: The final carbon distribution obtained after 2 × 10 1 0 Monte Carlo steps (MCS). Red (R), green (G), blue (B) colors refer to site types according to R = x, G = y, B = z. The above simulation results are obtained in Fe-2C (at.%) at 300 K in a 25 × 25 × 25 a 3 0 supercell.

Figure 4 :

 4 Figure 4: Predicted crystal structure of the equilibrium carbon configuration.

Figure 5 :

 5 Figure 5: a: Evolution of the short-range order (SRO) parameters. b: Evolution of the system enthalpy. The pair-interaction part of the system enthalpy is the sum over all C-C pair interaction energies. c: Evolution of the sublattice occupancy. d: Evolution of the tetragonality. ∆c is the variation of the lattice parameter along the tetragonal direction. e: Evolution of the carbon distribution. Red (R), green (G), blue (B) colors refer to site types according to R = x, G = y, B = z. The above simulation results are obtained in Fe-2C (at.%) at 300 K in a 25 × 25 × 25 a 3 0 supercell.

Figure 6 :

 6 Figure 6: a: Planar distributions of the average carbon concentration in different planes. b: Positions of the two representative clusters. Large spheres (orange and olive) represent carbon atoms within the control volumes of clusters. Small spheres (blue) represent carbon atoms outside these volumes. c: Evolution of the average carbon concentration in the two indicated clusters. These simulation results are obtained in Fe-2C (at.%) at 300 K in a 25 × 25 × 25 a 3 0 supercell after 4 × 10 7 simulation steps.

Figure 7 :

 7 Figure 7: Top: Comparison between the short-range order evolution of Fe-2C with and without applied stress. Bottom: The C distributions in Fe-2C with and without applied stress at the same physical time (2300 hours). The simulation temperature is set to 300 K.

Figure 8 :

 8 Figure 8: Evolution of the carbon distribution (a) and sublattice occupancy (b) in Fe-2C submitted to a pure shear stress such that σ 22 = -σ 33 = 200 MPa. Red (R), green (G), blue (B) colors refer to site types according to R = x, G = y, B = z. The simulation temperature is set to 300 K.

Figure 9 :

 9 Figure 9: Top: Evolution of short-range order (SRO) parameters in Fe-2C at different temperatures. Bottom: The final carbon distribution obtained in Fe-2C after 40 000 000 Monte Carlo steps (MCS) at different temperatures.

Figure 10 :

 10 Figure 10: a: Planar distributions of the average carbon concentration in different planes. b: Positions of the two representative clusters. Large spheres (orange and olive) represent carbon atoms within the control volumes of clusters. Small spheres (blue) represent carbon atoms outside these volumes. c: Evolution of the average carbon concentration in the two indicated clusters. These simulation results are obtained in Fe-2C (at.%) at 400 K in a 25 × 25 × 25 a 3 0

Figure 11 :

 11 Figure 11: Evolution of the carbon distribution obtained with a random single-sublattice initial state (a), a random completely-disordered initial state (b), and a random completely-disordered initial state with an applied uniaxial tension of 100 MPa along z direction (c). Red (R), green (G), blue (B) colors refer to site types according to R = x, G = y, B = z. These simulation results are obtained in Fe-2C (at.%) at 400 K in a 25 × 25 × 25 a 3 0 supercell.

Figure 12 :

 12 Figure 12: Evolution of the occupancy resulting from a random completelydisordered initial state without applied stress (left) and with an uniaxial tension of 100 MPa along z direction (right). The simulation temperature is set to 400 K.

Acknowledgments

This work was supported by the Agence Nationale de la Recherche (contract C-TRAM ANR-18-CE92-0021). Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources.

Data availability

The authors confirms that the data supporting the findings of this study are available within the article. Raw data that support the findings of this study are available from the corresponding author, upon reasonable request.