Dr Mark Gilbert

Alexander Zachariah

Mrs Suja

grandparents (Mr Alexander Kurian

Ms Mr Antony Edith Dinu

Ms Isabelle Molins

Mr Bourigaud

Rajan Evin

Keywords: Métaux, paysage énergétique, défauts, potentiels empirique, ab initio Metals, energy landscape, defects, empirical potentials, ab initio

Exploration du paysage énergétique de défauts dans les métaux cubiques centrés

Résumé

Les matériaux composants les réacteurs nucléaires subissent des conditions d'irradiation sévères, donnant lieu à des modifications de leurs propriétés mécaniques. Le vieillissement de ces matériaux en condition d'utilisation soulève des questions aussi importantes que celles liées à la sécurité des centrales existantes dans le cadre de l'allongement de leur durée de vie et celles liées aux systèmes énergétiques du futur tels que les réacteurs à fission de génération IV et le confinement de la fusion dans les tokamaks. Les réacteurs fonctionnant sur le principe de la fission nucléaire utilisent essentiellement les alliages de fer comme matériaux de structure. Les systèmes énergétiques futurs reposant sur la fusion nucléaire proposent d'utiliser les alliages à base tungstène. Le vanadium et le tantale sont également des matériaux envisagés pour des réacteurs à fusion. Le point commun entre ces différents matériaux est leur structure cristalline de symétrie cubique centrée CC.

Les collisions entre les particules irradiantes et les atomes formant le réseau cristallin des matériaux engendrent des défauts ponctuels dont la migration par activation thermique mène au regroupement de ces défauts sous forme d'amas responsables du vieillissement de matériaux. Mieux comprendre comment ces amas se forment est donc essentiel pour appréhender l'évolution des matériaux sous irradiation. Dans cette thèse nous avons étudié les propriétés du paysage énergétique de défauts ponctuels dans les métaux cubiques centrés fer, tungstène, vanadium et tantale. En s'appuyant sur la modélisation à l'échelle atomique, nous avons souhaité améliorer et réduire le caractère empirique des modèles phénoménologiques existants.

La modélisation des défauts ponctuels à l'échelle atomique peut être réalisée avec différentes méthodes se différenciant uniquement par la qualité de la description de l'interaction entre atomes. Les études utilisant des interactions atomiques exactes, type ab initio, nécessitent des calculs lourds, réduisant drastiquement les possibilités d'investigations systématiques et en particulier rendant impossible l'étude directe des amas de grandes tailles. A mesure que l'on réduit la qualité de la modélisation de l'interaction atomique, on réduit la fiabilité et le caractère prédictif du calcul. C'est le cas avec la modélisation des interactions atomiques via les potentiels semi-empiriques. Ceux-ci permettent toutefois de réaliser des études systématiques des amas en fonction de leur taille. Dans cette thèse nous avons pu dépasser le facteur limitant des calculs ab initio et éviter l'utilisation des potentiels semi-empiriques. Nous avons développé un modèle énergétique original pour les boucles de dislocation ainsi que pour les amas interstitiels tridimensionnels de type C15. Le modèle obtenu est sans limite de taille et peut être paramétré entièrement par les calculs ab initio. Afin de tester sa robustesse pour les grandes tailles d'amas nous avons également paramétré ce modèle par rapport à des calculs en potentiels semi-empiriques et comparé les prédictions du modèle aux simulations atomiques sur des amas de grandes tailles. Pour développer ce modèle,nous avons dû combiner le caractère discret des petit amas qui tient compte du voisinage chimique de chaque atome avec le caractère continu des grandes boucles de dislocation, traité suivant l'approximation de l'élasticité anisotrope. En utilisant cette approche il a été possible de prédire les énergies de formations des amas C15 et des boucles de dislocations pour toutes les tailles avec la précision des calculs ab initio. Les seules données nécessaires sont les énergies de formation ab initio d'une série de configurations de petits amas d'interstitiels dans les différents matériaux étudiés.

Grâce à notre développement nous avons pu déterminer : (i) la stabilité relative des boucles de dislocation d'interstitiels dans le cas du Fe d'après leur vecteur de Burgers. (ii) la stabilité des amas C15 par rapport aux amas de type boucle. Nous avons montré que les amas de type C15 étaient plus stables lorsqu'ils impliquent moins de 41 interstitiels dans le fer. Ce résultat est très important pour les expérimentateurs puisqu'il fournit une indication précise de la taille maximale de ces amas C15, qui correspond à environ 2 nm de diamètre dans Fe. (iii) dans le Ta nous avons pu mettre en évidence la même stabilité jusqu'à 20 interstitiels ouvrant une voie possible de détection expérimentale de ces amas correspondants à environ 1 nm de diamètre. (iv) pour l'ensemble des Right from the moment when I took the decision of pursuing my Ph.D., to this day when I am writing my note of Acknowledgement, I am indebted to a number of amazing people. People who have seen me go through ups and downs in this endeavour of mine and have stuck with me despite all of it; I owe this day to them and though, this is an insufficient response to all their goodness, I hope they may be able to sense my humble gratitude.

First and foremost, I would like to thank my supervisor, Dr. Mihai-Cosmin Marinica, for directing and supporting me throughout the course of my Ph.D. study. His constant guidance, patience, optimism and willingness to share his vast knowledge helped me immensely. From launching my first calculation on VASP to writing my thesis, he has always been willing to help me out. To say the least, I could not have imagined a better advisor and mentor for my Ph.D. study. I would also like to thank the director of my thesis, Dr. Laurent Proville, for his insightful advice, enriching comments and encouragement which provided me with valuable guidance. It has indeed been a pleasure to have worked with both of you.

I am also thankful to the jury: Dr. Laurent Pizzagalli, Dr. Pär Olsson, Dr. Charlotte Becquart, Dr. Robin Schäublin and Dr. Steve Fitzgerald, for taking time out from their busy schedule to review my thesis. I am simply honored to have such an inspiring group of scientists as my jury.

I am indebted to Dr. François Williame, for giving me this opportunity of carrying out research in the scientifically stimulating environment of SRMP, during my master's degree and my Ph.D. I would also like to extend my heartfelt gratitude to Dr. Jean-Luc Béchade, for his constant encouragement and support which have been a significant source of motivation for me.

Contents

I am also grateful to the scientists from SRMP, for the stimulating discussions, presentations and seminars that made me broaden my research perspective. Particularly, I thank Dr. Thomas Jourdan, Dr. Estelle Meslin, Dr. Anne Duchateau, Dr. Hélène Lefaix-Jeuland and fellow doctoral aspirant Sofi Gorondy Novak for EPIGRAPH-related discussions. I also thank, Dr. Fabien Bruneval, Dr. Emmanuel Clouet and Dr. Bernard Legrand for the useful seminars at SRMP. It has been an absolute privilege of mine to have undergone training at SRMP-CEA Saclay. I greatly appreciate the cultural and scientific diversity among the scholars that has indeed helped me grow as a person.

Introduction

Safety of a commercial nuclear fission plant is ensured by its four barriers: the nuclear fuel rod, the reactor cooling system, the containment building and the reactor building. Undoubtedly, nuclear power plants (both fission and planned fusion reactors) require a whole range of structural materials to maintain their integrity and to contain radiation under all circumstances. Since longevity and safe functioning of the nuclear power plant rely heavily on the adequacy of the structural materials, a thorough study of their behaviour under radiation damage is inevitable. Study of radiation damage in structural materials is thus central to the safety as well as smooth functioning of the nuclear power plants.

In numerous situations, the main materials of interest in the nuclear fission and fusion industries are the bcc transition metals, such as special steels or materials based on tungsten, tantalum and vanadium. According to Ref. (1 ), Molecular dynamics simulations and experimental studies have shown that the body centered cubic lattice demonstrates improved radiation resistance compared to the close-packed face centered cubic lattice due to reduced amount of vacancy and interstitial defect clustering that occurs directly within displacement cascades. This endorses continued development of ferritic/martensitic steels and vanadium alloys as promising candidates for fusion reactor structures. (1 ) Reduced activation Ferritic/Martensitic steels exhibit high strength at higher temperatures, as can be established from Refs. (2 -4 ) while vanadium alloys demonstrate a good combination of strength, ductility and radiation resistance (5 -7 ). Tungsten alloys are candidate materials for first wall and divertor components of future fusion reactors because of their high melting point and hence higher operating temperatures (8 -10 ). Tantalum is known for its high toughness, high-sputtering threshold energy, easy fabricability and low activation properties which make it a candidate for target material for spallation sources or even first wall material in fusion reactors (11 -14 ). As such, study Contents of defects in Fe, V, Ta and W serve as the basic foundation for future research regarding structural materials in the nuclear industry. It can also supplement the current industrial research and lead to better understanding of the behavior of potential structural materials.

To model the behaviour of materials in reactor theory needs multiscale approach. These multiscale models require some inputs like formation energies of the Self Interstitial Atom (SIA) defects. Although there exist many methods to determine the formation energy of SIAs, each has its limitation. Ab initio calculations, which will be discussed further in the chapters, incorporate some reasonable approximations and can be performed using standard ab initio simulation packages. Such calculations are based on Density functional theory and provide accurate formation energies for small-sized interstitial defects but are computationally costly at larger sizes (say, above 22 SIAs).

To address this size limitation, various empirical potentials have been developed for Fe and W with Embedded Atom Method (EAM) being the most popular. These EAM potentials have facilitated calculation of formation energy while being computationally feasible for large cluster sizes and fast. However, they do not yield accurate values for formation energy and are non-transferable, i.e. reliable description of material properties strongly depends on the relation of these properties with the parameters fitted. For instance, one important contribution of EAM potentials is allowing for testing scaling laws and theory but the formation energy predicted using different empirical potentials can have a difference of up to 400 eV for 1000 SIAs in Fe.

To summarize, existing accurate ab initio calculations to calculate formation energies are limited by size, calculations using EAM potentials are hindered by non-transferability and the easily accessible ad-hoc laws are empirical. That is what justifies the present study. My work aims to address the need for a reliable atomistic model that can provide reasonably accurate zero K formation energy or finite-temperature free energy without size-limitation and which can be applied with reasonable ease in further multi-scale studies. Although this study would be in principle applicable to all bcc metals, this thesis concentrates on interstitial defects in Fe, W, Ta and V.

The organization of chapters in this thesis is as follows.

The first chapter is dedicated to the literature survey of the theoretical background of radiation damage in the four transition metals of interest. It includes a brief summary of the radiation damage theory and a survey of the existing experimental and simulation-Contents based results. This chapter aims to emphasize the importance of this study and equip the reader with relevant background of the radiation damage in the concerned metals, as established before this study.

The second chapter presents a theoretical background of the existing atomistic techniques to calculate formation energies of defects in transition metals. The first section in this chapter traces the chronological evolution of the ab initio methods that are currently used to calculate accurate formation energies of 'limited number' (up to 22 selfinterstitial atoms) of interstitial defects. The second section involves description of the tight-binding approach while the third section is dedicated to 'approximate' empirical potentials with a special emphasis on EAM potentials for bcc metals.

The third chapter is devoted to the development of new EAM potentials for Fe. Firstly, existing empirical potentials and their shortcomings are discussed to highlight the need for a new empirical potential. Then, the lessons learned from the fitting strategies of existing empirical potentials are recapitulated to devise an effective strategy to develop new empirical potentials for studies involving irradiation-induced defects in Fe. This strategy is then implemented in the development of new empirical potentials with a description of the fitting procedure. Selection rules for the construction of stable C15 clusters are also presented. Finally, the relative stability of C15 clusters and dislocation loops is studied using existing and newly-developed empirical potentials.

The fourth chapter focusses on a systematic description of the steps involved in the development of ab initio accuracy model to predict formation energies of defect clusters with various geometries and sizes. The first section in this chapter recapitulates the elastic theory in the 'continuum limit' applicable for large two-dimensional (2D) loops. Then, the second section formulates integration of 'discrete' ab initio formation energies for small-sized defects with 'continuum' anisotropic elastic theory formulation for formation energy of large-sized defects. This new formulation provides a 'discretecontinuum' model to calculate formation energy of defects as a function of number of SIAs without size-limitation at 0K and no strain condition. The discrete part involves transformation of ab initio formation energies of small-sized defects into an energetic model based on geometry i.e. in terms of 1st and 2nd nearest neighbors. The third section outlines construction of a database of defect structures (2D and 3D) up to 22 SIAs, following some general rules from previous studies. The fourth section illustrates validation of the discrete-continuum model after tests using existing empirical potentials for Fe and W. After validation, the fifth section presents the ab initio based predictions Contents of the energy landscape of defects and also provides scaling laws for extrapolation which provide formation energies of defects (2D or 3D) of any size for the four metals solely as a function of number of interstitials. This chapter concludes with a summary of the main results obtained using the discrete-continuum model.

Chapter five is concerned with calculation of free energy of defects at finite temperatures. The first section emphasizes the need for free energy calculations. The second section explains extension of the discrete-continuum model, developed for 0K conditions, to calculate finite-temperature free energy. The temperature dependence of elastic constants and other variables of the four metals was derived from experimentally observed values of elastic constants and the corresponding variables at finite temperatures. This temperature dependence of elastic constants was included in the discrete-continuum model formulation for formation energy at 0K (explained in chapter 4) in order to develop an atomistic model with validity at finite temperatures. The third section supplies another possible method of free energy calculations involving the quasi-harmonic approximation. Then, the free energy calculations from the two models are compared to validate the finite-temperature discrete-continuum model. Two applications of this finite-temperature discrete-continuum model are also presented. The first application involves explanation of the mechanism of formation of 100 loops using the present finite temperature model while the second one involves use of this model in cluster dynamics for comparison with capillary law experimental results.

State of the Art Review

The practical realization of new and improved engineering technologies in extreme conditions, be it for nuclear reactors, medical application or space exploration, is directly related to the sustainability of available materials in the specified conditions. Even when a material sustains reasonably well in these specified extreme conditions, lifetime of such an installation is most often limited by the material in use. Thus, a thorough study of the sustainability and lifetime of materials are essential for any future application of novel engineering technologies that push the boundaries of extreme conditions. This study includes consideration of some possible changes and processes due to irradiation in material.

As mentioned in the introduction, we have chosen four body centered cubic (bcc) transition metals(Fe, W, V, Ta) that are promising candidates for future fusion reactors. When subjected to extreme radiation conditions, a chain of collisions give rise to a displacement cascade. This displacement cascade causes atom displacements and leads to formation of point defects. Since migration and reactions of these point defects influence the physical and mechanical properties of the material, it is worthwhile to carry out a study of these point defects. In this thesis, we restrict our study to interstitial defects in the four bcc metals of interest.

Objective: In this chapter, we present a brief summary of radiation damage and a thorough survey of the relevant published articles about radiation damage in iron, tungsten, vanadium and tantalum. The various sections and their contents are organized as follows:

• Section 1.1: The radiation damage effects in materials are briefly recapitulated.

• Section 1.2: An extensive survey of articles reporting direct observation of microstructure of irradiated samples of iron, tungsten, vanadium and tantalum, and 1 State of the Art Review related computation-based studies are summarized. These articles are systematically discussed in the following sections:

-Section 1.2.1: The different morphologies encountered in bcc transition metals, according to experiments and Density Functional Theory calculations, are discussed.

-Section 1.2.2: Overview of the microstructural studies, comprising experimental observations and simulation-based results, in irradiated iron are presented.

-Section 1.2.3: Overview of the microstructural studies, comprising experimental observations and simulation-based results, in irradiated tungsten are presented.

-Section 1.2.4: Overview of the microstructural studies, comprising experimental observations and simulation-based results, in irradiated vanadium are presented.

-Section 1.2.5: Overview of the microstructural studies, comprising experimental observations and simulation-based results, in irradiated tantalum are presented.

• Section 1.3: Finally, the main points of this chapter are recapitulated systematically to provide a clear idea of the experimental observations, the simulation-based results and the correspondence between them.

Radiation damage

Depending on the type of radiation and their interaction with metals, radiation can be classified (16 -18 ) into heavy charged particles, fast electrons, neutrons and photons. The interactions of these different radiations in material cause some radiation effects, which can typically be classified into four categorie (16 , 19 ):

• Impurity Production: This refers to products of nuclear reactions that may be initiated by neutrons or even helium and hydrogen production on neutralization 1 State of the Art Review of alpha and proton, respectively.

• Atom Displacement: This refers to the shift or relocation of atoms from their original position in the material. Atom displacement is observed in materials mainly in the form of interstitials, vacancies, Frenkel pairs, dislocations etc. Examples of point defects are shown in Fig. 1.1

• Ionization: Charged particulate radiations are directly ionizing because they result in immediate stripping of electrons in the materials giving ion pairs, though uncharged radiation are also capable of causing ionization indirectly.

• Energy Release: Since radiation continuously loses energy via interactions with matter, this leads to its energy loss and consequent heating along its path in the material.

In this thesis, we concentrate on defects caused by atom displacement specifically in the form of interstitials. Irradiation-induced displacement of even one atom can potentially lead to a chain of collisions in an otherwise (almost) perfect lattice of the target material. This chain of collisions gives rise to a displacement cascade. Interstitials and other point defects are formed mainly in the collisional phase of a displacement cascade. The collisional phase is where creation of a primary knock-on atom(PKA) leads to energy transfer from incident particle and all collisions of PKA take place. These collisions lead to displacement of lattice atoms and a surge in the vacancy and interstitial concentrations. A typical displacement cascade simulation is shown in Fig. 1.2(a,b) and the displacement cascade timescale is shown in Fig. 1.2(c). As seen in Fig. 1.2(c), the evolution of this displacement cascade can broadly be divided into three overlapping phases (20 -22 ): the Collisional Phase (includes the creation of PKA at about 10 -18 seconds and PKA-caused collisions at about 10 -13 seconds with temperature rise in cascade core of up to thousands of degrees), the Cascade Cooling Phase (lasts till about 10 -11 seconds with annihilations, energy dissipation and some clustering lowering of the cascade temperature to nearly equal that of the matrix) and the Diffusional Interaction phase (includes defect interactions, defect migration and thermally activated defect reactions depending on the temperature conditions).

The migration and reactions of the created point defects result in physical and mechanical changes in the material. A pictorial depiction of the classification of defect reactions is shown in Fig1.3. Since a point defect is an isolated defect occuring when [START_REF] Zinkle | [END_REF].

1 State of the Art Review a single atom is missing from its original position in the lattice (vacancy) or an extra atom is found to be at an interstitial position instead of a lattice position (interstitial), it can be thought of as the fundamental unit of defects. Based on geometry, a collection of point defects produced in a crystalline material lattice can be further divided into line defects, planar defects(e.g., Dislocation loops) and volume defects(e.g., Voids, Bubbles, Stacking fault tetrahedra).

Since a displacement cascade disrupts the perfect lattice of the target material and can influence physical and mechanical changes in materials, it is important to estimate the defects produced. Models based on different assumptions have been proposed to estimate the number of Frenkel pairs generated by a primary knock-on atom (PKA). The number of Frenkel pairs is directly proportional to the damage energy, and inversely proportional to the threshold displacement energy, which is the energy available to generate minimal atomic displacements.

Practical experimental situations might not satisfy all the criteria for an ideal system because of presence of impurities or special experimental conditions. Despite this drawback, the observed defect morphologies in materials provide information about the most stable defect configuration. Relative stability of the possible irradiation-induced defects can also be theoretically estimated based on their free energies of formation which mainly consists of the formation energy contribution. Since calculation of formation energies are easier to perform, a comparison of the formation energies of defects provides a generally acceptable scale to gauge the relative stability of defects in the target material. Nevertheless, formation energies in pure metals give a reasonably good starting point to understand effects of radiation damage in high purity metals or even alloys and to compare our theoretical understanding. This makes experimental and simulation-based studies of irradiation effects essential in understanding the relative stability of defects and so, the following section is dedicated to the literature survey comprising both experimental and simulation-based results.

Throughout this thesis, we assume Miller index notation. For the benefit of the reader, we recapitulate the Miller indices for a 3D cubic lattice: Notation [abc] denotes a direction in the direct lattice vector while the notation (abc) denotes a plane in the reciprocal lattice vectors. Due to symmetry of a cubic lattice, direction [abc] is perpendicular to the plane (abc). Notation abc refers to all directions that are symmetrically equivalent to [abc]. Notation {abc} refers to all planes that are symmetrically equivalent to (abc). Ref. [START_REF] Was | Fundamentals of Radiation, Materials Science, Metals and Alloys[END_REF].

1 State of the Art Review

Literature Survey

Irradiation of bcc transition metals, in general, leads to a large number of point defects in the form of self-interstitial atoms and vacancies. Interaction of these point defects with each other or with sinks like dislocation and grain boundaries or their migration can influence microstructural transformations. This microstructural evolution affects various properties of the metals which ultimately dictate their lifetime in a nuclear reactor or under any specified radiation condition. Throughout this thesis, we will concentrate on interstitial point defects or self-interstitial atoms (SIAs). Following an extensive survey of microstructural studies comprising of experimental and simulation-based results in irradiated samples of iron, tungsten, vanadium and tantalum, a brief overview of the different morphologies of radiation-induced defects is presented in Sec. 1.2.1. This is followed by a summary of the main results for each of the four transition metals of interest in this study, the results are listed chronologically for Fe in Sec. 1.2.2,for W in Sec. 1.2.3,for V in Sec. 1.2.4 and for Ta in Sec. 1.2.5. The experimental literature survey concerns only the limit of large loops, visible in experiments. In the limit of very small loop size, the old resistivity recovery or internal friction experiments are mostly in agreement with the theoretical findings.

Description of the different morphologies of defect clusters encountered in bcc metals

Transmission electron microscopy (TEM) observations show that in all irradiated bcc metals except Fe, dislocation loops with 1/2 111 Burgers vector are dominant which suggests that they are the most stable configurations for bundles of dumbbells. In Fe, observations of nanometric-sized clusters of SIAs by TEM techniques reveal the presence of planar loops, which can adopt either the 1/2 111 (highly mobile) or 100 (immobile) configurations, depending on temperature (24 -26 ). The relative stability of the two types of loops in Fe is influenced by magnetism, and it has been shown that 1/2 111 loops are more stable at low temperatures while 100 loops are more stable (27 , 28 ) at high temperatures (over 700 K). Recently, much progress has been made in the experimental field, enabling observation of small 100 loops in W under heavy ion irradiation at low temperatures, which vanish at high temperatures (29 , 30 ). The reason for the formation of the 100 loops in W is still under debate in literature. In the low size 1.2 Literature Survey limit, we can access information about the morphologies of SIA using resistivity recovery or internal friction experiments. In the intermediate defect cluster size range, spanning the interval between individual self-interstitial atoms and nanometric-sized dislocation loops, it is difficult to generate experimental data because of the high resolution of observations required to characterize such small objects.

Density Functional Theory (DFT) calculations provide quantitative insight into the nature of clusters containing a small number of defects. These DFT predictions broadly agree with experiment [START_REF] Ehrhart | Physics of Radiation Effects in Crystals[END_REF], which makes it desirable to extend predictions to clusters larger than a single SIA by packing dumbbells together in bundles, to form small dislocation loops. Fig. 1.4 shows formation energies of different mono-SIA interstitials relative to 111 dumbbell for various bcc transition metals [START_REF] Nguyen-Manh | Proceedings of the Twelfth International Conference on Fusion Reactor Materials (ICFRM-12)[END_REF]. As seen in this figure, DFT calculations show that the most stable single SIA in Fe adopts a configuration that corresponds to a 110 dumbbell, whereas in other bcc transition metals, a single SIA forms a defect aligned along the 111 direction, known as a crowdion (32 -36 ) (a few examples are shown in Fig. 1.5). According to recent DFT calculations (37 ), SIA clusters can also form three-dimensional structures in Fe with symmetry corresponding to the C15 Laves phase.

Being a newly proposed interstitial configuration, the geometry of C15 cluster is not very well known. Here, we discuss the structure of these clusters. The starting point of the C15 family of SIA clusters is the di-interstitial. The C15 interstitial configuration was initially obtained by combining the triangle and hexagonal ring building-blocks, as illustrated in Fig. 1.6. This construction illustrates the fact that making a closed-shell structure decreases the net number of additional atoms in the lattice from 3 to 2. In terms of SIAs and vacancies, this defect corresponds to 12 SIAs, placed at the edges of a truncated tetrahedron, surrounded by 10 vacancies, which indeed makes a total of two additional atoms in the bcc lattice. One can recognize that the highly symmetric structure formed around the central atom by its 12+4 nearest neighbours corresponds to the Z16 Frank-Kasper polyhedron [START_REF] Frank | [END_REF](see Fig. 1.7a). Another Z16 polyhedron with 6 neighbours in common can then be constructed by adding 6 SIAs and 4 vacancies. This makes a quadri-interstitial cluster, I 4 (see Fig. 1.7c). The two polyhedra have two different orientations and they are centred on two nearest neighbours of the bcc lattice. Other polyhedra can be added, forming I 6 and I 8 clusters as illustrated in Figs. 1.7g-h. By repeating this process, 3D I n clusters with a cubic periodic structure can be built. This crystallographic structure corresponds to the C15 Laves phase, or M gCu 2 structure (see Fig. 1.7i) (40 ). Detailed construction of bigger clusters will be presented in chapter In other bcc transition metals the bulk C15 structure indeed also has a low energy in particular in group VB metals where the energy relative to the bcc structure is comparable to that in Fe, but DFT calculations performed for quadri-interstitials in V, Nb, Ta, Cr, Mo and W show that the I C15 4 cluster is always significantly higher in energy than the configuration made by four parallel 111 crowdions. The energy difference is typically 2 to 3 eV, except in Ta where it is only 0.8 eV. Energy landscape of C15 clusters has been investigated in Fe, W, Ta and V using DFT calculations as part of this thesis and will be discussed further in chapter 4. The atypical behaviour of Fe with respect to the other bcc transition metals confirms the specificity of Fe regarding SIA properties (42 ). In Fe, these C15 aggregates are stable, immobile, and exhibit large antiferromagnetic moments. These C15 clusters have been found to form directly inside atomic displacement cascades, and are able to grow by capturing self-interstitial atoms from the surrounding material. without the vacancies and the cubic lattice. g, h: skeleton representation for the hexa-and octa interstitials. i: unit cell of the C15 Laves structure. For each cluster size, the configuration which is represented corresponds to the lowest energy one found within DFT.

1 State of the Art Review

Overview of radiation-induced defects in iron

Experimental observations

In irradiated iron, the presence of point defect clusters had been indicated by Mogford (43 ) who reported platelets of Carbon atoms aligned along {100} and Eyre (44 ) who observed point defect clusters (glissile dislocations sweeping up defects forming jogs).

The crystallography of dislocation loops in Fe was first observed by B. C. Masters (24 , 45 ). Masters irradiated 0.005 in foils of iron with Fe+ ions of energy 150 keV to a dose of 1.7×10 16 ions/cm 2 at specimen temperature of 823 K (550 • C). He observed 'interstitial' loops with 100 Burgers vector on {100} planes with absence of vacancy loops. The results of this paper were unexpected due to two reasons: Firstly, absence of vacancy loops. Masters reasoned with uncertain choice of nucleation (homo or hetero) for this absence and suggested further work using materials of high purity. The more surprising result at that point in time was the presence of only 100 SIA loops. Since interstitial loops with 1/2 111 Burges vector had already been observed in bcc Molybdenum (46 , 47 ), similar results were expected for iron, as in (48 ). Considering three modes of lowest energy as {111} 1/2 [111], {110} 1/2 [1 11 ] and {001} [001], he proposed a mechanism to explain formation of 100 by a reaction as follows:

a 2 [111] + a 2 [1 11 ] = a [100] . (1.1)
Further, he explained that 1/2 111 loops probably escape detection because they glide out of the thin foil during irradiation.

Little and Eyre (49 ) studied the geometry of dislocation loops in 0.25 mm iron foils using 1MeV electron irradiation at 823 K (550 • C). They observed rectilinear 100 loops of interstitial nature, endorsing the mechanism for formation of 100 explained in (50 ). They concluded that the formation of 100 loops are not just a consequence of excess interstitials as suggested by Masters because 100 loops were observed with electron irradiation as well.

Heavy ion damage in alpha iron was studied by English et al. (51 , 52 ). They irradiated iron foils at room temperature with ions ranging from Fe+ to W+ with atomic weights ranging from 56 to 184 and low dose of 5 × 10 12 ions/cm 2 . They did not observe any visible damage for self-ions. Defect yield increased with ion mass and both 1/2 111 and 100 loops of vacancy type were observed. Formation of 100 loops was surprising 1.2 Literature Survey since it was known to be energetically unfavourable.

A TEM study of low fluence neutron-irradiation with dpa=0.5-1 of 0.5 mm iron disk specimen (3 mm diameter) at temperatures of 455-1013 K was carried out by Horton et al. (53 ). Different microstructures were observed in the four temperature ranges. From 455 to 523 K, small defect clusters near pre-irradiation dislocation segments were observed but Burgers vector determination was inconclusive. In the range 548-573 K, clusters of small loops were seen and identified as 100 loops of interstitial type. In the 623-773 K range, 100 {100} loops of interstitial type were observed with a decrease in defect density and an increase in defect diameter. There was also a single 1/2 111 loop observed in this temperature window. There were no visible changes compared to unirradiated specimens in the temperature range 923-1013 K. According to them, presence of 100 loops was justified because of their less glissile nature and presence of Cottrell atmosphere which pins these loops while the behaviour at higher temperatures (>723 K) was attributed to small precipitates (presumably, carbides). This is one of the first papers to question the validity of the Eyre-Bullough mechanism for the formation of the 100 loops on grounds of observation of solely 100 loops at higher temperatures and total absence of 1/2 111 loops which are energetically favoured.

Robertson (54 ) investigated low dose neutron irradiation damage where thin foils of iron were irradiated to doses between 10 22 to 10 24 neutrons/m 2 by 1MeV neutrons at ambient temperature. Here, loops with both 100 and 1/2 111 Burgers vectors were identified with habit planes between the {110} nucleation and edge plane but nature of loops was not determined. Post-irradiation annealing at 643 K for 60 min did not show much change. Following the same irradiation and post-irradiation annealing with dose of 4 × 10 23 gave a majority of 1/2 111 and 33 percent of 100 while a dose of 1.2×10 24 gave higher percentage of 100 loops with all loops being of interstitial nature. In accordance with the Eyre and Bullough mechanism (50 ), an elastic continuum model by Little et al. (55 ) was used to calculate the relative probability of unfaulting of an a/2 110 loop into an a 100 loop or an a/2 111 loop at reactor ambient conditions. This value was found to be 6.5 × 10 -17 which could not explain results by Robertson, even if anisotropy of iron were to be included for calculation in this model. Robertson et al. (56 ) investigated damage structure in alpha-iron using 50 keV self ions up to a maximum dose of 8 × 10 13 ions/cm 2 at 40 K. They observed no damages below 8 × 10 13 ions/cm 2 and a non-linear increase in defect density at and above this dose. Of the loops observed after irradiation with self-ions of 50 keV, 77 percent were 1 State of the Art Review 1/2 111 while the remainder were 100 loops. The loop nature was not determined. These findings revealed that dislocation loops in iron are formed by spatial overlap of cascades and not from direct collapse of individual displacement cascades as for other metals at low temperatures. Arakawa et al. (57 ) reported observations of a new type of changing process in the Burgers vector of dislocations by simple heating or high-energy electron irradiation by in situ transmission electron microsopy. According to this article, small 1/2 111 interstitial loops of diameter < 50nm are capable of transformation into another 1/2 111 loop or an energetically unfavorable 100 loops while these 100 loops can further transform into a 1/2 111 loop. On simple heating, it was seen that small 1/2 111 loops (diameters <20-30 nm) exhibited 1D thermal motion parallel to its Burgers vector at temperatures above 450 K while 100 loops exhibited similar glide motion at higher temperatures (<770 K). Meanwhile, on irradiation, 1/2 111 loops (diameter <50 nm) showed 1D motion at low temperatures whereas 100 loops rarely showed any 1D motion.

Yao et al. studied heavy-ion irradiation of thin foils of Fe and FeCr alloys at lower doses (58 ). They irradiated Ultra-High Purity (UHP) Fe and FeCr alloys with 100 keV or 150 keV Fe+ or 100 keV Xe+ at room temperature(RT) and 573 K (300 • C). A threshold value of dose (2 × 10 16 ions/m 2 ) was noted, below which no defect was observed. Moreover, a low proportion of 1/2 111 loops for Fe indicated their loss from the foil due to high mobility and possibility of glide cylinders intersecting the surface of foil. This suggested that loops were more mobile in Fe than Fe-Cr alloys for same irradiation condition.

Hernandez et al. (59 ) continued the above study at higher doses i.e. up to 2 × 10 18 ions/m 2 i.e. 13 dpa. For UHP Fe at 573 K (300 • C), strings of loops were observed at a dose of 2 × 10 18 ions/m 2 . These dislocation loops grew to observable sizes with loop-types depending on thickness of the Fe foil. A dense dislocation network for high thicknesses (>100 nm), 1/2 111 interstitial loops were observed for intermediate thickness (between 50 and 100 nm) and 100 finger-shaped loops for low thickness (< 50nm). They endorse the model proposed by Wen et al.(60 ) where 100 loops are assumed to be formed by growth and coalescence of small loops. UHP Fe at room temperature showed similar results. At low doses (300 • C), Fe-Cr results in irradiation damage similar to Fe results but with lower loop number density. At higher doses, smaller maximum loop sizes and larger loop number densities are observed due to a lower mobility of small loops in Fe-Cr compared to Fe. So, the effect of Chromium is the reduction of mobility 1.2 Literature Survey of loops which is consistent with the results at lower doses (58 ) Yao et al. (28 ) studied temperature dependence of microstructure in Fe when irradiated with 150 keV Fe+ ions. At 673 K (400 • C), almost equal proportions of both immobile 100 and mobile 1/2 111 were observed. While 1/2 111 loops coalesced into larger loops of round or irregular shape, 100 loops were comparatively small in size. At 723 K (450 • C), rectilinear 100 loops dominated and 1/2 111 loops interacted elastically to form long chains. These long chains of 1/2 111 loops ceased to form above 723 K (450 • C) because small, mobile 1/2 111 loops were consumed by 100 loops such that only large 100 loops existed at 773 K (500 • C). They could not explain the abrupt disappearance of 1/2 111 loops above 738 K (465 • C). However, they did not report any direct conversion of 1/2 111 loops into 100 loops, contrary to Arakawa's observation (57 ). To give a plausible explanation, they did lend support to Bacon's MD simulations (61 ) which relied on small, self-interstitial clusters (bundles of 100 or 111 crowdions) formed after displacement phase of the cascade. Being mobile and able to execute thermally-activated one-dimensional glide along the crowdion direction, these small clusters were considered responsible for observable 111 and 100 loops by coalescence. They postulated that the relative proportions of the two loops depended on a number of factors such as loop stability at smaller sizes and mobility of 1/2 111 loops at higher temperatures to explain their experimental observations. Irradiation of iron with 150 keV Fe ions at 573 K (300 • C) gave a vast majority (90%) of 100 loops, compared to room temperature irradiations with 150 keV Fe ions(63%), 30 keV Ga ions(55%) or 100 keV Fe ions(86%). This paper was found to be consistent with theoretical papers by Dudarevet al. (27 , 62 ), discussed later in this section. Xu et al. (63 ) studied the effect of 150 keV Fe+ ion irradiation on Fe and Fe-Cr alloys at 573 K (300 • C). 92% of the loops were of type 1/2 111 in Fe while the percentages were 30%, 46% and 37% in Fe-5%Cr, Fe-8%Cr and Fe-11%Cr suggesting lack of any pattern with respect to the percentage of Cr in the alloy. All the loops were reportedly of interstitial nature. These experimental results were at variance with some earlier experiments, notably those by Masters(45 ) and Yao et al.(58 ) which were discussed in the previous paragraphs. Experimental observations by Gelles (64 , 65 ) and Porollo et al. (66 ) are also contradictory because they both reported increase in the fraction of 1/2 111 loops with increasing percentage of Cr in Fe-Cr alloys after neutron irradiation up to 15 dpa at 673-723 K (400 -450 • C) and up to 7 dpa at 673 K (400 • C), respectively.

Literature Survey

and FeCr alloys when irradiated with 100-150 keV Fe+ and Xe+ of dose up to 13 dpa at room temperature, 573 K (300 • C) and 773 K (500 • C). They observed both 100 and 1/2 111 loops at room temperature(RT) and at 573 K whereas only 100 loops existed at 773 K (500 • C). Although similar radiation damage was observed in FeCr alloys, defect yield did not show as much dependence on foil orientation as seen for Fe.

Nature of loops at low doses wasn't reported with certainty but they were reportedly of interstitial nature at higher doses.

Prokhodtseva et al. [START_REF] Prokhodtseva | [END_REF] studied bulk and thin foil of Fe after irradiation with Fe+ at room temperature with or without a second He beam. In bulk, single beam resulted in more than 75% of loops being 1/2 111 loops. In thin foil, about 96% loops were 1/2 100 loops and rest were 100 loops with single beam. The fraction reversed with double beam and about 99% were reported to be 1/2 111 loops along with 1% of 100 loops. These results prove that He stabilizes the 1/2 111 loops rendering them immobile, as predicted by MD simulations (69 ). Due to absence of free surfaces, there is no strong driving force to cause escape of 1/2 111 loops or their (then proposed) transformation into 100 loops in bulk. Prokhodtseva et al.(70 ) extended the above experiment to FeCr alloys as well, reiterating the role of He in stabilizing the mobile 1/2 111 loops and preventing their (proposed) transformation into 100 leading to higher proportions of 1/2 111 loops. Since the early loop population is generally dominated by 1/2 111 loops, they postulated that the 100 loops must be formed from addition/absorption reactions between them. It was concluded that Cr in Fe-Cr alloys reduce the mobility of the 1/2 111 loops as well, eventhough the combined effect of He beam and Cr content in alloys did not show any clear trend.

Simulation-based results

In 1965, Eyre and Bullough (50 ) proposed some explanations for the observation of 100 and 1/2 111 loops in experiments. They based their arguments on computer calculations by Erginsoy et al.(71 ) and Johnson et al.(72 ) which indicate that 110 single and double interstitials are stable. They postulated growth of 110 di-interstitials to ultimately form a platelet of interstitials in the {110} plane, further leading to a stacking fault in the {110} plane. This stacking fault can be eliminated by the two 1 State of the Art Review following dislocation reactions:

a 2 [110] + a 2 [ 111] = a [010] , (1.2) a 2 [110] + a 2 [00 1] = a 2 [11 1] , (1.3)
leading to formation of 100 and 1/2 111 loops, seen in experiments. They explained that the need for a higher energy 110 shear requires thermal energy assistance and hence leads to formation of 100 loops only at higher temperatures. According to experimental results from Masters (45 ), 100 loops were observed on irradiation of iron with Fe+ ions at 823 K (550 Meanwhile, developments in MD simulations progressed to explain experimental results obtained by radiation damage in iron. Some of the earliest MD simulations on radiation damage in bcc iron were those by Calder and Bacon (73 ), Phythian et al. (74 ), Stoller et al. (75 ) and Soneda et al. (76 ). Though various experimentalists continued to lend support to the Eyre and Bullough mechanism (50 ) for formation of 100 loops, this mechanism was considered unlikely on the basis of MD simulations from Calder and Bacon (73 ). However, these early MD simulations ignored the effect of glissile clusters whose 1D transport is essential to understand the damage accumulation. This led to the development of the production bias model used by Osetsky et al. (77 ,78 ).

With improvement in MD simulations, many mechanisms were proposed to explain the formation of 100 loops at high temperature in iron. Marian et al. (79 ) proposed formation of [001] loops from two 1/2 111 loops, in accordance with Masters's experimental results (45 ) using a modified Eyre-Bullough mechanism. They found using MD simulation the following two-step mechanism:

1 2 [111] + 1 2 [00 1] → 1 2 [110] , 1 2 [110] + 1 2 [1 10] → [100] ,
For such reactions to occur, they postulated that the intersecting loops were required to be of the same size, (preferably) same shape, should stabilize and grow 100 -type segments. However, they did not provide complete proof of this mechanism using MD 1.2 Literature Survey simulations. Further, they suggested absorption of smaller loops by larger loops using MD on nanometric-sized 1/2 111 loops, in accordance with Osetsky et al. (78 ) and Terentyev at al. (80 ). In 2008, Dudarev et al. (27 , 62 ) showed the importance of using the full anisotropic elasticity approximation at higher temperatures where the alphagamma phase transition occurs in iron. They calculated free energies of loops using Bacon's formalism [START_REF] Bacon | [END_REF] to explain increasing proportion of the 100 loops observed experimentally above 573 K (300 • C). However, mechanism of 100 loop formation remained unexplained.

As a consequence of development of new MD simulations, new interatomic potentials were formulated simultaneously. Commonly used many-body interatomic models of metals span from second moment of tight binding approximation (82 -86 ) to EAM (87 , 88 ) or Modified EAM (89 ) and higher order tight binding models (90 -95 ). The Force-Matching method proposed by Ercolessi and Adams in the late 90s [START_REF] Ercolessi | [END_REF] showed for the first time that the database used in the fit of the empirical potentials could play the same key role as the formalism used to mimic the electronic effects. The force-matching method was used to parametrize EAM potentials suitable for defects in bcc iron. Using the same fitting approach, Mendelev et al. (97 ,98 ) and Ackland et al.(99 ) proposed a parametrization of iron. Using a different database (point defect oriented), Marinica (37 , 100 ) obtained an improved parameterization for point defects. In 2005, Dudarev and Derlet designed an alternative approach (101 ) with a 'magnetic' potential which was also based on the EAM formalism. These empirical potentials permit simulation of cascades facilitating studies involving defects which is otherwise not possible when using more accurate means. However, the results may vary from one potential to other and also from DFT calculations. Recently,Ref.(102 ) showed that the threshold displacement energies, which serve as an important input for radiation damage studies, are significantly different when calculated using DFT-MD and from a widely-accepted interatomic potentials. In fact, the average value from the interatomic potential was found to be about 20% lower than the average standard value in literature. This raises concerns over the quality of empirical potentials available for MD simulations and the reliability of MD results. This means that the quantitative and predictive studies based on empirical potentials cannot be conclusive, though the results provide a good qualitative estimation. In chapter 3, we develop a new empirical potential to improve the existing empirical potentials and in chapter 4, we present a new model to overcome the shortcoming of empirical potentials.

Marian's mechanism (79 ) was at odds with Kirchhoff's law (103 ) which suggests formation of a [001] junction leading to a sessile dislocation-loop complex. In order 1 State of the Art Review to find the mechanism experimentally, Arakawa et al. (104 ) irradiated thick sheets of iron with electrons up to 3 × 10 2 5e -/m 2 at 110-190 K, followed by heating up to 290-800 K. They showed that collision of 1/2[111] and larger 1/2 [11 1] gives an even larger 1/2 [11 1] after absorption of the former, shown in Fig. 1.9. In general, collision of two loops with different Burgers vectors led to absorption of the smaller loop by the larger loop. This experiment provided evidence for absorption of smaller loops by larger loops but did not observe the formation of 100 loops from 1/2 111 loops, eventhough Xu and Terentyev et al. (105 ) showed MD simulations using two different potentials, A04 (99 ) and M10 (100 ). This still didn't provide conclusive evidence to prove this mechanism because experimental evidence of complete absence of 1/2 111 loops at higher temperatures goes against this mechanism (45 ). The main conclusion of the experimental and theoretical papers is that the proposed mechanism of formation of 100 loops can be realised only under some strict conditions on the size and the orientation of the loops. Another mechanism was proposed by Chen et al. (106 ) which demonstrated using MD simulations that rearrangement and reorientation could lead to change of small 1/2 111 loops into 100 loops.

DFT predicts that in Fe the orientation of these dumbbells changes from 110 to 111 depending on the number of SIA involved. The transition occurs around five SIAs (37 , 107 ). As mentioned earlier, Dudarev et al.(27 ) took into account the full anisotropic elasticity approximation to calculate the formation energies near the αγ transition. This treatment explained the reduction of free energy of 100 formation and their consequent stability at temperatures higher than about 573 K (300 • C). This paper will be further discussed in detail in the first part of the chapters 4 and 5. Although the stability of 100 was thus established, the mechanism for the formation of 100 loops has not been confirmed yet.

According to recent DFT calculations (37 ), SIA clusters can also form three-dimensional structures in Fe with symmetry corresponding to the C15 Laves phase. Prior to C15 clusters, the lowest energy structures according to DFT calculations published for SIA clusters in bcc-Fe were (36 , 108 ): the triangular configuration for the di-interstitial, the hexagonal ring configuration for the tri-interstitial, the parallel 110 dumbbell configuration for the quadri-interstitial, and compact arrangements of parallel 111 dumbbells for larger clusters. These low energy structures are shown in Fig. 1.5. The parallel dumbbell configurations have slightly larger energies for di-interstitial I 2 and tri-interstitial I 3 . The quadri-interstitial configuration made of a 111 dumbbell centred on an hexagonal ring had attracted attention because of its rather low energy and very low mobility (36 , 1.2 1 State of the Art Review 109 , 110 ). As shown in Ref. (37 ), these stable structures are actually particular cases of stable C15 clusters.

As reported in Ref. (37 ), the remarkable property of these C15 SIA clusters is that they appear to be particularly stable in bcc Fe according to DFT calculations. For a given n, the possible configurations for C15 clusters multiply which makes the construction of the most stable configuration challenging. The present work also includes the rules for construction of stable n-SIA C15 clusters which will be discussed in Sec. 3.4. DFT calculations were also done to understand dumbbell migration barrier in Fe-Cr as compared to pure Fe in Ref. (111 ). It was found that the pure Fe interstitials had a higher dumbbell migration barrier as compared to mixed interstitials in accordance with experiments.

As mentioned earlier, the mechanism of formation of 100 loops in Fe as a result of high temperature irradiation was not explained so far. Quite recently, Zhang et al. (112 ) proposed an explanation involving the nucleation of C15 clusters and their growth by trapping of single self-interstitials of 1/2 110 dumbbell structure. They used MD simulations to show that these C15 clusters can transform into both 100 and 1/2 111 loops after absorbing SIAs at small sizes. This mechanism will be further discussed in chapter 4. Although there is no experimental evidence to support the formation of these 3D C15 clusters or the theory proposed by Zhang et al.(112 ), it is worthwhile to discuss them in some detail considering their importance in understanding the energy landscape of defects in Fe and in reconciling experimental observations with the calculated energy landscape of defects in Fe. In preceding studies, the energy landscape of defects in Fe is predicted using empirical potentials only but it will be predicted on the basis of DFT calculation in chapter 4 of the current study.

Summary for iron:The number of experiments studying the irradiation-induced microstructural changes in iron is quite high compared to the other transition metals. All these experiments report a majority of mobile 1/2 111 interstitial loops at lower temperatures. Further, an increasing proportion of the comparatively immobile 100 is observed above 573 K which ultimately become the majority at around 773 K. The orientation of the sample is crucial in this observation, otherwise the 1/2 111 can migrate to surface, biasing the results of experiments. Moreover, the impurities can also bias the results impacting the mobilities of the different loop types. For example, He stabilizes the 1/2 111 loops rendering them immobile, as predicted by MD simulations. Cr in Fe-Cr alloys reduce the mobility of the 1/2 111 loops as well, eventhough the combined 1.2 Literature Survey effect of He beam and Cr content in alloys did not show any clear trend. Though the theoretical explanation was lacking till recently, most of the experiments have now been explained rather convincingly in the existing literature. According to DFT Calculations, the most stable mono-SIA in Fe is the 110 dumbbell. The 110 SIAs are the most stable in Fe up to around five SIAs and 1/2 111 loops are the most stable beyond the cluster size of around five interstitials. Taking into account elastic approximation to calculate the formation energies near the αγ, a reduction of free energy of 100 loops was demonstrated, thus making the 100 loops more stable compared to 1/2 111 loops at temperatures higher than about 800 K. Although the stability of 100 was thus established, the mechanism for the formation of 100 loops has not been confirmed yet. The theoretically proposed mechanism of the formation of 100 by collisions of two 1/2 111 loops can explain the formation of 100 loops in pure Fe, or Fe under He atmosphere. However, despite the fact that Cr reduce the mobility of 1/2 111 loops the formation of 100 loops is very active in F eCr alloys. Moreover, the absence of 100 interstitial dislocation loops at low temperatures remains unexplained, given that the 1/2 111 loops are very mobile at low temperature. Recent studies highlight the importance of C15 clusters in understanding and explaining the anomalous experimental observations of 100 loops in Fe, thus reconciling experiments with DFT calculations of energy landscape of point defects in bcc Fe. The sections 4.5.1 and 5.4 will provide a detailed discussion of the existing models and the present findings of this thesis.

Overview of radiation-induced defects in tungsten

Experimental observations A number of articles were published on neutron irradiation effects on tungsten and its alloys. However, these articles concentrated on effects of various irradiation conditions and transmutation elements on radiation hardening and electrical resistivity (113 -119 ). A summary of some experimental observations has been shown in Fig. 1.10, as reported by Hasegawa et al. (120 ), where the presence of loops and voids as a function of irradition temperature and displacement damage are marked in the graph. Eventhough these experiments did not involve a systematic study regarding the Burgers vector or nature of loops, they did provide a useful basis for future investigation in terms of irradiation conditions. In particular, the irradiation conditions for observation of loops and voids or only voids appeared to be quite well defined based on these experiments. Additionally, 1 State of the Art Review some studies also concentrated on interaction of tungsten with solute atoms (121 -125 ).

The experimental database studying the microstructural evolution of radiation damage in tungsten is not as vast as that for iron. In 1972, Häussermann et al. (126 -128 ) studied microstructural changes in thin foils of tungsten irradiated with 30-60 keV Au ions at room temperature. They reported a majority of small 1/2 111 loops of vacancy type along with a few 1/2 110 loops. This led them to presume the Eyre-Bullough mechanism (50 ), originally proposed for bcc iron (described in section 1.2.2). In keeping with this mechanism, they assumed aggregation of loops of 1/2 110 Burgers vector in the {110} plane till a certain critical size. Following which, they transform into 1/2 111 or 100 loops in order to remove stacking fault.

Jager and Wilkens (129 ) performed a similar experiment by irradiating tungsten foils with 60 keV Au ions at room temperature with dose of the order of 10 11 to 10 12 ions/cm 2 . Like Häussermann, they also suggested nucleation of 1/2 110 vacancy loops on {110} plane till they reach a critical radius. The critical radius for transformation into perfect Burgers vector loops by shear process was estimated in this article to be in the range 0.9 nm-2 nm. Moreover, 90% of the observed loops were either 1/2 111 or 110 , both in the {110} plane while the rest were assumed to be 100 loops as expected by the Eyre-Bullough mechanism.

Recently, Yi et al. (130 ) irradiated tungsten foils with 150 keV tungsten ions at 773 K (500 • C) with dose in the range 10 16 -10 18 W + m -2 at a rate of 6.25 × 10 14 W + /m 2 /s. They observed a majority(75%) of 1/2 111 loops, rest being 100 loops. Further, they reported various kinds of elastic interaction among loops, e.g. dragging to form strings, change of Burgers vectors of smaller loop while interacting with larger loops (104 ), absorption and coalescence. The nature of the loops was not solely vacancy, in contrast with earlier observations by Häussermann et al. (128 ) and Jager et al. (129 ). While 1/2 111 loops were found to be of interstitial and vacancy type in equal numbers, 100 loops were more of interstitial type. They suggested a possible loss of 1/2 111 loops by glide. Yi et al. (30 ) reported experimental characterisation of tungsten foils irradiated with 2 MeV self-ions at temperatures 573 K(with 3.3 × 10 17 , 10 18 and 2.5 × 10 19 W + /m 2 ), 773 K(with 3.3 × 10 17 , 10 18 and 3 × 10 18 W + /m 2 ) and 1023 K(with 10 18 W + /m 2 ) at 3 × 10 -4 dpa/s. This characterisation revealed a majority of 1/2 111 loops (> 60%) alongwith 100 loops for all irradiation conditions studied. These loops were mostly of Ferroni et al. (131 ) performed characterisation of tungsten foils irradiated with 2 MeV self-ions at temperatures ranging from 1073 to 1673 K (800 -1400 • C) with a dose of 10 14 ions/cm 2 at a rate of 2.73 × 10 10 ions/cm 2 /s. This article reported that loops were exclusively of 1/2 111 type and predominantly of interstitial type.

Simulation-based results

Since ab-initio calculations (132 ) in W showed that the formation energies of 1/2 111 loops are always less than those of 1/2 110 or 100 loops, the Eyre-Bullough mechanism didn't seem probable. Gilbert et al. (133 ) compared formation energies of various vacancy and interstitial defects in iron and tungsten using zero Kelvin molecular dynamics calculations. In tungsten, they concluded that spherical voids were the most stable energetically for a given number of vacancies. Assuming Arrhenius diffusion, they attributed difference in experimental observations to large characteristic equilibration time at room temperature. For iron, this characteristic time is around 8 s whereas it is of the order of 4 × 10 20 s for W. They postulated that extremely small value of vacancy diffusion coefficient and non-equilibrium conditions led to observation of interstitial or vacancy loops which are otherwise energetically unfavourable.

Different potentials for tungsten were published which varied in their development and also their accuracy or applicability in a particular study. For example, one of the first potentials developed for tungsten is the Finnis-Sinclair potential (83 ) which is a tight-binding second moment short-ranged potential and is fitted to fundamental bulk properties. This potential did not work well for studying self-interstitial atom defects, so Ackland and Thetford (85 ) enlarged the database accordingly to develop a new potential that is well adapted for studies ranging from bulk related properties to energy landscape of point defects and dislocations in tungsten. Another variant of this potential was proposed by Juslin and Wirth (134 ) which is developed mainly by a modification of the short range part of the potential by Ackland-Thetford and has been used to study helium-vacancy clustering in tungsten. A new parametrization of the Finnis-Sinclair type potential, fitted to include ab initio point defect formation energies, was given by Derlet et al. (33 ). This potential has been mainly found application in study of thermally acti-1.2 Literature Survey Figure 1.11: (a) shows formation energies of 100 loops in black and 111 loops in red as a function of cluster size, calculated using an embedded atom method(EAM) potential (135 ). A crossover is predicted here which is contradictory to calculations using some other potentials. (b) shows a plot of formation energy of 111 loops as a function of size using two different potentials, an EAM potential (135 ) in blue and potential by Ackland (85 ) in red.

vated migration of self-interstitial atoms but had certain drawbacks as well. Marinica et al. (135 ) used force-matching method to parametrize embedded atom method potentials. These potentials demonstrated qualitative agreement with DFT calculations on screw dislocations which wasn't achieved by other existing potentials but there remained some scope of improvement for properties like formation energies of self-interstitial atoms.

Although these potentials provide a feasible means to study radiation effects, reliability becomes an issue when predictions or calculations from different potentials contradict.

Consider the relative stability of point defect clusters some potentials (33 , 135 ) predict a crossover between 100 and 111 loops as shown in Fig. 1.11(a) while others didn't (85 ). Also, a significant difference in values is observed from different potentials, as shown in Fig. 1.11(b). These contradictions cause empirical potential calculations to be inconclusive.

Various DFT-based studies were carried out in order to furnish reliable calculations. Notably, Becquart et al. (121 ) calculated vacancy formation energy and migration 1 State of the Art Review energy which were found to be in good agreement with experimental values. They also concluded that the 111 loops were the most stable in tungsten and repulsive behaviour of di-vacancy which was later confirmed by Ventelon et al (136 ). Ab initio calculations have also been used as database to parametrize object kinetic Monte Carlo code to simulate isochronal annealing experiments of He desorption from W (137 ). Sand et al. (138 ) performed MD simulations of 150 keV collision cascades using PAR-CAS code (139 -141 ) to have a direct comparison with experimental observations of Yi et al. (130 ) using three different potentials. A modified EAM potential (33 ) having the repulsive part given by the universal Ziegler-Biersack-Littmark potential (142 ) and filled by Björkas et al. (143 ) was compared with Ackland-Thetford potential (85 ) and a Tersoff-type potential by Ahlgren et al. (144 ). In accordance with experimental results, they reported 100 and 1/2 111 loops of interstitial as well as vacancy type, the nature being sensitive to the potential used. As discussed recently in Ref.(145 ), DFT threshold displacement energy calculated with and without semi-core electrons vary significantly which affects the quality of the developed empirical potential and thus, the MD simulations. Mason et al. (146 ) reported experimental results and object Kinetic Monte Carlo simulations on irradiation of tungsten foils with self-ions of low fluence. The largest loops were found to be mainly composed of prismatic 1/2 111 loops of vacancy type. These experimental results were well reproduced by simulations on considering elastic interaction between loops. Setyawan et al. (147 ) carried out MD simulations using the Juslin-Wirth potential (134 ) to document the defects following displacement cascades at 300 K, 1025 K and 2050 K in tungsten. At 300 K and 1025 K, vacancy defects (>50) were predominantly 100 {100} loops or cavities while SIA clusters (>30) were 1/2 111 111 or 1/2 111 {110} loops. Rarely 100 {110} interstitial loops were reported at 1025 and 2050 K while 3D interstitial clusters dominated at 2050 K.

Summary for tungsten: While recent experimental observations reveal presence of a vast majority of 1/2 111 and a few 100 interstitial loops at low temperatures, these 100 interstitial loops disappear leaving only 1/2 111 interstitial loops at high temperatures. Empirical potentials seem to be inconclusive with some predicting 1/2 111 as most stable and others show a crossover between 1/2 111 and 100 loops suggesting competing stability of these loops as a function of number of interstitials. As such, MD 1.2 Literature Survey simulations using different potentials can be contradictory as well but DFT calculations have quite definitively established the stability of 1/2 111 interstitial loops over other possibilities.

Overview of radiation-induced defects in vanadium Experimental observations

In 1968, Rau et al (148 ) were among the earliest to report the effects of neutron irradiation on vanadium. Having irradiated vanadium disks with fast neutrons at temperatures from 573 to 1448 K (300 -1175 • C), they reported defect clusters that grew into resolvable dislocation loops (apparently) of interstitial nature. The microstructure was not reported due to impurities which were responsible for stability of loops at higher temperatures. These impurities, appearing as black spots decorating dislocation loops, were actually precipitates (most likely vanadium carbides or nitrides) formed by their interaction with vanadium. Additionally, they reported presence of damage halos around some particles (probably boron)which were supposed to occur due to recoiling products of 10 B(n, α) 7 Li. 149) irradiated ultra pure vanadium disks with fast neutrons and reported their observations at 693 K (420 • C), 903 K (630 • C) and 1023 K (750 • C). At 693 K, they observed three types of defects: aggregating small dislocation loops of Burgers vector 1/2 111 , inhomogeneous distribution of voids and radiation-induced planar precipitates involving impurities such as oxygen, carbon or nitrogen. Similar precipitates had been earlier reported by Wiffen et al. (150 ) and Rau et al. (148 ). At 903 K, extremely mobile and large dislocation loops were observed which interacted with voids to give dislocation tangles but no planar precipitates were seen. At 1023 K, only an inhomogeneous distribution of voids existed at sites previously having dislocation loops. Shiraishi et al. (151 ) studied effects of fast neutrons on 0.3 mm thick vanadium disks at 473 K (200 • C) which was followed by annealing at temperatures of 453 to 973 K (180 -700 • C). They observed growth of dislocation loops with temperature till around 873 K (600 • C) where 1/2 111 111 vacancy loops were found to have an average diameters of 300 Å, as seen in Fig. 1.12. These dislocation loops annealed out following 1 hr of heating at 973 K (700 • C). 1.2 Literature Survey Horton et al. (152 ) analyzed temperature dependence of microstructures in vanadium due to neutron irradiation. Vanadium disk of 0.5 mm thickness was irradiated with neutrons up to 1 dpa at temperatures between 455 and 925 K. They also observed three kinds of defects, as reported by Elen et al. (149 ). First type of defects were dislocation loops. Between 455 and 625 K, the dislocation loops were shown to grow in size. The Burgers vectors of these loops at 625 K was observed to be 1/2 111 which remained the same till atleast 773 K, in accordance with Elen et al (149 ). Secondly, homogeneously distributed (except at 775 K) cavities were observed which also grew in diameter with increasing temperature. Thirdly, platelet precipitates were observed at 725 and 775 K, lying along 100 on 012. These precipitates, previously observed by Wiffen et al. (150 ), Rau et al. (148 ) and Elen et al. (149 ), are expected to form due to interaction of vanadium self-interstitials with interstitial impurities. However, exact composition of these precipitates wasn't explained till then (153 ). Bradley et al. (154 ) examined microstructural change in vanadium wires and foils on irradiation with T(d,n) neutrons at 300 K, 475 K and 675 K. At 300 K, irradiation resulted in very small defect clusters which were assumed to be 1/2 111 loops based on earlier experimental observations. At 475 K, larger interstitial 100 loops on (100) planes were observed in addition to small 1/2 111 loops (nature unknown). At 675 K, a heterogeneous distribution of dislocation loops and planar precipitates was reported. Interestingly, precipitate distributions of similar kind were also observed by vacuum annealing of unirradiated samples for 300 h at 675 K. This led them to conclude that precipitates were thermally activated, contradicting previous studies that believed precipitates were radiation-induced (149 , 152 ). Observation of 100 interstitial loop is also another contradiction to previous studies (151 , 152 ). However, Bradley et al. attributed this contradiction to the difference in impurity levels among the specimens irradiated in various studies. Little,Bullough and Wood (155 ) had predicted presence of 100 loops in vanadium and niobium using their model that was originally proposed to explain presence of 100 dislocation loops in iron. One of the requirements of the model is a high preirradiation dislocation network density which is not seen in the case of annealed vanadium, so this model seemed an unlikely explanation for vanadium. Ohnuki et al. (156 ) performed microstructural studies of 3 mm vanadium disks (pure and impure) irradiated with neutrons up to 14 dpa at 870 K. After irradiation, pure vanadium was found to contain low-density large voids (dia = 50 nm) and higher-density small voids (dia = 5 nm); blocky precipitates (of 30 nm) and large platelet-shaped precipitates (several hundred nm long, 10 nm thick), and a network of dislocation lines. Meanwhile, irradiation of impure vanadium resulted in no void formation, under-developed network of dislocations and higher density of platelets with faulted features along 001 direction. Ochiai et al. (157 ) compared microstructural evolution in pure vanadium during ion (3 MeV Cu 3+ ) irradiation at constant and varying temperature (473/873, 673/873 and 873/473 K). As reported earlier (149 , 152 , 154 ), this study also reported interstitial dislocation loops, needle-like precipitates along 100 direction and voids. Fig. 1.13a shows the variation in density of these three types of defects as a function of temperature while Fig. 1.13b shows number density change at various temperatures and irradiation conditions. Nita et al. (158 ) performed stepwise temperature increase alongwith upward and downward temperature changes to study microstructural changes due to irradiation with 4 MeV nickel ions and reported similar trends of density. Watanabe et al. (159 ,160 ) studied change in microstructure of vanadium by comparing irradiation with positron at two different temperature regimes. One temperature regime involved irradiation of vanadium disks for 8 irradiation cycles consisting of irradiation at 498 K of about 0.05 dpa and at 613 K of about 0.45 dpa whereas the second regime involved continious irradiation at 613 K. The latter resulted in an ordered array 1.2 Literature Survey of voids while the former developed thin, long platelets of precipitates in 100 directions with a random distribution of small voids. These precipitates were identified as carbides which were also detected by conducting copper ion irradiation in varying temperature condition(473/873 K) and were found to be sensitive to temperature variation. There is one important conclusion that can be drawn from various temperature dependent irradiation studies (161 -164 ) of microstructure in vanadium: Pre-irradiation of specimens at lower temperature can suppress interstitial loops formation when temperature variation exceeded a certain characteristic temperature. This temperature was estimated at about 673 K for pure vanadium (157 ).

Elen et al. (

Simulation-based results

Ab initio calculations show that a 1/2 111 interstitial dumbbell is the most stable SIA in vanadium. Various interatomic potentials have been developed to perform better atomistic simulations and its implications (165 ) using Finnis-Sinclair (FS) formalism (166 , 167 ), using modified embedded atom method (MEAM) (168 ) by refitting of parameters developed by Baskes (89 ). 170 ) reported stability of 1/2 111 SIA loops over 1/2 110 and 100 in simulations using the FS potential. They also found easy rotation of these other dislocation loops into 1/2 111 loops by relaxation at low temperatures due to a significant difference in formation energies for 1/2 111 loops and 100 or 1/2 110 loops. Also, interaction of two 1/2 111 loops indicated a thermally unstable 100 junction that rotated into 1/2 111 direction at temperatures between 327 and 800 K (527 • C). A detailed study of diffusion of SIAs in vanadium (171 , 172 ) revealed that SIA jumps were correlated with the correlation factor being sensitive to temperature changes below 800 K.

Summary for vanadium: In general, dislocation loops of Burgers vector 1/2 111 , voids and planar precipitates are observed after neutron irradiation. However, there are contradictions among the experiments with one reporting the presence of 100 interstitial loops in addition to the above-mentioned defects and another stating the nature of 1/2 111 dislocation loops as vacancy-type. These contradictions have been attributed to the varying levels of impurities among the different experimental samples. Ab initio calculations have established that the 111 interstitial dumbbell is the most stable and even predictions using interatomic potentials concur regarding the stability of 1/2 111 interstitial loops. In fact, simulations using Finnis-Sinclair potential for V predicts easy 1 State of the Art Review rotation of 100 and 110 loops into 1/2 111 due to a significant difference in their formation energies.

Overview of radiation-induced defects in tantalum Experimental observations

The database for experimental study of irradiated tantalum is extremely scarce. Among these, only a handful of them actually report any details relevant to our field of interest (with respect to microstructural changes: observation of point defects, dislocation loops etc.). The first TEM study involving irradiated tantalum was reported in 1970 by Kulcinski et al. (173 ). They analyzed neutron-irradiated tantalum foils at 723 K (450 • C) and reported a large density of small loop and defect clusters with an average size of 50 Å and no voids of comparable size. Absence of large clusters of loops is attributed to lower interaction energy (174 ) due to smaller shear modulus for tantalum when compared to bcc metal like Molybdenum. Wiffen et al. (175 ) irradiated tantalum foils with neutrons at temperatures ranging from 698 to 1323 K (425 -1050 • C). Many screw dislocations with Burgers vector (presumably) 1/2 111 were observed. Electron microscopy revealed two different categories of irradiated specimens: one after irradiation at 698 K and another for irradiation after temperatures greater to equal to 858 K (585 • C). Irradiation at 698 K resulted in mainly small dislocation loops, presumed to be of mixed vacancy and interstitial types with a few larger loops and dislocation segments. They also observed slip of 1/2 111 dislocations on {110} planes. On the other hand, specimens irradiated at 858, 1063 and 1273 K (585, 790, 1000 • C) consisted of a void population as well as a dislocation structure. The void size was found to increase with increase in temperature of irradiation. Dislocation structures did not follow a general trend as for voids. After irradiation at 858 K, 1/2 111 screw dislocations were reported along with random segments. Meanwhile irradiation at 1063 K resulted in various areas free of dislocation, connected dislocation and even isolated loops in addition to voids pinning dislocation segments. At 1273 K, connected networks at smaller scale were observed as were a few voids without evidence of void-dislocation interaction. copper ion irradiation to 20 dpa between temperatures 773 and 1546 K. Below 1073 K, a high density of small dislocation loops and line dislocations are observed. At 973 and 1073 K, they were identified to be of vacancy type and so the loops at 773 and 873 K were also presumed to be of vacancy type. They also studied the evolution of microstructure at 973 K, as shown in Fig. 1.14. Interstitial type loops having a diameter of about 10 nm were seen to be uniformy distributed at 0.03 dpa. With further irradiation, these dislocation loops combined to form larger loops and appeared to tangle with line dislocations. Small vacancy type loops formed preferentially inside these large interstitial loops. Similar evolution studies at higher temperatures (1173 and 1408 K) revealed voids even at 0.03 dpa which grew in size with increasing dose. Basically, they concluded that the stable morphology of vacancy clusters changed from 2D loops to 3D voids and the transition temperature was between 973-1073 K.

Yasunaga et al. ( 177) also studied copper ion irradition of tantalum disks up to 3 dpa for temperatures ranging from room temperature to 1073 K. They reported presence of only vacancy loops and vacancy voids.

State of the Art Review

Simulation-based results

From first-principles, it has already been established that 1/2 111 loops are the most stable. As for the other transition metals, various potentials have been developed for tantalum. EAM potentials for tantalum include those fit to experimental properties (83 , 178 , 179 ), those fit on ab initio and experiments (180 , 181 ) along with MEAM potential as in Ref. ( 182) and an angular-dependent potential as in Ref. (183 ). These have been used in previous studies related to simulations of various dislocation processes (184 , 185 ).

Summary for tantalum:Scarcity of experiments and contrasting irradiation conditions among the various experiments prevent any conclusive evidence of finding a particular type of defect over another. While one experiment found screw dislocations with Burgers vector 1/2 111 , another reports growth of interstitial loops which is followed by a preferential formation of vacancy loops at higher temperatures. Nevertheless, it has already been established from first principles that 1/2 111 loops are the most stable.

Conclusions

In this chapter, we have presented a brief summary of radiation damage and a thorough survey of the relevant published articles about radiation damage in iron, tungsten, vanadium and tantalum. In this section we emphasize the main conclusions.

Experimental observations of irradiated bcc transition metals except Fe reveal the dominance of dislocation loops with 1/2 111 Burgers vector. Irradiation of Fe at low temperatures also results in 1/2 111 loops while irradiation at high temperatures results in 100 loops. Additionally, small 100 loops have been observed in W under heavy ion irradiation at low temperatures. Explanation of these results are still under debate.

In agreement with experiments in bcc transition metals except Fe, DFT calculations predict 1/2 111 dumbbell to be the most stable mono-SIA. Although there is no experimental evidence to support their formation, the 3D C15 SIA clusters appear to be particularly stable in bcc Fe according to DFT calculations (37 ). In iron, it was found theoretically that these C15 clusters are stable, immobile, and exhibit large antiferro-1.3 Conclusions magnetic moments. Theoretically, these C15 clusters have been found to form directly inside atomic displacement cascades and are able to grow by capturing self-interstitial atoms from the surrounding material. In other bcc transition metals the DFT calculations show that the energy of very small C15 cluster (up to 4 interstitial) is always significantly higher than that of parallel 111 crowdions. In Ta this difference is the lowest. A systematic description of the energy landscape of C15 clusters and dislocation loops (1/2 111 and 100 loops) will be the subject of Chapters 3, 4 and5. Below, we recapitulate the brief summary of each material treated in this thesis:

Summary for iron:The number of experiments studying the irradiation-induced microstructural changes in iron is quite high compared to the other transition metals. All these experiments report a majority of mobile 1/2 111 interstitial loops at lower temperatures. Further, an increasing proportion of the comparatively immobile 100 is observed above 573 K which ultimately become the majority at around 773 K. The orientation of the sample is crucial in this observation, otherwise the 1/2 111 can migrate to surface, biasing the results of experiments. Moreover, the impurities can also bias the results impacting the mobilities of the different loop types. For example, He stabilizes the 1/2 111 loops rendering them immobile, as predicted by MD simulations. Cr in Fe-Cr alloys reduce the mobility of the 1/2 111 loops as well, eventhough the combined effect of He beam and Cr content in alloys did not show any clear trend. Though the theoretical explanation was lacking till recently, most of the experiments have now been explained rather convincingly in the existing literature. According to DFT Calculations, the most stable mono-SIA in Fe is the 110 dumbbell. The 110 SIAs are the most stable in Fe up to around five SIAs and 1/2 111 loops are the most stable beyond the cluster size of around five interstitials. Taking into account elastic approximation to calculate the formation energies near the αγ, a reduction of free energy of 100 loops was demonstrated, thus making the 100 loops more stable compared to 1/2 111 loops at temperatures higher than about 800 K. Although the stability of 100 was thus established, the mechanism for the formation of 100 loops has not been confirmed yet. The theoretically proposed mechanism of the formation of 100 by collisions of two 1/2 111 loops can explain the formation of 100 loops in pure Fe, or Fe under He atmosphere. However, despite the fact that Cr reduce the mobility of 1/2 111 loops the formation of 100 loops is very active in F eCr alloys. Moreover, the absence of 100 interstitial dislocation loops at low temperatures remains unexplained, given that the 1/2 111 loops are very mobile at low temperature. Recent studies highlight the importance of C15 clusters in understanding and explaining the anomalous experimental 2 Methods: Theoretical background of atomistic methods

Atomistic methods, as the name suggests, involve consideration of each atom in the material in order to understand its response under given conditions. The various phenomena taking place at the atomic scale affect material behaviour at a macroscopic scale, which in turn influences the thermal and mechanical properties of this material. Applicability, lifetime and sustainability of a material under specific conditions can only be determined accurately when a complete understanding of material behaviour at atomistic and consequently, macroscopic level is attained. This serves mainly two purposes: Firstly, limitations of experiments are circumvented and secondly, atomistic methods provide the means to verify and validate our understanding of material. Numerous methods have been developed to effectively describe the atomic scale of materials, varying in computational cost, accuracy, approximations made etc.

Although there exist a number of numerical methods to characterize the energy landscape of point defects in metals, each has its limitation. Ab initio methods, based on Density functional theory, provide accurate formation energies for small-sized interstitial defects but are computationally costly at larger sizes. To overcome this size limitation, various empirical potentials have been developed, such as, the Embedded Atom Method(EAM) potentials for Fe and W which are widely used. These EAM potentials have facilitated calculation of formation energy of defects while being computationally feasible and fast. However, they do not yield accurate values for formation energy and are non-transferable giving a difference of up to 400 eV for clusters of 1000 interstitial atoms in Fe.

Objective: This chapter provides a brief insight into the atomistic methods undertaken during this thesis. Additionally, we discuss the advantages and the limitations of these atomistic methods which range from ab initio to interatomic potentials, through tight binding approximation.

The various sections and their contents are as follows:

• Section 2.1: In this section, a chronological account of the evolution of ab initio methods is presented with subsections on the Born-Oppenheimer approximation (Sec. 2.1.1), the Density Functional Theory (Sec. 2.1.2) andBloch's theorem (Sec. 2.1.3). • Section 2.2: A brief discussion of the principles of the tight-binding method are presented to describe the transition from ab initio methods to empirical potentials.

• Section 2.3: In this section, the formulation for Embedded Atom Method potential is briefly discussed and compared with tight binding empirical potentials. This section also presents a subsection on the existing empirical potentials for bcc metals(Sec.2.3.1).

• Section 2.4: A short conclusion is also provided at the end to highlight the main points drawn from this study of atomistic methods.

Ab initio methods

At the macroscopic level (except for cases dealing with relativistic velocities), classical mechanics provides the requisite physical laws to describe and predict the motion of objects by measurement of its position as a function of time. However, the accuracy of measurement deteriorates at atomic levels. This uncertainty doesn't stem from any limitation of measurement equipments but is an inherent feature of objects at an atomic scale in the realm of quantum mechanics. Quantitatively, Heisenberg's uncertainty principle gives a lower bound to this inherent uncertainty as ∆σ x ∆σ p ≥ h/4π, where ∆σ x and ∆σ p denote standard deviations in position and momentum while h is the Planck's constant. The equality holds if the wave-packet is a Gaussian.

In the quantum realm, a particle is described by a wave function ψ( r, t) which can be 2.1 Ab initio methods obtained by solving the time-dependent Schrödinger's equation:

i ∂ψ( r, t) ∂t = - 2 2m ∇ 2 ψ( r, t) + V ( r, t)ψ( r, t), (2.1)
where = h/2π, m is the mass of the particle and V( r,t) is the potential energy associated with the particle. The physical significance of ψ( r, t) is that probability of finding the particle at a position r at time t is proportional to ψ( r, t) * ψ( r, t). Given ψ( r 1 , t) for a particle, the probability of finding it at the position r 1 at time t is:

P ( r 1 , t) = ψ * ( r 1 , t)ψ( r 1 , t) ψ * ( r, t)ψ( r, t)d r , (2.2)
where the denominator signifies integration over all space.

If potential energy is assumed to be independent of time, the wave function can be split into a time-dependent part and a spatially-dependent part as ψ( r, t) = φ( r)ϕ(t) which leads to splitting of Eq. 2.2 into time-dependent and time-independent equations as follows:

i ∂ϕ(t) ∂t = Eϕ(t), (2.3)
- 2 2m ∇ 2 φ( r) + V ( r)φ( r) = Eφ( r), (2.4)
where E is a separation constant which was later established to be the energy associated with the particle of interest. The general solution for the time-dependent equation can be written as ϕ(t) = Ce -iEt/ . The LHS of time-independent is actually the Hamiltonian operator, Ĥ: Ĥ = -

2 2m ∇ 2 + V ( r). (2.5)
Substituting the Hamiltonian operator in Eq. 2.4 gives:

Ĥφ = Eφ. (2.6)
The Hamiltonian operator acts on a wavefunction to return the wavefunction along with a constant value for energy E. So, the time-independent Schrödinger's equation is an eigenvalue problem. This implies that for a set of possible eigenfunctions φ n , only certain discrete energy eigenvalues are feasible which satisfy Ĥφ n = E n φ n .

The time-independent Schrödinger's equation holds good, except for fast processes like laser excitation or the first stages of high energy cascades. Further on, only the time-independent Schrödinger's equation will be discussed which can be written for a general many-body system with N nuclei at positions R I and n electrons at positions r i as:

Ĥφ( r i , R I ) = Eφ( r i , R I ) (2.7)
where Hamiltonian for the system is given by:

Ĥ = N I=1 P I 2 2M I + n i=1 p i 2 2m - n i=1 N I=1 Z I | r i -R I | + 1 2 N I=1 N J=1 Z I Z J | R I -R J | + 1 2 n i=1 n j=1 1 | r i -r j | ,
(2.8) where P I , M I and Z I denote the momentum, mass and charge of the i th nucleus while p i denotes the momentum of the i th electron, m and e being mass and charge of an electron. The first two terms represent the kinetic energy of nuclei and electrons while the following three terms denote the Coulombic interction between electrons and nuclei, among nuclei and among electrons, respectively. For simplicity, Eq. 2.8 can be written as:

Ĥ = T z + T e + V ez ( r, R) + V zz ( R) + V ee ( r), (2.9)
where T z and T e denote the kinetic energy of nuclei and electrons whereas V ez ( r, R), V zz ( R) and V ee ( r) denote Coulombic interaction (potential energy) between electronsnuclei, among nuclei and among electrons, respectively. This time-independent Schrödinger's equation can be solved analytically only for the case of a hydrogen atom. For a system consisting of more than two interacting particles, no analytical solutions can be found. Also, satisfying Pauli's exclusion principle and indistinguishability of fermions further complicate the unattainable. So, the way forward can be paved only on the basis of some reasonable assumptions and simplifications which will be presented in the following sections

Born-Oppenheimer approximation

The Born-Oppenheimer approximation ( 186) is based on the fact that the nuclei are much heavier than electrons implying that the timescales of nuclear motion are much larger than those of electronic motion. In other words, the electrons are able to adjust quickly in response to any nuclear motion and this nullifies any effect that the nuclear motion would have had on the system.

Although Born-Oppenheimer approximation and the adiabatic approximation are of-

Ab initio methods

ten used interchangeably, there does exist a distinction in their implication (typically, a chemist's viewpoint). Strictly speaking, the Born-Oppenheimer approximation implies absence of correlation of electronic surface and ionic motion while adiabatic approximation implies no coupling between different electronic surfaces. Finally, it can be stated that the electrons evolve into an external potential provided by the nucleus, V ext = V ez . Following Born-Oppenheimer approximation, Eq. 2.9 can be simplified as:

Ĥ = T e + V ext ( r, R) + V ee ( r), (2.10)
where the symbols follow the previously-stated convention.

Density Functional Theory

Although the Born-Oppenheimer approximation resulted in some significant simplification, the actual problem of determination of many-body wavefunctions still involved a great amount of complexity. The quantum problem is well defined if the external potential and the number of electrons are given. In principle, the wavefunctions can be determined using Schrödinger's equation and these wavefunction give access to system properties like charge density, energy etc. However, groundstate properties of the system were not practically feasible.

This led Hohenberg and Kohn (187 ) to suggest an alternative which relied solely on the electron density ρ( r) without the main obstacle of determination of wavefunctions φ( r 1 , r 2 , ..., r n ) from many-body Schrödinger's equation. Electron density of n electrons can be written:

ρ( r) = n |φ( r, r 2 , ..., r n )| 2 d r 2 ...d r n , (2.11)
They provided a unique way to calculate groundstate energy of the system as a functional of the electron density:

E[ρ( r)].

Hohenberg-Kohn Theorems

The basic idea of their theorems (187 ) was acknowledgement of the significant role that the electron density, ρ( r), plays in defining the quantum system. Basically, this caused 2 Methods: Theoretical background of atomistic methods a shift in focus from the problem of determination of wavefunctions to electron charge density. Considering an arbitrary number of electrons subjected to an external potential V ext , the Hamiltonian for this system can be written, from Eq. 2.10:

Ĥ = T e + V ee + V ext , (2.12)
where, as stated, T e and V ee denote the kinetic energy and potential energy of electrons.

Theorem 1: The external potential V ext (and consequently, the ground state energy of a system of electrons) is uniquely determined by its ground state electron density, assuming that the groundstate is non-degenerate.

Proof: Let an external potential V (1) ext and Hamiltonian H (1) with groundstate wavefunction φ (1) lead to groundstate electron density ρ with a groundstate energy of E (1) =< φ (1) |H (1) |φ (1) >. Suppose an external potential V [START_REF] Zinkle | Proceedings of the 23rd Symposium of Fusion Technology[END_REF] ext and Hamiltonian H (2) with groundstate wavefunction φ (2) also give the same groundstate electron density but with a groundstate energy E (2) =< φ (2) |H (2) |φ (2) >. Using variational principle, we obtain:

E (1) < φ (2) |H (1) |φ (2) = φ (2) |H (2) + V (1) ext -V (2) ext |φ (2) , (2.13)
which implies,

E (1) < E (2) + [V (1)
ext -V

ext ]ρ( r)d r.

(2.14)

Similar treatment with interchanged ( 1) and ( 2) gives: ) ext -V

E (2) < E (1) + [V (2
ext ]ρ( r)d r.

(2.15)

Now, addition of Eq. 2.14 and Eq. 2.15 leads to the following absurd relation:

E (1) + E (2) < E (1) + E (2) (2.16)
Thus, by reductio ad absurdum, the external potential V ext is shown to be uniquely determined by the non-degenerate groundstate electron density ρ( r). Since V ext contributes to the complete Hamiltonian Ĥ, it can be concluded that the groundstate is itself a unique functional of the groundstate electron density. As stated, this theorem is valid only for groundstate and does not hold for a system with excited electrons.

Theorem 2: A universal functional F [ρ( r)] can be defined in terms of the density for any given external potential V ext . The exact groundstate is the global minimum value of this functional.

Proof: Firstly, we define the universal function F [ρ( r)] as follows:

F [ρ( r)] =< φ|T + V |φ > .
(2.17)

One can write energy of a system subjected to an external potential V ext as a functional of density ρ (1) as follows:

E[ρ (1) ] =< φ (1) |T + V |φ (1) > + < φ (1) |V ext |φ (1) >, (2.18) having a minimum value at the groundstate wavefunction φ and groundstate electron charge density ρ such that (variational principle): 1) ] = φ (1) |T + V |φ (1) + φ (1) |V ext |φ (1) > φ|T

