Atomic-scale modeling of twinning disconnections in zirconium

TL;DR

This study models atomic-scale disconnections in zirconium twins, revealing how core energies influence stability and correlating with experimental twin modes, while migration plays a lesser role.

Contribution

It introduces a combined atomistic and elasticity approach to quantify disconnection energies and stability in zirconium twin boundaries.

Findings

Most stable disconnections are not always those with smallest Burgers vector.

Crystallographic parameters match experimental twin modes.

Disconnection migration is not critical for twin mode selection.

Abstract

Twin growth in hexagonal close-packed zirconium is investigated at the atomic scale by modeling the various disconnections that can exist on twin boundaries. Thanks to a coupling with elasticity theory, core energies are extracted from atomistic simulations and the formation energy of isolated disconnection dipoles is defined. For twin systems where several disconnections can exist, because of this core contribution, the most stable disconnection is not always the one with the smallest Burgers vector. Crystallographic parameters of the disconnection with the lowest formation energy correlate well with twin modes observed experimentally. On the other hand, disconnection migration, characterized here by computing their migration energy and Peierls stress, does not appear critical for twin mode selection.

Figures

Tables

Table 1. Table 1: Perfect twin boundaries energies in mJ m -2 . The regime of strain corresponding to the experimental activation of the twinning system Yoo and Lee ( 1991 ) is mentioned by C (compression) or T (tension) and the temperature regime by HT (high temperature) and LT (low temperature).

		Twin	EAM	GGA
C	HT	\hkl10-11	150	96
	LT	\hkl11-22	209	355
T	HT	\hkl10-12	264	272
	LT	\hkl11-21	144	229

Table 2. Table 2: Crystallographic definition of the disconnections found stable in zirconium for different twinning systems. The four systems are defined by their twin plane K 1 subscript K 1 \mathrm{K_{1}} , their experimental twinning vector η → 1 subscript → 𝜂 1 \vec{\eta}_{1} , and their vector ζ → → 𝜁 \vec{\zeta} normal to the experimental shear plane. d K 1 subscript 𝑑 subscript 𝐾 1 d_{K_{1}} is the interplanar distance between two K 1 subscript K 1 \mathrm{K_{1}} planes. The disconnection Burgers vectors b → h subscript → 𝑏 ℎ \vec{b}_{h} are decomposed into edge and screw components, b e subscript 𝑏 𝑒 b_{e} and b s subscript 𝑏 𝑠 b_{s} , respectively along η → 1 subscript → 𝜂 1 \vec{\eta}_{1} and ζ → → 𝜁 \vec{\zeta} . The sign of the edge component is calculated for Zr ( γ 𝛾 \gamma =1.598). γ 𝛾 \gamma is the axial c / a 𝑐 𝑎 c/a ratio and every distance in this table is normalized by the lattice constant a 𝑎 a . Energy properties of the disconnections are also defined: core energy E core superscript 𝐸 core E^{\rm core} , elastic energy prefactor k 𝑘 k , core radius r c subscript 𝑟 c r_{\rm c} entering the definition of the formation energy (Eq. 4 ), migration energy E mig superscript 𝐸 mig E^{\rm mig} , and Peierls stress τ P subscript 𝜏 P \tau_{\rm P}

${\overrightarrow{b}}_h$	$b_e$	$b_s$	||${\overrightarrow{b}}_h$||	$E^{\mathrm{core}}$	$k$	$r_c/h$	$E^{\mathrm{mig}}$	${\tau }_P$
				(meV Å-1)	(meV Å-3)		(meV Å-1)	(GPa)
