Ferritic steels, made of a body-centered cubic (bcc) iron matrix with interstitial carbon solutes, are widely-used structural materials. However, the atomic-scale mechanisms which control their plasticity are still only partially understood. At low temperature, the plastic deformation of bcc metals is controlled by the mobility of the screw dislocations, which is hindered by both a strong resistance of the lattice itself, and the presence of other crystal defects, among which are solute atoms.

Atomic-scale models of dislocation mobility based on the Transition State Theory (TST) constitute a useful framework to model plastic flow in pure metals and in alloys. However, the approximations often used (harmonic approximation, constant activation entropy) yield poor predictions in iron. We used the recent projected average force integrator method to compute the activation free enthalpy for kink pair nucleation, including anharmonic effects.

The data show that the harmonic regime is limited to very low temperatures, below 20 K.

Non-linearities remain small below 100 K, allowing to compute an effective activation entropy, which increases when the activation enthalpy decreases, corresponding to an inverse Meyer-Neldel behavior. Integrating these effects in dislocation mobility models greatly improves the agreement with direct molecular dynamics (MD) simulations.

Extensions to Fe-C alloys are limited by the realism of the interatomic potentials available for this system. To address this issue, we combined two existing empirical potentials for Fe and Fe-C to reproduce both the Peierls mechanism and the carbon-induced screw dislocation core reconstruction found in ab initio calculations. Using this hybrid potential, MD simulations of the glide of screw dislocations in random solid solutions confirm a strong solute strengthening, caused by complex short-ranged interaction processes. We also considered an idealized geometry where a screw dislocation interacts with a row of carbon atoms. Combining MD simulations and saddle-point search methods, we unveil a very strong pinning when the solute separation is below about 100 Burgers vectors. This effect is due to the necessity to nucleate two consecutive kink pairs on the screw dislocation, with the second kink pair having a markedly increased activation enthalpy. We developed a harmonic TST model of this process that also integrates the entropic effects observed in pure iron, which yields a vii Contents good agreement with MD simulations conducted up to 300 K. This work provides elementary processes and parameters that will be useful for larger-scale models and in particular kinetic Monte Carlo simulations. viii

Sur les épaules de géants…

Avoir l'opportunité de rédiger une thèse sur un sujet scientifique pointu (et passionnant) est un privilège rare, rendu possible par un environnement exceptionnel. Je suis reconnaissant envers ces inconnus qui sans le savoir m'ont soutenu (et financé) dans un parcours vers le métier qui m'attirait, et m'ont fait confiance pendant trois ans pour tâcher de découvrir quelque chose d'utile. D'autres m'ont fait le privilège supplémentaire de m'y aider sciemment, en me donnant de leur temps, de leur savoir ou de leur soutien ; je veux les remercier dans ces quelques lignes.

En la matière, ma première pensée va à mes encadrants, qui m'ont formé, conseillé, encouragé -et plus encore qui m'ont aidé à tenir la barre face à la pandémie. David : sans surprise, cette thèse contient beaucoup de ce qui fait ton ADN et je suis fier d'avoir fait ce chemin avec toi. Merci pour ta patience et ton courage face à cette thèse qui ne marchait jamais comme on l'attendait, face à ces potentiels qui ne sont jamais là où on les attend. Merci aussi pour ta pédagogie grâce à laquelle j'ai tant appris, ta rigueur à laquelle n'échappe aucun détail et aucun indice, ton envie de toujours aller plus loin. Michel : celui par lequel tout a commencé, qui voit de vraies personnes là où d'autres ne voient que des étudiants un peu perdus. Il est clair que tu as été là à tous les moments clés depuis mes premiers pas dans le supérieur, et que tu as cru en moi dès le début. Tu m'as appris à prendre de la hauteur, à savoir lever la tête du guidon (sans venir aux sorties vélo), à toujours rester humble. Je te dois beaucoup. Fabienne : merci de nous avoir accompagnés tout au long de ce projet. Nos discussions ont été des jalons indispensables, l'occasion de prendre du recul, de se recentrer, de trouver de nouvelles pistes grâce à tes fameuses questions, exigeantes mais toujours pertinentes. Je me souviendrai aussi de la bienveillace dont tu as toujours fait preuve. Je tiens également à remercier Tom, pour ton aide sur une partie aussi imprévue que passionnante, que j'ai beaucoup aimé explorer. Viennent ensuite les compagnons de galère. Merci à mes co-bureaux ; Antoine, qui m'a appris à gravir les paysages énergétiques (comme l'alpiniste de la photo), Manon et Ce pour les souvenirs de notre repas de Noël, Maxime, qui fait toujours des bêtises : arigato pour le scooter et la fanfare asiatique ; Theophile, certifié résistant aux ouragans, merci pour les Jacobins et les paquitos. Merci Laurabelle d'avoir partagé un été de rédaction dans le confort des Algeco. Je dois aussi remercier une bande de docteurs partis avant moi : Alexis, pour les potins sur tu-sais-qui, Aurélien pour les posters, Gwenaëlle pour les galettes des rois et pour ses conseils charte graphique (que je n'ai pas suivi), Alexandre pour les invitations à l'escalade et Justine pour sa bonne humeur inébranlable. Merci à tous ceux qui font de MétAl une équipe si agréable, et pour ne citer qu'eux, à Pierre-Antoine, Patrice et à Jean-Yves (tes mails me manquent déjà), Éric, et aussi aux lutins de Noël. Merci aussi aux inombrables copains dans les autres équipes de MATÉIS , et en particulier à Mathilde qui a tant fait pour nous, et à Carmen, trois mois c'était trop court ! Merci à mes amis qui ne font pas de thèse d'avoir été là pour me rappeler que le monde ne s'est pas arrêté de tourner pendant ce temps. Merci à John pour les water-polo, à Numa pour les flams (et les blattes), à Thibaut pour le duo de sauces, à Camille pour les déguisements, à Alexandre pour la règle des tiers, à Majda pour les chocolats chauds, à Audrey pour les guêpes, et à Clarisse pour le phénix. Merci aussi à Pierre et Romain (sans qui je n'aurais jamais appris Python) d'être là depuis si longtemps, et d'avoir été mes encadrants dans ma 2 ème thèse, celle sur le Covid. Merci à Solen pour le harcèlement téléphonique, pour les recettes healthy, et d'être toujours de mon côté (même quand on n'est pas d'accord).

Pour finir, je remercie bien-sûr ma famille. Merci à mes parents d'avoir toujours cultivé notre curiosité, de l'avoir aidée à mûrir, et d'avoir toujours été là pour m'épauler. Je veux remercier aussi ma soeur, qui me montre la voie dans le monde des grands depuis le début et avec qui je continue de grandir (je te rattraperai !). Merci enfin Marco pour tes conseils ; sur la thèse, mais aussi sur l'extraction du café qui m'a permis de l'écrire.

Résumé

Les aciers ferritiques, constitués d'une matrice de fer cubique centré (CC) avec des solutés de carbone interstitiels, sont des matériaux de structure largement utilisés. Pourtant, les mécanismes qui contrôlent leur plasticité à l'échelle atomique ne sont encore que partiellement compris. A basse température, la déformation plastique des métaux CC est contrôlée par la mobilité des dislocations vis, qui est limitée à la fois par une forte résistance du réseau luimême, et par la présence d'autres défauts cristallins, dont notamment les atomes de soluté.

Les modèles de mobilité des dislocations à l'échelle atomique basés sur la théorie de l'état de transition (TET) fournissent un cadre pour modéliser la déformation plastique des métaux purs et de leurs alliages. Cependant, les approximations couramment utilisées pour appliquer ces modèles (approximation harmonique, entropie d'activation constante) donnent de médiocres résultats dans le fer. Nous avons calculé l'enthalpie libre d'activation pour la nucléation de paires de crans grâce à la récente méthode projected average force integrator, qui permet de prendre en compte les effets anharmoniques. Les données indiquent que le régime harmonique est limité à des températures très basses, inférieures à 20 K. Les non-linéarités restent faibles en dessous de 100 K, permettant de calculer une entropie d'activation effective, qui augmente lorsque l'enthalpie d'activation diminue, ce qui correspond à un comportement Meyer-Neldel inverse. L'intégration de ces effets dans les modèles de mobilité des dislocations améliore considérablement l'accord avec les simulations de dynamique moléculaire (DM).

L'application de ces approches aux alliages Fe-C demeure limitée par le réalisme des potentiels interatomiques disponibles pour ce système. Pour résoudre ce problème, nous avons combiné deux potentiels empiriques existants pour Fe et Fe-C afin de reproduire à la fois le mécanisme de Peierls et la reconstruction des coeurs de dislocation vis en présence de carbone mise en évidence par des calculs ab initio. Des simulations de DM du glissement d'une dislocation vis dans une solution solide aléatoire utilisant ce potentiel confirment un fort durcissement, causé par des processus complexes d'interaction à courte distance entre dislocation et solutés. Nous avons également étudié une géométrie idéalisée où une dislocation vis interagit avec une rangée d'atomes de carbone. En combinant des simulations de DM et des méthodes de recherche de point-col, nous dévoilons un très fort ancrage lorsque ix Contents la séparation des solutés est inférieure à environ 100 vecteurs de Burgers. Cet effet est dû à la nécessité de nucléer deux paires de cran consécutivement sur la dislocation, la deuxième paire ayant une enthalpie d'activation nettement augmentée. Nous avons développé un modèle de ce processus basé sur la TET harmonique, qui intègre également les effets entropiques observés dans le fer pur, en bon accord avec les simulations de DM menées jusqu'à 300 K.

Ce travail fournit des processus et des paramètres élémentaires qui seront utiles pour des modèles à plus grande échelle et en particulier des simulations Monte Carlo cinétique.

x 0 Introduction

Ageing of nuclear reactor vessels

Nuclear power is widely regarded as a low-carbon emission method of electricity generation [START_REF] Lenzen | Life cycle energy and greenhouse gas emissions of nuclear energy: A review[END_REF][START_REF] Warner | Life Cycle Greenhouse Gas Emissions of Nuclear Electricity Generation[END_REF]. Compared to fossile fuels, which are used for 60% of worldwide electricity production (see Fig. 0.1) and account for a large part of global greenhouse gas emissions, nuclear energy could contribute to the mitigation of global climate change and air pollution [START_REF] Kharecha | Prevented mortality and greenhouse gas emissions from historical and projected nuclear power[END_REF]. While the benefit of new investments in nuclear power compared to renewable energies is still debated [START_REF] Sovacool | Differences in carbon emissions reduction between countries pursuing renewable electricity versus nuclear power[END_REF], it has become clear that leveraging existing nuclear power capacity can help reach greenhouse gas reduction goals. However, this might be challenged by the limited lifetime of nuclear reactors, often estimated as 50-60 years [START_REF] Ballesteros | The role of pressure vessel embrittlement in the long term operation of nuclear power plants[END_REF]). In France for example, existing reactors have been exploited for 19 to 41 years, with a median age of 36 years (IAEA, 2019). Reactor closure would result in a reduction in electric production capacity, likely compensated by a massive use of fossile fuels -which would in turn increase greenhouse gas emissions. Instead, operation should be prolonged to more than 50 years to buy time as energetic transition and renewable deployment are ramping up [START_REF] Mainsant | Trajectoires d 'évolution du mix électrique 2020-2060[END_REF].

The possiblity of extending the lifetime of existing nuclear power plants has thus become a burning issue in the face of climate change emergency [START_REF] Pachauri | Climate change 2014: synthesis report. Contribution of Working Groups I, II and III to the fifth assessment report of the Intergovernmental Panel on Climate Change[END_REF] and ubiquitous calls to cut greenhouse gas emissions.

The main component of pressurized water reactors (PWR) is their vessel, made of 330-520 tons of ferritic steel (depending on the generation of the reactor) and housing the reactor core, which cannot be replaced at a reasonable cost. The integrity of the vessel can therefore be the limiting feature for the PWR operation lifetime. The vessel is also an essential safety component, as it acts as the second barrier of containment of radioactive material. One key aspect is thus to make sure that safety margins used for PWRs design remain valid, including in the conditions of an accident. Material properties can be computed experimentally and their evolution monitored, but investigating all relevant conditions is costly and difficult, es- pecially for accidental conditions. The use of physical models to extend empirical approaches is therefore growing (CEA-DEN, 2014).

Vessel steel is subject to thermal ageing and irradiation during operation, which cause strengthening, i.e. an increase of the yield stress, and an embrittlement that reduces its ability to deform plastically under a high mechanical load. The yield stress increase and loss of ductility that characterize irradiation hardening are illustrated in Fig. 0.2. The creation of irradiation defects involves a number of mechanisms at the atomic scale: the impact of an energetic particule results in a cascade of atomic displacements, and the formation of residual point defects (i.e. the primary damage) [START_REF] Becquart | Modelling the primary damage in Fe and W: Influence of the short range interactions on the cascade properties: Part 1 -Energy transfer[END_REF]; new defects interact with existing microstructural defects (e.g. dislocations, solutes), possibly affecting their mobility [START_REF] Zamzamian | A study on the mobility of 12111{01¯1} edge dislocation in low-carbon 𝛼-Fe and its interactions with damage cascade: On picosecond time scale using molecular dynamics simulations[END_REF]; and the accumulation or diffusion of these defects can create strong obstacles to dislocation motion (e.g. precipitates, dislocation loops) that cause hardening [START_REF] Tipping | Understanding and mitigating ageing in nuclear power plants: materials and operational aspects of plant life management[END_REF]. On the other hand, thermal ageing allows the diffusion of interstitial solutes, such as carbon in iron, to crystal defects and notably to dislocations, where they form Cottrell atmospheres [START_REF] Cottrell | Dislocation theory of yielding and strain ageing of iron[END_REF][START_REF] Wilde | Three-dimensional atomic-scale mapping of a Cottrell atmosphere around a dislocation in iron[END_REF] that limit dislocations mobility [START_REF] Berns | Ferrous Materials[END_REF]. While Cottrell atmosphere modelling has a long history [START_REF] Cottrell | Dislocation theory of yielding and strain ageing of iron[END_REF][START_REF] Zhao | Kinetics of Cottrell atmosphere formation during strain aging of ultra-low carbon steels[END_REF], and their formation in steels has been studied by several authors at the atomic scale [START_REF] Veiga | Atomistic modeling of carbon Cottrell atmospheres in bcc iron[END_REF][START_REF] Waseda | Formation of carbon Cottrell atmospheres and their effect on the stress field around an edge dislocation[END_REF][START_REF] Candela | Investigating the kinetics of the formation of a C Cottrell atmosphere around a screw dislocation in bcc iron: a mixed-lattice atomistic kinetic Monte-Carlo analysis[END_REF], an atomistically-informed model of strengthening by interstitial solutes is still lacking. Understanding the processes by which strengthening emerges is especially relevant in order to prevent failure in accidental conditions. Important modelling efforts are therefore devoted to their integration in plasticity models, which involves various techniques at different scales [START_REF] Lu | An Overview of Multiscale Simulations of Materials[END_REF].

(a) Yield stress augmentation associated with thermal ageing for different times and temperatures. Reproduced from [START_REF] De | Effect of dislocation density on the low temperature aging behavior of an ultra low carbon bake hardening steel[END_REF].

(b) Stress-strain curves for different irradiation conditions. Unirr. refers to the unirradiated material. Reproduced from [START_REF] Tipping | Understanding and mitigating ageing in nuclear power plants: materials and operational aspects of plant life management[END_REF] Figure 0.2: Effect of thermal ageing and irradiation on ferritic steels.

Multi-scale modelling approach of plasticity

The development of plasticity models aims at understanding the relations between the microstructure and mechanical properties of metals. Continuous models of crystal plasticity (CP), often based on the Finite-Element Method (FEM), are widely-used engineering tools, able to predict the mechanical properties in conditions that were not necessarily tested experimentally using general constitutive laws. But the accuracy of such approaches can be limited [START_REF] Pinna | Assessment of crystal plasticity finite element simulations of the hot deformation of metals from local strain and orientation measurements[END_REF][START_REF] Mello | Effect of microstructure on strain localization in a 7050 aluminum alloy: comparison of experiments and modeling for various textures[END_REF] and they remain phenomenological [START_REF] Tipping | Understanding and mitigating ageing in nuclear power plants: materials and operational aspects of plant life management[END_REF].

To improve the reliability of these models, connections should be drawn with lower-scale plasticity models [START_REF] Van Der Giessen | Roadmap on multiscale materials modeling[END_REF] which are based on the deformation mechanisms allowed by crystal defects -among which are dislocations, which control plastic deformation under most conditions [START_REF] Kubin | Dislocations, Mesoscale Simulations and Plastic Flow[END_REF]. As dislocation physics spans a wide range of length and timescales, there is no unique technique able to model dislocations from the scale of their core up to the polycrystalline microstructural scale, but a set of discrete and continuous numerical methods which together form a multi-scale scheme, presented in Fig. 0.3. It is to be noted that several of these methods can be coupled within a single simulation, e.g. using continuous methods in some regions to reduce the computational load and number of degrees of freedom associated with discrete methods. Prominent examples of coupling methods include the quasicontinuum method [START_REF] Shenoy | An adaptive finite element approach to atomic-scale mechanics -The quasicontinuum method[END_REF][START_REF] Miller | The Quasicontinuum Method: Overview, applications and current directions[END_REF][START_REF] Kochmann | The Quasicontinuum Method: Theory and Applications[END_REF], Coupled Atomistic/Discrete-Dislocation (CADD) models [START_REF] Shilkrot | Coupled Atomistic and Discrete Dislocation Plasticity[END_REF][START_REF] Shilkrot | Multiscale plasticity modeling: Coupled atomistics and discrete dislocation mechanics[END_REF] or the Macroscopic, Atomistic, Ab initio Dynamics (MAAD) method [START_REF] Abraham | Spanning the continuum to quantum length scales in a dynamic simulation of brittle fracture[END_REF].

At the mesoscopic scale, Discrete Dislocations Dynamics (DDD) methods represent dislocations as discrete lines, while Continuous Dislocations Dynamics (CDD) methods model a continuous dislocation density field. Both can be used to model dislocation networks, explaining collective phenomena such as strain hardening [START_REF] Sills | Dislocation networks and the microstructural origin of strain hardening[END_REF], and to propose constitutive laws usable in CP models. The velocity of dislocations in these simulations is classically given by laws derived from atomistic simulations of dislocations conducted with classical Molecular Dynamics (MD). As MD simulations operate at a much lower scale and with up to 10 10 times more degrees of freedom than DDD [START_REF] Bertin | Frontiers in the Simulation of Dislocations[END_REF], velocity laws are obtained from simulations of single dislocations -with an already large computational load-and the interaction between dislocations is thus not taken into account at the atomistic level. Yet, by leveraging important computational resources, it has become possible to use MD to simulate ∼ 10 3 times larger systems than single dislocations, to compare the predictions of DDD and MD at the same scale (≈ 1µm) [START_REF] Zepeda-Ruiz | Atomistic insights into metal hardening[END_REF], by identification of the dislocation network from atom positions using "in silico microscopy" techniques [START_REF] Stukowski | Extracting dislocations and non-dislocation crystal defects from atomistic simulation data[END_REF]. But beyond the scarcity of High Performance Computing (HPC) platforms suited for this kind of approach, the treatment of the much more complex dislocations-obstacles interaction mechanisms would require extravagant computational resources.

Instead, a more classical approach uses a limited set of atomistic simulations of dislocations in either pure metals or alloys and/or in the presence of other defects to parameterize lightweight analytical [START_REF] Rodney | Activation enthalpy for kink-pair nucleation on dislocations: Comparison between static and dynamic atomic-scale simulations[END_REF][START_REF] Domain | Simulation of Screw Dislocation Motion in Iron by Molecular Dynamics Simulations[END_REF][START_REF] Maresca | Theory of screw dislocation strengthening in random BCC alloys from dilute to "High-Entropy" alloys[END_REF] or Monte Carlo [START_REF] Stukowski | Thermally-activated non-Schmid glide of screw dislocations in W using atomistically-informed kinetic Monte Carlo simulations[END_REF][START_REF] Zhao | Direct prediction of the solute softening-to-hardening transition in W-Re alloys using stochastic simulations of screw dislocation motion[END_REF] mesoscale kinetic models of dislocation velocity. Putting the different mechanisms identified at the atomic scale in competition in mesocale models is especially useful to discriminate the most important ones, which should then be integrated in DDD or CP models [START_REF] Cereceda | Unraveling the temperature dependence of the yield strength in single-crystal tungsten using atomistically-informed crystal plasticity calculations[END_REF]. The accuracy of atomistically-informed models is fundamentally limited by the interatomic potential (IP) used to compute forces in atomistic simulations [START_REF] Bulatov | Computer Simulations of Dislocations[END_REF][START_REF] Van Der Giessen | Roadmap on multiscale materials modeling[END_REF], which can be only partially validated on experimental data [START_REF] Alexander | Interatomic potentials for irradiation-induced defects in iron[END_REF] and has to be built on large ab initio databases [START_REF] Cereceda | Assessment of interatomic potentials for atomistic analysis of static and dynamic properties of screw dislocations in W[END_REF] to achieve good accuracy and transferability.

The widely-used Density Function Theory method (DFT) is capable of modelling atomic interactions at the electronic level. It was used to model dislocations in a number of metals and semiconductors [START_REF] Rodney | Ab initio modeling of dislocation core properties in metals and semiconductors[END_REF], shedding light on small-scale processes in the region of the core that are crucial for dislocation mobility in pure metals [START_REF] Woodward | Flexible Ab Initio Boundary Conditions: Simulating Isolated Dislocations in bcc Mo and Ta[END_REF][START_REF] Itakura | First-principles study on the mobility of screw dislocations in bcc iron[END_REF][START_REF] Ventelon | Ab initio investigation of the Peierls potential of screw dislocations in bcc Fe and W[END_REF][START_REF] Dezerald | Modélisation ab initio des dislocations vis dans les métaux de transition cubiques centrés[END_REF][START_REF] Dezerald | Plastic anisotropy and dislocation trajectory in BCC metals[END_REF][START_REF] Kraych | Nonglide effects and dislocation core fields in BCC metals[END_REF] and in alloys, including the Fe-C system and other interstitials [START_REF] Ventelon | Dislocation core reconstruction induced by carbon segregation in bcc iron[END_REF][START_REF] Lüthi | First principles investigation of carbon-screw dislocation interactions in body-centered cubic metals[END_REF][START_REF] Lüthi | Attractive interaction between interstitial solutes and screw dislocations in bcc iron from first principles[END_REF][START_REF] Hachet | Screw dislocation-carbon interaction in BCC tungsten: an ab initio study[END_REF][START_REF] Romaner | Effect of Rhenium on the Dislocation Core Structure in Tungsten[END_REF], which have to be transferred to larger scales. Indeed, the scale accessible to DFT calculations lies under a thousand atoms in static simulations, and some elementary processes of dislocation motion -such as the kink pair nucleation in b.c.c. metals [START_REF] Proville | Modeling the thermally activated mobility of dislocations at the atomic scale[END_REF]-cannot be captured with this method; dynamical simulations (i.e. ab initio molecular dynamics) are also not applicable. One option is to parameterize higher scale models from DFT, such as line tension models [START_REF] Proville | Prediction of the kink-pair formation enthalpy on screw dislocations in 𝛼-iron by a line tension model parametrized on empirical potentials and first-principles calculations[END_REF]; [START_REF] Dezerald | First-principles prediction of kink-pair activation enthalpy on screw dislocations in bcc transition metals[END_REF]; [START_REF] Hachet | Screw dislocation-carbon interaction in BCC tungsten: an ab initio study[END_REF], velocity laws [START_REF] Hachet | Mobility of screw dislocation in BCC tungsten at high temperature in presence of carbon[END_REF], or kinetic Monte Carlo models [START_REF] Zhao | Simulating the mechanisms of serrated flow in interstitial alloys with atomic resolution over diffusive timescales[END_REF]. Classical MD simulations can also be used, while ensuring they are meaningfully connected to other scales, i.e. that the IP used has a reasonable accuracy compared to ab initio data, and that the obtained predictions are transferable to larger scales and other conditions (temperature, deformation rate, 0 Introduction defects concentration,…), for example with the help of an analytical model based on rate theory [START_REF] Pollak | Reaction rate theory: What it was, where is it today, and where is it going?[END_REF]. The connections between different simulation techniques used for modelling plasticity are illustrated in Fig. 0.4, with transfers of information between techniques.

Figure 0.4: Flowchart of the multi-scale modelling approach illustrating information transfer between techniques. Reproduced from Tipping (2010).

Towards a mobility law for dislocations in aged steel

Recent DFT results on the interaction between interstitials and dislocations mentioned above offer promising perspectives to better understand how solutes (either in a Cottrell atmosphere or a solid solution) can affect the mobility of dislocations, and ultimately the mechanical properties of steel. However, the results obtained in DFT simulations are limited to the case of high solute concentrations, while the random case of dilute interstitial alloys cannot be treated with DFT, notably because large simulation cells are needed to model kink pairs. The adaptation of existing models of solute strengthening developed for substitutional alloys [START_REF] Varvenne | Solute strengthening in random alloys[END_REF][START_REF] Maresca | Theory of screw dislocation strengthening in random BCC alloys from dilute to "High-Entropy" alloys[END_REF][START_REF] Rao | Theory of solid solution strengthening of BCC Chemically Complex Alloys[END_REF] to the case of interstitial bcc alloys with large solute-dislocation interaction energies also remains difficult, as an underlying assumption is that solutes do not significantly distort the dislocation 0.3 Towards a mobility law for dislocations in aged steel core [START_REF] Varvenne | Solute strengthening in random alloys[END_REF], while the opposite is suggested by DFT calculations in interstitial alloys [START_REF] Ventelon | Dislocation core reconstruction induced by carbon segregation in bcc iron[END_REF][START_REF] Lüthi | Attractive interaction between interstitial solutes and screw dislocations in bcc iron from first principles[END_REF]. The treatment of interstitials diffusion close to the dislocation core at high temperature also remains inaccessible to such models.