E[ρ (
+ V |φ + φ|V ext |φ = E[ρ], (2.19)
or equivalently, writing in terms of charge density:

E[ρ (1) ] = F [ρ (1) ] + V ext ρ (1) d r > E[ρ] = F [ρ] + V ext ρd r. (2.20)
Hence,

E 0 = E[ρ( r)] = minE[ρ (1) ( r)].
So, it can be concluded that the minimum of the energy functional of electron density indeed is the minimum of the energy of the system.

If F [ρ] were a sufficiently simple functional of ρ, a mere minimization would suffice to calculate the groundstate energy and density in a given external potential. The price to pay for this simplification is that the actual form of the functional E(ρ) is not known.

In mathematical terms, this means that a variational principle exists for E(ρ), which is minimized for the charge density of the ground state energy. However, determination of F [ρ] is fairly complex and this led to further development in the form of Kohn-Sham equations.

Kohn-Sham equations

Complexity of the quantum problem of determination of groundstate properties of a given system of electrons under a given external potential (using Schrödinger's equation) is mainly caused by the interaction of electrons among themselves i.e. the electronelectron repulsion. In order to simplify this complexity, 189 ) proposed

2 Methods: Theoretical background of atomistic methods mapping of a many-body interacting system with the actual external potential onto a non-interacting system of electrons in an effective potential such that the groundstate charge density remains the same as the original system of interacting potential. This effective potential is known as the Kohn-Sham potential. Kohn-Sham equations essentially provide a simpler means to carry out calculations because the Schrödinger-like equations just need to account for the external potential and not for all the possible interactions considered otherwise. The groundstate energy functional of interacting electrons can be written as:

E e [ρ( r)] = V ext ( r)ρ( r)d r + E KE [ρ( r)] + E H [ρ( r)] + E XC [ρ( r)], (2.21)
or equivalently,

E e [ρ( r)] = - N K=1 ρ( r) r -R K d r + n i=1 φ i ( r)( -∇ 2 2 )φ i ( r)d r + 1 2 ρ( r 1 )ρ( r 2 ) | r 1 -r 2 | d r 1 d r 2 + E XC [ρ( r)], (2.22)
where the first term gives the interaction of electrons with the external potentials (which should be determined), the second term is the kinetic energy of a system of noninteracting electrons, the third term gives the electron-electron repulsion (Hartree electrostatic energy of a charge distribution) and the fourth term is the exchange-correlation functional. This exchange-correlation function is a universal function whose form is not exactly known. Several approximations for this exchange-correlation term like Local Density Approximation (LDA) or the Generalized Gradient Approximation (GGA) are commonly used. These approximations will be discussed at the end of this section.

Minimization of the energy functional with respect to orbitals is performed to obtain the orbitals that give the groundstate energy. Although both give the same results, minimization is done with respect to φ * ( r) instead of φ( r). So, we get:

δE e δφ * i ( r) = δE KE δφ * i ( r) + δE ext δρ( r) + δE H δρ( r) + δE XC δρ( r) δρ( r) δφ * i ( r) = ε i φ i ( r), (2.23)
after functional derivatives the equation becomes:

-∇ 2 2 + V ext ( r) + d r ρ( r ) | r -r | + δE XC δρ φ i ( r) = ε i φ i ( r), (2.24)
The last equation resemble to Schrödinger's equation with a form such as:

[T + V KS ef f ]φ i ( r) = i φ i ( r), (2.25)
2.1 Ab initio methods and represents the Kohn-Sham equations. The effective self-consistent potential V KS ef f of non-interacting electrons in term of electron density reads:

V KS ef f = δ[E ext + E H + E XC ] δρ * i ( r) = V ext ( r) + d r ρ( r ) | r -r | + δE XC δρ (2.26)
In practice, we can evaluate the ground state energy of interacting electrons from the solutions of the Kohn-Sham equation Eq. 2.25 as:

E e (ρ) = i i f F D ( F -i ) + F [{R I }, ρ], (2.27)
where the electronic density can be deduced as

ρ( r) = f F D ( F -i )|φ i ( r)| 2 , f F D is
the Fermi-Dirac distribution function, F is the Fermi energy, and the residual term F is given by: .28) Finally in order to evaluate the total enbergy of the system we should add in the above equation the nucleus-nucleus interaction. The new F reads:

F [ρ] = E XC [ρ] -V KS ef f ( r)ρ( r)d r. ( 2
F [{R I }] = F [ρ] + 1 2 I =J Z I Z J | R I -R J | (2.29)
An initial guess for ρ( r) is used to calculate the effective potential V KS ef f in Eq. 2.26. This V KS ef f is, in turn, used to estimate the new density ρ( r) using Eq. 2.24 and ρ( r)

= i=1 f F D ( F -i )|φ i ( r)| 2
and so on. In this way, n individual electrons are used to find the self-consistent electron density of the original system of interacting particles. So far the only unknown part in the above Kohn-Sham self-consistent calculation are the exchange-correlation energy E XC .

Local Density Approximation(LDA): This is the simplest approximation to the unknown exchange-correlation term in the Kohn-Sham equations and considers only the local electron density contribution i.e. it depends only on the value of the electronic density at a single point. Assuming that the density of exchange-correlation energy of homogenous electron gas e xc (ρ) is known from the density of the electron gas ρ( r) at a given point r, the exchange-correlation energy can be written as a simple integration (189 ):

E LDA XC = e xc (ρ( r))ρ( r)d r (2.30)
2 Methods: Theoretical background of atomistic methods

The exact form of exchange-correlation of the homogeneous gas is unknown except at the asymptotic limits, i.e. very low and high density. Thus, it is usual to assume that xc can be written as the sum of exchange and correlation parts:

xc = c + x (2.31)
The exchange part x ( r) can be derived analytically and is given by the Dirac functional (190 ):

x (p) = -

3 4 3 π 1/3 ρ 4/3 (2.32)
In contrast, the correlation part is analytically unknown. However, accurate calculations using Quantum Monte Carlo integration have been done by Cerperley and Alder (191 ).

Many analytic interpolation of these numerical results have been proposed in the past, such as, or . Despite being derived from the homogeneous electron gas and the simplicity of LDA, this exchange correlation functional provides sufficient accuracy to address material science problems. In this thesis, we have not used the LDA approximation which is not very well adapted to metallic systems.

Generalized Gradient Approximation(GGA): A natural extension of the LDA approximation is to go beyond the local approximation in Eq. 2.30 by increasing the region of space which contributes to the exchange-correlation energy. The density ρ(r) can be formally expanded around the point r which is assumed to be the origin, as follows (189 ):

ρ(r) = ρ + ρ i ri + 1 2 ρ ij ri rj + ..., (2.33)
where ρ = ρ(0), ρ i = |∇ i ρ( r)| r=0 etc. The above development provides a series for exchange-correlation functional as:

E (0) XC = e xc (ρ( r))ρ( r)d rr (2.34) E (1) XC = f (1) [ρ( r), | ρ( r)|] ρ( r)d r (2.35) E (2) XC = f (2) [ρ( r), | ρ( r)|] 2 ρ( r)d r (2.36) E (0)
XC reproduces the LDA while E

(1)

XC gives the GGA. This approximation requires independent calculation of function of two variables, the density (ρ( r)) and the gradient of density(| ρ( r)|). One of the most reliable exhange-correlation functional of this 2.1 Ab initio methods type is the GGA-PBE or GGA- Perdew-Burke-Ernzerhof introduced in 1996 (194 ). This functional will be used for all DFT calculations performed throughout this thesis.

There is no feasible way to improve the exchange-correlation functional because the errors introduced by LDA or GGA are systematic and unavoidable. However, all other approximation presented below are controllable in the sense that some convergence parameters can be tuned in order to reduce the corresponding error.

Bloch's Theorem

Despite all the above approximations, complexity of a practical quantum problem arises from the high number of electrons that need to be considered. Bloch's theorem provides the requisite simplification for cases with periodic lattices. Bloch's theorem proves that if potential V ( r) is periodic such that V ( r + R) = V ( r) where R is a lattice vector of the periodic structure, then the wavefunction/eigenstate of an electron in this periodic potential is:

ψ j,k ( r) = u j ( r)e i k. r , (2.37)
where u j ( r) is also a periodic function with the same periodicity as V such that u j ( r + R) = u j ( r), j is the band index and k lies in the first Brillouin zone. First Brillouin zone represents a unique primitive cell in the reciprocal space. This implies that, in principle, our problem is reduced to determining the eigenstates corresponding to each possible k in the first Brillouin zone because the eigenstates are unique only in the first Brillouin zone.

Eigenvalues of the Hamiltonian vary smoothly over the Brillouin zone, so that in practice only finite number of points are needed to be treated. Consequently, any observable which should be integrated over the first Brillouin zone can be approximate as a sum of function values calculated at a finite number of points in the Brillouin zone. This finite grid is called the k-point mesh. Choosing a sufficiently dense mesh of integration points is crucial for the convergence of the results and is therefore one of the major objectives when performing convergence tests. The most used techniques to build a grid over the Brillouin zone is the which provides regular meshes to be centered at Γ for odd or even divisions, with possibility of an additional shift in both cases.

It is important to distinguish the materials which have band gap, such as insulators and semiconductors, from the metals investigated in this thesis. The latter requires denser k-meshes compared to the former. This is a consequence of the fact that the Fermi surface has a sharp interface in metals compared to the materials with gap. The convergence of the observables, such as forces or energies with the grid is an important step for any DFT calculations. Throughout this thesis, we will mention the k-mesh used each time in order to sample the Brillouin zone. For instance, in order to obtain the lattice parameter or the elastic constant of one of the materials of this thesis, the DFT calculations require only a primitive unit cell. The k-point mesh in the reciprocal space is set with a grid between 16 × 16 × 16 and 24 × 24 × 24 k points.

Basis set

A basis set is a set of known functions φ j ( r) that are used to represent orbitals of unknown analytic form ψ i ( r) using coefficients, quite similar to the role of orthogonal unit vectors in representation of any unknown vector using coefficients. Basically, the unknown solution-wavefunctions are expressed as a linear combination of known functions with expansion coefficients to facilitate computation. Typically, any unknown wavefunction can be written as:

ψ i ( r) = ∞ j c j φ j ( r).
(2.38)

Now, this basis set needs to be chosen such that they resemble the solution of oneelectron Schrödinger equation and facilitate computation. There can be various possible candidates for basis set but we mention the most commonly used ones.

Localized basis set: Localized basis set involves consideration of wavefunctions on each nucleus/atom. The wavefunctions of isolated atoms can be written with a form similar to the eigenfunctions for Hydrogen atom. Examples of localized basis set are atomic orbitals (AOs), Slater-type orbitals (STOs) and Gaussian-type orbitals (GTOs).

Plane-wave basis set: Plane wave basis set are the most suitable when dealing with periodic lattices which comply with Bloch's theorem. Major advantage is that the infinite sum can be restricted by a cut-off.

The one-electron wavefunction can be written in terms of a periodic function u j ( r). Periodic nature of u j ( r) allows us to expand as plane waves basis set in terms of the 2.2 Tight Binding Approach reciprocal lattice vectors G as

u j ( r) = G c j, G e i G. r , which gives: ψ j,k ( r) = G c j, G+ k e i( G+ k). r . (2.39)
For a periodic system, integrals in real space over the (infinitely extended) system are replaced by sums over the (finite) first Brillouin zone in reciprocal space, by virtue of Bloch's theorem. This finite number is chosen by imposing a condition on the reciprocal lattice vectors G. So, only those G are considered that allowed by the following cutoff criterion:

k 2 + G 2 2 ≤ E cut , (2.40)
where E cut is the cut-off energy used in DFT calculations and refers to the kinetic energy associated with a reciprocal vector G. Any k 1 > G will then be equivalent to k 1n G where n is an integer such that k 1n G lies in the first Brillouin zone. For the concerned metallic systems with valence orbitals as s, p and d, the minimal cut distance that we take is E cut = 350eV in our DFT calculations. In practical situations, like the k-points setting, fixing of the E cut must be done in preliminary convergence tests.

Even after considerable simplification, DFT calculations are limited by the size of the periodic box considered and are feasible only upto a few hundreds of atoms. Further approximations, like Tight-Binding Approach, can increase the accessible size-range and provide valuable information. Moreover, the Tight Binding approximation led to the formalism of empirical potential for metals.

Tight Binding Approach

The Tight Binding(TB) approximation is the intermediate approximation between ab initio methods and empirical potentials. A vast majority of interatomic potential proposed for transition metals have their roots in TB approximations. The main goal of this section is to describe the transition from ab initio method to empirical potentials. Here, I will present a brief description of the principles of the tight binding method. A detailed description of this method can be found elsewhere (92 , 94 ).

The main advantage of the tight binding approach is the accessibility to both electronic structure and total energy under the assumption that the core electrons are very localized due to tight binding. This approximation, is very similar to the quantum chemistry method of Hückel (196 ) or linear combination of atomic orbitals which turned out to be very useful in the study of molecules (197 ). The tight binding method is based on the fact that the core electrons in solids remain very localized, and thus, the discrete core levels of the atoms are only slightly broadened in the solid. The corresponding wave functions are not very different from the atomic wave functions in the vicinity of each atom. Tight binding method is also a very good approximation for valence electrons when the spatial extension of the atomic orbitals is smaller than the interatomic distance. Consequently, this method is suitable for the study of electrons arising from the narrow band (3d, 4d, 5d) of transition metals (92 , 94 ). Moreover, towards the end of 90' (93 , 198 , 199 ) it was proved that the delocalized s, p valence electrons and their hybridization with d electrons can be treated in tight binding approximation. This implied that only the valence electrons should be treated. This effectively reduces our problem from a manybody Hamiltonian to a one-electron Hamiltonian which experiences an effective potential due to all the atoms in the system. For N atoms, this one-electron Hamiltonian can be expressed as:

H T B = - 2 2m ∆ + N I=1 V ( r -R I ), (2.41)
H T B |φ α = α |φ α , (2.42)
where N I V ( r -R I ) denotes the effective potential experienced by valence electrons due to N atoms of the system. V r -R I = V r -R I ≡ V I can be interpreted as the effective spherically symmetric atomic potential centered on the I th atom. α and φ α are the eigenvalue and the eigenfunction, respectively, of the secular equation Eq.2.42. This one-electron problem -Eqs. 2.42 -can be solved in a space which is spanned by a set of atomic orbitals. Hence eigenfunction can be expressed as a linear combination of atomic orbitals localized on each atom:

|φ α = Iλ c α Iλ |Iλ , (2.43)
where |Iλ represents an atomic orbital. The angular dependence of the functions {|Iλ } has the symmetry of s, p x , p y , p z , d xy , d yz , d x 2 -y 2 and d 3z 2 -r 2 orbitals. For the transition metals, investigated in this thesis, the delocalized s and p electrons give a small contribution to the cohesive energy of the system and only the d band gives a good description 2.2 Tight Binding Approach of the system. The introduction of the sp valence orbitals becomes necessary when discussing noble metals for which the d band is filled. The s band is non-degenerate, the p band is three times degenerated and the d band is five times degenerated.

Using Eq.2.43, Hamiltonian in Eq.2.41 can be solved as:

I λ c α I λ Iλ|H T B |I λ = α I λ c α I λ Iλ|I λ , (2.44)
which gives:

I λ c α I λ Iλ|H T B |I λ = α c α Iλ .
(2.45)

The Hamiltonian matrix elements Iλ|H T B |I λ are parametrized using various approximations and are related to the hopping integrals β Iλ,I λ which reflects the ability of the electrons to jump from the orbital-site state |Iλ to orbital-site state |I λ . So, the performance of tight binding models is dependent on the Hamiltonian matrix elements in the atomic orbital basis. The most common way to express the above matrix was presented in 1954 by Slater and Koster who showed that all the off-diagonal elements of an spd base can be constructed from 10 independent parameters (200 ). The ability of the tight binding model to reproduce electronic effects as well as the quality of the transferability of the model are direct consequences of the quality of parametrization of the Slater-Koster parameters.

For the case of transition metals lying in the middle of series and having narrow d band, the tight binding model can be even more simpler. The narrow d band of transition metals can be explained by the small overlap of the d atomic orbitals. This d band is superposed upon an s band which is wide, due to the large overlap of the s atomic orbitals. These two bands slightly interact and this sd hybridization may be treated as a perturbation. This is what is called a pure d band model, without any contributions from the s electrons. The atomic orbital basis {|Iλ >} (Λ = 5 and λ = 1, ..., 5 labels respectively the |xy , |yz , |zx , |x 2y 2 and |3z 2r 2 orbitals) is, in general, considered orthonormal.

The origin of the density of states of energy is taken to be the center of gravity of the d-band. In this TB approximation, it is usual to divide the cohesive energy of the metal into an attractive and a repulsive contribution. Consequently, the energy of the system can be written: .46) In this model the band term in Eq. 2.46 is considered attractive (because of the delocalization of the valence d electrons when forming the solid) and the remainder of the total energy in Eq. 2.46, E rep , can be interpreted as repulsive core-core energy. The above form of the total energy is similar to the DFT energy given by Eq. 2.28. However, there are two points that need attention: first, the self-consistent character of the DFT method is lost in the tight binding method, and second, the form Eq. 2.46) of the total energy is not universal.

E T B R I = α α f F D ( F -α ) + E rep R I . ( 2
For a given material, a good set of parameters can be found by fitting the on-site terms, the hopping integrals and the phenomelogical repulsion potential to the experimental data or ab-initio calculations. Then, the transferability of the model is checked. Depending on the ingredients included in the form of the parameters, we can have major differences between the models. For example, the non self-consistent aspect of the tight binding method becomes visible in inhomogeneous systems such as metallic surfaces. When such systems are treated, an unphysical large charge transfer is present. This effect can be avoided by introducing an additional potential determined self-consistently in order to recover a local quasi-charge-neutrality (198 , 201 ).

The energetics of transition metals is dominated by the width, shape and occupation of the d-band. The widths and the shape of the electronic density of states of the system can be calculated using the low-order theory of moments that will allow us to find a simple analytic form for the energy components in Eq. (2.46). Friedel (202 ) proposed to focus on the global and local density of states, ρ I and ρ Iλ , respectively, of atom I and orbital λ, which can be easily written as:

ρ I (E) = 1 Λ λ ρ Iλ (E) = 1 Λ λ α |c α Iλ | 2 δ(E -E α ), (2.47)
where c α Iλ is defined by Eq.2.43 with ρ I (E) and ρ Iλ (E) being normalized to unity. He noted that all the physics of the system is revealed if the local densities of state is known. The direct solution of the secular equation by a full diagonalization, Eq. 2.45, can be replaced by a careful characterization of local density of states. The theory of moments provides a historic way to characterize the local density of state.

Let me recapitulate the definition and the interpretation of moments for the benefit of the reader. Given a function f (x), the quantity µ (n) is given by:

µ (n) = ∞ -∞ (x -x g ) n f (x)dx (2.48)
where n) is known as the centred moment of order n of the f (x) distribution and x g is the center of gravity of the function. From the mathematical point of view, the moments can be used to describe the "center of gravity" of the distribution, the width, the asymmetry, the tendency of a gap, and so on, of the function f (x). Precisely:

x g = 1 µ (0) ∞ -∞ xf (x)dx . (2.49) µ (
• µ (0) gives the area under f (x). This can be used as a normalization factor of the function f (x);

• µ (1) gives the "center of gravity" of f (x); the definition (2.48) assigns it to 0;

• µ (2) gives the momentum of "inertia" of f (x). Hence µ (2) is proportional to the width of the f (x) in the root mean square sense;

• µ (3) gives the asymmetry from 0, e.g. a large negative value of µ (3) corresponds to a long tail of f (x) in the region below µ (1) ;

• µ (4) measures the tendency for a gap in the middle of the band ... and one can continue the description of the "picture" of f (x) with higher order moments.

According to the Ducastelle-Cyrot-Lackmann theorem (203 , 204 ), the n th moment is equal to the sum over all possible path of n jumps with the same starting and ending orbital-site state |Iλ . The n th moment can be expressed as:

µ (n) I = 1 Λ λ Iλ|H n T B |Iλ (2.50) µ (n) I = 1 Λ λ I 1 λ 1 ... I n-1 λ n-1 β Iλ;I 1 λ 1 β I 1 λ 1 ;I 2 λ 2 β I 2 λ 2 ;I 3 λ 3 ....β I n-1 λ n-1 ;Iλ (2.51)
where β Iλ;I 1 λ 1 is the above introduced hopping parameter permitting the jump of an electron from the orbital-site state Iλ to orbital-site state I 1 λ 1 . If the second moment of the local density of states is evaluated in terms of the Hamiltonian matrix, then the 2 Methods: Theoretical background of atomistic methods second moment can be expressed as:

µ (2) = 1 Λ I 1 =I T r[B 2 ( R I 1 I )], (2.52)
where B( R I 1 I ) denotes the matrix for hopping integrals, deduced from the Eq.2.51.

Since the energy of occupied d-states E band is proportional to the width of the d-band µ

I (205 ), we obtain:

E band ∝ A I µ (2)
I , (2.53)

E band = A I I 1 =I g(R I 1 I ), (2.54)
where g(R

I 1 I ) = T r[B 2 (R I 1 I
)] and A is proportionality constant which depends on the occupation of the d-band. Within the rectangular d-band approximation it can be demonstrated that the proportionality constant is 0 for noble metals (in which the dband is filled). Consequently, there is no theoretical justification for using the above analytical form of the E band energy for noble metals. However, this type of semi-empirical potential was used successfully for noble metals, such as Cu, Ar or Au, by fitting the A constant and the form of the g function in order to match the experimental values of the cohesive energy, elastic constants etc.

For transition metals, the above formalism is the more suitable and simplified analytic form of the total energy that can be proposed. The total energy can be divided into an attractive contribution from valence d-electrons and a repulsive contribution from the nuclei similar to DFT, Eq.2.28, or Tight-binding form, Eq. 2.46, as:

E({R I }) = -A I I 1 =I g(R I 1 I ) + F [{R I }], (2.55)
which is the form of semi-empirical potentials using the second-moment tight binding approximation.

EAM potentials

The goal of any atomic potential is to represent the potential energy as a simple function of atomic coordinates. As a consequence, the electronic effects are integrated into

EAM potentials

a parametric form of the inter-atomic potential which is determined by a fit. This strategy greatly simplifies numerical implementation and drastically reduces the computation time for large systems compared to quantum chemical methods presented in the precedent sections. The tight binding empirical potentials given by the Eq. 2.55 answer to this requirement and the total energy can be decomposed into two terms: the first term is the sum of pair functions which accounts for the repulsion between nuclei and the second term is a many body function(interaction between three or more atoms) which mimics the cohesive d band energy.

In the begining of 1980, Daw and Baskes(206 ,207 ) proposed the embedded-atom method (EAM) wherein an atom is assumed to be embedded into the local electron density due to all the remaining atoms. They begin with the observation that the total energy of the system using DFT can be written as following from Eq. 2.29 and Eq. 2.21:

E [{R I }] = E KE [ρ] + E XC [ρ] + E H [ρ] + V ext ( r)ρ( r)d r + 1 2 I =J Z I Z J | R I -R J | . (2.56)
The first and the second terms of the above equation can be grouped together in a functional G[ρ] and the following assumption are made: (i) G(ρ) can be described by G[ρ] = g(ρ( r), ∇ρ( r), ∇ 2 ρ( r), ...)dr where g is assumed to be a function of the local electron density and its lower derivatives. (ii) The electron density ρ is assumed to be a linear superposition of the densities of individual atoms. Using these two assumptions, Daw and Baskes found the total energy of the system to be divided into a pairwise contribution which accounts for the electrostatic contribution and a many-body attractive local contribution:

U = N -1 I=1 N J=I+1 ϕ(R IJ ) + N I=1 Φ(ρ I ), (2.57)
where N is the total number of atoms in the system, R ij is the distance between atom I and atom J. Φ(ρ I ) is the embedding energy and ρ I is a function which represents the electronic density around the I th atom.

The second moment many-body potentials approximate the spd in tight-binding approximation while the EAM potentials are based on local approximation of atoms emerged in nearly free electron gas. It is interesting to note that even with evidently different approximations to start with, the cohesive term in both approaches involves many-body while the repulsive term involves a pair dependence. On comparison, it is found that the second moment potentials resemble a particular case of the EAM potentials.

2 Methods: Theoretical background of atomistic methods

In practical implementation, it is useful to give a simple form to the functions involved in Eq. 2.57. Now, ρ i = j ψ(R ij ) and the three functions ϕ(r), ψ(r) and Φ(r) can be written as summations over basis functions. These basis functions can be chosen to be cubic splines, exponentials or polynomials to obtain different variants of the EAM potentials. As an example, the widely-used EAM potentials for iron, developed by use cubic spline basis functions for ϕ(r) and ρ(r) whereas Φ is (generally) expanded as Φ = a 1 √ x + a 2 x 2 . We note that the only variance compared to second-moment tight binding approximation is due to the x 2 term.

Empirical potentials for bcc metals

As we have mentioned in the Introduction to this thesis, the main materials of interest in the nuclear fission industry are the body centered cubic (bcc) metals. This endorses continuous development of empirical potential for those metals in order to study of point defects, dislocation loops or extended defects such as dislocation or surfaces. As such, these studies serve as the basic foundation for future research regarding structural materials in the nuclear industry. E.g. molecular dynamics simulation up to microsecond with simulation box containing millions of atoms lead to better understanding of the behavior of structural materials exposed to radiation.

Success of molecular dynamics simulations depends significantly on the inter-atomic potential used. Hence, it becomes important to wisely consider the stakes before imposing approximations. Generally, this choice can made according to the system of study. The EAM formalism (as described in the previous section) works well for bcc metals (83 , 85 ) because of relatively simple shape of the electronic density of states (DOS). However, there exist other alternatives like Baskes' modified EAM (MEAM) (89 ), bond order potential (BOP) (91 , 93 ), or tight-binding method (94 )(described earlier in this chapter) for more accurate description of a possibly complicated DOS. Various BOP parametrizations have been developed for Fe and W (144 , 208 -210 ) but EAM potentials are the most commonly-used potentials due to their sufficiently good description coupled with ease of computational applicability.

Concerning the state of the art review of potentials for the materials investigates in this thesis we already have mentioned the situation in W in Sec. 1.2.3. In the next section we will focus on iron, the most important metal among those studies in this 2.4 Conclusions thesis due to its technological importance.

Conclusions

In this chapter we have reviewed the various atomistic methods which were used in the current study of defects, such as dislocation loops or C15 clusters, in metals.

Ideally, density functional theory (DFT) calculations would be utilised to capture the complete physics involved in a system with irradiation-induced defects. However, unachievable computational requirements of DFT calculations have fueled the search for alternatives considering reasonable approximations over the past decades. This search has led to development of a number of empirical potentials, ranging from pair potentials to EAM. Although these empirical potentials have been successful in making radiation damage studies feasible, inconsistency of results from different empirical potentials is a major shortcoming that hinders conclusive theoretical results.

In the following chapter, we present an application of the EAM potential by developing a new empirical potential for Fe, mainly to study irradiation-induced defects such as selfinterstitial atom clusters or dislocation loops.

Development of Embedded-Atom

Method potentials for defects in iron

The EAM potential formalism, described in the previous chapter, is an important tool to explore the mesoscopic regime between the two asymptotic limits of small clusters and observables dislocation loops. At this scale, the EAM formalism remains the only accessible atomistic method. This formalism is used in this chapter to develop a new EAM potential. Developing new EAM potentials using up-to-date ab initio and experimental data for all the metals of interest could potentially improve the quality of the accessible simulation-based radiation damage studies. This development is not a easy task, due to limitation of time, we develop an EAM potential only for for one metal, iron. For technological interest, we have chosen iron, the most important metal among those studied in this thesis.

As discussed in Sec. 1.2.1, C15 clusters have been reported recently using DFT calculations. Consideration of this new DFT information is important in development of a new EAM potential for Fe. Thus, we develop a new C15 oriented EAM potential in this chapter with the main objective of investigating the energy landscape of C15 clusters.

As pointed in Chapter 1, radiation damage studies in Fe have posed a perplexing scenario for a while now due to experimental observation of both 1/2 111 and 100 families of SIAs depending on temperature, unlike other bcc metals. In addition, DFT provides theoretical evidence (37 ) for the existence of highly stable, small-sized threedimensional SIA clusters with C15-type structure. In iron, the energy landscape of selfinterstitials atoms organized in loops is relatively known and has been widely investigated by various methods. However, the energy landscape of C15 clusters is not very well known and is very complicated due to an enormous number of possible configurations. Only a few recent studies investigate this problem. Systematic exploration of the energy landscape in search of the minimum energy C15 configurations gives good results for small clusters. Marinica et al. (37 ) used the Activation Relaxation Technique (211 ) for finding the lowest energy configurations for 2, 3 or 4 SIAs clusters. Nonetheless, the number of possible configurations grow exponentially with the size of the cluster, making a systematic search prohibitive at larger sizes. More advanced techniques using a genetic algorithm were proposed, making it possible to find the lowest energy configurations containing up to 10 SIAs (212 ). This section will solve this problem in a more pragmatic approach, we will deduce three selection rules which will be established from observation of the formation energies of several trial configurations generated using existing and newly-developed EAM potentials.

Objective: The goal of this chapter is to develop a potential for iron which includes the latest knowledge from the DFT scale of the energy landscape of C15 clusters. The motivation of such development is two-fold: (i) The newly-developed potential will be used to exhaustively explore the energetic landscape of C15 clusters up to clusters which contain hundreds of atoms. (ii) As discussed in Chapter 1, a mechanism has recently been proposed to explain the anomalous behaviour in Fe of forming 100 loops, through C15 clusters as intermediate reaction product (112 ). At very large sizes, we expect that the dislocation loops becomes more stable than the C15 cluster. The crossover between the two categories of clusters has practical importance. The confirmation of the Zhang's mechanism (112 ) depends on the crossover between the dislocation loops and the C15 clusters. A reliable potential for the energy landscape of C15 cluster is a key ingredient for finding the right crossover.

The various sections and their contents are as follows:

• Section 3.1: In this section, the existing empirical potentials for irradiation-induced defects in iron are discussed to emphasize the need for a new empirical potential.

• Section 3.2: A brief discussion of the fitting strategy of existing empirical potentials is used to strategize the development of new empirical potentials.

• Section 3.3: In this section, the fitting procedure for new Fe potentials is provided along with the database in Sec. 3.3.1 and the minimization of the cost function in Sec. 3.3.2 . On comparison with the existing Fe potentials, these new Fe potentials are found to be reasonably good. However, despite similar condition of fitting, some scatter is observed in the results which numerically demonstrate the limit of this fitting procedure. Nevertheless, the physics of the addressed system 3.1 Existing empirical potentials for defects in iron seems to be better described qualitatively.

• Section 3.4: Systematic exploration of the energy landscape in search of the minimum energy C15 configurations gives good results for small clusters. Nonetheless, the number of possible configurations grow exponentially with the size of the cluster, making a systematic search prohibitive at larger sizes. In this section, the C15 configuration is discussed as a continuation of Sec.1.2.1 and three selection rules for construction of C15 configurations are presented in Sec. 3.4.1.

• Section 3.5: The relative stability of C15 clusters and dislocation loops is analyzed in the light of the newly developed empirical potentials.

• Section 3.6: A conclusion is also provided at the end to highlight the main points drawn from this study of existing and newly-developed empirical potentials.

Existing empirical potentials for defects in iron

Performance and transferability of metal potential parameterization are obviously related to the underlying physical model. Commonly used many-body interatomic models of metals span from second moment of tight binding approximation (82 , 85 , 86 ) to EAM (87 , 88 ) or Modified EAM (89 ) and higher order tight binding models (90 -95 ). Some examples of different parametrizations for iron include a well known tight-binding second moment Finnis-Sinclair potential (83 ), the long range version given by Sutton and Chen ( 84) and the parametrization proposed in 1997 by Ackland et al. (213 ) (hereafter called A97). The force-matching method [START_REF] Ercolessi | [END_REF] was used to parametrize EAM potentials suitable for defects in bcc metals such as iron and tungsten. Using the same fitting approach, Mendelev et al (97 , 98 ) and Ackland et al (99 ) (A05) proposed a parametrization of iron. Using a different database (point defect oriented), Marinica (37 , 100 ) (M10) obtained an improved parameterization for point defects. Recently, an empirical potential potential was developed to study thermally activated glide of dislocations (214 ), this potential will be denoted as P12.

Our main goal in this chapter is to develop a potential which gives reasonably accurate formation energy values of C15 interstitial clusters relative to traditional planar clusters. Hereafter, these SIA clusters will be denoted as I conf n , where n is the net number of SIAs i.e. the number of additional atoms in the bcc lattice and conf denotes a particular configuration of interstitial cluster. As discussed in Sec. 1.2.1, the morphology of small interstitial clusters has been widely studied in the past and can be classified into two main classes: the parallel configurations formed by dumbbells packed together in bundles along 110 direction, I 110 n or 111 direction, I 111 n and secondly, the non-parallel configurations where dumbbells are not aligned. Considering the most stable parallel configurations, the DFT calculations predicted that the orientation of these dumbbells changes from 110 to 111 above about five SIAs in Fe. This holds for most empirical potential except for the potentials fit prior to Mendelev 2003 EAM potential (97 ) and a few exceptions such as the P 12 potential. The non-parallel configurations were observed for the first time in molecular dynamics simulation (36 ) using A05 potential. Due to very high vibrational entropy (100 , 110 ), these configurations are stabilized by increasing the temperature. DFT calculations showed that the triangular configuration for di-interstitial, denoted hereafter as I gao configuration. Using DFT calculations, the formation energy of triangle configuration is found to be even lower than the parallel di-interstitial I 110 2 with 0.11 eV (36 ). Moreover, the formation energy of these highly entropic non-parallel configurations, ring, tri-and tetra-interstitial configurations (for tetra-interstitial see the Figure 1.5) are very close to the formation energy of the parallel clusters. Quite recently, systematic exploration of the energy landscape using an Eigenvector following method, ARTn by Marinica et al. (110 ) revealed a very complex energy landscape of small interstitial cluster presenting a lot of non-parallel configurations. Coupling this systematic search to lattice dynamics free energy calculations, Marinica et al. (37 ) were able to reach regions of phase space which remained inaccessible by standard molecular techniques of investigation and subsequently, they provided evidence for the C15 clusters (37 ). Moreover, they showed that in the case of I 4 the ring configuration share the same energetic basin with the 3D C15 clusters using disconnectivity graph technique. These C15 clusters can be seen as generalizations of the oldest but low energy non-parallel configurations. The building block of the C15 clusters is the di-interstitial configuration, I C15 Fig. 3.1, the formation energies of the C15 clusters relative to that of the lowest energy parallel dumbbell configurations (37 , 41 , 107 ) are plotted as a function of cluster size. We have restricted the discussion to the most-widely used iron potentials in the last few years namely A97, A05, M10 and P12 potentials. However, the self-interstitial energy landscape provided by the P12 potential exhibits obvious anomalies making it unsuitable for studies concerning point defects. As seen in Fig. 3.1, the relative stability of C15 clusters with respect to the traditional clusters calculated using the M10 potential corresponds well while the A05 or older A97 potential do not agree very well. Despite its weakness with respect to point defect energy landscape, the Mendelev potentials (like A05) serve as reference for the study of extended defects, such as screw dislocations in iron. Until 2012, the 2003 Mendelev potential was the only existing EAM potential in literature that successfully predicted the compact non-degenerate core for the screw dislocation in agreement with the ab initio calculation and the glide plane of the screw dislocation in the plane in agreement with the experiments and the ab initio calculations.

However, all the Mendelev based potentials (97 -99 ) exhibit two maxima for Peierls barrier potential with a marked minimum mid-way whereas ab initio calculations indicate only one maximum, as seen in Fig 3 .2. The M10 potential, developed from the 2004 Ackland-Mendelev (A05) potential, gives very good energetics for point defects when compared with DFT calculations but the properties of screw dislocation energetic landscape are very poor in comparison. Despite the fact that the lowest configuration of the dislocation core of the screw dislocation is degenerate, the Peierls barrier potential displays an even more stable configuration mid-way between the degenerate configurations. This situation is unphysical and needs to be corrected. While the P12 potential alleviates this problem with a Peierls potential very close to the ab initio calculations of Ventelon et al. (216 ), point defect properties from P12 potential show poor correspondence and match older iron potentials like A97 better than modern potentials. So, each of the existing empirical potentials seem to have some drawback which prevents it from further use in the study of SIAs in iron.

Thus, it is necessary to develop a new empirical potential which works reasonably well to overcome the shortcomings posed by the existing empirical potentials and to facilitate the study of self-interstitial atom defects in iron. In the following section, we discuss the fitting strategies adopted by the most widely-used empirical potentials and plan a strategy for fitting a new empirical potential for iron in the sections to follow.

Lessons learned from evaluation of fitting strategy of existing interatomic potentials

As pointed out in Marinica et al.(135 ), the observables included in the database play a very important role in the parameter fitting. It is therefore recommended to simulate properties close to those that were fitted. The most straightforward approach to fit the potential only on experimental observables of the crystal such as: elastic constants, cohesive energy, Debye temperature etc. is insufficient. Some examples of different parametrizations for iron include a well known tight-binding second moment Finnis-Sinclair potential (83 ), the long range version given by Sutton and Chen( 84) and the parametrization proposed in 1997 by Ackland et al.(213 ). Although these potentials fit bulk properties very well, they fail to reproduce simple properties of point defects such as the most stable configuration of the single interstitial. In iron, all these parametrizations predict 111 configuration to be the most stable instead of 110 , as attested by DFT (32 , 34 , 107 ) and the experiment ( 25). This contradiction is not surprising since the goal of any atomic potential is to represent the potential energy as a function of atomic coordinates by integrating electronic effects into a parametric form which is determined by a fit. The environment of atoms which form the 110 mono interstitial dumbbell in iron is very different from the perfect bulk atom neighborhood: the length of the dumbbell is around 2 Å compared to 2.5 Å of the first nearest neighbor distance in bcc iron. Consequently, this change in the neighborhood results in a strong modification of the electronic density around interstitial atoms which is also reflected in the drastically changing of magnetism around the dumbbell. Such sharp effects cannot be caught by the physics injected in the tight binding second moment potentials. This lack of versatility of the physical model should be compensated by directly fitting a dumbbell environment.

More accurate potentials need to be developed using larger fitting data base which contain information beyond experimental quantities e.g. the total energy or forces acting on atoms in various configurations provided by ab initio methods. The Force-Matching method proposed by Ercolessi and Adams in the late 90s [START_REF] Ercolessi | [END_REF] showed for the first time that the database used in the fit of the empirical potentials play the same key role as the formalism used to mimic the electronic effects. The force-matching method was used to parametrize EAM potentials suitable for defects in bcc metals such as iron and tungsten. Using the same fitting approach, Mendelev et al (97 , 98 ) and Ackland et 3.3 Developing new Fe potentials al (99 ) proposed a parametrization of iron. Using a different database (point defect oriented), Marinica (37 , 100 ) obtained an improved parameterization for point defects. In 2005, Dudarev and Derlet designed an alternative approach (101 ) with a magnetic potential which was also based on the EAM formalism. By systematically enlarging the fitting database, these magnetic potential were continuously improved.

The main difference in the parametrization of the M10 and the A05 potentials is the database used. The M10 potential is parametrized on point defect related database including configurations of di-, tri-or tetra-interstitials provided by DFT calculations (36 , 107 ). Among these configurations, Marinica et al.( 100) also included the triangular configuration. In the procedure of fitting, only those parametrizations were selected for which the difference between the formation energy of triangle configuration, I gao 2 , and the parallel dumbbell configuration, I 110 2

, is low. This approach did not succeed in attaining triangular configurations with lower energy but managed to considerably reduce the difference from 0.31 eV using A05 potential to 0.07 eV using M10 potential (see the Table 3.1) (36 , 100 ). Perhaps, this low formation energy of the triangular configuration which serves as the building block for C15 clusters influenced the low formation energy of C15 clusters. In this new development, we will pay special attention to triangular configuration which will be included in the database in order to produce potentials with low formation energy for C15 clusters.

In the following section, we will develop new potentials for iron that would agree with the latest ab initio knowledge of energy landscape of radiation-induced point defects and dislocation loops, and Peierls barriers of screw dislocations. We follow a force matching fitting procedure similar to that applied by , Proville et al. (214 ) for iron and Marinica et al (135 ) for tungsten in order to produce a new set of empirical potentials for iron.

Developing new Fe potentials

Overall, the main fitting procedure for developing new potentials is similar to the one used in developing Fe (100 ) or W (135 ) potentials. As stated in (135 ), we combine the database containing three components in order to produce suitable potential parametrizations for radiation defects: experimental observables (elastic constants, cohesive energies, surface energies etc.), point defect related configurations and force match-3 Development of Embedded-Atom Method potentials for defects in iron ing matching method on the liquid or random iron configurations far from equilibrium. The details about production of database are described in the following subsection 3.3.1. The parametrization will be done using an on-the-fly iterative minimization of cost function, Sec. 3.3.2, until convergence to obtain several converged sets of parameters (the details are given also in (135 ). Details of the fitting and the minimization of the cost function can be found in the following subsection.

Fitting procedure: Database

The database contains the experimental observables (elastic constants, cohesive energies, surface energies etc.) and different minima configurations (self-interstitials and vacancies). We use force matching method on the liquid or random iron configurations far from equilibrium. The ab initio liquid information appears to play an important role in the fitting data base for extended defects.

The database used contains three components:

1. Experimental values of the solid state properties: lattice parameters and cohesive energies of the FCC/BCC structures and elastic constants C11, C12 and C44.

2. The ab initio formation energies of basic point defects in iron: the mono-interstitial with different orientations ( 110 , 111 , 100 , octahedral and tetrahedral) and the mono-vacancy. These ab initio calculations were performed within the Density Functional Theory (DFT) framework using SIESTA code, i.e. using the pseudopotential approximation and localized basis setsmade of 10 localized functions. The defect calculations were performed using the supercell approximation keeping the cell geometry fixed to the bulk equilibrium geometry and relaxing the atomic positions. The 6×6×6, 4×4×4 and 3×3×3 shifted k-pont grids were used in the 54, 128 and 250 atom cells, respectively. The Hermite-Gauss scheme for electronic density of state broadening was used with a smearing of 0.3 eV and the residual forces were smaller than 0.01eV/A. The standard Perdew-Burke-Ernzerhof Generalized Gradient Approximation (GGA) was used for exchange-correlation functional.

3. The ab initio forces acting on the atoms in the liquid or random state configurations. All the random configurations were generated using the Ackland-Mendelev potential for iron. From the initial perfect BCC supercell with 128 atoms, 15 atoms were removed in order to obtain the experimental liquid density of iron. Afterwards, this system state was propagated using molecular dynamics at temperature of 4000K over one million molecular dynamics steps (the integration step is set to 1fs). From this molecular dynamics trajectory, we extracted one liquid configuration each 50 ps. The time interval between snapshots was large enough to avoid correlation due to the molecular dynamics propagation. Finally, all 20 random configurations were used as input for the ab initio calculations. The positions of atoms were frozen and only the atomic forces acting on each atom in each liquid configuration are computed. The ab initio calculations were performed using SIESTA GGA with the same pseudo potential and basis sets used for defect calculations. Including different ab initio liquid configurations in the data base fitting, the parameterization of the potential contains information out of equilibrium (non-zero forces). Moreover, the objective function is sampled through the pair and embedded part of the EAM potential in continuum range of atom-atom distances in variance with the case when only minima are present in the database.

For example in the case of the vacancy, the radial distribution functions are not very different from the perfect bulk.

Fitting procedure: Minimization of the cost function

We chose as fitting tool, the ASSIMPOT code(217 ) where based on the principle of the variational assimilation and the adjoint model(218 , 219 ), one compares the results provided by the EAM model with the benchmark data (given experimental or ab initio) on all the configurations and then minimizes the deviation. This deviation is quantified using an objective function. For any set of parameters A the objective function is written as:

J(A) = obs w obs [Y obs -Y (A)] 2 + i w i [F i a-initio -F i (A)] 2
Y obs contains the observables to be fitted (such as, cohesive energy, elastic constants etc.) and Y (A) the corresponding values using the set of parameters A of the EAM potential. F i ab-initio is the target ab initio force and F i (A) is the corresponding force on the ith atom with A parameters. The w coefficients are the weights and are user defined, the choice of their values is to be discussed later.

Development of Embedded-Atom Method potentials for defects in iron

The cost function J(A) is minimized by use of an iterative algorithm. We start with a first set of parameters A and in each iteration step an improved vector of parameters is sought. The search direction is computed from the gradient of the cost function with respect to the parameters. The gradient needs the estimation of non linear target function and implies numerical instabilities. For this reason, the gradient is computed within the adjoint model. Alternatively, the gradient vector could be approximated by finite differences. The use of the adjoint model in the estimation of the gradient has two advantages over finite differences: high precision of the estimation and high computational efficiency.

In practice only one minimization of the cost function J(A) is not satisfactory. It turns out that it is impossible to fit energies and forces simultaneously: if one relaxes the positions with the obtained set of parameters A of EAM potential, relaxed energies deviate from their objective values. In other words, DFT-GGA and EAM force fields around defects differ significantly. This type of fitting is therefore used only in the first step. Then the atomic positions are no longer considered in the fit and for every minimization configuration, the energies are calculated using the atomic positions relaxed corresponding to the previous set of parameters. This procedure is iterated until convergence to obtain several converged sets of parameters. Of the many potentials developed, two potentials (EAM 2/3) seem to be agreeing well for most of the tests performed. These EAM potentials are quite good for SIAs (as seen in figure 3.4), bulk properties (see the table 3.1) and they exhibit one peak for Peierls barrier (figure 3.3). In particular, EAM 3 emerged as the most promising EAM potential for studying SIA defects.

Starting from the P 12 potential as an initial guess of parameters, we have developed new potentials for iron following exactly the same three stages described in (135 ). The difference here is that we selected only those parametrizations from the few hundreds of parametrizations which resulted in a physical Peierls potential at the end with qualitative agreement with DFT calculations and low disparity in the formation energy of I Two of the developed potentials, hereafter called EAM2 and EAM3, give satisfactory results for most tests. Firstly, correspondance with the bulk properties like the elastic constants is improved compared to the starting P 12 potential as it is seen from Tab. 3.1. Secondly, the energy landscape of small interstitial clusters is better matched. The mono SIA relative energies of 110 and 111 configurations of the new developed potentials are in agreement with DFT studies (32 , 34 , 107 ). The I configurations is fairly reduced. Thirdly, as seen in Fig. 3.3, the Peierls pathway of the screw dislocation exhibits one peak in qualitative agreement with the DFT calculations. However, it can be noted that having reasonable properties for the energetics of selfinterstitials impacts the quantitative agreement between DFT and the newly-developed potentials for the Peierls barrier. By comparison, the older P 12 potential is found to be closer to DFT calculation than the new potentials. Nevertheless, we recommend the new potentials because they integrate more physics than the older P 12 or M 10 potentials.

Finally, the new potentials are tested for the relative stability of C15 cluster compared to the traditional clusters. As presented in Fig. 3.4, the new potentials exhibit reasonable agreement with the DFT calculations. Although these new potentials don't compare as well as the M 10 potential, the energy landscape is quite accepable compared to the previous parametrizations.

In the following section, the C15 configuration is discussed as a continuation of Sec.1.2.1 and three selection rules for construction of C15 configurations are also presented in Sec.3.4.1.

Construction of the stable C15 clusters

Here we recall the construction of C15 clusters, briefly presented in Sec. 1.2.1. The building block of C15 clusters is a di-interstitial cluster. A simple way to insert a diinterstitial C15 cluster into a bcc matrix is to place a Z16 Frank-Kasper polyhedron having 12 atoms at the interstitial positions (see Fig. 3.5a) together with 10 vacancies around a given bcc atomic site. Larger C15 clusters can be described as sums of Z16 Frank-Kasper polyhedra having centers situated on a diamond network, which underlies the initial bcc structure (this network is shown in Fig. 3.5a-e). By repeatedly following the above-mentioned procedure of Z16 polyhedron addition, 3D clusters having a cubic periodic structure can be built. It is interesting to note that the subsequent structure of the center of Z16 polyhedra is same as the diamond structure obtained from the initial bcc structure by removing half of the initial bcc sites. This network is represented in the Fig. 3.6a. The final cubic unit cell of the crystallographic structure obtained after the 3D growth process is represented in the Figure 3 b (3nn) (eV) -0.02 -0.03 -0.03 -0.14 -0.03 -0.02 Table 3.1: Bulk, mono-and di-interstitial and vacancies properties for a few widelyused potential for iron, M10 (37 , 100 ), A05 (99 ) and the P12 (214 ) potential as well as the new developed potential marked as EAM2 and EAM3. The target value are computed from the ab-intio calculations, the same as were used in the development of the potential M10 and P12 (100 , 214 ). The a 0 and E coh denote the lattice parameter of the cubic corresponding structure and the cohesive energy, respectively. E f and E v is the formation energies of various orientation or configuration (the 1nn, 2nn and 3nn denotes the first, second and third nearest-neighbor configurations of the di-vacancy).

The ∆E

I tri 2 -I <110>
3.4 Construction of the stable C15 clusters phase or C15 structure. In this homo-atomic type of C15 structure, interstitials occupy the Cu sites and half of the original bcc sites are empty while the others are occupied and correspond to the Mg sites.

Selection rules for C15 construction

The energy landscape of C15 clusters is not very well known and is very complicated due to an enormous number of possible configurations. The number of possible configurations of n SIA clusters in vaccuum (as metallic clusters, molecules, proteins) varies as exp(n) (220 ). The situation is even more complicated for the SIA-clusters which are embedded into bcc matrix because the interaction between SIAs and the continuum of bcc states gives rise to many more configurations. As a result, the full investigation of

Construction of the stable C15 clusters

the energy landscape for large SIA-clusters (more than 10 interstitial atoms) is a herculean task. So, there is a need to develop a new strategy in order to search the global minimum. For the sake of simplicity, we use the subsequent diamond network formed by the center of Z16 polyhedra in order to generate the C15 clusters by pointing/marking only the center of the Z16 polyhedron. Similar to the configurations represented in the Fig. 3.5, we can construct larger clusters. In spite of this simplified representation, the number of possible configurations remains large. After a careful study of the various configurations and their corresponding formation energies, some selection rules are proposed in order to limit these choices:

Rule 1: the generated configuration will contain only those Z16 centers which are connected to its nearest neighbors. This rule prevents the construction of configurations formed by two (or more) disconnected clusters, e.g. we eliminate those n SIA cluster configurations which consist of two separate clusters of p and q self-interstitials atoms, where p + q = n. Let us take the case of I C15 4 which has two Z16 centers located in random lattice positions. According to this rule, only the configurations with centers in nearest neighbor position will be treated (as in Fig. 3.5b). All other possibilities will generate two disconnected clusters and are ineligible. This rule leads to only one possible configuration for cluster with three Z16 centers (as in Fig. 3.5c).

In order to prevent the construction of open chains of Z16 centers, we impose the second rule which allows only those clusters which have a special topology. Let's take the example of I C15 10 , presented in the Fig. 3.7, where many configurations for n=10. All the EAM potentials confirm the same trend: the lowest energy configuration forms a closed ring with 5 Z16 centers. This observation leads to the second rule.

Rule 2: Closed hexagonal path made of 6 Z16 centers are favored whenever possible. The smallest structure having 6 Z16 connected centers is the I C15 11 cluster which is shown in the Fig. 3.5e or in the Fig. 3.8, using the network of centers of Z16 representation. Loop closure then occurs for specific sizes, referred to as magic numbers. The next magic numbers are observed for 17 and 23 SIAs. These structures indeed have very low formation energies. E. g. the formation energy for 11 SIAs is particularly low due to the completion of one closed loop. Such low formation energies can be observed for some specific number of atoms and these are referred to as magic numbers. Next magic numbers are observed for I C15 cluster which is created by a planar geometry of 2 closed loops. Owing to different number of interstitial atoms in the two configurations, conclusive interpretations can't be made. Therefore, different configurations of 22 interstitial atoms are used to compare energies of planar and compact 3D forms. Careful observation reveals that the closed loops in compact form have lower formation energies as compared to closed loops in planar form. Hence: Rule 3: the C15 clusters must be constructed in the most compact 3D way respecting the first 2 rules.

Having developed new potentials in the previous section and having laid down selection rules for construction of stable C15 configurations, the formation energies of C15 clusters 1/2 111 dislocation loops will be compared using different empirical potentials. 993 Development of Embedded-Atom Method potentials for defects in iron

Relative stability of C15 clusters and dislocation loops

As repeatedly emphasized, precise knowledge of the energy landscape of interstitial defects is central to carrying out reliable simulations. Here we compare the formation energies of the dislocation loops 1/2 111 and C15 clusters. Using the pragmatic approach presented in the previous section expressed by the three "selection rules" we constructed configurations of C15 clusters up to hundreds of SIAs. Fig. 3.9 shows the formation energies of these two defect-types with cluster size for different empirical potentials. All the presented potentials show a crossover between 1/2 111 and C15 clusters at sizes lower than 100 SIAs. However, as evident from Fig. 3.9, these crossover values are dependent on the empirical potential used. It is thus impossible to give a precise value for this crossover from the presented data. As mentioned in the introduction: these empirical approaches give a good basis for description of the physics but the quantitative prediction are dependent on the potential and thus, peculiar. This characteristic remains even if great effort are made in order to fit a reliable set of parameters, as in the present development. These inconsistent crossover values, when used in simulations, can result in conflicting predictions.

Having reviewed the various atomistic methods, we conclude that this conflict can be resolved only by going beyond the realm of atomistic methods. In the following chapter, continuum methods are combined with accurate atomistic methods to assess their ability to address the shortcomings of the empirical potentials.

Conclusions

In this chapter we have reviewed the various empirical potentials used in the study of defects, such as dislocation loops or C15 clusters, in iron. Unachievable computational requirements of DFT calculations have fueled the search for alternatives considering reasonable approximations over the past decades. This search has led to the development of a number of empirical potentials, ranging from pair potentials to EAM potentials. Although these empirical potentials have been successful in making radiation damage 100 3 Development of Embedded-Atom Method potentials for defects in iron studies feasible, inconsistency of results from different empirical potentials is a major shortcoming that hinders conclusive theoretical results.

Conclusions

In this chapter, new empirical potentials are developed for Fe, mainly to study irradiationinduced defects such as self-interstitial atom clusters or dislocation loops. These empirical potentials are developed using embedded atom method formalism and are fitted on experimental values of solid-state properties, ab initio formation energies of basic point defects and ab initio forces acting on the atoms in the liquid or random state configurations. Various bulk and defect properties are compared to validate the transferability of the new potential.

We have also explored the energy landscape of the three-dimensional C15 self-interstitial atom clusters using the present potentials. The complex energy landscape of C15 clusters is explored and we establish three selection rules which facilitate the construction of lowest energy cluster configurations. These rules have practical importance enabling the construction of C15 clusters with hundreds of interstitials with minimal effort.

The newly-developed potentials are compared to selected existing potentials and are used to address the relative stability of 1/2 111 dislocation loops and C15 clusters. This analysis gives a crossover between C15 and 111 cluster below 100 SIAs. However, a precise value is missing despite the effort put in the fitting of the new empirical potentials. The shortcomings of the existing Fe potentials highlight the lack of any universally applicable potential.

Ironically, even developments and improvements in the empirical potentials fail to reach a consensus on the energy landscape of radiation-induced defects in bcc metals. This provides motivation to continue development of novel methods in order to eliminate this ambiguity based on choice of potential.

The Discrete-Continuum model

The energetics of interstitial clusters with nanometer size plays an important role, being a key ingredient that enables the connection between the asymptotic limits: isolated point defects that can be modeled using ab initio methods, and large observable dislocation loops. Limitation on the size of a DFT simulation cell in transition metals does not permit the exploration of clusters containing more than a few tens of SIAs. This technical problem can be overcome in part through the development of inter-atomic potentials based on the embedded atom method (EAM) as discussed in Section 2.3 and Chapter 3, but the unavoidable approximations result in the loss of accuracy and transferability. As explained in Section 2.3 and Chapter 3, most of the EAM potentials developed to study clusters of SIAs are built to fit the energetics of small clusters of SIAs provided by ab initio methods (33 , 37 , 98 , 99 , 101 , 134 ). Because of this, all the potentials provide similar results in the small cluster size limit but there is significant scatter in the predicted formation energies over the nanometer size range for loops (100 , 221 ) and C15 clusters (37 , 112 ). For instance, the EAM potentials proposed in Refs. (98 -101 , 213 ) can be used to compute the formation energy of nanometer-sized clusters in Fe containing 1000 SIAs in the form of 1/2 111 dislocation loops but the results span a fairly broad interval from 400 eV to 700 eV. Similar scatter is observed for other bcc elements, such as W, and for different types of clusters ( 100 or C15).

One way of circumventing this discrepancy is to establish scaling laws from elasticity theory and then to use these laws to extrapolate DFT calculations from small clusters to larger scales. One model, proposed by Soneda et al. (76 ) two decades ago, postulated an ad hoc function for the formation energy in terms of the number n of SIAs forming the cluster. The formation energy takes the form E f (n) = P 0 + P 1 n 2/3 , where P 0 and P 1 are adjustable parameters. This popular model, sometimes referred to as capillary model, has been widely used in the literature and over a hundred studies have used this simple law in order to parameterize kinetic Monte Carlo or cluster dynamics simulations for the time evolution of a distribution of clusters (see for instance Refs. (107 , 146 , 222 , 103 4 The Discrete-Continuum model 223 ) and Refs. therein). However, as we shall see further in this section, this model yields large uncertainties at large sizes when its parameters are fitted to properties of small clusters.

Since the elastic theory is known to work well for large interstitial defects but not at smaller sizes, we proceed to develop an atomistic model incorporating a 'clusterexpansion' type discrete contribution that will enable us to calculate the formation energy of interstitial defects at 0 K as a function of no. of SIAs. This involves integration of ab initio formation energies for small-sized defects with anisotropic elastic theory formulation for formation energy of large-sized defects and re-formulation of the general formula (without size-limitation) to calculate formation energy using the full anisotropic elastic theory and discrete contribution.

Objective: In this chapter, we develop a discrete-continuum model for the formation energy of clusters, which combines cluster expansion and elasticity for crystalline solids, enabling us to predict the formation energies for large SIA clusters directly from ab initio calculations performed on small clusters. The various sections, subsections, and their contents are as follows:

• Section 4.1: A systematic description of the key aspects of the continuum methods adapted to dislocation loops is presented. In section 4.1.1, the continuum part of formation energy formulation is developed starting from the anisotropic elastic theory for 2D loops.

• Section 4.2: The development of the discrete-continuum formulations for 2D loops and 3D C15 clusters are presented in the following subsections:

-Section 4.2.1 Formation energy formulation using discrete-continuum model for 2D loops.

-Section 4.2.2: Formation energy formulation for 3D clusters is developed by considering the 3D clusters as Eshelby's inclusions in an anisotropic cubic matrix.

• Section 4.3: Construction of database of small-sized interstitial (database) defects up to 22 SIAs is explained in this section. In particular, the motivation and simplified methodology of construction of both 2D loops and 3D clusters are discussed. ter, and the stacking fault energy. The latter is very high in bcc metals, and as a result stacking faults do not form, and are hence neglected in the following.

Other SIA clusters in bcc metals, such as C15 clusters, are 3D objects and have different elastic constants than the anisotropic bcc medium. These C15 clusters will be considered as Eshelby inclusions in anisotropic bcc matrix and the elastic energy of these C15 clusters will also be treated accordingly, as given in Sec. A1.1.

Under the influence of an external force, the material can permanently deform (plastic deformation) or regain its original shape (elastic deformation). (Note: Here, we restrict to elastic deformations.) Determination of elastic stresses and strains as well as their interdependence are key inputs to the study of their influence on materials. This shall be the focus in the following pages of this section. Based on the elastic theory, elastic formation energy of loops can be calculated as a function of the elastic constants and the distortion field components. This calculation has been carried out by two approaches, as illustrated by Bacon(81 ) and Stroh(224 ). Bacon's approach will be briefly summarized in the following subsection.

Formation energy formulation using anisotropic elastic theory for 2D loops

Here, derivation of the total linear strain energy per unit length without core contribution is presented in accordance with the elastic theory, as explained in ( 81),( 224). If the displacement field u i ( x) is known, then the elastic energy density can be expressed as a function of elastic stiffness tensor C ijkl and distortion components u i,j ( x) ≡ ∂u i ∂x j :

∆E = 1 2 C ijkl u i,j ( x)u k,l ( x), (4.1)
where i,j,k and l can vary from 1 to 3 for a three-dimensional body and summation over repeated indices is assumed. Due to high distortions at and near the dislocation, this definition only holds till a certain distance from the dislocation. This distance is approximated to be given by magnitude of Burgers vector and is called the core radius. The total elastic energy of an arbitrary dislocation is obtained by integrating the elastic density in Eq. 4.1 over its region of validity i.e. the entire volume V of the body except a cylindrical region at the dislocation with core radius. Total linear elastic strain energy is:

E = 1 2 V C ijkl u i,j ( x)u k,l ( x)dV = 1 2 S C ijkl u k,l ( x)u i ( x)dS j . (4.2)
In the previous equation, the volume integral is converted into closed surface integral using Gauss theorem. The surface which enclose the volume V should be carefully chosen. In the immediate neighbourhood of the core, the strain is too high to be described by elastic theory. Consequently, elastic field diverges and hence, the dislocation core region should be excluded from the elastic description. Excluding the core region, the closed surface integral in Eq. 4.2 has three different contributions: S δ , S R and S + , S -. S δ is the curved surface of the cylindrical core region surrounding the dislocation at a distance equal to core radius, S R is the external surface of the body and S + , S -are the straight surfaces connecting the two curved surfaces (S δ and S R ), as it presented in the Fig. 4.1. Using this contour surface the singularity due to the dislocation is removed. Thus, the total elastic strain energy can be split into various surface contributions as:

E = E S + + E S δ + E S -+ E S R . (4.3)
In order to simplify the previous equation we combine E S + and E S -into E S , i.e. S + and S -are collapsed into S (S represents the limit when two surfaces approach infinity close each to other). From the definition of Burgers vector b of the dislocation line, we can set a leap in the solution of the elastic field u i (S -)u i (S + ) = b i . We will denote the outward unit normal vector from V as n S on S -andn S on S + . We obtain the following equation which is valid for an infinite dislocation ( 81):

E = 1 2 S C ijkl b i n S j ( x)u k,l ( x)dS + 1 2 S δ C ijkl u i ( x)u k,l ( x)dS j - 1 2 S R C ijkl u i ( x)u k,l ( x)dS j . (4.4)
Static equilibrium of the body imposes the condition that the total force on the body is zero. Hence, using body force free displacement equations for static equilibrium condition,

C ijkl u k,lj ( x) + f i = 0, (4.5)
where f i denotes body force density. It can be shown that the integrals over S δ and S R cancel each other because the resultant of forces in the core region or the entire body is null. So, the final equation becomes:

E = E S . (4.6)
Using isotropic elasticity the Eq. 4.6 can be easy integrated into analytic form. E.g. following, the textbook deduction of Hirth and Lothe (103 ) the energy per unit line of 4.1 Continuum Methods for 2D dislocation loops edge dislocation can be written as (Eq. 3.55 of ( 103)):

E = µb 2 4π(1 -ν) ln R δ (4.7)
where µ is the shear modulus of the isotropic material, ν is its Poisson ratio. R is the effective range of the elastic field of the loop and δ is the radius of the dislocation core.

In anisotropic elasticity is impossible to integrate Eq. 4.6 in analytic form. Historically, two formulations of anisotropic elasticity are used in order to evaluate this equation: Bacon's theory ( 81) or Stroh's sextic formalism (224 , 225 ). Wide numbers of studies use these formalisms. In the study of defects, we cite only a few, such as the work of Dudarev (27 , 133 ) using Bacon formalism or the work of based on Stroh formalism. In this study, both formulations have been tested and yielded very similar results. Here, we briefly describe the Bacon's formalism.

Bacon's formalism utilized some results provided by Willis in 1970. Willis (231 ) derived displacement and distortion fields due to a straight dislocation at a point x with R | x| as follows:

u i ( x) = 3 η=1 F i ( χ η )ln χ η . x R , (4.8)
u i,p ( x) = 3 η=1 χ η p F i ( χ η )ln 1 χ η . x , (4.9)
where the function F i ( χ η ) is defined as:

F i ( χ η ) = 1 π b j C lmjk n k χ η m N il ( χ η ) n r ∂D ∂χr ( χ η ) , (4.10)
where denotes the imaginary part. Given that l is a unit vector along the dislocation line, m and n are unit vectors such that m, n and l form a right-handed set, or m = n× l.

For an edge dislocation, n = b/b and χ η = m + nω η . ω (1) , ω (2) and ω (3) are the three complex roots of the sextic equation S(ω) = det [C ijkl (m j + n j ω) × (m l + n l ω)] = 0, which are situated in the upper half of the complex plane ω = ω + i ω, where ω > 0, and N ik ( χ) is the matrix adjoint to L ik ( χ) = C ijkl χ j χ l , and D( χ) = det L ik ( χ).

Integrating Eq. 4.4 between the limits δ and R with the field expressed by Eqs.4.8 and 4.9, the total linear elastic strain energy per unit length, without the core contribution 4The Discrete-Continuum model can be written as: 4.11) where R is gain the effective range of the elastic field of the loop and δ is the radius of the dislocation core.

E = 1 4π b i b m n j n q C ijkl C npmq 3 η=1 χ η p χ η l N kn ( χ η ) n s ∂D ∂χs ( χ η ) ln R δ , (
Comparing the edge dislocation energy per unit line using anisotropic elastic theory, given by the the above equation, and the isotropic form given by Eq.4.7 it can be noted that the results have similar structure. Both approximations, isotropic and anisotropic, can be written into compact form:

E = Kb 2 4π ln R δ (4.12)
K is known as energy coefficient or prefactor or prelogarithmic factor of dislocation line and for edge dislocation which has the values of µ/(1ν) in isotropic elasticity or complicated form, which can very easily be deduced from Eq. 4.11, in the anisotropic elasticity case.

Although the major contribution to elastic energy comes from the total linear elastic strain energy without core contribution, the core contribution is required to accurately determine the elastic energy and to further use it as a state function. Thus, the complete expression to calculate total elastic energy can be expressed as ( 27): (4.13) where E is the total linear elastic strain energy per unit length without core contribution as described by Eq. 4.11, E δ is the core-traction energy per unit length and E c is the nonlinear core-traction energy per unit length.

E total = Ed l + E δ d l + E c d l,
To derive the formation energy of dislocation loops directly from anisotropic elasticity is very difficult. However, we can obtain an analytic expression using isotropic elasticity. Moreover, using the above mentioned similarity between isotropic and anisotropic formulation, we can provide an analytic formula in both approximations. Considering circular prismatic dislocation loop in an isotropic medium with its Burgers vector perpendicular to the loop's plane, the total energy can be derived to be (as in Section 6.4 of Hirth and Lothe (103 )): where µ is the shear modulus of the isotropic material, ν is its Poisson ratio and R is the radius of the equivalent circular loop. As mentioned, this equation is only applicable to isotropic materials and is inadequate when dealing with anisotropic materials. Replacing µ/(1ν) by K, as energy prefactor of dislocation line, we obtain the common formula for formation energy of prismatic dislocation loops using isotropic and anisotropic elastic theory:

E total = 2πR µb 2 4π(1 -ν) ln 4R eδ + E δ + E c , ( 4
E total = 2πR Kb 2 4π ln 4R eδ + E δ + E c , (4.15)
Using the above formulation of the formation energy, Dudarev et al. ( 27) brought forth two main ideas to predict relative stability of different morphology of loops with change in temperature for different bcc metals. Firstly, the full anisotropic elasticity approximation should be considered to treat dislocation loops in bcc metals due to strong anisotropy of some bcc metals, like iron. As a consequence, ratio of prelogarithmic term for different Burgers vectors in Eq. 4.11 was shown to vary depending on the anisotropy of the material considered(seen in Fig. 4.2). Secondly, an approximation to study temperature dependence of the linear elastic energy involves consideration of temperature dependence of elastic constants in the prelogarithmic energy factor in Eqs.4.11. Using this approximation, they demonstrated the stability of 100 loops in iron at high temperatures which was not explained till then.

A comparison of the 1st terms in Eq. 4.13 and Eq. 4.15 relates the parameter K to the integration of the total linear strain energy per length over the perimeter of the circular loop. Except for special cases, K can not be determined analytically and is generally calculated numerically. Since Kb 2 value in Eq. 4.15 represents the dominant contribution to total elastic energy calculation, the ratio of k b = Kb 2 of 100 and 1/2 111 loops gives a reliable way to estimate the relative stability of loops. Clearly, when this ratio is greater than 1, 111 dislocation loops are more stable than those of 100 . Observing Fig. 4.3, it can be concluded that the 1/2 111 dislocation loops remain more stable than 100 for all bcc metals of this thesis, except iron. There appears to be a crossover of this ratio from greater than 1 to less than 1 indicating enhanced stability of 100 dislocation loops at higher temperatures.

Although E. Clouet (232 ) provides an expression to calculate the core traction contribution in an anisotropic material, isolated dislocation and isolated dislocation dipole were treated. As such, core traction contribution for dislocation loops cannot be calculated using this expression. Since the non-linear remaining terms, E δ and E c , cannot be determined experimentally or analytically. It will be shown in the next chapter that these terms can be fit on feasible DFT calculations to obtain reasonable estimates on extrapolation.

The formula for formation energy of prismatic dislocation loop, given by 4.15 can be related to the number of interstitials which form that loop. Imposing that the surface area of the loop is equal to n times the average surface area per SIA, we obtain that R * is f b a bcc √ n, where a bcc is the lattice parameter of the bcc structure and f b is a factor that depends on the Burgers vector and the habit plane of the loop, i.e. f b = 1/(2π) and √ 2/(2π) for the 100 {001} and 1/2 111 {110} loops, respectively. The latter consideration allows us to rewrite Eq. 4.15 as a function of three unknown parameters, T, P 0 and P 1 : On the basis of this graph, the 1/2 111 dislocation loops are more stable than 100 for bcc metals, except for iron at higher temperatures (greater than 830 K). The k b -values are plotted up to melting point for the bcc metals, except iron where the plot is limited till its Curie temperature.

E elastic (n, b; δ, E c-δ ) = T √ n ln (n) + P 1 √ n + P 0 , ( 4
where the P 0 term is introduced to match the atomic data in the limit of small SIAs cluster size, e.g., n = 1 or n = 2, for which the concepts of perimeter and surface are not well defined.Using the sextic formalism (27 , 81 ), for the case of pure prismatic loop, the term T can be written as: .17) with all the definitions givrn above, in Eq 4.11.

T = 1 2 f b a bcc ln f b a bcc δ b i b m n j n q C ijkl C npmq 3 η=1 χ η p χ η l N kn (χ η ) n s ∂D(χ η ) ∂χs dθ. ( 4
The elastic theory, detailed above, is adapted to the treatment of large clusters. However, we note from Eq. 4.16 that the elastic energy varies as the square root of n, so that for small clusters, different contributions, either from the shape of the loops or from the internal structure of the loop, are expected to become dominant below a certain value of n. Additionally, the values of δ and E δ-c in Eq. 4.15 cannot be determined solely from elastic theory but they must be determined from atomistic calculations. This implies that parameters P 1 and P 0 can not be calculated analytically and are unknown.

An example of parameterization of Eq. 4.16 is shown in Fig. 4.4a. The best set of parameters for this model has been obtained using a database of clusters smaller than 22 SIAs, which are accessible to a DFT computation. The atomistic formation energies are computed using an EAM inter-atomic potential (99 ) in order to check the validity of our parametrization for large SIAs clusters through a comparison between predictions and direct atomic scale simulations. Two strategies have been tested for fitting. In the first case, all the three parameters, P 0 , P 1 and T were fitted. For that case, not described here, the predictions made from Eq. 4.16 for large clusters stringently diverge from the atomistic values, with some relative errors up to 60% for the two families of loops that have been examined. In a second method, T was computed from the elastic tensor as shown in Eq. 4.17 and only P 0 and P 1 were adjusted with respect to the formation energies computed at the atomic scale for clusters with n < 22. Using this approach, as can be seen in Fig. 4.4a, the predictions are much better since the error is around 10% for the 1/2 111 and less than 4% for the 100 loops. The error is smaller than ad-hoc laws proposed by Soneda (76 ) but still important in absolute terms for 1/2 111 loops. Even worse, the accuracy of the elastic model depends strongly on the choice of the database used for the fit. For the same number of clusters involved in the database, we can arbitrarily change the error by choosing various shapes of clusters. In order to reduce the variability of the results due to differences in the shapes of small clusters included in the database, and to reduce the relative error below 3% for the two types of 4.2 Formation energy formulations of defects using discrete-continuum model loops, we must increase the maximum size of loops included in the database to 53 SIAs. However, this is not accessible to DFT simulations in transition metals with conventional computers because to obtain the formation energies for clusters larger than 53 SIAs with reasonable accuracy, the total number of atoms needed in the simulation cell is larger than 5000. The fact that a pure elastic model parametrized on small clusters fails to correctly predict the formation energies for large clusters can be ascribed to the two following reasons: (i) the description of dislocation loops with finite core extensions is inappropriate for small clusters where the enclosing dislocation core is comparable in size to the radius of the cluster; (ii) the perimeter of the enclosing loops is fixed by a function indexed on integer values i.e. the number of interstitial atoms. In order to emphasize the latter point, we have reported in Fig. 4.4b the ratio of the convex hull perimeter of loops to the perimeter deduced from n using the criterion described above. This ratio converges very slowly to 1 and for sizes included in the fit (n between 3 to 21 SIAs) the value ranges from 0.56 to 0.85. Even if the ratio of convex hull perimeter and perimeter deduced as function of n using √ n criterion can be improved, the ambiguity in the definition of perimeter of small and large loops remains. Therefore the two points noted above imply that the parameters fitted to the data derived for small loop sizes are not representative of larger loop sizes generating large errors in the adjustment/extrapolation procedure.

Formation energy formulations of defects using discrete-continuum model

As the name suggests, the discrete-continuum model combines discrete contribution consisting of DFT calculations with anisotropic elastic theory to obtain the formation energy formulation for 2D loops and with Eshelby's inclusion method to obtain formation energy formulation for 3D C15 clusters. The development of these two discrete-continuum formulations are presented in the following subsections.

4

The Discrete-Continuum model

Formation energy formulation of 2D loops using discrete-continuum model

Hence, due to the size limitation of ab initio calculations, it is impossible at present -or in the near future -to parameterize an elastic model using Eq. 4.16 for the formation energies of SIA nanometric clusters. In order to overcome this difficulty, we add a cluster-expansion like term to the elastic model, which takes into account the discrete structure of small dislocation loops:

E f ormation (n) = E discrete (n, n 1 , n 2 ..; {P j }) + E elastic (n, b; δ, E c-δ ) . (4.18)
The discrete term depends on a set of parameters {P j }, and we impose a requirement that it vanishes in the asymptotic limit n -→ ∞ i.e.E discrete (n,

n 1 , n 2 ..; {P } j ) -→ 0 as [n → ∞].
The discrete nature and the geometric structure of clusters are accounted for in the term E discrete through a topological mapping to the local neighborhood of each dumbbell which is defined by the number of first (n 1 ), second (n 2 ) or higher nearest neighbor pairs of dumbbells. The distance between dumbbells is defined as the distance between their centers. For example in the case of a 1/2 111 dislocation loop with a {110} habit plane, the first and second nearest neighbor shells each have 4 nearest neighboring dumbbells, situated at distances √ 3a bcc /2 and a bcc , respectively.

The discrete part of the energy for a dislocation loop containing n SIAs is written as the sum of contributions from all dumbbells:

E discrete = n i=1 E i = n i=1 f (n; n i 1 , n i 2 )E(n i 1 , n i 2 ). (4.19)
The local energy associated with the i th dumbbell of the cluster is expressed as

E i = f (n; n i 1 , n i 2 )E(n i 1 , n i 2 )
, where the function E(n i 1 , n i 2 ) fully determines how E i depends on the dumbbell neighborhood, i.e. on the number of the 1 st and 2 nd nearest neighbor dumbbells in the habit plane, denoted n i 1 and n i 2 , respectively. Function f (n; n i 1 , n i 2 ) fixes the weight for the i th dumbbell energy E(n i 1 , n i 2 ). In order to define the latter function, we note that various atomic scale studies (100 , 233 , 234 ) have confirmed that the interatomic distance between two atoms that form the dumbbells situated far from cluster's edges recovers bulk coordination. For example, relaxation of dumbbells recovers perfect bulk 1 st nearest neighbor distance √ 3a bcc /2 in the center of clusters. As a result, the dumbbells that are close to the center of loops, with full nearest neighbor shells, make no contribution to the energy of the system other than twice the cohesive bulk energy. This means that in terms of the formation energy, these dumbbells give no contribution to the discrete energy. Therefore the function f (n, n i 1 , n i 2 ) should be zero for the dumbbells with their full nearest neighbour shell. A second constraint on this function is given by the asymptotic limit at large n, i.e. f (n; n 1 , n 2 ) -→ 0 as [n -→ ∞. Hence we consider the following product form: f (n; n 1 , n 2 ) = g(n)h(n 1 , n 2 ), where h(n 1 , n 2 ) equals unity for the atoms which do not have full nearest neighborhood, and zero otherwise, and g(n) -→ 0 for large clusters. In order to ensure the condition E discrete (n, n 1 , n 2 ..; {P } j ) -→ 0 as n -→ ∞ we have tried many monotonically decreasing functions of the form proportional to 1/n α for g(n), with α in the interval from 0.5 to 1. The best choice for α was found to be 0.55.

To reduce the sum in Eq. 4.19, we can rewrite the discrete energy contribution by introducing the number of dumbbells having n 1 and n 2 first and second neighbors, (i n 1 ,n 2 ): (4.20) where N b 1 and N b 2 are the bulk numbers of first and second neighbors and P 2 is a constant. A pair formulation is neither a necessity nor a constraint in this approach, since the model 4.2 Formation energy formulations of defects using discrete-continuum model can be readily extended to more complex types of interaction. In order to exemplify this energetic model, let's take the example of a small 100 {001} loop containing 15 SIAs, which is sketched in Fig. 4.5. In this case, the discrete part of the energy can be written as: 55 [2E(1,2) + 4E(2,2) + 3E(2,3) + 2E(4,2) + 3E(4,3)] + P 2 . (4.21) In the above equation h (4, 4) = 0 andh(1, 2) = h(2, 2) = h(2, 3) = h(4, 2) = 1. Combining the elastic and discrete parts, we find that the formation energy of a loop with n SIAs is:

E discrete = N b 1 n 1 =0 N b 2 n 2 =0 i n 1 ,n 2 f (n; n 1 , n 2 )E(n 1 , n 2 ) + P 2 ,
E discrete = 1 15 0.
E f ormation (n) = T √ n ln (n) + P 1 √ n + P 0 + N b 1 n 1 =0 N b 2 n 2 =0 i n 1 ,n 2 f (n; n 1 , n 2 )E(n 1 , n 2 ). (4.22)
We note that parameter P 0 derived from elasticity is combined with P 2 , deduced from the discrete model, to give just one constant, denoted by P 0 in Eq. 4.22. Our goal here is to produce an analytical model that defines a general functional form of the scaling law that describes the formation energy E f ormation (n) of clusters as a function of their size n. The advantage of this formulation, Eq. 4.22 is that a full set of parameters E(n 1 , n 2 ) and P 0,1 can be obtained from ab initio formation energies derived using a training series of configurations of small interstitial clusters. This new model combines a discrete contribution to the energy, evaluated using a cluster expansion formalism, with a term derived using a treatment of prismatic loops based on elasticity theory as described in the chapter.

Formation energy formulation for 3D clusters using discrete-continuum model

The strategy described above has also been adopted to develop a model for C15 clusters. C15 inclusions have different elastic properties in comparison to the host bcc matrix (41 ), and the corresponding energy is treated using the formalism of isotropic Eshelby inclusion (225 ,235 ,236 ). The discrete contribution to the formation energy takes into account the particular structure of C15 clusters. In comparison with the case of dislocation loops, they have an additional contribution from the atoms having perfect C15 bulk coordination situated inside the clusters. Thence the core region of a C15 cluster does contribute non-negligibly to the formation energy. In addition, the bcc 4The Discrete-Continuum model bulk atoms of the perfect lattice are replaced by C15 bulk atoms with different cohesive energies and consequently this difference must be accounted for. The remaining atoms of the SIA cluster, which do not have the perfect C15 bulk coordination, correspond to the interface between the C15 cluster and the bcc matrix. These interfacial atoms also contribute significantly to the formation energy. The present energetic model is close to the Zhang et al. (112 ) model used for interpolating the formation energy of C15 clusters provided by EAM inter-atomic potential calculations. In the Zhang et al. model, the number of atoms situated both at the interface and in the perfect C15 bulk, are deduced from the asymptotic limit of large clusters. The present model is used for predicting the formation energy of C15 clusters from DFT calculations. The difference between the convex hull surface and the surface computed from the number of SIAs in a cluster is fairly large for small sizes (even larger than that for perimeter of loops, Fig. 4.4b). Therefore, the number of atoms with perfect C15 bulk coordination N C15 , the number of interfacial atoms N i , as well as the surface area S C15 and volume V C15 of C15 clusters are deduced directly from the geometry of the cluster. The formation energy expression used in our model is written as follows:

E f ormation (n) = S C15 γ + 6V C15 µε 2 α + N C15 (E bcc coh -E C15 coh ) + N i ∆E i , (4.23)
where E bcc coh and E C15 coh are the cohesive energies of the perfect bulk bcc and of the perfect C15 structures, respectively. Coefficient ∆E i is the average energy of atoms at the interface and γ is the interface energy per unit area between the bcc matrix and the C15 inclusion. The second term in Eq.4.23 is the energy of isotropic Eshelby's inclusion, with a C15 cluster treated as an inclusion in the otherwise isotropic bcc matrix. Eshelby discovered an elegant way of calculating the stress, strain and displacement fields, both in the inclusion as well as in the matrix, by using a superposition of linear elasticity and Green's function formalism (235 , 236 ). The same approach allows the computation of the strain energy contribution in the presence of a C15 cluster. Here, µ is the isotropic shear modulus of the matrix, while α = 1 + 3µ/(4B C15 ) where B C15 is the bulk modulus of C15 clusters and ε is the misfit strain, which can be computed directly from atomic scale calculations. ∆E i and γ are adjusted with respect to atomic scale simulations while the cohesive energies of bcc and C15 clusters are determined from atomistic calculations (41 ). In the limit of large spherical C15 clusters, the previous equation can be written as a function of the number n of SIAs as in the Zhang model ( 112):

E f ormation (n) = 2γ s 9πΩ 2 1/3 n 2/3 + 12Ωµε 2 α n + 3n(E bcc coh -E C15 coh ), (4.24)
where Ω is the atomic volume of bcc iron. In the large limit, a cluster with n interstitials is obtained by replacing 2n bcc atoms by 3n C15 atoms which gives the volume of the database, are given in the following section.

Parametrization of the discrete part of the present energetic model relies on the topology of cluster configurations. The number of first or second nearest neighbors of each dumbbell is an essential ingredient for the formation energy calculations in Eq. 4.22. In this section, we present the constraints that we have imposed in the construction of cluster geometries in order to set-up correctly, from the mathematical point of view, the fitting procedure of E(n 1 , n 2 ) parameters.

Database construction of dislocation loops

The construction of extrapolation database is based on previous studies of various authors (27 , 100 , 110 , 133 , 240 ) which assert that the closed-loop configurations such as rectangles, squares or circles are more stable than possible elongated configurations for the same number of SIAs. Adhering to this requirement, hundreds of configurations of 100 and 1/2 111 dislocation loop types were constructed for SIAs ranging from 2 to 1500. For all these configurations, we investigated the local environment of each dumbbells. To facilitate better understanding of the neighborhood behavior in this set of configurations, occurrence of each (n 1 , n 2 ) pair was plotted where n 1 refers to the number of first nearest neighbors and n 2 refers to the number of second nearest neighbors. The plotted histogram revealed that certain (n 1 , n 2 ) pairs do not occur while other pairs are overrepresented as shown in Figure 4.6.

The configurations contained in the training database of the discrete-continuum model were built in keeping with the selective dumbbell neighborhood behavior of large clusters expressed above. These configurations were constructed such that all the occurring (n 1 , n 2 ) pairs appear in the database as well so that each occurring pair is considered for fitting. Due to the compact form of clusters, there are some extra (n 1 , n 2 ) pairs of neighbors that occur in database. However, this is not expected to pose a problem because they do not contribute significantly in the extrapolation. The histogram of the neighborhood of the training database is shown in All energies are normalized to the asymptotic limit, taken as the formation energy in a simulation cell containing 207646 + n atoms. Lines are guides for the eyes, obtained using a fit to a 4 th order polynomial. Formation energies have been corrected using the elastic dipole correction method to account for the finite size of the simulation cells (230 ).

Construction of database

Validation of the Discrete-Continuum model

Database construction of C15 clusters

The present database for C15 clusters should include the lowest energy configurations.

In order to form a cluster with a given number of SIAs, the number of possible choices of Z16 centers is quite large and guessing the configurations with the lowest formation energies is problematic. The database for C15 clusters contains up to 20 configurations of SIAs. This limit is fixed, as in the case of loops, by the accuracy in the formation energies derived from DFT calculations with a cell containing 1024+n atoms.The number of possible configurations grow exponentially with the size of the cluster, making a systematic search prohibitive at larger sizes. We use the pragmatic approach, configurations for the present database were generated using the three selection rules presented in detail in the Sec. 3.4.

Validation of the Discrete-Continuum model

Discrete-continuum model can be parameterized through simulations performed using different EAM inter-atomic potentials, which allows us to test model predictions for large clusters, using large simulations cells. Several EAM potentials for Fe (99 -101 , 213 ) and for W (33 , 85 , 134 , 135 ) were used for our tests. Note that in the Fe EAM potential published in Ref.(100 ), a typo was corrected in Ref. (37 ).

The set of cluster geometries used for training the discrete-continuum model with EAM energies is the same as the one that will be used later for parameterizing the model from DFT data. The tests were performed for 100 {100} and 1/2 111 {110} dislocation loops containing up to 1200 SIAs. Three types of shapes were considered to construct configurations: rectangular, circular, and hexagonal, where the sides of the polygon correspond to the dense directions of the habit planes. The database of C15 clusters contains configurations with sizes up to 110 SIAs. The C15 configurations were mostly generated in accordance with the three rules mentioned in the previous section. The few configurations which do not obey these rules will be discussed later.

The atomistic formation energies of clusters of SIAs were computed using zero Kelvin atomic relaxation simulations. The asymptotic values of the formation energies were obtained by introducing interstitial clusters in a constant volume simulation cell with millions of atoms, sufficient to remove any residual size effect. The system was relaxed using a conjugate gradient technique with a convergence criterion on the maximal force per atom of lower than 0.02 eV/ Å. We have also performed calculations where the criterion was 0.001 eV/ Å, resulting in minor changes in the formation energies, less than 0.001 eV.

The formation energies calculated with EAM potentials were compared with predictions made using the discrete-continuum model (see Fig.4.10). For dislocation loops the difference is less than 2% (see Fig. 4.10a and 4.10b for 1/2 111 and 100 loops, respectively) using the EAM potential from Ackland-Mendelev for Fe (99 ). Similar results are obtained for all the EAM potentials tested for Fe and W. For the C15 clusters the error is slightly larger, i.e. 3% (see Fig. 4.10c). The error may reach 5% for some clusters e.g. clusters containing 62, 64, 66 and 67 SIAs. The main reason for this discrepancy is that such clusters don't have compact geometries, hence they break the spherical symmetry assumed in the model. These configurations were created by infringing the third rule given in section 4.3. We estimate that such configurations are not significant for the purpose of this study, being far from the lowest energy configurations.

Calculation of formation energy for interstitial defects at 0K

Having validated the discrete-continuum model, we can now proceed to its parameterization using the DFT formation energies of configurations included in our database. For calculation of DFT formation energies of database configurations, the DFT simulation cell for n SIAs has been chosen to contain between 250+n and 1024+n atoms in such a way that the relative error in the formation energy is lower than 2% from the converged values (see Fig. 4.9). The DFT calculations were performed using VASP within the projector augmented wave (PAW) framework (241 -244 ). The plane wave energy cutoff is 350 eV and the Hermite-Gaussian broadening-width for Brillouin zone integration is 0.2 eV. The calculations are performed including the p semi-core states. The exchange correlation energy is evaluated using the Perdew-Burke-Ernzerhof (PBE) Generalized Gradient Approximation (GGA). The k-point grid mesh was chosen from 3 3 for 250+n cell up to (1 or 2 3 ) for the 1024+n cell. W and V are non-magnetic materials and iron is treated in the ferromagnetic state, which is a reasonable approximation in the low temperature limit Each configuration is relaxed using the conjugate gradient technique with a convergence criterion on the force on each atom of 0.02 eV/ Å. The size of the supercell remains fixed in order to ensure constant volume-per-atom simulations. All the formation energies were adjusted using the dipole correction (230 ). In Varvenne et al (230 ), it has been proved that dipole correction from constant volume, or zero strain, yield a good correction value for the formation energy of defects in ferromagnetic iron, such as dislocation loops or C15 clusters.

Ab initio based predictions of SIA cluster formation energies

The formation energies for the 100 and 1/2 111 loops, as well as C15 clusters, were computed using the discret continuum model for Fe, W, V (except 100 for V) and Ta, and the results are shown in Figs. 4.11, 4.12, 4.14 and 4.13, respectively. As mentioned earlier, for all bcc metals, the experimental evidence, within the limit of detection in TEM, for instance 1-2 nm radius for loops, confirms that the most frequently observed morphology at low temperature corresponds to 1/2 111 loops. The formation energies predicted by the discrete-continuum model for large dislocation loops 129 4 The Discrete-Continuum model (sizes larger than hundreds of SIAs) are in agreement with this observation since the 1/2 111 loops are predicted to be energetically more stable. The model also shows that, for dislocation loops from 10 SIAs to sizes visible in TEM, 1/2 111 loops always have smaller formation energy than the 100 loops in both Fe and W. The present study is at odds with some EAM potentials for Fe (37 , 100 ) and for W (33 , 133 , 135 ), which predict a crossover in the relative stability of two families of loops around 200 SIAs. Below this critical size the 100 loops would be more stable in W whereas they would be more stable above the critical size in Fe. The origin of this inversion in the relative stability of loops is still unclear. It is worth noting that the discrete-continuum model is able to reproduce the crossover predicted by the EAM potentials if the model was calibrated using the database corresponding to the same potential. When the model is calibrated to the database derived from DFT, the model predicts no crossover between the loop formation energies.

In contrast, the DFT-based predictions show crossovers between C15 clusters and loops. In Fe, one crossover appears with 1/2 111 loops at clusters around 51 SIAs in size, corresponding to a 1.5 nm diameter C15 cluster. There is also a crossover with 100 loops, this time both in W and Fe at 21 and 91 SIAs, respectively. In V there is no crossover, 1/2 111 loops are the most stable configurations for all defect cluster sizes.

In Fe, the present results reconcile the theoretical predictions with experiments, where only the 1/2 111 loops were observed under irradiation at low temperature, by giving some support to a mechanism recently identified as a possible route of formation of the 1/2 111 and 100 loops involving the collapse of larger C15 clusters (112 ). The possible formation mechanisms of 100 loops in Fe were addressed in the past in several studies, some examples are given in Refs. (79 , 105 , 106 , 245 ). In particular, Refs. (79 , 105 , 245 ) propose a mechanism based on the reaction between two 1/2 111 loops having appropriate size and specific orientations. The mechanisms proposed by Marian et al. (79 ) and Xu et al. (105 ) are similar, the only difference being that Xu et al. showed that this reaction holds for larger clusters and has stochastic components. Another scenario by Chen et al. proposed transformation of 1/2 111 loops into 100 by correlated translation-rotation of SIAs forming the loop (106 ). All these mechanisms involve a certain number of stringent conditions, such as the direction of loop migration and the size of the loops, which make the corresponding events highly infrequent. Zhang et al. (112 ) proposed an alternative idea involving the nucleation of C15 clusters and their growth by trapping of single self-interstitials of 1/2 111 dumbbell structure. In Fe, small C15 clusters are energetically very stable and act as traps for small mobile SIAs. Moreover, they are kinetically trapped, meaning that the lowest energy reaction pathway that allows C15 clusters to transform into planar loops corresponds to very large energy barriers resulting in highly improbable transitions. In Ref. (37 ), it was shown that the lowest energy pathway that transforms a 4-SIA C15 cluster into a planar loop is of the order of a few electron-Volts. Under irradiation, small mobile interstitial clusters, such as 1/2 111 or 110 loops, are continuously produced, facilitating the growth of C15 clusters which can reach very large sizes, even larger than the crossover between C15 and traditional loops because of their kinetic trapping.

At large sizes the transformation of C15 clusters into dislocation loops with 1/2 111 or 100 orientation becomes very likely. This transformation is demonstrated even on the time scale of molecular dynamics simulations, by Zhang et al. (112 ). Therefore, the 4

The Discrete-Continuum model frequency of formation of 100 loops is definitely larger than in any other mechanism proposed in the past. The only condition is that the C15 clusters should be formed at small sizes, which is confirmed by DFT calculation of ( 37) and present study for small (up to 8 SIAs) and large clusters (nanometric sized), respectively. Although the mechanism proposed by Zhang et al. is rather convincing to explain the formation of 100 loops at high temperature in Fe, it doesn't explain why these loops are not observed at low temperature (28 ). Present work resolves this contradiction by revealing the DFT relative energy of large clusters. All the past interpretations of Zhang et al mechanism were based on EAM potentials energetic landscape which is different from the presents DFT findings. As shown in Fig. 4.11, the crossover between the C15 clusters with 111 and 100 loops occur at 51 and 91 SIAs, respectively, in Fe. It means that the C15 clusters which could form under irradiation and have sizes larger than 51 and smaller than 91 SIAs can decay only into the 1/2 111 clusters. This could explain the absence of 100 loops because the C15 clusters are more stable in this size range (between 51 and 91 SIAs). We expect that C15 clusters should have sizes much larger than 91 SIAs in order to have non-zero probability to transform into 100 loops, which further increases the size range where 100 cannot appear. Even though our interpretation does not exclude the possibility of 100 loop formation directly under irradiation at low temperature, it drastically reduces such probability in agreement with experimental observations (28 , 70 ).

For small cluster sizes in W and V, the formation energies of C15 clusters are much higher than for 1/2 111 loops. In W for small sizes, between 7 and 21 SIAs, the C15 clusters have slightly lower formation energies than 100 loops as shown in Fig. 4.12a and 100 clusters become energetically more favorable than the C15 clusters containing more than 21 SIAs. In V, 100 loops have the formation energies that are between those of 1/2 111 loops and C15 clusters at all sizes. We used a relatively restricted set of calculations to parameterize an energetic model for 100 loops in V, and so this conclusion is given on the basis of calculations for intermediate cluster sizes performed for 2, 4, 10 and 20-SIA clusters.

The energy landcape for Ta seems to be a special intermediate case between Fe and W, as shown in Fig. 4.13. Like Fe, a crossover between 111 and C15 clusters does exist but at a much smaller cluster size ( 17 SIAs). Beyond this crossover, the energy landscape resembles that of W with 111 being the most stable and formation energies of C15 clusters higher loops. Although conclusive experimental evidence of the irradiation-induced defects in Ta is missing, required resolution might pose a major obstacle. Nevertheless, 134 4.5 Calculation of formation energy for interstitial defects at 0K

Ab initio scaling laws for formation energy

One of the goals of this paper is to provide simple analytical scaling law formula for the formation energy of self-interstitial clusters. The interest of such a formulation is the practical application in multi-scale techniques including kinetic Monte Carlo simulations and cluster dynamics or dislocation dynamics studies. Using the present analytical scaling law, we restrict the input required for parametrization of defect energetics to the number of interstitial atoms and their type. These new simple scaling law provide reliable formation energies over a very broad range of defect sizes for any subsequent multi-scale study.

Therefore, based on the dependence of the elastic contribution of the formation energy of loops, in Eq. 4.16,and C15 clusters,in Eq. 4.24, on the number of interstitials, we propose a simple analytical expression in order to fit the formation energies predicted by the discrete-continuum model, for the loops,

E f (n) = a 0 √ n ln (n) + a 1 √ n + a 2 , (4.25)
and of C15 clusters,

E f (n) = a 0 n 2/3 + a 1 n + a 2 . (4.26)
It is worth noting that for sizes larger than 15 SIAs these two laws are a very good fit to the formation energies, with an absolute error lower than 1 eV. Parameters of Eq.4.25 and 4.24 for Fe, W, V and Ta are given in table

Conclusions

In this chapter, the formation energies of SIA clusters are investigated for four bcc metals, namely Fe, W, V and Ta. The main result was the development and validation of a discrete-continuum model that makes it possible to perform ab initio-level accurate calculations for clusters without any size limitation. The model allows us to treat various cases of interstitial dislocation loops and C15 clusters from clusters containing a few SIAs to nanometer size.

From the interpretation of the present results it can be concluded that above ∼ 100 SIAs 1/2 111 loops are always the most stable family of SIA-clusters -in agreement with experimental observations of irradiation defects at low temperature in bcc metals. However, these results are at odds with calculations made using various EAM interatomic potentials, which yield spurious predictions concerning the relative stability of 100 and 1/2 111 loops (100 ). Future developments of such potentials should consider the information provided in the present paper, and include the appropriate additional fitting conditions on the potential parameters.

Our study shows that in Fe, C15 clusters are the most stable clusters of defects for sizes lower than 51 SIAs, which is a size range not accessible to direct TEM observations. Our model also supports the theory of formation for 100 loops proposed by Zhang et al. (112 ). In the present work, we do not include thermal effects and magnetic excitations. As a consequence, our results are comparable only to low temperature experiments. The results obtained, sheds some light on the absence of 100 loops in low temperature experiments, and reconciles the Zhang mechanism with the experimental evidence. However, in order to validate entirely our expectations, further analysis is required.

Finally, our work makes it possible to establish scaling laws for the formation energies of various types of clusters in various materials, which is significant for multi-scale simulations such as kinetic Monte-Carlo simulations (107 , 223 , 246 -248 ), cluster dynamics studies (249 , 250 ), or mean field approximations (251 ), where simple analytic laws are needed to model the energy of large clusters. However, to enable the use of scaling laws in multi-scale simulations, the effects of temperature must be accounted for. The present formulation of the discrete-continuum model can be extended to address the formation free energies e.g. by including the temperature dependence of elastic constants. These 5 Finite-temperature extension of discrete-continuum model

In the previous chapter, a discrete continuum model based on atomistic parametrization was developed to calculate nearly accurate formation energies of interstitial defects at 0 K without size-limitation. In order to apply the discrete-continuum model to realistic environments, there is a need for extension of this model in terms of temperature. In this chapter, we aim to extend the developed atomistic model beyond 0 K, permitting calculation of nearly accurate formation energies of interstitial defects without sizelimitation at any given finite temperature.

Objective: The objective of this chapter is to develop a finite temperature extension of the discrete continuum approach. This new extension will be compared to atomistic free energy calculations for validation. When possible, the present development is also directly compared to the corresponding experiment.

The various sections, subsections, and their contents are as follows:

• Section 5.1: In this section, the importance and relevance of free energy calculations is explained.

• Section 5.2: This section deals with the extension of the previously-developed Discrete-Continuum model to obtain free energy at finite temperatures. Since study of temperature dependence of elastic constants of the chosen metals is central to this approach, this section also includes relevant discussion about adiabatic and isothermal elastic constants in Section 5.2.1.

• Section 5.3: Firstly, a brief summary of the canonical ensemble is presented. The goal of this section is to perform free energy atomistic calculations of interstitial 5 Finite-temperature extension of discrete-continuum model defects at finite temperature in order to validate the finite extension of the discrete continuum approach. This section includes the following subsection:

-Section 5.3.1: Harmonic and quasi-harmonic approximations are considered to obtain the temperature dependence of the observables required by the finite temperature discrete continuum model, such as lattice parameters, elastic constants or cohesive energies.

Apart from the obvious energy landscape, there are many possible applications of the current developments of the discrete-continuum model. The following two sections will treat two applications:

• Section 5.4 investigates the influence of temperature on the relative stability of various cluster in iron. Also, the relevance of the mechanism of formation of 100 loops in bcc iron at higher temperature is discussed.

• Section 5.5: Here, free energy values predicted by the finite temperature discrete continuum model are used as input in cluster dynamics simulations in order to compare the simulation results with time evolution of loops density in post irradiation experiments.

Free Energy

The objective of this chapter is to present a preliminary analysis of thermal impact on point defects properties, which influence the atomistic free energy landscape at finite temperatures. Unarguably, finite temperature excitations play an important role in the physics of defects. Consider, for example, self-diffusion in a bcc material at a finite temperature T, the concentration of defects as well as the diffusion rates are controlled by the free energy landscape of vacancies. As evident, study of any temperature-dependent process would directly or indirectly involve free energy calculation. For such calculations, ideally, Gibbs free energy definition may be used due to its general formalism. However, when adequate information is available about the process/system, the definition of free energy can be tailored depending on the thermodynamic variables of the process/system. In this section, we assume a zero-pressure condition and consequently, the thermodynamic functions like Gibbs (G = E + P V -T S) and enthalpy (H = E + P V )

Free Energy

are replaced by Helmholtz free energy(F = E -T S) and energy(E) of the system, respectively. The Helmholtz free energy surface F (T, V, N ), a function of no. of moles of the substance N , volume V and temperature T , is a central quantity in thermodynamics of solids because it fully determines all other thermodynamic quantities such as heat capacities or elastic properties at finite temperatures. From the differential of the Helmholtz function (dF = -P dV -SdT ), a simple thermodynamics relation in the form of partial differentiations can provide access to the entropy surface:

S(V, T ) = - ∂F ∂T V , (5.1)
or the pressure surface P (V, T ):

P (V, T ) = - ∂F ∂V T . (5.2)
Further, the isothermal bulk modulus at a given temperature, defined as change of pressure per unit fractional change in volume at constant temperature T, can also be calculated in terms of the Helmholtz free energy:

B(T ) = -V ∂P ∂V T = V ∂ 2 F ∂V 2 T . (5.3)
As such, the Helmholtz free energy is a quantity of great importance in understanding chemical and biochemical processes. This claim will be further supported when treatment of the canonical ensemble will be briefly explained in section5.3.

If we consider the concentration of defects, it is found to be related to the free energy landscape through the formation free energy landscape. The isobaric equilibrium concentration of defects is given by:

C eq (T ) = exp - F f (T ) kT , (5.4)
where the formation free energy of the defect F f (T ) is computed in the same way as the formation energy i.e. by taking the difference of the corresponding function of a perfect crystal and a crystal with defect. For example, the formation free energy in the case of n intertitials (or vacancies) in a crystal with N atoms and volume V :

F f (T, V, n) = F (T, V, N ± n) - N ± n N F bulk (T, V, N ) , (5.5)
using +/-for interstitials/vacancies, respectively.

where the parameter T is evaluated as:

T = 1 2 2πf b a bcc K (5.8) K = 1 2π ln f b a bcc δ b i b m n j n q C ijkl C npmq 3 η=1 χ η p χ η l N kn ( χ η ) n s ∂D( χ η ) ∂ χs dθ, (5.9)
and where all the symbols in these equations hold the same meaning as defined in Sec. 4.1.1 of Ch.4. While this formulation provides the energy landscape at 0 K, it can be extended to obtain the corresponding energy landscape at finite temperatures by means of a simple manipulation which facilitates inclusion of temperature dependence in the current 0 K formulation. The scheme proposed here is similar to the seminal idea of Dudarev et al (27 ). In the present formulation, more consistent schemes are proposed. Also, the limits of the accuracy of the model are provided by comparative studies with atomistic free energy calculations.

At the end of the previous chapter, extrapolations of the formation energies at different sizes provided by the discrete-continuum model were fitted with scaling laws. These scaling laws were simple analytical functions in the frame of pure elastic theory and were dependent only on the size of the interstitial defect clusters of a particular defect-type (Section 4.5.2). In order to include temperature dependence, we assume that the elastic functional form doesn't change at finite temperatures and that only the parameters of the scaling laws will vary. For the case of loops, the proposed function for the formation free energy F (n, T ) is similar to that of formation energy given by pure elastic formula Eq. 4.16:

F (n, T ) = T (T ) √ n ln n + P 1 (T ) √ n + P 0 (5.10)

In the above equation, temperature dependence is integrated in the definition of parameter T and P 1 (T ) by considering temperature dependence of the elastic or atomistic properties of the system. At 0 K, the two parameters are well defined being equal to:

T(0) = 1 2 2πf b a bcc (0)K (5.11) P 1 (0) = 2πf b a bcc (0) E δ-c (0) + K ln 4f b a bcc (0) eδ (5.12)
At 0 K, the formation free energy becomes the formation energy. By simple fitting procedure on the formation energies predicted by the discrete continuum model the two parameters P 0,1 can be deduced. The T parameter is imposed to the value numerically 5 Finite-temperature extension of discrete-continuum model given by Eq. 5.9. Following this fit, from the value of P 1 (0) at 0 K, the values of core-traction energy E δ-c can be deduced:

E δ-c (0) = P 1 (0) 2πf b a bcc (0) -K(0) ln 4f b a bcc (0) eδ (5.13)
The value of the core-traction term is an important ingredient of the scaling law. Moreover, the value of E δ-c can be used to test the reliability of the empirical potential and to predict the correct values of formation energies in the mesoscopic scale. At small sizes of loops, the core energy contributes significantly to the formation energy. The discrete-continuum model parametrized by the present DFT calculations predicts the values given in Table 5.1. It is interesting to note that these values have previously been computed using empirical potential, only. In the case of Fe, Dudarev (27 ) et al. reported the core-traction energy for 1/2 111 {11 2} ( 100 {001}) dislocation using two empirical potentials and molecular dynamics calculations as 0.82 eV/ Å (0.97 eV/ Å) and 0.71 eV/ Å ( 0.74 eV/ Å) for EAM potentials of and Marinica (100 ), respectively. It is noteworthy that the core-traction energy of 1/2 111 {11 2} dislocation is not far from the DFT values while the dislocation 100 {001} is largely underestimated with a value roughly half of the DFT value. This difference is found to be even more significant in the case of W when core-traction energy for the same dislocations is computed using two EAM potentials. The values are 1.24 eV/ Å (0.76 eV/ Å) and 2.51 eV/ Å ( 2.13 eV/ Å ) for and EAM4 of Marinica et al. (135 ), respectively. None of these investigated potentials are close to the DFT values and also predict an unphysical behaviour with a lower core energy in the case of 100 dislocations compared to 111 . This contradiction in values of E δ-c for the two types of loops 100 and 1/2 111 predicts unreliable results related to the relative stability of 1/2 111 loops compared to 100 loops, as mentioned in the first and fourth chapters for Fe and W. Furthermore, as we will discuss later, this term of scaling law is crucial in the extrapolation of the discrete-continuum model to finite temperature.

Having the core-traction energy as key ingredient, the discrete model can be extended to finite temperature. Certain assumptions should be made. Firstly, the P 0 parameters are assumed to be independent of temperature. The temperature dependence of the other two parameters P 1 (T ) and T(T ) can be considered in different ways, thus resulting in the following three models:

• In the first model, temperature dependence is ensured only through the elastic

E δ-c eV/ Å Fe W Ta V 1/
T(T ) = 1 2 2πf b a bcc (0)K(T ) P 1 (T ) = 2πf b a bcc (0) E δ-c (0) + K(T ) ln 4f b a bcc (0) eδ (5.14)
• In the second model, temperature dependence of the linear thermal expansion is also considered:

T(T ) = 1 2 2πf b a bcc (T )K(T ) P 1 (T ) = 2πf b a bcc (T ) E δ-c (0) + K(T ) ln 4f b a bcc (T ) eδ (5.15)
• In the third model, temperature dependence of the core-traction energy is also included with the temperature as: 5.16) In all the three models, the parameter K(T ) is calculated, using Eq. 5.9, through the elastic constants at the temperature T .

T(T ) = 1 2 2πf b a bcc (T )K(T ) P 1 (T ) = 2πf b a bcc (T ) E δ-c (T ) + K(T ) ln 4f b a bcc (T ) eδ (
For the model 3, an additional temperature dependence is given by the temperature dependence of the core-traction energy. However, the complete first principles deduction of core-traction free energy is inaccessible. Even if empirical potentials are used, the direct determination of the core-traction free energy is challenging and requires the free energy of very large simulation boxes, a task difficult to achieve even with the latest, fast and parallel computers. Hence, some assumptions about this term will be made, the accuracy of these approximations will be numerically demonstrated later in this chapter. The core-traction free energy is associated with the free energy of a tube of radius δ enclosing the edge of dislocation loop. The magnitude of the core-traction free energy stored in this tube is fixed by the zero K energy while the temperature dependence is controlled by an effective vibrational entropic term. This dependence is given by the usual form of free energy:

E δ-c (T ) ≡ F δ-c (T ) = E δ-c -T S δ-c
( 5.17) where, S δ-c denotes the effective formation vibrational entropy of the tube, per unit length.

In order to estimate the effective formation energy, S δ-c three assumptions have been made: (i) Firstly, the main contribution to the formation entropy comes only from the atoms of tube which have an atomic environment different from the bulk i.e. only the atoms which are in the border of the loop. This assumption is entirely consistent with the derivation of discrete term of the discrete continuum model presented in the Sec. 4.2.1. (ii) Secondly, the formation entropy given by the self-interstitials located in the border of the loops are treated in Einstein approximation, i.e. only local contributions are dominant. (iii) Thirdly, the value of the formation entropy for each interstitial in the border of loops is considered as the formation entropy of single self-interstitial in bcc matrix, 100 or 111 , for 100 and 1/2 111 loops, respectively. Using the last assumption, the finite temperature extension of the discrete continuum model remain entirely parametrized using first principles or experimental observables. Also, the last assumption implies that even the self-interstitials located in the corner of loops contribute with the same vibrational entropy as those located in the edges. In iron, the formation entropy S 111 f , of 111 self-interstitial has been already computed in DFT by Lucas and Schäublin (69 ) as 4.2 k B . For 100 self-interstitial, the value of S 100 f is 1.8 k B , as computed in this thesis. Using these values for the perimeters depending on the geometrical structure of the loops and using the same neighbourhood analysis, as in the previous chapter, it can be easily estimated that the effective vibration entropies S δ-c are 1.3 k B / Å and 2.1 k B / Å for 100 and 1/2 111 loops, respectively. The Fig. 5.1 sketches the free energies obtained using the three models of free energy discussed above. Even at elevated temperatures, the differences between the three models are less than 10 % which attests the hypothesis that the main influence of temperature is included in the temperature dependence of the dominant term T. The source 5.2 Finite Temperature Extension of Discrete-Continuum model of the comparatively small difference between the models 2 and 3 is the core-traction free energy. This denotes the importance of the parametrization of the function proposed for the E δ-c (T ) to detail the description of the free energy landscape of dislocation loops. Similarly, formation energy for the C15 clusters is derived in Sec.4.2.2 as:

E C15 f ormation (n) = 2γ s 9πΩ 2 1/3 n 2/3 + 12Ωµε 2 α n + 3n(E bcc coh -E C15 coh ), (5.18)
which can be manipulated in the same way, as demonstrated for 2D loops, to include temperature dependence of the parameters as:

F C15 (n, T ) = 2γ s 9πΩ(T ) 2 1/3 n 2/3 + 12Ω(T )µ[C ijkl (T )]ε(T ) 2 α[C ijkl (T ), B C15 (T )] n + 3n(E bcc coh (T ) -E C15 coh (T )). (5.19)
The finite temperature extension of discrete-continuum model for C15 clusters is more robust. The only variable which is fitted from the zero K calculation is the γ s term.

5 Finite-temperature extension of discrete-continuum model

The interfacial term γ s is assumed independent of temperature change. The required parametrization of all the variables in Eq. 5.19 can be obtained from experiment (bcc elastic constants or the thermal expansion) or can be deduced from first principles calculations (variation with the temperature of the elastic constants of the C15 phase, cohesive energies and the mismatch ). As indicated in the previous chapter, the above equation utilizes the isotropic implementation. We have implemented also the anisotropic variant by replacing the contribution of the second term of Eq. 5.19, by the energy given by spherical anisotropic Eshelby inclusion.

The key ingredient for the finite temperature extension of discrete continuum model is consideration of the temperature dependence of elastic constants for each bcc metal. An extensive literature survey provides us with the experimental values of the elastic constants at 0K and/or higher, finite temperatures. These experimental values are fit with appropriate polynomials of temperature to obtain an expression for calculation of elastic constants of each bcc metals. The experimental values and the curves obtained by fitting a polynomial in temperature are shown in Fig. 5.2. Using this temperaturepolynomial expression, Eqs.5.10 and 5.19 are evaluated to obtain the energy landscape at finite temperatures. However, we should note that all the elastic constants which come from the experiment are adiabatic elastic constants while the present free energy models assume a thermal equilibrium between the loops and bulk matrix which requires the iso-thermal elastic constants. The difference between the two types of elastic constants are presented in the following subsection.

Adiabatic and iso-thermal elastic constants

In order to define the elastic constants at finite temperatures the free energy can be written as Taylor series in strain components: .20) From the Helmholtz free energy, the isothermal elastic constants can be defined as follows: while the adiabatic elastic constants are defined in the same manner but at constant entropy. Using the definition of free energy, we obtain: 5.22) where U is the internal energy of the system. Alternatively, the adiabatic elastic constant can be defined as the Taylor development of the internal energy of the system: .23) In the same manner the isothermal bulk modulus B T and the adiabatic bulk modulus B S are defined as derivatives of the free energy with volume at constant temperature and entropy, respectively:

F ( , T, V ) = F (T, V ) + V ij σ ij ij + 1 2 ijkl ∂ 2 F ∂ ij kl T ij kl + . . . ( 5
C T ijkl = 1 V ∂ 2 F ∂ ij kl T , ( 5
C S ijkl = 1 V ∂ 2 F ∂ ij kl S = 1 V ∂ 2 U ∂ ij kl S , (
U ( , T, V ) = U (T, V ) + V ij σ ij ij + 1 2 ijkl ∂ 2 U ∂ ij kl S ij kl + . . . ( 5
B T = V ∂ 2 F ∂V 2 T (5.24) B S = V ∂ 2 F ∂V 2 S (5.25)
While the adiabatic bulk modulus B S is typically measured in experiments, it is more straightforward to calculate the isothermal bulk modulus in theory. The two variants of the bulk modules coincide at T=0 K and deviate at higher temperatures. Using simple equilibrium thermodynamic relations, the difference between adiabatic and isothermal elastic constants can be related to this deviation of the bulk modulus: .26) or in Voigt notations:

C S ijkl -C T ijkl = δ ij δ kl (B S -B T ) ( 5
C S 11 -C T 11 = B S -B T (5.27) C S 12 -C T 12 = B S -B T (5.28) C S 44 -C T 44 = 0 (5.29)
Using thermodynamics relations, generalized to include elastic strain energy, Hearmon (265 , 266 ) have shown that the difference between isothermal and adiabatic elastic compliances(stiffnesses) decreases(increases) linearly with temperature. For most mechanical purposes the difference between the two classes of elastic constants are negligible, being few percent or less for most of the cases. However, this difference can be important at high temperatures.

Finite Temperature Extension of Discrete-Continuum model

The distinction between isothermal and adiabatic elastic constants arises naturally from measurement methods followed in experiments: slow or static loading experiments give isothermal elastic constants while rapid measurements, such as ultrasonic wave determination of elastic constants, yield adiabatic constants because the system does not get enough time to achieve thermal equilibrium. In order to have access to the free energy using the finite temperature version of discrete-continuum-model, it is suitable to use isothermal elastic constants. These isothermal elastic constants are related to the experimental adiabatic elastic constants via thermodynamic relations such as (267 ): 5.30) where

B S B T = C P C V = 1 + αT γ = 1 + α 2 T V B T C V (
C P,V = T ∂S ∂T P,V
, is the specific heat at constant pressure/volume (5.31)

α = 1 V
∂V ∂T P , is the volumetric thermal expansion coefficient (5.32)

γ = V ∂P ∂V V = αV B T C V
, is the Grüneisen parameter (5.33) This relation in Eq. 5.30 enables calculation of isothermal bulk modulus from experimental adiabatic bulk modulus, thermal expansion coefficient and the Grüneisen parameter.

However a simpler method for calculating the isothermal bulk modulus is to use the Maxwell relation: (5.34) and substitute

C P -C V = T V α 2 B T ,
C V = C P B T B S
from Eq. 5.30 to directly obtain B T as: (5.35) where C P is the isobaric heat capacity which is very well known experimentally.

B T = C P B S C p + T V α 2 B S ,
Consequently, using Eq. 5.26 along with the experimental isobaric heat capacities C P , the volumetric thermal expansion α and the adiabatic bulk modulus B S , we can access the isothermal bulk modulus (as seen in Fig. 5.3) and other elastic constants. The first three rows of graphs in Fig. 5.3 show the variation of the volume, volumetric expansion coefficient α and the experimental isobaric heat capacities C P with temperature for Fe, W, Ta and V. The last row represents the ratio B S /B T calculated using Eq. 5.35. Here, empty symbols denote the experimental values while the lines represent an extrapolation. Units: Volume is in Å3 , α in 10 -6 K -1 , specific heat in Jmol -1 K -1 and temperature in Kelvin.

5.3 Atomistic models for interstitial defects at finite temperature

Atomistic models for interstitial defects at finite temperature

Canonical ensemble

The first step towards an appropriate statistical treatment of any system is construction of a statistical ensemble in accordance with the physical conditions involved. Depending on the known thermodynamic quantities, three principle ensembles are microcanonical, canonical and grand-canonical ensembles (270 ). This classification of ensembles is not exhaustive and other ensembles can be defined based on the system of interest. In our study, canonical ensemble is the most appropriate choice because a finite, specified temperature is assumed instead of a constant, known energy. Following is a brief recapitulation of the main results for canonical ensemble.

By the second law of thermodynamics, a system in the canonical ensemble is at equilibrium when the entropy is maximized or the Helmholtz free energy is minimized. In general, the probability function of a canonical ensemble or the canonical distribution can be expressed as:

P α = e -βEα Z , (5.36)
where the exponential term e -βEα denotes the Boltzmann factor with β = 1 k B T , k B is the Boltzmann constant, T is the temperature of the system and the denominator Z = Σ α e -βEα is the partition function which denotes sum over all states. Now, the probability of a system in a canonical ensemble having energy in the interval E to E +dE is P (E) = Σ α P α where α refers to the states with energy lying in the concerned energy interval or equivalently, P (E) = cΩ(E)e -βEα , where Ω(E) is the number of states with energy is the concerned energy interval. In principle, if the canonical partition function Z or the number of states Ω(E) can be calculated, then it becomes possible to derive the other thermodynamic quantities and to obtain the required thermodynamic quantities of interest for that system. It can be shown that choosing canonical partition function Z provides mathematical convenience and eliminates cumbersome counting of states for Ω(E) (271 ), especially when dealing with macroscopic systems.

Under certain conditions, the problem can be treated classically by replacing summation over ( r i , p i ) states by integration over d 3 rd 3 p phase space. When considering N indistinguishable particles in a continuous phase space, the canonical partition function 5.3 Atomistic models for interstitial defects at finite temperature less than the dissociation energy. Being the quantum analog of the classical harmonic oscillator problem, quantum mechanical solutions can be derived by reinterpretation of the classical solutions. This analogy somewhat simplifies the problem at hand. Let the adiabatic potential energy of a system of N atoms be denoted by U (r, r 0 ), where r = ( r 1 , r 2 , ..., r N ) refers to the set of positions of N atoms and r 0 = ( r 1,0 , r 2,0 , ..., r N,0 ) refers to their equilibrium positions. Given the position r i of the i th atom, the displacement from equilibrium of the i th atom is given by u i = r ir i,0 , where r i,0 is the equilibrium position of the corresponding atom. To consider sum of N vectors over the three (x, y, z) coordinates, we consider i to vary from 1 to 3N such that u 1x , u 1y , u 1z , ..., u N x , u N y , u N z = {u i } i=1,3N = u 1 , u 2 , u 3 , ..., u 3N -2 , u 3N -1 , u 3N . The potential energy of N atoms can be expanded into a Taylor series about its equilibrium position as:

U (r, r 0 ) = U (r 0 ) + i C i u i + 1 2 i,j C ij u i u j (5.42) + 1 6 i,j,k C ijk u i u j u k + • • • , (5.43)
where

C i = ∂ U (r, r 0 ) ∂ r i r=r 0 , (5.44
)

C ij = ∂ 2 U (r, r 0 ) ∂ r i ∂ r j r=r 0 , (5.45)
C ijk = ∂ 3 U (r, r 0 ) ∂ r i ∂ r j ∂ r k r=r 0 , etc.
The coefficients C ij , C ijk ,... are the atomic force constants of the second order, third order,..., respectively.

The harmonic approximation (HA) holds for small displacements compared to the nearest neighbor distance and it implies that the third and higher order terms in the Taylor expansion can be neglected from Eq. 5.43. The first term in Eq. 5.42 is a constant 157 5 Finite-temperature extension of discrete-continuum model and defines the equilibrium potential energy of the system. Since the Taylor expansion is done at the minimum energy of the system, the first derivative of the potential at equilibrium has to be zero i.e. C i = 0 for all i. Hence, the harmonic potential U (r, r 0 ) will depend quadratically on the equilibrium position r 0 due to the non-zero third term in Eq. 5.42. At high temperatures, when the displacements u i are too large, it becomes necessary to consider anharmonic terms ∼ r n (n > 2) too.

Using harmonic approximation, the classical Hamiltonian of N atoms of a system can be written in the following form: (5.46) and the equations of motion for these N atoms are

H = 1 2 i m i u2 i + 1 2 ij C ij u i u j ,
m i üi = - 1 2 j [C ji + C ij ] u j .
(5.47)

The above linear equation relates the force m i ri to the displacements u j through the force constant C ij . These constants can be: i) calculated directly from the derivatives of the interatomic potential, empirical potential or ab initio, using Eq. (5.45), or ii) fitted from experimental elastic constants and phonon frequencies.

The Hamiltonian of Eq. ( 5.46) is a quadratic in displacement {u i } i=1,3N . A feasible method to solve such coupled differential equations of motion is to assume a linear transformation of displacement of the kind: (5.48) in order to diagonalise the quadratic form of Hamiltonian in Eq. 5.46 into the form: .49) This new set of coordinates U p is known as normal coordinates. In phase-space, the component U p defines the motion of a simple harmonic oscillator of frequency ω p and mass 1. Substituting the solution from Eq. 5.48 in Eq. 5.46, one gets:

u i = 1 √ m i 3N p=1 L ip U p .
H = p 1 2 Up 2 + 1 2 ω 2 p U 2 p . ( 5
H = 1 2 i p L 2 ip U2 p + 1 2 j,i p,p 1 √ m i m j C ji L jp L ip U p U p = 1 2 i p L 2 ip U2 p + 1 2 p,p L -1 D L p,p U p U p .
(5.50)

5.3 Atomistic models for interstitial defects at finite temperature where the D matrix defined by D ij = C ij / √ m i m j is known as a dynamic matrix. Now, Eq. 5.48 can represent the solution only if Eq. 5.50 can be converted into the required, diagonalised form in Eq. 5.49 by imposing the necessary conditions on the unknowns of the supposed solution. By comparison of the two equations, it can be concluded that L -1 D L should be diagonal and L should be a unitary matrix which diagonalizes D.

Formally, one has: (5.51) where the eigenvalues ω p of the dynamical matrix represent the normal frequencies associated with the normal modes. The eigenvectors of the dynamical matrix D can be taken to satisfy the usual orthogonality conditions:

L -1 D L p,p = ω 2 p δ p,p ⇔ j D ij L jp = ω 2 p L ip ,
L ip L ip = δ p,p (5.52)
The orthogonal nature of the matrix L is imposed by the quadratic form of the kinetic energy and the symmetric form of the matrix C. If L satisfies the Eqs.5.51 and 5.52, only then does its form preserve kinetic energy and potential energy forms simultaneously.

If the system has a periodical boundary condition, it can be shown that the three frequencies are exactly zero and correspond to the modes which represent a translation of the overall system. Taking into account the oscillatory behavior of the normal coordinates, the small displacements u i can now be rewritten as:

u i = 1 √ m i p L ip U p (0)e -iωpt , (5.53)
The coordinate U p (0) is the amplitude of the p th normal mode. Note finally that the symmetry considerations can drastically reduce the present 3N × 3N dimensions of the dynamical matrix D.

Following (5.45), the C matrix is symmetric. The equilibrium conditions of the system ensure that it represents a convex function due to minimization condition or, equivalently, that the matrix D is positively defined. This implies that ω 2 p ≥ 0 or that the frequency ω p is real.

Beyond normal modes analysis, we should pay attention to the anharmonic terms in Eq. (5.43). The anharmonic terms introduce couplings among normal modes which can be treated as correction terms. For instance, these terms allow the thermal equilibrium to be reached, a process which could not occur if the modes were completely independent.

5 Finite-temperature extension of discrete-continuum model

Harmonic free energy

Using harmonic approximation, the partition function can be easily deduced using the energy of independent harmonic oscillators. For a single harmonic oscillator with the frequency ω p , the canonical partition function at temperature β reads:

Z p = n exp -β ω p n + 1 2 = exp -1 2 β ω p 1 -exp (-β ω p ) = 1 2 sinh β ωp 2 .
(5.54)

The partition function for a harmonic solid with 3N independent oscillators is the product of these individual oscillator partition functions. Consequently, the free energy of the harmonic solid can be written as: 5.55) or

F vib = -β -1 ln p Z p = 1 2 p ω p + β -1 p ln [1 -exp (-β ω p )], (
F vib = β -1 p ln 2 sinh β ω p 2 .
(5.56)

The first term of the previous equation is known as zero-point energy i.e. the canonical free energy of the system in the limit of T → 0. The vibrational entropy becomes:

S vib = - ∂F vib ∂T = k B p -ln [1 -exp (-β ω p )] + β ω p exp (β ω p ) -1 , (5.57)
or .58) and the vibration internal energy is calculated from the expression:

S vib = k B p -ln 2 sinh β ω p 2 + βω p 2 coth β ω p 2 , ( 5
U vib = - ∂ ln p Z p ∂β = p ω p 2 coth β ω p 2 .
(5.59)

From the above-mentioned results for F vib , U vib and S vib , it can be verified that the thermodynamic relation F vib = U vib -T S vib holds good. It is worth noting that the entropy at temperature β of phonon with energy ω p can be rewritten as:

S p = k B [(1 + n BE ) ln (1 + n BE ) -n BE ln(n BE )] , (5.60)
where n = n BE (β, ω p ) gives the number of particles with the corresponding energy and at the given temperature in Bose-Einstein distribution:

n (β, ω p ) = 1 e β ωp -1 .
(5.61)

Atomistic models for interstitial defects at finite temperature

In the case of fermions, entropy has a form similar to Eq. 5.60 but with a negative sign before k B and number of particles from Fermi-Dirac distribution for fermions instead of Bose-Einstein distribution.

Classical harmonic free energy

In the limit of high temperatures, the Bose-Einstein distribution function becomes Boltzmann distribution and the phonons behave as classical oscillators system. In this limit, the free energy and entropy of a classical solid become:

F vib = β -1 p ln (β ω p ) (5.62) S vib = -k B p ln k B T ω p + 1 (5.63)
In the case of solids, the limit of high temperature is given by Debye temperature at which all phonons of the solids are excited.

Quasiharmonic free energy surface

A pure harmonic limit assumes that the vibrational frequencies ω p do not depend on the interatomic distances. Consequently, the vibrational free energy F vib depends only on T without any dependence on V . To overcome this drawback, the harmonic approximation can be generalized to the quasi-harmonic approximation where the pulsations ω p are assumed to depend on volume. This situation is consistent with any strain or deformation applied to the solid. Let a constraint, say X, be applied to the system. This constraint can be volume V or the strain tensor or any externally applied load. Now, the pulsations ω p will depend on the constraint X and the Helmholtz free energy can be written:

F QHA (X, T ) = U (X) + 1 2 p ω p (X) + β -1 p ln {1 -exp [-β ω p (X)]} . (5.64)
The U (X) is the energy of the crystal or the free energy at 0 K while the second and the third terms denote the zero-point energy and the thermal contribution, respectively. If X = V , then the above expression of the QHA free energy gives the equation of state of the crystal: 5.4: The values of lattice parameter (a 0 ), isothermal bulk modulus (B) and the free energy (F ) of the bcc phase of iron ) from QHA approximation implemented in this thesis and force brute calculation (as it is explained in the text). The solid line and circles emphasize the QHA and brute force calculations, respectively. The phonon gas is quantized (blue) or treated in classical approximation (orange).

P (T, V ) = - ∂F ∂V = - ∂U ∂V + 1 V p ω p (V )γ p 1 2 + 1 exp [β ω p (V )] -1 , ( 5
where γ p 's are known as the Grüneisen mode parameters and are defined as:

γ p = - V ω p (V ) ∂ω p (V ) ∂V (5.66)
The Grüneisen mode parameters γ p are assumed to be independent of the volume V . These parameters can be computed by a series of ω p (V ) calculations, using harmonic approximation at various volume around equilibrium. In general, 2-4 calculations are necessary to provide a good estimation of the first partial derivatives in Eq. 5.66.

Once the Grüneisen mode parameters γ p are computed, any thermodynamics property of the system which depends on the volume can be calculated. For example, thermal expansion becomes accessible using one of Maxwell's thermodynamics relations as: 5.67) or by using the derivation in Eq. 5.65, the thermal expansion can also be written as:

α(T ) = 1 V ∂V ∂T P = - (∂P/∂T ) V (∂P/∂V ) T = 1 B T ∂P ∂T V , (
α(T ) = 1 B T p ω p γ p ∂ ∂T 1 exp [β ω p (V )] -1 . (5.68)
A more precise, though not more efficient way to take into account the volume dependence of the harmonic free energy is to fit the full free energy surface F (T, V ) using a formula of the kind: 5.69) This recipe was proposed by the group of J. Neugebauer et al.(272 ,273 ) and calculations involving 7 or 10 phonons at different volumes are needed to achieve high accuracy c i (T ) coefficients.

F QHA (T, V ) = 3 i=0 c i (T )V i (
Moreover, using Eq. 5.2, the equilibrium volume of a solid at given temperature T and pressure P can be deduced by solving the following equation for V :

P = P (V, T ) = - ∂F ∂V T (5.70)
If zero pressure condition is imposed in the previous equation, the volume corresponds to the equilibrium volume V eq (T ) at the temperature T . Then, computing this equilibrium volume at many temperatures we can have access to volume expansion coefficient. Also the corresponding free energy at zero pressure becomes F (T, V eq (T )).

5.3 Atomistic models for interstitial defects at finite temperature to the artefacts of empirical potential and have their origin in the mentioned differences between quantum and classical phonon gas. The situation is very interesting for the bulk modulus: most of the studies in literature (e.g. DFT calculation) report the bulk modulus which are fitted from the energy-volume curve. Using the free-energy curves we show a difference of up to 2 GPa between the two approximations, the correct value being the quantum value. Thus, a potential (empirical or pseudo-for DFT calculation) which perfectly fit the experimental value does not necessarily translate into a very good potential. Probably, the future comparison of experiment-DFT calculations should take this difference into account.

In order to access the isothermal and adiabatic elastic properties, quantities such as heat capacities C P,V of Eq. 5.31 should be provided. This can be achieved by computing the entropy of system and its partial derivatives at constant pressure or volume, respectively. In harmonic approximation, it is very simple to access entropy using classical, Eq. 5.63, or quantum approximation, Eq. 5.57. In order to achieve the volume dependence, we fit the entropy computed at 7-10 volumes around equilibrium with the same function in volume as Eq. 5.69. The equilibrium value of the lattice parameter is already deduced for each temperature of interest, hence, the function S(V ) and its derivative (∂S/∂T ) V ) should be evaluated at one point for deducing C V (T ) at the equilibrium volume. C P (T ) is then obtained from C V (T ) using Maxwell relation, Eq. 5.34. Using this procedure, we calculate a0, C P,V for pure bcc and C15 structure of iron using M10 potential (100 ), as shown in Fig. 5.5.

Through the assumed volume dependence of vibrational frequencies of the crystal in QHA, all equilibrium thermal properties of the system become immediately accessible. This approach takes into account the anharmonic effects implicitly through the volume dependence of the vibrational frequencies.

Elastic constants at finite temperatures using QHA approximation Determination of elastic constant at finite temperature is quite straightforward, using the same procedure as well established for 0 K (277 ). Using QHA approximation, we determine the equilibrium volume of the solid for a specific temperature. For this volume, we distort the primitive vectors using a well-defined strain tensor. Then using the development of the free energy in term of strain given by Eq. 5.20 we obtain the The phonon gas is quantized (blue) or treated in classical approximation (orange).

5 Finite-temperature extension of discrete-continuum model elastic constants. Let's consider a general deformation with the strain . The crystal primitive vectors a i are transformed into vectors a i of deformed crystal as:

   a 1 a 2 a 3    =    a 1 a 2 a 3    (I + ) (5.71)
where the I is 3×3 identity matrix and is the symmetric strain tensor. We consider only non-rotating strain which also preserve the constant volume. If we take the following two deformations:

11 =    δ/2 0 0 0 -δ/2 0 0 0 δ 2 /(1 -δ 2 )    , 44 =    0 δ/2 0 δ/2 0 0 0 0 δ 2 /(4 -δ 2 )    ,
then the difference in free energy at a given volume V for both deformations becomes an even function in δ which can be written as:

∆F 11 (δ) = V (C 11 -C 12 ) δ 2 + O[δ 4 ], (5.72) ∆F 44 (δ) = 1 2 V C 44 δ 2 + O[δ 4 ]. (5.73)
By performing a series of deformations around the equilibrium volume given by QHA approximation, the above equation can be fitted in order to obtain the curvature. For the Taylor expansion of free energy in strain to hold good (Eq.5.20),the value of maximal strain should be small. Here, we use the maximal strain less than 1%. Using also the value of the bulk constant, deduced previously in QHA approximation, we have all the elastic constants, C 11 , C 12 and C 44 at all temperatures of interest.

As described above, we show the elastic constant of bcc phase computed using QHA approximation in Fig. 5.6 with empirical potentials for Fe (99 ) and W (135 ). At low temperatures, we find the intrinsic differences between classical and quantum approximation due to the zero point energy. We should point the anomalous behavior of W, for which, using the present EAM potential, the isothermal elastic constants increase with the temperature. This behavior is probably an artefact of the empirical potential. The experiments, summarized in Fig 5 .2 suggest this deficiency of the present EAM potential for tungsten. In Fig. 5.7 we have compared the results from the three free energy model, entirely parametrized using potential, to the atomistic calculations using the same potential. In the case of loops the model 1 and the model 2 are very close and give large discrepancies compared to atomistic calculations. However, the 3 rd model is very close to the atomistic calculation. The precision is high with error less than about 3 % in the range of temperature 0 -1000 K. As expected, the form of the function proposed for the dependence in the temperature of the core-traction energy is crucial in the prediction on the discrete-continuum model. The same agreement holds for the case of anisotropic version of discrete-continuum for the C15 clusters, as it shown in 5.8.

The comparison of both models, discrete-continuum and atomistic, indicates that the discrete-continuum model is able to predict, with striking precision and without any adjustable parameter, the free energies given directly from atomistic calculation. The accuracy of predictions of discrete continuum is remarkable, even at higher temperatures. The present validation opens up many possibilities in the investigation of the free energy landscape of mesoscopic dislocation loops. In the following section of this chapter, we will present only two of the many possible applications of the extended version of the discrete-continuum model.

Application: mechanism of formation of 100 loops in bcc iron

At zero K, using discrete-continuum model, we have shown in the previous chapter that the 1/2 111 loops are always the most stable family of SIA-clusters above ∼ 50 SIAs. At lower sizes < 50 SIAs, the C15 clusters are the most stable clusters. In the intermediate regime the C15 clusters have energies lying between 1/2 111 and 100 loops. Above ∼ 100 SIAs the traditional dislocation loops, 100 or 1/2 111 are more stable than 3D C15 cluster. The energy landscape obtained, sheds some light on the absence or low concentration of 100 loops in low temperature experiments, and reconciles the Zhang mechanism (112 ) with the experimental evidence. Our findings support the theory of formation for 100 loops proposed by Zhang et al. (112 ) which completes the panel of mechanisms proposed by Marian (79 ) and Xu and Terentyev (105 ). These 5.4 Application: mechanism of formation of 100 loops in bcc iron mechanism are described in detail in Sec. 1.2.2 and Sec.4.5.1. However,our The free energy predictions are needed to explore high temperature limit. These free energy calculations should include thermal effects and magnetic excitations, explicitly. Our finite temperature discrete continuum model does not include the magnetic excitation explicitly but these excitation are included through the temperature dependence of the experimental elastic constants. The finite temperature discrete continuum model for loops include these excitation entirely. However, for C15 clusters, the observables needed by the finite temperature discrete continuum model are not available from experiments. These observables should be computed using DFT calculations.

Treating the magnetic excitation, even in ab initio method is very challenging and this field is subject to continuous development. Consequently, the bulk properties of the C15 5 Finite-temperature extension of discrete-continuum model 0 300 600 900 6.64 isothermal elastic constants of the ferromagnetic C15 phase of iron. The circles denote the temperatures chosen for QHA approximation, the black line being a interpolation of those points.

5.4 Application: mechanism of formation of 100 loops in bcc iron phase will be computed only using the ferromagnetic approximation. The thermal effects will be treated using the previous developed QHA approximation. As we have pointed out in the previous sections the quantity needed to complete the experimental data are temperature dependence of the lattice parameter, cohesive energy and the elastic constants of the C15 phase. In order to validate this approach, we have computed the corresponding quantities also for the bcc phase which can be directly compared to the experiment.

In order to perform phonons calculations we have chosen PWSCF code as part of the Quantum-Espresso package (278 ). This choice is motivated by the possibility of computation of the second derivatives of the potentials energy surfaces Density-Functional Perturbation Theory. Using this approach, the force constants can be deduced directly in the reciprocal space. By Fourier transform and using a uniform grid mesh in the Brillouin zone (4×4×4 in these calculations), the force constants can be recovered in the direct space. The major advantage of this technique is that the simulation cell can be restricted only to the unit cell i.e. one atom per cell and six atoms per cell for bcc and C15 phase, respectively. DFT calculations were performed using the Generalized Gradient Approximation and the Projector Augmented-Wave scheme with semicore states in the electronic configuration of iron. Energy cutoff of plane-waves was set to 90 Ry. The Hermite Gaussian broadening was set to 0.04 eV. The choice of the K-point grid is extremely important because lower grid could yield results for an order of magnitude higher or lower than the converged results. 24×24×24 and 14×14×14 mesh grids are used for bcc and C15 respectively. In addition to vibrational free energy, we have also included the electronic free energy, which is very easy to compute once the Kohn-Sham orbitals are computed.

For QHA approximation, we have used 10 volumes in order to fit the volumetric free energy surface. Once the equilibrium volumes are deduced, we have computed the elastic constants using the procedure described in Sec. 5.3.1: 14 temperatures have been chosen, the deformation consists of 6 strain values at each temperature.

The results are displayed in the Figs. 5.9 and 5.10, for ferromagnetic bcc and C15 phase of iron. In bcc iron the brute DFT results seems to be shifted compared to experiment. This type of behaviour is expected due to systematic errors of the exchangecorrelation functional in the estimation of the lattice parameter and the elastic constants at zero K. However by scaling DFT results with the experimental values, the DFT prediction get closer to the experimental results for B and C 44 . The discrepancies for 5.4 Application: mechanism of formation of 100 loops in bcc iron lattice parameter and the C have their origin, in the anharmonic vibrational contribution and the magnetic excitation. Here, bcc Fe is treated in the ferromagnetic state of iron and the present QHA approximations, so these DFT calculations cannot predict the abrupt decrease of C' near the α → γ transition. The same errors are expected to impact the calculations of the corresponding observables for the C15 phase also. Consequently, some errors are expected in the prediction of the discrete continuum approach for the C15 phase.

To resume for loops we have used experimental values for lattice parameters and elastic constants as well as the DFT values of formation entropies of 111 and 100 mono-interstitial as inputs. In the case of C15 the parametrization is done by ab initio only. All the inputs of the finite temperature discrete continuum models being set, the prediction about the crossover between the free energies of 100 and 1/2 111 loops and the C15 clusters are presented in the Fig. 5.11 by full line.

We have tested the sensitivity of the model to variation with the temperature of the elastic constants of C15 phase. Beyond the full ab initio prediction and anisotropic elasticity, emphasized by full line in Fig. 5.11 we have tested different temperature dependence of the elastic constant of C15 phase. The experimental isothermal elastic constant of the bcc phase is rescaled to zero K DFT values of C15 elastic constants and used as input in the isotropic or anisotropic model. The obtained values are compared with the three energy models for the free energy of the dislocations loops. All these predictions are included in the errors bars presented in the Fig 5 .11. However, even after this error analysis the question of stability in temperature of the C15 phase, remains. Theoretical or experimental answers for this question is a challenge and probably should be addressed in the future.

The main information from the Fig 5 .11 is that there is a crossover between loops and C15 clusters for the entire range of temperature between 0 K and the temperature of αγ transition. This implies two important conclusions:

• Firstly, the C15 clusters can be observed for higher temperatures as they remain the most stable clusters up to 900 K, near the transitions αγ

• Secondly, due to this thermal stability of C15 clusters, the mechanism of formation of 100 proposed by Zhang can operate in any range of temperature. Moreover, the difference in crossover sizes of C15 clusters with the two types of dislocation loops decreases with increasing the temperature impacting the possibility of formation of clusters 100 . Due to this, the concentration of 100 becomes higher and higher with increasing temperature.

As we already mentioned we do not claim that the present findings will replace the existing accepted mechanism of formations of 100 loops proposed by Marian (79 ) and Xu and Terentyev (105 ). We claim that the present mechanism can occur with much higher probability than the previous mechanism. For instance, let us consider the experiments mentioned in section 1.2.2. According to Marian (79 ) and Xu and Terentyev (105 ), migration of 1/2 111 loops should produce 100 loops under some specific conditions. Experiments proved that Cr can slow down the migration of 1/2 111 loops in irradiated FeCr. This implies that the 100 loops should not be formed. However, 100 dislocation loops exist in F eCr alloys. The scenario proposed in this thesis can explain the above experimental results. Preliminary calculations of the C15 clusters in iron suggest increased stabilization of C15 clusters in the presence of Cr (279 , 280 ).

Application: concentration of dislocation loops post irradiation

In this section, finite-temperature free energy derived from discrete-continuum model is used to perform cluster dynamics modeling. These results, along with cluster dynamics results from the widely-used capillary law, are compared with experimental results.

In 2015-16, A. Duchateau and E. Meslin performed an irradiation experiment with an iron sample in helium atmosphere to study the density and size of dislocation loops. They used a polycrystalline ultra-high purity α-iron discs of 100 µm thickness and 3 mm diameter. These discs were irradiated with 60 keV He ions at room temperature, producing an estimated dose of 0.5 dpa. This He irradiation of iron samples is done to prevent the non-conservative Ostwald ripening process (249 , 281 ) which is caused by presence of powerful vacancy sources, like vacancy clusters and surfaces of thin samples, capable of annihilating interstitial loops. After irradiation, at room temperature, isochronal annealing at 3 K/min is carried out up to 830 K which is followed by isothermal annealing at 830 K, as shown in Fig. 5.12. Standard diffraction contrast methods using bright field technique were used to determine the morphology and nature of point defect clusters. The dislocation loops were found to be interstitial 1/2 111 loops. The mean radius and surface of these dislocation loops were also measured with time. The blue dots in Fig. 5.13 show the experimental data points, after annealing during 50 min. The same distributions were recorded at 72 min and 100 min. The goal of this section is to compare the cluster dynamics results, based on the free energy model proposed in this chapter, with the experimental distribution of sizes of dislocation loops at 72 min and 100 min.

Cluster dynamics

Cluster dynamics modeling provides an efficient means to study time evolution of clusters representing the real system by using chemical rate theory. Whenever interstitial loops are present, vacancy clusters are also considered because vacancies can be created during annealing of interstitial loops (249 ). If concentration of clusters containing n interstitials (n ≥ 1) is denoted by C n , equation for immobile cluster (250 ) is given by: 5.74) where M is the set of mobile species and J n,n+m is a net reactive flux between cluster n and cluster n + m, due to the mobility of cluster m. This flux is defined by: (5.75) where absorption rate β n,m is determined as:

dC n dt = m∈M J n-m,n - m∈M J n,n+m , (
J n,n+m = β n,m C n C m -α n+m,m C n+m ,
β n,m = 4πZ n,m (r n + r m )D m , (5.76)
D n is the diffusion coefficient of n (D n = D 0,n exp (-F m n /kT )), r n is an effective radius of cluster n and Z n,m is an efficiency factor which accounts for long-range elastic interactions between defects n and m. The emission rate can be calculated as: .77) In this expression V at is the atomic volume and F b n+m,m is the free binding energy of clusters n and m, defined in terms of the cluster free formation energies F f n as:

α n+m,n = β n,m V at exp - F b n+m,m kT . ( 5
F b n+m,m = F f n + F f m -F f n+m . (5.78)
For mobile species, Eq. ( 5.74) becomes: .79) where Ω is the set of all species.

dC n dt = m∈M J n-m,n - m∈M J n,n+m - m∈Ω J m,n+m . ( 5
From the cluster dynamics equations described above, it can be seen that a number of parameters are required to carry out cluster dynamics modeling. These parameters are mentioned in the Appendix A1. Additionally, free formation energies are needed as input. Since the capillary model for cavities gives a good asymptotic value as well as a smooth transition for small vacancy clusters, it is widely used (107 , 250 ). For loops, the four models derived in this thesis (Model 1, 2 and 3 are same as those introduced in Sec. 5.2 while Model 0 is another parametrization where elastic constants considering temperature dependence are adjusted according to DFT calculations) are considered along with the popular capillary law which is very similar to the one for cavities, so the binding energy of an interstitial to a loop is given by:

E b n,1 = E f 1 - E f 1 -E b 2,1
2 2/3 -1 n 2/3 -(n -1) 2/3 .

(5.80)

5 Finite-temperature extension of discrete-continuum model

Experimentally, loops are shown to evolve due to an Ostwald ripening by emission and capture of vacancies, since emission of interstitials is not possible owing to their high formation energy (249 ) (Fig. 5.14). Therefore, the relevant parameter for this study is the free binding energy of a vacancy to an interstitial dislocation loop, which reads (n ≥ 1): For the capillary law, this binding energy can be deduced from Eq. ( 5.80), since:

F b n,-1 = F b (n+1)-1,-1 = F f -1 + F f n+1 -F f n . ( 5
F b n,-1 = (F f -1 + F f 1 ) -F f 1 + F f n+1 -F f n (5.82) = (F f -1 + F f 1 ) -F b n+1,1 (5.83)
For the four other models, it is computed directly from the values of F f n .

Binding energies of vacancies and interstitial loops given by the 5 models are compared in Fig. 5.15. Then it is transformed into a distribution as a function of the number of interstitials, according to

C n = C(r) dr dn , (5.88)
and assuming that loops are 111 spherical loops, so that πr 2 b = nV at with b = 0.165 nm.

The initial distribution is shown in Fig. 5.13.

Results

The cluster distributions at different times, t = 72 min and t = 100 min are compared in Fig 5 .16. We see that the capillary law leads to a very different result from the other models and deviates significantly from the experimental results. On the contrary, the other four models lead to results in good agreement with experimental results.

It is known that these results are highly sensitive to the sum of the vacancy formation and migration free energies. It is interesting to see what would be the value of this sum to obtain a good agreement with experimental results, when a capillary law is used. A fair agreement can be obtained if, for example, the migration free energy is lowered from 0.51 eV to 0.45 eV (Fig. 5.17). In turn, this leads to a self-diffusion coefficient which is clearly overestimated with respect to measurements (Fig. 5.18). This tends to invalidate the capillary law. It remains to check that results are indeed sensitive only to the mono-vacancy properties and to the energetics of loops. Sensitivity analysis is performed using model 0 as a reference. It is found that the results are not sensitive to individual values of F f -1 and F m -1 but the sum F f -1 + F m -1 . Moreover, mobility of mono-to tetra-interstitials and that of di-to penta-vacancies do not modify the results, as shown in Appendix A1.3. Thus, the robustness of our results is confirmed. We have shown that the finite temperature discrete continuum model predict the variation of the post irradiation loops concentration with the temperature compared to the direct experiments very well. Moreover, the precision of the discrete continuum model does not leave much arbitrariness in the choice of the parameter. In fact, only the vacancy activation energy can be considered as an arbitrary variable. On fitting vacancy property with the finite discrete-continuum model, the value is found to be the value predicted by the self-diffusion experiments but the value predicted by using capillary law shows some deviation. The present finite temperature discrete continuum is perfectly compatible with the self-diffusion experiments. Beyond this fact, the precision of our model, open up new experimental possibility such as the experimental determination of the activation energy of the vacancy from the present experiment.

Conclusions

Conclusions

Our work makes it possible to extend the zero K scaling laws for the formation energies of various types of clusters in various materials to finite temperature. This extension was done to preserve the possibility to be completely parametrized from first principles cal-185 5 Finite-temperature extension of discrete-continuum model culations. The present finite temperature extension of the discrete continuum approach has been validated by atomistic calculations using quasi-harmonic approximation.

The discrete continuum model is able to predict, with striking precision and without any adjustable parameter, the free energies given directly the atomistic calculation. The accuracy of predictions of discrete continuum is remarkable, even at higher temperatures. The present validation opens up many possibilities in the investigation of the free energy landscape of mesoscopic dislocation loops. Finally, we have presented two of the many possible applications.

As the results obtained from the zero K discrete-continuum model relate the low temperature experiments in Fe with Zhang mechanism, the finite temperature extension of this model supports the same mechanism for explanation of high temperature experiments in Fe. The experimental observation of 100 dislocation loops at high temperatures is justified using the finite-temperature extension of the discrete-continuum model and the Zhang mechanism. As a second application, free energy values predicted by the finite temperature discrete continuum model are used as input in cluster dynamics simulations in order to compare the simulation results with time evolution of loops density in post irradiation experiments. The comparison with our model shows better results than the commonly-used capillary law.

Validation of the 0 K discrete-continuum model and its extension for finite temperatures coupled with successful applications of both the versions of the model show promise for more possible implementations.

Conclusions and perspectives

A fairly-detailed literature survey of irradiation-induced effects on Fe, W, V and Ta has been compiled, recapitulating most of the relevant experimental results along with significant simulation as well as computation results. The conclusions corresponding to each transition metal is presented below:

• Fe: The number of experiments studying the irradiation-induced microstructural changes in iron is quite high compared to the other transition metals. All these experiments report a majority of mobile 1/2 111 interstitial loops at lower temperatures. Further, an increasing proportion of the comparatively immobile 100 is observed above 573 K which ultimately become the majority at around 773 K. The orientation of the sample is crucial in this observation, otherwise the 1/2 111 can migrate to surface, biasing the results of experiments. Moreover, the impurities can also bias the results impacting the mobilities of the different loop types. For example, He stabilizes the 1/2 111 loops rendering them immobile, as predicted by MD simulations. Cr in Fe-Cr alloys reduce the mobility of the 1/2 111 loops as well, eventhough the combined effect of He beam and Cr content in alloys did not show any clear trend. Though the theoretical explanation was lacking till recently, most of the experiments have now been explained rather convincingly in the existing literature. According to DFT Calculations, the most stable mono-SIA in Fe is the 110 dumbbell. The 110 SIAs are the most stable in Fe up to around five SIAs and 1/2 111 loops are the most stable beyond the cluster size of around five interstitials. Taking into account elastic approximation to calculate the formation energies near the αγ transition, a reduction of free energy of 100 loops was demonstrated, thus making the 100 loops more stable compared to 1/2 111 loops at temperatures higher than about 800 K. Although the stability of 100 was thus established, the mechanism for the formation of 100 loops has not been confirmed yet. The theoretically proposed mechanism of the formation of 100 by collisions of two 1/2 111 loops can explain the formation of 100 loops in 6 Conclusions and perspectives pure Fe, or Fe under He atmosphere. However, despite the fact that Cr reduce the mobility of 1/2 111 loops the formation of 100 loops is very active in F eCr alloys. Moreover, the absence of 100 interstitial dislocation loops at low temperatures remains unexplained, given that the 1/2 111 loops are very mobile at low temperature. Recent studies highlight the importance of C15 clusters in understanding and explaining the anomalous experimental observations of 100 loops in Fe, thus reconciling experiments with DFT calculations of energy landscape of point defects in bcc Fe.

• W: While recent experimental observations reveal presence of a vast majority of 1/2 111 and a few 100 interstitial loops at low temperatures, these 100 interstitial loops disappear leaving only 1/2 111 interstitial loops at high temperatures. Empirical potentials seem to be inconclusive with some predicting 1/2 111 as most stable and others show a crossover between 1/2 111 and 100 loops suggesting competing stability of these loops as a function of the number of interstitials. As such, MD simulations using different potentials can be contradictory as well but DFT calculations have quite definitively established the stability of 1/2 111 interstitial loops amongst other possibilities.

• V: In general, dislocation loops of Burgers vector 1/2 111 , voids and planar precipitates are observed after neutron irradiation. However, there are contradictions among the experiments with one reporting the presence of < 100 > interstitial loops in addition to the above-mentioned defects and another stating the nature of 1/2 < 111 > dislocation loops as vacancy-type. These contradictions have been ascribed to the varying levels of impurities among the different experimental samples. Ab initio calculations have established that the 111 interstitial dumbbell is the most stable and even predictions using interatomic potentials concur regarding the stability of 1/2 111 interstitial loops. In fact, simulations using Finnis-Sinclair potential for V predicts easy rotation of 100 and 110 loops into 1/2 111 due to a significant difference in their formation energies.

• Ta: Scarcity of experiments and contrasting irradiation conditions among the various experiments prevent any conclusive evidence about the type of defect encountered in Ta. While one experiment found screw dislocations with Burgers vector 1/2 111 , another reports growth of interstitial loops which is followed by a preferential formation of vacancy loops at higher temperatures. Nevertheless, from first principles calculations we can expect that 1/2 111 loops are the most stable.

The existing theoretical background of the atomistic methods (ranging from ab initio to EAM potentials via Tight Binding Approach) are revisited. Also, the various EAM potentials for bcc metals are discussed. Ideally, density functional theory (DFT) calculations would be employed to capture the complete physics involved in a system with irradiation-induced defects. However, unachievable computational requirements of DFT calculations have fueled the search for alternatives considering reasonable approximations over the past decades. This search has led to development of a number of empirical potentials, ranging from pair potentials to EAM.

Although many empirical potentials have been successful in making radiation damage studies feasible, inconsistency of results from different empirical potentials is a major shortcoming that hinders conclusive theoretical results. Thus, a new potential is developed for Fe specifically for point-defect simulation studies. A comparative study of the results from empirical potentials is also presented, justifying the need for a new EAM potential. Further, the construction of C15 clusters is illustrated with the selection rules. On comparison of the energy landscape of defects in Fe using existing and newly-developed empirical potentials, it is established that better results can only be obtained by going beyond empirical potentials and thus, a new model is proposed.

A new model is developed from 'discrete' ab initio formation energies for small-sized defects and 'continuum' anisotropic elastic theory formulation for formation energy of large-sized defects. The so-called discrete-continuum model makes possible ab initiolevel accurate calculations for clusters without any size limitation. The scaling laws predicted by the discrete-continuum model about the size dependence of the clusters energy are validated by comparison with empirical potentials spanning the very same scales within extensive atomistic simulations. The model allows us to treat various cases of interstitial dislocation loops and C15 clusters from clusters containing a few SIAs to nanometer size.

From the interpretation of the present results it can be concluded that above ∼ 100 SIAs 1/2 111 loops are always the most stable family of SIA-clusters -in agreement with experimental observations of irradiation defects at low temperature in bcc metals. However, these results are at odds with calculations made using various EAM interatomic potentials, which yield spurious predictions concerning the relative stability of 100 and 1/2 111 loops (100 ). Future developments of such potentials should consider the information provided in the present paper, and include the appropriate additional fitting conditions on the potential parameters. Our study shows that in Fe, C15 clusters are the 6 Conclusions and perspectives most stable clusters of defects for sizes lower than 51 SIAs, which is a size not accessible to direct TEM observations. Our model also supports the theory of the formation for 100 loops proposed by Zhang et al. (112 ). In this zero K discrete-continuum model, we do not include thermal effects and magnetic excitations. As a consequence, our results are comparable only to low temperature experiments. The results obtained, shed some light on the absence of 100 loops in low temperature experiments, and reconcile the Zhang mechanism with the experimental evidence. However, in order to validate entirely our expectations, further analysis is required.

Finally, our work makes it possible to establish scaling laws for the formation energies of various types of clusters in various materials, which is significant for multi-scale simulations such as kinetic Monte-Carlo simulations (107 , 223 , 246 -248 ), cluster dynamics studies (249 , 250 ), or mean field approximations (251 ), where simple analytic laws are needed to model the energy of large clusters. However, to enable the use of scaling laws in multi-scale simulations, the effects of temperature must be accounted for.

Our work makes it possible to extend the zero temperature scaling laws for the formation energies of various types of clusters in various materials to finite temperature. This extension was done to preserve the possibility to be completely parametrized from first principles calculations. The present formulation of the discrete-continuum model is extended to address the formation free energies, e.g. by including the temperature dependence of elastic constants. The present finite temperature extension of the discrete continuum approach has been validated by atomistic calculations using quasi-harmonic approximation. The comparison of both models indicates that the discrete continuum model is able to predict, with striking precision and without any adjustable parameter, the free energies given directly the atomistic calculation. The accuracy of predictions of discrete continuum is remarkable, even at higher temperatures. The present validation opens up many possibilities in the investigation of the free energy landscape of mesoscopic dislocation loops. Finally, we have presented two of the many possible applications. As the results obtained from the zero K discrete-continuum model relate the low temperature experiments in Fe with Zhang mechanism, the finite temperature extension of this model supports the same mechanism for explanation of high temperature experiments in Fe. The experimental observation of 100 dislocation loops at high temperatures is justified using the finite-temperature extension of the discrete-continuum model and the Zhang mechanism. As a second application, free energy values predicted by the finite temperature discrete continuum model are used as input in cluster dynamics simulations in order to compare the simulation results with time evolution of loops density in post irradiation experiments. The comparison with our model shows better results than the commonly-used capillary law.

The key ingredient useful for the extension of the discrete-continuum model at zero K to finite temperatures is the temperature dependence of the core-traction free energy. For the first time, we give the ab initio-based value of the energy of the core-traction for various dislocations at zero K in this thesis. At the present, none of the tested empirical potentials can reproduce the DFT results. The future development of empirical potentials should include this knowledge. At finite temperatures, the results and predictions can be significantly influenced by this term and further studies are needed to better understand the temperature-dependent core-traction energy.

Validation of the 0 K discrete-continuum model and its extension for finite temperatures coupled with successful applications of both the versions of the model show promise for more possible implementations. Our findings can influence even the experimental field. Because of the precision of our energetic model, it can be used as a basic tool by experimentalists in order to reinterprete some experimental data. As an example, we have shown in this thesis how the experimental determination of density of loops in post-irradiation isothermal experiments can give access even to determination of vacancy activation energy for any given temperature when coupled with cluster dynamics and our energetic model. Moreover, the cluster expansion form of the discrete-continuum model permits further extension to bcc metals with impurities, such as Fe-Cr, Fe-Mn which could be useful for radiation damage studies related to nuclear fission reactors and W-H, W-He for future fusion reactors.

A1.1 Eshelby's inclusion and Eshelby's spherical inhomogeneity transformation strain is assumed e * ij = δ ij . This seems reasonable because it is usually so in diffusional solid-state transformations. Kroner(238 ) provided approximate solutions to calculate the strain energy of an incoherent ellipsoidal precipitate in both isotropic and anisotropic matrices as a function of the ellipsoidal aspect ratio but Kroner's approximation was shown to have significant errors when used for anisotropic matrix, as explained in (237 ). The system of spherical precipitates in anisotropic matrix dealt by Lee matches well with our scenario of spherical C15 clusters in anisotropic Fe matrix with one exception that C15 clusters are coherent in BCC Fe matrix while precipitates considered in the paper are incoherent. This difference is reconciled by a paper by Christian(239 ) which establishes that there is no distinction between coherent and incoherent precipitates if only a pure dilatational transformation strain is assumed. In light of our earlier assumption of purely dilatational transformation strain, Lee's system fully corresponds with our system of interest.

In order to use the results from (237 ), we define two tensors: However, this integral cannot be solved analytically and a numerical solution using Simpson's rule is carried out which ultimately provides the total strain energy for spherical C15 clusters in anisotropic iron matrix. Eshelby's method and the above evolution is a fairly suitable way to treat the larger size limit of C15 clusters but does not work for smaller C15 clusters. However, this is fit numerically due to limitation of DFT calculations and the next chapter will describe the method using the feasible DFT calculations. The values of formation and migration energies and entropies of the mono-vacancy are fitted to reproduce self-diffusion experimental values (247 , 282 , 283 ), which only depend on the sum of the formation and migration free energies (Fig. A1.1).

A1.2 Parameters for cluster dynamics

A1 Appendix

A1.3 Sensitivity to other parameters in cluster dynamics

Sensitivity analysis is performed using model 0 as a reference (Fig. A1.2).

• Sensitivity to F f -1 + F m -1 , whatever the values of F f -1 and F m -1 (Fig. A1.3).

• No role of mobility of di-, tri-, tetra-and penta-vacancy (Fig. A1.4)

• No role of mobility of di-, tri-and tetra-interstitials (Fig. A1.5).

• No role of mobility of mono-interstitial (Fig. A1.6).

Contents 5 . 6

 56 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Conclusions and perspectives A1 Appendix A1.1 Eshelby's inclusion and Eshelby's spherical inhomogeneity . . . . . . . . A1.2 Parameters for cluster dynamics . . . . . . . . . . . . . . . . . . . . . . . A1.3 Sensitivity to other parameters in cluster dynamics . . . . . . . . . . . .

1. 1 Figure 1 . 1 :

 111 Figure 1.1: Schematic representation of some possible atom displacements in a material with face centered cubic structure resulting in different positions of the selfinterstitial atoms(SIA) and vacancies. Figure taken from Ref. (15 )

Figure 1 . 2 :

 12 Figure 1.2: Interstitial (green) and vacancy (red) configurations in 20 keV displacement cascade are depicted at (a) peak damage state at 0.3 ps and (b) after 20 ps. Evolution timescale of defect concentration in typical displacement cascade in shown in (c). Figures (a) and (b) are taken from Ref. (23 ) and figure (c) from Ref. (20 ).

1. 1

 1 Figure 1.3: Various kinds of defect reactions are summarized. Classification is based on Ref. (22 ). 23

Figure 1 . 4 :

 14 Figure 1.4: (a,b)This graph shows formation energies of different interstitial configurations relative to 111 interstitial dumbbell for various bcc transition metals.(31 , 38 ) (c)Various SIA configurations are shown(left to right): 111 , 110 , SIA at tetrahedral position, 100 , SIA at octahedral position. Figure taken from Refs.(31 , 38 ).

Figure 1 . 5 :Figure 1 . 6 :

 1516 Figure 1.5: The structure of small self-interstitial clusters in bcc iron. The lowest energy structures for clusters made of parallel dumbbells obtained within DFT are represented for clusters with 2 to 8 SIAs, with a transition from 110 to 111 orientations between 4 and 5 SIAs. The lowest energy structures within DFT known for the di-and tri-interstititals, namely the triangular and ring configurations, are also shown, as well as the related configuration for quadriinterstitial. SIAs are represented by pink or light blue spheres and vacancies by dark blue cubes.

Figure 1 . 7 :

 17 Figure 1.7: Structure of small C15 interstitial clusters in a bcc lattice. a-c: Representation by vacancies (blue cubes) and interstitials (orange spheres) of the di-, tri-and quadri-interstitial clusters. For the di-interstitial, the atoms of the bcc lattice at the center and at the edges of the Z16 Frank-Kasper polyhedron are also represented. d-f: same as a-c in a skeleton representation, i.e.without the vacancies and the cubic lattice. g, h: skeleton representation for the hexa-and octa interstitials. i: unit cell of the C15 Laves structure. For each cluster size, the configuration which is represented corresponds to the lowest energy one found within DFT.

  Figure 1.8: Fraction of 1/2 111 and 100 loops in Fe as a function of temperature. This graph shows the temperature dependent change of microstructure in iron when irradiated with 150 keV Fe ions. Figure taken from Ref.(28 ).

Figure 1 . 9 :

 19 Figure 1.9: (a)-(d) Schematical representation of the experimentally observed absorption of a smaller loop by a larger loop with a change in Burgers vector of the smaller loop while (e) Plot showing the diameter of absorbing and absorbed loops(in nm). This graph illustrates that absorbed loop size is always found to be smaller than the absorbing loop size. Figure taken from Ref.(104 )

Figure 1 . 10 :

 110 Figure 1.10: As a result of irradiation in JMTR (Japan Materials Testing Reactor of JAEA-Oarai), JOYO (Sodium cooled High Flux Reactor of JAEA-Oarai) and HFIR (High Flux Isotope Reactor of Oak Ridge National Laboratory of the USA), the observed microstructure are plotted in the graph as a function of irradiation temperature and displacement damage, as reported in the Refs. (115 , 117 -119 ). The irradiation conditions for observation of loops and voids or only voids appeared to be quite well defined based on these experiments. Figure reproduced from Ref. (120 ).

Figure 1 . 12 :

 112 Figure 1.12: Microstructure after irradiation and annealing of Vanadium specimen for 1 hr at 873 K (600 • C), showing 1/2 111 loops of vacancy type. Figure taken from Ref. (151 ).

Figure 1 . 13 :

 113 Figure 1.13: Defect density evolution in Vanadium, irradiated using 3 MeV Cu 3+ ions, is studied at various temperature variations a) shows the variation in density of three types of defects(loops, precipitates and voids) as a function of temperature b) shows number density of interstitial dislocation loops at various temperatures and irradiation conditions. Figure taken from Ref. (157 )

  Figure 1.14: Microstructural evolution in Tantalum disks after Cu ion irradiation with increasing dose at 973 K. Dose is marked above each image. Figure taken from Ref. (176 )

2 or I triangle 2 (

 22 being first reported byGao et al. in (215 )), corresponding to three atoms in interstitial positions and one vacancy sharing the same bcc lattice site (see the Figure1.5) is even more stable than the traditional I 110 2

2 2 inFigure 3 . 1 :

 231 Figure 3.1: Formation energies of the C15 SIA clusters in bcc Fe calculated with respect to the lowest energy parallel-dumbbell configurations i.e. with a 110 orientation upto 4 SIAs and a 1/2 111 orientation at larger sizes. The DFT results are compared to those of selected empirical potentials.

Figure 3 . 2 :

 32 Figure 3.2: Comparison of Peierls barrier for selected potentials and SIESTA DFT-GGA.

Figure 3 . 3 :

 33 Figure 3.3: Comparison of Peierls barrier for selected existing potentials, SIESTA DFT-GGA and the newly-developed empirical potentials.

Figure 3 . 4 :

 34 Figure 3.4: Formation energies of the C15 SIA clusters in bcc Fe calculated with respect to the lowest energy parallel-dumbbell configurations i.e. with a 110 orientation upto 4 SIAs and a 1/2 111 orientation at larger sizes. The DFT results are compared to those of selected empirical potentials and the newlydeveloped empirical potentials.

Figure 3 . 5 :

 35 Figure 3.5: (a-d) Top: Structure of small C15 interstitial clusters in a bcc lattice of the di-, tetra-, hexa-and octo-interstitial clusters, in a skeleton representation, i.e., only SIAs are represented as orange spheres without any representation of vacancies and cubic lattice sites. (a-d) bottom: centers of the Z16 Frank-Kasper polyhedron corresponding to the top C15 skeletons are represented by green spheres. (e) The 11 SIA C15 cluster, the lowest size which forms a closed ring with the centers of Z16 Frank-Kaspers polyhedra. This ring is emphasized by blue bonds connecting the centers of Z16 polyhedra.

Figure 3 . 6 :

 36 Figure 3.6: (a) The spheres in green represent the possible centers of the Z16 polyhedra. We can easily recognize the diamond structure. The plotted cube is the cubic unit cell of the original bcc lattice. (b) The unit cell of the cubic phase of Laves C15 or M gCu 2 , with the Mg atoms in green and the Cu atoms in orange. In our convention the green atoms correspond to the atoms which are the center of Z16 polyhedra and the orange atoms are the atoms which are in the interstitial positions.

97 3

 97 three closed loops are formed and packed in the most compact way possible. I C15 17 is presented in the Fig. 3.8. Development of Embedded-Atom Method potentials for defects in iron

Figure 3 . 7 :

 37 Figure 3.7: 4 configurations of I C15 10 clusters. The green and yellow spheres are the possible centers of Z16 polyhedra. Green corresponds of one orientation of the Z16 polyhedra and yellow to the orientation rotated at 90 degrees along Z axis. The blue bonds link the centers of trial configurations. The configuration (a) has the lowest formation energy for all three potentials M10, EAM2 and EAM3. The formation energy of (b), (c) and (d) configurations is higher for the potential (i) M10 with 1.15 eV, 1.41 eV and 1.70 eV, respectively (ii) EAM2 with 0.77 eV, 1.05 eV and 1.21 eV, respectively and (iii) EAM3 with 0.79 eV, 1.02 eV and 1.15 eV, respectively.

3. 4

 4 Construction of the stable C15 clusters

Figure 3 . 8 :

 38 Figure 3.8: The configuration of I C15 11 and I C15 17 which are formed by one and three closed loops of Z16 centers, respectively. The I C15 18 cluster contains 2 plain closed loops yielding a planar shape. For representing the C15 clusters the same convention is used as in .Fig3.7.

Figure 3 . 9 :

 39 Figure 3.9: Formation energy of the 1/2 111 loops and C15 clusters, calculated for selected potentials, are plotted as a function of cluster size. For each potential, the yellow empty squares and line represent the formation energies of C15 clusters while the purple circles and line represent the formation energies of the 1/2 111 loops. The crossover between these two families of defects is seen to vary considerably depending on the potential.

Figure 4 . 1 :

 41 Figure 4.1: Surfaces considered for dislocation elastic energy calculations are shown: S δ is the curved surface surrounding the dislocation, S R is the curved surface encompassing the body and S + and S -are the straight lines connecting the two curved surfaces to remove the singularity due to the dislocation. (81 , 224 ) Here, X represents a straight dislocation.

. 14 )Figure 4 . 2 :

 1442 Figure 4.2: Variation of the ratio of pre-logarithmic term in Eq. 4.11 of two different Burgers vector dislocation is shown as a function of the rotation in terms of angle of vector, l.Figure is reproduced from (27 ).

  Figure 4.2: Variation of the ratio of pre-logarithmic term in Eq. 4.11 of two different Burgers vector dislocation is shown as a function of the rotation in terms of angle of vector, l.Figure is reproduced from (27 ).

. 16 )Figure 4 . 3 :

 1643 Figure 4.3: Ratio of k b = Kb 2 of 100 dislocation loops to k b value of 1/2 111 dislocation loops is plotted as a function of temperature. This ratio provides information regarding the relative stability of the compared dislocation loops.On the basis of this graph, the 1/2 111 dislocation loops are more stable than 100 for bcc metals, except for iron at higher temperatures (greater than 830 K). The k b -values are plotted up to melting point for the bcc metals, except iron where the plot is limited till its Curie temperature.

Figure 4 . 4 :

 44 Figure 4.4: (a)Formation energies of the 1/2 111 and 100 dislocation loops in Fe against the number of SIAs, computed: (i) from atomic scale simulations using an EAM inter-atomic potential (99 ); (ii) using predictions based on the law P 0 +P 1 n 2/3 (see Ref. (76 )); (iii) using the anisotropic elasticity theory (see Eq. 4.16) which is parametrized with atomic scale calculations up to n=21 SIAs. P 0 and P 1 are fitted while term T is computed directly from the elastic tensor associated to the EAM potential. (b) The ratio of the convex hull perimeter of the dislocation loops to the perimeter deduced from the discrete number n of SIAs contained in the cluster. The full line curve was fitted using the function (1 -(1/(x a 1 + a 2 ))), with values of 0.70 and 0.88 for the exponent a 1 of 1/2 111 and 100 loops, respectively.

Figure 4 . 5 :

 45 Figure 4.5: Structure of a 100 loop in the {001} habit plane containing 15 SIAs showing the number of 1 st and 2 nd nearest neighbors of all dumbbells. Note: The loops were constructed such that each dumbbell has at least one dumbbell in 1 st nearest neighbor position. As a result, possible number of 1 st nearest neighbors varies from 1 to 4 while the number of 2 nd nearest neighbors varies from 0 to 4. There exists just one exception which is the case of mono-SIA.

Figure 4 . 6 .

 46 Database configurations were limited to 20 SIAs in 50 configurations and 22 SIAs in 31 configurations for 100 and 1/2 111 , respectively. All the training database configurations for 100 and 1/2 111 clusters are shown in Figure 4.7 and Figure 4.8.

Figure 4 . 6 :

 46 Figure 4.6: Histograms showing number of occurrence of pairs of dumbbells with respect to the type (n 1 , n 2 ) for all (a) training (small clusters) and (b) validation (small up to large clusters) 100 configurations where n 1 is the number of first nearest neighbors and n 2 refers to the number of second nearest neighbors. Possible (n 1 , n 2 ) dumbbell pairs of the type (0, n) where n = 0 to 4 were not included due to absence of such pairs in both training and validation configuration sets.

Figure 4 . 7 :

 47 Figure 4.7: Database configurations of 100 loop type. The cluster dumbbells are projected (and represented) in the 001 habit plane. The color of each dumbbell is assigned according to the number of first and second neighbors, (n 1 , n 2 ). The color assignment map is shown in Figure 4.5 of the paper.

Figure 4 . 8 :

 48 Figure 4.8: Database configurations of 1/2 111 loop type in 110 habit plane, projected in the 110 plane. The same color convenction is applied as in the Figure 4.7.

Figure 4 . 9 :

 49 Figure 4.9: Formation energies of n SIAs of (a)1/2 111 and (b) 100 dislocation loops in Fe, computed with an EAM potential (99 ) with different configurations in cubic simulation cells containing 250 + n (red circle), 432 +n (blue square), 686 + n (orange triangle up) and 1024+n (green rhombus) atoms.All energies are normalized to the asymptotic limit, taken as the formation energy in a simulation cell containing 207646 + n atoms. Lines are guides for the eyes, obtained using a fit to a 4 th order polynomial. Formation energies have been corrected using the elastic dipole correction method to account for the finite size of the simulation cells (230 ).

Figure 4 . 10 :

 410 Figure 4.10: Cluster formation energies as a function of cluster size for (a) 1/2 111 and (b) 100 loops and for (c) C15 SIA clusters in Fe calculated using the Ackland-Mendelev potential for Fe (99 ). Open diamonds represent the direct EAM results derived from simulations using large cells while the blue full diamonds represent values predicted by the discrete-continuum model.The relative errors are plotted as insets. Note that for the nanometric clusters the relative error is less than 3%.

4. 5

 5 Calculation of formation energy for interstitial defects at 0K

Figure 4 .

 4 Figure 4.11: (a) DFT formation energies of 1/2 111 , 100 and C15 clusters in Fe (empty circles, squares and diamonds, respectively) and the DFT based predictions made using the discrete-continuum model (full circles, squares and diamonds, respectively) (b) Extrapolation of the formation energies at large sizes for the 1/2 111 loops, 100 loops and C15 clusters in Fe -empty symbols. Full lines represent the elastic model (Eq. 4.16) parameterized using the points predicted by the present discrete-continuum model. This extrapolation can be done without size limitation. Note the crossover between 1/2 111 loops and the C15 clusters at 51 SIAs, and between 100 loops and C15 clusters at 91 SIAs.

Figure 4 .

 4 Figure 4.12: (a) DFT formation energies of 1/2 111 , 100 and C15 clusters in W (empty circles, squares and diamonds, respectively) and the DFT based extrapolation from discrete-continuum model (coloured lines). (b) Extrapolation of formation energies at large sizes for the 1/2 111 loops, 100 loops and C15 clusters in W -empty symbols. Full lines represent the elastic model (Eq. 5) parameterized on the points predicted by the present discrete-continuum model.

Figure 4 .

 4 Figure 4.13: (a) DFT formation energies of 1/2 111 , 100 and C15 clusters in Ta (empty circles, squares and diamonds, respectively) and the DFT based extrapolation from discrete-continuum model (coloured lines). (b) Extrapolation of formation energies at large sizes for the 1/2 111 loops, 100 loops and C15 clusters in Ta -empty symbols. Full lines represent the elastic model (Eq. 5) parameterized on the points predicted by the present discrete-continuum model.

Figure 4 .

 4 Figure 4.14: DFT and discrete-continuum model predictions for the formation energies of 1/2 111 interstitial loops and C15 SIA clusters in V. The same conventions as in Fig.4.11 and Fig.4.12 are used.

Figure 5 . 1 :

 51 Figure 5.1: The formation free energy (eV) of 1/2 111 dislocation loops at 800 K as a function of size of the cluster in the three free energy models, 1, 2 and 3, mentioned in the text through the Eqs.5.14, 5.15 and 5.16, respectively.

. 21 )Figure 5 . 2 :

 2152 Figure 5.2: Various adiabatic elastic constants i.e. Bulk modulus B, shear modulus C', elastic constants C 44 and anisotropy A are plotted as a function of temperature for Fe, W, Ta and V. Different references were used for each element to obtain the variance of the elastic constants in the temperature range of interest: Fe (254 -256 ), W (257 -259 ), Ta (258 , 260 , 261 ) and V (262 -264 ). Here, empty symbols denote the experimental values while the lines represent an extrapolation. In graphs where experimental values are missing e.g. B for Ta and V, experimental bulk modulus values were not provided and the extrapolation was obtained by deriving B from the other experimental elastic constants. Unit of elastic constants is GPa and temperature is in Kelvin.

