New Method for Calculation of Radiation Defect Dipole Tensor and Its Application to Di-Interstitials in Copper

: The effect of external and internal elastic strain ﬁelds on the anisotropic diffusion of radiation defects (RDs) can be taken into account if one knows the dipole tensor of saddle-point conﬁgurations of the diffusing RDs. It is usually calculated by molecular statics, since the insufﬁcient accuracy of the available experimental techniques makes determining it experimentally difﬁcult. However, for an RD with multiple crystallographically non-equivalent metastable and saddle-point conﬁgurations (as in the case of di-interstitials), the problem becomes practically unsolvable due to its complexity. In this paper, we used a different approach to solving this problem. The molecular dynamics (MD) method was used to calculate the strain dependences of the RD diffusion tensor for various types of strain states. These dependences were used to determine the dipole tensor of the effective RD saddle-point conﬁguration, which takes into account the contributions of all real saddle-point conﬁgurations. The proposed approach was used for studying the diffusion characteristics of RDs, such as di-interstitials in FCC copper (used in plasma-facing components of fusion reactors under development). The effect of the external elastic ﬁeld on the MD-calculated normalized diffusion tensor (ratio of the diffusion tensor to a third of its trace) of di-interstitials was fully consistent with analytical predictions based on the kinetic theory, the parameters of which were the components of the dipole tensors, including the range of non-linear dependence of the diffusion tensor on strains. The results obtained allowed for one to simulate the anisotropic diffusion of di-interstitials in external and internal elastic ﬁelds, and to take into account the contribution of di-interstitials to the radiation deformation of crystals. This contribution can be signiﬁcant, as MD data on the primary radiation damage in copper shows that ~20% of self-interstitial atoms produced by cascades of atomic collisions are combined into di-interstitials.