$K_1$=\hkl10-11; $d_{K_1}$=$\frac{\sqrt{3}\gamma }{\sqrt{4{\gamma }^2+3}}$; ${\overrightarrow{\eta }}_1$=$\left[10\overline{1}\overline{2}\right]$; $\overrightarrow{\zeta }$=$\left[1\overline{2}10\right]$
${\overrightarrow{b}}_1$	$\frac{3}{\sqrt{4{\gamma }^2+3}}>0$	0	0.826	-87	742	0.5	10	0.30
${\overrightarrow{b}}_2$		$\pm \frac{1}{2}$	0.529	-4	544	0.28	6	0.17
${\overrightarrow{b}}_{3e}^{-/+}$		0	1.08	236	742	0.5	120 / 208	0.21 / $>0.21$
${\overrightarrow{b}}_{3m}^{-/+}$	$\frac{15-4{\gamma }^2}{2\sqrt{4{\gamma }^2+3}}>0$	$\pm \frac{1}{2}$	0.826	289	648	0.37	25 / 20	$>0.39$ / 0.39
${\overrightarrow{b}}_4$		0	0.334	117	742	0.5	116	$>2$
$K_1$=\hkl11-22; $d_{K_1}$=$\frac{\gamma }{2\sqrt{{\gamma }^2+1}}$; ${\overrightarrow{\eta }}_1$=$\left[11\overline{2}\overline{3}\right]$; $\overrightarrow{\zeta }$=$\left[1\overline{1}00\right]$
${\overrightarrow{b}}_1$	$\frac{1}{\sqrt{{\gamma }^2+1}}>0$	0	0.529	-74	750	0.5	0.2	0.02
${\overrightarrow{b}}_3^{}$		0	0.294	119	750	0.5	22	1.87
${\overrightarrow{b}}_3^{-/+}$		0	0.294	237	750	0.5	40 / 10	$>1.47$ / 1.47
${\overrightarrow{b}}_4^{}$	$\frac{3-{\gamma }^2}{\sqrt{{\gamma }^2+1}}>0$	0	0.238	254	750	0.5	32	1.82
${\overrightarrow{b}}_4^{-/+}$	$\frac{3-{\gamma }^2}{\sqrt{{\gamma }^2+1}}>0$	0	0.238	250	750	0.5	22 / 14	2.38 / $>2.38$
$K_1$=\hkl10-12; $d_{K_1}$=$\frac{\sqrt{3}\gamma }{2\sqrt{3+{\gamma }^2}}$; ${\overrightarrow{\eta }}_1$=$\left[10\overline{1}\overline{1}\right]$; $\overrightarrow{\zeta }$=$\left[1\overline{2}10\right]$
${\overrightarrow{b}}_2$	$\frac{3-{\gamma }^2}{\sqrt{3+{\gamma }^2}}>0$	0	0.189	5	753	0.5	2	0.15
$K_1$=\hkl11-21; $d_{K_1}$=$\frac{\gamma }{\sqrt{4{\gamma }^2+1}}$; ${\overrightarrow{\eta }}_1$=$\left[11\overline{2}\overline{6}\right]$; $\overrightarrow{\zeta }$=$\left[1\overline{1}00\right]$
${\overrightarrow{b}}_{1/2}$	$\frac{1}{2\sqrt{4{\gamma }^2+1}}>0$	0	0.149	-25	758	0.5	< 0.1	< 0.01

Equations

$${\epsilon }_{k3}^0={\epsilon }_{3k}^0=\frac{b_k\ell }{2HL}.$$

$$E^{\mathrm{form}}(\ell ,H,L)=E^{\mathrm{elas}}(\ell ,H,L)+2E^{\mathrm{core}}.$$

$$E^{\mathrm{elas}}(\ell )=\frac{1}{4\pi }k{b}^2\ln \left(\frac{\ell }{r_{\mathrm{c}}}\right).$$

$$E^{\mathrm{form}}(\ell )=\frac{1}{4\pi }k{b}^2\ln \left(\frac{\ell }{r_{\mathrm{c}}}\right)+2E^{\mathrm{core}}.$$

Full text

I Introduction

Twinning is a necessary deformation mode in crystals where only a limited number of dislocation slip systems can be activated. This is the case in hexagonal close-packed (hcp) metals Partridge (1967); Christian and Mahajan (1995) where dislocations with $1/3\,\hkl<1-210>$ Burgers vectors are the main carriers of plastic deformationHirth and Lothe (1982) but cannot accommodate any strain along the \hkl<c> axis. This can be done only through the activation of twinning or the glide of dislocations with $1/3\,\hkl<1-213>$ Burgers vectors. Whereas $1/3\,\hkl<1-213>$ dislocations account for plastic deformation at fairly high temperatures, twinning becomes the dominant mechanism at low temperature or high strain rate.Yoo (1981) While the crystallography of twinning is well asserted,Christian and Mahajan (1995) the mechanisms controlling twin formation and growth are still the object of active research.Beyerlein et al. (2014) Twinning can be divided in three steps: twin nucleation, lateral propagation and thickening.

This article focuses on the last step, where the motion of the twin boundary leads to thickening. This motion occurs thanks to disconnections,Christian and Mahajan (1995); Hirth et al. (2016) i.e. interface dislocations with a step character gliding along the twin boundary to propagate locally the twin in the parent crystal. The presence of disconnections on the twin boundary can result from the interaction of the boundary with dislocations coming from the bulk crystal Serra and Bacon (1995, 1996); Wang et al. (2012); Kadiri et al. (2015); Fan et al. (2015, 2016); Wang and Agnew (2016) or from nucleation under the action of the applied stress and temperature. Ghazisaeidi et al. (2014); Luque et al. (2014) Whether twin thickening is controlled by the formation or the migration of the disconnections is not clearly known. In hcp metals, although several disconnections with different Burgers vectors and / or different step heights are usually possible on a twin boundary,Serra et al. (1988, 1991); Bacon and Vitek (2002) only a single twinning mode, corresponding to a given shear direction and intensity, is reported for each twin system.Partridge (1967); Christian and Mahajan (1995) One can then wonder if the selection of this mode is driven by a competition between the formation of the different possible disconnections or between their migration. Answering such a question requires a description at the atomic scale. Although crystallography and elasticity theory can rationalize many properties of disconnections, the core region of the line defects controls their migration and has a non negligible contribution to the defect formation energy, sometimes counterbalancing the elastic contribution as shown below.

