First principle calculations of pressure dependent yielding in solute strengthened aluminium alloys

The pressure dependence of the yield stress in solute strengthened aluminium alloys is investigated by first principle calculations. The solute elements studied are magnesium, silicon and copper. A fixed boundary cluster model is employed to calculate the interaction energies between the edge dislocation and the solutes, while simultaneously controlling the hydrostatic pressure in the system. The results show a systematic increase in yield stress with increasing hydrostatic pressure for all solute elements. The calculated pressure dependence is in qualitative agreement with experiments, but underestimated quantitatively. It is suggested that the experimentally observed pressure dependence is caused by both the static and the transient dilatancy of dislocations. In contrast to magnesium and copper atoms, silicon atoms are found to interact non-elastically with dislocations within the core field, indicating that the favourable position for the silicon atoms is in the distorted sites in the


Introduction
Atomistic scale investigations are essential tools in understanding the behaviour of dislocations. Dislocation motions and interactions are the cause of plastic flow and work hardening in metals and alloys. Simulations of these phenomena can help develop material models for properties and structural analysis in both academia and industry, and aid in calibrating such models. Further development of an experimentally validated framework for multi-scale material modelling would reduce the need of time-consuming and expensive testing. An accurate material model must be able to account for all relevant phenomena in order to be robust and reliably used in industry. One particular effect often neglected in higher scale constitutive models for aluminium alloys is the strength-differential (SD) effect. The SD effect signifies a difference in flow stress depending on whether the material is subjected to tensile or compressive loading. This effect has been studied since the 70s, when it was first observed in high-strength steel [1][2][3]. Later, it was shown to also be relevant for aluminium alloys, both for the 1xxx-series [4] and the age-hardenable 6xxx-series [5].
Several hypotheses for the mechanism behind the SD effect were tested by early pioneers on the topic [1,6]. The SD effect was initially believed to be an artefact caused by the experimental setup, such as friction in compression testing. This was eliminated through thorough investigations, which resulted in various other hypotheses [7], namely microcracking due to quenching and residual stresses from prior processing of the material. The former was ruled out after the same effect was observed in materials resistant to microcracking [1]. The latter failed on the account that the SD effect should disappear after cyclic loading or at plastic strain of a few percent [8]. Two strong hypotheses remained: non-linear elastic effects, as discussed by Hirth and Cohen [1], and a permanent volume change during plastic deformation studied by Drucker [6]. The latter was investigated by Spitzig and Richmond [4,9] and disregarded on the basis that this hypothesis predicts a significantly higher volume expansion than measured.
Spitzig and Richmond [4] also investigated if the SD effect was caused by the pressure sensitivity of dislocation motion, which is a nonlinear elastic effect. Two dislocation models were compared, one by Shmatov [10] and another by Jung [11]. The study yielded good agreement between theories and predictions of both the volume change and the influence of the hydrostatic pressure on the shear modulus and the shear strength. The volume change was predicted through a change in the dislocation density. The effect on the shear modulus was argued to be a result of the influence of the hydrostatic pressure on the dislocation motion. Based on their results, Spitzig and Richmond [4] proposed that the flow stress, expressed in terms of the von Mises equivalent stress eq , depends linearly on the hydrostatic pressure P according to where 0 is the yield stress at = P 0 and is a constant governing the pressure dependence. Bulatov et al. [12] investigated the dislocation motion in pure aluminium by molecular dynamics (MD), and, in particular, the effect of hydrostatic pressure on the transient dilatancy caused by the movement of a dislocation. The MD calculations yielded -values in good agreement with the value reported by Spitzig and Richmond [4]. The SD effect in age-hardenable aluminium alloys was investigated by Holmen et al. [5]. They conducted finite element modelling in combination with experiments, and found good agreement between simulations and experiments when adopting the -value reported in [4] for aluminium. The study by Holmen et al. illustrates the importance of the SD effect in high strength aluminium alloys, and further strengthen the hypothesis proposed by Spitzig and Richmond [4].
However, the fundamental mechanism of the SD effect is still under discussion. It is not known if there exists a pressure dependence of the interaction between dislocations and other defects, and whether it is important for alloy systems. In this study, density functional theory (DFT) is applied to investigate the effect of hydrostatic pressure on the solute-dislocation interaction energy and the yield stress of solutestrengthened aluminium alloys. The influence of hydrostatic pressure on the static dilatation of an edge dislocation in the presence of solute atoms has been studied, and a difference in yield stress during compression and tension has been observed, in qualitative agreement with literature.

