Céline Varvenne

Olivier Mackain

Emmanuel Clouet
email: emmanuel.clouet@cea.fr

Vacancy clustering in zirconium: an atomic scale study ✩

Keywords: Ab initio calculations, Point defects, Stacking faults, Dislocation loops, Zirconium

The stability properties of vacancy clusters in hexagonal close-packed Zr, cavities and dislocation loops, are investigated at the atomic scale, with a modeling approach based on density functional theory and empirical potentials. Considering the vacancy-vacancy interactions and the stability of small vacancy clusters, we establish how to build the larger clusters. The study of extended vacancy clusters is then performed using continuous laws for defect energetics. Once validated with an empirical potential, these laws are parameterized with ab initio data. Our work shows that the easy formation of a loops can be explained by their thermodynamic properties.

Introduction

Zirconium alloys are widely used in the nuclear industry as a cladding material. In nuclear reactors, they are subjected to a fast neutron flux, leading to the creation of a large amount of point defects, both vacancies and self-interstitials. These point defects then diffuse and can be trapped by the different sinks of the system, or can cluster to form larger defects, like dislocation loops and cavities [START_REF] Onimus | Radiation effects in zirconium alloys[END_REF]. Vacancy clusters can also appear in quenched zirconium alloys [START_REF] Carpenter | Vacancy precipitation in zirconium alloys[END_REF].

Extensive experimental studies have been carried out in the past to determine the structure of these defect clusters in hexagonal close-packed (hcp) Zr and its alloys (see Ref. [START_REF] Onimus | Radiation effects in zirconium alloys[END_REF] for a recent review). At low irradiation doses, perfect dislocation loops with a = 1/3 11 20 Burgers vector are observed [START_REF] Northwood | Characterization of neutron irradiation damage in zirconium alloys -an international "round-robin" experiment[END_REF][START_REF] Griffiths | Neutron damage in zirconium alloys irradiated at 644 to 710 K[END_REF][START_REF] Griffiths | A review of microstructure evolution in zirconium alloys during irradiation[END_REF]. These a loops are both of interstitial and of vacancy type. Their habit plane is close to the prismatic plane of the hcp lattice. The same perfect loops, all of vacancy type, are obtained in quenched Zr alloys [START_REF] Carpenter | Vacancy precipitation in zirconium alloys[END_REF]. Under irradiation, loops lying in the basal planes are also observed for the highest irradiation doses [START_REF] Griffiths | Neutron damage in zirconium alloys irradiated at 644 to 710 K[END_REF]. These loops are faulted with a Burgers vector 1/6 20 23 , thus with a c component. They are all of vacancy type. Finally, cavities are observed in only a very few specific cases [START_REF] Griffiths | Hvem study of the effects of alloying elements and impurities on radiation damage in Zr-alloys[END_REF][START_REF] Griffiths | Study of point defect mobilities in zirconium during electron irradiation in a high-voltage electron microscope[END_REF]. These vacancy and interstitial clusters have important consequences on the macroscopic behavior of zirconium. Like in other metals, strong hardening is associated with the presence of these defects [START_REF] Onimus | Radiation effects in zirconium alloys[END_REF][START_REF] Boyne | Numerical simulation of irradiation hardening in zirconium[END_REF]. Irradiation also leads in hcp Zr to dimensional changes without any applied stress: a Zr single crystal undergoes under irradiation an elongation along the a axis of the hcp lattice and a shortening along the c axis, with no significant volume change [START_REF] Carpenter | Irradiation growth of zirconium single crystals: A review[END_REF]. The growth strain remains small at low fluence, but a breakaway growth is observed at higher fluence [START_REF] Carpenter | Irradiation growth of zirconium single crystals: A review[END_REF][START_REF] Fidleris | The irradiation creep and growth phenomena[END_REF]. This breakaway is correlated with the appearance of the c type vacancy loops [START_REF] Griffiths | Neutron damage in zirconium alloys irradiated at 644 to 710 K[END_REF][START_REF] Griffiths | The formation of c-component defects in zirconium alloys during neutron irradiation[END_REF].

Understanding the formation of these clusters is of prime importance to be able to model the kinetic evolution of the microstructure under irradiation and of the associated macroscopic behavior. This requires first to know the relative stability of these clusters, in particular vacancy clusters for which different types coexist. Atomistic simulations appear as the suitable tool for such a study, as they can provide information on cluster sizes which are not accessible by other techniques and which are necessary to build higher level models or theories. Several works already attempted to address this question [START_REF] Kapinos | Computer simulation of vacancy loops and stacking faults in zirconium[END_REF][START_REF] Kulikov | Vacancy dislocation loops in zirconium and their interaction with self-interstitial atoms[END_REF][START_REF] De Diego | Structure and properties of vacancy and interstitial clusters in α-zirconium[END_REF][START_REF] De Diego | On the structure and mobility of point defect clusters in alpha-zirconium: a comparison for two interatomic potential models[END_REF]. They either showed that the most stable vacancy clusters are cavities [START_REF] Kulikov | Vacancy dislocation loops in zirconium and their interaction with self-interstitial atoms[END_REF] or basal loops [START_REF] Kapinos | Computer simulation of vacancy loops and stacking faults in zirconium[END_REF][START_REF] De Diego | On the structure and mobility of point defect clusters in alpha-zirconium: a comparison for two interatomic potential models[END_REF]. This appears in contradiction with experimental observation indicating that the easiest clusters to form are vacancy loops lying in the prismatic planes. But all these simulations relied on empirical potentials, either long ranged pair potential [START_REF] Kapinos | Computer simulation of vacancy loops and stacking faults in zirconium[END_REF] or Embedded Atom Method (EAM) [START_REF] Kulikov | Vacancy dislocation loops in zirconium and their interaction with self-interstitial atoms[END_REF][START_REF] De Diego | Structure and properties of vacancy and interstitial clusters in α-zirconium[END_REF][START_REF] De Diego | On the structure and mobility of point defect clusters in alpha-zirconium: a comparison for two interatomic potential models[END_REF]. These central forces empirical potentials are known to poorly model stacking faults in hcp transition metals. Legrand [START_REF] Legrand | Relations entre la structure électronique et la facilité de glissement dans les métaux hexagonaux compacts[END_REF] showed that one needs to correctly account for the electronic filling of the valence d band, and thus to consider the angular dependence of the atomic bonding, in order to obtain a good description of these stacking faults. As vacancy loops, at least the smallest ones, are faulted, it is worth looking at the stability of the vacancy clusters in hcp Zr with a better modeling of the atomic bonding than with the previously used empirical potentials. Ab initio calculations represent a nice alternative but they can only be used to study small clusters containing a few vacancies. We therefore propose to use an hybrid approach relying both on ab initio calculations and empirical potentials to model these clusters.