In this perspective, disconnections are modeled here at the atomic scale in Zr, an hcp transition metal of upmost technological interest in the nuclear industry.Lemaignan (2012) The purpose is to extract from these atomistic simulations the key quantities describing disconnection formation and migration, and to confront these results with the twinning modes reported experimentally. The developed method proposes to extract the disconnection core energy from atomistic simulations using small simulation cells and to calculate the elastic energy thanks to an analogy between a disconnection dipole and an Eshelby inclusion, taking full account of the elastic anisotropy and inhomogeneity. Small computation cells are used in order to later perform *ab initio *calculations but for the moment, this approach is validated using a many-body potential of the embedded-atom method (EAM) type. The adequacy of the EAM potential is first assessed by comparison with *ab initio *calculations of perfect twin boundaries. All possible disconnections that have been proposed for the four different twin systems reported in ZrYoo and Lee (1991) are then modeled. The formation energies of disconnection dipoles are then defined considering both core and elastic contributions. The disconnection migration is finally studied by calculating their migration energy and Peierls stress.

II Perfect twin boundaries

The empirical interatomic potential is validated by comparing the structures and energies of perfect twin boundaries in hcp Zr with *ab initio *calculations. The four twin systems that can be activated in Zr,Yoo and Lee (1991) corresponding to \hkl10-11, \hkl11-22, \hkl10-12 and \hkl11-21 twin boundaries, are considered.

II.1 Methods

The atomic interactions are described using the EAM potential developed by Mendelev and Ackland for zirconium and referred to as #3 in reference Mendelev and Ackland, 2007. This potential, developed for bulk hcp Zr with a special emphasis on stacking faults controlling dislocation dissociation, has already been shown to give a reliable description of Zr plasticity.Khater and Bacon (2010); Clouet (2012); Chaari et al. (2014a, b); Lu et al. (2015a, b); Szewc et al. (2016)

Ab initio calculations were performed with the VASP code, Kresse and Furthmüller (1996) using the Perdew-Becke-Erzenhof Perdew et al. (1996) parameterization for the GGA exchange and correlation functional. The interactions between the core and outer electrons are modeled through the projector augmented wave approximation, considering valence (5s2 4d2) and semicore electrons (4s2 4p6) in the outer shell. The plane wave cutoff energy is set at 460 eV. The Brillouin zone is sampled with a regular $\Gamma$-centered mesh corresponding to $20\times 20\times 9$ $k$-points for the primitive hcp crystal cell. The electronic density of states is integrated using the Methfessel-Paxton broadening function with a smearing parameter of 0.1 eV. Atomic positions are relaxed until the force on each atom is less than $5\ 10^{-3}$ eV Å-1.

Twin boundaries are modeled with periodic boundary conditions in all directions, thus introducing two twin planes in the simulation cell. A minimum number of atomic planes has to be kept between the two twin planes to limit their interactions. Relaxation perpendicular to the twin boundary is considered by allowing an increase or a decrease of the periodicity vector in this direction, with the corresponding displacement $\delta_{z}$ localized in the planes adjacent to the twin boundaries.Kumar et al. (2015) Such a relaxation improves the convergence of the twin energy with the number of atomic planes in the simulation cell. For ab initio calculations, cells containing nine atomic planes between twin planes are used leading to a convergence better than 5 mJ m-2 on twin energies.

II.2 Atomic structure

Twinning operations can either correspond to a reflection in the twinning plane $K_{1}$ (Type I), or a rotation of $\pi$ along the twinning direction $\eta_{1}$ (Type II).Christian and Mahajan (1995) \hkl10-11 and \hkl10-12 twin systems are simultaneously of type I and II because in these cases, both operations are equivalent. On the other hand, the atomic structures obtained by type I or type II twinning are not the same for the \hkl11-21 and \hkl11-22 twin systems. Structures corresponding to both types have been relaxed with ab initio and the EAM potential. It was found that type II is unstable for \hkl11-22 and relaxes to the type I twin structure in both ab initio and EAM. For \hkl11-21 twins however, both types are stable but the type I structure has a much lower energy, 144 mJ m-2 (229 mJ m-2 in ab initio) compared to 610 mJ m*-2* for the type II structure (749 mJ m*-2* in ab initio). This is in agreement with previous *ab initio *calculations, not only in Zr Kumar et al. (2015); Morris et al. (1994, 1995) but also in other hcp elements,de Jong et al. (2015); Ni et al. (2015) and with the structures observed experimentally with high resolution electron microscopy in various hcp metals such as Ti,Kasukabe et al. (1993); Pond et al. (1995); Braisaz et al. (1996) Zn,Lay and Nouet (1994); Braisaz et al. (1997) Co,Zhang et al. (2012) or Mg.Sun et al. (2014) The structures predicted by the EAM potential (figure 1) agree with *ab initio *results with a slight exception for the \hkl11-21 twin boundary (figure 1d). The stable structure obtained with both the EAM potential and by ab initio is of type I. With the EAM potential however, there is an additional translation parallel to the interface. The amplitude of this displacement, 0.07 $a$ is small ($a=3.234$ Å is the lattice parameter). The interface without the translation (structure drawn in grey in fig. 1d) has an energy only 3 mJ m-2 higher than the fully relaxed structure but is unstable. This artifact of the EAM potential is therefore considered harmless for the study of the \hkl11-21 twinning system. A similar stable structure for this perfect twin boundary has been observed by Bacon and Serra Bacon and Serra (1991) in Ti with a many-body potential of the Finnis-Sinclair form.