This motivates a detailed investigation of solute-dislocation interactions at the atomic scale in the Fe-C system, to include ab initio insights into mesoscale dislocation dynamics models.

The aim is to develop a dislocation mobility law in the presence of C, and to extract local rules of interaction with the defects from atomistic simulations (see Fig. 0.4). To achieve this, a number of challenges must be overcome: 1 Methods

Atomic interactions

Atomistic simulations

In this project, our goal is to use molecular dynamics (MD) to simulate the motion of screw dislocations in bcc Fe and Fe-C alloys, with a ⃗ 𝑏 = 1 2 [111] Burgers vector, gliding in a ( 1 10) plane. Due to the usual tradeoff between model size and computational time in atomistic simulations, the simulation of dislocations, which are extended line defects able to interact at long range, involves many technical concerns. General issues regarding the construction of the dislocation model are discussed in this chapter. Chap. 2 is devoted to a discussion of interatomic potentials.

All atomic system visualizations were produced using the Ovito software [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITO-the Open Visualization Tool[END_REF]. Dislocation line positions were obtained from atomic positions using the DXA algorithm available in Ovito [START_REF] Stukowski | Extracting dislocations and non-dislocation crystal defects from atomistic simulation data[END_REF].

Molecular dynamics simulations

To run MD simulations, we use the widely-adopted LAMMPS package [START_REF] Plimpton | Fast Parallel Algorithms for Short-Range Molecular Dynamics[END_REF], as it offers numerous capabilities, has great parallel performance on CPU and GPU hardwares and is a standard in the community.

All dynamical simulations are conducted in the NVE ensemble, i.e. using the fix nve LAMMPS command. During dislocation motion, a temperature rise of a few kelvins is therefore observed due to the plastic work, and is included in our reported data. Dynamical simulations are run using the GPU package of LAMMPS [START_REF] Brown | Implementing molecular dynamics on hybrid high performance computers-short range forces[END_REF] to accelerate force evaluations. For dynamical simulation, LAMMPS was built based on the 19 Sep. 2019 release, with MANYBODY, REPLICA and GPU packages enabled and support for the GZIP library and openmpi. For static benchmarking of potentials (Chap. 2), we built the 4 Jan. 2019 release with the REPLICA, MANYBODY and USER-MEAMC packages.

1 Methods

Static simulation of dislocations

As the DFT method is highly computationally expensive, the number of simulated atoms must not exceed a few hundred. Different methods that allow static simulations of dislocations in small systems can be used [START_REF] Rodney | Ab initio modeling of dislocation core properties in metals and semiconductors[END_REF].

The most straightforward approach, called the cluster approach, consists in defining a cylinder of atoms, oriented along the [111] direction of the crystal. A screw dislocation is inserted in the same direction at the centre of the crystal, with periodic boundaries along the line direction. The dimension in this direction can be set to 1 Burgers vector for an infinite straight dislocation. This method has several issues due to the limited size of the simulated domain, notably regarding the elastic field of the core, the change of boundary conditions seen by the dislocation when it moves, and the numerical efficiency due to the presence of vacuum around the cylinder. More sophisticated approches use a combination of elasticity and interatomic potentials to implement flexible boundary conditions, where atoms simulated by DFT are surrounded by a wider domain where other simulation methods are applied. But this introduces several technical complications, and does not allow to extract dislocation energies (see Ref. [START_REF] Rodney | Ab initio modeling of dislocation core properties in metals and semiconductors[END_REF] for a detailed discussion).

Instead, to avoid external boundary issues, a method allowing the use of a tri-periodic cell was developed [START_REF] Bigger | Atomic and electronic structures of the 90°partial dislocation in silicon[END_REF][START_REF] Ismail-Beigi | Ab Initio Study of Screw Dislocations in Mo and Ta: A New Picture of Plasticity in bcc Transition Metals[END_REF]. To have fully periodic conditions, a dislocation dipole of opposed Burgers vectors is inserted to obtain a null total Burgers vector. By periodicity, it is representative of a 2D square lattice of dislocations quadrupoles. This arrangement has the advantage to cancel the Peach-Koehler image force on the dislocations. The quadrupolar simulation cell is illustrated in Fig. 1.1. To construct this square arrangement of dislocations of opposed burgers vectors, inserted in two equivalent crystallographic positions (denoted as easy and described in Chap. 3), periodicity vectors should respect geometrical constraints as described in Refs. [START_REF] Ventelon | Core structure and Peierls potential of screw dislocations in 𝛼-Fe from first principles: cluster versus dipole approaches[END_REF][START_REF] Ventelon | Ab initio investigation of the Peierls potential of screw dislocations in bcc Fe and W[END_REF]. A setup that matches these constraints and typically used in ab initio simulations is a cell of length 𝑛 × 𝑏 in the [111] direction, which contains 𝑛 × 135 atoms, with 𝑛 being an integer. By construction, the dimensions in [ 1 1 2] and [1 1 0] directions are respectively 15 × √ 2/3𝑎 0 and 9 × √2/2𝑎 0 . If a dislocation of the dipole moves (e.g. in simulations of dislocation glide), the elastic energy is affected by both the elastic interaction of the dipole, and the deformation associated with the motion of the dislocation. Elastic corrections based on anisotropic linear elasticity can be used to deduce these effects from the total energy, and can be applied using the Babel package1 . The interested reader is referred to [START_REF] Clouet | Ab Initio Models of Dislocations[END_REF] 1.1 Atomistic simulations for a review on ab initio methods for the modelling of dislocations. This simulation setup will be used to compare properties computed with interatomic potentials to ab initio results in representative conditions. [START_REF] Ventelon | Simulation ab initio des coeurs de dislocation vis dans le fer[END_REF].

Dynamical simulation of dislocations

For dynamical simulations, it is preferable to avoid fully periodic simulation cells, as the strong attraction between the dislocations of the dipole may result in cross-slip and dislocation annihilation.

Instead, we used a cell represented in Fig. 1.2, where we introduced a single screw dislocation along the ⃗ 𝑥 = [111] direction using the isotropic elastic displacement solution [START_REF] Hirth | Theory of dislocations[END_REF] followed by an energy minimization. Periodic boundary conditions are applied in the 𝑥𝑦 ( 1 10) glide plane. Free surfaces are created by extending the cell dimension in the ⃗ 𝑧 = [ 1 10] direction by more than the cutoff radius of the interatomic potential in use, such that atoms in opposite boundaries do not interact. An additional tilt of 𝑏/2 is applied in the ⃗ 𝑦 direction, to account for the plastic strain of the dislocation, resulting in a triclinic cell. These conditions generate a periodic array of dislocations (PAD) [START_REF] Rodney | Molecular dynamics simulation of screw dislocations interacting with interstitial frank loops in a model FCC crystal[END_REF][START_REF] Bacon | Chapter 88 Dislocation-Obstacle Interactions at the Atomic Level[END_REF].

Methods

Different dimensions in the ⃗ 𝑥 direction 𝐿 𝑥 were tested, while dimensions along ⃗ 𝑦 and ⃗ 𝑧 and refered as 𝐿 𝑦 and 𝐿 𝑧 were set to (138 Å, 159 Å) in dynamic simulations (corresponding to 192000 atoms if 𝐿 𝑥 = 40𝑏) and (138 Å, 80 Å) at zero kelvins (96000 atoms). At finite temperature, a larger 𝑧 dimension is necessary to lower the attractive effect of free 𝑧 surfaces on the dislocation, which induces cross-slip. Dynamical simulations were performed under constant strain rate, using the flexible boundary conditions as presented in [START_REF] Rodney | Activation enthalpy for kink-pair nucleation on dislocations: Comparison between static and dynamic atomic-scale simulations[END_REF]. For this purpose, two slabs of atoms denoted 𝑆 + and 𝑆 -are defined underneath the top and bottom 𝑧-surfaces, with a thickness larger than the interatomic potential cut-off radius (approximately 6 Å for classical EAM potentials). An initial velocity is applied to all atoms of each slab, depending on the desired 𝑥𝑧 shear rate γ :

𝐿 𝑦 𝐿 𝑥 𝐿 𝑧 ⃗ 𝑦 = [ 1 1 2] ⃗ 𝑥 = [111] ⃗ 𝑧 = [1 1 0] S + S -
𝑣 ± 0,𝑥 = ± γ 𝐿 𝑧 2 . (1.1)
At every time step, the force in the ⃗ 𝑥 = [111] direction on each atom 𝑖 in S ± is corrected

1.1 Atomistic simulations in order to ensure that the total force in S + and S -is zero:

𝐹 𝑖 𝑥 ← 𝐹 𝑖 𝑥 -⟨𝐹 𝑥 ⟩ ± , (1.2)
where ⟨𝐹 𝑥 ⟩ ± refers to the average force in S ± in the ⃗ 𝑥 direction. In this way, the centres of mass of the slabs move at a constant velocity, thus imposing a constant strain rate, while allowing atoms in S ± to adapt to the motion of the dislocation [START_REF] Rodney | Activation enthalpy for kink-pair nucleation on dislocations: Comparison between static and dynamic atomic-scale simulations[END_REF] A caveat of this type of simulation cell is the use of different boundary conditions in the 𝑥, 𝑦 (both periodic) and 𝑧 (free surfaces) directions. As the crystal is allowed to change its 𝐿 𝑧 dimension while 𝐿 𝑥 and 𝐿 𝑦 are constrained, thermal expansion is possible along 𝑧 and blocked along other directions, resulting in an anisotropic residual stress.

Remark The standard procedure presented above might be improved by imposing the 𝑧 coordinate of the centre of mass of 𝑆 ± , enabling an equivalent treatment of the thermal expansion in the 𝑥, 𝑦, and 𝑧 directions. This can be done by susbstracting the average force in the atoms group to each atom of the group:

𝐹 𝑖 𝑧 ← 𝐹 𝑖 𝑧 -⟨𝐹 ± 𝑧 ⟩. (1.3)
Which corresponds to the following LAMMPS snippet: compute average_fz_s1 atoms_slab_1 reduce ave fz variable fz_to_add_s1 equal -c_average_fz_s1 3 fix freeze_s1_z atoms_slab_1 addforce 0 0 v_fz_to_add_s1

Listing 1.2: LAMMPS code snipppet used to prevent atoms movement in z direction.

Computational techniques

The rate at which a system transits from one minimum energy state to the next by thermal activation (e.g. a dislocation moving through a bcc lattice) is the key component for a thermally activated velocity law, but its computation is notoriously difficult and requires assumptions [START_REF] Ásgeirsson | Exploring Potential Energy Surfaces with Saddle Point Searches[END_REF].

The transition state theory (TST) [START_REF] Wigner | The transition state method[END_REF] allows to express the transition rate 𝑝 between minimum energy states (MES). In its most common form, it writes:

𝑝 ∝ exp ( -Δ𝐺 𝑘𝑇 ) , (1.4)
where the Gibbs energy of activation Δ𝐺 is the difference between the Gibbs energy of the initial state 𝐺 0 and the transition state 𝐺 † .

Within the widely used harmonic approximation of atom vibrations [START_REF] Vineyard | Frequency factors and isotope effects in solid state rate processes[END_REF], the temperature dependence of the enthalpy Δ𝐻 and entropy Δ𝑆 is neglected. Δ𝐺 can then be decomposed as Δ𝐺 = Δ𝐻 -𝑇 Δ𝑆, yielding:

𝑝 ∝ exp ( Δ𝑆 𝑘 ) exp (- Δ𝐻 𝑘𝑇 ) , (1.5)
Δ𝐻 and Δ𝑆 can be computed separately at zero Kelvin using the techniques described in sections 1.2.1 and 1.3 respectively.

While the harmonic approximation is well verified at low temperature, the behaviour at higher temperature can strongly differ from the harmonic prediction due to anharmonic effects [START_REF] Swinburne | Uncertainty and anharmonicity in thermally activated dynamics[END_REF][START_REF] Sato | Anharmonic effect on the thermally activated migration of {1012} twin interfaces in magnesium[END_REF]. A recently developed technique [START_REF] Swinburne | Unsupervised calculation of free energy barriers in large crystalline systems[END_REF] allows the direct computation of the Gibbs energy of activation and is usable on systems containing millions of atoms, with no assumption on thermal vibrations. This method is described in Section 1.3.1 and applied to the case of dislocation glide in iron in Chap 3.

Harmonic enthalpy calculation

Saddle-point search techniques are a class of numerical methods used to localise TS and compute their activation enthalpy. In the case where the TS is located between two known MES, double-ended search techniques can be used to find the minimum energy paths (MEP) connecting both2 . This is the case for the study of dislocation glide in the ( 1 01) plane, where the known MES are adjacent Peierls valleys.

The CI-NEB method

The Nudged Elastic Band (NEB) method [START_REF] Henkelman | A climbing image nudged elastic band method for finding saddle points and minimum energy paths[END_REF] allows the convergence to the MEP that is closest to an initial guess. This initial path is discretized in a series of successive atomic states called replicas, which are usually obtained by linear interpolation of atoms positions in cartesian coordinates between two known MES, designated as initial and final. This simple approach is sufficient to provide an initial guess that is close to a MEP in many cases, but is not optimal for the kink-pair mechanism. In this case, MES are infinite straight dislocation lines in two adjacent Peierls valleys, symmetrical in the direction of the line (i.e. the [111] direction). A linear interpolation between the two states preserves this symmetry, which then has to be broken during the NEB calculation, thus increasing the number of steps.

A more efficient approach is to explicitly provide an initial path where the transition goes through the formation of a kink pair, as shown in Fig. 1.3. To do so, we slice the cell in three regions, defined by cutting planes that are perpendicular to the ⃗ 𝑥 ≡ [111] direction. The first region has a thickness 𝑖/𝑅 × 𝐿, where 𝑖 is the index of the replica starting at 0, R the total number of replicas and L the cell dimension along ⃗ 𝑥. In this region, atomic positions are taken from the final state, corresponding to the screw segment that has crossed the barrier.

A buffer region of thickness 𝑙 0 makes a smooth transition to the rest of the cell where the dislocation is in its initial position. In the presence of carbon atoms along the dislocation line, a linear interpolation between the initial and final state is sufficient, as kinks are present in either of both states (see Chap. 4 and especially Fig. 4.10 for more details).

After the initial path is defined, the replicas are relaxed iteratively in the direction of the energy gradient projected in the hyperplane perpendicular to the path, allowing the convergence to the MEP. In order to control the distribution of replicas along the MEP, an harmonic bias spring force is added between adjacent images [START_REF] Henkelman | A climbing image nudged elastic band method for finding saddle points and minimum energy paths[END_REF]. This procedure,

1
𝐹 𝑐 = -∇𝐸 + 2(∇𝐸 ⋅ τ 𝑐 ) ⋅ τ 𝑐 (1.6)
where ∇𝐸 is the energy gradient with respect to atoms coordinates and τ 𝑐 a normalized vector tangent to the path. A classical optimization algorithm is then used until the force has converged. Some of the limitations of the CI-NEB method are the sensitivity of the final MEP to the strength of the bias spring force [START_REF] Kolsbjerg | An automated nudged elastic band method[END_REF], the high computational cost of increasing the number of replicas along the path, or the addition of the bias force, which makes the force convergence ill-defined. In the case of dislocation unpinning from carbon presented in Chap. 4, complications arise for defining the position of carbon final state. We therefore preferred to use CI-NEB in combination of a single-end search method, the ARTn technique 1.2 Computational techniques described below, to automatically determine a final state which is connected by a simple barrier to the initial state. The correct determination of the final state is especially important in this case as we used it as the initial state of another calculation (see Chap. 4 for details).

NEB under an applied stress

The NEB method can be used to compute the enthalpy barrier between two positions of a dislocation, under an applied stress. For this, the LAMMPS command fix addforce can be used to add an arbitrary force on a group of atoms. Plastic work is produced as the dislocation moves along the 𝑦 coordinate, which must be included in the total energy to obtain the enthalpy. This is done by using the fix_modify energy keyword. A snippet illustrating how to impose an external stress on a slab of atoms (as of LAMMPS ver. When performing such calculations, we faced difficulties due to a bug in the NEB routine of LAMMPS. When loading pre-defined atomic replicas (neb ... each command), the periodic boundary conditions were not enforced, resulting in energy jumps when atoms crossed 1 Methods the cell boundaries. After we identified the bug, it was corrected in the main LAMMPS distribution in its 30 Nov. 2020 version3 . As a workaround before the bug was corrected, we used an in-house MD code ("MERLIN") which provides an alternative implementation of standard MD techniques, including the NEB method. The code was previously used in several studies on dislocations [START_REF] Rodney | Activation enthalpy for kink-pair nucleation on dislocations: Comparison between static and dynamic atomic-scale simulations[END_REF][START_REF] Proville | Quantum effect on thermally activated glide of dislocations[END_REF]. In the case of enthalpy barriers computed under an applied stress, we present the results obtained with this implementation to avoid re-running computationally intensive calculations. To allow the reproduction of our results, the equivalent LAMMPS directives presented in Listing 1.3 can be used in the patched version of the package. In our calculations using the MERLIN code, we used the ARTn technique (presented below) in addition to the CI-NEB method, to avoid difficulties in the definition of the final state.

The ARTn method

The activation relaxation technique nouveau (ARTn) technique [START_REF] Barkema | Event-based relaxation of continuous disordered systems[END_REF][START_REF] Malek | Dynamics of Lennard-Jones clusters: A characterization of the activation-relaxation technique[END_REF] is often used to explore the energy landscape around a known MES [START_REF] Mousseau | The Activation-Relaxation Technique: ART Nouveau and Kinetic ART[END_REF], searching for saddle-points that connect to unknown MES. It thus requires the knowledge of only one MES, which is highly valuable when final states cannot be easily guessed. This method consists in three main steps:

1. Starting from the MES, the system is iteratively pushed in a random direction in order to leave the harmonic bassin. After each iteration, a relaxation in the hyperplane perpendicular to the random direction is performed, and the lowest eigenvalue 𝜆 0 of the Hessian matrix is computed using the Lanczos method [START_REF] Malek | Dynamics of Lennard-Jones clusters: A characterization of the activation-relaxation technique[END_REF].

The system is considered to have left the harmonic basin when 𝜆 0 is negative (or inferior to a negative threshold set by the user).

2. The activation step, or minimum mode following step, is then performed to converge to the activated state. The system is pushed in the direction of the saddle point following the direction of the eigenvector ê 0 associated with 𝜆 0 , and then relaxed in the hyperplane perpendicular to ê 0 . This operation is repeated until the total force is below a preset force tolerance.

3. Starting from the found saddle state, the system is pushed in the direction of ê 0 and then fully relaxed to obtain the final state. In order to verify that the saddle state is 1.3 Harmonic entropy calculation directly connected to the initial state, the saddle state system is also pushed along -ê 0 and relaxed.

When using the MERLIN code (i.e. only for enthalpy calculations under an applied stress), we used the ARTn technique in combination to the NEB method, as follows. CI-NEB was first performed under an applied stress to obtain an estimate of the direction of the saddle state and the overall shape of the barrier, using 250 NEB integration steps. The estimated direction of the saddle state is used to push the system out of the harmonic bassin in the first step of the ARTn technique, instead of a random vector. The ARTn method is then used as described above to determine the saddle point with a force tolerance set to 2.5 × 10 -5 eV Å -1 .

Harmonic entropy calculation

The effect of entropy is taken into account by TST models, but computing its value is often difficult. Within the harmonic approximation, the vibrational entropy can be linked to the frequency of each vibration mode 𝑖 in the initial 𝜈 0 𝑖 and the activated state 𝜈 † 𝑖 , obtained by diagonalizing the Hessian matrix in the initial and activated states. The transition rate is then written [START_REF] Vineyard | Frequency factors and isotope effects in solid state rate processes[END_REF]:

𝑝 = 𝜈 ⋆ ∏ 𝐷 𝜈 0 𝑖 ∏ 𝐷-1 𝜈 † 𝑖 exp ( -Δ𝐻 𝑘𝑇 ) , (1.7)
where 𝜈 ⋆ is a constant that accounts for the attempt frequency (of the order of the Debye frequency), and the configurational entropy for the considered transition event. It can be a function of other variables, such as the number of possible nucleation sites, the temperature or an applied stress, but these dependencies are neglected here. The activated state is characterized by a unique unstable mode which has a negative eigenvalue that should be excluded from the calculation, as the system does not oscillate along this mode. There is therefore one less vibration mode at the denominator.

Evaluating the terms 𝜈 𝑖 requires the diagonalization of the Hessian matrix, i.e. the matrix of the second-order derivatives of atoms cartesian coordinates scales by the atomic masses. For each atom, the Hessian matrix has a 3 × 3 dimension, as there are three degrees of freedom for atoms movement, resulting in a matrix of dimension 3𝑛 × 3𝑛 for a system of 𝑛 atoms.

After diagonalization, the eigenmodes frequencies are obtained from the real strictly positive eigenvalues 𝜉 𝑖 as:

𝜈 𝑖 = 1 2𝜋 √ 𝜉 𝑖 , (1.8)
and the associated eigenvectors are also obtained. The diagonalization must be conducted at both the initial state and the activated state, which can be obtained by the CI-NEB or ARTn methods.

The main drawback of this approach is the high computational cost of diagonalizing the Hessian matrix that scales cubically with the number of atoms, often leading authors to neglect the effect of vibrational entropy in large systems. To reduce the size of the matrix, we performed the calculation in a cylinder around the dislocation in a PAD cell (as used by [START_REF] Barvinschi | Quantum Peierls stress of straight and kinked dislocations and effect of non-glide stresses[END_REF]), with a radius of ∼ 28 Å containing about 𝑁 𝑑 = 17 × 10 3 atoms. The negligible effect of the finite size of the cylinder was checked by performing a few full diagonalizations in collaboration with M.C. Marinica, which required considerable CPU time and advanced shared-memory linear algebra routines [START_REF] Proville | Quantum effect on thermally activated glide of dislocations[END_REF].

Free energy calculation

The use of harmonic enthalpy and entropy calculations has been ubiquitous over the last two decades, and helped solve many problems in materials science. But the harmonic assumption only holds at low temperatures [START_REF] Lesar | Finite-temperature defect properties from free-energy minimization[END_REF][START_REF] Gilbert | Free energy generalization of the Peierls potential in iron[END_REF][START_REF] Sato | Anharmonic effect on the thermally activated migration of {1012} twin interfaces in magnesium[END_REF] and a discrepancy can appear for temperatures as low as 10% of the melting temperature [START_REF] Swinburne | Uncertainty and anharmonicity in thermally activated dynamics[END_REF]. To account for anharmonic effects, several free energy calculation methods are available. Among them, collective variable-based methods are able to reduce the dimensionality of phase space [START_REF] Laio | Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science[END_REF], and have found enormous successes in studies of molecular systems. But collective variables suited for the study of dislocation glide cannot be identified [START_REF] Swinburne | Unsupervised calculation of free energy barriers in large crystalline systems[END_REF], making this method inapplicable for this problem. A path-based thermodynamic integration method from [START_REF] Cheng | Computing the absolute Gibbs free energy in atomistic simulations: Applications to defects in solids[END_REF] requires the full diagonalization of the Hessian matrix, resulting in a prohibitive computational cost for large systems. Another popular method is the finite temperature string method (E et al., 2005), which features a linear scaling of the computational time, but requires an ad-hoc parameterization [START_REF] Swinburne | Unsupervised calculation of free energy barriers in large crystalline systems[END_REF].

The projected average force integrator (PAFI) method recently introduced by Swinburne and Marinica (2018) enables free energy calculations with a linear scaling of the computation time and no assumption on the nature of thermal vibrations. In this method, a MEP is obtained by the CI-NEB method, and a series of hyperplanes perpendicular to the path are defined, as represented in Fig. 1.4. Then, constrained molecular dynamics is performed in each hyperplane, consisting in a thermalization followed by a sampling by overdamped Langevin dynamics, which allows the system to move towards a minimum free energy posi-tion. Note that the sampling can be repeated independently (i.e. with perfect parallelism) in order to localize the minimum free energy position more accurately: the number of repetitions therefore provides a simple way of adjusting the precision of PAFI calculations. Under the locality assumption -i.e. assuming the temperature pathway remains close to the MEP, a common condition in path-based methods (E et al., 2005)-the free energy gradient can be calculated [START_REF] Swinburne | Unsupervised calculation of free energy barriers in large crystalline systems[END_REF]. By integrating the free energy gradient along the path, the true free energy profile is finally obtained. Due to the sampling at finite temperature and the numerical treatment to extract results, free energy data computed using PAFI have an uncertainty4 , which increases with the temperature, as systems are then able to visit states that are far from the minimum free energy position, and which can be directly reduced by performing more sampling.

Figure 1.4: Illustration of the sampling performed by PAFI. Starting from an initial minimum energy path (MEP) 𝑋 0 (𝑟), the system is thermalized and constrained MD is performed in hyperplanes that are perpendicular to the MEP (shown in orange), resulting in a global displacement towards a minimum free energy position. An estimate of the minimum free energy position is obtained by repeating the thermalization and constrained MD several times. Following the same approach in a series of hyperplanes along the MEP, the minimum free energy path is determined. (Courtesy of TD Swinburne.)

A C++/MPI implementation of PAFI is freely available to the community and maintained by T.D. Swinburne 5 It uses LAMMPS as its backend, with the fix pafi directive 6 , and is thus fully compatible with all LAMMPS pair styles or settings, allowing e.g. Gibbs energy calculations under an external applied stress.

Gibbs energy barrier calculations under an applied stress with PAFI raise several technical issues and require a specific methodology described in the next paragraph. PAFI supports the execution of arbitrary LAMMPS code before the thermalization/sampling starts.

To apply an external stress, we use the LAMMPS directives used to perform NEB under an applied stress (see Listing 1.3) in the PreRun section of the PAFI configuration file. As multiple LAMMPS runs are performed by PAFI, the different fix, compute and region defined in PreRun must all be discarded in the PostRun section using the unfix, uncompute and region delete commands, respectively. The PAFI configuration used in our simulations is available at https://github.com/arn-all/phd-thesis.

Construction of the MEP

To start with, a candidate for the MEP between two Peierls valleys is computed using the CI-NEB method under an applied stress 𝜏 (as presented in Sec. 1.2.1), with a 1 × 10 -5 eV Å -1

force tolerance and 16 replicas (or knots). The obtained enthalpy profile is characterized by an energy difference between the initial and final states due to the plastic work: 𝜏 𝑏𝑑𝐿 (𝐿 is the dislocation length, 𝑑 the distance between Peierls valleys). The activation enthalpy Δ𝐻 is also affected by the applied stress, and decreases as the stress increases. As a result, the ratio between 𝜏 𝑏𝑑𝐿 and Δ𝐻 rapidly increases with the stress (it is close to 10 when 500 MPa are applied) and the NEB knots distribution between the initial and activated states tends to be sparse (see Fig. 1.5 (a)).

While this MEP candidate could in principle be used as a reference pathway in PAFI, a spline interpolation between knots with a large energy difference is likely to result in both poor integration accuracy, as the free energy gradient would experience large variations, and a poor estimate of the MEP tangent, which is used for the projection. This would in turn lead to erroneous projections of the reaction coordinate of thermalized systems, and in a large error in the integration of the Gibbs energy. In addition, performing the PAFI sampling in hyperplanes located after the activated state is useful to identify the maximum of the free energy profile, but it represents an additional computational cost that can be important if many knots are located after the saddle state. To increase the knots density in the zone of interest -i.e. between the initial and activated states-and remove unnecessary knots, a remeshing is performed7 , as represented in Fig. 1.5. b) is relaxed to a MEP using the "free-end NEB" method (see List. 1.4). Flat segments are then eliminated to avoid spline issues when interpolating the MEP in PAFI.

After a series of knots is constructed (Fig. 1.5 (b)), a "free-end NEB" relaxation is performed, using options of the LAMMPS fix neb command as shown in List. 1.4. The difference with a regular CI-NEB is that the last knot is not necessarily an equilibrium position. On the last knot, a spring force is applied perpendicularly to the energy path, preventing it from falling in neighbouring energy bassins. The spline interpolation of the MEP between knots that is performed by PAFI can be inaccurate on flat segments, where splined functions tend to oscillate. To avoid flat segments, we remove knots that have an absolute energy difference to the previous knot that is less than 1% of the activation enthalpy. If the knot that should be removed is the one with the maximum energy, the previous knot is removed instead. To avoid interpolation errors between knots, the PAFI Rediscretize option is disabled, such that hyperplanes used by PAFI are located at the position of the knots, instead of allowing arbitrary positions along the interpolated MEP.

The number of sampled hyperplanes is therefore the same as the number of knots along the MEP.

tions.

Methods

Following this procedure, we obtain a candidate MEP whose stability can be tested using a tool provided in PAFI (pafi-test-path), which produces a report based on a 0 K PAFI run which allows to identify flat energy segments, large energy jumps or large deviations between the MEP and the MFEP caused by interpolation issues.

Computational load and parallelism

A single PAFI run requires about 1 × 10 6 force evaluations, and supports parallelism. 2 Assessment of interatomic potentials for Fe and Fe-C

The accuracy of MD simulations is highly dependent on the quality of the interatomic potential used, and a careful choice of the potential is thus crucial.

Among the vast range of properties that can be tested, we chose to focus on the most relevant for our study, which are dislocation properties in pure iron, and dislocation-carbon interactions in dilute solutions. In iron, two main features of screw dislocations are their core structure and their Peierls barrier, both known from DFT calculations [START_REF] Ventelon | Core structure and Peierls potential of screw dislocations in 𝛼-Fe from first principles: cluster versus dipole approaches[END_REF][START_REF] Rodney | Ab initio modeling of dislocation core properties in metals and semiconductors[END_REF]. Even if it has been widely studied in the last decades, reproducing these properties with interatomic potentials while keeping an accurate description of other crystal properties remains difficult [START_REF] Barvinschi | Quantum Peierls stress of straight and kinked dislocations and effect of non-glide stresses[END_REF][START_REF] Alexander | Interatomic potentials for irradiation-induced defects in iron[END_REF][START_REF] Starikov | Angular-dependent interatomic potential for large-scale atomistic simulation of iron: Development and comprehensive comparison with existing interatomic models[END_REF]. In the next section, a litterature survey allows us to compare both properties for the most used potentials for Fe. In the opposite, dislocation properties of Fe-C potentials were often overlooked, and had to be computed.

We tested interatomic potentials for Fe-C alloys available on the NIST interatomic potentials repository [START_REF] Becker | Considerations for choosing and using force fields and interatomic potentials in materials science and engineering[END_REF][START_REF] Hale | Evaluating variability with atomistic simulations: The effect of potential and calculation methodology on the modeling of lattice and elastic constants[END_REF], listed in Table 2.1. The potentials are based on different formalisms, namely Embedded Atom Method (EAM), Modified Embedded Atom Method (MEAM) and Tersoff, all implemented in the LAMMPS MD package [START_REF] Plimpton | Fast Parallel Algorithms for Short-Range Molecular Dynamics[END_REF]. The Fe-C potentials from Becquart, Veiga and Hepburn use the Fe interaction proposed by Ackland, and therefore share the same properties in pure iron. While dislocation properties of the Ackland potential are known [START_REF] Gordon | Screw dislocation mobility in BCC metals: the role of the compact core on double-kink nucleation[END_REF], the Fe part of Liyanage, Lee and Henriksson potentials, which were respectively developed in Refs. [START_REF] Lee | Atomistic modeling of thermodynamic equilibrium and polymorphism of iron[END_REF], [START_REF] Lee | Second nearest-neighbor modified embedded atom method potentials for bcc transition metals[END_REF] and (Müller et al., 2007a) were never tested for dislocation properties. Table 2.1 also includes information about two EAM potentials for pure iron as well as ab initio DFT data, used for benchmarking.

Following our tests, the poor dislocation properties of Fe-C potentials motivated the construction of a new interatomic potential, denoted as "Hybrid", which is an hybridization of the potentials of Proville (Fe) and Becquart (Fe-C) aimed at improving the dislocation properties of the Becquart potential. Details on the construction of the Hybrid potential can be found in Sec. 2.4, while its properties are presented in comparison to other potentials throughout this chapter.

For clarity, we refer to the interatomic potentials by the name of the first author of the associated publication, e.g. "Ackland potential". A comprehensive evaluation of the performance of interatomic potentials (i.e. beyond dislocation properties) is a heavy task that is out of the scope of the present work, and rather calls for a collective effort. Several online collaborative projects [START_REF] Becker | Considerations for choosing and using force fields and interatomic potentials in materials science and engineering[END_REF][START_REF] Hale | Evaluating variability with atomistic simulations: The effect of potential and calculation methodology on the modeling of lattice and elastic constants[END_REF][START_REF] Elliott | Knowledgebase of Interatomic Models (KIM) Application Programming Interface (API)[END_REF] aim at giving a comprehensive and reproductible assessment of all available potentials, and provide readily-available properties to help potential selection.

But at the time of writing, only some classical potentials were available on these platforms; tested properties were limited to basic crystal properties such as the linear thermal expansion coefficient or the elastic constants; and the vacancy migration energy [START_REF] Li | Vacancy formation and migration energy for bcc Fe[END_REF] was the only implemented test related to defects using the NEB method.

2.1 Dislocation core structure

Pure iron

The structure of dislocation cores drives many features of dislocation mobility in bcc metals [START_REF] Rodney | Ab initio modeling of dislocation core properties in metals and semiconductors[END_REF]. Dislocation core structures can be visualized using "differential dis-2.1 Dislocation core structure placement" (DD) maps [START_REF] Vítek | The core structure of ½(111) screw dislocations in b.c.c. crystals[END_REF] which represent atomic displacements between the perfect crystal taken as a reference, (𝑠 0 ) and the relaxed dislocation (𝑠 1 ). Maps are obtained by projecting the position of atoms in cells with a thickness of 1 𝑏 along the [111] direction (i.e.

the 𝑥 axis), so each atom on the map represents an atomic column by periodicity, and their po- Using this representation, the easy and hard core structures of screw dislocations present a characteristic pattern and can be easily recognized. Burgers circuits are also visible, as they are represented by closed circuits of head-to-tails arrows.

The compact, non-degenerate easy core structure, broadly accepted as the stable configuration in iron [START_REF] Domain | Simulation of Screw Dislocation Motion in Iron by Molecular Dynamics Simulations[END_REF][START_REF] Ventelon | Core structure and Peierls potential of screw dislocations in 𝛼-Fe from first principles: cluster versus dipole approaches[END_REF][START_REF] Rodney | Ab initio modeling of dislocation core properties in metals and semiconductors[END_REF], is shown in Fig. 2.1 (a). All potentials presented in table 2.1 were tested on this property by inserting a dislocation in an easy position in iron, and minimizing the atomic forces with a 1 × 10 -6 eV Å -1 tolerance. Almost all tested potentials predict a non-degenerate, easy compact core (shown in Fig. 2.1 (a)), in agreement with ab initio results. But the potential of Lee predicts a degenerate core shown in Fig. 2.1 (b), and the Henriksson potential favours a structure that is alike the hard structure, both results being in disagreement with ab initio data and making these potentials unsuited for the simulation of screw dislocations.

Carbon-induced core reconstruction

In the presence of carbon, the structure of dislocation cores in iron can be greatly modified, according to DFT calculations [START_REF] Ventelon | Dislocation core reconstruction induced by carbon segregation in bcc iron[END_REF].

A dislocation relaxed with an array of carbon atoms placed every ⃗ 𝑏 in the octahedral position nearest to the dislocation core (denoted 𝑂 (1) and shown in Fig. 2.2) spontaneously shifts to a hard position and adopts a prismatic configuration, where the three atomic columns of the core have the same atomic coordinates projected along the [111] direction. In this configuration, which is similar to the unstable hard configuration in pure iron [START_REF] Dezerald | Modélisation ab initio des dislocations vis dans les métaux de transition cubiques centrés[END_REF], carbon atoms sit at the centre of the prisms and stabilize the hard core. The initial and force-minimized structures are illustrated in Fig. 2.3. This arrangement has similarities with the structure of cementite (Fe 3 C) and octahedral sites. As seen in Table 2.2, the first nearest neighbour (1NN) Fe-C distances in the prismatic structure are close to the 1NN distance in cementite, and to the 2NN distance in octahedral sites in bulk iron 1 . The core reconstruction is not limited to the FeC system, but was also identified in other bcc metals [START_REF] Lüthi | First principles investigation of carbon-screw dislocation interactions in body-centered cubic metals[END_REF] and with other interstitial solutes [START_REF] Lüthi | Attractive interaction between interstitial solutes and screw dislocations in bcc iron from first principles[END_REF].

For all tested potentials, we observed a conversion from an easy to a distorted hard core structure as the dislocation was relaxed with C in the 𝑂 (1) site. The observed structure is close to the ab initio structure, although distorted and asymmetrical (see Fig. 2.1 (d) for the 2.1 Dislocation core structure [START_REF] Lüthi | Modélisation ab initio des interactions entre dislocations vis et solutés interstitiels dans les métaux de transition cubiques centrés[END_REF].

O (1) O (2) O (3) O (4) O (5) O (6) O (7)
DFT calculations also evidenced a similar reconstruction upon relaxation when C atoms are inserted in octahedral sites around the core, for distances as large as 3.5 Å [START_REF] Lüthi | Modélisation ab initio des interactions entre dislocations vis et solutés interstitiels dans les métaux de transition cubiques centrés[END_REF], corresponding to the 𝑂 (7) position. In our tests with interatomic potentials, the core reconstruction was observed when C is inserted in 𝑂 (2) for the potentials of Lee and Liyanage only. In the case of Becquart, Veiga and Hybrid potential, the core is converted to a distorted degenerate core shown in Fig. 2.5. With other interatomic potentials, the presence of C in 𝑂 (2) position did not involve a change of the core structure or a movement of the dislocation.

Peierls barrier

The potential proposed by [START_REF] Mendelev | Development of new interatomic potentials appropriate for crystalline and liquid iron[END_REF] was the first empirical potential to predict a non-degenerate compact easy core [START_REF] Mrovec | Magnetic bond-order potential for iron[END_REF] 2021) also predicts a metastable intermediate core structure, and is thus not usable for the simulation of screw dislocations.

The Peierls barrier predicted by the Liyanage potential was not available in the litterature, and was thus computed. We used the LAMMPS implementation of the NEB method to determine the Peierls barrier in the ( 1 01) glide plane. NEB calculations are conducted using a force criterion for convergence of 1 × 10 -4 eV Å -1 , a spring force set to 0.1 eV Å -1 , with 48 replicas. The result is shown in Fig. 2.7. This potential, provides a poor description of the Peierls barrier, with a stable intermediate state which has a lower energy than the easy core.

A potential proposed by [START_REF] Proville | Quantum effect on thermally activated glide of dislocations[END_REF] succeeded at providing the correct shape and height of the barrier, making it suitable for the simulation of screw dislocations. A reparameterization of this potential was proposed more recently to make it suitable for simulations of C15 self-interstitial atom clusters [START_REF] Alexander | Interatomic potentials for irradiation-induced defects in iron[END_REF], but this was at the expense Mori and Ozaki (2020) (36 meV/𝑏). As an example, the Peierls barrier and corresponding dislocation trajectory computed with Mori potential is represented on Fig. 2.6. Even though they achieve a much better transferability than the Proville potential, these models are 3 to 5 orders of magnitude more computationally expensive [START_REF] Starikov | Angular-dependent interatomic potential for large-scale atomistic simulation of iron: Development and comprehensive comparison with existing interatomic models[END_REF], making their use for simulations at the scale of millions of atoms prohibitive, as seen in Fig. 2.8.

At the time of writing, Proville potential was therefore the only potential to achieve quantitative agreement with the ab initio Peierls barrier with a low computational load, and was thus adopted by other authors [START_REF] Narayanan | Crystal plasticity model for BCC iron atomistically informed by kinetics of correlated kinkpair nucleation on screw dislocation[END_REF][START_REF] Shinzato | An atomistically informed kinetic Monte Carlo model for predicting solid solution strengthening of body-centered cubic alloys[END_REF][START_REF] Ghafarollahi | Solute/screw dislocation interaction energy parameter for strengthening in bcc dilute to high entropy alloys[END_REF]. A complementary discussion on Proville potential can be found in Chap. 3.

10 0 10 1 10 2 10 3 10 4

Normalized computational time Time is normalized by 7.9 × 10 -7 s/atom/step. Data extracted from [START_REF] Starikov | Angular-dependent interatomic potential for large-scale atomistic simulation of iron: Development and comprehensive comparison with existing interatomic models[END_REF], where we added our timing for the Mori NN potential [START_REF] Mori | Neural network atomic potential to investigate the dislocation dynamics in bcc iron[END_REF], which is ∼ 10 3 times slower than Mendelev. Potentials names follow the convention adopted in Ref. [START_REF] Starikov | Angular-dependent interatomic potential for large-scale atomistic simulation of iron: Development and comprehensive comparison with existing interatomic models[END_REF].

Ac97 Ol09 Chi11 Zh04 Ch06 Men03 MCM11 Mar07 MEAM15 ADP BOP11 Mori20 GAP18

Alloy properties

The task of developing interatomic potentials for Fe-C alloys is very complex, as describing pure species accurately is already challenging. As an example, the great diversity of phases and phase transformations found in iron [START_REF] Ou | Molecular dynamics simulations of fcc-to-bcc transformation in pure iron: a review[END_REF][START_REF] Meyer | Martensite-austenite transition and phonon dispersion curves of studied by molecular-dynamics simulations[END_REF][START_REF] Müller | Analytic bond-order potential for bcc and fcc iron -Comparison with established embedded-atom method potentials[END_REF][START_REF] Song | A molecular dynamics simulation study of the velocities, mobility and activation energy of an austenite-ferrite interface in pure Fe[END_REF]Hoyt, 2012), carbon (de Tomas et al., 2019), and in the Fe-C phase diagram [START_REF] Okamoto | The C-Fe (carbon-iron) system[END_REF] are hard to represent with potentials, especially considering the difficulty to fit potentials on free energy data [START_REF] Ackland | Interatomic Potential Development[END_REF][START_REF] Ackland | Calculation of free energies from ab initio calculation[END_REF], instead of zero-temperature ab initio results. Available potentials can thus only provide a limited description of Fe-C alloys.

