Emmanuel Clouet
email: emmanuel.clouet@cea.fr

Screw dislocation in zirconium: an ab initio study

Keywords: numbers: 61, 72, Lk, 61, 72, Bb

Plasticity in zirconium is controlled by 1/3 1 210 screw dislocations gliding in the prism planes of the hexagonal close-packed structure. This prismatic and not basal glide is observed for a given set of transition metals like zirconium and is known to be related to the number of valence electrons in the d band. We use ab initio calculations based on the density functional theory to study the core structure of screw dislocations in zirconium. Dislocations are found to dissociate in the prism plane in two partial dislocations, each with a pure screw character. Ab initio calculations also show that the dissociation in the basal plane is unstable. We calculate then the Peierls barrier for a screw dislocation gliding in the prism plane and obtain a small barrier. The Peierls stress deduced from this barrier is lower than 21 MPa, which is in agreement with experimental data. The ability of an empirical potential relying on the embedded atom method (EAM) to model dislocations in zirconium is also tested against these ab initio calculations.

I. INTRODUCTION

Plasticity in α-zirconium is controlled by dislocations with a 1/3 1 210 Burgers vector gliding in the prism planes of the hexagonal compact (hcp) lattice. [1][2][3][START_REF] Caillard | Thermally Activated Mechanisms in Crystal Plasticity[END_REF] The relative ease of prismatic glide compared to basal glide has been shown to be linked to the ratio of the corresponding stacking fault energies, which in turn is controlled by the electronic structure. [START_REF] Legrand | [END_REF] In particular, Legrand [START_REF] Legrand | [END_REF] used a tight binding model to show that prismatic slip in transition metals of the IV B column (Zr, Ti, Hf) originates from the electronic filling of the valence d band. As a consequence, it appears necessary to take into account the anisotropy of the d orbital, and hence the angular dependence of the atomic bonding, to model dislocations in these transition metals 6 . One cannot therefore rely on central forces empirical potential and needs to take account of the electronic structure. Tight binding models 7,8 or ab initio calculations [START_REF] Ferrer | Zirconium in the nuclear industry: thirteenth international symposium[END_REF][START_REF] Domain | Mat. Res. Soc. Symp. Proc[END_REF][START_REF] Domain | [END_REF][12][13] show indeed that a 1/3 1 210 screw dislocation, either in Zr or in Ti, spreads in prism planes, in agreement with the prismatic glide observed experimentally. But none of these atomistic simulations calculate the Peierls stress of a screw dislocation. According to some authors, 7,14 its core structure is not completely planar, which may be the cause of a high Peierls stress.

It is true, experimentally, that screw dislocations glide with difficulty compared to other dislocation characters in zirconium or titanium alloys: characteristic microstructures, with long and straight dislocations aligned along their screw orientation, are observed at low temperature, 3,[START_REF] Ferrer | Zirconium in the nuclear industry: thirteenth international symposium[END_REF][START_REF] Bailey | [END_REF][16][17][18] and the flow stress is strongly temperature dependent, 2,3,16,[19][20][21][22][23] in agreement with the assumption of a high Peierls barrier which must be overcome by the nucleation of double kinks. But experiments also show that the yield stress in zirconium or in titanium strongly decreases with a decreasing amount of interstitial impurities like oxygen. 2,3,16,19,22,23 It is therefore probable that the high Peierls stress of screw dis-locations has an extrinsic cause. In pure zirconium or pure titanium, this Peierls stress may not be as high and screw dislocations are probably gliding as easily as other orientations.

Recently, Mendelev and Ackland 24 developed an empirical potential for Zr in the Embedded Atom Method (EAM) formalism. Using this potential, Khater and Bacon 25 showed that it leads to a screw dislocation that spontaneously spreads in the prism plane, the configuration dissociated in the basal plane being metastable. They also showed that the Peierls stress of a screw dislocation gliding in the prism plane is not so different from the Peierls stress of an edge dislocation and that this stress is small (22 MPa for the screw and 16 MPa for the edge). These results therefore support experimental findings stating that screw dislocations are gliding in prism planes with a low Peierls stress in pure zirconium. But these simulations rely on a central forces potential, which is not well suited to describe dislocations in a hcp transition metal like Zr, as described above. More reliable atomistic simulations, incorporating a description of the electronic structure, are therefore needed to confirm this easy glide of screw dislocations in pure zirconium. This article aims to use ab initio calculations so as to fully characterize the core structure of a 1/3 1 210 screw dislocation in zirconium and estimate its Peierls stress. We also examine generalized stacking faults as dislocation core structures are closely related to them. Two different ab initio methods, Siesta 26 and Pwscf 27 , have been used. Siesta offers the advantage of efficiency, allowing simulating more atoms than with a standard ab initio code, whereas Pwscf offers the advantage of robustness. All calculations are also performed with Mendelev and Ackland EAM potential, 24 so as to identify its ability to model dislocations in zirconium. In addition, this empirical potential is used to study the convergence of our results with the size of the simulation cell.

II. ATOMIC INTERACTION MODELING

Atomistic simulations have been performed both with an empirical interatomic potential and with ab initio calculations. The empirical potential that we used is the EAM potential developed by Mendelev and Ackland. 24 This potential is labeled #3 in Ref. 24. It is supposed to be well suited to model dislocations, as ab initio values 28 of the stacking fault energies in the basal and prism planes have been included in the fitting procedure. Using this potential, Khater and Bacon 25 showed that a 1/3 1 210 screw dislocation spontaneously dissociates in the prism plane and that a metastable configuration dissociated in the basal plane also exists.

The ab initio calculations are relying on the Density Functional Theory (DFT) in the Generalized Gradient Approximation (GGA) with the functional proposed by Perdew, Burke, and Ernzerhof (PBE) and the pseudopotential approximation. Two different ab initio codes are used, Siesta 26 and Pwscf. 27 In the Siesta code, 26 valence electrons are described by a localized basis set corresponding to a linear combination of pseudoatomic orbitals with 13 functions per atom. We used a norm conserving pseudopotential of Troulliers-Martins type with 4p electrons included as semicore. Electronic density of state is broaden with the Methfessel-Paxton function with a broadening of 0.3 eV and the integration is performed on a regular grid of 14 × 14 × 8 k-points for the primitive hcp unit cell and an equivalent density of k-points for the supercells used for the defect calculations. The charge density is represented on a real space grid with a spacing of 0.08 Å (Mesh cutoff: 450 Ry) that is reduced to 0.06 Å (800 Ry) for dislocation calculations. This approach, i.e. the basis and the pseudopotential, has already been used to study vacancy diffusion in zirconium 29 and comparison with plane waves DFT calculations has led to a reasonable agreement.

In the Pwscf code, 27 valence electrons are described with plane waves using a cutoff energy of 28 Ry. The pseudopotential is ultrasoft of Vanderbilt type with 4s and 4p electrons included as semicore. [START_REF]Pseudopotential file Zr.pbe-nsp-van[END_REF] The same kpoint grid and the same electronic broadening are used as with Siesta code.

To validate these different atomic interaction models, it is worth comparing their results to available experimental data for some bulk properties of Zr. All models lead to an equilibrium lattice parameter and a c/a ratio in good agreement with experimental data (Tab. I). In particular, a ratio lower than the ideal 8/3 ∼ 1.633 value is obtained in all cases.