II.3 Energies of twin boundaries

The energies of perfect twin boundaries are gathered in table 1. Different values are obtained with the EAM potential and the *ab initio *calculations. However, if one looks separately to the twin systems activated under compression or tension, Yoo and Lee (1991) the same relative stability is obtained with both methods. For compression twins, the \hkl10-11 boundary is more stable than \hkl11-22, and for tension, \hkl11-21 is more stable than \hkl10-12. A qualitative agreement is obtained between the EAM potential and *ab initio *calculations, similar to the agreement previously obtained on stacking fault energies.Mendelev and Ackland (2007); Chaari et al. (2014a)

One can conclude that this empirical potential leads to a reasonable description of the various twinning systems that can be activated in Zr, thus justifying its use to study twin growth at an atomic scale.

III Disconnection formation

Twin boundaries move thanks to the migration of disconnections.Christian and Mahajan (1995); Hirth et al. (2016) Disconnections are formed either by absorption of dislocations coming from the bulk Serra and Bacon (1995); Fan et al. (2015, 2016); Wang et al. (2012) or by nucleation on the twin boundary.Luque et al. (2014) Once formed, disconnections migrate along the twin boundary and propagate the twinned crystal. The following section describes the structure of the different disconnections that can appear on the four twinning systems activated in Zr and define their formation energy, before addressing their migration in the next section.

III.1 Crystallography

Disconnections are defined as steps along the twin boundary with a dislocation character.Hirth and Pond (1996) The step height is called $h$ and for purpose of simplicity is hereafter divided by the interplanar spacing $d_{K_{1}}$, so that $h$ is an integer. The dislocation content is characterized by the Burgers vector $\vec{b}_{h}$ corresponding to the shear needed to make the twin grow by $h$ planes at the expense of the parent crystal. Several disconnections corresponding to different $h$ and $\vec{b}_{h}$ can exist for a given twinning system. Keeping only the possibilities leading to the smallest $h$ and $\vec{b}_{h}$, as in previous atomistic studies, Serra et al. (1988, 1991); Serra and Bacon (1995, 1996); Wang et al. (2011); Li et al. (2012); Khater et al. (2013) table 2 gathers the crystallographic definitions of the different disconnections considered in this study (see also the dichromatic patterns in figure 2). All these disconnections have been found stable. In this table, the Burgers vector is decomposed into an edge and a screw component, with the line direction $\vec{\zeta}$ of the disconnection taken normal to the experimental twinning direction $\vec{\eta}_{1}$.Christian and Mahajan (1995) Due to this choice, the twin modes observed experimentally correspond to pure edge disconnections. The orientation has also been chosen so that compression (respectively tension) along the \hkl<c> axis activates twinning through the glide of disconnections with negative (respectively positive) edge component.

From a purely elastic point of view, the disconnections with the smallest Burgers vectors are expected to have the lowest formation energies. However, the steps also induce atomic rearrangements that represent a core energy cost. This contribution can be intuitively expected to be minimal for the smallest step heights. Since the smallest Burgers vectors are usually not obtained for the smallest step heights (Tab. II), the competition between the core and elastic energies prevents to predict which disconnection is the most stable for a given twinning system.

In the particular case of the \hkl11-21 twin system, the $\vec{b}_{1}$ disconnection is not stable and spontaneously dissociates into two $\vec{b}_{1/2}$ disconnections.Serra et al. (1988, 1991); Khater et al. (2013) Although the height of this disconnection is not an integer, the $\vec{b}_{h}$ notation will be kept.

A shuffling is imposed to the atoms between the former and the new twin planes following the method described by Serra et al.Serra et al. (1988, 1991) Several initial shuffling may be possible, leading in some cases to different disconnection cores after relaxation, that is to say disconnections with the same $\vec{b}_{h}$ and $h$ but different atomic structures and thus different formation and migration energies. This happens for $\vec{b}_{3}$ and $\vec{b}_{4}$ disconnections on \hkl11-22 twin planes, where two different configurations of the disconnection dipole can be stabilized, depending on the initial shuffling.