At zero Kelvin, most FeC potentials predict that carbon occupies octahedral sites of the bcc lattice -in agreement with ab initio calculations [START_REF] Domain | Ab initio study of foreign interstitial atom (C, N) interactions with intrinsic point defects in 𝛼 -Fe[END_REF], which predict a higher energy for tetrahedral interstitials (2.14 eV) than octahedral (1.25 eV) [START_REF] Liyanage | Structural, elastic, and thermal properties of cementite (Fe3 C) calculated using a modified embedded atom method[END_REF]. The migration energy for carbon across octahedral sites 𝐸 𝑚 is thus 0.90 eV [START_REF] Becquart | Atomistic modeling of an Fe system with a small concentration of C[END_REF]. Potentials favoring tetrahedral occupation [START_REF] Ruda | Interatomic potentials for carbon interstitials in metals and intermetallics[END_REF] [START_REF] Hepburn | Metallic-covalent interatomic potential for carbon in iron[END_REF][START_REF] Becquart | Atomistic modeling of an Fe system with a small concentration of C[END_REF][START_REF] Lee | A modified embedded-atom method interatomic potential for the Fe-C system[END_REF][START_REF] Liyanage | Structural, elastic, and thermal properties of cementite (Fe3 C) calculated using a modified embedded atom method[END_REF] compared to the DFT target value [START_REF] Becquart | Atomistic modeling of an Fe system with a small concentration of C[END_REF].

Carbon-dislocation binding energy

The binding energy of carbon atoms to dislocations is of great importance in the mechanism of solute strengthening [START_REF] Labusch | A Statistical Theory of Solid Solution Hardening[END_REF][START_REF] Maresca | Theory of screw dislocation strengthening in random BCC alloys from dilute to "High-Entropy" alloys[END_REF], and is close to -0.7 eV according to ab initio methods in the Fe-C system [START_REF] Lüthi | Modélisation ab initio des interactions entre dislocations vis et solutés interstitiels dans les métaux de transition cubiques centrés[END_REF] -denoting an highly attractive interaction-. It was calculated as the energy difference per Burgers vector between a dislocation decorated by carbon atoms (𝐸 𝑑+𝐶 ), evenly spaced at a distance 𝑑 𝑐𝑐 along the dislocation, and a dislocation in pure iron (𝐸 𝑑 ), minus the energy associated with the introduction of a carbon atom in bulk iron (𝐸 𝑏+𝐶 -𝐸 𝑏 ) [START_REF] Lüthi | Attractive interaction between interstitial solutes and screw dislocations in bcc iron from first principles[END_REF]:

𝐸 𝑏𝑖𝑛𝑑𝑖𝑛𝑔 = 1 2 [𝐸 𝑑+𝐶 (𝑑 𝑐𝑐 ) -𝐸 𝑑 ] -[𝐸 𝑏+𝐶 -𝐸 𝑏 ]. (2.1)
Here, 𝐸 𝑏+𝐶 is the energy per 𝑏 of a cell that contains one carbon atom in an octahedral position, and 𝐸 𝑏 is the energy of the same cell with no solute atom (subscript b stands for bulk).

The factor 1/2 accounts for the fact that there are two dislocations in the cell.

As seen in Fig. 2.10, the obtained binding energies are negative, meaning an attraction of the carbon atoms to the dislocations, for all potentials except Hepburn potential. The magnitude of the energy varies widely between potentials, with Henriksson and Lee potentials being closest to ab initio results [START_REF] Lüthi | Modélisation ab initio des interactions entre dislocations vis et solutés interstitiels dans les métaux de transition cubiques centrés[END_REF]. Moreover, ab initio calculations showed that the binding energy is minimum for a distance between carbon atoms in the core of 3𝑏, which represents an optimum between the gain of adding a carbon atom in the core and the cost of placing the dislocation in a hard configuration between carbon atoms. Most potentials presented in Fig. 2.10 reproduce this feature, although with a minimum at 𝑑 𝑐𝑐 = 2𝑏 in most cases. The Hybrid potential reproduces the C-dislocation binding energy quantitatively for 𝑑 𝑐𝑐 larger than 1b.

Carbon-carbon binding energy in bulk Fe

In bulk iron, the interaction energy of carbon atoms sitting in octahedral sites has been calculated by ab initio methods [START_REF] Becquart | Atomistic modeling of an Fe system with a small concentration of C[END_REF], as a function of the separation distance and crystallographic direction of the C-C pair. Here, we use the same simulation cell as in Ref. [START_REF] Becquart | Atomistic modeling of an Fe system with a small concentration of C[END_REF] in order to determine the binding energy in the [111] direction predicted by interatomic potentials. The cell is fully periodic and contains 128 atoms of bulk iron. The binding energy is calculated as:

𝐸 𝑏𝑖𝑛𝑑𝑖𝑛𝑔 (𝑑) = 𝐸 𝐶-𝐶 (𝑑) -𝐸 𝑏 -2[𝐸 𝑏+𝐶 -𝐸 𝑏 ], (2.2)
where 𝐸 𝐶-𝐶 (𝑑) is the energy of a cell that contains two carbon atoms in octahedral positions and separated by a distance 𝑑 in the [111] direction.

In DFT simulations, there is a repulsive interaction (0.09 eV) for carbon atoms separated by 1𝑏 along the [111] direction, and an attraction (-0.16 eV) for a separation of 2𝑏 [START_REF] Becquart | Atomistic modeling of an Fe system with a small concentration of C[END_REF]. The potentials presented in Fig. 2.11 show at best a partial agreement with these [START_REF] Lüthi | Modélisation ab initio des interactions entre dislocations vis et solutés interstitiels dans les métaux de transition cubiques centrés[END_REF]. Horizontal lines serve as guides for the eye.

Alloy properties

data. In particular, the second neighbor repulsion is systematically underestimated.

Regarding the hybrid potential, the interaction energies are close to Becquart and Veiga potentials at a distance of 1b, but close to zero at 2b, instead of being attractive. The carboncarbon interaction energy we obtained is thus not fully satisfactory compared to DFT. Though, it is negligible compared to dislocation core effects in this work, as we consider large C-C distances in bulk, due to the relatively low carbon concentration of the cells.

Dipolar tensor

In bulk iron, interstitial carbon atoms sitting in octahedral positions cause a local anisotropic strain due to the anisotropy of the octahedral site. Within the framework of the linear elasticity theory, this strain induces in a cell of fixed volume, 𝑉, a stress increment given by the dipolar tensor 𝑃 𝑖𝑗 [START_REF] Veiga | Comparison of atomistic and elasticity approaches for carbon diffusion near line defects in 𝛼-iron[END_REF]:

𝜎 𝑖𝑗 = - 1 𝑉 𝑃 𝑖𝑗 . (2.3)
One method to determine 𝑃 𝑖𝑗 is to measure 𝜎 𝑖𝑗 in relaxed, triperiodic, simulation cells of varying volumes, in which a carbon atom is inserted in an octahedral position. We compared the dipolar tensor obtained with the hybrid potential and with Becquart potential. The former matched experimental values better, as seen in Fig. 2.12.

Figure 2.12: Stress induced by the insertion of a carbon atom in an octahedral site as a function of the inverse volume of the simulation cell. Experimental data are from [START_REF] Schultz | Point defects in body-centered cubic transition metals[END_REF]. Symbols correspond to our MD simulations. Dashed lines are linear regressions.

Construction of an hybrid interatomic potential