We also compared the theoretical elastic constants with experimental data (Tab. I): a good agreement is also obtained. The computed elastic constants are the relaxed ones: 35 as the hcp lattice contains two atoms in its primitive unit cell, some internal degrees of freedom may exist when applying a homogeneous strain. One needs to allow for atomic relaxations when computing C 11 , C 12 or C 66 constants. [START_REF] Cousins | [END_REF] It is also possible to calculate inner elastic constants to characterize these internal degrees of freedom. These are also given in Tab. I using the notations introduced by Cousins. [START_REF] Cousins | [END_REF] Two of these inner elastic constants, e 11 and e 33 , are related to the phonon frequencies of the optical branches at the Γ point 31

ω i = 2 Ωe ii m , (1)
where Ω = a 2 c √ 3/4 is the atomic volume and m the atomic mass. The last inner elastic constants d 21 couples the internal degrees of freedom to the homogeneous strain. It leads to a contribution to the elastic constants 31

δC 12 = d 21 2 e 11 . (2)
If C 0 ij are the unrelaxed elastic constants, i.e. the elastic constants calculated by imposing a homogeneous strain to the hcp lattice without letting atoms relax from their initial positions, the true elastic constants are given by 31 C 11 = C 0 11 -δC 12 , C 12 = C 0 12 + δC 12 , and C 66 = C 0 66 -δC 12 , all other elastic constants being unchanged.

III. STACKING FAULT ENERGIES

Dislocation dissociation is controlled by the existence of a metastable stacking fault for the corresponding plane. According to the results obtained by Khater and Bacon 25 with Mendelev and Ackland EAM potential, 24 TABLE II. Stacking fault energies in the basal plane, γ b , and in the prism plane, γp, calculated with different atomic interaction models, including Vasp calculations of Domain et al.. 28 and of Udagawa et al.. 37 R is the ratio defined by Legrand, 5 and d eq b and d eq p are the dissociation lengths of a screw dislocation respectively in the basal (Eq. 5) and in the prism plane (Eq. 6).

EAM

Siesta Pwscf Vasp a 1/3 1 210 screw dislocation can dissociate either in a basal or in a prism plane. To characterize these eventual dissociations, we compute generalized stacking fault energies [START_REF] Vitek | Dislocations in Solids[END_REF]36 -or γ-surfaces -for both the basal and the prism planes.

A. Methodology γ-surfaces describe the energy variation when two parts of a crystal are rigidly shifted for different fault vectors lying in a given crystallographic plane. Atoms are allowed to relax in the direction perpendicular to the fault plane. We calculate these γ-surfaces for both the basal and prism planes using full periodic boundary conditions. To introduce only one fault in the simulation cell, the same shift as the one corresponding to the fault vector is applied to the periodicity vector perpendicular to the fault plane. No free surfaces is therefore introduced in the simulation cell, which allows a fast convergence of the fault energies with the number of stacked planes. A periodic stacking of at least 16 {0001} planes is used for the basal fault and 12 {10 10} planes for the prismatic fault. This corresponds to a distance between fault planes h 0001 = 8c ∼ 41 Å and h 10 10 = 6 √ 3a ∼ 34 Å. Increasing the number of planes in the stacking modifies the energies by less than 1 mJ m -2 . Generalized stacking fault energies are calculated on a regular grid of 10 × 10 fault vectors and are then interpolated with Fourier series.

(a) Pwscf -basal plane

3/2 1 1/2 0 a/3 [12 -10] 1 2/3 1/3 0 a/2 [101 -0] 0 0.25 0.5 0.75 γ (J m -2 ) (b) Siesta -basal plane 3/2 1 1/2 0 a/3 [12 -10] 1 2/3 1/3 0 a/2 [101 -0] 0 0.25 0.5 0.75 γ (J m -2 ) (c) EAM -basal plane 3/2 1 1/2 0 a/3 [12 -10] 1 2/3 1/3 0 a/

B. γ-surfaces

1. Basal plane γ-surfaces corresponding to the basal plane are shown in Fig. 1 for the three interaction models we used. In all cases, a local minimum can be found at 1/3[1 100] which corresponds to the intrinsic I 2 fault. [START_REF] Hull | Introduction to Dislocations[END_REF] This minimum does not vary when full atomic relaxations are allowed . This appears clearly in Fig. 2 where the fault energy predicted by the different interaction models are compared along the [1 100] direction. We will see latter that this has consequences on the stability of a screw dislocation dissociated in the basal plane. The γ-surface calculated with the EAM empirical potential has another minimum located in 2/3[1 100]. This is an artifact of the potential: one expects instead a maximum as this fault vector transforms the original BABABA stacking of the basal planes in a BABBCB stacking. Ab initio calculations confirm that this fault vector gives a maximum. Finally, it is worth pointing that both ab initio techniques give a very similar γsurface: the shapes are identical and the amplitudes do not differ by more than 10%.

Prism plane

γ-surfaces for the prism plane are shown in Fig. 3. Both ab initio methods lead to a similar γ-surface [Figs. 3(a) and (b)]. In particular, both Pwscf and Siesta predicts the existence of a minimum at halfway of the Burgers vector, i.e. in 1/6 [1 210]. Like for the basal fault, this minimum does not vary when full atomic relaxations are allowed. As can be seen on the projection of these γ-surface along the [1 210] direction (Fig. 4), this minimum is a little more pronounced with Pwscf than with Siesta. The same minimum was also present in the Vasp calculations of Domain et al., 28 but they obtained a lower value γ p of the fault energy in this point (Tab. II). On the other hand, Udagawa et al. 37 a value of about 200 mJ m -2 . The γ-surface calculated with the EAM potential is quite different [Fig. 3(c)] as the point located in 1/6 [1 210] is indeed a maximum and not a minimum like with ab initio calculations. A minimum is found for a point located in a/6 [1 210] + 0.14c [0001]. One therefore expects that a 1/3 1 210 {10 10} dislocation dissociates into two partial dislocations with a Burgers vector component orthogonal to the one of the perfect dislocation. In particular, a screw dislocation should dissociate in two partial dislocations with a small edge character. Khater and Bacon 25 showed that the empirical potential of Ackland et al. [START_REF] Ackland | [END_REF] suffers from the same artefact. As noted by Bacon and Vitek, 41 all central forces potentials stabilize indeed such a stacking fault a/6 [1 210] + αc [0001] with α = 0, a minimum also predicted by a simple hard sphere model. 42 This minimum either disappears or is located exactly in a/6 [1 210] (α = 0) only when the angular dependence of the atomic interactions is taken into account. The value γ p of the stacking fault energy obtained with Mendelev and Ackland EAM potential 24 for this minimum is much lower than our ab initio value (Tab. II). This is quite normal as Mendelev and Ackland used the ab initio value of Domain et al. to fit their potential.

C. Dislocation dissociation

Before using atomistic simulations to obtain dislocation core structures and their associated Peierls stress, it is worth looking what can be learned from these γ-surface calculations. Legrand 5 proposed a criterion based on the elastic constants and the metastable stacking fault energies to determine if glide occurs in the base or prism plane in a hcp crystal. According to this criterion, prismatic glide is favored if the ratio R = C 66 γ b /C 44 γ p is larger than 1. Tab II shows that this is the case for the EAM potential and the Pwscf calculation, as well as for the Vasp calculation of Domain et al. 28 and of Udagawa et al.. 37 On the other hand, Siesta leads to a value too close to 1 to be able to decide between basal and prismatic glide.