Most disconnection dipoles are symmetrical: both disconnections are then equivalent. But some asymmetric dipoles have also been found. The two different core regions are then characterized by a different atomic density. The disconnection with the higher (respectively the lower) density will be labeled $\vec{b}_{h}^{\mathrm{+}}$ (respectively $\vec{b}_{h}^{\mathrm{-}}$).

III.2 Atomistic simulations

The atomistic simulations of disconnections are performed with the EAM potential described in the previous section. The twin plane lies in the $(xy)$ plane with the $x$-direction along the $\vec{\eta}_{1}$ twinning direction. The unitary cell is repeated from thirty up to two hundred times in this direction, leading to a cell length $L$ along $x$ (Fig. 3). The cell dimension $H$ in the $z$-direction normal to the twin plane is chosen large enough for the twin boundary not to interact with its periodic images. In practice, cells with up to sixty planes between the two twin boundaries have been used.

A disconnection dipole of length $\ell$, with a line vector $\vec{\zeta}$ along the $y$-direction, is introduced in the simulation box by displacing all atoms according to the elastic field created by a dislocation dipole of Burgers vector $\vec{b}_{h}$ located at a distance $h$ above the twin boundary, taking full account of periodic boundary conditions.Cai et al. (2003)

An homogeneous strain $\varepsilon^{0}$, corresponding to the plastic strain introduced in the simulation cell by the creation of the disconnection dipole, is also applied to cancel the stress in the cell.Rodney et al. (2017) Its non-zero components are given by:

[TABLE]

Atomic positions are then relaxed thanks to a conjugate gradient algorithm. As this relaxation can change the length $\ell$ of the dipole, a residual stress may be observed after relaxation. The applied homogeneous strain is adjusted to cancel this stress. Further relaxation does not change the atomic structure of the disconnections and leads to zero stress. In the following, the actual length $\ell$ of the disconnection dipole will be defined through equation 1, using the strain $\varepsilon^{0}_{13}$ that has been applied to the simulation box to cancel the stress. In cases where the applied strain $\varepsilon^{0}_{23}$ is not zero, i.e. when the disconnection has a screw component, it has been found that the lengths $\ell$ deduced either from $\varepsilon^{0}_{13}$ or $\varepsilon^{0}_{23}$ do not differ by more than 1 Å.

Figure 4 represents the relaxed atomic structure of a $\vec{b}_{2}$ disconnection dipole on the \hkl10-12 twin boundary. Two steps can be clearly seen on the interface separating the parent and twinned crystals. Away from these steps, one recovers the structure of the perfect twin boundary. Nye tensor can be extracted from the relaxed configurations Hartley and Mishin (2005) to visualize the Burgers vector density created by both disconnections. The density is located at the steps and vanishes over a distance approximately equal to five times the lattice parameter $a$, confirming the localized character of the shear incompatibility. The dipole length $\ell$ calculated with equation 1 is reported on this figure. This definition of $\ell$ agrees with the dipole length that can be inferred from the relaxed atomic structure. The same behavior is observed for the different disconnection dipoles studied in the present work.

III.3 Formation energy

The formation energy $E^{\mathrm{form}}$ of a disconnection dipole is defined as the energy difference between a simulation cell containing a disconnection dipole along one of its twin boundary and the same cell with two perfect twin boundaries, both cells being under zero stress. The formation energy depends not only on the length $\ell$ of the disconnection dipole, but also on the dimensions of the simulation cell. This is illustrated in figure 5 where the formation energy of the $\vec{b}_{2}$ disconnection dipole on the \hkl10-12 twin boundary is displayed as a function of $\ell$ for various cell heights $H$ and a fixed cell width $L$. The energy variations are caused by the interactions between the two disconnections composing the dipole and their periodic images. As these interactions are elastic, one can decompose the formation energy in two parts, an elastic and a core contribution:

[TABLE]

The elastic energy $E^{\mathrm{elas}}$, which contains all the interactions with the surrounding microstructure, varies with the cell dimensions and with the dipole length. On the other hand, the core contribution $E^{\mathrm{core}}$ accounts for the cost of atomic disorder at the very core of the disconnections, which cannot be described by linear elasticity. This contribution should be an intrinsic property of the disconnections and thus should not depend on the surrounding microstructure. For symmetric dipoles, both disconnections have the same core energy, explaining the factor 2 in equation 2. For asymmetric dipoles, $E^{\rm core}$ is then the average core energy between the $\vec{b}_{h}^{\mathrm{+}}$ and $\vec{b}_{h}^{\mathrm{-}}$ disconnections.