Here, we showed that Proville potential features a Peierls barrier in good agreement with DFT data (see Fig. trell atmospheres formation and locking mechanisms [START_REF] Waseda | Formation of carbon Cottrell atmospheres and their effect on the stress field around an edge dislocation[END_REF][START_REF] Veiga | Monte Carlo and molecular dynamics simulations of screw dislocation locking by Cottrell atmospheres in low carbon Fe-C alloys[END_REF][START_REF] Tchitchekova | A novel method for calculating the energy barriers for carbon diffusion in ferrite under heterogeneous stress[END_REF][START_REF] Candela | Investigating the kinetics of the formation of a C Cottrell atmosphere around a screw dislocation in bcc iron: a mixed-lattice atomistic kinetic Monte-Carlo analysis[END_REF]. Both Proville and Veiga potentials are EAM potentials derived from the seminal Mendelev potential [START_REF] Mendelev | Development of new interatomic potentials appropriate for crystalline and liquid iron[END_REF]. We thus propose here an hybridization of these potentials to combine their advantages, while keeping the properties of Proville potential in pure Fe.

In EAM potentials, the energy 𝐸 𝑖 of an atom 𝑖, of species 𝛼, is the sum of two terms within a cutoff distance: a pairwise interaction 𝜙 𝛼𝛽 summed over all neighbors 𝑗 (species 𝛽); and the function 𝐹 𝛼 that represents the embedding of atom 𝑖 in the electron cloud formed by its neighbors. It is written as:

𝐸 𝑖 = 𝐹 𝛼 (∑ 𝑖≠𝑗 𝜌 𝛽 (𝑟 𝑖𝑗 )) + 1 2 ∑ 𝑖≠𝑗 𝜙 𝛼𝛽 (𝑟 𝑖𝑗 ), (2.4)
where 𝑟 𝑖𝑗 is the distance between atoms 𝑖 and 𝑗, and 𝜌 𝛽 is the electron density function.

In order to combine Proville and Veiga potentials, we used a classical Gauge transformation [START_REF] Ackland | Interatomic Potential Development[END_REF], which consists in rescaling the electron density and embedding function of iron in a way that does not affect the Fe-Fe interaction. In iron, we define the electron density function ρ 𝐹 𝑒 and the embedding function F 𝐹 𝑒 (𝜌) of the new hybrid potential as follows, based on the corresponding functions of the Proville potential which are identified by a ⋆ symbol:

ρ 𝐹 𝑒 = 𝐴 × 𝜌 ⋆ 𝐹 𝑒
(2.5)

F 𝐹 𝑒 (𝜌 𝐹 𝑒 ) = 𝐹 ⋆ 𝐹 𝑒 ( 𝜌 𝐹 𝑒 𝐴 ) , (2.6)
where 𝐴 is a constant. The function 𝜙 𝐹 𝑒-𝐹 𝑒 (𝑟) was taken from Proville potential [START_REF] Proville | Quantum effect on thermally activated glide of dislocations[END_REF] with no modification. As a consequence, the interaction between iron atoms within the cutoff radius is the same as with the Proville potential:

F 𝐹 𝑒 ( ρ 𝐹 𝑒 ) = 𝐹 ⋆ 𝐹 𝑒 (𝜌 ⋆ 𝐹 𝑒 ) (2.7)
This way, we emphasize that the hybrid potential has the same materials properties in pure iron as Proville potential. However, it must be noted that in the case of alloys, the electron density is modified, as all neighbors contribute to the term ∑ 𝑖≠𝑗 𝜌 𝛽 (𝑟 𝑖𝑗 ), regardless of their species. All functions involving carbon, namely 𝜌 𝐶 (𝑟), 𝐹 𝐶 (𝜌), 𝜙 𝐹 𝑒-𝐶 (𝑟) and 𝜙 𝐶-𝐶 (𝑟) were taken from the Veiga potential [START_REF] Veiga | Comments on "atomistic modeling of an Fe system with a small concentration of C[END_REF] with no modification.

The introduction of parameter 𝐴 in Eq. 2.5 provides a degree of freedom, which was used to tune the hybrid potential in order to obtain the best agreement with the carbon-dislocation interaction energies computed with DFT. The result is shown in Fig. 2.10 and was obtained with 𝐴 = 0.55.

In the original tabulated Proville potential, the electronic density function 𝜌 ⋆ becomes abruptly constant for 𝑟 <2 Å, resulting in a discontinuity in the forces around 2.0 Å. This discontinuity is confirmed when computing the energy of a dimmer of Fe atoms and its derivatives at a varying distance with LAMMPS. This result is available publicly in the online OpenKIM archive [START_REF] Tadmor | Verification Check of Dimer C1 Continuity v004[END_REF][START_REF] Elliott | Knowledgebase of Interatomic Models (KIM) Application Programming Interface (API)[END_REF]. In bcc iron, the discontinuity is usually harmless as the shortest equilibrium interatomic distance, which corresponds to 𝑏 = 1/2[111], is close to 2.4 Å. But in the Fe-C system, octahedral sites have 1NN distances below 2 Å (see Table 2.2), possibly causing large errors in force integration and a shift of the total energy. As a simple way to correct for this discontinuity, we replace the original tabulated values between 1.95 and 2.05 Å by a cubic spline using the UnivariateSpline method of the SciPy 1.5.2 package.

The Hybrid potential is used throughout the rest of this work, with its parameter 𝐴 set to 0.55. In additional tests, we obtained a migration barrier for C interstitials of 0.96 eV, which is close to the ab initio value of 0.90 eV [START_REF] Becquart | Atomistic modeling of an Fe system with a small concentration of C[END_REF], but with a rough energy landscape around the migration path between octahedral sites. We thus do not recommend using this potential for studies on diffusion mechanisms, or with large C concentrations where solutes occupy adjacent lattice sites. Following these results, it appears that Proville potential is the only potential to enable large scale simulations of screw dislocations in iron in agreement with the Peierls barrier obtained with DFT, and with an acceptable computational load. But as any simplified interatomic model, it also has several downsides [START_REF] Starikov | Angular-dependent interatomic potential for large-scale atomistic simulation of iron: Development and comprehensive comparison with existing interatomic models[END_REF] and does not describe "real-Fe": it poorly reproduces the phonon spectra for bcc Fe, the formation energy of selfinterstitial atoms, or the energy of symmetric grain boundaries, and the kink pair nucleation enthalpy is underestimated (this value is compared to other methods in 3). As new fitting strategies based on machine learning give promising results [START_REF] Mishin | Machine-learning interatomic potentials for materials science[END_REF], novel potentials for Fe that outperform Proville, especially in terms of transferability, are expected to emerge over the coming years.

Conclusion and recommendations

For the Fe-C system, there was no available potential with a single-humped Peierls barrier, leading us to recommend our hybridization of Proville and Veiga potential. While it also has drawbacks and is likely to find a limited use, it unlocks a major limitation of current potentials, allowing for the first time the simulation of carbon interacting with screw dislocations that move with the correct glide mechanism.

3 Thermally-activated glide of screw dislocations in pure bcc iron Low-temperature plasticity of bcc metals is controlled by the mobility of screw dislocations [START_REF] Caillard | Thermally Activated Mechanisms in Crystal Plasticity[END_REF][START_REF] Kubin | Dislocations, Mesoscale Simulations and Plastic Flow[END_REF][START_REF] Proville | Modeling the thermally activated mobility of dislocations at the atomic scale[END_REF], which is limited by both the strong resistance of the bcc atomic lattice, and interactions with other defects [START_REF] Bacon | Chapter 88 Dislocation-Obstacle Interactions at the Atomic Level[END_REF], among which are dislocations, grain boundaries, precipitates or solute atoms. In this chapter, the effect of the bcc lattice itself on dislocation motion is studied, while interactions between dislocations and interstitial carbon atoms will be discussed in Chap. 4.

Lattice resistance in bcc metals

In the bcc lattice projected in the (111) plane, the centers of the triangles formed by [111] atomic rows are high-symmetry positions for1 2 ⟨111⟩ screw dislocations [START_REF] Kubin | Dislocations, Mesoscale Simulations and Plastic Flow[END_REF]. The center of upward triangles (in the convention adopted in Fig. 3.1 (b)) is a position of minimum energy, denoted as easy, where the helicity is reversed compared to the perfect crystal. In downward triangles, the helicity of the screw dislocation cancels with the helicity of the crystal, and atoms of the core have the same coordinate along the [111] direction [START_REF] Dezerald | Plastic anisotropy and dislocation trajectory in BCC metals[END_REF]. This position is an energy maximum according to ab initio calculations [START_REF] Ventelon | Core structure and Peierls potential of screw dislocations in 𝛼-Fe from first principles: cluster versus dipole approaches[END_REF]. Another high-symmetry position, called split, is obtained when the dislocation is located in the close vicinity of one [111] atomic row 1 . This position has a high energy according to DFT results, but it is a common artefact of interatomic potentials to feature a metastable split structure (see Fig. 2.6).

To move from one easy position to the next, the dislocation has to pass through a highenergy state located in the dividing plane defined by the hard and split positions (see Fig. where sits the split core, and antitwinning (AT) region where sits the hard core, illustrated in green and pink in Fig. 3.1 [START_REF] Dezerald | Ab initio modeling of the two-dimensional energy landscape of screw dislocations in bcc transition metals[END_REF]. While this deviation is low in the case of iron [START_REF] Dezerald | Plastic anisotropy and dislocation trajectory in BCC metals[END_REF]) (𝛼 = -1.7°in the absence of an applied stress), another common artefact of interatomic potentials for iron is to predict a large deviation (see Fig. 3.1 Lattice resistance in bcc metals the barrier 2 . At low temperature and applied stress, the kink pair mechanism is the preferred mechanism for dislocation glide [START_REF] Kubin | Dislocations, Mesoscale Simulations and Plastic Flow[END_REF]. Conversely, for levels of stress that are above the Peierls stress, the Peierls barrier vanishes (i.e. no thermal activation is required for dislocations to move) and dislocations move while remaining straight. Under an applied shear stress, a plastic work 𝜏 𝑏𝑥𝐿 is dissipated when a dislocation of length 𝐿 moves by a distance 𝑥 in the glide direction, and we note the enthalpy for kink pair formation Δ𝐻 𝑘𝑝 . In Fig. 3.3 (a), two examples of minimum enthalpy paths are represented, one with a low applied stress, and the other with a higher stress. In Fig. 3.3 (b), the stress-dependence of Δ𝐻 𝑘𝑝 for Proville potential is presented up to 𝜏 𝑃 . With increasing levels of stress, the difference between the equilibrium and saddle states for kink pair formation decreases and kinks become smaller, as shown in Fig. 3.4.

It is common to fit the evolution of Δ𝐻 (𝜏 ) by a Kocks law [START_REF] Kocks | Thermodynamics and kinetics of slip[END_REF] of the form:

Δ𝐻 (𝜏 ) = Δ𝐻 0 (1 -( 𝜏 𝜏 𝑃 ) 𝑝 ) 𝑞 , (3.1)
where 0 < 𝑝 ≤ 1 and 1 ≤ 𝑞 ≤ 2 are fitting constants and Δ𝐻 0 is the zero-stress barrier for kink pair nucleation. The low stress limit of the kink pair nucleation energy Δ𝐻 (𝜏 = 0), also 2 In pure bcc metals, the secondary Peierls barrier, which is the energetic cost of propagating the kinks along the line, is negligible -while in alloys, kinks can be pinned by solute atoms (as described in Chap. 4). Kink-pair nucleation and subsequent propagation in pure metals is thus often simply designed as kink pair nucleation or formation.

noted 𝐸 𝑘𝑝 is ill-defined as with no applied stress, the kinks are only stabilized by their interaction with periodic images: the activated kink pair is therefore symmetrical, with a width that is one half of the simulation cell dimension 𝐿 in

Glide of screw dislocations

Atomistic simulations

We simulated the glide of a screw dislocation in a PAD simulation cell (see Fig. 1.2) under an imposed shear rate. Examples of stress strain curves obtained at different temperatures with the Proville potential are shown in Fig. 3.5. In dynamical simulations, the stress can be computed using the relation:

𝜏 𝑥𝑧 = 𝐹 + 𝑥 -𝐹 - 𝑥 2𝐴 (3.2)
where 𝐹 ± 𝑥 is the total force in the upper or lower slab (+ andrespectively), projected in the ⃗ 𝑥 direction (i.e. the glide direction) and A is the area of the free surfaces of the crystal, i.e. 𝐴 = 𝐿 𝑥 × 𝐿 𝑦 . The stress computed using Eq. 3.2 undergoes large oscillations (Rodney, 2007), which are eliminated using a Savitzky-Golay filter3 of order 1 [START_REF] Savitzky | Smoothing and Differentiation of Data by Simplified Least Squares Procedures[END_REF], with a window of 51 points (the time difference between successive points is 0.5 ps).

The shear stress can also be calculated as the internal stress of the cell 𝜏 𝑥𝑧 (named 𝑝𝑥𝑧 in LAMMPS), multiplied by the fraction of the cell volume occupied by the crystal, resulting in less oscillations. We verified that depending on the method used to calculate the shear stress and the chosen average window length, variations of the average stress remained below 20 MPa for the typical stress levels of our MD simulations.

The flow stress is then extracted as the average of the successive peaks of stress in the plastic regime4 (marked by round symbols). It can be seen that the flow stress decreases with the temperature, as glide is a thermally activated process, and that the peaks height can strongly fluctuate, as dislocation jump is a stochastic process. Five independent MD simulations with different initial atomic velocities were performed to obtain an accurate measure of the flow stress.

Due to the computational load associated with the resolution of classical dynamics equations, our MD simulations are limited to ∼ 10 10 -10 13 steps × atoms, corresponding e.g. to the ∼ 10 5 atoms needed to model dislocation glide, simulated over ∼ 10 6 time steps of 1 fs, for an overall simulated time of the order of 1 ns. In order to reach stress levels that are large enough to observe plastic events in such limited simulated time, high strain rate and/or temperature values are necessary. At low temperatures, typical MD simulations are therefore conducted at strain rates ∼ 10 10 times larger than experimental conditions (see e.g. Proville and Rodney (2020) for a review on atomistic simulations of dislocations), leading to much higher flow stress values.

Kinetic models

Kinetic models based on the Transition State Theory (see Chap.1.2) can help compare MD simulations conducted at high shear rates to experimental conditions. In this section, we consider different models of dislocation glide based on the TST, which are then discussed in comparison to MD simulations. where Δ𝐺 𝑘𝑝 is the Gibbs free energy of nucleation for a kink-pair, 𝜈 is an attempt frequency assumed constant, 𝛽 = 1/(𝑘𝑇 ), and 𝐿 is the dislocation segment length available for kink pair nucleation. Although non-Schmid effects would be important to include in a comprehensive study [START_REF] Christian | Some surprising features of the plastic deformation of body-centered cubic metals and alloys[END_REF][START_REF] Dezerald | Plastic anisotropy and dislocation trajectory in BCC metals[END_REF][START_REF] Kraych | Nonglide effects and dislocation core fields in BCC metals[END_REF], we here only consider the effect of the resolved shear stress 𝜏 in the {110} glide plane of the dislocation, instead of the whole stress tensor [𝜎].

Computing Δ𝐺 𝑘𝑝 remains notoriously difficult [START_REF] Swinburne | Uncertainty and anharmonicity in thermally activated dynamics[END_REF], especially given the dimensions of the system that need to be considered to model the nucleation of kink-pairs, which is beyond the limits of conventional methods (see Sec. 1.3.1). Several assumptions are classically used to simplify the calculation of Δ𝐺 𝑘𝑝 , which will be discussed in regard to full Gibbs energy calculations allowed by the recent PAFI method and presented in this chapter.

The connection between the dislocation velocity given in Eq. 3.3 and the flow stress computed in strain-controlled simulations is obtained using Orowan's equation in the steadystate regime of deformation:

𝓋 = γ 𝜌𝑏 (3.4)
where 𝜌 = 1/(𝐿 𝑦 𝐿 𝑧 ) is the dislocation density.

Neglecting stress fluctuations between jump events and equating Eq. 3.3 and 3.4, the flow stress 𝜏 * is solution of:

γ 𝜌𝑏 = 𝜈𝐿𝑒 -𝛽Δ𝐺 𝑘𝑝 (𝜏 * ,𝑇 ) , (3.5)
which can be written:

Δ𝐺 𝑘𝑝 (𝜏 * , 𝑇 ) = 𝑘𝑇 log ( 𝜌𝑏𝜈𝐿 γ ) . (3.6)
For given simulation conditions, the dimensionless ratio Γ = 𝜌𝑏𝜈𝐿 γ can be considered fixed as well as 𝑇, and the value of 𝜏 * can be determined numerically if Δ𝐺 𝑘𝑝 is a known function of 𝜏 and 𝑇. Within this approach, it is customary to adjust the term 𝜈 to best match the MD data [START_REF] Rodney | Activation enthalpy for kink-pair nucleation on dislocations: Comparison between static and dynamic atomic-scale simulations[END_REF], which is usually obtained with a frequency on the order of the Debye frequency, i.e. 𝜈 ∼ 1 × 10 13 s -1 .

While it gives a simple way of modelling dislocation glide knowing only Δ𝐺 𝑘𝑝 (or an approximation), this model has issues at low stress levels, as it only considers the forward motion of the dislocation, and assumes the kink pair migration time negligible compared to nucleation time. In the limit of a low applied stress, backward jumps become more likely and must be accounted in the average velocity. In addition, kink propagation along the line is driven by the applied stress, and becomes limited under a low stress. For both reasons, the dislocation velocity should tend to zero for low stress levels, while in the above model, 𝓋 = 𝜈𝐿 exp(-𝛽Δ𝐺(0, 𝑇 )) ≠ 0 at zero stress, which is unphysical.

Thermally-activated forward and backward-jump model

To account for the possibility of reverse-jumps at low stress [START_REF] Nabarro | One-dimensional models of thermal activation under shear stress[END_REF], an additional term can be added to Eq. 3.3. To move backwards, the reverse energy barrier is Δ𝐺 𝑘𝑝 plus an additional energy cost due to the plastic work between the initial and final states 𝜏 𝑏𝑑𝐿, which gives:

𝓋 = 𝜈𝐿 [𝑒 -𝛽Δ𝐺 𝑘𝑝 -𝑒 -𝛽 (Δ𝐺 𝑘𝑝 +𝜏 𝑏𝑑𝐿) ] (3.7)
When 𝜏 𝑏𝑑𝐿 ≫ Δ𝐺 𝑘𝑝 , this expression is equivalent to Eq. 3.3. One consequence is that reverse jump can remain unlikely down to low stress levels in the case of long segments, thus not effectively reducing the dislocation velocity. Values of Δ𝐺 𝑘𝑝 will be given later in this chapter (e.g. in Fig. 3.12), but it can already be noted that its low-stress, low-temperature values are close to Δ𝐻 (𝜏 = 0) = 0.54 eV (according to Fig. 3.3), and decrease with T and 𝜏. In Writing 𝑥 = 𝛽𝜏 𝑏𝑑𝐿/2, Eq. 3.7 is equivalent to:

𝓋 = 𝜈𝐿𝑒 -𝛽Δ𝐺 𝑘𝑝 𝑒 -𝑥 (𝑒 𝑥 -𝑒 -𝑥 ) (3.8)
Using 2 sinh(𝑥) = 𝑒 𝑥 -𝑒 -𝑥 , we obtain the following formulation:

𝓋 = 2𝜈𝐿 sinh ( 𝛽𝜏 𝑏𝑑𝐿 2 ) 𝑒 -𝛽(Δ𝐺 𝑘𝑝 +𝜏 𝑏𝑑𝐿/2) (3.9)
Because of the sinh term, this model allows the velocity to tend to zero when 𝜏 tends to zero, as forward and backward jumps are equally probable under zero stress. Here, the flow stress 𝜏 * verifies:

Δ𝐺(𝜏 * , 𝑇 ) = 𝑘𝑇 log ( 2𝜌𝑏𝜈𝐿 γ sinh ( 𝛽𝜏 * 𝑏𝑑𝐿 2 ) 𝑒 -𝛽𝜏 * 𝑏𝑑𝐿/2 ) , (3.10)
which can be solved numerically. The 𝑒 𝑥 term in sinh can cause numerical overflow, and better stability is therefore achieved by rewriting it as:

Δ𝐺(𝜏 * , 𝑇 ) = 𝑘𝑇 log ( 𝜌𝑏𝜈𝐿 γ (1 -𝑒 -𝛽𝜏 * 𝑏𝑑𝐿 )) , (3.11)
to obtain values of 𝜏 * knowing Δ𝐺(𝜏 * , 𝑇 ).

Dislocation glide model by Po et al. (2016)

The low-stress domain can also be understood as a low kink migration velocity 𝓋 𝑘 (𝜏 , 𝑇 ) domain, where kink pair expansion is the limiting feature, even if there are many kink nucleation events along the line e.g. at high temperature or for long segments. Since lattice resistance on kinks is negligible in pure bcc metals, kink velocity is limited by phonon drag and can be written as:

𝓋 𝑘 = 𝜏 𝑏 𝐵(𝑇 ) (3.12)
where the effective drag coefficient is 𝐵(𝑇 ) = 𝐵 0 + 𝐵 1 (𝑇 ). 𝐵(𝑇 ) depends on the temperature in the general case, but is approximately constant in the case of screw dislocations, i.e. 𝐵 1 = 0 [START_REF] Po | A phenomenological dislocation mobility law for bcc metals[END_REF]. The value 𝐵 0 = 9.8 × 10 -4 Pa s was determined for screw dislocations in tungsten based on the work of [START_REF] Stukowski | Thermally-activated non-Schmid glide of screw dislocations in W using atomistically-informed kinetic Monte Carlo simulations[END_REF] and the Frenkel-Kontorova model from [START_REF] Swinburne | Theory and simulation of the diffusion of kinks on dislocations in bcc metals[END_REF]. In this work, we use the same value as a rough estimate for our calculations in iron with satisfactory results, even though a precise determination of 𝐵 0 would be valuable. Reducing the value of 𝐵 0 would result in a higher kink velocity, and a softening at high temperatures where kink migration is dominant (see below), with no significant change of the model predictions at low temperature. The dislocation velocity is expressed as:

𝓋 = ℎ𝐽 𝑘𝑝 𝑋 (3.13)
where 𝐽 𝑘𝑝 is the kink pair nucleation rate obtained from the classical kink-diffusion theory of [START_REF] Hirth | Theory of dislocations[END_REF] and is proportional to 𝓋 𝑘 , and ℎ is the kink height. 𝑋 is the length available for kink pair nucleation, which is the total segment length 𝐿 if it can only support one kink pair (which is typical of MD simulations). For longer segments, it corresponds to the average distance between existing kink-pairs 𝜆, which is obtained from the equilibrium concentration of kink-pairs [START_REF] Hirth | Theory of dislocations[END_REF], as a function of Δ𝐺 𝑘𝑝 :

𝜆 = 2𝑎𝑒 𝛽Δ𝐺 𝑘𝑝 /2 (3.14)
where 𝑎 is the crystal lattice parameter. An effective measure of 𝑋 is therefore obtained by the geometrical average 𝑋 = 𝜆𝐿/(𝜆 + 𝐿).

The main result of the velocity law proposed by Po et al. is an expression of dislocation velocity that accounts for both the low and high stress regime, written:

𝓋 = 𝓋 𝑘 2ℎ 𝑎 𝐿 2𝑎𝑒 𝛽Δ𝐺 𝑘𝑝 /2 + 𝐿 𝑒 -𝛽Δ𝐺 𝑘𝑝 /2 . (3.15)
Note that 𝓋 is proportional to 𝜏 through the term 𝓋 𝑘 , giving a null velocity when the applied resolved shear stress is zero. Isolating Δ𝐺 𝑘𝑝 and introducing Orowan's equation as above, we obtain:

Δ𝐺 𝑘𝑝 (𝜏 * , 𝑇 ) = 2𝑘𝑇 log ( 1 4𝑎 ( √ 𝐿 2 + 16ℎ𝐿𝜌𝑏 2 𝜏 * γ 𝐵(𝑇 ) -𝐿)) .
(3.16)

The models presented above all require the knowledge of Δ𝐺 𝑘𝑝 , which is usually approximated using hypothesis presented below. A comparison of the models introduced in this section and a discussion of their limits will then be presented in § 3.2.4.

Classical assumptions in Harmonic TST

While the harmonic assumption on atomic vibrations considerably simplifies the evaluation of the Gibbs energy [START_REF] Vineyard | Frequency factors and isotope effects in solid state rate processes[END_REF] -which can be decomposed in an enthalpic and entropic contribution that are independent of the temperature-, further approximations are often necessary to evaluate the entropic term. For large systems, harmonic entropy calculations are indeed limited by both the computational time and shared-memory cost associated with the diagonalization of the full Hessian matrix [START_REF] Proville | Quantum effect on thermally activated glide of dislocations[END_REF]. [START_REF] Gilbert | Stress and temperature dependence of screw dislocation mobility in 𝛼-Fe by molecular dynamics[END_REF] proposed to fit Δ𝐻 𝑘𝑝 based on MD simulations of dislocation glide, while neglecting the effect of entropy. While it results in accurate velocity predictions, such approach uses a large number of numerical fitting parameters and does not provide an understanding of the physical mechanisms.

Constant activation entropy

The simplest approach is to assimilate the activation entropy as a constant noted Δ𝑆 0 , yielding Δ𝐺 = Δ𝐻 -𝑇 Δ𝑆 0 . In the case of a thermally activated velocity model (Eq. 3.3), the constant Δ𝑆 0 can be included in the prefactor as: (3.17) where the constant 𝜈 * is an effective attempt frequency which can be fitted to best match the MD data. However, as the model of Po et al. does not have an attempt frequency term (see Eq. 3.15), we set 𝜈 = 2 × 10 13 s -1 and only change the value of Δ𝑆 0 for clarity in the comparison between models.

𝓋 = 𝜈𝐿𝑒 -𝛽(Δ𝐻 (𝜏 )-𝑇 Δ𝑆 0 ) = 𝜈𝐿𝑒 Δ𝑆 0 /𝑘 𝑒 -𝛽Δ𝐻 (𝜏 ) = 𝜈 * 𝐿𝑒 -𝛽Δ𝐻 (𝜏 ) ,
If this approximation has been validated for a specific dislocation by comparison with MD simulations of the glide of Lomer dislocations in Al [START_REF] Rodney | Activation enthalpy for kink-pair nucleation on dislocations: Comparison between static and dynamic atomic-scale simulations[END_REF], there are evidence of temperature dependence of energy barriers for other dislocation processes, including dislocation nucleation [START_REF] Warner | Origins and implications of temperature-dependent activation energy barriers for dislocation nucleation in face-centered cubic metals[END_REF][START_REF] Ryu | Entropic effect on the rate of dislocation nucleation[END_REF], obstacle by-pass [START_REF] Saroukhani | Harnessing atomistic simulations to predict the rate at which dislocations overcome obstacles[END_REF] and cross-slip [START_REF] Santos-Güemes | Influence of the stress state on the cross-slip free energy barrier in Al: An atomistic investigation[END_REF]. Moreover, the Gibbs free energy of a straight screw dislocation in BCC iron has also been shown strongly dependent on temperature [START_REF] Gilbert | Free energy generalization of the Peierls potential in iron[END_REF].

Meyer-Neldel law

Another classical approximation to include the evolution of the entropy without explicitly computing it is to invoke the Meyer-Neldel (MN) compensation rule [START_REF] Yelon | Multi-excitation entropy: its role in thermodynamics and kinetics[END_REF], which states that the activation entropy of a process is proportional to its activation enthalpy: (3.18) with 𝑇 𝑀𝑁 being a characteristic temperature and Δ𝑆 0 a constant. If 𝑇 𝑀𝑁 is positive, this means that an increase of Δ𝐻 leads to an increase of Δ𝑆, leading to a partial compensation in the Gibbs energy Δ𝐺 = Δ𝐻 -𝑇 Δ𝑆. The Gibbs energy can then be written:

Δ𝑆 = Δ𝐻 /𝑇 𝑀𝑁 + Δ𝑆 0 ,
Δ𝐺 𝑘𝑝 (𝜏 * , 𝑇 ) = Δ𝐻 (1 - 𝑇 𝑇 𝑀𝑁 ) + Δ𝑆 0 . (3.19)
Similarily to Eq. 3.17, Δ𝑆 0 can be recasted in the 𝜈 * term for thermally activated velocity laws.

Usually observed when varying a material parameter, as in the diffusion of different solutes in a given host matrix [START_REF] Gelin | Enthalpy-entropy compensation of atomic diffusion originates from softening of low frequency phonons[END_REF] or magnetization reversals in memory elements of different composition [START_REF] Desplat | Entropy-reduced retention times in magnetic memory elements: A case of the Meyer-Neldel Compensation Rule[END_REF], the compensation law has also been reported

for stress-driven processes in a given material, i.e. :

Δ𝑆([𝜎]) ∝ Δ𝐻([𝜎]) 𝑇 𝑀𝑁 , (3.20)
for instance for dislocation nucleation [START_REF] Hara | Adaptive strain-boost hyperdynamics simulations of stress-driven atomic processes[END_REF], creep [START_REF] Wang | Entropic effect on creep in nanocrystalline metals[END_REF]), cross-slip (Esteban-Manzanares et al., 2020) and twin motion [START_REF] Sato | Anharmonic effect on the thermally activated migration of {1012} twin interfaces in magnesium[END_REF] -and could similarily be used for kink pair nucleation.

Applications of the MN law in TST models include dislocation nucleation [START_REF] Zhu | Temperature and strain-rate dependence of surface dislocation nucleation[END_REF][START_REF] Amodeo | Site dependence of surface dislocation nucleation in ceramic nanoparticles[END_REF] and glide [START_REF] Clouet | Screw dislocations in BCC transition metals: From ab initio modeling to yield criterion[END_REF]. [START_REF] Gelin | Enthalpy-entropy compensation of atomic diffusion originates from softening of low frequency phonons[END_REF] recently tracked down the origin of the compensation law in the harmonic regime to a systematic softening of lowfrequency vibrational modes between the initial and activated states.

Comparison and discussion

In this section, we compare the flow stress predicted by each model introduced in Sec. 3.2.2, using values of Δ𝐺 determined using different classical assumptions, with flow stress values obtained from strain-rate controlled MD simulations. In MD simulations, 𝜌 = 4 × 10 15 m -2 , 𝐿 = 40𝑏, we take 𝜈 = 2 × 10 13 s -1 and γ is set to 1 × 10 6 s -1 or 1 × 10 7 s -1 , yielding respectively model, the flow stress decreases until reaching a critical temperature 𝑇 𝐶 which can be seen as the limit of the thermally-activated regime, before glide becomes athermal. 𝑇 𝐶 satisfies Eq. 3.6 :

Γ = 𝜌𝑏𝜈𝐿 γ = 2 ×
𝑘𝑇 𝐶 log (Γ) = Δ𝐺(𝜏 * = 0, 𝑇 𝐶 ), (3.21)
and corresponds to a null flow stress. We have here a well-defined athermal temperature, but plastic flow occurs due to the imposed shear rate (with a finite dislocation velocity 𝓋 = γ /(𝜌𝑏)), while the stress remains zero on average. This behavior is unphysical and is an artefact of the forward jump model, as dislocation glide by the kink pair mechanism can only occur when a non-zero stress 𝜏 allows the expansion of kinks. This issue is addressed by including the possibilty of backward jumps under a low stress. As a result, the athermal temperature 𝑇 𝐶 is no longer rigorously defined. However, we see in that enforces 𝓋 = 0 if 𝜏 = 0. With this model also, the limit of a null flow stress is never reached, and the stress decreases asymptotically at high temperature. At high temperature, the decrease of the flow stress is slower than with the thermally activated model, making it more difficuly to extrapolate an athermal temperature from the low flow stress domain.

From the expression of the velocity in the model of Po et al. (Eq. 3.15), different regimes emerge, where the velocity is limited either by kink migration or kink pair nucleation. For long segments, where nucleation is readily achieved thanks to the great amount of available nucleation sites, or at high temperature where Δ𝐺 𝑘𝑝 is lowered, the average distance between kinks pairs is 𝜆. Kink pair migration then consists in kinks crossing the distance 𝜆 to annihilate on one another, with an average migration time 𝑡 𝑚𝑖𝑔 = 𝜆/𝓋 𝓀 which does not depend on 𝐿. In this domain, the flow stress is limited by kink migration, and is higher compared to the thermally activated model -which only accounts for nucleation-as seen in Fig. 3.7. On the other hand, if the segment length 𝐿 ≪ 𝜆, e.g. in the case of a simulation cell with a short dislocation segment, kink migration will be faster, as it will stop after crossing the distance 𝐿.

In addition, short segments provide only a low number of nucleation sites, which limits the kink pair nucleation rate. In this regime, kink nucleation becomes the limiting feature and 𝓋 depends on 𝐿 as both the nucleation time and the migration time 𝑡 𝑚𝑖𝑔 = 𝐿/𝓋 𝓀 depend on L.

These two regimes can be identified in the expression of the dislocation velocity (Eq. 3.15).

Using the approximation 𝐿 ≫ 2𝑎𝑒 𝛽Δ𝐺 𝑘𝑝 /2 it is simplified as:

𝓋 = 𝓋 𝑘 2ℎ 𝑎 𝐿 2𝑎𝑒 𝛽Δ𝐺 𝑘𝑝 /2 + 𝐿 𝑒 -𝛽Δ𝐺 𝑘𝑝 /2 ∼ 𝓋 𝑘 2ℎ 𝑎 𝑒 -𝛽Δ𝐺 𝑘𝑝 /2 , (3.22)
which is dependent on 𝐿 and corresponds to the regime where kink migration is limiting. If the contrary approximation 𝐿 ≪ 2𝑎𝑒 𝛽Δ𝐺 𝑘𝑝 /2 is used to simplify Eq. 3.15, we obtain:

𝓋 = 𝓋 𝑘 2ℎ 𝑎 𝐿 2𝑎𝑒 𝛽Δ𝐺 𝑘𝑝 /2 + 𝐿 𝑒 -𝛽Δ𝐺 𝑘𝑝 /2 ∼ 𝓋 = 𝓋 𝑘 2ℎ 𝑎 𝐿𝑒 -𝛽Δ𝐺 𝑘𝑝 . (3.23)
which is close to the velocity of the thermally-activated model (Eq. 3.3).

The ratio 𝐿/𝜆 which determines the limiting feature for dislocation velocity is plotted against Δ𝐺 𝑘𝑝 in Fig. 3 , 𝐿 = 40𝑏, γ = 1 × 10 7 s -1 ), while red dots correspond to typical experimental conditions (𝜌 = 1 × 10 12 m -2 , 𝐿 = 1 × 10 -6 m, γ = 8.5 × 10 -4 s -1 ).