In this work we focus on the stability of vacancy clusters in hcp Zr. Stability of small clusters, as well as stacking-fault and surface energies, are investigated with ab initio calculations and then compared to predictions obtained with two recent EAM potentials developed by Mendelev and Ackland [START_REF] Mendelev | Development of an interatomic potential for the simulation of phase transformations in zirconium[END_REF]. As a result of this comparison, one EAM potential is selected to study larger vacancy clusters. This allows us to validate analytical laws based on continuous models able to describe their formation energies. These analytical laws are finally parameterized on ab initio calculations so as to conclude on the relative stability of the different vacancy clusters in hcp Zr.

Details of atomistic simulations

Our ab initio calculations are based on the Density Functional Theory (DFT), using the Pwscf code of the Quantum Espresso package [START_REF]QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials[END_REF]. All calculations are performed in the Generalized Gradient Approximation with the exchange-correlation functional of Perdew-Burke-Ernzerhof [START_REF] Perdew | Generalized gradient approximation made simple[END_REF]. Valence electrons are described with plane waves, using a cutoff of 28 Ry. The core electrons are replaced by an ultrasoft pseudo-potential of Vanderbilt type, including 4s and 4p electrons as semicore. The electronic density of state is broadened with the Methfessel-Paxton function, with a broadening of 0.3 eV. The integration is performed on a regular grid of 14 × 14 × 8 k-points for the primitive cell and an equivalent density of k-points for the supercells used in defect calculations. This ab initio modeling approach has been already validated on Zr bulk properties in a previous study [START_REF] Clouet | Screw dislocation in zirconium: An ab initio study[END_REF].

Ab initio calculations of vacancy clusters, including the single vacancy, are performed in a periodic supercell corresponding to 5 × 5 × 4 hcp primitive unit cells and containing 200 atomic lattice sites. Only the atomic positions are relaxed while the periodicity vectors are kept fixed (constant volume calculations). The elastic correction described in Ref. [START_REF] Varvenne | Point defect modeling in materials: coupling ab initio and elasticity approaches[END_REF] is applied so as to remove the elastic interaction of the vacancy cluster with its periodic images.

The two EAM potentials we used were developed by Mendelev and Ackland [START_REF] Mendelev | Development of an interatomic potential for the simulation of phase transformations in zirconium[END_REF]. They are labeled #2 and #3 in Ref. [START_REF] Mendelev | Development of an interatomic potential for the simulation of phase transformations in zirconium[END_REF]. Both of them give a reasonable description of the bulk properties of hcp Zr. EAM #3 potential has already been used to calculate the properties of small vacancy and interstitial clusters by De Diego et al. [START_REF] De Diego | On the structure and mobility of point defect clusters in alpha-zirconium: a comparison for two interatomic potential models[END_REF]. It is believed to be better to study defect properties in hcp Zr, as some stacking fault energies in the basal and prism planes have been adjusted on ab initio values. EAM #2 potential is particularly designed to describe the hcp-bcc transition, but it also gives a reasonable description of defects in hcp Zr. We will see in the following that it is actually better suited than the EAM #3 potential to study vacancy clustering. Atomistic simulations with these empirical potentials are performed with a 100 × 100 × 50 supercell containing 1 million of atomic lattice sites. This size gives well converged energies for all the investigated defects.

Before using these different atomic models to study vacancy clustering, it is worth comparing their results for the single vacancy properties. The obtained vacancy formation and migration energies are compared with experimental data in Table 1. All three models lead to a vacancy formation energy which is compatible with the lower-bound value given by positron annihilation spectroscopy (PAS) [START_REF] Hood | Diffusion and vacancy properties of α-Zr[END_REF]. The vacancy can migrate along two non-equivalent pathways: one inside the basal plane (E mig bas ) and the other one out of the basal plane , respectively in the basal plane and along the c axis, non-null components of the elastic dipole tensor, and relaxation volumes. Energies and dipoles are given in eV and the relaxation volume is normalized by the atomic volume Ω. [START_REF] Neely | Damage rate and recovery measurements on zirconium after electron irradiation at low temperatures[END_REF] 0.58 ± 0.04 Expt. (PAS) [START_REF] Hood | The recovery of single crystal α-Zr from low temperature electron irradiation -a positron annihilation spectroscopy study[END_REF][START_REF] Hood | Diffusion and vacancy properties of α-Zr[END_REF] ≥ 1.5 0.65 ± 0.05 Expt. (growth kinetics) [START_REF] Buckley | The direct observation of irradiation damage in zirconium and its alloys[END_REF] 0.65 Expt. (Huang) [START_REF] Ehrhart | Self-interstitial-atoms and their agglomerates in hcp metals[END_REF][START_REF] Ehrhart | Basic Defects in Metals[END_REF] -0.1

(E mig c
). The migration energies obtained with ab initio calculations are 0.54 eV in the basal plane and 0.65 eV out of the basal plane. This suggests a significant anisotropy of vacancy diffusion, with a fast diffusion inside the basal plane. Previous ab initio calculations in Zr [START_REF] Vérité | Anisotropy of the vacancy migration in Ti, Zr and Hf hexagonal closepacked metals from first principles[END_REF] have already shown such an anisotropy, and the same anisotropy is obtained with the EAM potentials. This is in agreement with the experimental characterization of selfdiffusion performed by Hood et al. [START_REF] Hood | α-Zr selfdiffusion anisotropy[END_REF][START_REF] Hood | Self-and Hf diffusion in alpha-Zr and in dilute, Fe-free, Zr(Ti) and Zr(Nb) alloys, Defect Diff[END_REF], who obtained a ratio of 0.6 ± 0.2 between the diffusion coefficient along the c axis and in the basal plane. In addition, the average migration energy given by ab initio calculations is in very good agreement with the experimental ones, deduced either from resistivity recovery [START_REF] Neely | Damage rate and recovery measurements on zirconium after electron irradiation at low temperatures[END_REF], positron annihilation spectroscopy (PAS) [START_REF] Hood | The recovery of single crystal α-Zr from low temperature electron irradiation -a positron annihilation spectroscopy study[END_REF][START_REF] Hood | Diffusion and vacancy properties of α-Zr[END_REF], or TEM characterization of irradiation growth [START_REF] Buckley | The direct observation of irradiation damage in zirconium and its alloys[END_REF]. Whereas the average migration energy is also good with EAM #3, EAM #2 overestimates this energy.

We also use the methodology of Ref. [START_REF] Varvenne | Point defect modeling in materials: coupling ab initio and elasticity approaches[END_REF] to deduce from our atomistic simulations the elastic dipole when the vacancy is in its stable configuration. The ab initio values (Tab. 1) indicates a contraction more important along the c axis than in the basal plane. This is consistent with the anisotropic displacement field evidenced by Ehrhart et al. [START_REF] Ehrhart | Self-interstitial-atoms and their agglomerates in hcp metals[END_REF][START_REF] Ehrhart | Basic Defects in Metals[END_REF], using Huang X-ray diffuse scattering experiments. The relaxation volume of the vacancy, δV rel , can be deduced from this elastic dipole. 2 Ab initio calculations lead to a higher relaxation , where C ij are the elastic constants of hcp Zr.