To check the validity of the energy decomposition in equation 2, the elastic energy of the disconnection dipole is calculated separately. In this elastic calculation, the dipole is modeled as an Eshelby inclusion Mura (1987) of length $\ell$ and height $h$ with an eigenstrain $\varepsilon^{*}_{i3}=b_{i}/2h$ (Fig. 6). The elastic energy is computed within linear elasticity theory thanks to the Fast Fourier Transform (FFT) approach of Moulinec and Suquet,Moulinec and Suquet (1998) taking full account of elastic anisotropy, elastic inhomogeneity, and periodic boundary conditions.

As one can see in figure 5, the calculated elastic energy shows the same variations with dipole length and cell dimensions as the formation energy extracted from the atomistic simulations, confirming that these variations are elastic in nature. Withdrawing the elastic contribution from the formation energy, one obtains a constant energy contribution corresponding to the core energy (Fig. 5). The variation of core energies obtained with such a procedure does not exceed 2 meV Å*-1*, confirming that the decomposition of the formation energy into elastic and core contributions (Eq. 2) is meaningful.

The core energies obtained for all disconnection dipoles are given in table 2. As the elastic energy may be larger than the formation energy, negative core energies are obtained for some of the disconnections. We will show in the next paragraph that this is simply a consequence of the energy decomposition assumed when modeling the disconnection dipole as an Eshelby inclusion. One has to remember that this decomposition is arbitrary and that despite those negative core energies, the formation energy of all dipoles is positive.

III.4 Isolated disconnection dipoles

Now that the decomposition of the formation energy into an elastic and a core contribution (Eq. 2) is validated, the elastic energy of an isolated dipole is calculated in order to withdraw the interactions with the periodic images.

The elastic energy of an isolated dipole is calculated using the same FFT approach as in the previous section. In practice, the cell used for this calculation remains periodic, but the cell dimensions $H$ and $L$ are chosen much larger than the dipole dimensions $h$ and $\ell$, so that the result becomes independent of these cell dimensions. Figure 7 shows that the elastic energy varies as the logarithm of the dipole length following the same expression as for a dislocation dipole:

[TABLE]

The FFT calculations differ from this analytical expression only at small dipole lengths. However, as soon as the dipole length becomes larger than its height ($\ell\gg h$), the elastic energy can be expressed analytically through Eq. 3. $k$ is a constant that depends only on the elastic constants and the disconnection orientation Dupeux and Bonnet (1980) while $r_{\rm c}$ is an effective core radius. These two parameters are obtained by fitting equation 3 to the results of the FFT calculations at the large dipole lengths and results are displayed in table 2. The main variations of $k$ are caused by the disconnection character, with $k$ being larger for pure edge than for mixed disconnections. Elastic anisotropy also leads to a slight variation of this parameter. On the other hand, the ratio $r_{\rm c}/h$ only depends on the disconnection character, with $r_{\rm c}\approx h/2$ for pure edge and smaller values for mixed orientations. The effective core radius entering the definition of the elastic energy (Eq. 3) is therefore much smaller than the actual extension of the disconnection core, as can be inferred for instance from the spreading of the Burgers vector density obtained from the Nye tensor (Fig. 4). As a consequence, the region where linear elasticity breaks down is more extended than the cylinder of radius $r_{\rm c}$ around the disconnection line.

Combining equations 2 and 3 gives the formation energy of an isolated disconnection dipole.

[TABLE]

The core radius inferred from FFT calculations being smaller than the actual extent of the disconnections is the reason why negative core energies are obtained for some of the disconnections (Tab. 2).

Figure 8 shows the formation energy (Eq. 4) as a function of the disconnection dipole length for all isolated disconnection dipoles that exist in zirconium. Since the atomic simulations are relying on a EAM potential, one can also use these simulations with large enough cells to obtain a disconnection dipole isolated from its periodic images. This offers a way to validate the analytical model. The points on figures 8a and 8c are the results of those atomistic calculations and show a quantitative agreement with the analytical expression. As the size of the simulation cells used to obtain the core energy are compatible with *ab initio *calculations, one can conceive using the same approach based on the partition of the disconnection energy into a core and an elastic contributions to obtain an *ab initio *prediction of the formation energy of the disconnection dipole. This is however left for future work.