Introduction
The efficiency of the development of structural and functional materials depends on the level of available knowledge of the laws and mechanisms for the formation of radiation microstructures, defects, and properties of the materials. High-dose reactor tests of materials for determination of the degree of radiation change (degradation) of their structure, and properties are very expensive and time-consuming. Moreover, it is difficult to obtain a neutron spectrum in a fission reactor relevant to a fusion reactor environment. To support the development of advanced materials for nuclear fusion technology, theoretical, modelling, and simulation research is necessary to determine parameters of materials and to build physically based models of material property changes under irradiation in conditions specific to fusion reactors.
To build physical models of crystal deformation under damaging radiation (e.g., radiation creep, radiation void swelling), it is necessary to determine the effect of external and internal elastic strain fields on the anisotropic diffusion of radiation defects [1][2][3][4]. This Symmetry 2021, 13, 1154 2 of 11 effect can be taken into account if one knows the dipole tensor of saddle-point configurations of the diffusing RDs [5][6][7]. The dipole tensor P ij (i, j = 1, 2, 3) is usually calculated by molecular statics [7][8][9][10][11][12][13], since the insufficient accuracy of the available experimental techniques makes determining it experimentally difficult. Different techniques within molecular statics have been developed for calculation of the dipole tensor. The first one is to calculate this quantity by summing the first moments of the forces acting on the atoms, fixed at their perfect lattice positions, in the region surrounding the region of free atoms in which the defect is located [7][8][9][10]. The second technique uses a direct relation between the dipole tensor and the strain derivative of the defect formation energy, calculated at constant strain [11]. The third technique uses the proportionality of the dipole tensor to the volume-averaged stress caused by a defect in the simulation box with periodic boundary conditions [12,13]. All of these techniques give the same result. However, for an RD with multiple crystallographically non-equivalent metastable and saddle-point configurations, the problem of anisotropic diffusion modelling becomes practically unsolvable, due to the need to use physically unreasonable simplifications or assumptions.
In this paper, we used a different method for solving this problem. It was demonstrated in [14] that the molecular dynamics (MD) method can be successfully used to study anisotropic diffusion of point defects in metals under external stress fields. We used the MD method to calculate the strain dependences of the RD diffusion tensor D ij for various types of strain states given by the strain tensor ε kl (k, l = 1, 2, 3). These dependences were used to determine the RD elastodiffusion tensor d ijkl , the components of which are the coefficients for linear terms of the expansions in the degrees of strain of the calculated strain dependences. The crystal symmetry and the defect diffusion mechanism determine the symmetry of the tensor, which corresponds to the dipole tensor of the effective saddlepoint configuration of the defect, which takes into account the contributions of all real saddle-point configurations. Figure 1 schematically illustrates the idea of the proposed method for calculation of P ij . and internal elastic strain fields on the anisotropic diffusion of radiation defects [1][2][3][4]. This effect can be taken into account if one knows the dipole tensor of saddle-point configurations of the diffusing RDs [5][6][7]. The dipole tensor Pij (i, j = 1, 2, 3) is usually calculated by molecular statics [7][8][9][10][11][12][13], since the insufficient accuracy of the available experimental techniques makes determining it experimentally difficult. Different techniques within molecular statics have been developed for calculation of the dipole tensor. The first one is to calculate this quantity by summing the first moments of the forces acting on the atoms, fixed at their perfect lattice positions, in the region surrounding the region of free atoms in which the defect is located [7][8][9][10]. The second technique uses a direct relation between the dipole tensor and the strain derivative of the defect formation energy, calculated at constant strain [11]. The third technique uses the proportionality of the dipole tensor to the volume-averaged stress caused by a defect in the simulation box with periodic boundary conditions [12,13]. All of these techniques give the same result. However, for an RD with multiple crystallographically non-equivalent metastable and saddle-point configurations, the problem of anisotropic diffusion modelling becomes practically unsolvable, due to the need to use physically unreasonable simplifications or assumptions.
In this paper, we used a different method for solving this problem. It was demonstrated in [14] that the molecular dynamics (MD) method can be successfully used to study anisotropic diffusion of point defects in metals under external stress fields. We used the MD method to calculate the strain dependences of the RD diffusion tensor Dij for various types of strain states given by the strain tensor εkl (k, l = 1, 2, 3). These dependences were used to determine the RD elastodiffusion tensor dijkl, the components of which are the coefficients for linear terms of the expansions in the degrees of strain of the calculated strain dependences. The crystal symmetry and the defect diffusion mechanism determine the symmetry of the tensor, which corresponds to the dipole tensor of the effective saddlepoint configuration of the defect, which takes into account the contributions of all real saddle-point configurations. Figure 1 schematically illustrates the idea of the proposed method for calculation of Pij. For crystals of a cubic crystal system, the diffusion tensor for small strains can be written as [1]: where D(0) = 1/3 Tr Dij(0) is a third of the trace of the tensor Dij, in the absence of external strains, and δij is the Kronecker symbol (here and hereafter, the Einstein summation convention for a set of indexed terms is used). The tensor dijkl has the same symmetry as the tensor of elastic constants of the crystal cijkl, which means that it has three independent components, d11, d12, d44 (in the Voigt notation), for cubic crystal system crystals. As the three independent values defining the tensor dijkl, it is more convenient to choose its three eigenvalues: d (1) = d11 + 2d12, d (2) = d11 − d12, d (4) = 2d44. For crystals of a cubic crystal system, the diffusion tensor for small strains can be written as [1]: where D(0) = 1/3 Tr D ij (0) is a third of the trace of the tensor D ij , in the absence of external strains, and δ ij is the Kronecker symbol (here and hereafter, the Einstein summation convention for a set of indexed terms is used). The tensor d ijkl has the same symmetry as the tensor of elastic constants of the crystal c ijkl , which means that it has three independent components, d 11 , d 12 , d 44 (in the Voigt notation), for cubic crystal system crystals. As the three independent values defining the tensor d ijkl , it is more convenient to choose its three eigenvalues: d (1) = d 11 + 2d 12 , d (2) = d 11 − d 12 , d (4) = 2d 44 . The strain tensor can be represented in the form [1] (here and hereafter, a crystallographic coordinate system (CCS) with axes along 100 is used, unless otherwise specified): where Rewriting Equation (1), taking into account (2, 3), gives: Thus, for a full determination of d ijkl in the cubic crystal system crystals, it is sufficient to consider only three specific strain states, for example, V 1 , V 2 , V 4 . Knowing d ijkl , we could determine the dipole tensor of the effective RD saddle-point configuration [1], allowing us to simulate the diffusion of defects in elastic fields of any arbitrary type, including inhomogeneous ones (e.g., dislocation fields).
In this paper, the proposed method was used to determine the dipole tensor of the effective saddle-point configuration of di-interstitials in an FCC copper crystal. Copper and copper-based alloys (e.g., CuCrZr alloy) are planned to be used as functional materials of fusion reactors [15][16][17], since these materials possess excellent thermal and electrical conductivity. Di-interstitials are one of the most frequently formed types of self-interstitial atomic clusters under damaging irradiation (cluster analysis of MD data, spatial distributions of self-interstitial atoms (SIAs) produced by atomic collision cascades with damaging energies from 1 to 50 keV in Cu, suggested that~20% of SIAs are combined into di-interstitials [18][19][20]). They have multiple different metastable and saddle-point configurations, which makes them convenient to demonstrate the physical advantages of the method used on this type of RDs. The components P ij for di-interstitials are approximately twice as large as for single SIAs [10]. This leads to a stronger elastic interaction of the di-interstitials with external elastic fields compared to single SIAs, since in the linear approximation, the interaction energy is equal to [5,11].
Therefore, the contribution of di-interstitials to crystal radiation deformation under cascade-forming irradiation becomes significant, despite the fact that less of them are generated as compared to single SIAs.