Figure 5 . 3 :

 53 Figure 5.3: Volume, volumetric expansion coefficient α (alpha in graph), C p and the ratio B S /B T are plotted as a function of temperature for Fe, W, Ta and V. Experimental values showing temperature dependence of volume and α are taken from Ref. (268 ) while those for specific heat are taken from Ref. (269 ).Here, empty symbols denote the experimental values while the lines represent an extrapolation. Units: Volume is in Å3 , α in 10 -6 K -1 , specific heat in Jmol -1 K -1 and temperature in Kelvin.

Figure

  Figure 5.4: The values of lattice parameter (a 0 ), isothermal bulk modulus (B) and the free energy (F ) of the bcc phase of iron ) from QHA approximation implemented in this thesis and force brute calculation (as it is explained in the text). The solid line and circles emphasize the QHA and brute force calculations, respectively. The phonon gas is quantized (blue) or treated in classical approximation (orange).

Figure 5 . 5 :

 55 Figure 5.5: The values from QHA approximation of lattice parameter (a 0 ), C P,V of the bcc (left) and C15 (right) phase of iron using M10 potential (100 ). The phonon gas is quantized (blue) or treated in classical approximation (orange).

Figure 5 . 6 :

 56 Figure 5.6: The values from QHA approximation of the isothermal elastic constant of bcc phace using two empirical potentials for iron (99 ) (left) and W (135 ) (right). The phonon gas is quantized (blue) or treated in classical approximation (orange).