Because of the two contributions entering the formation energy, the most stable disconnection on a twin boundary is not always the one with the smallest Burgers vector. This is the case for the \hkl11-22 twin boundary where the symmetric configuration of the $\vec{b}_{3}$ disconnection dipole has a low core energy that compensates its elastic energy. As a consequence, in the considered range of dipole lengths ($10\leq\ell\leq 300$ Å), this disconnection dipole has a lower formation energy than the $\vec{b}_{4}$ dipole (Fig. 8), despite its higher Burgers vector (Tab. 3). The $\vec{b}_{4}$ disconnection dipole will become more stable only for lengths $\ell$ larger than 1.5 106 Å, a length too large to be meaningful for disconnection nucleation. Both $\vec{b}_{3}$ and $\vec{b}_{4}$ disconnections are pure edge with a Burgers vector collinear to the twinning direction \hkl<11-2-3>, which has been experimentally reported for this \hkl11-22 twin boundary.Rapperport (1959); Rapperport and Hartley (1960); Partridge (1967); Paton and Backofen (1970); Yoo (1981); Christian and Mahajan (1995) But because their Burgers vectors have edge components of opposites signs (Tab. 2), the $\vec{b}_{3}$ disconnection dipole propagates the twin under compression of the \hkl<c> axis, whereas the $\vec{b}_{4}$ is activated under tension. The lower formation energy obtained for the $\vec{b}_{3}$ disconnection dipole is therefore in agreement with the \hkl11-22 twinning system being active in compression.Akhtar (1973); Yoo (1981) The shear magnitude $s=-b_{3}/3d_{\hkl{11-22}}=0.231$ induced by the glide of this disconnection also corresponds to the magnitude $s=0.225$ determined experimentally.Rapperport (1959); Rapperport and Hartley (1960)

For the \hkl10-11 twin boundary, the $\vec{b}_{4}$ disconnection dipole is the most stable for dipole lengths larger than $20$ Å (Fig. 8). This disconnection has a Burgers vector in agreement with the twinning direction experimentally reported, Paton and Backofen (1970); Yoo (1981); Christian and Mahajan (1995) with a sign corresponding to the \hkl10-11 twin system active under compression, Akhtar (1973); Yoo (1981) and a shear amplitude $s=-b_{4}/4d_{\hkl{10-11}}=0.109$ corresponding to the amplitude $s=0.104$ experimentally assessed in Zr.Jensen and Backofen (1972) But for small dipole lengths, the $\vec{b}_{2}$ dipole becomes more stable because of its low core energy. This disconnection dipole, which will also be active under compression, is not pure edge and thus corresponds to a different twinning direction. As noted by Serra et al. Serra et al. (1991) this second \hkl10-11 twinning mode has been proposed Christian and Mahajan (1995) but not observed experimentally.

For the \hkl11-21 twin, our model gives a negative formation energy below 25 Å. This is unphysical and illustrates the limits of our approach for this very low energy disconnection: because of the large spreading of this disconnection, one cannot assume that the two disconnections composing the dipole are interacting only through their elastic field for such small separation distances.

Finally, it is worth noting that the formation energies of the most stable disconnections vary strongly between the four different twin systems (Fig. 8). These energies are higher for the two systems active under compression, \hkl10-11 and \hkl11-22, and lower for the tension systems, \hkl10-12 and \hkl11-21. This is mainly due to the low core energies of disconnections activated in tension.

IV Disconnection migration

After modeling disconnection dipole formation, we now focus on their migration. Both the migration energy and Peierls stress of a disconnection are calculated, i.e. respectively the energy barrier to be overcome by the disconnection to glide without any applied stress and the resolved shear stress which cancels this barrier, thus allowing the disconnection to glide without the help of thermal activation.

IV.1 Migration energy

Migration energies are calculated with the nudged elastic band (NEB) method.Henkelman and Jonsson (2000) The initial and final states of the NEB chain correspond to the same disconnection dipole with a length $\ell$ differing by one periodicity vector. The dipoles are chosen large enough so that this length variation leads to a negligible variation of the elastic energy along the path. The energy variation given by the NEB calculation then corresponds directly to the disconnection migration energy. When dipoles are asymmetrical, we performed NEB calculations for both disconnections separately in order to obtain the two different migration barriers. Calculations are performed under zero applied stress. The NEB calculations were performed with the same periodicity vectors for each replica. It implies, through equation 1, that the applied stress along the path is not rigorously constant. The variations are however small and we chose in practice the periodicity vectors canceling the stress for a disconnection dipole halfway between the initial and final states.

Finally, one should point out that the length of the disconnection dipoles along their line direction are minimal and thus does not allow for the formation of kinks.Combe et al. (2016) The obtained energy barriers therefore correspond to the 1D migration of the disconnections. They constitute an upper limit of the 3D energy barriers when disconnection migration proceeds via the nucleation of double kinks and are a necessary input for line tension or elastic models predicting kink pair formation.Proville et al. (2013); Kraych et al. (2016)

Results displayed in figure 9 show the barriers for the different disconnections on the four twin boundaries. It can be observed that the energetic landscapes are different for all disconnection cores. On the other hand, the migration energies are in general much lower than the formation energies (Fig. 8).

For the \hkl10-11 twin system, only the $\vec{b}_{4}$ disconnection is compatible with the experimental twinning elements (twinning direction, intensity and activity under \hkl<c> compression). This disconnection has nevertheless a migration energy one order magnitude higher than the $\vec{b}_{1}$, $\vec{b}_{2}$ and $\vec{b}_{3m}$ disconnections, which can also exist on this twinning plane. The migration energy therefor does not appear as a key factor in the selection of the disconnection responsible for twin growth. The same conclusion is reached for the \hkl11-22 twin system where the $\vec{b}_{3}$ disconnections, which are the only ones that can account for twin growth under compression, do not have the lowest migration energy.