ments. For the thermally activated model in the other hand, the flow stress slowly tends to zero as longer segments increase the nucleation probability without augmenting the kink migration time, which is neglected.

We emphasise that the kink migration or nucleation-limited regimes are not solely determined by the value of 𝐿, as 𝜆 is a function of Δ𝐺(𝜏 , 𝑇 ). In Fig. 3.8, it can be seen that the motion of a long segment of 1 µm can be limited by kink nucleation if Δ𝐺 is large( i.e. in the case of a low shear rate) or if the temperature is low5 .

Gibbs energy approximations In Fig. 3.9, we compare our MD data to the flow stress obtained based on both harmonic assumptions of Δ𝐺 𝑘𝑝 presented in Sec. We also tested a Meyer Neldel law to compute Δ𝐺 𝑘𝑝 , which is presented in the second row of Fig. 3.9 with 𝑇 𝑀𝑁 = 500 K. We can see that the agreement with MD data is made worse by introducing the MN law compared to the constant entropy assumption, as it makes the decrease of the flow stress much stiffer. The poor agreement strongly suggests that this system doesn't obey a Meyer-Neldel compensation rule.

These models are therefore not sufficient to explain the temperature evolution of the flow stress observed in the MD simulations. In addition, they are based on unchecked assumptions on Δ𝑆 whose effect is difficult to decorelate from the effect of the underlying velocity model.

The temperature and stress dependence of Δ𝐺 𝑘𝑝 has been estimated by [START_REF] Gilbert | Free energy generalization of the Peierls potential in iron[END_REF] using thermodynamic integration with the potential of Mendelev, showing a strong stress and temperature dependence and marked anharmonicity. These results suggest that an improved agreement with MD data could be obtained by determining Δ𝐺 𝑘𝑝 more accurately. In the following, we use direct Gibbs energy calculations to determine Δ𝐺 𝑘𝑝 , and implement our results in the well-established velocity models discussed in this section to verify the validity of commonly used expressions of Δ𝐺 𝑘𝑝 , using a minimum of adjustable parameters. 3.3 Beyond the harmonic approximation

Beyond the harmonic approximation

Approaches developed in the previous section are unable to accurately predict the flow stress at high temperature. The harmonic assumption which sets the basis of the harmonic TST is known to hold only at low temperatures, and anharmonic effects are expected to appear at higher temperature [START_REF] Lesar | Finite-temperature defect properties from free-energy minimization[END_REF][START_REF] Gilbert | Free energy generalization of the Peierls potential in iron[END_REF][START_REF] Sato | Anharmonic effect on the thermally activated migration of {1012} twin interfaces in magnesium[END_REF][START_REF] Swinburne | Uncertainty and anharmonicity in thermally activated dynamics[END_REF].

The PAFI method introduced in Sec. 1.3.1, which allows the evaluation of free enthalpy barriers associated with kink pair nucleation, is used in this section to explore a broad range of stress and temperature. Technical details regarding the PAFI method and our calculation methodology were already presented in Sec. 1.3.1. In our calculations, we use simulation cells that are twice smaller in the 𝑧 direction (i.e. the direction with no boundary conditions) compared to the MD simulations, in order to reduce the computational load. The cell dimensions are given by the following numbers of unit cells along 𝑥, 𝑦, 𝑧: (𝑁 𝑥 , 𝑁 𝑦 , 𝑁 𝑧 ) = (40,20,20). Additional tests conducted with (𝑁 𝑥 , 𝑁 𝑦 , 𝑁 𝑧 ) = (40, 15, 15) and (40, 30, 20) showed no significant effect of the cell dimension in the 𝑦, and 𝑧 directions. In Fig. 3.10, we present values of the free enthalpy of activation for kink pair nucleation Δ𝐺 𝑘𝑝 computed with applied stresses 𝜏 = 200 MPa and 500 MPa as a function of temperature.

Δ𝐺 𝑘𝑝 is calculated as the difference between the maximum and initial Gibbs energy along barriers such as those presented in Fig. 3.11 (a). In Fig. 3.10, the values of Δ𝐺 𝑘𝑝 computed with the PAFI method are also compared to the harmonic prediction obtained by diagonalization of the Hessian matrix (see Sec. 1.3). It appears that both methods are in agreement close to zero Kelvin, but quickly diverge as Δ𝐺 𝑘𝑝 becomes non-linear even at low temperature (𝑇 ≳ 20 K). This unexpected behavior means that harmonic entropy calculations, a wellestablished method for computing the low-temperature entropy [START_REF] Vineyard | Frequency factors and isotope effects in solid state rate processes[END_REF] would lead to a strong overestimation of the activation entropy in this case, even at low temperature.

The curvature of the energy landscape at the initial and activated states can therefore only be used very locally, limiting the harmonic regime to very low temperatures. When temperature allows the system to visit neighboring states around the equilibrium position, it experiences a curvature that differs from the one of an idealized harmonic surface. This overestimation of the entropy is likely due to a roughness of the potential energy landscape predicted by Proville potential and to locally flat regions in the energy landscape of the interatomic potential, which could be an artefact of the interatomic potential.