volume than the experimental value reported by Ehrhart et al. [START_REF] Ehrhart | Self-interstitial-atoms and their agglomerates in hcp metals[END_REF][START_REF] Ehrhart | Basic Defects in Metals[END_REF]. The relaxation volume given by EAM #3 potential is in good agreement with ab initio results, whereas EAM #2 potential leads to a very small relaxation volume.

Stability of small vacancy clusters

The stability of small vacancy clusters is studied by calculating their binding energies. After building the configurations corresponding to the chosen vacancy clusters, atomic positions are relaxed with a conjugate gradient algorithm. We define then the binding energy of a cluster containing n vacancies as the difference between the formation energies of n isolated vacancies and the formation energy of the cluster:

E b (V n ) =nE f (V 1 ) -E f (V n ) =nE(V 1 ) -E(V n ) -(n -1)E(bulk),
where E(V 1 ), E(V n ) and E(bulk) are the energies of the same simulation cell containing respectively one isolated vacancy, the vacancy cluster and no defect. A positive value of the binding energy indicates that the interaction between the vacancies is attractive and that the cluster is stable.

Divacancy

In order to understand how to build the vacancy clusters, we first analyse the interaction between two vacancies. Different configurations are investigated: the second vacancy is placed on the successive neighboring shells of the first vacancy, at distances lower than two lattice parameters (see Fig. for the detailed configurations). The corresponding binding energies are presented in Fig. 1.

Ab initio calculations show that the interaction is attractive only when the vacancies are first nearest neighbors. The first nearest-neighbor configuration with the divacancy lying in the 1/6 [02 23] direction (a configuration) is twice more attractive than the b configuration lying in the basal plane. All other configurations, corresponding to vacancies separated by more than one first nearest-neighbor distance, are repulsive. In particular, the d configuration, which lies along the c axis, shows a strongly negative binding energy (-0.26 eV). For the most distant e to g configurations, the magnitude of the binding energy decreases, but the interaction still remains repulsive. Similar results are found for divacancies in hcp titanium [START_REF] Raji | Ab initio pseudopotential study of vacancies and self-interstitials in hcp titanium[END_REF][START_REF] Connétable | Firstprinciples study of diffusion and interactions of vacancies and hydrogen in hcp-titanium[END_REF]: only first nearest-neighbor configurations are attractive, with however almost equal values for both a and b configurations, and the d configuration lying along the c axis is strongly repulsive.

The same divacancy configurations are investigated with the EAM potentials #2 and #3. The binding energies obtained with EAM #3 are close to zero for all the configurations (Fig. 1). The vacancies do not interact, even when the vacancies are first nearest-neighbors. This is incompatible with the DFT results and with the vacancy clustering observed in experiments. This potential is therefore not well suited to describe vacancy clustering in hcp Zr, and we will mainly ignore it in the following. On the other hand, EAM #2 potential leads to an attractive interaction when vacancies are first nearest-neighbors, and to zero binding when the vacancies are further. This is qualitatively consistent with ab initio results, although the empirical potential overestimates the divacancy stability and does not account for the repulsive character of the c to g configurations. It is also worth pointing out that the potential fails to discriminate between the a and b configurations. This is a direct consequence of the central force approximation used by EAM potentials, where no angular dependence is included. As a consequence, these potentials cannot catch the difference between configurations a and b where the two vacancies composing the divacancy have the same environment and are separated by an almost equal distance. The same limitation of EAM potentials will impact the relative stability of vacancy clusters predicted by EAM #2 potential, as it will be seen later.

Compactness of stable clusters: the tri-vacancy

We now look at how to build larger vacancy clusters. Based on the results obtained for the divacancy, only clusters formed by vacancies which are first nearest-neighbors are considered. This leads for the tri-vacancy to nine different clusters (Table 2). Both ab initio calculations and the EAM #2 potential predict positive binding energies for all these nine clusters. The most stable clusters are the compact ones, which involve the largest number of first nearest-neighbor interactions (clusters 1, 2 and 5). Two different compact configurations can be formed in a basal plane, which are crystallographically not equivalent (clusters 1 and 2). With both interaction models, the most stable one is the configuration 2. The compact tri-vacancy lying in a prismatic plane (configuration 5) has the same formation energy as the most stable basal configuration.

Relative stability of compact clusters

As compact clusters are the most stable ones, we can now separate these clusters in different groups Table 2: Total binding energy for different configurations of a tri-vacancy (in eV). The vacancies composing the cluster are sketched by squares on a basal projection of the hcp lattice, whereas spheres correspond to atoms. White and grey symbols are lying respectively in the z = 0 and z = c/2 basal planes. When necessary, squares of different sizes are used for vacancies lying in different basal planes separated by a distance c.

Configuration

Ab initio EAM #2 so as to compare their relative stability. We consider three different groups:

V 3 V 4 V 5 V 6 V 7
• basal clusters, where all vacancies are lying in the same basal plane. These clusters can be seen as precursors of c loops.

• prismatic clusters, where all vacancies are lying in the same prismatic corrugated plane (plane denoted Aα in Fig. 6a). These clusters can be seen as precursors of a loops.

• 3D clusters maximizing the number of vacancies in position of first nearest-neighbors (precursors of cavities).

We investigate different configurations for each group and retain only the most stable ones. The configurations obtained for clusters containing between 3 and 7 vacancies are shown in Fig. 2. The most stable configuration in a given group is always the one for which the number of first nearestneighbor vacancies is maximum. The binding energies for the different types of clusters containing between 2 and 7 vacancies are shown in Fig. 3. Both ab initio calculations and the EAM #2 potential show that the 3D clusters are the most stable. This is not surprising as these clusters are the most compact and therefore maximize the number of attractive interaction between vacancies in first nearest-neighbor positions. Like for the di-vacancy, ab initio calculations show that the prismatic clusters are slightly more stable than the basal ones, for clusters containing at least 5 vacancies. The empirical potential, on the other hand, predicts the same stability for both types of plane clusters. Despite this limitation, and an overestimation of cluster stability, this empirical potential manages to give a reasonable description of vacancy clusters. In particular, the configurations of the most stable clusters predicted by this potential are the same as the ab initio ones for each cluster type, and the defect structures after atomic relaxation is also equivalent. The EAM #2 potential appears therefore well-suited to study vacancy clustering in hcp Zr.

Stacking faults and surfaces

As vacancy clustering leads to faulted dislocation loops, it is worth looking at stacking fault energies before studying the stability of large vacancy clusters. We also study energies of different plane surfaces, as these surface energies will be used then to model cavities. Comparison between ab initio calculations and results obtained with empirical potentials will give insights on the ability of these potentials to model large vacancy clusters.

Basal stacking faults