One can also use dislocation elasticity theory [START_REF] Hirth | Theory of Dislocations[END_REF] to compute the dissociation distance of a dislocation both in the basal and prism planes. According to elasticity theory, the energy variation caused by a dissociation of length d is

∆E diss (d) = -b (1) i K ij b (2) j ln d r c + γd, (3)
where b (1) and b (2) are the Burgers vectors of each partial dislocation, γ the corresponding stacking fault energy, K the Stroh matrix [START_REF] Stroh | [END_REF]45 controlling dislocation elastic energy, and r c the core radius. Minimization of this energy leads to the equilibrium dissociation length

d eq = b (1) i K ij b (2) j γ . ( 4
)
When the hcp crystal is oriented with the x, y, and z axis respectively along the [10 10], [0001], and [1 210] directions, for a dislocation lying along the z direction, the K matrix is diagonal with its components given by 46-4849

K 11 = 1 2π C11 + C 13 C 44 C11 -C 13 C 33 C11 + C 13 + 2C 44 , K 22 = C 33 C 11 K 11 , K 33 = 1 2π
d eq b = (3K 33 -K 11 ) a 2 12γ b . (5)
According to the minimum of the γ-surface in the prism plane (Fig. 3), a 1/3[1 210] dislocation can dissociate in this plane in two partial dislocations with Burgers vectors 1/6[1 210] ± αc/a[0001]. The parameter α controls the position of the stacking fault minimum along the [0001] direction, i.e. α = 0 for Pwscf and Siesta, and α = 0.14 for the EAM potential. The dissociation length of a screw dislocation in the prism plane is then

d eq p = K 33 a 2 -4α 2 K 22 c 2 4γ p . (6)
The dissociation lengths d eq b and d eq p calculated from the elastic constants and the stacking fault energies are given in table II. Elastic constants predicted by the atomic interaction models are used in each case. For all energy models, one expects a larger dissociation in the prism than in the basal plane. We will compare in the following section these dissociation lengths predicted by elasticity theory with the ones observed in our atomistic simulations of the dislocation core structure. Our atomistic simulations of the core structure of a screw dislocation are based on full periodic boundary conditions. [START_REF] Bulatov | Computer Simulations of Dislocations[END_REF] This requires introducing a dislocation dipole in the simulation cell. Two periodic arrangements of the dislocations have been used (Fig. 5). In the O arrangement, dislocations with opposite Burgers vectors are located on the same prism and basal planes, i.e. the two foreseen glide planes. On the other hand, only dislocation with the same Burgers vectors can be found on a given prism or basal plane in the S arrangement.

Periodicity Both dislocation arrays are quadrupolar: the vector joining the two dislocations composing the primitive dipole is D = 1/2( U 1 + U 2 ). Because of the centrosymmetry of this arrangement and the symmetry properties of the Volterra elastic field, this ensures that the stress created by other dislocations is minimal at each dislocation position. The cut vector A defining the dislocation dipole is obtained by a π/2 rotation of D.

The dislocation dipole is introduced in the simulation cells by applying to all atoms the elastic displacement predicted by anisotropic elasticity theory [START_REF] Stroh | [END_REF]45 taking full account of the periodic boundary conditions. [START_REF] Cai | [END_REF] A homogeneous strain is also applied to the simulation cell so as to cancel the plastic strain introduced by the dislocation dipole and minimize the elastic energy. This strain is given by

ε 0 ij = - b i A j + b j A i 2S ,
where

S = |( U 1 ∧ U 2 )
• e z | is the surface of the simulation cell perpendicular to the dislocation lines. This homogeneous strain adds some tilt components to the periodicity vectors. Atoms are then relaxed until all components of the atomic forces are smaller than 5 meV Å-1 for Siesta, 2 meV Å-1 for Pwscf, and 0.1 meV Å-1 for the EAM potential. Starting from perfect dislocations, atom relaxation leads to dislocations spread in the prism plane, whatever the interaction model (EAM, Siesta, or Pwscf) and whatever the simulation cell used. This can be clearly seen by plotting differential displacement maps as introduced by Vitek. 52 These maps (Fig. 6) show that the strain created by the screw dislocation spreads out in the (10 10) prism plane and that displacements outside this plane are much smaller.

To characterize this spreading in the (10 10) prism plane, we extract from our atomistic simulations the disregistry D(x) created by the dislocation. This is de- fined as the displacement difference between the atoms in the plane just above and those just below the dislocation glide plane. The derivative of this function, ρ(x) = ∂D/∂x corresponds to the dislocation density. Fig. 7 shows the disregistry obtained for the three interaction models. In all cases, the b discontinuity created by the screw dislocation does not show a sharp interface, but spreads on a distance ∼ 10 Å.

Peierls [START_REF] Peierls | Proc. Phys. Soc[END_REF] and Nabarro [START_REF] Nabarro | Proc. Phys. Soc[END_REF] built a model that leads to a simple expression of the disregistry. The analytical expression they obtained [START_REF] Hirth | Theory of Dislocations[END_REF] can be extended to a dissociated dislocation. As suggested by the prismatic γ-surface (Fig. 3), we assume that the screw dislocation dissociates in two equivalent partial dislocation separated by a distance d. Based on the Peierls-Nabarro model, we write the disregistry created by a single dislocation

D dislo (x) = b 2π arctan x -x 0 -d/2 ζ + arctan x -x 0 + d/2 ζ + π/2 , (7)
where x 0 is the dislocation position, d its dissociation length, and ζ the spreading of each partial dislocation. We need then to take into account that we do not have only one dislocation on a given prism plane but a periodic array (Fig. 5). The disregistry created by an array of period L is

D L (x) = ∞ n=-∞ D dislo (x -nL) = b 2π    arctan   tan π L [x -x 0 -d/2] tanh πζ L   +π x -x 0 -d/2 L + 1 2 + arctan   tan π L [x -x 0 + d/2] tanh πζ L   +π x -x 0 + d/2 L + 1 2 , ( 8
)
where ⌊•⌋ is the floor function. For the O arrangement (Fig. 5a), the disregistry in the prism plane should be given by D(x) = D L (x) -D L (x -L/2) with L = 2mc, whereas it should be D(x) = D L (x) with L = mc for the S arrangement (Fig. 5b).

We fit this analytical expression of the dislocation disregistry to the data coming from the atomistic simulations. Fig. 7 shows a good agreement between atomistic simulations and the model, using only three fitting parameters: the dislocation position x 0 , the dissociation length d and the spreading ζ. This procedure therefore allows us to determine the location of the dislocation center. For all interaction models, we find that this center lies in between two (0001) atomic planes. One can see in Fig. 6 that this position corresponds to a local symmetry axis of the differential displacement map. This is different from the result obtained by Ghazisaeidi and Trinkle 13 in Ti where the center of the screw dislocation was found to lie exactly in one (0001) atomic plane, a position that corresponds in Zr to the saddle point between two Peierls valleys, as it will be shown below.

The dissociation length of the screw dislocation obtained through this fitting procedures are shown in Fig. 8 for both periodic arrangements. We observe variations with the dislocation periodic arrangement used in the simulation, as well as with the size of the simulation cell. The results are nevertheless close to the predictions of elasticity theory ( §III C) for all the three interaction models. We could even see with the EAM potential that the dissociation lengths extracted from atomistic simulations converge to the value given by elasticity theory for large enough unit cells. It is thus relevant to describe the screw dislocation as dissociated in two partial dislocations linked by a stacking fault, although the dissociation length remains small.