Nevertheless, the observed non-linearities remain small compared to the variations caused by the stress or the temperature. In Fig 3 .11, it can be seen that the overall curves remain close to a linear behavior at low temperature. We have therefore linearized the Δ𝐺 𝑘𝑝 data at low-temperature linear regressions is reported in Fig. 3.13 as a function of both the activation enthalpy Δ𝐻 𝑘𝑝 and the stress. It appears that Δ𝑆 𝑟𝑒𝑔 varies by several 𝑘 𝐵 , which contrasts with the approximation of a constant activation entropy. In addition, a clear proportionality is observed in both cases, which can be seen as two manifestations of a same effect, as Δ𝐻 𝑘𝑝 varies close to linearly with the stress (see Fig. In Fig. 3.14, we implemented an inverse Meyer Neldel law in the glide models of Sec. 3.2.2 using the parameterized value of 𝑇 𝑀𝑁 and Δ𝑆 0 , respectively 484 K and 12.3𝑘 𝐵 . Again, both models presented side-by-side use the same expression of Δ𝐺 𝑘𝑝 . With the inverse MN law, both models yield a satisfactory agreement with MD data. The thermally activated model appears slightly more accurate, and its predictions go through the error bar of most data points. We emphasize that this result is obtained with no adjustable parameter in the case of the model of Po et al., and with 𝜈 as the only adjustable parameter in the case of the thermally activated model (which was set to the same value 𝜈 = 5 × 10 10 s -1 throughout this chapter). Both models are reasonably accurate within the entire range of temperature and shear rates presented here. These results suggest that the temperature dependence of the flow stress seen in low temperature MD simulations can be explained within the framework of the harmonic TST, using the inverse MN law. At higher temperature, anharmonicity is expected, as a deviation from linearity of the Gibbs energy is seen in Fig. 3.12. Beyond the linear regime, a deviation from linearity is indeed observed in Fig. 3.12, corresponding to an anharmonic regime where Δ𝐺 𝑘𝑝 decreases non-linearly. We have checked that for the moderate stresses and temperatures considered here, both usually discussed anharmonic effects, i.e. thermal dilatation and weakening of the elastic moduli [START_REF] Ryu | Entropic effect on the rate of dislocation nucleation[END_REF], do not affect the energy barriers. The main effect is probably a widening of the activated kink pair [START_REF] Swinburne | Unsupervised calculation of free energy barriers in large crystalline systems[END_REF], which results in a deviation of the MGEP at finite temperature from the zero-Kelvin MEP. Interestingly, we see in Fig. 3.12 that at low applied stresses, Δ𝐺 𝑘𝑝 decreases faster than the linear regression while it is the opposite at high applied stresses. This implies that using a linear approximation at high temperatures may result in either over-or underestimating the dislocation velocity depending on the applied stress, resulting in erroneous curvature of flow stress curves in Fig. PAFI data (symbols) is fitted to Eq. 3.24, giving a set of parameters given in Tab 3.1. Each temperature data is represented in a different subplot, with the corresponding fit shown as a dark solid line. All curves are reported in the background of other sublots as light gray lines for easier comparison.

temperatures (including the anharmonic regime), the data can be accurately fitted by a phenomenological approach, using a modified Kocks law [START_REF] Kocks | Thermodynamics and kinetics of slip[END_REF]:

Δ𝐺 𝑘𝑝 (𝜏 , 𝑇 ) = Δ𝐺 0 (𝑇 ) (1 -( 𝜏 𝜏 𝑃 (𝑇 ) ) 𝑝 ) 𝑞 , (3.24)
with:

𝜏 𝑃 (𝑇 ) = 𝑎 𝑇 + 𝑇 0 + 𝑏, (3.25)
and Δ𝐺 0 (𝑇 ) = Δ𝐺 0 -𝑇 × 𝐺 𝑇 , with Δ𝐺 0 = 0.6 eV and 𝐺 𝑇 = 2 × 10 -4 eV K -1 . Here, 𝜏 𝑃 (𝑇 ) is analogous to the zero-temperature Peierls stress for a finite temperature 𝑇, since it verifies

Δ𝐺(𝜏 𝑃 (𝑇 ), 𝑇) = 0. (3.26)
At low temperature (i.e. below ∼ 100 K), 𝜏 𝑃 (𝑇 ) is strongly temperature-dependent, and becomes constant at higher temperature, as reflected by Eq. 3.25. The obtained set of fitting parameters 𝑎, 𝑏, 𝑇 0 , 𝑝, 𝑞 in Fe are given in Tab. 3.1. For temperatures below 100 K, a linear fit gives the best match to the data. The corresponding fit is shown as a function of the stress in Fig. 3.15. It can be seen that the proposed expression is able to fit Gibbs energy data within a wide range of temperature and stress in the anharmonic domain. With Eq. 3.24 and 3.25, we propose a phenomenological expression of Δ𝐺 𝑘𝑝 (𝜏 , 𝑇 ) that accounts for anharmonic effects encountered at high temperature.

We used Eq. 3.24 to compute Δ𝐺 𝑘𝑝 and determine the flow stress of iron using the two models previously discussed, with the results presented in Fig. 3.16. Following the approach of Fig. 3.14, models predictions are compared with MD data. As can be seen, the agreement of the anharmonic models with MD data is very satisfactory, and data points are well predicted by the models. We note that the agreement with MD is only very slightly improved compared to the harmonic inverse MN law, which appears accurate enough, especially considering its great simplicity. However, at higher temperature, we expect the anharmonic model to better predict the flow stress. But MD simulations could not be conducted at higher temperature to verify this assumption, due to the instability of the dislocation, which tends to exit the simulation cell through the free surfaces when the temperature reaches 400 K in our simulation setup.

The values of Δ𝐺 𝑘𝑝 obtained with PAFI (Eq. 3.24) can also be used to compare the predictions of the different models to experimental flow stress data, e.g. from [START_REF] Brunner | Temperature and Strain-Rate Dependence of the Tensile Flow Stress of High-Purity 𝛼-Iron below 250 K. I. Stress/Temperature Regime III[END_REF]. It must be noted that in the low temperature limit, a well-known discrepancy appears Th. activ. jump (Eq. 3.11) Po et al. (Eq. 3.16) Experiments (secondary axis)

Figure 3.17: Experimental flow stress compared to models predictions. Experimental data from [START_REF] Brunner | Temperature and Strain-Rate Dependence of the Tensile Flow Stress of High-Purity 𝛼-Iron below 250 K. I. Stress/Temperature Regime III[END_REF] is compared to the predictions of the thermally activated model, which uses our parameterization of Eq. 3.24. Note that experimental data is plotted on a separate axis to account for the usual lowtemperature discrepancy between TST models and experiments.

long segments.

When comparing to MD simulations, the accuracy of the predictions depends for a larger part on the assumptions used for the calculation of Δ𝐺 𝑘𝑝 , often more than on the velocity model itself. Comparison to MD also shows that common assumptions, such as a constant entropy, or the Meyer-Neldel law, fail to reproduce the flow stress of iron at high temperatures. We used a novel free energy integration method (PAFI) to evaluate Δ𝐺 𝑘𝑝 beyond the harmonic domain, showing strong anharmonic effects, and non-linearity at temperatures as low as 20 K that quickly diverge from the predictions of HTST methods based on a diagonalization of the Hessian matrix. Using the PAFI method, the Gibbs energy was computed in a large range of stresses and temperatures and fitted with simple laws to be implemented in kinetic models. The flow stress predicted using this fitted law is in excellent agreement with MD simulations up to temperatures that are above the room temperature.

PAFI data also show that the Gibbs activation energy tends to follow an inverse Meyer-Neldel behavior at temperatures below 100 K. Implementing an inverse Meyer-Neldel law parameterized on our Gibbs energy data yields a satisfactory agreement with flow stress values obtained from MD simulations up to the high temperatures. Complementary analysis of atomic vibration modes obtained by diagonalization of the Hessian matrix was performed following the approach of [START_REF] Gelin | Enthalpy-entropy compensation of atomic diffusion originates from softening of low frequency phonons[END_REF] as an attempt to understand the physical origins of this inverse Meyer-Neldel behavior. The results show that when the applied stress (and conversely, the enthalpy) changes, the vibrational modes of the iron matrix (i.e. delocalized in the simulation cell) are affected. We observe that in contrast to the observations of Gelin et al. ( 2020), the softening peak seen in the difference of vibrational density of states (VDOS)

between the initial and activated states is increased when the enthalpy decreases, resulting in an increase of the activation entropy. However, the strong deviation of the Gibbs energy from the HTST prediction starting at 20 K raises questions on the validity of these observations above very low temperatures. As the distance between the initial and activated states in configurational space is reduced when the stress is close to the Peierls stress (as the critical kink is then only a small buldge, as seen in Fig. 3.2), we expect the activation entropy to decrease in the high stress domain. We were able to confirm this prediction with complementary HTST calculations using the EAM4 potential for tungsten [START_REF] Marinica | Interatomic potentials for modelling radiation defects and dislocations in tungsten[END_REF].

But the same calculations conducted in iron revealed that Δ𝑆 increases even more stiffly under a high stress, corresponding to a softening of vibration modes, which is likely a local artefact of the interatomic potential.

At this point, this work is limited by the realism of the Proville interatomic potential. The strong anharmonicity of the potential at moderate temperature, as well as the non-linearity that appears at low temperature raise questions on the realism of the potential. These results should thus be confirmed with other interatomic potentials for iron that are valid for the simulation of dislocations, such as the recently developed machine learning-based potentials from [START_REF] Goryaeva | Efficient and transferable machine learning potentials for the simulation of crystal defects in bcc Fe and W[END_REF]. These potentials show an excellent agreement with the ab initio Peierls barrier and kink-pair nucleation energies derived from ab initio-parameterized line tension models.

Perspectives for the PAFI method

The PAFI method has shown a great efficiency for free-energy calculations, compared either to classical harmonic approaches or to recent developments in free energy calculation methods [START_REF] Grabowski | Thermally-activated dislocation mobility in bcc metals: An accelerated molecular dynamics study[END_REF]. In the future, this method could be generalized to a broad range of thermally activated processes, allowing a better understanding of these processes as well as more accurate predictions. Using this type of method, new interatomic potentials could be automatically compared to free energy data obtained e.g. with ab initio methods [START_REF] Ackland | Calculation of free energies from ab initio calculation[END_REF][START_REF] Zhang | Calculating free energies of point defects from ab initio[END_REF][START_REF] Grabowski | Ab initio vibrational free energies including anharmonicity for multicomponent alloys[END_REF], while it is common to benchmark potentials on zero-kelvin properties, which offers no guarantee that predictions remain valid at high temperature.

Several improvements to the PAFI method could be proposed to make it more easily usable by the community and faster. A first limitation of the method in its current state is its sensitivity to the density of knots along the initial minimum energy path (MEP), which must be sufficiently high in regions where the free energy varies strongly, to avoid integration error that would be caused by discontinuities of the free energy gradient. Here, we avoided this issue by using a large number (20-30) of NEB replicas and PAFI hyperplanes to compute each barrier, but this approach is highly inefficient and a similar result could be obained with less hyperplanes placed at optimal positions. One difficulty is that the position of the regions where the free energy varies strongly cannot be anticipated from 0 K calulations, as the position of the initial and transition states are typically shifted with increasing temperature (see Fig. 3.18). This issue could be addressed by adding knots on-the-fly at optimal positions where large free energy gaps remain, e.g. by adapting approaches used for adaptative NEB path construction [START_REF] Kolsbjerg | An automated nudged elastic band method[END_REF]. The hyperplane associated with a new knot is then defined by spline interpolation of the MEP, or using the free end NEB method for a more precise result, which allows to find a minimum energy path between non-equilibrium positions.

Reac. coord.

Free Energy NEB path

Free Energy gap increasing temperature shift of the intial state Another current limitation of the method is the time nedded to perform the several hundreds of independent samplings that are necessary to compute a barrier with a good precision. The parallelization scheme of PAFI, shown in Fig. 1.6, allows to run several samplings on the same hyperplane simultaneously, each on a different worker (group of CPUs), connected by the Message Passing Interface (MPI) with a perfect parallel scaling (i.e. no time is lost in communication between workers, as they are fully independent). While offering a simple way of starting a large number of simulations simultaneously and to collect their results, this approach uses a fixed number of workers, meaning that all CPUs have to be requested at once in shared high performance computing (HPC) platforms. On a busy computing platform, this results in a long queue time until all requested CPUs are available. On the other hand, if more CPUs are available on the plateform, only the requested number is used, thus not fully exploiting the available resources. These limitations could be lifted by submitting an independent job for each worker (or groups of workers) of all hyperplanes to the schedueler at once, which will then be allocated depending on the available resources.

Screw dislocation pinning by interstitial carbon in iron

The presence of interstitial carbon atoms can affect the thermally-activated nucleation and subsequent propagation of kink-pairs, even with low concentrations of solutes [START_REF] Caillard | An in situ study of hardening and softening of iron by carbon interstitials[END_REF]. While theories of solute hardening are well established in substitutional alloys [START_REF] Pachauri | Climate change 2014: synthesis report. Contribution of Working Groups I, II and III to the fifth assessment report of the Intergovernmental Panel on Climate Change[END_REF][START_REF] Varvenne | Solute strengthening in random alloys[END_REF], their application to interstitials remains difficult due to strong solute-dislocation interactions in the core region. In this chapter, we investigate the effect of interstitial carbon on dislocation mobility in Fe-C with atomistic simulations, confirming a powerful locking effect of carbon. We propose a model of dislocation pinning based on our model of dislocation glide in Fe (Chap. 3), which explains the origin of the strengthening effect induced by carbon.

Introduction

The addition of carbon in iron to improve mechanical properties is at the basis of steel metallurgy [START_REF] Berns | Ferrous Materials[END_REF]. For dilute alloys, the change of mechanical properties at low temperature is related to the interaction between solute atoms and dislocations1 .

In substitutional alloys, solutes have a low mobility and dislocations can be considered as moving through a field of fixed obstacles [START_REF] Zhao | Direct prediction of the solute softening-to-hardening transition in W-Re alloys using stochastic simulations of screw dislocation motion[END_REF]. With interstitials on the other hand, like carbon in iron, solutes are mobile and able to diffuse. The elastic field of dislocations favors their segregation around dislocation lines, forming Cottrell atomospheres [START_REF] Cottrell | Dislocation theory of yielding and strain ageing of iron[END_REF]) that were evidenced experimentally using atom probe techniques [START_REF] Blavette | Three-dimensional atomicscale imaging of impurity segregation to line defects[END_REF][START_REF] Wilde | Three-dimensional atomic-scale mapping of a Cottrell atmosphere around a dislocation in iron[END_REF]. Cottrell atmospheres can greatly reduce the mobility of dislocations, which necessitate a high stress to escape their atmosphere and be set into motion [START_REF] Schoeck | The flow stress of iron and its dependence on impurities[END_REF]. This effect is known as static strain ageing, and is typically associated with an important yield drop at the beginning of the plastic domain [START_REF] Kubin | On static strain ageing[END_REF], caused by the collective depinning of dislocations from carbon. In a domain of deformation rate and temperature where dislocations and carbon move at a comparable velocity, a regime of serrated flow is observed (dynamic strain ageing) where atmospheres are able to catch up with dislocations before they can move any further [START_REF] Zhao | Simulating the mechanisms of serrated flow in interstitial alloys with atomic resolution over diffusive timescales[END_REF] (see e.g. stress more than 10-fold compared to a solid solution of low concentration [START_REF] Veiga | Monte Carlo and molecular dynamics simulations of screw dislocation locking by Cottrell atmospheres in low carbon Fe-C alloys[END_REF].

When solutes are homogeneously distributed in solid solution, a hardening and a softening are observed experimentally depending on the temperature domain [START_REF] Caillard | An in situ study of hardening and softening of iron by carbon interstitials[END_REF], as represented schematically in Fig. 4.3.

Several solute strengthening theories were proposed to explain the strengthening in random solid solutions of substitutional atoms, based on interactions between dislocations and solutes at the atomic scale (Patinet and [START_REF] Pachauri | Climate change 2014: synthesis report. Contribution of Working Groups I, II and III to the fifth assessment report of the Intergovernmental Panel on Climate Change[END_REF][START_REF] Varvenne | Solute strengthening in random alloys[END_REF]. A classical conception of hardening, which is typical of fcc metals, is attributed to the bowing of elastic dislocation lines between pinning points under the applied stress [START_REF] Foreman | Dislocation movement through random arrays of obstacles[END_REF], which increases until it reaches either a critical angle 𝜙 𝐶 for weak obstacles, or until segments on each side of the obstacle collide, resulting in their annihilation, in the case of strong obstacles. In bcc metals, due a large Peierls barrier, this mechanism does not apply as straight dislocation segments are favored and glide is allowed by the kink pair mechanism (see Sec. 3.1).

In this case, strengthening can be explained by an attractive interaction between the solutes and the dislocation, impeding dislocation glide. Several authors implemented the solutedislocation interaction in analytical models of substitutional alloys [START_REF] Labusch | A Statistical Theory of Solid Solution Hardening[END_REF][START_REF] Suzuki | Dislocation Dynamics and Plasticity[END_REF][START_REF] Pachauri | Climate change 2014: synthesis report. Contribution of Working Groups I, II and III to the fifth assessment report of the Intergovernmental Panel on Climate Change[END_REF][START_REF] Varvenne | Solute strengthening in random alloys[END_REF][START_REF] Larosa | Solid solution strengthening theories of high-entropy alloys[END_REF][START_REF] Maresca | Theory of screw dislocation strengthening in random BCC alloys from dilute to "High-Entropy" alloys[END_REF][START_REF] Rao | Theory of solid solution strengthening of BCC Chemically Complex Alloys[END_REF][START_REF] Ghafarollahi | Solute/screw dislocation interaction energy parameter for strengthening in bcc dilute to high entropy alloys[END_REF], or simulations based on the Kinetic Monte Carlo (KMC) method [START_REF] Zhao | Direct prediction of the solute softening-to-hardening transition in W-Re alloys using stochastic simulations of screw dislocation motion[END_REF]. Recently, [START_REF] Zhao | Simulating the mechanisms of serrated flow in interstitial alloys with atomic resolution over diffusive timescales[END_REF] proposed KMC simulations of interstitial W-O alloys, parameterized on DFT simulations [START_REF] Zhao | Electronic Structure Calculations of Oxygen Atom Transport Energetics in the Presence of Screw Dislocations in Tungsten[END_REF], that were able to reproduce dynamical strain ageing behavior. A central parameter in With interstitials, local effects of solutes on the dislocation core were evidenced in DFT calculations [START_REF] Ventelon | Dislocation core reconstruction induced by carbon segregation in bcc iron[END_REF][START_REF] Lüthi | Attractive interaction between interstitial solutes and screw dislocations in bcc iron from first principles[END_REF], presented in Fig. 2.3. It consists in the stabilization of a prismatic core structure unstable in pure Fe and other bcc metals [START_REF] Lüthi | First principles investigation of carbon-screw dislocation interactions in body-centered cubic metals[END_REF], evidenced with C and with other interstitial solutes in iron [START_REF] Lüthi | Attractive interaction between interstitial solutes and screw dislocations in bcc iron from first principles[END_REF]. Due to its low mobility [START_REF] Lüthi | Modélisation ab initio des interactions entre dislocations vis et solutés interstitiels dans les métaux de transition cubiques centrés[END_REF], this reconstructed, prismatic core is expected to control the plasticity of bcc interstitial alloys when solutes are mobile and decorate the dislocations [START_REF] Hachet | Screw dislocation-carbon interaction in BCC tungsten: an ab initio study[END_REF]. But DFT simulations remain limited by their high computational cost as studies at a larger scale are necessary.

MD simulations were also used to study the effect of carbon on non-screw dislocations in showing a dramatic reduction of dislocation mobility, linked to a higher kink-pair nucleation enthalpy [START_REF] Wang | Atomistic simulations of carbon effect on kink-pair energetics of bcc iron screw dislocations[END_REF]. However, the interatomic potentials used in these studies show a poor agreement with dislocation properties in iron computed with DFT -including the shape and height of the Peierls barrier-which makes their conclusions open to question. In chapter 2, we compared Fe-C potentials from the literature and proposed a novel potential which allows the simulation of screw dislocations interacting with carbon. In Chap. 3, we investigated the temperature dependance of the flow stress in iron, and determined that using this potential, the harmonic approximation remains valid up to the ambiant temperature in typical MD simulations conditions. At higher temperatures in the Fe-C system, a different regime of deformation is seen experimentally, attributed to the effect of edge dislocations that become less mobile than screws [START_REF] Caillard | An in situ study of hardening and softening of iron by carbon interstitials[END_REF].

Based on the hybrid interatomic potential developed in Chap. Solute atoms outside this volume are left unchanged between the different simulations. In Fig. 4.5 (a), we compare the stress-strain curves obtained with a full solid solution (𝑛 𝑝𝑙𝑎𝑛𝑒𝑠 = 0), and a carbon-depleted zone around the dislocation glide plane comprised of 𝑛 𝑝𝑙𝑎𝑛𝑒𝑠 atomic planes on both sides. In the case 𝑛 𝑝𝑙𝑎𝑛𝑒𝑠 = 2, the vicinity of carbon atoms attracts the dislocation, which cross slips and gets pinned to the carbon, resulting in a strenghtening that is comparable to the full solid solution. When 𝑛 𝑝𝑙𝑎𝑛𝑒𝑠 is larger than 4 (corresponding to approximately 8 Å), the stress level reveals no significant strengthening in comparison to pure iron (green curve). This advocates for a pinning mainly controlled by short-range interactions.

Unpinning from an array of C atoms

Given the complexity of dislocation glide in a random solid solution seen in previous section, we defined a model system when a screw dislocation detaches from an array of carbon atoms directly located in its core. This allows to exclude cross-kink pinning, and the effect of mid-and long-range interactions, to focus on the elementary unpinning process. For these simulations, two C atoms were inserted in second-nearest octahedral positions along the dislocation line (denoted 𝑂 (2) , see Fig.

Effect of solutes spacing

We investigated the effect of solute spacing by varying ℓ and 𝐿 on dislocations of varying ℓ+𝐿 lengths. However, we have seen that the length of the shorter ℓ segment has a negligible effect, such that all results can be presented in Fig. 4.9 as a function of 𝐿 only. The figure also contains the flow stress computed in pure iron, which does not show any significant dependence on the dislocation length 𝐿. This is coherent with the fact that without obstacles, the flow stress depends only logarithmically on 𝐿. The reason is that the flow stress is solution of Orowan equation γ = 𝜌𝑏𝑣, with 𝜌 the density of mobile dislocations (here 𝜌 = 1/𝐿 𝑦 𝐿𝑧)

and 𝑣, the dislocation velocity. As introduced in Sec. 3.2.2, the velocity 𝓋 of a dislocation of length 𝐿 can be expressed by a thermally activated forward-jump model: By way of contrast, the average unpinning stress in presence of C atoms shows a very strong 𝐿-dependence when 𝐿 100 𝑏. 𝜏 * reaches more than 1 GPa when solutes are closest, meaning that the stress required to overcome them is more than twice the flow stress and is actually even higher than the 0 K Peierls stress in pure Fe (∼ 1 GPa). This demonstrates the powerful locking capability of an array of C atoms, with a spacing in the range 5 to ∼ 50 𝑏 when the temperature is small enough that the solute atoms cannot follow the dislocation in its glide. This very strong pinning for closely spaced C atoms is consistent with DFT calculations [START_REF] Lüthi | Attractive interaction between interstitial solutes and screw dislocations in bcc iron from first principles[END_REF]. For larger C spacings, the unpinning stress reduces to the flow stress in pure iron and no strengthening or weakening is observed when 𝐿 100 𝑏.

In addition, as mentioned above, the data show only a small difference in 𝜏 * whether the dislocation is pinned by evenly (ℓ = 𝐿, red squares) or unevenly spaced C atoms (ℓ ≤ 𝐿, blue circles). The latter case was investigated by repeated simulations with different values of ℓ, ranging between 0 and 𝐿, which is why different data points share the same 𝐿 in Fig. 4.9.

The reason for a negligible influence of the ℓ segment is that unpinning is controlled by the nucleation of the second kink pair (step 5 of Fig. 4.7), which has a higher probability to occur on the longer 𝐿 segment.

Effect of C atoms on kink pair nucleation enthalpy NEB calculations

In order to better understand the origin of the very strong pinning induced by the C atoms at short distance, we used the NEB method to identify the minimum energy paths that connect the different stages of unpinning identified in MD simulations (i.e stages presented in Note that the enthalpy globally decreases due to the work of the applied stress, which tilts the potential energy landscape.

The unpinning process involves the successive nucleation of two kink-pairs, between 𝑆 1 and 𝑆 2 and between 𝑆 2 and 𝑆 3 . As can be seen in Fig. 4.10 (a), the kink-pairs in the activated states (𝑆 ⋆ 1 and 𝑆 ⋆ 2 ) have not formed at the C atom, but rather in-between the C atom and one of its periodic images. The configuration 𝑆 2 contains two kinks on either side of the C atom. It is metastable. This confirms the MD observation that kinks of opposite signs cannot spontaneously annihilate when there is a C atom in-between them. The enthalpy barrier to form this first kink pair is noted Δ𝐻 1 . A second barrier of height Δ𝐻 2 has to be overcome to form a second kink pair (state 𝑆 ⋆ 2 ) and release the dislocation. This also confirms the MD observation that a kink crossing two Peierls valleys can overcome a C atom and unpin the dislocation. The next valleys (not shown here) are crossed with activation barriers close to that in pure Fe, Δ𝐻 Fe . Fig. 4.10 (c,d) show examples of activation barriers for dislocation lengths 𝐿 = 20 and 40 𝑏 as well as in pure Fe. When 𝐿 is small (20 𝑏), we observe a significant increase of Δ𝐻 2 compared to Δ𝐻 1 , which on the other hand remains close to Δ𝐻 Fe . Therefore, the C atom does not affect significantly the nucleation of the first kink-pair, but the presence of this kink pair which cannot annihilate because of the C atom makes nucleation of the second kink pair much harder. This increase is at the origin of the strengthening seen in the MD simulation.

Conversely, larger values of 𝐿 (60 𝑏) yield almost no change compared to pure iron. When the C atoms are widely separated, the kink-pairs, which form away from the C atoms, no longer feel the presence of the solute atoms. We note a slight decrease of Δ𝐻 1 , implying that the C atoms help somewhat the nucleation of the first kink-pair. But since Δ𝐻 2 , which controls the unpinning, matches Δ𝐻 Fe , no softening is produced, as confirmed by the MD simulations. This is different from substitutional solid solutions, where solutes can help the nucleation of kink-pairs and induce a softening effect [START_REF] Ghafarollahi | Solute/screw dislocation interaction energy parameter for strengthening in bcc dilute to high entropy alloys[END_REF].

Statistical model

In order to connect the energy barriers with the MD simulations, we used a statistical model based on the TST to predict the average unpinning stress. This model is adapted from Ref. [START_REF] Rodney | Activation enthalpy for kink-pair nucleation on dislocations: Comparison between static and dynamic atomic-scale simulations[END_REF], which considered thermally-activated glide in absence of obstacles.

The central quantity of the statistical model is the survival probability, 𝑊 (𝑡, 𝑡 ′ ), the probability that the system in a metastable state at time 𝑡 remains in that state until 𝑡 ′ .

𝑊 𝑖 (𝑡, 𝑡 ′ ) = exp (-∫ 𝑡 ′ 𝑡 𝑝 𝑖 ( τ 0 𝑢) 𝑑𝑢) (4.3)
In the present case, we need to consider two survival probabilities: 𝑊 1 in 𝑆 1 and 𝑊 2 in 𝑆 2 .

The probability that, starting at 𝑡 = 0 in an unstrained system, the dislocation is still pinned at time 𝑡 is expressed as:

𝑊 (𝑡) = 𝑊 1 (0, 𝑡) -∫ 𝑡 0 𝜕𝑊 1 𝜕𝑡 ′ (0, 𝑡 ⋆ )𝑊 2 (𝑡 ⋆ , 𝑡)𝑑𝑡 ⋆ , (4.4)
where the first term is the probability that the dislocation is still in state 𝑆 1 at time 𝑡 and the 4.3 Unpinning from an array of C atoms second term is the probability that the dislocation jumps from 𝑆 1 to 𝑆 2 at time 𝑡 ⋆ ∈ [0, 𝑡] but then remains in 𝑆 2 until 𝑡. The survival probabilities depend on the transition rate, i.e. the probabilities per unit time, 𝑝 𝑖=1,2 (𝜏 , 𝐿, 𝑇 ), that the dislocation jumps from 𝑆 1 to 𝑆 2 (𝑖 = 1) and from 𝑆 2 to 𝑆 3 (𝑖 = 2). These rates are functions of the dislocation length 𝐿 and the temperature 𝑇 as well as the internal stress 𝜏. Since the MD simulations are stain-rate controlled, the stress increases linearly with time, at a rate τ 0 = 𝜇 γ , where 𝜇 is the shear modulus. The survival probability is then expressed as:

𝑊 (𝑡) = exp(-∫ 𝑡 0 𝑝 1 ( τ 0 𝑢) 𝑑𝑢) + ∫ 𝑡 0 [𝑝 1 ( τ 0 𝑡 * ) exp(-∫ 𝑡 * 0 𝑝 1 ( τ 0 𝑢) 𝑑𝑢) exp(-∫ 𝑡 𝑡 * 𝑝 2 ( τ 0 𝑢 -Δ𝜏 ) 𝑑𝑢)] 𝑑𝑡 * .
(4.5)

The probability to unpin at time 𝑡 is given by 𝜔(𝑡) = -𝑑𝑊 /𝑑𝑡 and the average unpinning stress is thus:

⟨𝜏 ⋆ ⟩ = 𝜇 γ ∫ ∞ 0 𝑡𝜔(𝑡)𝑑𝑡 (4.6)
To compute ⟨𝜏 ⋆ ⟩, we thus need to express the rates, 𝑝 𝑖=1,2 (𝜏 , 𝐿, 𝑇 ). Within the TST, they are written as: 𝑝 𝑖 (𝜏 , 𝐿, 𝑇 ) = 𝜈𝐿 exp ( -Δ𝐺 𝑖 (𝜏 , 𝐿, 𝑇 ) 𝑘𝑇 ) (4.7)

where 𝜈 is the attempt frequency and Δ𝐺 𝑖 is the Gibbs activation energy of the process.

As extensively discussed in Chap. 3, Δ𝐺 is notoriously difficult to compute and approximations are often used. While we were able to directly compute Δ𝐺 in the case of a dislocation gliding in pure iron, the present case is more complex and the use of the PAFI method remains challenging. In our tests, difficulties arise from the competition between carbon migration and kink pair formation, which both can happend during the constrained MD sampling done by PAFI. Technical concerns on how to inhibit carbon migration while not affecting crystal vibration modes remain difficult to address.

We therefore use approximations to compute Δ𝐺 𝑖 , based on the results obtained in pure Fe (see Chap. 3). We consider that Δ𝐻 and Δ𝑆 do not depend on 𝑇 and Δ𝐻 (𝜏 , 𝐿) is directly obtained from the zero Kelvin MEP computed in previous section. At low temperature, we derived a stress/entropy linear relation from our Gibbs energy calculations in Fe:

Δ𝑆 ≃ 𝑘 𝐵 𝛼 𝐹 𝑒 𝜏 (4.8)

with 𝛼 𝐹 𝑒 = 1.6 × 10 -2 MPa -1 , which compared well with MD simulations. Here, we also use the attempt frequency 𝜈 = 5 × 10 10 s -1 already used in Chap. 3.

Computing Eq. 4.7 requires a continuous expression of Δ𝐻 𝑖 (𝜏 , 𝐿). For this purpose, we used again a classical Kocks law [START_REF] Kocks | Thermodynamics and kinetics of slip[END_REF]: The unpinning stresses obtained from Eq. 4.5 with the continuous approximation of Δ𝐻 𝑖 in Eq. 4.9 and the attempt frequency and activation entropy computed in pure iron are shown as a dashed line in Fig. 4.9. The model reproduces the strong pinning effect at short solute spacing, which is caused by the increased activation enthalpy, and the convergence to the flow stress in pure iron at large spacing. In-between, the model tends to underestimate the pinning stress.

Δ𝐻
Length-dependent entropy We have made a strong simplification by assuming that the activation entropy Δ𝑆 was not affected by the solute atoms. We can reasonably expect that the presence of solute atoms along the dislocation will affect the vibrational modes of the system and thus modify the activation entropy. We attempted to perform free energy barrier calculations with the PAFI method in the case of Fe-C, but faced difficulties as carbon tends to migrate to neighbouring sites during the sampling at finite temperature. In the limit of

Conclusions

We used atomistic simulations to study the interaction of screw dislocations with interstitial carbon atoms in iron. Simulations of full solid solutions showed a strong locking caused by the short-range attraction between the dislocation and the solutes. The contributions of all carbon atoms close to the dislocation line result in a collective effect that amplifies the strengthening.

In order to clarify the origin of this strengthening and to overpass the difficulty to analyse such collective effect in MD simulation, we focus on elementary interaction processes between a screw dislocation and an array of carbon atoms. Based on enthalpy calculations, we built a statistical model which is able to reproduce MD results with good agreement.

Though generally neglected, the effect of entropy discussed in Chap. 3 was again significant in the case of the Fe-C system. To take it into account, we incorporated a stress-and lengthdependent harmonic activation entropy law using a phenomenological approach, to avoid explicitly computing the harmonic entropy. Indeed, diagonalizing the Hessian matrix for the many possible configurations would result in a high computational cost, and instabilities around the MEP already mentioned could flaw the results. A more detailed study, based on free energy calculations in the FeC system for different solute spacings and levels of stress, would be very valuable and is an important perspective of this work.

Conclusion and perspectives

The addition of carbon in steels to improve mechanical properties make these alloys prone to ageing, which can be detrimental over long timescales, especially in the case of nuclear pressurized water reactors where it could limit the operation lifetime. Steel ageing is caused by the locking of dislocations by carbon solutes, which strongly impedes plastic deformation. In this work, we investigated atomic-scale mechanisms that enable the mobility of dislocations in iron and carbon steels, in order to develop a mobility law of dislocations in aged steels, transferable to larger-scale models.

Chapter 2: Assessment of interatomic potentials A first barrier for atomistic simulations of dislocations in the Fe-C system is the accuracy of available interatomic potentials. In Chap. 2, we tested all Fe-C empirical potentials available from the reference NIST database [START_REF] Becker | Considerations for choosing and using force fields and interatomic potentials in materials science and engineering[END_REF][START_REF] Hale | Evaluating variability with atomistic simulations: The effect of potential and calculation methodology on the modeling of lattice and elastic constants[END_REF] on screw dislocations and alloy properties, with a special focus on the Peierls mechanism, the carbon-induced screw dislocation core reconstruction, and the carbon-dislocation binding energy. Our calculations show that all tested potentials fail to reproduce some or all of these properties accurately compared to reference ab initio data. We therefore implemented a simple combination of two EAM potentials for Fe and FeC with complementary properties, the Fe potential from [START_REF] Proville | Quantum effect on thermally activated glide of dislocations[END_REF] and the Fe-C from [START_REF] Becquart | Atomistic modeling of an Fe system with a small concentration of C[END_REF] modified in [START_REF] Veiga | Comments on "atomistic modeling of an Fe system with a small concentration of C[END_REF], in a way that Fe-Fe interactions are unchanged compared to the potential of [START_REF] Proville | Quantum effect on thermally activated glide of dislocations[END_REF]. The resulting hybrid potential shows a satisfactory agreement with ab initio data, with the limited computational footprint of an EAM potential -while avoiding the heavy task of fitting a new interatomic potential from an ab initio dataset. This hybrid potential opens the way for atomistic simulations of carbon-dislocation interaction mechanisms, such as solute strenghthening.

However, one caveat of the hybrid potential is the roughness of its energy landscape along the carbon migration path between octahedral positions, which makes this potential unsuited for simulations involving diffusion mechanisms, where the potential of [START_REF] Becquart | Atomistic modeling of an Fe system with a small concentration of C[END_REF] remains preferred. At the time of writing, it appears that there is still no available poten-tial suited for both screw dislocation glide and carbon diffusion simulations. The increasing use of machine learning methods for interatomic potentials development will likely help develop new potentials with better accuracy and transferability, thus overcoming the limitations discussed here, perhaps with a computational cost several orders of magnitude larger than classical EAM (see e.g. [START_REF] Mishin | Machine-learning interatomic potentials for materials science[END_REF] for a review).

Potentials testing remains a crucial methodological step, which is often overlooked and cannot be comprehensive, due to the heavy task of setting up calculations of many different properties, and multiple methods can often be used to determine the same property, possibly leading to different results. Projects such as the iprPy tool [START_REF] Hale | Evaluating variability with atomistic simulations: The effect of potential and calculation methodology on the modeling of lattice and elastic constants[END_REF] or OpenKIM database [START_REF] Elliott | Knowledgebase of Interatomic Models (KIM) Application Programming Interface (API)[END_REF] address this issue by automating potential properties calculations, making useful results freely available, but they still lack many important material properties, and need to be extended by contributions from the community.

Chapter 3: dislocation glide in Fe In pure bcc iron, plastic deformation is controlled by the mobility of dislocations, which is limited by the lattice resistance at low temperature. In Chap. 3, we computed the flow stress of iron using shear-rate controlled MD simulations of isolated screw dislocations based on the potential of [START_REF] Proville | Quantum effect on thermally activated glide of dislocations[END_REF]. We tested several classical dislocation glide models based on the harmonic transition state theory as well as often-used assumptions for the calculation of the Gibbs energy of activation for kink pair nucleation (constant entropy, Meyer-Neldel law), which fail to reproduce our MD data in the entire range of temperature.

The recent projected average force integrator (PAFI) method enables direct Gibbs activation energy Δ𝐺 calculations for thermally-activated processes in large systems, with no assumption on the nature of atomic vibrations, i.e. unlocking the usual limitations of entropy calculations based on harmonic TST, which scale cubically with the number of atoms in the system, and become inaccurate at high temperature due to anharmonic effects [START_REF] Swinburne | Uncertainty and anharmonicity in thermally activated dynamics[END_REF]. Calculations in a broad range of stresses and temperatures show a strong effect on Δ𝐺, with marked anharmonic effects. A deviation from harmonic TST calculations is observed starting at temperatures as low as 20 K, suggesting that harmonic calculations do not always provide an accurate picture of the Gibbs energy landscape in the vicinity of the initial and transition states, even at low temperature.

Since non-linearities remain small below ∼100 K, the Gibbs activation energy was linearized to extract an effective activation enthalpy, which interestingly increases when the enthalpy decreases, corresponding to an inverse Meyer-Neldel behavior. Integrating the results of PAFI simulations into kinetic models of dislocation glide, we obtained a much improved agreement with direct MD simulations compared to classical approaches. Using the inverse Meyer-Neldel law, we were are able to reproduce MD simulations up to 200-300 K. These results advocate for a better integration of entropic and anharmonic effects that are often overlooked in atomistic studies of crystal defects. In the present case, the low-temperature non-linearity could be caused by local irregularities of the Gibbs energy landscape of the potential of [START_REF] Proville | Quantum effect on thermally activated glide of dislocations[END_REF], and should be verified using other atomic interaction models, such as the newly-developed machine-learning potential of [START_REF] Goryaeva | Efficient and transferable machine learning potentials for the simulation of crystal defects in bcc Fe and W[END_REF].

PAFI is an efficient and versatile free-energy calculation method which can be very useful for the study of a broad range of thermally activated processes. Yet, this recent method still has to overcome some challenges to be more easily usable by the community and more efficient, especially regarding its sensitivity to the initial MEP or its parallelization scheme.

Chapter 4: dislocation glide and pinning in Fe-C In steels, dislocation mobility is limited not only by the lattice resistance but also by their interactions with interstitial carbon atoms. In Chap 4, we use the hybrid interatomic potential developed in Chap. 2 to simulate dislocations interacting with a random solid solution of carbon atoms.

Our simulations confirm a local reorganization of the dislocation core as previously evidenced in ab initio simulations, and a strong dislocation pinning caused by complex shortranged carbon-dislocation interactions. In order to investigate the elementary pinning mechanisms, we studied the idealized case of a dislocation interacting with a row of carbon atoms, using a combination of MD simulations, saddle-point search techniques, and a harmonic TST kinetic model based on the calculations presented in Chap. 3. Our results indicate that carbon stabilizes the first kink pair by preventing kink annihilation, likely due to the local reconstruction of the dislocation core. As a result, a second kink pair has to form for the dislocation to unpin, with a markedly increased activation enthalpy. This effect causes a high strenghtening effect, with an unpinning stress that can be more than twice the flow stress of pure iron for closely-spaced solutes.

A more comprehensive investigation of entropic effects would be necessary, as these are affected by the presence of solutes, but remains challenging. Preliminary tests with the PAFI method revealed that under finite temperature and applied stress, carbon atoms tend to jump to other sites during the sampling. We presented a simple phenomenological approach to illustrate how a segment length-dependent activation entropy could improve the agreement with MD simulations, but Gibbs energy calculations in this system remain challenging, and might require adaptations in the PAFI code. Yet, investigating anharmonicity in alloys would have a great interest for studies above room temperature e.g. on serrated flow [START_REF] Zhao | Simulating the mechanisms of serrated flow in interstitial alloys with atomic resolution over diffusive timescales[END_REF], where we expect anharmonic effects to be dominant. The present work remains limited by the accuracy of the interatomic potential we used, and further investigation, especially concerning entropic effects, is difficult without a more accurate potential for Fe-C.

We have shown here that the elementary interaction process between a screw dislocation and carbon atoms involves the successive nucleation of two kink pairs, since a carbon atom is able to block a single kink and prevent its annihilation with a kink on the other side at least on MD timescale.

Our results in Fe and Fe-C could be combined and implemented in mesoscale models, which constitutes the next step towards a constitutive law of aged steel. In Fig. 1, the flowchart for multiscale simulation of plasticity proposed by Tipping ( 2010) is adapted to highlight the contribution of the present work, as well as connections to larger scale models.

In particular, integrating both the effect of the lattice resistance at high temperature and the effect of interstitial solute (which could be treated using local rules) in kinetic Monte Carlo

This work

Mesoscale models The contribution of the present work is emphasized (in green), as well as its main perspectives (in blue).

  Merci d'être venu enrichir cette thèse avec tes idées (et ton humour !) et d'avoir donné de ton temps pour m'aider avec PAFI ! Je souhaite ensuite remercier les membres de mon jury, dont les travaux ont m'ont beaucoup inspiré et qui m'ont fait l'honneur d'évaluer mon travail. Je remercie également mes collègues de l'IRSN, de toutes disciplines, pour les échanges iii Contents très enrichissants (et trop rares !). Je veux aussi remercier Sébastien et ses collègues du CCIN2P3 pour leur aide et leurs conseils très précieux.

Figure 0 . 1 :

 01 Figure 0.1: Share of the world total electricity production from fossil fuels, nuclear and renewables. Source: Our World in Data based on BP Statistical Review of World Energy & Ember (2021)

  Figure 0.3: Multi-scale mapping of the main techniques used in metals plasticity modelling.

  The need of an accurate IP for the transfer of information from ab initio to atomistic simulations constitutes a first barrier. A survey of existing IPs for Fe and Fe-C presented in Chap. 2 exposes the lack of a Fe-C potential with satisfying dislocation properties. We propose a novel Fe-C EAM potential, based on two potentials of the litterature, which best reproduces dislocation properties in iron compared to DFT, and unlocks the simulation of dislocation pinning by C.Dislocation mobility in pure iron In pure bcc metals, dislocation mobility is strongly limited by the resistance of the lattice, which has to be well understood before integrating the effect of solutes. In Chap. 3, we show that existing models of dislocation glide fail to reproduce the flow stress obtained from low-temperature MD simulations in pure iron. We compare different models and classical assumptions used within the harmonic transition state theory to anharmonic Gibbs energy calculations by leveraging the recent projected average force integrator method. Our results unveil strong entropic effects, usually overlooked in dislocation glide models, and anharmonicity especially at high temperature. We propose a model of dislocation glide in Fe in excellent agreement with MD observations, setting the fundation for a model of dislocation glide in the Fe-C system.Effect of pinning by interstitial solutes In Chap.4, the effect of carbon atoms on dislocation mobility is studied atomistically. MD simulations of dislocation glide in a random solid solution confirm a powerful solute strengthening caused by complex processes. We analyse it using a model geometry, where a row of aligned carbon atoms is inserted in the dislocation core. We use a combination of MD simulations, minimum-energy path calculations and extend a stochastic model validated in iron to explain the strengthening induced by carbon. We unveil that carbon disrupts the glide process, as unpinning necessitates the successive nucleation of two kink-pairs, with an activation enthalpy 0 Introduction markedly increased compared to iron at low C-C separations.

Figure 1 . 1 :

 11 Figure 1.1: Quadrupolar arrangement of dislocations typically used in DFT simulations. The simulation cell (grey shape) contains two dislocations of burgers + ⃗ 𝑏 and -⃗ 𝑏. Figure extracted from (Ventelon, 2008).

Figure 1 . 2 :

 12 Figure 1.2: Simulation cell used for dynamical simulations. A screw dislocation (indicated by a green line) is inserted along the [111] direction in the middle of a crystal of bcc iron. Two regions denoted 𝑆 ± are used to shear the crystal (in blue). Atoms visualization was obtained with Ovito.

Figure 1 . 5 :

 15 Figure 1.5: Strategy used to prepare a MEP for PAFI, based on a NEB path with sparse knots density along the barrier. (a) Knots located after the saddle state (𝑟 > 𝑟 𝑚𝑎𝑥 ) and below an energy treshold (𝐸 𝑡ℎ ) are eliminated. (b) A series of evenly spaced knots (in blue) is obtained by linear interpolation of atoms position between the selected knots (in green). (c) The path created in (b) is relaxed to a MEP using the "free-end NEB" method (see List. 1.4). Flat segments are then eliminated to avoid spline issues when interpolating the MEP in PAFI.

  neb ${spring} perp ${spring} end last ${spring} neb 0.0 1e-5 10000 10000 50 each neb_pos_file.$u Listing 1.4: LAMMPS code snipppet used to perform "free-end NEB".

Figure 1 . 6 :

 16 Figure 1.6: Parallel scheme as implemented in the PAFI code. Here, the first two hyperplanes are sampled by 2 workers, with 2 repetitions, with 4 CPUs per worker.

Figure 2 . 1 :

 21 Figure 2.1: Differential displacements of atoms in the dislocation core region. See main text for interpretation of. Different structures were obtained: (a) compact core (DFT, all potentials except Lee and Henriksson); (b) degenerate core (Lee); (c) symmetrical core (DFT); (d) asymmetrical core obtained with hybrid potential and other potentials (see text). Circles represent the relaxed position of the atoms. Red circles represent carbon atoms. Head-up triangles correspond to easy positions, and head-down triangles correspond to hard positions.

Figure 2 . 2 :

 22 Figure 2.2: Naming convention for octahedral sites around a screw dislocation (green cross). Reproduced from (Lüthi, 2017).

Figure 2 . 3 :

 23 Figure 2.3: Reconstruction of the dislocation core in the presence of C obtained with DFT. (a) Dislocation core, with an array of C atoms in 𝑂 (1) site every 𝑏, before (left) and after (right) performing an energy relaxation. (b) and (c) are the corresponding DD map for each system. Reproduced from[START_REF] Lüthi | Modélisation ab initio des interactions entre dislocations vis et solutés interstitiels dans les métaux de transition cubiques centrés[END_REF].

Figure 2

 2 Figure 2.6: Peierls barrier and corresponding core trajectory. DFT data (Ventelon and Willaime, 2007) is compared to different potentials. Figure adapted from Ref. (Mori and Ozaki, 2020).

Figure 2 . 8 :

 28 Figure 2.8: Normalized computational time for different Fe interatomic potentials.Time is normalized by 7.9 × 10 -7 s/atom/step. Data extracted from[START_REF] Starikov | Angular-dependent interatomic potential for large-scale atomistic simulation of iron: Development and comprehensive comparison with existing interatomic models[END_REF], where we added our timing for the Mori NN potential[START_REF] Mori | Neural network atomic potential to investigate the dislocation dynamics in bcc iron[END_REF], which is ∼ 10 3 times slower than Mendelev. Potentials names follow the convention adopted in Ref.[START_REF] Starikov | Angular-dependent interatomic potential for large-scale atomistic simulation of iron: Development and comprehensive comparison with existing interatomic models[END_REF].

Figure 2

 2 Figure 2.10: Carbon-dislocation binding energy, as a function of carbon spacing 𝑑 𝑐𝑐 along the dislocation line.Values obtained with different interatomic potentials are compared to ab initio data[START_REF] Lüthi | Modélisation ab initio des interactions entre dislocations vis et solutés interstitiels dans les métaux de transition cubiques centrés[END_REF]. Horizontal lines serve as guides for the eye.

Figure 2

 2 Figure 2.11: Carbon-carbon binding energy in bulk iron as a function of carbon spacing 𝑑 𝑐𝑐 in the [111] direction.Values obtained with different interatomic potentials are compared to ab initio data[START_REF] Becquart | Atomistic modeling of an Fe system with a small concentration of C[END_REF].

  2.7), along with the correct dislocation core structure, but this potential is limited to pure iron. For most Fe-C properties, there is a relatively large gap between interatomic potentials and DFT values, and none of the potentials stand out. Yet, Veiga potential qualitatively agrees with DFT on carbon-carbon and carbon-dislocation binding energies. It was used in previous studies to investigate carbon-dislocation interactions, Cot-

  Interatomic potentials were tested in regard to different Fe and Fe-C properties, yielding the following conclusions: Performance For the simulation of dislocations, only the most computationally efficient potentials can be used. They are potentials based on the EAM, MEAM or ADP formalisms. All tested Fe-C potentials follow these formalisms and are thus relatively fast. Dislocation core Most tested potentials feature a compact non-degenerate core in iron, except the potentials of Lee and Henriksson. The core reconstruction induced by C is reproduced by most potentials, with artefacts. Peierls barrier Proville potential provides the best quantitative agreement with DFT over all reviewed potentials. Excluding the most computationally expensive potentials, Proville is the only potential to agree qualitatively with DFT (single maximum). It is thus the only available potential for Fe to provide a physical representation of screw dislocations with good efficiency. C-dislocation binding energy The potentials of Henriksson, Lee, and the Hybrid potential are the only ones to accurately reproduce the DFT value of the C-dislocation binding energy. Solute strengthening effects simulated with other potentials are expected to differ from DFT predictions. Energy landscape The Hybrid potential features a rough energy landscape around the carbon migration path, and should therefore be used with caution for simulations involving the displacement of C. In this case, computational techniques that are highly 2.5 Conclusion and recommendations dependent on the local topology of the energy landscape (e.g. Hessian matrix diagonalization) might provide incorrect results.

  Figure 3.1: Screw dislocation motion in the bcc lattice. (a) Schematic representation of the (1 1 0) glide plane and maximum resolved shear stress plane (MRSSP). (b) Projection of the BCC lattice in the (111) plane showing the high-symmetry dislocation positions (E, easy core; H, hard core; S, split core). (c) Example of dislocation trajectory (green line) calculated with DFT in Mo and. (d) Energy profile along the hard-split line calculated with DFT in Mo. Reproduced from Dezerald et al. (2016).

  2.6): in the case of the Proville potential for iron, 𝛼 ∼ 24°[START_REF] Dezerald | Modélisation ab initio des dislocations vis dans les métaux de transition cubiques centrés[END_REF].The kink pair mechanism At finite temperature, thermal activation can help a short segment of the dislocation to cross the barrier, while it remains connected to the rest of the line by two non-screw segments called kinks, forming a kink pair. Under an applied stress, kinks experience a Peach-Koehler force that tends to move them in opposite directions along the line, thus expanding the kink pair as shown in Fig 3.2, and allowing the entire line to cross

Figure 3 . 2 :

 32 Figure 3.2: Sketch of the nucleation and migration of a kink-pair. The red line is a dislocation of burgers vector ⃗ 𝑏 = 1 2 [111]. The periodic energy landscape on which the dislocation is gliding is represented by a sinusoidal gray surface. Green arrows represent kinks migration due to the Peach-Koehler force.

  [111] direction, making Δ𝐻 (𝜏 = 0) L-dependent. The fitting parameter Δ𝐻 0 of the Kocks law provides an estimate of Δ𝐻 (𝜏 = 0) by extrapolation. The value obtained with Proville potential, close to 0.60 eV[START_REF] Proville | Quantum effect on thermally activated glide of dislocations[END_REF][START_REF] Ghafarollahi | Theory of double-kink nucleation in dilute BCC alloys[END_REF] can be compared to predictions of line tension models parameterized on ab initio data. Depending on the line tension model and ab initio parameterization, reported values for 𝐸 𝑘𝑝 are 0.86 eV[START_REF] Proville | Prediction of the kink-pair formation enthalpy on screw dislocations in 𝛼-iron by a line tension model parametrized on empirical potentials and first-principles calculations[END_REF], 𝐸 𝑘𝑝 = 0.91 eV (Dezerald et al., 2016) (using the same tension model but different ab initio parameterization) and 𝐸 𝑘𝑝 = 0.73 eV[START_REF] Itakura | First-principles study on the mobility of screw dislocations in bcc iron[END_REF] with a different model and ab initio code. A precise determination of 𝐸 𝑘𝑝 from experiments remains difficult[START_REF] Dezerald | First-principles prediction of kink-pair activation enthalpy on screw dislocations in bcc transition metals[END_REF] and leads to different estimates of 𝐸 𝑘𝑝 depending on the stress/temperature domain where experiments are conducted and on the analytical treatment use to extract 𝐸 𝑘𝑝 . To avoid these issues,[START_REF] Dezerald | First-principles prediction of kink-pair activation enthalpy on screw dislocations in bcc transition metals[END_REF] retained the lowest stress-high temperature values from[START_REF] Brunner | Temperature and Strain-Rate Dependence of the Tensile Flow Stress of High-Purity 𝛼-Iron below 250 K. I. Stress/Temperature Regime III[END_REF] as a reference: 0.83 eV, even though it gives an underestimation of 𝐸 𝑘𝑝 . Two recent potentials for Fe trained on ab initio data using machine learning methods yield 𝐸 𝑘𝑝 = 0.77 eV and 0.84 eV[START_REF] Goryaeva | Efficient and transferable machine learning potentials for the simulation of crystal defects in bcc Fe and W[END_REF]. The value obtained with Proville potential is therefore underestimated compared to other methods, with the effect of underestimating the lattice resistance to dislocation motion, and to possibly overestimate the relative effect of other concurrent thermally activated processes (e.g. interactions with defects).

Figure 3 . 3 :

 33 Figure 3.3: Kink pair nucleation enthalpy. (a) Minimum enthalpy path associated with the nucleation of a kink-pair, at two different levels of stress. (b) kink pair nucleation enthalpy as a function of the applied stress. All curves were obtained with Proville potential. Reproduced from Proville et al. (2012) (Supp. Materials).

Figure 3 . 4 :

 34 Figure 3.4: Shape of the transition state for kink pair formation as a function of the applied stress, with Proville interatomic potential. The equilibrium position of the initial straight dislocation 𝑔 𝑒𝑞 (𝜏 ) is indicated on the lower left. 𝑧/𝑏 is the position along the [111] direction normalized by the burgers vector length and 𝑥/𝑎 the coordinate in the glide direction normalized by the distance between Peierls valleys. The stress 𝜏 is normalized by the Peierls stress 𝜏 𝑃 .Reproduced from[START_REF] Ghafarollahi | Theory of double-kink nucleation in dilute BCC alloys[END_REF].

Figure 3 . 5 :

 35 Figure 3.5: Typical stress-strain curves from MD simulations, at different temperatures. Simulations conducted with Proville potential with an imposed shear rate γ = 1 × 10 7 s -1 . Stress peaks are represented by circles. A dashed line represents the average of the peaks for each temperature (i.e. the computed flow stress).

Fig. 3 Figure 3 . 6 :

 336 Fig. 3.6, we present values of 𝜏 𝑏𝑑𝐿 as a function of the segment length 𝐿, showing that for long segments, the plastic work remains dominant even for stress levels as low as 5-10 MPa, making reverse jumps highly unlikely. Experimental flow stress measurements conducted at high-temperature typically feature long dislocation segments (𝐿 ∼ 1 µm), and low flow stress values, corresponding to a domain where Δ𝐺 𝑘𝑝 is reduced due to the temperature, and possibly negligible compared to 𝜏 𝑏𝑑𝐿. An estimation of experimental dislocation density from the work of Brunner and Diehl (1991) yields an estimate of 𝐿 ≈ 1 √ 𝜌 = 10 3 to 10 4 Burgers vectors. In our MD simulations in the other hand, where 𝐿 = 40 b, the plastic work is much lower and backward jumps are more likely to happen under low stresses.

Figure 3 . 7 :

 37 Figure 3.7: Flow stress predicted using the different models of dislocation glide, as a function of (a) temperature and (b) the segment length (L). In (a), 𝐿 = 40 b (Γ = 10 4 ). In (b), the same conditions are simulated, but the value of 𝐿 varies.

  Fig. 3.7 (a) that the difference between models is limited to only very low stresses, allowing to define 𝑇 𝐶 by extrapolation, yielding the same 𝑇 𝐶 as the forward model. As the two thermally activated models show no difference except for very low stress values, we do not use the simple forward jump model in the rest of this chapter, and refer to the forward-and-backward jump model as "thermally activated model". With the model of Po et al., the velocity is proportional to a kink velocity term: 𝓋 ∝ 𝓋 𝑘 ∝ 𝜏

Figure 3 . 8 :

 38 Figure 3.8: Comparison of the terms 𝐿 and 𝜆 = 2𝑎𝑒 𝛽Δ𝐺 𝑘𝑝 /2 as a function of Δ𝐺 𝑘𝑝 . The term 𝐿 is dominant in the upper half of the plot, i.e. for low values of Δ𝐺 𝑘𝑝 . Orange circles correspond to an estimation of 𝐿/𝜆 in typical MD simulations conditions ( 𝜌 = 4 × 10 15 m -2, 𝐿 = 40𝑏, γ = 1 × 10 7 s -1 ), while red dots correspond to typical experimental conditions (𝜌 = 1 × 10 12 m -2 , 𝐿 = 1 × 10 -6 m, γ = 8.5 × 10 -4 s -1 ).

  Figure 3.9: MD flow stress data compared to model predictions assuming either a constant entropy or a Meyer Neldel law. Each subplot presents the predictions of a model (thermally activated on the left, Po et al. on the right), for two different simulation conditions (continuous curves), and compared to corresponding MD simulations (symbols, same color corresponds to same conditions). When assuming a constant entropy, Δ𝑆 0 is set to 8.5𝑘 𝐵 . For the MN law Δ𝑆 0 = 6𝑘 𝐵 , 𝑇 𝑀𝑁 = 500 K.

Figure 3 . 10 :

 310 Figure 3.10: Gibbs activation energy at low temperature, at 200 and 500 MPa. PAFI calculations were conducted with a high level of accuracy, resulting in a narrow uncertainty domain (gray area). The results of harmonic entropy calculations conducted at 200 and 500 MPa are presented as dashed lines.

Figure 3

 3 Figure 3.11: Minimum Gibbs free energy paths associated with the nucleation of a kink pair under 200 MPa in Fe, at different temperatures. (a) MGEP profiles: PAFI free energy data (symbols) and corresponding splines (lines). (b) Temperature evolution of the Gibbs energy of activation. A linear regression in the domain below 100 K is shown as a dashed line. Both plots share the same 𝑦 axis.

Figure 3 . 13 :

 313 Figure 3.13: Effective activation entropy Δ𝑆 as a function of (a) the activation enthalpy Δ𝐻 𝑘𝑝 and (b) the shear stress 𝜏. The slopes of linear regressions of Fig. 3.12 are plotted as a function of the activation enthalpy in (a) and stress in (b) (square symbols), and fitted to: (a) an inverse Meyer-Neldel entropy-enthalpy relation (Eq. 3.18), (b) a linear stress-activation entropy relation.

  3.3). At low temperature, Δ𝑆 can therefore be expressed either as Δ𝑆 = -Δ𝐻 𝑘𝑝 𝑇 𝑀𝑁 + Δ𝑆 0 with 𝑇 𝑀𝑁 positive, corresponding to an inverse Meyer Neldel behavior, or as a function of the stress Δ𝑆 = 𝛼𝜏 + Δ𝑆 0 . The parameterization of both equations is presented in Fig. 3.13. Interestingly, the inverse Meyer Neldel behavior seen at low temperature corresponds to a reinforcement effect between enthalpy and entropy, in contrast to the classical compensation effect.

Figure 3 .

 3 Figure 3.14: MD flow stress data compared to models predictions using an inverse Meyer Neldel law. Each subplot presents the predictions of a model (thermally activated on the left, Po et al. on the right), for two different simulation conditions (continuous curves), and compared to corresponding MD simulations (symbols, same color corresponds to same conditions).

Figure 3 . 15 :

 315 Figure 3.15: Gibbs free energy of nucleation of a kink pair as a function of stress.PAFI data (symbols) is fitted to Eq. 3.24, giving a set of parameters given in Tab 3.1. Each temperature data is represented in a different subplot, with the corresponding fit shown as a dark solid line. All curves are reported in the background of other sublots as light gray lines for easier comparison.

Figure 3 .

 3 Figure 3.16: MD flow stress data compared to models predictions using our Gibbs free energy calculations. Each subplot presents the predictions of a model (thermally activated on the left, Po et al. on the right), for two different simulation conditions (continuous curves), and compared to corresponding MD simulations (symbols, same color corresponds to same conditions).

Figure 3 . 18 :

 318 Figure3.18: Sketch of the distribution of knots along a free energy barrier at different temperatures. As the temperature increases, the MFEP is more and more deviated from the MEP and the difference between successive knots can become large.

  Fig. 4.1). Cottrell atmosphere formation kinetics and mobility are therefore important to predict the mechanical properties of aged steel and the plastic flow in the serrated flow regime. Atomistic simulations of their formation using atomistic Monte Carlo approaches were conducted by Waseda et al. (2017) and Candela et al. (2020) (Fig. 4.2), providing insights on the filling of lattice sites near the dislocation. Molecular dynamics simulations of decorated dislocations expose a powerful pinning effect of carbon, which increases the unpinning

  these approaches is the dislocation-solute binding energy which depends on the interatomic potential used. Predictions of different interatomic potentials are compared to ab initio values in Sec. 2.3.1.