Background
Solute strengthening is the increase in strength due to the impeding of dislocation motion by solute elements. The mechanism is governed mainly by the interactions between the elastic mismatch of solute elements and the pressure field of dislocations. This interaction increases the energy barrier for dislocation mobility. Solute strengthening can be treated by classical models through elasticity theory, in which case the region nearest to the dislocation is omitted as to avoid the divergency in the stress field. Leyson et al. [13] showed that the atoms close to the dislocation have a large contribution to the strengthening, and thus they are crucial for an accurate model. An approach combining elastic and atomistic models is capable of taking both the classical and quantum effects into account. The numerical model employed in this study is similar to the model of Leyson et al. [13], which is discussed and explained in detail in the review article by Varvenne et al. [14]. The next portion of this section will serve as a brief summary of the key components and properties of the model.
The model is based on the idea that a statistical distribution of solutes interacts with the pressure field of a dislocation. An infinitely long straight edge dislocation, placed in a statistical solute configuration, will yield a total energy of the system above the equilibrium value. To minimise the energy of the system, the dislocation is allowed to bow out to gain a more favourable solute configuration at the cost of the energy gained through excess line tension. This is illustrated in Fig. 1, where connected line segments of length have relaxed in energy wells. The energy wells are spaced by a bow-out amplitude w. The change in the total energy as a function of the length of the line segments and the bow-out amplitude can be expressed by [14] where is the line tension, L is the length of the dislocation line, and b is the length of the Burgers vector. The ratio L/2 represents the number of pinned segments of length . By assuming a dilute binary alloy, the key mesoscale quantity for solute strengthening E w ( ) p can be approximated by where c is the solute concentration and U x y ( , ) i j is the solute-dislocation interaction energy for a solute placed at position x y ( , ) i j . The coordinates x i and y j are defined to be orthogonal to the dislocation line and to each other, and the origin is placed at the centre of the dislocation. The positive direction of x i is defined to be parallel to the glide direction of the edge dislocation. The sum in Eq. (3) is taken over the whole space, but scales as 2 . An energy minimisation of Eq.
(2), with respect to w and , will yield the characteristic values c and w c of the segment length and bow-out amplitude that are specific for the system. Note that the minimisation of Eq. (2), with respect to , yields a characteristic value w ( ) c , which depends on the bow-out amplitude. The energy barrier to unpin a 2 c -long dislocation segment is expressed as [14] where the prefactor 1.22 is a result of the numerical constants emerging during the derivation. For the dislocation to glide, it has to be affected by an applied resolved shear stress. The stress matching the energy barrier will be the zero Kelvin yield stress, which can be calculated as To compare with experiments, the yield stress at higher temperatures is required. Dislocation motion can be thermally activated, and thus motion can be initiated at lower shear stress at higher temperatures. From transition state theory, assuming a quasi-static loading, the finite temperature yield stress can be approximated by [14] = k T E 1 ln , where k b is the Boltzmann constant, T is the temperature, and 0 and are the reference and macroscopic strain rates, respectively. Assuming a polycrystalline sample during a tension or compression test, the shear yield stress must be multiplied by a Taylor factor to account for the bulk behaviour. The Taylor factor is dependent on the texture and loading direction, and a Taylor factor of 3 is generally accepted for random texture, and enables the calculated stress to be comparable to the yield stress found experimentally. Fig. 1. The solid red line represents a dislocation in a favourable energy configuration after bow out. The dashed line represents the straight dislocation before bow out. The red circles represent solutes that have a positive interaction energy with the dislocation, while blue circles represent solutes with a negative interaction energy. w is the bow-out amplitude, and is the length of a line segment.
An expression for the solute-dislocation interaction energy is needed in order to calculate the characteristic bow-out amplitude from Eq. (2) and (3). The interaction energy is well defined by elasticity theory far away from the dislocation core, and can be approximated by where V is the misfit volume associated with the solute and p x y ( , ) i j is the pressure field of the dislocation. The pressure field of a straight dislocation is given by the Volterra solution through anisotropic elasticity [15], where is the Poisson ratio, and µ is shear modulus. It is evident from Eq. (8) that the pressure field is divergent when approaching the origin, i.e. the centre of the dislocation. The elasticity approach is inaccurate close to the dislocation, which is why the solutedislocation interaction energy in this region must be calculated using an atomistic approach. In an atomistic approach, different atomistic models are used to isolate the interaction. Considering a system consisting of aluminium matrix, a dislocation at the origin and a solute at position x y ( , ) i j , the calculation is done by where E tot is the total energy of the full system, E dis is the energy of having a dislocation in the matrix including the free surface, and E sol is the energy difference of substituting a solute in the matrix, see Fig. 2.
The atomistic models used in the DFT calculations are discussed in Section 3.