C. Core energy

We obtain the dislocation core energy by subtracting the elastic energy from the excess energy given by the atomistic simulations. This excess energy is the energy difference between the simulation cell containing the dislocation dipole and the same cell without any defect. The elastic energy calculation takes into account the elastic anisotropy [START_REF] Stroh | [END_REF]45 and the effect of periodic boundary conditions. [START_REF] Cai | [END_REF] The obtained core energy are shown in Fig. 9. Like for the dissociation length, some variations with the size of the simulation cell and the periodic arrangement can be observed. We did not manage to link both quantities, i.e. the core energy and the dissociation length. They are not simply related by the expression of the energy variation with the dissociation length (Eq. 3), and the change of the elastic interaction between dislocations caused by their dissociation could not fully explain the variation of the core energy either. As the spread-ing of partial dislocations also depends on the size of the simulation cells, the variation of the core energy probably have a more complex origin than simply a variation of the dissociation length. We could use, with the EAM potential, a simulation cell large enough to obtain a converged value of the core energy, 163 meV Å-1 for a core radius r c = b. It is worth pointing that this value does not depend on the periodic arrangement used in the simulation (Fig. 9). This can be achieved thanks to a proper account of the core traction contribution to the elastic energy. [START_REF] Clouet | [END_REF] A difference of ∼ 10 meV Å-1 would have been observed between the S and O periodic arrangement without this contribution. Ab initio calculations lead to a dislocation core energy of 145 ± 5 meV Å-1 for Siesta and 125 ± 5 meV Å-1 for Pwscf (r c = b in both cases).

D. Dissociation in the basal plane

Basal slip is observed experimentally only at high temperature (above 850 K) and for a higher resolved shear stress than the one needed to activate prismatic slip. 56 At lower temperatures, no basal slip could be observed, 3,57 even when the monocrystal was oriented so as to favor basal slip. Our atomistic simulations lead to a screw dislocation configuration dissociated in the prism plane, which clearly cannot glide easily in the basal plane. It is worth looking if another configuration, which could glide in this basal plane, also exists. Using the same EAM potential, Khater and Bacon 25 showed that a screw dislocation can also dissociate in the basal plane. The basal configuration is obtained by introducing in the same basal plane two partial dislocations of Burgers vector 1/3[1 100] and 1/3[0 110], in agreement with the location of the minimum on the basal γsurface (Fig. 1). After relaxation of the atomic positions, the dislocation remains dissociated in the basal plane, as can be seen from the corresponding differential displacement map (Fig. 10). A fit of the screw component of the disregistry created by the dislocation in the basal plane leads to a dissociation length d = 6.0 Å, a higher value than predicted by elasticity theory (Tab. II). This configuration is metastable. It has indeed an energy higher by 62 meV Å-1 than the configuration dissociated in the prism plane.

We check if ab initio calculations also lead to such a metastable configuration. Starting from a screw dislocation dissociated in two partial dislocations in the basal plane, both Siesta and Pwscf lead after relaxation of the atomic positions to the stable configuration dissociated in the prism plane. This is true both with the S (n = 5) and the O (n = 4) periodic arrangements. These ab initio calculations show then that such a dissociation of the screw dislocation in the basal plane is unstable. The metastable configuration observed with the EAM potential appears to be an artifact of the empirical potential. This is not specific to the Mendelev and Ackland potential 24 as a configuration dissociated in the basal plane is stabilized by any central forces potential. 25,41,58 V. PEIERLS BARRIER Before calculating ab initio the Peierls barrier of the screw dislocation, we use the EAM potential to assess the validity of the method and to check the convergence of the Peierls barrier with the size of the simulation cell.

A. Methodology

We determine the Peierls barrier for the screw dislocation gliding in the prism plane. This is done using a constrained minimization between two adjacent equilibrium configurations of the dislocations. We move both dislocations in the same direction by one Peierls distance λ P = c/2 between the initial and final states, so as to keep constant the distance between dislocations.

Two different algorithms are used to perform the constrained minimization, the simple drag method and the more robust nudged elastic band (NEB) method. [START_REF] Henkelman | Progress in Theoretical Chemistry and Physics[END_REF] Intermediate configurations are built by linearly interpolating the atomic coordinates between the initial and final states. We define the corresponding reaction coor-

dinate ζ = ( X -X I ) • ( X F -X I )/ X F -X I 2
, where X, X I , and X F are the 3N vectors defining atomic positions for respectively the intermediate, initial, and final configurations. In the drag method, the minimization is performed on all atomic coordinates with the constraint that ζ remains fixed for each of the nine intermediate images. Nine intermediate images are also used in NEB method with a spring constant k = 0.1 eV Å-2 .

So as to obtain the variation along the path of the dislocation energy with its position, we need to determine the dislocation position x Dislo (ζ) for each intermediate image, once it has been relaxed. This is done thanks to a fit of the disregistry in the prism plane, as described in §IV B. This allows us to check that both dislocations in the simulation cell are moving in a coordinated way: the distance between them remains fixed. As a consequence, there is no variation of the elastic energy along the path and the energy variation ∆E P (ζ) obtained by the constrained minimization corresponds to a variation of the core energy, i.e. the Peierls energy. We deduce the Peierls energy ∆E P (x Dislo ) by eliminating the reaction coordinate ζ between ∆E P (ζ) and x Dislo (ζ). Fig. 11 illustrates the whole procedure for a given dislocation periodic arrangement using both constrained minimization techniques. The variation ∆E P (ζ) slightly differs between both techniques: for a given reaction coordinate ζ, the drag method leads to a state of lower energy than NEB method. This corresponds to a small difference in the dislocation position: this position deviates more with drag than with NEB method from a linear variation. Nevertheless, one obtains at the end the same Peierls barrier ∆E P (x Dislo ), whatever the method used. These differences observed for the functions ∆E P (ζ) and x Dislo (ζ) between drag and NEB methods increase with the size of the simulation cell, i.e. with the number of degrees of liberty. For large simulation cells (containing more than 3600 atoms in the S periodic arrangement for instance), the drag method sometimes fails to find a continuous path between the initial and final states: one has to use the NEB method in theses cases. Only much smaller simulation cells can be studied ab initio. For these sizes, the drag and the NEB methods always lead to the same result. We will therefore only use the drag method in the ab initio calculations, as it costs much less CPU time. Finally, we interpolate with Fourier series the periodic functions ∆E P (ζ) and x Dislo (ζ) -λ P ζ. This leads to a smooth function ∆E P (x Dislo ) that can be derived. The Peierls stress σ P is deduced from the maximal slope of this function,

σ P = 1 b Max ∂∆E P ∂x Dislo . (9)
We obtain a Peierls stress σ P = 24 ± 1 MPa for the EAM potential. Khater and Bacon 25 determined, for the same empirical potential, a Peierls stress σ P = 22 MPa using molecular statics simulations under applied stress. As the agreement between both methods is good, we see that the Peierls stress can be defined either from the slope of the Peierls barrier or from the minimal applied stress under which the dislocation glides indefinitely.