Condensation of vacancies in a basal plane results in the creation of a dislocation loop of Burgers vector b 1 = 1/2 [0001]. This corresponds to the removal of a platelet of one atomic layer in the perfect stacking BABABA of basal planes and leads to the formation of a highly energetic stacking sequence, BAB.BABA, denoted BB in the following. The stacking can then evolve so as to lower the energy of the vacancy loop by creating two different stacking faults [START_REF] Hull | Introduction to Dislocations[END_REF]: an extrinsic fault E, which corresponds to the stacking BABCABA and does not change the Burgers vector of the dislocation loop, or an intrinsic fault I 1 , which corresponds to the stacking BABCBCB and leads to a dislocation loop with Burgers vector b 2 = 1/6 20 23 . To better understand the formation and stability of these different stacking faults we use the concept of generalized stacking faults [START_REF] Vitek | Intrinsic stacking faults in body-centred cubic crystals[END_REF][START_REF] Vitek | Non-planar dislocation cores: A ubiquitous phenomenon affecting mechanical properties of crystalline materials[END_REF].

Extrinsic stacking fault

The extrinsic fault E is formed from the BAB.BABA stacking by the glide of one atomic plane (Fig. 4a). We compute the stacking fault energy for different glide vectors F lying in the basal plane. In these calculations, atoms are allowed to relax only in the direction perpendicular to the fault plane. We used a stacking of 15 (0001) planes in the ab initio calculations, which corresponds to a distance h 0001 = 15c/2 between fault planes and is high enough to prevent any interaction between the fault plane and its periodic images. Generalized stacking fault energies are calculated on a regular 10 × 10 grid and are then interpolated with Fourier series.

The obtained energy as a function of the fault vector, or γ-surface (Fig. 4b, c andd), shows a minimum for a fault vector 2/3 [1 100] which corresponds to the metastable extrinsic stacking fault E. This is the only minimum which exists on the ab initio γ-surfaces. In particular, the BB stacking, corresponding to a fault vector 0 or 1/3 [1 100], is unstable. This is more clearly seen on Fig. 4e which corresponds to a plot of the fault energy along the [1 100] direction. On the other hand, both empirical potentials EAM #2 and #3 predict that the BB stacking is an energy minimum. This artifact of empirical potentials leads them to stabilize the BB stacking for small vacancy loops, whereas one expects from the ab initio results that such a BB stacking will naturally relax to create an extrinsic E fault. This may be the reason why special relaxation techniques had to be used in Ref. [START_REF] De Diego | On the structure and mobility of point defect clusters in alpha-zirconium: a comparison for two interatomic potential models[END_REF] to obtain the stable structure of vacancy clusters lying in the basal planes.

The minimum energy, corresponding to the extrinsic stacking fault, is γ E = 274 mJ.m -2 , with ab initio calculations, in good agreement with already published values [START_REF] Domain | Atomic-scale ab initio study of the Zr-H system: II. interaction of H with plane defects and mechanical properties[END_REF][START_REF] Poty | Classification of the critical resolved shear stress in the hexagonal-closepacked materials by atomic simulation: Application to α-zirconium and α-titanium[END_REF]. The empirical potential EAM #2 underestimates this fault energy whereas a good agreement is obtained with EAM #3 (Tab. 3).

Intrinsic stacking fault

An intrinsic stacking fault of type I 1 is created from the BAB.BABA stacking by shifting one part of the crystal by a fault vector F lying in the basal plane (Fig. 5a). The I 1 fault corresponds to a fault

Ab initio

EAM Pwscf Ref. [START_REF] Domain | Atomic-scale ab initio study of the Zr-H system: II. interaction of H with plane defects and mechanical properties[END_REF] Ref. [START_REF] Udagawa | Ab initio study on plane defects in zirconium-hydrogen solid solution and zirconium hydride[END_REF] Ref. [ 5b) and that the BB stacking is an energy maximum. The EAM #2 potential leads to a γ surface in qualitative agreement (Fig. 5c). On the other hand, EAM #3 potential predicts that both the BB stacking and the I 1 fault are energy minima (Fig. 5d). Like for the extrinsic generalized stacking fault, this empirical potential leads to a stabilization of the BB stacking versus the formation of an intrinsic I 1 fault. It is worth pointing out that a similar artifact had already been mentioned for the basal I 2 γ-surface with this potential [START_REF] Clouet | Screw dislocation in zirconium: An ab initio study[END_REF].

Ab initio calculations lead to an energy minimum γ I1 = 147 mJ.m -2 , still in good agreement with previous ab initio calculations [START_REF] Domain | Atomic-scale ab initio study of the Zr-H system: II. interaction of H with plane defects and mechanical properties[END_REF][START_REF] Poty | Classification of the critical resolved shear stress in the hexagonal-closepacked materials by atomic simulation: Application to α-zirconium and α-titanium[END_REF] (Tab. 3). Both EAM #2 and #3 potentials underestimates this fault energy, with the larger error for EAM #2 (Fig. 5e). All energy models lead to the following order between the energies of the different basal stacking faults: γ I1 < γ I2 < γ E , in agreement with predictions based on an analysis of broken bonds between pairs of atoms [START_REF] Hirth | Theory of Dislocations, 2nd Edition[END_REF][START_REF] Hull | Introduction to Dislocations[END_REF]. Such an ordering of the fault energies was not retrieved by the long ranged pair potential used in Ref. [START_REF] Kapinos | Computer simulation of vacancy loops and stacking faults in zirconium[END_REF]. As EAM #2 and #3 are short ranged potentials relying on a central force approximation, they naturally lead to the relation 3γ I1 = 2γ I2 = γ E [START_REF] Hirth | Theory of Dislocations, 2nd Edition[END_REF]. Our ab initio calculations show that such a relation is only approximate (Tab. 3) and that the angular contribution of the atomic interaction causes deviations from this idealized picture. As first pointed out by Legrand [START_REF] Legrand | Relations entre la structure électronique et la facilité de glissement dans les métaux hexagonaux compacts[END_REF], a fully predictive modeling of atomic interactions in hcp transition metals like Zr needs a proper account of these angular contributions, and thus to go beyond simple empirical potentials of the EAM type.

Prism stacking faults

When removing a vacancy platelet in a corrugated {10 10} prismatic plane, a prismatic stacking fault is formed (Fig. 6a). This stacking fault, which is associated with the creation of a dislocation loop of Burgers vector b 3 = 1/2 10 10 , is the same as the one involved in the dissociation in a prismatic plane of a 1/3 1 210 dislocation. The unfaulting of the vacancy loop occurs by a 1/6 1 210 shearing of the fault plane, which leads to a perfect vacancy loop with Burgers vector b 4 = 1/3 2 11 0 (Fig. 6b). The γ-surface associated with this unfaulting mechanism has already been described in Ref. [START_REF] Clouet | Screw dislocation in zirconium: An ab initio study[END_REF], both for ab initio calculations and EAM #3 potential. Ab initio leads to an energy γ (10 10) = 211 mJ.m -2 for the metastable stacking fault in this prismatic plane. EAM #2 and #3 respectively overestimates and underestimates this fault energy (Tab. 3).

Surface energies