Figure 5 . 7 :

 57 Figure 5.7: Formation free energy of dislocation loops of 100 and 1/2 111 with various sizes using atomistic QHA calculations and the three free energy models derived from discrete continuum model. All models are based on the interatomic potential Ackland-Mendelev (99 ). Inset is the shape of the loop corresponding to each graph (the same convention is applied as in Figs 4.6 and 4.7).

Figure 5 . 8 :

 58 Figure 5.8: Formation free energy of C15 containing 36 SIAs given by atomistic QHAcalculations and the finite temperature discrete continuum model (using isotropic and anisotropic elasticity). All models are based on the interatomic potential. The cluster is drawn using the centers of Z16 polyhedron lying on the bcc matrix which are linked with blue bonds (the same convention is applied as in Fig.3.7)

Figure 5 . 9 :

 59 Figure 5.9: The lattice parameter, adiabatic bulk modulus, C' and C 44 adiabatic elastic constants of α-iron. The experimental results are emphasized by full triangle while the DFT using QHA approximation results by black circles or full black lines. The circles denote the temperatures chosen for QHA approximation, the black line being an interpolation of those points. The orange curves denote the DFT results rescaled to the experimental value of the corresponding observable at zero K.

Figure 5 . 10 :

 510 Figure 5.10: The DFT values of lattice parameter, isothermal bulk modulus, C' and C 44isothermal elastic constants of the ferromagnetic C15 phase of iron. The circles denote the temperatures chosen for QHA approximation, the black line being a interpolation of those points.