We now examine, still with the EAM potential, how this Peierls barrier varies when the size of the simulation cell decreases up to reaching a number of atoms that can be handled in ab initio calculations (Fig. 12). Both the S and O periodic arrangements give the same Peierls barrier, and hence the same Peierls stress, for a large enough simulation cell ( 1000 atoms). The Peierls barrier increases when the size of the simulation cell decreases with the S periodic arrangement, whereas it decreases with the O periodic arrangement. As a consequence, the S and O periodic arrangement should respectively lead to an upper and lower limits of the Peierls stress for small simulation cells. We will therefore use the S periodic arrangement to calculate ab initio this Peierls barrier. This will allow us to confirm that the Peierls stress is as low as indicated by experiments and this EAM potential. Moreover, it is worth pointing that the expected Peierls barrier should be small: ∆E P = 0.4 meV Å-1 at the saddle point according to the EAM potential. This corresponds to a difference of energy 2b∆E P = 2.6 meV for a simulation cell of minimal height containing a dislocation dipole. Such a small energy difference may be problematic because of the precision of ab initio calculations. Looking for an upper limit of this value is easier as it minimizes the problems associated with this precision.

Finally, it is worth pointing that the dissociation length varies during the dislocation migration, and this variation is more pronounced (∼ 10%) for the smallest simulation cells. But, like for the core energy, we did not manage to relate the size dependence of the Peierls barrier to this variation of the dissociation length.

B. Ab initio barriers

The Peierls barriers obtained by ab initio calculations are shown on Fig. 13(a) for Pwscf and Fig. 13(b) for Siesta. In both cases, the height of the barrier decreases when the number of atoms increases, in agreement with what has been observed with the EAM potential for the same S periodic arrangement. For a given number of atoms, the ab initio barriers are a little bit lower than the ones obtained with the EAM potential. Results obtained with Siesta are noisy: the energy barrier that we want to calculate is so small that it needs a really strict convergence criterion on atomic forces for the relaxation. We did not manage to reach such a precision with Siesta. On the other hand, the barriers obtained with Pwscf are smooth. We can therefore interpolate the ab initio results with Fourier series, and estimate then the Peierls stress (Eq. 9). This leads to σ P = 36 MPa for the simulation cell containing 100 atoms and σ P = 21 MPa for 144 atoms. Considering that these values, obtained in small simulation cells, are upper limits, ab initio calculations predict a Peierls stress smaller than 21 MPa for a screw dislocation in zirconium gliding in a prism plane.

The comparison of this ab initio estimate of the Peierls stress with experimental data 2,3,19,57 is quite challenging. As pointed out in the introduction, the yield stress of zirconium strongly depends on its oxygen content. No experimental data exists for a purity better than 0.07% O (in atomic fraction). Moreover, all measure- ments have been performed at temperatures higher than 77 K. The determination of pure zirconium flow stress at 0 K therefore needs some extrapolation of experimental data. Without a clear understanding of the mechanisms responsible of zirconium hardening by oxygen impurities, such an extrapolation is quite hazardous. Nevertheless, a graphical comparison (Fig. 14) of experimental data with our ab initio Peierls stress shows a reasonable agreement, keeping in mind that the value 21 MPa has to be considered as an upper limit.

VI. CONCLUSIONS

Using two ab initio approaches (Pwscf and Siesta), both in the DFT-GGA approximation, we have shown that a 1/3 1 210 screw dislocation in hcp zirconium dissociates in two partial dislocations with a pure screw character. This is in agreement with the minimum in 1/6 1 210 found for the generalized stacking fault energy in the prism plane. We could extract the dissociation length from our atomistic simulations. Although this dissociation length is small (d ∼6 Å for Pwscf and d ∼ 4 Å for Siesta), it is in reasonable agreement with the one predicted by elasticity theory.

The EAM potential of Mendelev and Ackland 24 leads to the same structure of the dislocation dissociated in the prism plane. The metastable configuration dissociated in the basal plane predicted by this potential is not stable in ab initio calculations, both with Pwscf and Siesta. This configuration is therefore an artifact of this potential, probably induced by the deep minimum in 1/3 1 100 found with this empirical potential for the generalized stacking fault in the basal plane. This minimum is much more shallow in ab initio calculations.

We could also obtain an ab initio estimate of the Peierls stress of the screw dislocation gliding in the prism plane. Calculations with Pwscf lead to an upper limit of 21 MPa for this Peierls stress. This small value shows that screw dislocations can glide quite easily in pure zirconium, thus confirming what had been obtained with Mendelev and Ackland EAM potential. Such a small Peierls stress is in agreement with experimental data, once the hardening of oxygen impurities has been considered.

FIG. 1 .

 1 FIG. 1. Generalized stacking fault energy in the basal plane calculated with (a) Pwscf, (b) Siesta and (c) EAM potential. The arrows indicate Burgers vectors of the partial dislocations corresponding to a dissociation in the basal plane. The dashed line is the [1 100] direction used in Fig. 2. Contour lines are drawn at the base every 50 mJ m -2 .

FIG. 2 .

 2 FIG.2. Generalized staking fault energy in the basal plane along the [1 100] direction (cf. Fig.1) calculated with Pwscf, Siesta and EAM potential.

  instead of being constrained to the direction perpendicular to the fault plane. All methods give a close value for the fault energy γ b in this minimum (Tab. II). A good agreement is also obtained with the values calculated by Domain et al.28 and by Udagawa et al.37 using Vasp ab initio code. The depth of this local minimum is more pronounced with the empirical EAM potential [Fig. 1(c)] than in ab initio calculations [Fig. 1(a) and (b)]

FIG. 3 .

 3 FIG. 3. Generalized stacking fault energy in the prism plane calculated with (a) Pwscf, (b) Siesta and (c) EAM potential. The arrows indicate Burgers vectors of the partial dislocations corresponding to a dissociation in the prism plane. Contour lines are drawn at the base every 75 mJ m -2 .

FIG. 4 .

 4 FIG. 4. Generalized staking fault energy in the prism plane along the [1 210] direction calculated with Pwscf, Siesta and EAM potential.

1 2 C

 2 44 (C 11 -C 12 ), where C11 = √ C 11 C 33 . The γ-surface of the basal plane (Fig. 1) indicates a possible dissociation 1/3[1 210] → 1/3[1 100] + 1/3[0 110]. The dissociation length in the basal plane is then, for a 1/3[1 210] screw dislocation,

FIG. 5 .

 5 FIG. 5. Screw dislocation periodic arrangements used for atomistic simulations. U1 and U2 are the periodicity vectors of the arrangement, and A the cut vector of the dislocation dipole. The thin vertical line corresponds to the prism glide plane.

  vectors of the O arrangement are, before introducing the dislocations, U 1 = na 1 2 [10 10]-mc [0001], U 2 = na 1 2 [10 10] + mc [0001], and U 3 = a 1 3 [1 210]. The integers n and m are taken equal to keep an aspect ratio close to a square. This simulation cell has been used for ab initio calculations with n = 4 (128 atoms), n = 5 (160 atoms), and n = 6 (288 atoms). For the S arrangement, U 1 = na 1 2 [10 10], U 2 = mc [0001], and U 3 = a 1 3 [1 210]. Ab initio calculations have been performed with n = m and varying between n = 5 (100 atoms) and n = 8 (256 atoms).

FIG. 6 .

 6 FIG.6. Differential displacement maps around one of the two 1/3[1 210] screw dislocations composing the dipole for the S periodic arrangement with n = m = 7 (196 atoms). Atoms are sketched by circles with a color depending on the (1 210) plane to which they belong. The arrow between two atomic columns is proportional to the [1 210] component of the differential displacement between the two atoms. Displacements smaller than 0.1b are not shown. Crosses × correspond to the positions of the two partial dislocations, and + to their middle, i.e. the position of the total dislocation.