Large cavities observed by TEM in Zr irradiated with electrons show facets in the basal (0001), prismatic {10 10} and pyramidal {10 11} planes [START_REF] Griffiths | Hvem study of the effects of alloying elements and impurities on radiation damage in Zr-alloys[END_REF][START_REF] Griffiths | Study of point defect mobilities in zirconium during electron irradiation in a high-voltage electron microscope[END_REF]. We now consider the surface energies for these three different planes. For each surface of interest, a crystal block cut in the suitable planes is placed into vacuum. For DFT calculations, supercells contain ∼ 12 atomic layers, and the vacuum slab is ∼ 10 Å thick. This ensures the convergence of the surface energies.

Table 3 displays the basal, prismatic and pyramidal surface energies, calculated with the different interaction models. Our ab initio calculations show that the pyramidal surface has the lowest energy, with the basal surface being the next ones. This agrees with TEM observations showing that large cavities have facets mainly in the pyramidal and basal planes [START_REF] Griffiths | Hvem study of the effects of alloying elements and impurities on radiation damage in Zr-alloys[END_REF][START_REF] Griffiths | Study of point defect mobilities in zirconium during electron irradiation in a high-voltage electron microscope[END_REF]. Our ab initio results for the basal and prismatic surfaces match well those of previous ab initio studies [START_REF] Domain | Atomic-scale ab initio study of the Zr-H system: II. interaction of H with plane defects and mechanical properties[END_REF][START_REF] Udagawa | Ab initio study on plane defects in zirconium-hydrogen solid solution and zirconium hydride[END_REF] (Tab. 3).

The EAM #2 potential underestimates these surface energies and predicts that the basal surface is the most stable one, instead of the pyramidal surface. Values predicted by EAM #3 are closer to ab initio results, but this potential does not really discriminate between the different plane surfaces.

Stability of large vacancy clusters: dislocation loops and cavities

We now study the stability of larger vacancy clusters, cavities and dislocation loops lying either in the basal or prismatic planes. As pointed out in the introduction, cavities are hardly observed in irradiated zirconium [START_REF] Griffiths | Hvem study of the effects of alloying elements and impurities on radiation damage in Zr-alloys[END_REF][START_REF] Griffiths | Study of point defect mobilities in zirconium during electron irradiation in a high-voltage electron microscope[END_REF] and vacancies mainly condensate in the form of dislocation loops [START_REF] Northwood | Characterization of neutron irradiation damage in zirconium alloys -an international "round-robin" experiment[END_REF][START_REF] Griffiths | Neutron damage in zirconium alloys irradiated at 644 to 710 K[END_REF][START_REF] Griffiths | A review of microstructure evolution in zirconium alloys during irradiation[END_REF][START_REF] Onimus | Radiation effects in zirconium alloys[END_REF]. The formation of a loops lying in the prismatic planes seems to be more favorable than the formation of c loops in the basal planes. a loops already appear at low irradiation doses, whereas c loops only appear for larger irradiation doses, when the irradiation growth of the crystal accelerates. We propose to examine if these experimental observations can be understood through stability arguments.

Ab initio calculations cannot be used to study such large clusters. On the other hand, the two previous sections have shown that empirical potentials suffer from limitations. We therefore propose to use an hybrid approach to model large vacancy clusters.

We calculate the formation energy of each type of defect for various cluster sizes (up to ∼ 380 vacancies) with the EAM #2 empirical potential. We choose this potential because it accounts for the binding between vacancies and it reasonably describes the relative stability of the vacancy clusters. (cf. §3). Results are then used to validate energy models based on a continuous description of vacancy clusters. We finally parameterize these continuous models with quantities deduced from ab initio calculations. It allows us to extrapolate the DFT results to larger sizes and to discuss the relative stability of the different vacancy defects.

Introduction of vacancy loops in atomistic simulations

Vacancy loops are introduced in our atomistic simulations by first removing the atoms inside the vacancy platelets and then applying to all atoms in the simulation box the displacement field predicted by elasticity theory for the corresponding dislocation loop. The displacement created by a dislocation loop of Burgers vector b is given by the Burgers formula [START_REF] Hirth | Theory of Dislocations, 2nd Edition[END_REF]:

u i ( x) = - b i Ω( x) 4π + L ... (1)
where Ω( x) is the solid angle subtended by the loop area at x. It corresponds to the plastic displacement created by the loop and is a purely geometrical term. The second term is a closed line integral which accounts for the elastic relaxation. It can be evaluated using either isotropic [START_REF] Barnett | The displacement field of a triangular dislocation loop[END_REF][START_REF] Barnett | The displacement field of a triangular dislocation loop -a correction with commentary[END_REF] or anisotropic [START_REF] Lazar | Dislocation loops in anisotropic elasticity: displacement field, stress function tensor and interaction energy[END_REF] elasticity theory. As Eq. ( 1) is only used to generate the initial configuration, which is then relaxed with the empirical potential, we only retain the plastic part of the displacement field. The solid angle is calculated with the closed-form expression given by Van Oosterom [START_REF] Van Oosterom | The solid angle of a plane triangle[END_REF].

Experimentally [START_REF] Griffiths | A review of microstructure evolution in zirconium alloys during irradiation[END_REF], a loops are circular for radius below 40 nm and elliptic above. No precise information is available for the shape of the c loops. As these c loops are formed with a background of numerous a loops, one needs to choose image conditions where a loops are invisible to see c loops in TEM. As a consequence, these c loops are usually imaged on their edge. Dislocation loops are introduced in our simulation boxes as hexagonal loops. This morphology is reasonable with regards The symbols correspond to the results of atomistic simulations and the lines to continuous laws: Eq. 5 for dislocation loops and Eq. 6 for cavities.

to the experimental data. Previous atomistic studies have also shown that the formation energy of the loops only slightly depends on their shape [START_REF] De Diego | Structure and properties of vacancy and interstitial clusters in α-zirconium[END_REF].

In previous works [START_REF] Kapinos | Computer simulation of vacancy loops and stacking faults in zirconium[END_REF][START_REF] Kulikov | Vacancy dislocation loops in zirconium and their interaction with self-interstitial atoms[END_REF][START_REF] De Diego | Structure and properties of vacancy and interstitial clusters in α-zirconium[END_REF][START_REF] De Diego | On the structure and mobility of point defect clusters in alpha-zirconium: a comparison for two interatomic potential models[END_REF], the vacancy loops were obtained by removing vacancy platelets in the relevant planes. Atomic relaxations, eventually followed by annealing sequences, were used to find the stable configurations. With the procedure used here, based on the Burgers formula (Eq. 1), it is possible to introduce separately each kind of loop and to control the stacking fault created by the loop by choosing the corresponding Burgers vector.

Basal dislocation loops

We now examine the stability of the different vacancy loops lying in the basal planes. Loops with a BB stacking fault are formed when choosing a Burgers vector b 1 = 1/2 [0001]. Loops with an intrinsic fault I 1 are formed with b 2 = 1/6 20 23 . To build loops with an extrinsic fault E, we use the recipe of Hull and Bacon [START_REF] Hull | Introduction to Dislocations[END_REF]. Two loops of the same size separated by one atomic layer are formed on top of the other with Burgers vectors 1/12 40 43 and 1/12 4043 . This results in a loop with Burgers vector b 1 = 1/2 [0001], but with an extrinsic fault.