Method and model
The DFT calculations were performed using the Vienna Ab Initio Simulation Package (VASP) [16,17]. The functional used for the calculations is the Perdew-Burke-Ernzerhof generalized gradient approximation [18], which has shown to be accurate for metals and is relatively fast computationally. The plane-wave energy cut-off was 400 eV and a gamma sampling of 0.18 k-points per Å was used to model the Brillouin zone.
The bulk fcc aluminium lattice constant has been calculated for multiple levels of hydrostatic pressure. A change in the hydrostatic pressure P will result in a volumetric strain, represented by a minuscule change in the lattice constant. This is resulting in a linear relation between the hydrostatic pressure and the lattice constant, which was calculated in a preliminary study, see Fig. 3. In these calculations, two different sized supercells were used to minimise the numerical error.
Modelling of dislocations with periodic plane wave DFT has an inherit challenge due to the asymmetric nature of the long-range stress field generated by the dislocation. With periodic boundary conditions, the system has to be large enough for the stress field to be insignificant at the boundaries, or the model must be designed in such a way that the stresses at the boundaries are cancelled. Alternatively, the system could be surrounded by vacuum in two directions, creating a so-called cluster model. The surface of the cluster model must be corrected by a fixed or a flexible boundary condition to represent an infinite matrix [19]. Such a model would circumvent the problems of the long range displacement field of the dislocation, but at the cost of introducing an artificial surface. Applying a fixed boundary condition adds the benefit of having a controllable hydrostatic pressure in the system. The drawback of a fixed boundary condition is a constrained displacement field in the centre of the dislocation, referred to as the core field. The approximation is justified by a faster decay of the non-elastic part of the displacement field compared to the elastic part.
The cluster model adopted in this work is illustrated atomistically in Fig. 4. The cluster is divided into an outer and inner region, represented by purple and grey colours, respectively. The radius of the inner region is R 1 , and the outer region is between R 1 and R 2 . The model is periodic in the [112] direction, parallel to the dislocation line.
The centre of the cluster is referred to as the centre-of-symmetry (COS), represented by a green large in i j tot is the total energy of the full system, E dis is the energy introduced by the dislocation and matrix, while E sol is the energy introduced by the solute.  The displacement field is given by the Volterra solution from anisotropic elasticity theory [15], which was imposed to all atoms in the model. The atoms in the outer region of the model were constrained during the relaxations, under the assumption that these atoms were sufficiently far away from the partial dislocations so that the non-elastic part of the displacement field is negligible. Relaxations with different R 1 values, R [10,22] 1 Å, showed a maximum displacement of any atom of less than 0.05 Å near the interface between the fixed and the relaxed region. The value chosen for R 1 puts a limit to how many positions there are for solutes in the relaxed region, as the strain field of solutes are excluded in the fixed region. The inner and outer radii were set to = R 14 1 Å and = R 20 2 Å for systems containing Mg solutes, while for Si and Cu containing systems both radii were increased by 6 Å due to the higher selection space for solute placement needed in the analysis. In the periodic [112] direction, a slab thickness of one unit cell length, 4.96 Å, was used. The full atomistic model contains 371 and 640 atoms for the small and large version, respectively. The minimum width of the surrounding vacuum was set to be 10 Å, making the interactions with periodic images negligible.
The initial relaxation of the inner region was performed for systems with different hydrostatic pressures P ranging from −600 MPa to 600 MPa. The pressure was superimposed by changing the lattice constant in accordance to Fig. 3.
The complete system consists of two partial dislocations constituting a pure edge dislocation at the COS, and one solute substituted at different aluminium sites for different calculations. The various sites used for substitution are indicated in the figures in Section 4. Due to the mirror symmetry, only one side of the COS was sampled. This assumption was confirmed in preliminary calculations. The interaction energy was calculated as described by Eq. (9) and Fig. 2.
E sol can be calculated using two methods. The first method relies on an assumption that the solute-dislocation interaction energy is antisymmetric over the tension-compression boundary of the dislocation stress field. As a consequence, the interaction energies of systems with an even number of solutes, sampled at both sides of the boundary, will be cancelled out. E sol can then be found by subtracting the energy introduced by the dislocation from the total energy of the systems and dividing by the number of solute sites sampled. The second method is finding the energy difference between two systems, one of pure aluminium and one where an aluminium atom is substituted with another element. The latter method should be performed with a thin simulation cell to reproduce the system geometry in Fig. 4. The numerical implications of the two methods are discussed in Section 4.
The parameters used to calculate the yield stress for systems at different hydrostatic pressure are summarised in Table 1. The numerical value of the solute concentration has been chosen in accordance with previous experimental and numerical studies. The results are presented in the next section.