Figure 4 . 1 :

 41 Figure 4.1: Snapshots from KMC simulations of a screw dislocation in tungsten interacting with interstitial oxygen atoms, showing different dynamical regimes: (a) with fast solute diffusion, the dislocation and solutes motion are uncorrelated, (b) for comparable solutes and dislocation velocities, solutes are able to segregate to the core and pin the dislocation before it can unpin and move further, (c) when solutes diffusion is slow, the dislocation moves in a field of static obstacles.Reproduced from[START_REF] Zhao | Simulating the mechanisms of serrated flow in interstitial alloys with atomic resolution over diffusive timescales[END_REF].

Figure 4 . 2 :

 42 Figure 4.2: Initial (a) and final (b) configurations of a mixed-lattice kinetic Monte Carlo simulation containing a screw dislocation (the core is represented by three black dots) and 25 C atoms (spheres) at 300 K. Black arrows in (a) represent the displacement of atoms between (a) and (b). Reproduced from (Candela et al., 2020).

Figure 4 . 3 :

 43 Figure 4.3: Softening and hardening effects due to carbon addition.Figure reproduced from Caillard (2011).

  Figure 4.3: Softening and hardening effects due to carbon addition.Figure reproduced from Caillard (2011).

iron

  [START_REF] Tapasa | Computer simulation of the interaction of carbon atoms with self-interstitial clusters in 𝛼-iron[END_REF];[START_REF] Veiga | Comparison of atomistic and elasticity approaches for carbon diffusion near line defects in 𝛼-iron[END_REF];[START_REF] Waseda | Formation of carbon Cottrell atmospheres and their effect on the stress field around an edge dislocation[END_REF];[START_REF] Zamzamian | A study on the mobility of 12111{01¯1} edge dislocation in low-carbon 𝛼-Fe and its interactions with damage cascade: On picosecond time scale using molecular dynamics simulations[END_REF][START_REF] Ishii | Conjugate Channeling" Effect in Dislocation Core Diffusion: Carbon Transport in Dislocated BCC Iron[END_REF], but simulating its effect on screw dislocations remains challenging. Screw dislocation locking by carbon was simulated usingMD Veiga et al. (2015);[START_REF] Khater | Dislocation glide in Fe-carbon solid solution: From atomistic to continuum level description[END_REF] and energetic computations[START_REF] Becquart | Atomistic modeling of an Fe system with a small concentration of C[END_REF];[START_REF] Clouet | Dislocation interaction with C in 𝛼 Fe: A comparison between atomic simulations and elasticity theory[END_REF];[START_REF] Hanlumyuang | Interactions between carbon solutes and dislocations in bcc iron[END_REF]

  2, it becomes possible to simulate by MD the glide of screw dislocations in random solid solutions of C atoms, at low temperature. Simulations of a full solid solution are used to assess the importance of short-range core interactions compared to long-range elastic ones. A detailed analysis of an elementary pinning configuration, where a screw dislocation interacts with C atoms of varying separation is also performed. For this, we use a combination of MD simulations, saddle-point search methods and a statistical model of the unpinning stress based on the harmonic transition state theory (HTST).