With the EAM #2 potential, the BB loops are stable for small sizes, but they are the less energetically favorable defects (Fig. 7). They become unstable for clusters containing more than 160 vacancies. The E loops have a lower formation energy than the I 1 loops for the small sizes and a greater energy for the largest sizes (Fig. 7). This is in agreement with a competition between the stacking fault energy (γ E > γ I1 ) and the elastic energy ( b 1 2 < b 2 2 ), as it will be seen below. The stability crossover between both types of loops occurs around 200 vacancies, which corresponds to a loop radius R = 2.4 nm.

In order to extrapolate the data for loops containing a larger number of vacancies, we consider a line tension model [START_REF] Bacon | Anisotropic continuum theory of lattice defects[END_REF]. Within this model, the energy of a dislocation loop of radius R is given by:

E f loop (R) = πR 2 γ + L K( t)ds ln R r c , (2)
where γ is the stacking fault energy and r c the core radius of the dislocation loop. The coefficient K( t) is the factor appearing in front of the logarithm when defining the elastic energy of a straight dislocation. It only depends on the bulk elastic constants, the Burgers vector of the loop, and the direction of the loop element ds through its local tangent t. We calculate it according to anisotropic elasticity using Stroh sextic formalism [START_REF] Stroh | Dislocations and cracks in anisotropic elasticity[END_REF][START_REF] Stroh | Steady state problems in anisotropic elasticity[END_REF][START_REF] Clouet | Elastic energy of a straight dislocation and contribution from core tractions[END_REF]. We assume a circular shape to calculate the closed line integral and define an average value of this coefficient:

K = 1 2π 2π 0 K( t)dθ, (3)
Analytical expressions are available for these dislocation loops lying in the basal plane (Appendix A.2), but the numerical evaluation has to be done in the more general case, in particular for the loops lying in the prismatic plane that will be considered below. Approximate expressions for this coefficient can also be obtained if one assumes isotropic instead of anisotropic elasticity (Appendix A.1). The exact shape of the loop is considered below through a shape factor f used as a fitting parameter. f = 1 for a circular loop. Deviations from this ideal value occur for non circular loops because both the perimeter, as a function of the number of vacancies, and the average energy coefficient (Eq.

3) depend on the exact shape of the loop. The link between the loop radius R and the number n of vacancies constituting the loop is established considering that n vacancies occupy the same volume as a disk of radius R and of thickness b e , the edge component of the loop Burgers vector normal to its habit plane. For the three basal loops, we have b e = c/2. Considering that the volume of each vacancy is one atomic volume Ω = √ 3/4 a 2 c, the equality πR 2 b e = nΩ leads to the relation:

R = a √ 3n 2π 1/2 = R 1 √ n, (4)
with R 1 = a( √ 3/2π) 1/2 . The continuous expression of the basal dislocation loops energy can therefore be expressed as a function of their number n of vacancies:

E f loop (n) = πR 1 2 γn + 2πf R 1 K√ n ln R 1 √ n r c .
(5) As expected, for large loops (large n), the stacking fault energy represents the dominant contribution to the formation energy. The shape factor f and the core radius r c are used here as fitting parameters, in order to obtain the best agreement between the continuous expression 5 and the results of atomistic simulations. In this fitting procedure, the stacking fault energies γ and the elastic coefficients K are fixed to their values calculated with EAM #2 for the corresponding basal loop (Tables 3 and4). Despite the simplicity of the line tension model, the agreement with the atomistic results is good (Fig. 7). The line tension model perfectly fits atomistic simulations for clusters containing at least 50 vacancies. Some discrepancies appear for smaller clusters (cf. inset in Fig. 7), but the predictions of the line tension model are still reasonable. The fitted parameters are given in Table 4: the shape factor only slightly deviates from its ideal value (f = 1), and the core radius r c is close to the norm of a Burgers vector (b e = c/2), as expected from elasticity theory. This shows that the integration of the exact shape of the loops through these two fitting parameters is a reasonable procedure. No attempt was made to calculate them exactly, as it would require a more complex treatment within elasticity theory [START_REF] Bullough | The elastic energy of a rhombus-shaped dislocation loop[END_REF][START_REF] Bacon | The anisotropic elastic energy of a rhombus-shaped dislocation loop[END_REF][START_REF] Schoeck | The elastic energy of dislocation loops in anisotropic media[END_REF] than the simple line tension model used here. We also note that variations of these parameters between the different loops are small.

Prismatic dislocation loops

We now look at vacancy loops lying in the prismatic {10 10} planes. Faulted loops are created with a Burgers vector b 3 = 1/2 10 10 and perfect loops with b 4 = 1/3 2 11 0 .

The stability of these prismatic loops, as predicted by atomistic simulations using the EAM # Table 4: Parameters defining the formation energy of the different vacancy loops (Eq. 5). The core radii rc (normalized by the lattice parameter a) and the shape factors f have been obtained by fitting the atomistic results of EAM #2 potential. The elastic coefficients K (in eV/ Å) are deduced from the elastic constants, corresponding either to EAM #2 potential or ab initio calculations, using anisotropic elasticity and Eq. 3. 7. Perfect loops are unstable towards the faulted loops when they contain less than 150 vacancies (R = 2 nm). In the range where both types of loops are stable, an inversion of stability is observed when increasing the loop size, at a size corresponding to ∼ 250 vacancies (R = 2.7 nm). For small clusters, faulted loops are the most stable ones because they have a smaller Burgers vector, whereas perfect loops are more stable for large defects, the stacking fault becoming too costly. Both types of prismatic loops are more stable than the basal loops. This differs from what was obtained in previous atomistic simulations [START_REF] Kapinos | Computer simulation of vacancy loops and stacking faults in zirconium[END_REF][START_REF] De Diego | Structure and properties of vacancy and interstitial clusters in α-zirconium[END_REF][START_REF] De Diego | On the structure and mobility of point defect clusters in alpha-zirconium: a comparison for two interatomic potential models[END_REF], using different empirical potentials.

r c f K EAM Pwscf Basal: fault BB 0.
We again compare the atomistic results with a continuous law for the dislocation loop energy. The expression is still given by Eq. 5, but with R 1 = ac/2π now. Using the core radius r c and the shape factors f as fitting parameters, we obtain a perfect agreement between the continuous laws and the EAM # 2 energies (Fig. 7). The parameters obtained through this fitting procedure are given in Table 4: like for the basal loops, these parameters have reasonable values.

Cavities

We now use the EAM #2 potential to study the stability of cavities. They are introduced in simulation boxes as spherical vacancy clusters of increasing size. The formation energies are shown in Fig 7. We find that cavities are always more stable than the vacancy-loops, whatever their nature.