Figure 5 . 11 :

 511 Figure 5.11: The crossover (expressed in number of interstitials) in the relative stability of 100 (blue) and 1/2 111 (orange) loops with respect to C15 clusters in α-iron using finite temperature discrete continuum approach. The shaded area denotes the estimated errors. The light blue shaded regions represents the error arising due to temperature dependence of elastic constants and the inherent error in the zero K discrete-continuum model. Using the best parametrization for zero K, the darker shaded region around both lines denotes the error arising from only the use of finite-temperature discretecontinuum method.

Figure 5 . 13 :

 513 Figure 5.13: Initial loop distribution (after annealing during 50 min).

Figure 5 . 14 :

 514 Figure 5.14: Schematic representation of the Ostwald ripening of interstitial loops by vacancy emission. The thickness of the arrows is proportional to the vacancy flux

Figure 5 . 15 :

 515 Figure 5.15: Binding energies of vacancies with interstitial loops

Figure 5 . 16 :

 516 Figure 5.16: Comparison of loop distributions obtained with different loop energy models at (a) t = 72 min and (b) t = 100 min.

Figure 5 . 17 :

 517 Figure 5.17: Loop distributions obtained at different times with a capillary law, for a migration free energy of the vacancy equal to (a) 0.51 eV and (b) 0.45 eV.