Molecular Dynamic Model
The MD techniques used in this work closely follow the ones described in [21], which contains their detailed explanations. Di-interstitial diffusion was simulated using interatomic interaction potential [22], since it describes the set of experimental data on the bulk properties of copper crystals and their self-point defects well [22,23]. The model crystallite was a microcanonical ensemble with periodic boundary conditions. The crystal lattice constant a, before the onset of strains, was chosen in such a way that at a given temperature T, the pressure P in the crystallite was equal to zero, with an accuracy of 0.1 eV/nm 3 . The time step value was chosen in such a way that the mean atom displacement per iteration was~0.005a. Verlet integration was used to integrate the equations of motion [24]. The locations of the two SIAs making up the di-interstitial were determined by analyzing the number of atoms in the Wigner-Seitz cells (WSC). The position of the di-interstitial was taken to be that of one of the two SIAs.
The model crystallites were rectangular parallelepipeds. Their characteristics (edge directions, lengths, number of atoms in the crystallite) are summarized in Table 1. The model crystallites were strained by giving the model crystallite atoms displacements u i = ε ij x j , where x j are the atomic coordinates before straining.

Calculated Diffusion Characteristics
To calculate the RD diffusion tensor D ij , each simulated diffusion trajectory of the defect was divided into several isochronous intervals of duration τ. The value of τ for each temperature and strain was set from a compromise between the need for a large number of isochronous intervals to increase the calculation accuracy and the need for a sufficient number of jumps to be contained within the isochronous intervals to reproduce the main correlations and physical features of the RD diffusion [21]. For each isochronous interval, the di-interstitial displacement vector R nm was calculated; here, n is the trajectory number, and m is the isochronous interval number. Next, the quantity D = R nm ⊗R nm /2τ was calculated; here, D is the matrix defining the tensor D ij in the CCS, and the operation . . . denotes averaging over all n and m. The diffusivity D was defined as Tr D ij /3. The errors of D and the tensor components D ij were estimated as standard errors for the specified set of n and m. A total of 100 diffusion trajectories were simulated for each given set of strain state and temperature. The total physical simulation time of 100 trajectories ranged from 8.49 µs to 6.61 ns, with a change in T from 300 K to 1100 K (2.06 µs at 500 K). The data set obtained allowed for calculating the diffusion characteristics with such high accuracy that the effect of crystal strains on them was made significant, even at the 0.1% strain.
The tracer diffusivity D tr (self-diffusivity per one SIA) was calculated by the Einstein-Smoluchowski relation [25], using the initial and end positions of the model crystallite atoms for each of the 100 program runs simulating the di-interstitial diffusion. For each dataset, the self-diffusivity intermediate values were calculated. The final D tr value was obtained by their averaging. The D tr error was estimated as the standard error for the specified set of 100 intermediate values.
The correlation factor f tr was calculated as the ratio of D tr to D. Since the value of f tr is sensitive to the defect diffusion mechanism, it can be used to indicate changes in the di-interstitial diffusion mechanism with changes in temperature or cluster size. For example, in the case of a one-dimensional diffusion of crowdions, f tr = 0 [26], whereas in the case of a three-dimensional diffusion of 100 dumbbells under Johnson's mechanism, mboxemphf tr = 0.44 [27]. The di-interstitial diffusion mechanism can be mixed when the di-interstitial makes a specific number of jumps along one of the directions 110 , subsequently changing this migration direction.