4. 2

 2 Dislocation mobility in a solid solutionWe start by simulating the glide of a screw dislocation interacting with a random solid solution of C atoms. In order to observe several pinning-unpinning events in the course of a MD simulation, and the combined effect of several C atoms along the dislocation, we set the concentration to 0.52 at%, i.e. 1000 solute atoms in a box containing 192000 Fe atoms. It was created by inserting carbon atoms in octahedral sites randomly selected in the PAD simulation cell. Considering the short simulated time (a few nanoseconds) and the investigated temperatures, carbon atoms are not able to diffuse and are considered immobile in the rest of this chapter.

Figure 4 . 4 :

 44 Figure 4.4: Dislocation line (green) adopting a complex three-dimensional kinked structure in presence of carbon atoms (red spheres). Only the solutes closest to the dislocation are represented for clarity. The dislocation line was identified using the DXA algorithm (Stukowski et al., 2012).

Figure 4 . 5 :

 45 Figure 4.5: Constant shear rate simulations with a carbon-depleted zone around the glide plane. (a) Stress-strain curves during constant shear-rate MD simulations of random solid solutions at 300 K. Blue lines represent simulations in which a carbon depleted zone of 𝑛 𝑝𝑙𝑎𝑛𝑒𝑠 atomic planes above and below the glide plane was created, as illustrated in (b). The carbon concentration is 0.52 at when 𝑛 𝑝𝑙𝑎𝑛𝑒𝑠 = 0. The gray line in (a) represents the same deformation test conducted in pure Fe.

Figure 4 . 6 :

 46 Figure 4.6: Local conversion of the screw dislocation core during a dynamical simulation. The smaller red atom is a carbon atom, while all others are iron atoms. Iron atoms that adopt a distorted prismatic structure near the carbon atom, are highlighted in blue. The dislocation line obtained using the DXA algorithm is shown in green. For clarity, only core atoms are shown.

  2.2), aligned in the[111] direction inside the dislocation glide plane to avoid cross-slip. Due to the periodic boundary conditions in the[111] direction, the carbon atoms form an infinite one-dimensional array with a spacing between carbon atoms defined by two lengths: ℓ and 𝐿, with ℓ ≤ 𝐿 as illustrated in Fig.4.7 (a).

Fig. 4 .Figure 4 . 7 :Figure 4 . 8 :

 44748 Fig. 4.7 details the full unpinning processes of a dislocation, with the corresponding stressstrain curve, simulated with ℓ = 10 𝑏, 𝐿 = 30 𝑏 at 40 K. A step-by-step visualisation of the MD simulation can be decomposed in 6 successive images shown in Fig. 4.7 (a). The dislocation is initially in the Peierls valley before the C atoms (step 1). It transits to the valley just after the C atoms (step 2) at a relatively low stress of about 250 MPa, which is marked by a slight stress drop highlighted by a dashed circle in Fig. 4.7 (b). The magnitude of the stress drop is Δ𝜏 = 𝜇𝜌𝑏𝑑 (𝑑 is the distance Peierls valleys) and correspond to the plastic work associated with the first jump. This configuration is the true pinning configuration, which resists dislocation motion until a high applied stress of about 950 MPa. At this point, a kink

Figure 4 . 9 :

 49 Figure 4.9: Dependence of the average unpinning stress on the largest dislocation segment length 𝐿. MD at (a) 100 K and (b) 300 K. Blue circles are when ℓ < 𝐿, red squares when ℓ = 𝐿. Half-filled orange symbols refer to the flow stress in pure iron computed by MD in cells of different lengths 𝐿. The dashed and solid lines are predictions of a statistical model (see text for details).

Fig. 4

 4 Fig. 4.7 (a)), followed with an ART convergence to identify precisely the saddle configurations. We consider here the case 𝐿 = ℓ where only one C atom needs to be inserted in the simulation cell. The reference state for the NEB calculations, 𝑆 0 , is constructed in the same way as in the dynamic simulations, by inserting a C atom in an 𝑂 (2) site, and relaxing the system. Final states are constructed by introducing the dislocation in different Peierls valleys in the glide direction. The initial NEB path was obtained by linear interpolation between the initial and final state. Equilibrium and activated states for a dislocation with 𝐿 = 30 𝑏 at an applied stress of 600 MPa are shown in Fig. 4.10 (a), with the corresponding energy barriers in Fig. 4.10 (b).

FeFigure 4

 4 Figure 4.10: NEB calculations of the two-step unpinning process of a dislocation from a row of C atoms. (a) Successive metastable and activated states during the interaction process. Atoms of the core are identified using the centrosymmetry parameter (Kelchner et al., 1998). (b) Sketch of the corresponding energy profile. (c) and (d) show the activation enthalpies of both transitions (circles) for two different solute separations. The solid lines are fits using Kocks' law (Eq. 4.9). (e) and (f) show the fitted parameters Δ𝐻 0,𝑖 and 𝜏 𝑃,𝑖 of Kocks' law. Solid lines in (e) and (f) represent the average of Δ𝐻 0 and a fit of 𝜏 𝑃 using Eq. 4.10 for both transitions. Dashed lines in (c) to (f) represent the corresponding quantities in pure iron.

  𝑖 (𝜏 , 𝐿) = Δ𝐻 0, 𝑖 (𝐿) (sufficient to set the exponents 𝑝 = 0.6 and 𝑞 = 1.1 to the values obtained by fitting the kink pair formation enthalpy in iron with Proville potential and to only adjust Δ𝐻 0 and 𝜏 𝑃 . The fits are shown as continuous lines inFig. 4.10 (c, d) while the values of Δ𝐻 0,𝑖 (𝐿) and 𝜏 𝑃,𝑖 (𝐿) are reported in Fig.4.10 (e, f). They represent the effect of the solute spacing on the activation enthalpy. We found Δ𝐻 0,1 ≈ Δ𝐻 0,2 , independent of 𝐿 and only slightly larger than the value in pure iron (see Fig.4.10 (e)). For the sake of simplicity, we used their average value, Δ𝐻 0,𝑖 , in the model. The effective Peierls stresses 𝜏 𝑃,𝑖 on the other hand strongly depends on the solute spacing and are different for the first and second barrier.As seen in Fig.4.10(f), they can be accurately fitted as:𝜏 𝑃,𝑖 (𝐿) = 𝑎 𝑖 (𝐿/𝑏) 𝑏 𝑖 + 𝑐 𝑖 (4.10) with 𝑎 1 = 7.89 × 10 5 MPa, 𝑏 1 = 3.0, 𝑐 1 = 959 MPa, 𝑎 2 = 7.10 × 10 5 MPa, 𝑏 2 = 2.7, 𝑐 2 = 1.01 × 10 3 MPa. The strong dependence of Δ𝐻 2 on L seen in Fig.4.10 (c) and (d) implies that short segments have a much reduced nucleation probability compared to larger segments, resulting in the negligible effect of ℓ observed in Sec. 4.3.1.

  models of dislocation pinning as developed in[START_REF] Zhao | Direct prediction of the solute softening-to-hardening transition in W-Re alloys using stochastic simulations of screw dislocation motion[END_REF];[START_REF] Zhao | Simulating the mechanisms of serrated flow in interstitial alloys with atomic resolution over diffusive timescales[END_REF];[START_REF] Shinzato | An atomistically informed kinetic Monte Carlo model for predicting solid solution strengthening of body-centered cubic alloys[END_REF] is the main perspective of this work to better understand the effect of carbon solutes. Such a model would allow the study of the complex and collective interplay between solute-dislocation interactions and stress-mediated solute diffusion leading to the formation of Cottrell atmospheres, over length-and time-scales unaccessible to MD simulations, more accurately than approaches fully based on harmonic theories.

Figure 1 :

 1 Figure 1: Flowchart of multiscale modelling of plasticity adapted from Tipping (2010).The contribution of the present work is emphasized (in green), as well as its main perspectives (in blue).

  Methods

	Direct interpolation		Kink pair construction
	r i Final	[111]	Final	i/R⨉L	r i
	Initial		Initial	l 0
	Figure 1.3: Creation of an initial "guessed" NEB path that explicitly contains a kink pair. Direct interpolation (left) is compared to our approach (right). With direct interpolation, atoms positions 𝑋 are linearly interpolated between the initial and final state. In our approach, for replica 𝑟 𝑖 (with 1 ≤ 𝑖 ≤ 𝑅 -1), a zone of length 𝑖/𝑅 × 𝐿 in [111] direction has the same atoms coordinate as the final state. A buffer zone of size 𝑙 0 = 5 Å makes a linear transition to the rest of the cell, where atoms positions are those of the initial state.
	reffered to as the classical NEB method, necessitates a lot of iterations to reach an accurate
	estimate of the energy of a state of interest, which is usually the activated state (i.e. the local
	energy maximum). As forces are evaluated on all replicas at each iteration, this results in a
	high computational workload to obtain the value of this state of interest. This issue is ad-
	dressed in the climbing-image NEB method (CI-NEB), available in LAMMPS, which allows
	to converge more efficiently to the saddle. After a limited number of classical NEB itera-
	tions, the replica with the highest energy is selected as the best activated state candidate. It
	is then forced to climb to the local maximum by reversing the component of the force in the
	direction of the path tangent. The climbing force is thus:

Table 2

 2

	Notation

.1: List of interatomic potentials and ab initio data used in this study.

  Degenerate core structure induced by C, with the Hybrid potential. Differential displacements of atoms in the dislocation core region after relaxation with an array of C atoms inserted in the 𝑂 (2) every 𝑏.features a metastable intermediate split-core structure[START_REF] Gordon | Screw dislocation mobility in BCC metals: the role of the compact core on double-kink nucleation[END_REF] associated with a bent trajectory between easy positions, as shown in Fig.2.6. Due to this artefact, glide

	2.2 Peierls barrier
	Figure 2.5:
	in agreement with DFT. Though, it

occurs by the nucleation of half-kinks with a low barrier, at odds with ab initio data both qualitatively and quantitatively. Attempts to artificially enforce the nucleation of full kinkpairs in NEB calculations of

[START_REF] Wang | Atomistic simulations of carbon effect on kink-pair energetics of bcc iron screw dislocations[END_REF]

using the Becquart potential still resulted in ill-shaped barriers, with a metastable intermediate state. The potential developed by Gordon et al. (2011) aims at correcting the barrier shape, but still has a metastable split core (though less deep), and a barrier height close to 12 meV/𝑏. A recent potential based on the Angular-Dependent Potential (ADP) formalism

[START_REF] Mishin | Phase stability in the Fe-Ni system: Investigation by first-principles calculations and atomistic simulations[END_REF]

, fitted with the help of machine learning techniques

[START_REF] Mishin | Machine-learning interatomic potentials for materials science[END_REF]

by

Starikov et al. (

  are to be avoided, especially to study the interaction of C and dislocation cores. Values of 𝐸 𝑚 for potentials available in the litterature are presented in Fig.2.9. The potentials of Liyanage and Henriksson strongly underestimate this value, resulting in a large overestimation of diffusion kinetics.

	Hybrid				0.95
	DFT				0.90
	Hepburn				0.89
	Becquart				0.85
	Lee				0.82
	Henriksson		0.32
	Liyanage	0.21
	0.0	0.2	0.4	0.6	0.8
			Energy (eV)

Figure 2.9: Migration barrier for carbon between octahedral sites. Data extracted from Refs.

  3.2.3, i.e. either a constant entropy, or a Meyer-Neldel law. Models presented side-by-side use the same expression of Δ𝐺 𝑘𝑝 . It can be seen that assuming a constant entropy Δ𝑆 0 = 7𝑘 𝐵 yields good results at low temperature. The thermally activated model gives a discrepancy at higher temperature, which is larger with the model of Po et al. that predicts a slower decrease of the flow stress.

  This would in addition allow to use computing nodes that are not designed for distributed parallel calculations, i.e. which are not connected by a high-performance network such as an Infiniband network or do not have a MPI library installed. Workflow managers (see e.g.[START_REF] Rocklin | Dask: Parallel Computation with Blocked algorithms and Task Scheduling[END_REF];[START_REF] Huber | AiiDA 1.0, a scalable computational infrastructure for automated reproducible workflows and data provenance[END_REF];[START_REF] Köster | Sustainable data analysis with Snakemake[END_REF]) could be used to automatically submit and track jobs on HPC platfoms with minimum code addings to PAFI.

  The dashed and solid lines are predictions of a statistical model (see text for details). When used with correct values of Δ𝐺 𝑘𝑝 , this simple model showed very good agreement with MD data in Chap. 3. One caveat of this model is its low-stress limit, which does not prevent the flow stress from going to zero as neither reverse jump or kink migration is taken into account. However, this issue is irrelevant in the highstress, low-temperature domain considered in the following, where this model yields good results.

	𝓋 = 𝜈𝐿 exp ( -	Δ𝐺 𝑘𝑝 (𝜏 * , 𝑡) 𝑘 𝐵 𝑇	)	(4.1)
	which depends only logarithmically on 𝐿.	(	𝜈𝜌𝑏𝐿 γ ),	(4.2)

where 𝜈 is an attempt frequency and Δ𝐺 𝑘𝑝 is the Gibbs activation energy. The flow stress 𝜏 * is thus solution of Δ𝐺 𝑘𝑝 (𝜏 * , 𝑡) = 𝑘 𝐵 𝑇 log

Available at http://emmanuel.clouet.free.fr/Programs/Babel/index.html

When only the initial state is known, single-ended techniques such as the Activation Relaxation Technique presented in this section can be used to sample the accessible saddle states.

See the relevant pull-request on LAMMPS public development repository: https://github.com/lammps/ lammps/pull/2486#pullrequestreview-540994609

For a discussion on uncertainty in PAFI, see: https://github.com/tomswinburne/pafi/blob/master/ error_analysis.pdf

Online access at: https://github.com/tomswinburne/pafi. The code version used in this work and presented in this section is available at: https://github.com/tomswinburne/pafi/tree/fa82b82c73.

Online documentation at: https://docs.lammps.org/fix_pafi.html

See https://github.com/arn-all/phd-thesis for the Python implementation used for our calcula-

As EAM potentials take interatomic distances as their only input, this overlapping brings the difficulty of accurately reproducing different environments which have similar interatomic distances with the same potential. This issue is partly addressed in the case of MEAM potentials, which also include an angular component[START_REF] Baskes | Modified embedded-atom potentials for cubic materials and impurities[END_REF].

As each row is at the intersection of three different triangles, each split position exists in three different variants.

This type of smoothing filter uses a least-squares fits to polynomials on a moving window along the curve to smooth its variations while preserving the height of peaks better than a simple rolling average.

Peaks are detected using the scipy.signal.find_peaks Python module with a prominence set to 15. Increasing this parameter has the effect of ignoring smaller peaks.

This is precisely the case encountered at the end of this chapter, when testing glide models on typical experimental conditions, based on Gibbs energies computed with the Proville potential.

In more concentrated alloys, other phases can appear depending on their stability as predicted by the phase diagram. Here, we focus on low-carbon ferritic steel.

This effect is observed even though the annihilation of two opposite kinks would lead to a large energy gain, which can be estimated from the drop of enthalpy in the last stage of the Peierls barrier under stress (for 0.6 GPa the Fig.

3.3 (a)), and which is likely higher than the energy of the carbon in the dislocation. This result suggests that there is a high energy barrier for the kinks to cross the carbon, which is larger than the barrier of nucleation on top of the nucleated kink pair.

temperatures below 100 K, giving an approximation of Δ𝑆 = -𝜕𝐺 𝑘𝑝 /𝜕𝑇 assimilated to a constant. This close-to-linear domain is alike a harmonic regime, with the difference that it does not directly derive from the vibrational modes around the initial and activated states. We also emphasise how the small deviations from linearity at low temperatures affect the computed activation entropy, making such calculations very challenging. To be able to measure these small variations in Fig. 3.10, we improved the precision of the Gibbs energy calculations by increasing the number of independent samplings, resulting in a narrow uncertainty domain (represented by a light grey region in Fig. 3.12) that allows to capture all non-linearities, but with a high computational load. The linearization that we propose can thus also be seen as a pragmatic approach to save computational resources by neglecting small non-linearities. In Fig. 3.12, we present values of Δ𝐺 𝑘𝑝 computed over a large range of stresses and temperatures. As expected, Δ𝐺 𝑘𝑝 decreases when the temperature or the applied stress increase.

Also note that the sampling error is larger than in Fig 3 .10, as less sampling repetitions were realized to save computation time and explore more stress/temperature conditions. For all tested values of the stress, the variations of Δ𝐺 𝑘𝑝 remain close to linear below 100 K, as highlighted by the blue dotted lines. The effective activation entropy Δ𝑆 𝑟𝑒𝑔 obtained from between the 0 K extrapolation of experimental flow stress data, which is close to 0.4 GPa in iron [START_REF] Kubamoto | Thermally activated slip deformation between 0.7 and 77 K in high-purity iron single crystals[END_REF], and predictions of classical TST models, which are close to 1 GPa (Fig. 3.3). Several authors proposed to correct this discrepancy either phenomenologically [START_REF] Domain | Simulation of Screw Dislocation Motion in Iron by Molecular Dynamics Simulations[END_REF] or by accounting for the effect of dislocations junctions [START_REF] Bulatov | Nodal Effects in Dislocation Mobility[END_REF], collective dislocation dynamics effects [START_REF] Gröger | Explanation of the discrepancy between the measured and atomistically calculated yield stresses in body-centred cubic metals[END_REF] or quantum effects [START_REF] Proville | Quantum effect on thermally activated glide of dislocations[END_REF]. While these processes are of great importance, their integration is beyond the scope of the present work. In order to compare relative variations, we plot experimental and simulated data on separate scales [START_REF] Clouet | Screw dislocations in BCC transition metals: From ab initio modeling to yield criterion[END_REF]. Another caveat of comparison to experimental data is the estimation of the segment length 𝐿 and dislocation density 𝜌 [START_REF] Brunner | Temperature and Strain-Rate Dependence of the Tensile Flow Stress of High-Purity 𝛼-Iron below 250 K. I. Stress/Temperature Regime III[END_REF]. Here, we take 𝐿 = 1 µm, 𝜌 = 1 × 10 12 m -1 , γ = 8.5 × 10 -4 s -1 and use value of 𝜈 = 5 × 10 10 s -1 used previously.

In Fig 3.17, we compare experimental flow stress data to the predictions of both dislocation glide models, based on our parameterization of Eq. 3.24. Both models yield similar predictions in this case, as we have 𝐿/𝜆 ∼ 10 -2 , i.e. glide is limited by kink pair nucleation. It can be seen that this parameterization allows the models to reproduce well the overall relative evolution of the experimental flow stress. The athermal temperature 𝑇 𝐶 which marks the end of the thermally activated regime is well defined with both models with this set of conditions, even though the definition of 𝑇 𝐶 has been discussed in this chapter (see Sec. 3.2.4). The difference between model predictions and experiments is mainly due to an underestimation of the athermal temperature, leading to a discrepancy that gets larger with increasing temperature.

While this result is encouraging for the integration of entropic and anharmonic effects in kinetic models of dislocation glide, a more realistic comparison to experimental data would benefit from the integration of the low-temperature effects mentioned above, as well as the use of a more realistic interatomic potential. In addition, neither our implementation of the thermally activated model nor of the model of Po et al. include non-schmid effects which are well-known to play a role in experimental conditions [START_REF] Spitzig | The effect of orientation and temperature on the plastic flow properties of iron single crystals[END_REF][START_REF] Hale | Insights on activation enthalpy for non-Schmid slip in body-centered cubic metals[END_REF][START_REF] Stukowski | Thermally-activated non-Schmid glide of screw dislocations in W using atomistically-informed kinetic Monte Carlo simulations[END_REF][START_REF] Kraych | Nonglide effects and dislocation core fields in BCC metals[END_REF][START_REF] Clouet | Screw dislocations in BCC transition metals: From ab initio modeling to yield criterion[END_REF].

Conclusions

In this chapter, we discussed different dislocation velocity models and tested their validity using MD simulations. The tested models notably differ by their complexity, ranging from a simple forward-jump model to the more comprehensive approach of [START_REF] Po | A phenomenological dislocation mobility law for bcc metals[END_REF] which takes into account the kink migration velocity and the equilibrium density of kink pairs along short segments (i.e. a few Burgers vectors), the available length is small enough for the line to advances while remaining straight instead of forming a double kink. For this process, a much reduced activation entropy is expected compared to the kink pair mechanism, as (i) both the initial and activated states have a unique position due to the invariance in [111] direction, leading to a null configuational entropy and (ii) the vibration modes of two straight segments are likely to be less different than those of a straight line and a kink pair. For longer segments on the other hand, the behaviour of pure iron is recovered, and Δ𝑆 ∼ 𝛼 𝐹 𝑒 𝑘 𝐵 𝜏 at low temperature, which is independent of L. In order to assess the possible effect of an activation entropy which increases with 𝐿 and smoothly tends to the value of pure iron, we tested a simple phenomenological law 𝛼(𝐿) = 𝛼 𝐹 𝑒 tanh ( 𝐿 100𝑏 ) for both barriers. This is again a simplification since we can expect the activation entropy to vary differently for both barriers. However, we can see in Fig. 4.9 where the prediction of the model with the length-dependent activation entropy is shown as a continuous line, that the simple phenomenological law greatly improves the agreement between the statistical model and the MD simulations. Fig. 4.11 allows to compare the prediction of the two models plotted on the same scale. This suggests that the activation entropy is reduced for closely-spaced defects, which is left to investigate in future work.