Results and discussion
The contour plots of the solute-dislocation interaction energy for Mg and Cu, at P 0, are presented in Fig. 5. The most energetically favourable position for the Mg solutes is in the lowest density regions of the system, i.e. the tensile area of the dislocation. With Mg solutes, the highest energy is found when the solutes are placed above the partial dislocation cores, in the compressed area. The energetically favourable positions of the Mg solutes can be predicted by elasticity theory, since they are influenced by the difference in misfit volume of the elements. Mg has a positive misfit energy, see Table 1, resulting in a compressive strain in the matrix, while Cu and Si have a negative misfit volume. Thus, the interaction energy is expected to be qualitatively reversed over the tension-compression boundary for Cu and Si.
Although this is the case for Cu, see Fig. 5b, the Si solutes seem to be expressing a non-elastic effect manifested through an increase in the interaction energy in the stacking fault, see Fig. 6a. The preferred position of a Si atom is close to the partial dislocation cores, where the lattice is most displaced. There is also a localised energy plateau approximately five sites above the centre of symmetry, where the interaction energy turn positive. For Cu solutes a non-elastic effect is seen on the tensile side through an increased interaction energy in the stacking fault, but a change from negative to positive interaction energy on the compressive side is not observed. The change of sign in the interaction energy for Si was first suspected to be an artefact caused by the small size of the system, = R 14 1 Å. However, the non-elastic behaviour was as prominent when repeating the calculations with an increased radius = R 20 1 Å, indicating that the observation was caused by other factors which will be discussed later. Based on the investigation of the Si system, the R 1 value was increased to 20 Å in the production runs for Si and Cu as a precaution, and for the added benefit of a larger selection space for solute placements valuable for further analysis. Another possible reason for the localised energy plateau in the interaction-energy map could be the slab thickness of one unit cell that results in a continuous column of Si atoms. To check if this is the case, a calculation was performed with the system thickness doubled. In agreement with Layson et al. [13], this results in an insignificant change in the interaction energy of less than 0.01 eV, indicating that the thickness of the slab is not the source of the observed energy plateau.
From the observations, Si behaves differently depending on the position of the nearest neighbours, indicating a chemical effect not captured by elasticity. One likely cause of the plateau is the covalent bonds of Si. The covalent bonds are directional and do not match with the positions of the Al sites in the fcc matrix, nor in the stacking fault. Owing to the displaced nature of the dislocation core, it can potentially be a more suitable position for the Si atoms, thus making them more tightly bonded. Charge-density difference maps of different Si sites were investigated, but did not yield any conclusive results. The observations in the interaction energy map motivate further studies on the interaction between solutes and dislocations, and the interaction between solutes and distorted matrices. The increased interaction energy away from the core of the partial dislocations may affect the local Table 1 Numerical parameters used in the analysis of solute-strengthening in aluminium binary alloy: c is the solute concentration, V is the misfit volume of the respective solute, where values are calculated by Leyson et al. [20], is the dislocation line tension energy, calculated by Dong et al. [21] and used by Leyson et al. [20], is the chosen strain rate, and and µ are Poisson's ratio and the shear modulus, respectively. kinetics of solutes. This can in turn have an influence on other phenomena, such as dynamic strain ageing where cross-core diffusion is a possible mechanism [22]. The non-elastic effect on the interaction energy in Fig. 6a had an effect on the value of the key mesoscale quantity E p , and thus on the value of the characteristic bow-out amplitude w c . To explore the magnitude of this effect, different samplings of solute placement were used when calculating w c . Fig. 6b shows an example of an alternative solute sampling, where the rows nearest the dislocation are sampled. As seen in the figure, the inelastic behaviour is not well captured due to the low sampling set. However, it was found that the position of the peaks in the interaction energy on the nearest rows, caused by the dissociated dislocation, plays a more critical role on w c .
In addition to the solute sampling near the core region, the choice of interpolation of the interaction energy also influenced the calculations of w c . Two different methods of interpolation were specifically analysed. The interaction-energy map was interpolated by linear interpolation of all atoms, or in a two-step method. The two-step interpolation was first done by interpolating the interaction energy on the tensile and compression sides separately, and then the two interpolated fields were merged. The idea is that interpolation over the discrete jump from a negative to a positive interaction energy is problematic. The difference in the resulting w c of the two interpolation procedures was noticeable. That said, the solute-strengthening model by Leyson et al. [20] has shown to correctly take the interaction energy in the near vicinity of the dislocation core into account, and these solute atoms The black and white crosses are the partial dislocations and the centre of symmetry, respectively. The interaction energy inside the black circle, enclosing the relaxed region of the system, was found by linear interpolation between the interaction energy outside the black circle, calculated from elasticity theory, and the white dots. Fig. 6. Solute-dislocation interaction energy map for Si, at approximately zero hydrostatic pressure, with different solute sampling, where (a) has a more complete sampling than (b). The white dots represent solute placement where individual DFT calculations have been conducted. The black and white crosses are the partial dislocations and the centre of symmetry, respectively. The interaction energy inside the black circle, enclosing the relaxed region of the system, was found by linear interpolation between the interaction energy outside the black circle, calculated from elasticity theory, and the white dots.
provided the most significant contribution to w c as indicated in their paper. In brevity, a decision was made to stay consistent throughout all the different pressure states of the system, and the results shown here are thus calculated using linear interpolation over the entire system.
The inelastic behaviour of Si brought questions about the calculation of the solute energy, E sol , entering into Eq. (7). A change in the solute energy is seen as a constant change in the interaction energy at all positions, and this could explain the change in sign that can be seen in Fig. 6a. The initial calculations of this quantity were carried out under the assumption that the interaction energy is antisymmetric over the tension-compression boundary of the dislocation stress field, as explained in Section 3. This assumption might not be valid as the solute energy could be biased by a solute tightly bonded, or strongly repelled by the dislocation. The alternative method in which E sol is calculated separately, avoids interfering with a dislocation. However, in this method the fluctuations due to the constrained core field are not corrected. The error of the core field should be relatively low, based on the convergence of the system size. The solute-dislocation interaction energy, with the solute energy calculated separately, can be seen in Fig. 7.
These results indicate that Si should be tightly bonded to dislocations, and not strongly repelled by the compression side to the extent predicted by elasticity theory. This suggests that the core region is energetically favourable for Si atoms to reside, due to a distorted matrix. Neither Mg or Cu showed similar effects, and the solute energy of Mg and Cu did not change considerably based on the calculation method. The interaction energy near a dislocation could be relevant when studying nucleation of precipitates at dislocation lines, where the interaction energy will affect both the local concentrations of solutes and the kinetics. This subject will not be expanded on in this study.
The interaction energy found in this study is in good agreement with previous studies [20], with the exception of the energy map of Si. The solute-strengthening model was used to calculate the yield stress, and the zero pressure yield stress is used to compare with experiments for Mg and Cu. The results are shown in Table 2. Note that Eq. (6) is used to convert the predicted values at 0 K to corresponding values at 78 K. The low yield stress for Cu compared to experimental results is consistent with results obtained by Leyson et al. [20]. The higher values in the experimental results were attributed to the concentration of iron in the Al-Cu alloys under experimental testing [23]. The temperature is set to 78 K for all the presented results regarding yield stress in this study for simplicity and consistency. Fig. 8 and 9 show the calculated yield stress for solute strengthening as a function of hydrostatic pressure for Mg, Cu and Si together with linear interpolations. The yield stress is undoubtedly higher for superimposed compression than for tension. This trend is evident for all cases and in qualitative agreement with experimental studies [4,5]. The uncertainties in the calculated yield stresses are difficult to analyse, since the model relies on the relative values between interaction energy maps with different atom positions. The maximum value of the interaction energy was higher in compression than in tension, although the results are not conclusive due to limited statistics. However, the strength difference found is less pronounced compared to experimental observations. As mentioned, Spitzig and Richmond [4] proposed that the flow stress is linearly dependent on pressure, with a parameter determining the pressure dependence, see Eq. (1). This parameter was said to be material dependent [4], i.e. it should be dependent on the elastic constants of the matrix material, meaning that similar values of are expected for all solute elements. The -parameter values predicted in this study were 14, 20 and 26 TPa −1 for Mg, Si, and Cu, respectively. The results indicate a pressure dependence that is affected by the solute elements. However, there is not enough data to firmly conclude such a pressure dependence, and a more thorough investigation on the subject could yield further insight. One possible explanation is that the hydrostatic pressure affects the solute-dislocation interaction energy differently depending on the solute, because the strain field of each solute atom is unique with different associated volumetric strain.
The predicted values of are consistently lower than what has been found in previous studies, both experimentally [4,5] and in a numerical study by Bulatov et al. [12]. The experimental study by Spitzig and Richmond reported 58.7 TPa −1 [4]. In their MD calculations, Bulatov et al. [12] investigated the influence of hydrostatic pressure on the transient dilatation generated by moving dislocations, but they did not consider any static dilatation, such as solute-dislocation interaction, dislocation-dislocation interaction or interaction with precipitates. In the present work, solute-dislocation interactions are included, and the results indicate that these interactions have a considerable contribution to the strength-differential effect for solute-strengthened alloys. In light of these observations, we propose that the solute strengthening effect is caused by the pressure dependence of both the static and transient dilatation of dislocations. Both effects must be considered when studying the mechanisms behind the SD effect. This study strengthens the hypothesis that the strength differential effect is a non-linear elastic effect as proposed by Spitzig and Richmond [4]. The results presented give a direct indication that pressure has a non-linear elastic effect on the interaction between solutes and dislocations and thus the energy barrier experienced by the dislocation.