FIG. 7 .

 7 FIG. 7. Disregistry D(x) created by the screw dislocation in the prism plane and corresponding dislocation density ρ(x) = ∂D(x)/∂x for the S periodic arrangement with 196 atoms (n = m = 7). Symbols correspond to atomistic simulations and lines to the fit of the Peierls-Nabarro model to these data. For clarity, disregistries D(x) have been shifted by 0.25 between the different interaction models.

FIG. 8 .FIG. 9 .

 89 FIG.8. Dissociation length of the screw dislocation. Symbols correspond to data extracted from atomistic simulations through the fit of the disregistry (Eq. 8) for both the O and S periodic arrangements and for different sizes of the simulation cell. The solid lines are the predictions of elasticity theory based on the stacking faults (Tab. II). On the right vertical axis, the experimental c lattice parameter has been used to normalized the dissociation length by the distance λP = c/2 between Peierls valleys.

FIG. 10 .

 10 FIG.10. Differential displacement map of the metastable configuration of a 1/3[1 210] screw dislocation dissociated in the basal plane, as obtained with the EAM potential for the O periodic arrangement with n = 6 (288 atoms). Crosses × correspond to the positions of the two partial dislocations, and + to their middle, i.e. the position of the total dislocation.

FIG. 11 .

 11 FIG.11. Peierls barrier of a screw dislocation calculated with the EAM potential for the S periodic arrangement with 1600 atoms (n = m = 20). Two different methods for finding the minimum energy path have been used: the NEB method and a simple constrained minimization (drag). Symbols correspond to the results of atomistic simulations and lines to their interpolation with Fourier series.

FIG. 12 .

 12 FIG.12. Variation of the Peierls barrier with the size of the simulation cell calculated with the EAM potential for both periodic arrangements.

FIG. 13 .FIG. 14 .

 1314 FIG. 13. Peierls barrier for a screw dislocation gliding in its prism plane calculated ab initio with (a) Pwscf and (b) Siesta for the S periodic arrangement. Symbols correspond to ab initio results and lines to their interpolation by Fourier series.

TABLE I .

 I Bulk properties of hcp Zr calculated with different atomic interaction models and compared to experimental data: lattice parameter a, c/a ratio of the hexagonal lattice, relaxed elastic constants Cij , inner elastic constants 31 eij and dij, phonon frequencies ω1 and ω3 of the optical branches at the Γ point (Eq. 1), inner elasticity contribution to elastic constant δC12 (Eq. 2).

		Expt.	EAM Siesta Pwscf
	a ( Å)	3.232 32	3.234	3.237	3.230
	c/a	1.603 32	1.598	1.613	1.601
	C11 (GPa)	155.4 a	142.	140.	140.
	C33 (GPa)	172.5 a	168.	168.	168.
	C12 (GPa)	67.2 a	75.	86.	70.
	C13 (GPa)	64.6 a	76.	68.	65.
	C44 (GPa)	36.3 a	44.	24.	26.
	C66 (GPa)	44.1 a	33.5	27.	35.
	e11 (meV Å-5 )		27.0	17.0	18.6
	e33 (meV Å-5 )		118.	122.	101.
	d21 (meV Å-4 )		30.0	38.9	36.6
	ω1 (THz)	2.66 ± 0.02 33	2.60	2.08	2.16
	ω3 (THz)	4.23 ± 0.15 33	5.43	5.57	5.03
	δC12 (GPa)		5.33	14.3	11.5

a Experimental elastic constants

34

have been measured at 4 K.

ACKNOWLEDGMENTS

This work was performed using HPC resources from GENCI-CINES and GENCI-CCRT (Grant Nos 2011-096020 and 2012-096847).


C Simenel
email: cedric.simenel@anu.edu.au

A Wakhle

B Avez

D J Hinde

R Du Rietz

M Dasgupta

M Evers

C J Lin

D H Luong

Effects of Nuclear Structure on Quasi-fission

The quasi-fission mechanism hinders fusion of heavy systems because of a mass flow between the reactants, leading to a re-separation of more symmetric fragments in the exit channel. A good understanding of the competition between fusion and quasi-fission mechanisms is expected to be of great help to optimize the formation and study of heavy and superheavy nuclei. Quantum microscopic models, such as the time-dependent Hartree-Fock approach, allow for a treatment of all degrees of freedom associated to the dynamics of each nucleon. This provides a description of the complex reaction mechanisms, such as quasi-fission, with no parameter adjusted on reaction mechanisms. In particular, the role of the deformation and orientation of a heavy target, as well as the entrance channel magicity and isospin are investigated with theoretical and experimental approaches.

I. INTRODUCTION

The formation of the heaviest nuclei usually involves fusion-evaporation reactions [1][2][3][4]. The latter are strongly hindered in the case of heavy ion reactions by two competing mechanisms: (i) the quasi-fission (QF) process, and (ii) the statistical fission of the compound nucleus (CN).

Quasi-fission occurs in the early stage of the collision [5][6][7], when the two reactants form a di-nuclear system, that is, two fragments linked by a neck. An important nucleon transfer usually occurs from the heavy fragment toward the light one. The two fragments then re-separate with more mass symmetry than the entrance channel, without forming a compound nucleus.

Typical QF times are shorter than 10 -20 s [5][6][7][8]. These times have to be compared with fusion-fission times which can be longer than 10 -16 s [9]. This shows that the two mechanisms are of very different nature. In fact, the QF process is a dynamical mechanism depending on the characteristics of the entrance channel, while CN fission is pure-ly statistical and is determined by temperature and angular momentum only.

The QF mechanism is responsible for the fusion hindrance observed in heavy systems [10]. In these reactions, an additional energy above the Coulomb barrier, sometimes called "extra-push" energy [11], is needed for the system to fuse and form a CN. Note that lighter systems may also exhibit quasifission, although with a smaller probability. For instance, quasi-fission has been observed in 16 O, 32 S+ 238 U [12][13][14], and in 32 S+ 208 Pb [14].