The same result was obtained in Ref. [START_REF] Kulikov | Vacancy dislocation loops in zirconium and their interaction with self-interstitial atoms[END_REF] with a different atomic potential.

To interpolate these results of atomistic simulations, we consider the formation energy of a spherical cavity, taking into account only its surface energy. This leads for a cavity containing n vacancies to

E f cav (n) = 4π a 2 c 3 √ 3 16π 2/3 f σ n 2/3 . (6)
f is a geometrical factor which is equal to 1 for a spherical cavity. It will be taken as a fitting parameter. The surface energy σ appearing in this expression is an average energy. It can be defined from the energies of plane surfaces, σ 0001 , σ 10 10 and σ 10 11 , using the Wulff construction (cf. Appendix B). This construction ensures that the ideal spherical cavity considered in Eq. 6 has the same surface energy as the real faceted cavity. Using the values predicted by EAM #2 potential for plane surfaces (Tab. 3), we obtain σ = 1420 mJ.m -2 . The results of atomistic simulations are then perfectly reproduced by Eq. 6 with a shape factor f = 1.03 (Fig. 7). The value obtained for this fitting parameter is close to its ideal value (f = 1). This shows the validity of our modeling, despite its simplicity. The continuous laws for the defects energetics, even if their expressions are simple, fit well the atomistic results obtained with the EAM #2 potential. We therefore use the same laws, but now with parameters deduced from ab initio calculations. The elastic coefficients K and the stacking fault energies are fixed to their DFT values. The surface energy needed for cavities is derived from the same Wulff construction, using ab initio results for plane surfaces. This leads to σ = 1690 mJ.m -2 , a value slightly lower than the experimental estimate at 0K, 1900 mJ.m -2 [START_REF] Tyson | Surface free energies of solid metals: Estimation from liquid surface tension measurements[END_REF]. The remaining parameters of the continuous models, i.e. the core radii r c and the shape factors f , are kept fixed to their values deduced from the fitting of EAM #2 results, as these parameters could not be calculated ab initio. The resulting stability curves are displayed in Fig. 8.

Ab initio modeling

This ab initio based modeling predicts that the loops lying in the prismatic planes are the most stable defects for the whole range of defect sizes. This is in agreement with experiments showing that these prismatic a loops are the easiest vacancy clusters to create [START_REF] Onimus | Radiation effects in zirconium alloys[END_REF]. For the smallest sizes, the faulted prismatic loops have the lowest formation energy and when the loop radius becomes larger than 5.2 nm, perfect prismatic loops become more stable (cf. inset in Fig. 8).

Contrary to what was obtained with the EAM #2 potential, cavities are not the most stable defects. This disagreement of the empirical potential with ab initio calculations arises both from an underestimation of the surface energy (σ = 1420 mJ.m -2 with EAM #2 instead of 1690 mJ.m -2 with ab initio) and from an overestimation of the prismatic stacking fault energy (γ 10 10 = 357 mJ.m -2 with EAM #2 instead of 211 mJ.m -2 ). As cavities are observed only in a few cases in pure Zr [START_REF] Griffiths | A review of microstructure evolution in zirconium alloys during irradiation[END_REF][START_REF] Griffiths | Hvem study of the effects of alloying elements and impurities on radiation damage in Zr-alloys[END_REF][START_REF] Griffiths | Study of point defect mobilities in zirconium during electron irradiation in a high-voltage electron microscope[END_REF], ab initio predictions appear more reliable than results of EAM #2 potential. According to this ab initio modeling, cavities are the less stable defects also at very small sizes. This contradicts our DFT calculations for small clusters containing up to 7 vacancies, as these calculations indicate that volume aggregate are more stable than plane clusters (Fig. 3). This illustrates the limitations of the continuous laws used to describe cluster energy. Such simple laws are valid only for large enough clusters. At small sizes, a full atomic description is needed.

We finally comment on loops lying in the basal planes. As ab initio calculations show that the BB fault is unstable, BB loops are not considered. Like with the EAM #2 potential, the basal loops are less stable than loops lying in the prismatic planes. The ab initio model predicts a stability inversion between the E and I 1 basal loops for R = 1.4 nm.

Conclusion

The stability of vacancy clusters in pure hcp zirconium has been studied using an atomistic modeling approach based on both ab initio calculations and empirical potential. DFT calculations performed for small vacancy clusters show that the interaction between vacancies is attractive only when they are first nearest neighbors. Such an interaction is partly additive, leading to a higher stability for compact clusters. The empirical potential EAM #2 derived by Mendelev and Ackland [START_REF] Mendelev | Development of an interatomic potential for the simulation of phase transformations in zirconium[END_REF] gives a reasonable description of vacancy clustering. This potential allows us to study clusters containing up to 380 vacancies, so as to validate simple analytic laws giving a continuous description of the formation energy for the different configurations. A good agreement has been obtained between these two descriptions, at an atomic and at a continuous scale. We could then parameterize such continuous laws with ab initio calculations. This shows that the most stable vacancy clusters are dislocation loops, either perfect or faulted, lying in the prismatic planes. This is in agreement with experiments where such a loops are usually the easier to form.

The continuous laws, which were used to describe the formation energy of the different vacancy clusters, are classical models which can be found in any metallurgy textbook. In particular, Eq. 5 for planar clusters directly derives from dislocation theory [START_REF] Hirth | Theory of Dislocations, 2nd Edition[END_REF], within an anisotropic line tension approximation. We emphasize the usefulness of such laws, which allow discussing the relative stability of different loops, like we did, and like it was previously done in Ref. [START_REF] Yoshida | Dislocation loops with stacking fault in quenched aluminum[END_REF] for vacancy loops in quenched aluminum or Ref. [START_REF] Dudarev | Effect of the α -γ phase transition on the stability of dislocation loops in bcc iron[END_REF] for interstitial loops in irradiated iron.

These laws then offer a convenient framework to model the kinetics of point defect clustering, using for instance cluster dynamics simulations [START_REF] Kiritani | Analysis of the clustering process of supersaturated lattice vacancies[END_REF]. As the variation of the energy with the size of the loop differs from the one corresponding to 3D clusters, like precipitates, the long time evolution leads to a coarsening regime different from the usual LSW theory. The corresponding coarsening model, developed by Kirchner, Burton and Speight [START_REF] Kirchner | Size distribution of dislocation loops[END_REF][START_REF] Burton | The coarsening and annihilation kinetics of dislocation loop[END_REF], has been shown to agree with experimental data [START_REF] Powell | A determination of the elastic energy of dislocation loops from coarsening kinetics[END_REF][START_REF] Liu | Evolution of dislocation loops in silicon in an inert ambient-i[END_REF][START_REF] Bonafos | Ostwald ripening of end-of-range defects in silicon[END_REF]. To develop such a modeling of the kinetic evolution of irradiated Zr, it will be necessary however to parameterize the same type of continuous laws for interstitial clusters. Finally, these laws are also a good way to study the influence of impurities and alloying elements on defect stability, through the modification of stacking fault or surface energies [START_REF] Domain | Atomic-scale ab initio study of the Zr-H system: II. interaction of H with plane defects and mechanical properties[END_REF][START_REF] Udagawa | Effect of sn and nb on generalized stacking fault energy surfaces in zirconium and gamma hydride habit planes[END_REF].