Figure 5 .

 5 Figure 5.18: Self-diffusion coefficient in α-iron. The green dot represents the value obtained if the properties of the vacancy are fitted to reproduce experimental loop measurements with a capillary law. Yellow and blue dots are from Ref. (282 ) and Ref. (283 ), respectively

1 .

 1 Auxiliary Tensor, D (286 ): It is defined such thatu c i,l (x) = -σ * kj D ijkl (x) ,(A1.7)where u c i,l is the constrained displacement gradient and σ * kj is the eigenstress inside the inclusion. Within a spherical inclusion, D is shown to be: ij z k z l sin φdΘdΦ.(A1.8)Using the definition of Eshelby tensor, it can also be shown that:S ijkl (x) = -1 2C lkmn (D iklj (x) + D jkli (x)) .

Figure

  Figure A1.1: Self-diffusion coefficient in α-iron. The black line corresponds to the values of the formation and migration free energies adopted. The red dot indicates the temperature at which the simulations are performed (T = 833 K).

  Figure A1.2: Reference calculation for sensitivity analysis, model 0

  Figure A1.4: Calculation without the mobility of small vacancy clusters

Figure A1. 5 :

 5 Figure A1.5: Calculation without the mobility of small interstitial clusters

  • C) but not if neutron-irradiated iron at 333 K (60 • C) is subjected to post-annealing heating up to 773 K (500 • C). Based on this, Eyre and Bullough ruled out the possibility of formation of 100 loops by combination of two closely spaced 1/2 111 loops. In addition, they acknowledged presence of square 100 loops and circular 1/2 111 .