The QF process is known to be affected by nuclear deformation and orientation at energies close to the fusion barrier [12,[15][16][17][18][19][20]. The role of isospin has also been investigated theoretically [21]. Recently, the influence of entrancechannel magicity and isospin on quasi-fission has been investigated [22]. Magic shells are indeed expected to generate "cold valleys" in the potential energy surface, favouring the formation of a compact CN [23][24][25]. In addition, magic nuclei are difficult to excite, reducing energy dissipation, and, then, allowing more compact di-nuclear systems [26,27].

Many experimental data on QF are now available for comparison with theoretical models in order to test their predictive power. This is indeed crucial to have reliable theoretical models in order to drive future experiments on heavy elements formation. Macroscopic approaches have been thoroughly used in the past [28,29]. In addition, the recent increase of computational power allowed micros-copic descriptions of nuclear dynamics. For instance, systems as heavy as actinide collisions have been studied with microscopic approaches [30][31][32][33][34].

In particular, the time-dependent Hartree-Fock (TDHF) theory [35] is particularly well suited at low energies where a proper treatment of the interplay between reaction mechanisms and nuclear structure is crucial. Indeed, TDHF calculations treat both dynamics and ground-state structure on the same footing, i.e., with the same energy density functional (EDF) as the only phenomenological input. As a result, the TDHF approach has been successful in describing several reaction mechanisms, such as fusion, nucleon transfer, and deep-inelastic collisions (see Ref. [36] for a review).

We first illustrate the fusion hindrance with TDHF calculations of fusion thresholds in heavy systems. Then we present an experimental study of the role of magicity and isospin on the quasi-fission mechanism. Finally, we discuss results of a recent study of quasi-fission with TDHF calculations.

II. FUSION HINDRANCE IN HEAVY SYSTEMS

Recent TDHF calculations have been performed to investigate collisions of heavy systems leading to quasi-fission and formation of super-heavy compound nuclei [31,33,34,[37][38][START_REF] Simenel | proceedings of the 11th International Conference on Nucleus-Nucleus Collisions (NN2012)[END_REF][START_REF] Guo | proceeding of the NSRT12 conference[END_REF]. In particular, fusion hindrance in heavy systems was predicted [START_REF] Simenel | proceedings of the 11th International Conference on Nucleus-Nucleus Collisions (NN2012)[END_REF][START_REF] Guo | proceeding of the NSRT12 conference[END_REF]. The additional energy needed for fusion to occur has been found to be of the same order of magnitude as the extra-push energy [START_REF] Simenel | proceedings of the 11th International Conference on Nucleus-Nucleus Collisions (NN2012)[END_REF] determined with the phenomenological approach of Swiatecki [11].

Fig. 1 shows the evolution of contact time as a function of energy for different central collisions obtained with the tdhf3d code [START_REF] Kim | [END_REF]. Here, the contact time is determined arbitrarily as the time the system spend with a distance between the centers of mass of the fragments smaller than 15 fm. The most asymmetric reaction, 48 Ti+ 208 Pb with Z 1 Z 2 = 1804, produces long contact times at energies above the proximity barrier [42] which could be assimilated to fusion. In this case, no (or little) extra-push energy is needed to fuse.

Increasing the mass of the projectile, as in the 70 Zn+ 208 Pb reaction with Z 1 Z 2 = 2460, changes drastically this behaviour. In this case, we observe a slow increase of the contact time above the barrier, with a maximum of ∼ 4 zs at ∼ 1.15B prox. . These calculations do not show any fusion for this system. Instead, the di-nuclear system always encounters quasi-fission.

The 90 Zr+ 124 Sn reaction is more symmetric, but with an intermediate charge product Z 1 Z 2 = 2000. At and up to 10% above the barrier, the collision time increases slowly, similarly to the 70 Zn+ 208 Pb case, indicating that quasi-fis-sion occurs within this energy range. At higher energies, the TDHF calculations predict long contact times (greater than 20 zs) which may be associated to fusion reactions. As a result, an extrapush E T DHF X ∼ 23 MeV is predicted for this system, which is of the same order of magnitude, but slightly larger, than the value given by the extra-push model E S wi. X ∼ 15 MeV [11]. Similar agreements have been obtained for other systems in Ref. [START_REF] Simenel | proceedings of the 11th International Conference on Nucleus-Nucleus Collisions (NN2012)[END_REF].

Detailed investigations on the quasi-fission mechanism are mandatory to understand fusion hindrance in heavy systems. In particular, a good understanding of how these mechanisms are affected by entrance channel properties is of utmost importance to optimise the formation of super-heavy elements.

III. QUASI-FISSION: INTERPLAY BETWEEN SHELLS AND ISOSPIN

A series of experiments to investigate the quasi-fission process have been performed recently at the Australian National University in Canberra [8,22,[START_REF] Hinde | AIP Conf. Proc[END_REF][START_REF] Lin | [END_REF][45]. In Ref. [22], we showed that shell effects in the entrance channel may affect the quasi-fission process. This was done by measuring massangle distributions (MAD) of the fragments in several reactions at sub-barrier energies. The experimental setup is shown in Fig. 2. The beams were produced by the 14UD electrostatic accelerator. Two Multi-Wire Proportional Chambers (MWPC) were used to measure time and positions of both fission fragments in coincidence (see Ref. [22] for more details on the geometry of the setup). Time of flight (ToF) and positions were converted into fragment masses and angles using two-body kinematics.

Fig. 3 shows the resulting MAD (up) and the projections on the mass-ratio axis (bottom). The mass ratio is defined as

M R = m 2 /(m 1 + m 2 )
where m 2 (resp. m 1 ) is the mass of the fragment in the back (front) detector. Fission and quasi-fission fragments are located between the bands at extreme M R corresponding to (quasi-)elastic and deep-inelastic events. Fig. 4 sketches the "trajectory" of the fragments in the MAD. In particular, short scission times induce correlations, i.e. the distribution forms a finite angle with the M R = 0.5 axis. Note that this simple picture neglects possible shell effects in the exit channel which may enhance the production of magic nuclei such as 208 Pb and "delay" the mass drift toward symmetry. In addition, long-time fission of super-heavy systems may occur via asymmetric channels [46,47], i.e., fusion-fission does not necessarily populate the region of the MAD around the M R = 0.5 axis.

Correlations between mass and angle of the fission fragments can be seen on some spectra of Fig. 3 (see,e.g., the case of 44 Ca+ 208 Pb which has the highest statistics). Such correlations increase the width of the mass ratio distribution (lower panels). We then use the width of the fission fragment mass distribution to quantify the amount of "fast reseparation" which could be associated to quasi-fission events. In particular, the 16 O+ 238 U reaction, which is known to exhibit only a small amount of quasi-fission [12], gives a width which can be considered as an upper limit for pure fusionfission events. Fig. 5 shows the width of the fission-like fragment mass distributions as a function of the number of magic numbers in the entrance channel. We see that, apart for the 40 Ca+ 208 Pb case which is discussed below, there is a clear link between the two quantities, i.e., the more the magicity, the smaller the width. This is interpreted as a hindrance of the transfer toward symmetry process due to shell-effects in the entrance channel.

The effect of entrance channel magicity on fusion as been investigated in Ref. [48]. Fig. 6 shows a comparison of fusionevaporation cross-sections for reactions forming 220 Th * [48]. We can see that fusion-evaporation cross-sec-tions decrease when the mass asymmetry increases. On the contrary, these cross-sections increase when shell effects are present in the entrance channel. Indeed, 48 Ca is doubly magic and 124 Sn has a magic proton number. Assuming that capture cross-sections, i.e., the sum of QF and fusion cross-sections, are not sensitive to shell-effects, we conclude that shell effects in the entrance channel hinder quasi-fission while it favours fusion of the fragments.

We see in Fig. 5 that the 40 Ca+ 208 Pb does not lie on the global trend. In Ref. [22], we interpreted this apparent discrepancy as an effect of charge equilibration occurring in the early stage of the collision. This equilibration process is indeed very rapid as shown by the TDHF calculations reported in Fig. 7. The latter are performed at the experimental energy E c.m. = 179.1 MeV, and varying L to obtain different collision times. ∆(N/Z) f is the difference between the N/Z ratio of the fragments after the collision. The initial value of ∆(N/Z) is 0.54 for 40 Ca+ 208 Pb and 0.14 for 48 Ca+ 208 Pb. We observe a (partial) equilibration (i.e., a reduction of ∆(N/Z) f ) within ∼ 2 zs in the 40 Ca+ 208 Pb system. This time is short as compared to the quasi-fission time which, according to simple simulations based on the model of Ref. [8], is longer than 10 zs. This means that the collision partners change their N and Z at contact and, as far as quasi-fission is concerned, behave like a non-magical system. This effect is not present in 48 Ca+ 208 Pb which is less N/Z asymmetric and does not encounter a charge equilibration (see the red dashed line in Fig. 7). Fusion is then more probable in this reaction which preserves its magic nature in the di-nuclear system. We conclude that shell effects in the entrance channel hinder quasi-fission (and then, favour fusion) only for systems with small N/Z asymmetry.

IV. QUASI-FISSION WITHIN THE TDHF APPROACH

The previous studies showed the importance of the quasifission process as a mechanism in competition with fusion, hindering the formation of heavy systems. A quantum and microscopic theoretical framework able to describe properly the quasi-fission properties would be of great importance to get a deep insight into the interplay between structure properties and the QF mechanism.

We saw in Section II that quasi-fission was responsible for fusion hindrance in heavy systems. We now use the TDHF approach to investigate the quasi-fission process. Fig. 1 indicates that fusion hindrance appears in systems with larger charge products than in the 48 Ti+ 208 Pb system (see also Ref. [START_REF] Simenel | proceedings of the 11th International Conference on Nucleus-Nucleus Collisions (NN2012)[END_REF]). In the following we study the 40 Ca+ 238 U reaction (Z 1 Z 2 = 1840) around the barrier which, indeed, leads to quasi-fission at the mean-field level [START_REF] Simenel | proceedings of the 11th International Conference on Nucleus-Nucleus Collisions (NN2012)[END_REF]49].

Examples of density evolutions obtained with the tdhf3d code are shown in Fig. 8 for this system at E c.m. = 206 MeV. The left column shows the density for a collision with a tip of the 238 U nucleus at L = 80 , leading to an average exit channel 55 Mn+ 218 Fr. The middle and right columns show a collision with the side of 238 U at L = 20 . In this case the exit channel is more symmetric: 111 Rh+ 167 Ho in average. These reactions correspond to quasi-fission, i.e., an important multi-nucleon transfer from the heavy fragment toward the light one with several zs life-time of the di-nuclear system which is typical for QF [5][6][7][8]. We also see that the mass equilibration (i.e., the formation of two fragments with symmetric masses) is not complete and may depend on the initial conditions.

To get a better insight into the mass transfer mechanism, we plot in Fig. 9 the average mass of the light fragment in the outgoing channel of 40 Ca+ 238 U at E c.m. = 206 MeV as a function of the initial angular momentum. Note that the mass of the fragments are determined as expectation values of one-body operators. Distribution of probabilities around these values could be extracted at the TDHF level thanks to particlenumber projection techniques [50]. However, the calculation of such distributions is going beyond the scope of this paper. We see in Fig. 9 that all collisions with the tip lead to quasifission with partial mass equilibration. In particular, this orientation never leads to fusion, while the other orientation produces long contact time which may lead to fusion at L ≤ 10 . In fact, collisions with the tip favour the production of a heavy fragment in the 208 Pb region up to L = 70 which indicates that these QF reactions are strongly affected by shell effects. On the contrary, collisions with the side does not seem to produce an excess of 208 Pb. This may be due to the fact that more compact configurations are reached with this orientation, favouring more mass symmetric exit channels for the most central collisions. In addition, QF is observed for collisions with the side up to L = 30 only. For larger L, the overlap between the fragments is too small to allow the formation of a di-nuclear system, leading essentially to quasi-elastic reactions. Usual experimental observables to investigate the quasifission process include, in addition to the fragment mass distribution, the kinetic energy of the fragments [13] and the scattering angle [6,51]. Fig. 10 shows the evolution of the scattering angle in 40 Ca+ 238 U at E c.m. = 206 MeV for collisions with the tip of 238 U as a function of angular momentum. For symmetry reasons, central collisions (L = 0 ) induce backward scattering at 180 • , while non-central collisions produce a rotation of the di-nuclear system. The higher the angular momentum, the higher the angular velocity. As a result, increasing L induces quasi-fission with more forward angle down to ∼ 17 • at L = 80 (see Fig. 8 left column for the associated density evolution). This indicates almost half a rotation before emission of the fragments at L = 80 . For L ≥ 100, quasi-elastic scattering occurs with no orbiting, i.e., the nuclei follow a Coulomb trajectory.

The study of correlations between mass transfer and scattering angle, and comparisons with experimental data taken at ANU are ongoing [49].

V. CONCLUSIONS AND OUTLOOKS

The quasi-fission process, which is responsible for the fusion hindrance in heavy system, strongly depends on the entrance channel properties. Measured mass angle distributions in several Ca,Ti+Pb,Hg systems show narrower fragment mass distributions compatible with fusion when the colliding partners have strong shell effects and similar N/Z ratio, such as in the 48 Ca+ 208 Pb system. Microscopic mean-field calculations with the TDHF approach show that magic shells also affect the outgoing channel by favouring the production of 208 Pb-like fragments in the 40 Ca+ 238 U reaction. It is also shown that the amount of mass transfer depends strongly on the orientation of the actinide.

TDHF calculations of heavy-ion collisions provide observables, such as the mass of the fragments and the scattering angle, which can be directly compared to experimental data.

However, a complete determination of mass-angle distributions requires models which include beyond mean-field fluctuations. For instance, extensions of TDHF calculations including fluctuations at the time-dependent RPA level are now available [52]. Alternatively, a stochastic mean-field approach could also be considered [53]. These approaches should be applied to investigate quasi-fission in the near-future.

FIG. 1 :

 1 FIG.1: Contact time for heavy-ion central collisions as a function of center of mass energy normalised to the proximity barrier[42]. The arrow indicates a lower limit.

FIG. 2 :

 2 FIG. 2: Experimental setup for the measurement of MAD.

FIG. 3 :

 3 FIG. 3: (upper panels) Measured MAD. (lower panels) Projected mass ratio spectra. Gaussian fits to the region around M R =0.5 are shown (turquoise lines). Gaussian functions with σ M R = 0.07 (thin red lines) are shown for reference.

FIG. 5 :

 5 FIG.5: Widths of the mass distributions of the fission-like events as a function of the magicity (quantified by the number of magic numbers N m ) in the entrance channel[22].

FIG. 7 :

 7 FIG.7: TDHF calculations of charge equilibration: the difference between final and initial N/Z of the fragments is shown as a function of the contact time. Adapted from[22].

FIG. 8 :

 8 FIG. 8: Snapshots of the TDHF isodensity at half the saturation density (ρ 0 /2 = 0.08 fm -3 ) in the 40 Ca+ 238 U system at E c.m. = 206 MeV. (left) Collision with a tip of 238 U at L = 80 . (middle and right) Collision with the side of 238 U at L = 20 . Snapshots are shown every 1.5 zs.

FIG. 9 :

 9 FIG. 9: Mass of the light final fragment in 40 Ca+ 238 U collisions at E c.m. = 206 MeV as a function of the angular momentum. Collisions with the tip (red circles, solid line) and with the side (blue squares, dashed line) of 238 U are considered.

FIG. 10 :

 10 FIG.10: Scattering angle as a function of angular momentum in collisions of a40 Ca with the tip of a 238 U at E c.m. = 206 MeV.

Acknowledgements

The authors are grateful to N. Rowley for motivating discussions. N. Lobanov and D. C. Weisser are thanked for intensive ion source development. The calculations have been performed on the Centre de Calcul Recherche et Technologie of the Commissariat à l' Énergie Atomique, France, and on the NCI National Facility in Canberra, Australia, which is supported by the Australian Commonwealth Government.

Partial support from ARC Discovery grants DP0879679, DP110102858, and DP110102879 is acknowledged.