Specific Features of the Diffusion Tensor in FCC Crystals
During the diffusion of a di-interstitial in an FCC Cu crystal, its constituent SIAs usually make jumps to the first nearest neighbors (NN). Since the share of jumps to the second NN is a hundredths of a percent of the total number at the considered temperatures and strains, it was further assumed that the di-interstitial makes jumps only to the first NN. If the di-interstitial makes a random walk with probability n k of making a jump in the direction k (k = 1, 2 . . . 6 [011], respectively, n 1 + n 2 + . . . + n 6 = 1), the normalized diffusion tensor D ij in the CCS takes the form: One can see from (6) that Tr D ij = 3, while the absolute values of the components cannot exceed 3/2. Figure 2 shows the MD-calculated dependences D(T) and f tr (T) for di-interstitials in Cu in the absence of external elastic fields (here and hereafter, the data inaccuracy in the figures is not indicated, since it was smaller than the size of the symbols presented on the graphs). The dependence D(T) was non-linear in Arrhenius coordinates (Figure 2a): the parameters of Arrhenius approximations at low and high temperatures differed ( Table 2). Despite the low activation energy, which is usually inherent in one-dimensional diffusion of self-point defects in metals, di-interstitial diffusion is essentially three-dimensional: f tr values are markedly different from zero and increase with temperature ( Figure 2b) from 0.17 at 300 K to 0.36 at 700 K. a specific number of jumps along one of the directions 110, subsequently changing this migration direction.

Specific Features of the Diffusion Tensor in FCC Crystals
During the diffusion of a di-interstitial in an FCC Cu crystal, its constituent SIAs usually make jumps to the first nearest neighbors (NN). Since the share of jumps to the second NN is a hundredths of a percent of the total number at the considered temperatures and strains, it was further assumed that the di-interstitial makes jumps only to the first NN. If the di-interstitial makes a random walk with probability nk of making a jump in the direction k (k = One can see from (6) Figure 2 shows the MD-calculated dependences D(T) and f tr (T) for di-interstitials in Cu in the absence of external elastic fields (here and hereafter, the data inaccuracy in the figures is not indicated, since it was smaller than the size of the symbols presented on the graphs). The dependence D(T) was non-linear in Arrhenius coordinates (Figure 2a): the parameters of Arrhenius approximations at low and high temperatures differed ( Table 2). Despite the low activation energy, which is usually inherent in one-dimensional diffusion of self-point defects in metals, di-interstitial diffusion is essentially three-dimensional: f tr values are markedly different from zero and increase with temperature ( Figure 2b) from 0.17 at 300 K to 0.36 at 700 K.  Table 2).  Table 2).  Figure 3 shows the MD-calculated dependences D(ε α ), where α = 1, 2, 4, for the diinterstitials in Cu at 500 K. These dependences were approximated well by the expressions:

Diffusivity
where D(0) = 1.33·10 -4 cm 2 /s. The dependences D(ε 2 ) and D(ε 4 ) are even, since due to the crystal's symmetry, its rotation by 90 • around the axis [001] is equivalent to a change in the sign of the strain state to the opposite, while the crystal rotation cannot change the properties of its defects.   Figure 3 shows the MD-calculated dependences D(εα), where α = 1, 2, 4, for the diinterstitials in Cu at 500 K. These dependences were approximated well by the expressions: where D(0) = 1.33·10 -4 cm 2 /s. The dependences D(ε2) and D(ε4) are even, since due to the crystal's symmetry, its rotation by 90° around the axis [001] is equivalent to a change in the sign of the strain state to the opposite, while the crystal rotation cannot change the properties of its defects.  (Table 3). For the MD-calculated dependences

Normalized Diffusion Tensor D ij
For ε ij = 0, D ij = δ ij within the calculation inaccuracy for the considered T. Strain state V 1 does not change the symmetry of the crystal, so D ij = δ ij must also be true. Strain states V 2 and V 4 lower the symmetry of the crystal, which should lead to a change in D ij , imposing certain relations for its components (Table 3). For the MD-calculated dependences D ij (ε α ), these relations are actually fulfilled within the framework of the calculation inaccuracy (Figures 4 and 5). Table 3. Relations for the components of the normalized diffusion tensor, D ij (i, j = 1, 2, 3), under strain states V α (α = 1, 2, 4) following from the crystal symmetry.

Elastodiffusion and Dipole Tensors
The eigenvalues of the elastodiffusion tensor d (α) can be obtained by adjusting their values in such a way that expression (4) would describe the calculated dependences