Table 4 .

 4 1: The anisotropy factor (A) and the elastic constants, bulk modulus (B), shear modulus (C'), C 11 , C 12 and C 44 (in GPa) of different bcc metals from experiments and DFT calculations.

	4 The Discrete-Continuum model
		Fe		W		V		Ta
		Expt. DFT Expt. DFT Expt. DFT Expt. DFT
	B	173.0 150.7 314.2 303.2 160.7 184.2 194.2 193.1
	C' 53.0	63.0	163.8 162.3 58.1	74.2	54.1	40.6
	C 11 243.7 234.6 532.6 519.6 238.2 283.3 266.3 247.4
	C 12 137.7 108.7 204.0 190.0 122.1 134.8 158.2 166.1
	C 44 122.0 116.0 163.1 141.2 38.3	32.1	87.4	68.3
	A	2.30	1.84	0.99	0.87	0.81	0.43	1.61	1.68

Table 4 .

 4 4.2. For any subsequent use, it is recommended to compute the formation energies using best-fit parameters from Table4.2 in Eqs. 4.25 and 4.26 for size of clusters higher than n = 15. For lower values of n the formation energies of various configurations are provided in the Tab. 4.3. The presented configurations are not necessarily those with the lowest formation energy. Moreover, it should be noted that in the case of iron, as is mentioned earlier in Section 4.3, the lowest energy configurations of small parallel clusters are 110 and not 111 or 100 . 2: Best fit parameters (Eq.4.25 and Eq. 4.26) for the formation energies extrapolated using the discrete-continuum model for the three types of clusters in three different bcc crystals. a 0 , a 1 and a 2 are expressed in eV.

			4.5 Calculation of formation energy for interstitial defects at 0K
	Element		a 0			a 1				a 2
		1/2 111	100	C15	1/2 111	100	C15		1/2 111	100	C15
	Fe	1.6049	1.7768 0.45330	5.3523	7.1595 1.3547 -0.14732 -5.8180 9.3542
	W	3.9290	4.8488 1.0967	7.9242	13.6984 3.3295	6.2009	-8.2585 44.6802
	V Ta	1.0303 2.2649	-2.3573 3.3221 1.5365	0.8541 1.0118	-4.4935 2.1860 10.7638 -0.1153 -20.7005 1.3499 3.8621 -7.3536 -
		Fe			W		V			Ta
	size C15 1/2 111	100	C15 1/2 111	100	C15 1/2 111	C15	1/2 111	100
		4.90	5.31		10.48	12.96		2.78	4.740
	8.03	8.63	8.33 23.71	18.40	21.31 7.99	4.80	11.609	8.439
	11.31	11.54	12.73 27.33	25.01	28.76 9.78	6.86	12.484
	11.28	14.48	14.32 34.35	30.79	34.81 10.75	8.00	14.870 14.858 16.889
	14.39	17.28	17.22 42.11	36.39	41.23 13.22	9.63	18.709 17.959 21.533
	15.85	19.80	19.70 45.34	41.52	46.59 13.71	10.86	19.469 20.480
	17.96	21.65	21.79 53.04	45.25	53.33 16.15	11.34	23.280 21.940
	18.76	24.22	23.74 56.29	50.44	57.26 16.59	12.87	24.094 24.696 28.447
	20.77	26.61	27.01 62.20	55.09	62.86 18.37	13.97	26.804 27.011 31.159
	21.70	28.29	28.96 65.96	58.49	67.89 19.04	14.49	27.848 28.362 34.411
	22.71	30.27	30.64 69.88	62.95	72.84 19.84	15.69	29.115 30.587 36.628
	25.55	31.77	32.70 76.84	66.23	76.83 22.04	16.20	32.596 31.858 39.490
	26.67	35.08	35.85 80.82	73.85	84.80 22.84	18.85	33.803 36.271 42.359
	28.55	35.09	38.37 86.13	73.64	89.21 24.16	17.84	36.071 35.376 44.177
	29.18	37.28	39.46 90.12	77.81	92.93 25.74	18.97	37.608 45.607

Table 4 .

 4 3: DFT formation energies (expressed in eV) of smallest interstitial clusters, up to 15 interstitials, as a function of size. In this table, only the lowest formation energies are reported, at the same number of self-interstitial atoms, for the configurations used in the present database.

	4 The Discrete-Continuum model

  zero K results are comparable only to low temperature experiments.

	2.87
	2.86
	2.85
	2.84	0 300 600 900

f is the difference between the formation energy of the diinterstitial triangle and 110 configurations (positive/negative values indicate that 110 /triangle is the most stable.)

Finite-temperature extension of discrete-continuum model

Acknowledgements

1 State of the Art Review observations of 100 loops in Fe, thus reconciling experiments with DFT calculations of energy landscape of point defects in bcc Fe. The sections 4.5.1 and 5.4 will provide a detailed discussion of the existing models and the present findings of this thesis.

Summary for tungsten: While recent experimental observations reveal presence of 1/2 111 and 100 interstitial(mainly) loops at higher temperatures (>573 K), previous experiments observed 110 and 1/2 111 vacancy loops at room temperature. Empirical potentials seem to be inconclusive with some predicting 1/2 111 as most stable and others show a crossover between 1/2 111 and 100 loops suggesting competing stability of these loops as a function of number of interstitials. As such, MD simulations using different potentials can be contradictory as well but DFT calculations have quite definitively established the stability of 1/2 111 interstitial loops over other possibilities.

Summary for vanadium: In general, dislocation loops of Burgers vector 1/2 111 , voids and planar precipitates are observed after neutron irradiation. However, there are contradictions among the experiments with one reporting the presence of 100 interstitial loops in addition to the above-mentioned defects and another stating the nature of 1/2 111 dislocation loops as vacancy-type. These contradictions have been attributed to the varying levels of impurities among the different experimental samples. Ab initio calculations have established that the 111 interstitial dumbbell is the most stable and even predictions using interatomic potentials concur regarding the stability of 1/2 111 interstitial loops. In fact, simulations using Finnis-Sinclair potential for V predicts easy rotation of 100 and 110 loops into 1/2 111 due to a significant difference in their formation energies.

Summary for tantalum:Scarcity of experiments and contrasting irradiation conditions among the various experiments prevent any conclusive evidence of finding a particular type of defect over another. While one experiment found screw dislocations with Burgers vector 1/2 111 , another reports growth of interstitial loops which is followed by a preferential formation of vacancy loops at higher temperatures. Nevertheless, it has already been established from first principles that 1/2 111 loops are the most stable.

Continuum Methods for 2D dislocation loops

• Section 4.4: In this section, validity of the discrete-continuum model is established using existing empirical potentials for iron and tungsten.

• Section 4.5: Following validation of discrete continuum model, we proceed to its parameterization using the DFT formation energies of configurations included in our database. With this ab initio parametrization the formation energies of interstitial defects without size-limitation at 0 K are provided. The subsections include the following:

-Section 4.5.1: Ab initio based predictions of formation energy landscape of SIA clusters are presented and discussed.

-Section 4.5.2: A useful application of the discrete-continuum model are the ab initio scaling laws for formation energy calculations. These scaling laws are supplemented by the formation energy values corresponding to various defect-types are given and should be utilised when number of SIAs is less than 15 for each bcc metal.

• Section 4.6: Finally, the main points of this chapter are recapitulated systematically to provide a clear picture of the energetic model developed.

Continuum Methods for 2D dislocation loops

Continuum methods involve treatment of materials as continuous mass by assuming the constituents to be infinitesimal elements with identical properties as the bulk. This treatment is fairly valid at macroscopic levels and is commonly used for elasticity studies.

Here, we deal with the problem of calculation of total elastic formation energy of a 2D dislocation loop and a 3D cluster in the bcc metal medium. It is noteworthy that most of the bcc metals (except W) considered in this thesis work, i.e. Fe, V and Ta, are anisotropic (See Tab. 4.1).

The 2D dislocation loops are treated using the full anisotropic elasticity approximation, as explained in following paragraphs. According to the elastic theory of dislocations (103 ), the formation energy of a 2D plate-like clusters of SIAs is related essentially to two quantities: the line energy density of the edge dislocation which encloses the clus-4.3 Construction of database C15 cluster as: V C15 = 2nΩ. The convex hull surface and the atomic interface energies, i.e., the first and last terms of Eq.4.23, have been combined into a single term that is the first term of Eq.4.24. The prefactor γ s in the first term, which gives the dependence in n 2/3 , plays the role of interface energy. Because of the fact that S C15 and N i terms of the discrete formulation in Eq. 4.24 take the form n 2/3 in the infinite limit, the new interface energy γ s combines γ and ∆E i .

The elastic energy of C15 clusters can more accurately be considered by treating the second term of Eq.4.23, as Eshelby inhomogeneities in anisotropic bcc matrix and the elastic energy of these C15 clusters be treated accordingly, as given in Appendix A1.1. As explained in the appendix, this treatment considers C15 clusters as spherical inhomogeneities in an anisotropic cubic matrix in accordance with . The results obtained by this alternative treatment have not been presented because they don't vary much qualitatively from the current results.

Construction of database

The energy landscape of dislocation loops has been widely studied at the atomic scale by various authors (27 , 100 , 110 , 133 , 240 ) and the lowest energy configurations of various clusters have been already reported in the literature. 1/2 111 dislocation loops were generated by inserting 111 dumbbells in the {110} habit plane so as to form compact clusters. 100 dislocation loops were generated by inserting atomic 100 dumbbells in the {001} habit planes. In addition to these configurations, we have found that some specific cluster geometries needed to be included as 'database' in order to have a good parameterization for E(n 1 , n 2 ). As explained later in this section, the choice of these configurations is made in such a way that the neighborhood of large loops is reproduced, i.e., those containing thousands of SIAs. Here, the size of clusters included in the database is limited by the feasibility of ab initio calculations. The size of clusters in the database ranges from 2 to 20 SIAs with a total of about 50 configurations for 100 and from 2 to 22 SIAs with 31 configurations for 1/2 111 . The choice of the largest cluster size, i.e. 22 SIAs for 1/2 111 , is justified by the low relative error in the formation energy in simulation cells with 1024 + n atoms (8a × 8a × 8a cells). As shown in Fig. 4.9, the relative errors in the formation energies due to this size limit are lower than 2%. More details, along with the exhaustive list of configurations included in the 5 Finite-temperature extension of discrete-continuum model

For a system at volume V , the free energy of N atoms can be expressed as:

where F el is the Mermin free energy of electrons (252 , 253 ) which takes into account the electronic internal energy and the entropy of electrons at the temperature T . The vibration free energy F vib contains the internal energy of phonons and the vibrational entropy of ions. The above entropic quantities, electronic S el and well as vibrational S vib are difficult to determine. The general tendency in the community is to neglect the comparatively smaller electronic contribution and replace the vibrational entropy effects by empirical laws. The main large scale simulations which fill the gap between theory and experiments, such as the kinetic Monte Carlo, cluster dynamics, dislocation dynamics simulations, use these empirical laws to calculate the contribution from the vibrational part.

Having established the physical significance and relevance of calculating free energy, calculation of free energy is illustrated in Sec. 5.2 by the finite temperature extension of discrete continuum model. Furthermore, in the Sec. 5.3.1 is presented the quasi harmonic approximation which provide us the atomistic free energy results for loops or C15 clusters used to validate the finite temperature discrete continuum model.

Finite Temperature Extension of

Discrete-Continuum model

Free energy calculations using discrete-continuum model As explained in Sec.4.1.1 of Ch.4, the formation energy for 2D loops at 0K using the discrete-continuum model is:

5 Finite-temperature extension of discrete-continuum model for a system with Hamilton H( r 1 , r 2 , ..., r N ; p 1 , p 2 , ..., p N ) can be written as:

(5.37)

Once the canonical partition function is calculated for a system, all thermodynamic observables can be derived. For example, the Helmholtz free energy is found to be:

(5.38)

Differentiating the Helmholtz free energy, thermodynamic properties can also be calculated as follows:

)