Conclusion
The pressure dependence of the yield stress in solute strengthened Al alloys has been investigated by first principle calculations using density functional theory. The studied solute elements were Mg, Cu and Si. The results show a clear increase in yield stress as the hydrostatic pressure is amplified, in agreement with published experimental studies Fig. 7. Solute-dislocation interaction energy map for Si, at approximately zero hydrostatic pressure, with the energy cost of a solute calculated separately. The white dots represent solute placement, where individual DFT calculations have been conducted. The black and white crosses are the partial dislocations and the centre of symmetry, respectively. The interaction energy inside the black circle, enclosing the relaxed region of the system, was found by linear interpolation between the interaction energy outside the black circle, calculated from elasticity theory, and the white dots.  [4]. A cluster model with fixed boundary conditions has successfully been utilised for calculating the interaction energy between solutes and an effective edge dislocation under varying superimposed hydrostatic pressure. The hydrostatic pressure has a visible effect on the solutedislocation interaction energy, manifesting a pressure dependent yield stress for all solute elements studied. The calculated interaction energy between Si and the two Shockley dislocations in aluminium alloys is more complicated than elasticity theory would suggest. In contrast to Mg and Cu solutes, the solute-dislocation interaction energy map for Si atoms indicates that Si is more attracted to the dislocation than repelled. This attraction is not predicted by elasticity theory, and suggested here to be a chemical effect. A more thorough investigation could yield valuable information about the local kinetics in systems containing solutes and other defects.

Data availability
The raw data used to reproduce the presented results are available in the Zenodo repository, https://doi.org/10.5281/zenodo.3686913 [24].

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.