1

 1

3 -

 3

Figure 1 :

 1 Figure 1: Binding energies of a divacancy calculated with the different energy models. The corresponding divacancy configurations are sketched on the hcp lattice, where V 1 denotes the position of the first vacancy and the letters a to g the position of the second one.

Figure 2 :

 2 Figure 2: Most stable configurations obtained for basal, prismatic and 3D clusters containing between 3 and 7 vacancies. The same conventions as in Tab. 2 are used to represent atoms and vacancies.

Figure 3 :

 3 Figure 3: Total binding energies for the different types of vacancy clusters (basal, prismatic and 3D) calculated with different energy models: (a) ab initio calculations and (b) EAM #2.

Figure 4 :

 4 Figure 4: Extrinsic E generalized stacking fault. (a) Formation mechanism of the stacking fault. (b) Ab initio, (c) EAM #2 and (d) EAM #3 γ-surfaces. (e) Comparison of the fault energies obtained with the different energy models along the [1 100] direction.

Figure 5 :

 5 Figure 5: Intrinsic I 1 generalized stacking fault. (a) Formation mechanism of the stacking fault. (b) Ab initio, (c) EAM #2 and (d) EAM #3 γ-surfaces. (e) Comparison of the fault energies obtained with the different energy models along the [1 100] direction.

Figure 6 :

 6 Figure 6: Unfaulting mechanism in prismatic planes. (a) The removal of one corrugated (10 10) plane creates a prismatic stacking fault. (b) The prismatic fault is removed by a shear of amplitude 1/6 [1 210] in the (10 10) plane.

Figure 7 :

 7 Figure 7: Formation energy of large vacancy clusters in Zrpredicted by EAM #2 potential. Hexagonal loops lying either in the basal (0001) or the prismatic {10 10} planes are considered with different stacking faults, as well as spherical cavities. The symbols correspond to the results of atomistic simulations and the lines to continuous laws: Eq. 5 for dislocation loops and Eq. 6 for cavities.

Figure 8 :

 8 Figure 8: Formation energies of large vacancy clusters predicted by continuous laws parameterized on DFT results. The inset shows the stability inversion between faulted and perfect loops lying in the prismatic planes.

Table 1 :

 1 Vacancy properties in hcp Zr: formation energies E for , migration energies E mig bas and E mig c

  C 33 -C 13 )(P 11 +P 22 )+(C 11 +C 12 -2C 13 )P 33 (C 11 +C 12 )C 33 -2C 13 2

	2 δV rel	=

(

Table 3 :

 3 Surface energies, σ 0001 , σ 10 10 , and σ 10 11 , and stacking fault energies, γ E , γ I 2 , γ I 1 , and γ 10 10 , in hcp Zr. The results of the present work, calculated either with Pwscf or the empirical potentials EAM #2 and #3 are compared to other ab initio values from the literature obtained with Vasp[START_REF] Domain | Atomic-scale ab initio study of the Zr-H system: II. interaction of H with plane defects and mechanical properties[END_REF][START_REF] Udagawa | Ab initio study on plane defects in zirconium-hydrogen solid solution and zirconium hydride[END_REF][START_REF] Poty | Classification of the critical resolved shear stress in the hexagonal-closepacked materials by atomic simulation: Application to α-zirconium and α-titanium[END_REF]. All energies are given in mJ.m -2 .

  Ab initio calculations show that this is the only minimum which exists on the corresponding γ surface (Fig.

					35]	#2	#3
	σ 0001	1600	1560	1600	-	1270 1540
	σ 10 10	1670	1640	1660	-	1340 1540
	σ 10 11	1550	-	-	-	1340 1550
	γ E	274	249	-	300	164	297
	γ I2	213	200	227	228	110	198
	γ I1	147	124	-	168	55	99
	γ 10 10	211	145	197	-	357	135
	vector 1/3 [1 100] or 2/3 [1 100].

Acknowledgments

This work was performed using HPC resources from GENCI-[CINES/CCRT/IDRIS] (Grant 2013-096847). AREVA is acknowledged for financial support.

Appendix A. [START_REF] Onimus | Radiation effects in zirconium alloys[END_REF]

. Isotropic elasticity

The exact calculation of the elastic coefficient K appearing in the definition of the loop energy (Eq. 5) can be tough. We use in this appendix isotropic elasticity to obtain simple analytical expressions which can be used as a first approximation. Within isotropic elasticity, the factors defining the elastic energy of an edge or a screw dislocation are respectively

where µ is the shear modulus and ν is Poisson's ratio. For a specific material, a proper choice has to be done to obtain these average elastic constants.

Basal loops with an extrinsic fault E are pure prismatic loops with b e = c/2. One simply gets in this case

, whatever the shape of the loop. Components of the Burgers vector for basal loops with an intrinsic fault I 1 are varying along the loop. Using an angle θ to define the direction of the loop tangent, one can write b e = c/2 + a sin (θ) √ 3/3 and b s = a cos (θ) √ 3/3. The average defined by Eq. 3 leads for a circular loop to

Faulted loops lying in the prismatic planes are pure prismatic loops with b e = a √ 3/2, and thus

Finally, for perfect loops lying in the prismatic planes b e = a[ √ 3 + sin (θ)]/2, and b s = a sin (θ)/2. This leads to

Appendix A.2. Basal loops Thanks to the transverse isotropy of the hexagonal crystal, one can take full account of the elastic anisotropy and derive exact expressions of the elastic energy for dislocation loops lying in the basal planes [START_REF] Chou | The energy and line tension of a dislocation in a hexagonal crystal[END_REF]. When the hcp crystal is oriented with the x, y, and z axis respectively along the [10 10], [0001], and [1 210] directions, the Stroh matrix defining the elastic energy of a dislocation lying along the z direction is diagonal with

where C11 = √ C 11 C 33 . The elastic coefficient of a basal loop with an extrinsic fault E is then

For a basal loop with an intrinsic fault I 1 , one gets

No analytical expression is available for loops lying in the prismatic planes. For these loops, one needs either to perform a numerical evaluation, like the one of the present study ( §5.2), or to use the previous approximations based on isotropic elasticity.

Appendix B. Wulff construction

We use the Wulff construction [START_REF] Christian | The Theory of Transformations in Metals and Alloys -Part I: Equilibrium and General Kinetic Theory[END_REF][START_REF] Porter | Phase Transformations in Metals and Alloys[END_REF] to define an isotropic surface energy σ from the surface energies σ 0001 , σ 10 10 , and σ 10 11 corresponding respectively to the basal, prismatic and pyramidal planes. Such a construction predicts that the equilibrium shape of cavities is faceted. Considering facets only in the (0001), {10 10}, and {10 11} planes, the surface of each facet type is proportional to