Thus, calculation of the canonical partition function of the system by classical treatment provides a convenient means to obtain the required thermodynamic information, provided conditions of validity hold (271 ). Unfortunately the integration over the highly dimensional phase space in order to deduce the partition function Z is possible in only few simple cases. One of them, the harmonic approximation, will be discussed further ahead in this section. Formation free energy landscape of defects requires calculation of free energy which in turn can be derived from the partition function of the system by considering an appropriate ensemble. Free energy calculation using quasi-harmonic approximations will be discussed in Sec.5.3.1.

Quasi-Harmonic approximation

Before proceeding to the description of quasi-harmonic approximations, results from treatment of harmonic approximation are recapitulated.

Harmonic Approximation

Any arbitrary potential surface can be approximated by a harmonic potential in the close vicinity of its equilibrium position, provided the corresponding energy of the system is 5 Finite-temperature extension of discrete-continuum model

We have tested the present variant of the quasi-harmonic approximation against the brute force calculations for observables such as the free energy, lattice parameter and the bulk modulus of the bcc structure at various temperatures. This test stresses the flexibility of the polynomial form proposed in Eq. 5.69 for the free energy (V, T ) surface as well as the first derivative (lattice parameter) and the second derivative (bulk modulus) with respect to the volume. The brute force calculations are performed on a grid with temperatures {T i } ranging from 0 to 1800 K in steps of 100 K. These calculations were performed as follows:

• At a temperature T i , the corresponding lattice parameter of the bcc lattice a T i 0 are deduced. A first guess a T i guess is provided which is the lattice parameter deduced from the previous temperature T i-1 , a

(For the lowest temperature of grid we started with the lattice parameter deduced from 0 K static calculation).

• Around the guessed lattice parameter, a T i guess , we deformed the bcc solid with a volumetric strain ± 0.8 % and 7 bcc systems are built in this interval.

• The free energy for each of the 7 systems is computed.

• From the free energy-volume plot and using a fit with Birch-Murnagham-Vinet equation of state (274 -276 ) we got the value of the bcc lattice parameter a T i 0 as well as the free energy and the bulk modulus at the temperature T i .

• The procedure is repeated for the temperature T i+1 .

The comparison between the two approaches is emphasized in the figure 5.4. The agreement between the two approaches is excellent up to 10 -5 eV, 10 -4 Å and 10 -2 GPa values of the free energy, lattice parameter and bulk modulus, respectively. This agreement holds for both treatments of the free energy in Eq. 5.69, classical using Eq. 5.62 and quantum using Eq. 5.55.

It is interesting to note the differences between the classical and quantum free energy of the phonon gas. As expected, for temperatures higher than Debye temperature (347 K for Ackland-Mendelev potential for iron) both approaches, classical and quantum are identical. Below Debye temperature down to very low temperature, near 0 K, there are small differences in all three observables due to the quantum zero point energy which has no counterpart in the classical approximation. These discrepancies are not related 5.3 Atomistic models for interstitial defects at finite temperature

Validation of finite temperature discrete-continuum approach using atomistic calculations

In this section we will compare the predictions of the three free energy models proposed as extension of the discrete continuum model and the direct atomistic calculations. This comparison is made for all types of clusters 1/2 111 , 100 and C15. The shape of clusters is chosen with various sizes in order to test the Einstein approximation for the core-traction free energy, as explained in the Sec. 5.2. The size of the clusters are large, up to few hundreds of SIAs. In order to have reliable results of the free energy for this size of clusters the bcc matrix should have 128 000 atoms. Computing the frequency of phonons of these boxes is a challenging task. Consequently, this type of comparison is impossible to perform using interatomic interactions other than the EAM potentials.

We will use the Ackland-Mendelev potential ( 99), one of the potentials already used in order to validate the zero K discrete-continuum model in Sec.4.4

Firstly, using Ackland-Mendelev potential ( 99) and QHA approximation, we numerically obtain the quantities needed in order to parametrize the finite temperature discrete continuum model. The three models for loops given by Eqs. 5.14, 5.15 5.16, as well as the model proposed for C15 clusters Eq.5.19 are investigated. To parametrize the previous equations, the required quantities are a bcc (T ), C 11,12,44 (T) as well as the mono-interstitial formation entropies, S . Additional quantities to complete the input list of discrete continuum model for C15 clusters is ∆ C15,bcc (T ) = E C15 coh (T ) -E bcc coh (T ). All these calculations can be rapidly done using simulation boxes containing maximum 128 (bcc matrix) +1 (for mono-interstitial entropy) atoms. The prediction of discrete continuum models for two different sizes of loops for each type are presented in the Fig. 5.7. The prediction of free energy for C15 clusters are reported in Fig. 5.8

The atomistic calculations are done using the present implementation of QHA approximation. The number of atoms is 128000 ± n (where n is the number of SIAs in cluster). 12 volumes are generated for each system in order to accurate interpolate the free energy landscape using Eq.5.69. Each of 12 boxes is relaxed and the boxes are considered as relaxed when all the forces are smaller than 10 -4 eV/ Å. Normal modes at each volume are computed by diagonalizing the Hessian matrix. Exact diagonalization of the Hessian was performed using diagonX package, already used in (214 ). The free 5.5 Application: concentration of dislocation loops post irradiation A preliminary analysis of C15 clusters in Fe presents a scenario similar to the one solved by Eshelby using the popular Eshelby's inclusion method (235 , 236 ) which is illustrated in (284 ). However, our scenario of interest is complicated by two facts. Firstly, the elastic constants of the matrix and the inclusion are different. Technically speaking, C15 clusters are inhomogeneities. Secondly, the matrix is not isotropic as for Eshelby's case but anisotropic (Zener anisotropy factor, A = 2.4 for Fe (285 )), as mentioned in Table4.1 of Sec.4.1.

Originally, Eshelby provided elastic energy formulation only for a uniform permanent (inelastic) deformation called inclusion in a homogeneous linear elastic isotropic solid matrix. Here, the elastic constants of the inclusion and the matrix were assumed equal but this formulation can be readily extended to inhomogeneities using the equivalent inclusion method by assuming ellipsoidal inhomogeneities for simplicity. This means that in order to calculate strain energy contribution to formation energy of C15 clusters by the equivalent inclusion method, the inhomogeneity is essentially replaced by an equivalent inclusion such that the stress and strain remain the same as that inside the original inhomogeneity. These results will be presented. Having established the complexity of the problem, a further simplification of the case of C15 clusters is carried out by reasonably approximating them as spherical inhomogeneities in anisotropic matrix. This approximation and other assumptions used to treat C15 clusters will be substantiated A1 Appendix further ahead in this section.

All primed-symbols (C ijkl , e * ij , etc.) are used for the inhomogeneity i.e. C15 clusters while the same symbols without prime are for the matrix phase or equivalent inclusion, depending on the context. Let σ I ij and e c ij denote the stress and constrained strain inside an equivalent inclusion while σ I ij and e c ij denote the stress and constrained strain inside the inhomogeneity. By assumption of equivalent inclusion method for an elliptical inhomogeneity, we have the following equations:

A fourth order tensor (popularly known as Eshelby's tensor S) is defined as:

where e * ij is the eigenstrain in the equivalent inclusion i.e. strain in the absence of any external stress. Using Hooke's law, we can write stresses inside the inhomogeneity and the equivalent inclusion as follows:

where σ c ij , σ * ij and C ijkl denote the constrained stress, eigenstress (solely caused by eigentrain) and elastic constants in the inhomogeneity, respectively while σ c ij , σ * ij and C ijkl refer to those in the equivalent inclusion. e * kl is the eigenstrain in the inhomogeneity. Using the eq. A1.4 to satisfy the conditions in eq. A1.1 and A1.2, we obtain an expression to calculate eigenstrain of the equivalent inclusion as:

Once the eigenstrain for the equivalent strain is calculated, the total strain energy for the inhomogeneity in an anisotropic medium can be calculated as:

Although theoretical solutions for inhomogeneities in anisotropic matrix are well formulated, there is no general analytical solution for anisotropic matrix.

In order to numerically solve this problem, we make certain reasonable assumptions. Firstly, C15 clusters are assumed to have spherical shape since the most stable configurations of C15 clusters are closely packed 3D structures. Secondly, a pure dilatational where M ij = C ijkl Z k Z l and Z = (sin φ cos θ, sin φ sin θ, cos φ). Now, re-writing eq.A1.3 with the newly defined tensor, e c ij = T ijkl C klmn e * mn . (A1.12)

In eq.A1.6, we substitute the value for stress inside equivalent inclusion from eq.A1.4 and using eqs.A1.1-A1.2, we get the following form to calculate the strain energy:

For a spherical precipitate in cubic matrix, there exist only three independent components of T: T 1111 , T 1122 and T 1212 which simplifies the problem at hand. Summarising results using Kroner's notation(237 ), we have .14) where

In order to calculate the sum T 1111 + 2T 1122 , an integral of the following form is to be solved: .16) where Title : Energy landscape of defects in body-centered cubic metals.

Keywords :Metals, energy landscape, defects, empirical potentials, ab initio Abstract : The structural materials of nuclear reactors are subjected to severe irradiation conditions, leading to changes in their mechanical properties. In this thesis, we studied the properties of the energy landscape of point defects in bcc metals, i.e. Fe, W, V and Ta. We overcome the size-limitation of ab initio calculations and unreliability of semi-empirical potentials by development of a unique energetic model for calculation of formation energies of dislocation loops as well as for three-dimensional interstitial cluster type C15.We also expanded the scope of our energy model to the free energy of defect formation, thereby establishing predictions for different temperatures. Additionally, by using this energy model, we showed significant improvements over older energy laws, such as the capillary law widely-used in multiscale computation cluster dynamics.