Elastodiffusion and Dipole Tensors
The eigenvalues of the elastodiffusion tensor d (α) can be obtained by adjusting their values in such a way that expression (4) would describe the calculated dependences D ij (ε α ) in the region of the small ε α , in which the linear nature of such dependences is observed. For the specified region, the expressions following from (4) are valid: Using Equations (10)- (12) to adjust the MD data presented in Figures 3 and 5 in the region |ε α | ≤ 0.1% gave, for d (α) /(D(0)β), the values 4.945 eV, -2.693 eV, 14.085 eV for α = 1, 2, 4, respectively.
In [1], it was analytically shown that the eigenvalues of the elastodiffusion tensor in cubic crystal system crystals could be determined if the dipole tensors of the defect in the stable P e ij and saddle-point P s ij configurations are known. Table 4 shows the relations obtained in [1] connecting d (α) with P e ij and P s ij for different symmetries of the defect saddle-point configurations. Table 4. Eigenvalues of the elastodiffusion tensor for various symmetries of the defect saddle-point configurations.

Symmetry d (1) /(D(0)β) d (2) /(D(0)β) d (4) /(D(0)β)
Cubic Di-interstitial diffusion is characterized by a combination of different diffusion mechanisms; hence, a di-interstitial has multiple different saddle-point configurations. However, using the relations in Table 4, it is possible to obtain the dipole tensor of the effective saddle-point configuration to which all real saddle-point configurations contribute (each with its own weight).
All three calculated eigenvalues d (α) were nonzero for Cu, which corresponded, according to Table 4, to the orthorhombic symmetry of the effective saddle-point configuration of the di-interstitial. Therefore, for the components P s ij , the following expressions are true:

Discussion
Knowledge of the dipole tensor of the defect saddle-point configuration allows calculation of its diffusion tensor for any strain state and simulation of the RD anisotropic diffusion in mechanical fields, generated by external crystal loads and/or internal sources (dislocations, sub-boundaries, etc.), by the object kinetic Monte Carlo (OKMC) method, taking into account both the effects of elastic crystal anisotropy and crystal symmetry [28]. The probabilities n k used in the OKMC method can be represented as [28]: where E int k is the interaction energy of the defect saddle-point configuration for the migration direction k, with the strain field ε ij calculated by Equation (5). Applying Equation (14) to (6), and introducing the designation x = (P s 33 − P s 11 )βε 2 , y = 2P s 12 βε 4 , we obtain for non-zero components D ij under strain states V 2 and V 4 : Figure 5 shows the dependences (15) and (16) and the corresponding MD dependences. Equations (15) and (16) that use the dipole tensor as a parameter described the MD data with high accuracy for all the considered strain values, including in the region of nonlinear dependence D ij on strains, which favorably distinguishes them from expression (4) that uses the elastodiffusion tensor as a parameter, which accurately described the MD data for |ε 4 | ≤ 0.001 only. Therefore, the dipole tensor of the effective saddle-point configuration of the di-interstitial for FCC Cu obtained in this study allowed for simulating of di-interstitial anisotropic diffusion in different elastic fields with high accuracy.

Conclusions
A new method for calculation of the dipole tensor of the saddle-point configuration of radiation defects was proposed. In this method, the molecular dynamics method was used to calculate the strain dependences of the RD diffusion tensor for homogeneous elastic fields of different types. These dependences were used to determine the RD elastodiffusion tensor. The crystal symmetry and the RD diffusion mechanism determined the tensor symmetry. This corresponded to the dipole tensor of the effective RD saddle-point configuration, which takes into account the contributions of all real saddle-point configurations.
The proposed approach was used for studying di-interstitial diffusion characteristics in copper. The symmetry of the calculated dipole tensor of the effective saddle-point configuration was orthorhombic. For all the strain states considered, the effect of the external elastic strain field on the MD-calculated normalized di-interstitial diffusion tensor (ratio of diffusion tensor to a third of its trace) was fully consistent with theoretical expressions, the parameters of which were the components of the dipole tensor, including those outside the range of linear dependence of the diffusion tensor on strains.
The results obtained allowed for simulation of di-interstitial anisotropic diffusion in external and internal elastic fields and for taking into account the di-interstitial contribution to crystal radiation deformation. This contribution is significant, since, as MD data on the primary radiation damage in copper showed,~20% of self-interstitial atoms produced by cascades of atomic collisions are combined into di-interstitials.
The equations given in the present work are ready to use to study the defect diffusion characteristics in FCC and BCC crystals. However, one can apply the proposed approach to crystals of systems other than the cubic one.