IV.2 Peierls stress

In addition to migration energies, the Peierls stress $\tau_{\rm P}$ has been determined by applying gradually a strain producing a resolved shear stress on the disconnections in their glide direction (see Tab. 2). The Peierls stress corresponds to the critical resolved shear stress for which a disconnection of the dipole moves by at least one Peierls valley without the help of thermal fluctuations. We checked by performing NEB calculations under applied stresses that the Peierls stress thus obtained corresponds to the applied stress which cancels the disconnection migration barrier. An example is shown in figure 10 for the $\vec{b}_{3}^{\mathrm{S}}$ disconnection of the \hkl11-22 twin system. The direct calculation leads to a Peierls stress of 1.87 GPa (Tab. 2), in agreement with the NEB calculation showing almost no energy barrier for an applied stress of 1.8 GPa (Fig. 10). We note that for asymmetric dipoles, only the lowest Peierls stress is obtained.

As for the migration energy, the obtained Peierls stresses (Tab. 2) exhibit a large variability amongst disconnections. For the twinning systems where several different disconnections exist, the ones corresponding to the experimental twinning elements, $\vec{b}_{4}$ on \hkl10-11 and $\vec{b}_{3}$ on \hkl11-22, are not the ones with the lowest Peierls stress.

The $\vec{b}_{4}$ disconnection, which is responsible for \hkl10-11 twin growth, shows a particular behavior. It was not possible to make this disconnection glide, even for an applied stress as large as 2 GPa. This is confirmed by NEB calculations under applied stresses (Fig. 11), showing that an energy barrier still remains with an applied stress of 2 GPa. The motion of the $\vec{b}_{4}$ disconnection along the \hkl10-11 plane takes place in four successive steps, with three intermediate metastable configurations. Despite a different sensitivity to the applied stress, none of the energy barriers associated with these successive steps disappears under applied stress. In particular, the first barrier appears uncoupled with the applied stress. Migration of the $\vec{b}_{4}$ disconnection will therefore require thermal activation which may be the reason why \hkl10-11 twinning is active only at high temperature, Yoo and Lee (1991) typically above $400\,^{\circ}\mathrm{C}$ in Ti.

V Conclusions

The two stages of disconnection nucleation relevant for twin thickening, i.e. the formation of disconnection dipoles and migration of disconnections along twin boundaries, have been modeled in hcp zirconium. Using an EAM potential validated on *ab initio *calculations, four different twin systems were studied at an atomic scale.

A precise definition of the disconnection dipole energy is allowed thanks to an approach coupling elasticity and atomistic developed here. It was shown that the formation energy of a disconnection dipole is composed of an elastic contribution, which depends on the surrounding microstructure, and a core energy, which is an intrinsic property of the disconnections. Computing separately the elastic contribution, thanks to an analogy with an Eshelby inclusion, core energies could be extracted from atomistic simulations to define the formation energy of an isolated dipole. This method was shown to reproduce well the behavior of large systems despite being parameterized on simulations of small sizes. Of course, once the core energy is known, it becomes possible to compute formation energies in more complex environments, like stress concentrations, or for more complex geometries, like disconnections loops.Luque et al. (2014) In the latter case, the variation of the core energy with the disconnection character will need to be known. The present study was restricted to disconnections normal to the twinning direction, that usually have a pure edge character. But the same approach can be applied to various line orientations of the disconnection dipole, to obtain the variation of the core energy with the disconnection character.

Considering the different disconnections that can exist on twin boundaries, it has been shown that the most stable is not always the one with the smallest Burgers vector, i.e. the one with the lowest elastic energy. Because of the core contribution, it sometimes happens that a disconnection with a larger Burgers vector is more stable over a wide range of dipole lengths. It appears that the energy scale controlling disconnection formation and migration in Zr is much higher for the \hkl10-11 and \hkl11-22 twin systems active under compression than for the \hkl10-12 and \hkl11-21 tensile twin systems. Considering both the core and elastic contributions, atomistic simulations predict that the twinning modes experimentally observed in hcp Zr correspond to the disconnections with the lowest formation energy. On the other hand, these disconnections are not necessarily the ones migrating the most easily. They generally neither have the lowest migration energy nor the lowest Peierls stress. As a consequence the nucleation of disconnection dipoles appears more critical than their migration for twin mode selection in Zr.

Acknowledgements.

This work was performed using HPC resources from GENCI-CINES and GENCI-CCRT (Grants 2016-096847). DR acknowledges support from LABEX iMUST (ANR-10-LABX-0064) of Université de Lyon (program “Investissements d’Avenir”, ANR-11-IDEX-0007)