A First-Principles Study of the Cu-Containing β″ Precipitates in Al-Mg-Si-Cu Alloy

The nanostructured β″ precipitates are critical for the strength of Al-Mg-Si-(Cu) aluminum alloys. However, there are still controversial reports about the composition of Cu-containing β″ phases. In this work, first-principles calculations based on density functional theory were used to investigate the composition, mechanical properties, and electronic structure of Cu-containing β″ phases. The results predict that the Cu-containing β″ precipitates with a stoichiometry of Mg4+xAl2−xCuSi4 (x = 0, 1) are energetically favorable. As the concentration of Cu atoms increases, Cu-containing β″ phases with different compositions will appear, such as Mg4AlCu2Si4 and Mg4Cu3Si4. The replacement order of Cu atoms in β″ phases can be summarized as one Si3/Al site → two Si3/Al sites → two Si3/Al sites and one Mg1 site. The calculated elastic constants of the considered β″ phases suggest that they are all mechanically stable, and all β″ phases are ductile. When Cu atoms replace Al atoms at Si3/Al sites in β″ phases, the values of bulk modulus (B), shear modulus (G), and Young’s modulus (E) all increase. The calculation of the phonon spectrum shows that Mg4+xAl2−xCuSi4 (x = 0, 1) are also dynamically stable. The electronic structure analysis shows that the bond between the Si atom and the Cu atom has a covalent like property. The incorporation of the Cu atom enhances the electron interaction between the Mg2 and the Si3 atom so that the Mg2 atom also joins the Si network, which may be one of the reasons why Cu atoms increase the structure stability of the β″ phases.


Introduction
Heat treatable Al-Mg-Si(-Cu) alloys in the 6xxx series are a common category of structural materials used in the construction and transportation industries. These alloys can be customized to have a desirable combination of properties, such as good formability, high specific strength, and corrosion resistance [1][2][3]. After proper aging treatment, the strength of the alloy can be greatly improved. This is mainly due to the precipitates that can contribute to the strengthening mechanisms by hindering the dislocation movement [4,5], particle strengthening σ p [6], and coherency of the particles [7]. The mechanical properties of these alloys can be greatly influenced by the composition, morphology, scale, and distribution of these solute atom nanostructures [8]. The precipitation sequence for Al-Mg-Si alloys is generally considered to be [9,10]: Among the precipitates [4,[13][14][15] formed in the aged Al-Mg-Si alloys, needle-like β precipitate is the most effective strengthening phase [14] responsible for the peak-hardening effect [16]. The β phase is a metastable precipitate phase, which is semi-coherent with the Al matrix in the needle cross-section, the space group is C2/m, a = 15.16 Å, b = 4.05 Å, c = 6.74 Å, and β = 105.3 • [17,18]. The monoclinic β phase was originally proposed to have the composition of Mg 5 Si 6 [17]. However, according to recent experimental and theoretical studies, the composition of β would fluctuate around Mg 5 Al 2 Si 4 [19][20][21][22]. Furthermore, the most recent density functional theory (DFT) calculations inferred very minor formation enthalpy differences for β -Mg 5+x Al 2−x Si 4 (−1 < x < 1) [21]. These results indicate that the composition of β phase in Al matrix may change under certain conditions. For example, the dispersed nano-precipitates can be affected by the addition of Mg and/or Si, as well as other elements like Cu [5,[23][24][25][26][27]. The addition of Cu is demonstrated to increase the age-hardening response, and it promotes the generation of higher number density and smaller size precipitates [14,[28][29][30][31][32]. Therefore, a certain amount of Cu is usually added into Al-Mg-Si alloys. The addition of Cu increases the complexity of the precipitation sequence [32,33]. The precipitation sequence of Al-Mg-Si-Cu alloys is reported as [34]: SSSS → solute clusters → GP − zones → β , L/S/C, QP, QC → β , Q → Q, Si Previous work used various experimental and theoretical methods to study the incorporation of Cu in β , and analyzed the Cu atoms as foreign solute atoms in the phases [20]. Cu addition could further enhance the positive effect of pre-aging on bake hardening for Al-Mg-Si alloys [35]. It has been demonstrated by high-angle annular dark-field scanning transmission electron microscopy (HAADF-STEM) that Cu is mainly confined to the Si3/Al sites (Si or Al atoms completely occupy) of the β structure [26,[35][36][37], which as mentioned was also supported from DFT-based calculations [38]. The β precipitates in Al-Mg-Si-Cu alloy were detected with an average composition of 28.6Al-38.7Mg-26.5Si-5.17Cu (at. %) using atom probe tomography (APT) and high-resolution energy-dispersive X-ray (EDX) mapping [36]. Furthermore, the addition of Cu has no effect on the type of β precipitate, Cu atoms incorporate in β and some of Mg, Si and Al in β unit cell are substituted by Cu atoms [39].
As mentioned above, the β precipitation behavior in Al-Mg-Si-Cu alloys has been investigated using various characterization methods. However, the detailed structures and stabilities are still unclear of Cu-containing β phases in these alloys, and these structural refinements could be supported by first-principles results [40]. In addition, we predict energy-lowering site occupations and stoichiometries of the β phases, where experimental information is incomplete. Understanding the structure of Cu-containing β precipitates is essential to elucidate the precipitation sequence in heat-treatable Al-Mg-Si (-Cu) alloys.
In the present work, first-principles calculations based on density functional theory (DFT) [41] were used to study the Cu-containing β phases. Based on the structural information obtained by experimental methods, first-principles atomistic calculations can provide structural, chemical, and energetic information [40]. A large number of Cucontaining β structures were constructed searching for possible stable configurations and structural stability, kinetic stability, and mechanical stability were also considered. Finally, the characteristics of Cu atoms occupying sites were analyzed through the electronic structure.

Atomic Model
For the β phases, the formation enthalpies and lattice parameters of Mg 4 Al 3 Si 4 , Mg 5 Al 2 Si 4 , Mg 6 AlSi 4 , and Mg 5 Si 6 were computed for each of the models of the crystal structures available in the literature [17,18,21], allowing a critical assessment of the validity of the models. Figure 1 shows four atomic models of the β without Cu. The Wyckoff site information of the energetically most favorable β -Mg 5 Al 2 Si 4 is shown in Table 1 [18,19].

Atomic Model
For the β″ phases, the formation enthalpies and lattice parameters of Mg4Al3Si4, Mg5Al2Si4, Mg6AlSi4, and Mg5Si6 were computed for each of the models of the crystal structures available in the literature [17,18,21], allowing a critical assessment of the validity of the models. Figure 1 shows four atomic models of the β″ without Cu. The Wyckoff site information of the energetically most favorable β″-Mg5Al2Si4 is shown in Table 1 [18,19].  [17,21]. (a) Mg5Si6 from Zandbergen [17]; (b) Mg4Al3Si4 and (c) the Mg5Al2Si4 from Hasting [19]; (d) Mg6AlSi4 from Ehlers [21]. The relative location of each site is marked in (c). Table 1. Wyckoff site information (x, y, z) in the β″-Mg5Al2Si4 phase [18,19]; atomic configuration is shown schematically in Figure 1c.

Computational Details
The first-principles calculations were performed utilize the plane wave pseudopotential method, as implemented in the highly efficient Vienna ab initio simulation package (VASP) [42,43], The electron-ion interactions were described through projector augmented wave (PAW) [44,45]. The exchange-correlation function were constructed by the generalized gradient approximation (GGA) of Perdew-Burke-Ernzerhof (PBE) [46]. All structures were fully relaxed with respect to atomic positions as well as all lattice parameters in order to find the lowest-energy structure. The electron wave function was expanded in plane waves up to a cutoff energy of 450 eV. The β′′phase was represented by a conventional cell with 22 atoms according to the experimental results, and 3 × 12 × 8 Γcentered k-point meshes were employed in the Brilluion zone sampling and generated automatically by following the Monkhorst-Pack sampling scheme [47], while the 3 × 3 × 8 Γ-centered k-point meshes and 1 × 4 × 1 supercells were employed for calculation of "replacement energy" (the detailed definition is explained below). Atoms were relaxed until  [17,21]. (a) Mg 5 Si 6 from Zandbergen [17]; (b) Mg 4 Al 3 Si 4 and (c) the Mg 5 Al 2 Si 4 from Hasting [19]; (d) Mg 6 AlSi 4 from Ehlers [21]. The relative location of each site is marked in (c). Table 1. Wyckoff site information (x, y, z) in the β -Mg 5 Al 2 Si 4 phase [18,19]; atomic configuration is shown schematically in Figure 1c.

Computational Details
The first-principles calculations were performed utilize the plane wave pseudopotential method, as implemented in the highly efficient Vienna ab initio simulation package (VASP) [42,43], The electron-ion interactions were described through projector augmented wave (PAW) [44,45]. The exchange-correlation function were constructed by the generalized gradient approximation (GGA) of Perdew-Burke-Ernzerhof (PBE) [46]. All structures were fully relaxed with respect to atomic positions as well as all lattice parameters in order to find the lowest-energy structure. The electron wave function was expanded in plane waves up to a cutoff energy of 450 eV. The β phase was represented by a conventional cell with 22 atoms according to the experimental results, and 3 × 12 × 8 Γ-centered k-point meshes were employed in the Brilluion zone sampling and generated automatically by following the Monkhorst-Pack sampling scheme [47], while the 3 × 3 × 8 Γ-centered k-point meshes and 1 × 4 × 1 supercells were employed for calculation of "replacement energy" (the detailed definition is explained below). Atoms were relaxed until their residual forces converged to 0.01 eV/Å. The phonon spectra were obtained using the Phonopy package [48].
The four-parameter Birch-Murnaghan equation of state with its linear form [49] is employed to estimate the equilibrium total energy (E 0 ), volume (V 0 ), where a, b, c, and d are fitting parameters. More details can be found in our previous work [50].
Compared with the energy of solid solution containing a Cu atom, the energy gain of the Cu atoms incorporated in β is referred to as "replacement energy". In order to construct the Cu-containing β phases, it is necessary to determine the possible occupation sites of Cu in β phases. Additionally, computing the replacement energy (see Ref. [51]) can be used as a criterion for the possible occupation sites of solute atoms. There have been previous studies addressing the first-principles calculations for describing replacement energies of different sides. Since the replacement energy of Cu atoms at Mg2 and Mg3 sites were not shown in Saito's work [38], one Cu atom was introduced into a 1 × 4 × 1 supercell and the preference of Cu atoms for each non-equivalent site in β was evaluated using the method described by Saito et al. [51], but with higher calculation precision.
To solve the compositional uncertainty preliminarily, the reported C2/m symmetries [18] were deliberately reduced to the level where only pairs of atoms (e.g., the two Cu atoms) were regarded as equivalent. This implies that space group P2/m was used throughout and there are 11 different sites within the unit cell. Besides, no partial occupancies were considered and vacancies were ignored. The replacement energy for Cu incorporation in β can be described as follows: where H are the calculated enthalpy of the system, β 0 are the Cu-free structure, Ξ are the sides in β 0 , X are the solute atoms incorporated in the precipitates, and S are substitutional sites in the Al matrix. A certain atom X incorporates on site Ξ is referred to as "{X → Ξ}". The formation enthalpy of solid solution (SS), ∆H form SS , was used to find out the most energetically favorable configurations in the atomic models. Since there is no stable fcc structure for Mg and Si, their formation energies in relation to SS were determined as follows: where E sub (Mg) and E sub (Si) are the enthalpies of substituting Al atoms by Mg and Si atoms, respectively. E sub (Mg) and E sub (Si) were calculated in a 3 × 3 × 3 Al supercell with one Mg/Si atom and 107 Al atoms with a k-point meshes of 5 × 5 × 5. The enthalpy of substituting a Mg atom was defined as: where E (Al) is the enthalpy of a 3 × 3 × 3 Al supercell. The definition of E sub (Mg) was also feasible for E sub (Si). Finally, in order to compare the structures with different Al content, the formation enthalpy can also be expressed in kJ/mol of solute atoms, instead of kJ/mol [52]. This transformation is achieved as follows: . This is a common definition of formation enthalpy in the literature [9,21,52].
The elastic constant can be represented by a 6 × 6 matrix. Based on the symmetry of the crystal structure, the independent elastic constants of the monoclinic crystal are reduced to 13, as shown in Formula (5): The stress-strain method based on the generalized Hooke's theorem is used to calculate the elastic constants of each crystal [53]. For more detailed stress-strain method description, please refer to [54]. The relationship between elastic constant C ijkl , stress tensor δ kl , and strain tensor δ kl can be expressed as: The Hill model [55] is used to further obtain the bulk modulus (B), shear modulus (G), and Youngs modulus (E) of the crystal through the elastic constant. The Hill model takes into account that the calculation results of the Voigt model and the Reuss model will be high and low, respectively, and take the arithmetic mean of the values of the Voigt model and the Reuss model. For monoclinic crystal structure, the formula for calculating the bulk modulus (B) and shear modulus (G) of monoclinic crystals using Voigt model and Reuss model are [56]: wherein: The formula for calculating the bulk modulus (B), shear modulus (G), and elastic modulus (E) of monoclinic crystal by Hill model [55] is:

Structure Stability
The replacement energy is shown in Figure 2 and alternative solute atoms Mg/Si were incorporated for comparison with Cu at different sites. In order to more intuitively express the competitive occupation sites of Cu atoms, the variable ∆ is introduced and the ∆ values of Cu atoms at different sites in different β configurations are shown in Figure 3. The ∆ represents the "competitiveness" between Cu atoms and other solute atoms at each site, it is the difference between the lowest replacement energy of Mg/Si solute atoms and the replacement energy of Cu atoms. The larger the value of ∆, the more likely the Cu atom will occupy the site. Consequently, due to the low Cu occupancy in β , only three designated Cu sites (Si1, Si3, Mg1, see Figure 1c) were allowed to host Cu atoms according to the relative value of replacement energy (refer to Figure 2). This conclusion is consistent with previous research [38].

Structure Stability
The replacement energy is shown in Figure 2 and alternative solute atoms Mg/Si were incorporated for comparison with Cu at different sites. In order to more intuitively express the competitive occupation sites of Cu atoms, the variable Δ is introduced and the Δ values of Cu atoms at different sites in different β″ configurations are shown in Figure 3. The Δ represents the "competitiveness" between Cu atoms and other solute atoms at each site, it is the difference between the lowest replacement energy of Mg/Si solute atoms and the replacement energy of Cu atoms. The larger the value of Δ, the more likely the Cu atom will occupy the site. Consequently, due to the low Cu occupancy in β″, only three designated Cu sites (Si1, Si3, Mg1, see Figure 1c) were allowed to host Cu atoms according to the relative value of replacement energy (refer to Figure 2). This conclusion is consistent with previous research [38].   For checking the reliability of the calculations, Table 2 displays the structural parameters for selected β″ configuration without Cu atom, along with the results of earlier theoretical and experimental studies of β″. Available calculation results of formation enthalpies are shown in Table 3. The formation enthalpies of the 33 possible unit cells have been plotted in Figure 4, including the configuration without Cu atom. Since the given formation enthalpy of per solute atom (eV/solute atom) essentially presents the solute chemical potentials, the zero-temperature convex hull can be constructed to deduce the precipitation order of the system [57]. It can be seen that Cu occupying one column of each Si3 column pair is found to be the energetically most favorable option for the set of Mg4Al2CuSi4 compositions. While the formation enthalpy of Mg4Al2CuSi4 is −0.337 eV/solute atom, the formation enthalpy of Mg5AlCuSi4 is −0.335 eV/solute atom, which is similar For checking the reliability of the calculations, Table 2 displays the structural parameters for selected β configuration without Cu atom, along with the results of earlier theoretical and experimental studies of β . Available calculation results of formation enthalpies are shown in Table 3. The formation enthalpies of the 33 possible unit cells have been plotted in Figure 4, including the configuration without Cu atom. Since the given formation enthalpy of per solute atom (eV/solute atom) essentially presents the solute chemical potentials, the zero-temperature convex hull can be constructed to de-duce the precipitation order of the system [57]. It can be seen that Cu occupying one column of each Si3 column pair is found to be the energetically most favorable option for the set of Mg 4 Al 2 CuSi 4 compositions. While the formation enthalpy of Mg 4 Al 2 CuSi 4 is −0.337 eV/solute atom, the formation enthalpy of Mg 5 AlCuSi 4 is −0.335 eV/solute atom, which is similar to that of Mg 4 Al 2 CuSi 4 . This is consistent with the observed in previous experiments that Cu atoms mainly occupy Si3 sites [36]. The energy gained when replacing Mg/Si/Al with at the Wyckoff sites is clearly varying with x. When Cu atoms occupy two sites (that is, x = 2), Mg 4 AlCu 2 Si 4 is the energetically most favorable phase, and Cu atoms occupy two Si3 columns. When Cu atoms occupy three sites, Mg 4 Cu 3 Si 4 is the most stable structure, in which Cu atoms occupy one Mg1 site and two Si3 sites, which is consistent with experimental observations [36]. The results show that stoichiometry of Cu-containing β phase is suggested as Mg 4 Al 3−x Cu x Si 4 (1 ≤ x ≤ 3). Since the formation enthalpy of Mg 5 AlCuSi 4 is very close to that of Mg 4 Al 2 CuSi 4 , it can also be taken into account. This result emphasizes the possibility of fluctuations between various compositions as a function of the local alloying element concentration for the physical system during precipitated phases growth. Then the structural parameters of low energy configurations from Figure 4 also have been displayed in Table 2. As discussed above, sole minimization of the β phase formation enthalpy supports the well-defined Mg 4+x Al 2−x CuSi 4 (x = 0, 1) unit cell shown in Figure 5.

Elastic Properties
Here, we compare the mechanical properties of β″ with or without Cu atoms. The                  Table 3. Formation enthalpies of β phases with different configurations in this work. The Cu occupied sites and its number are also listed in detail.  . Formation enthalpies of calculated β″ phases. When Cu atoms occupy x sites the number of Cu atoms in unit cell is corresponding to 2x. The lower energy configuration is marked with a specific shape, and Mg5Si6 is also marked as a reference. represents the Mg5Si6, represents the Mg5Al2Si4, represents the Mg5AlCuSi4 where Cu atoms occupy a Si3 site, represents the Mg4Al2CuSi4 where Cu atoms occupy a Si3 site, represents the Mg4AlCu2Si4 where Cu atoms occupy two Si3 sites, and represents the Mg4Cu3Si4 where Cu atoms occupy a Mg1 site and two Si3 sites.

Elastic Properties
Here, we compare the mechanical properties of β″ with or without Cu atoms. The elastic constants of key β″ phases that are most likely to precipitate during aging were calculated by using fully relaxed crystal structures, and the results are listed in Table 4. According to the Born stability criterion [61], the elastic constants of Mg4Al2CuSi4 and Mg5Al-CuSi4 all meet the stability criteria of monoclinic crystals. This further supports the stability of Mg4+xAl2−xCuSi4 (x = 0, 1) obtained from the formation enthalpy. The elastic constants C11, C22, and C33 are much greater than the other elastic constants in all calculated β″ phases, resulting in an obvious elastic anisotropy. In order to understand the anisotropic characteristics of these precipitation phases, the Young's modulus anisotropies are evaluated by three-dimensional map as shown in Figure 6. Comparing Figure 6a and c, it can be seen that after Cu atoms substituted Al atoms on the Si3/Al sites, the Young's modulus

Elastic Properties
Here, we compare the mechanical properties of β with or without Cu atoms. The elastic constants of key β phases that are most likely to precipitate during aging were calculated by using fully relaxed crystal structures, and the results are listed in Table 4. According to the Born stability criterion [61], the elastic constants of Mg 4 Al 2 CuSi 4 and Mg 5 AlCuSi 4 all meet the stability criteria of monoclinic crystals. This further supports the stability of Mg 4+x Al 2−x CuSi 4 (x = 0, 1) obtained from the formation enthalpy. The elastic constants C 11 , C 22 , and C 33 are much greater than the other elastic constants in all calculated β phases, resulting in an obvious elastic anisotropy. In order to understand the anisotropic characteristics of these precipitation phases, the Young's modulus anisotropies are evaluated by three-dimensional map as shown in Figure 6. Comparing Figure 6a and c, it can be seen that after Cu atoms substituted Al atoms on the Si3/Al sites, the Young's modulus (E) anisotropy increases significantly; similar results are also shown in Figure 6b,d. This phenomenon indicates that the growth rate of the Cu-containing β phases may be faster than that of the β without Cu. It is consistent with the previous study that Cu can accelerates the age-hardening response [14,28,30]. (E) anisotropy increases significantly; similar results are also shown in Figure 6b,d. This phenomenon indicates that the growth rate of the Cu-containing β″ phases may be faster than that of the β″ without Cu. It is consistent with the previous study that Cu can accelerates the age-hardening response [14,28,30].  Based on the elastic constants in Table 4, the bulk modulus (B), shear modulus (G), and Young's modulus (E) of polycrystalline are calculated by the Hill model [55], and the results are listed in Table 5  Based on the elastic constants in Table 4, the bulk modulus (B), shear modulus (G), and Young's modulus (E) of polycrystalline are calculated by the Hill model [55], and the results are listed in Table 5. Comparing the values of E, G, and B of Mg 4 Al 3 Si 4 and Mg 4 Al 2 CuSi 4 , it can be seen that the values of E, G, and B of β with Cu atoms are higher than that of β without Cu atoms. This relationship is also shown between Mg 5 Al 2 Si 4 and Mg 5 AlCuSi 4 . In general, the Young's modulus (E) can be used to measure the stiffness of the material. The stiffness of the material is greater with the increasing of Young's modulus (E) [64]. It is obvious that the stiffness is enhanced after Cu incorporate into Si3 sites. Pugh [65] proposes using the ratio of the bulk and shear modulus, B/G, to predict brittle or ductile behavior of materials. According to the Pugh criterion, if B/G is more than 1.75, ductile behavior is expected; otherwise, the material would be brittle. From Table 4, the B/G values of calculated β phases are all larger than 1.75, therefore, all the compounds of β phase are ductile with or without Cu atoms and the ductility decreases after Cu atoms incorporate into β . In addition, Poisson's ratio v has been used to measure the shear stability of the lattice, which usually ranges from −1 to 0.5. The smaller the value, the stronger the ability of the crystal to maintain stability during shear deformation [66]. The value of Poisson's ratio v > 0.26 means the ductility of the materials, and the Poisson's ratio of metals is usually 0.25< v < 0.35 [67]. As one can see, all β configurations show ductility with minor differences. It is consistent with the conclusion based on Pugh criterion.

Phonon Spectra
In addition, the dynamic stability is also taken into account. The phonon spectra of Mg 4 Al 2 CuSi 4 and Mg 5 AlCuSi 4 are shown in Figure 7. From Figure 7, one can see that there is no virtual frequency of configuration Mg 4 Al 2 CuSi 4 and Mg 5 AlCuSi 4 , which is generally considered to be dynamically stable.

Phonon Spectra
In addition, the dynamic stability is also taken into account. The phonon spectra of Mg4Al2CuSi4 and Mg5AlCuSi4 are shown in Figure 7. From Figure 7, one can see that there is no virtual frequency of configuration Mg4Al2CuSi4 and Mg5AlCuSi4, which is generally considered to be dynamically stable.

Electronic Structure
The total and partial electronic density of states (TDOSs and PDOSs) for four types of β configurations are calculated to explore the influence mechanism of electronic interaction on structural stability and mechanical properties, as shown in Figure 8, with the Fermi level set to zero. It is evident that incorporating Cu does not change the metallic characteristic of the β phase due to the finite DOS at the Fermi level. At the Fermi level, the TDOS for four types of β configurations at the Fermi level varies. The greatest n (E f ) is 7.41 states/eV/cell in Mg 5 Al 2 Si 4 , followed by 6.40 states/eV/cell in Mg 4 Al 3 Si 4 , 5.27 states/eV/cell in Mg 4 Al 2 CuSi 4 , and 4.18 states/eV/cell in Mg 5 AlCuSi 4 . This indicates that the Cu-containing β phases have a smaller n (E f ). In general, a smaller pseudo gap value n (E f ) corresponds to a more stable structure [68]. This indicates that Mg 4+x Al 2−x CuSi 4 (x = 0, 1) are more stable than the β phases without Cu. The Si-s (range from around 11 eV to 7 eV) and Si-p states (from around 7 eV to the Fermi level) dominate the TDOS of Mg 4 Al 3 Si 4 and Mg 5 Al 2 Si 4 below the Fermi level. In between (ranging from about −7 eV up to −4 eV) regimes, a mixture of s and p character exists, indicating strong hybridization. Especially from −7 eV to −5 eV, the shapes of Si-s and Si-p are very similar, indicating that there is a strong interaction between Si atoms. This may be the origin for the formation of the Si-network; the Si-network acts as a stable skeleton of these phases [32,69]. One can see that Mg-s/Al-s and Si-p in the range from −7 to −4 eV, originating mainly from the s-p hybridization of Si atoms and Mg/Al atoms. The s-states and p-states of Al, Mg, and Si are strongly hybridized above the Fermi level. From Figure 8, it should be noted that, below the Fermi level, the Cu-d state is formed. The s/p orbitals of Mg, Si, and Al all interact with the Cu-d state, and there is obvious electron transfer. The Si-p orbital and the Cu-d orbital are hybridized to form a covalent like bonding, and more electrons are transferred to the new orbital formed by the p-d hybridization. In order to gain a better understanding of the electronic structure of the studied system, the charge density distributions were used as an additional method. The charge-density difference between the (DFT) converged charge density and the isolated atomic charge densities were employed. Figure 9 shows the charge density difference contour plot for the (010) plane to analyze the interaction between Al, Mg, Si, and Cu atoms for the β″ phases. Here we clearly see that there has indeed been a transfer of charge to all the In order to gain a better understanding of the electronic structure of the studied system, the charge density distributions were used as an additional method. The charge-density difference between the (DFT) converged charge density and the isolated atomic charge densities were employed. Figure 9 shows the charge density difference contour plot for the (010) plane to analyze the interaction between Al, Mg, Si, and Cu atoms for the β phases. Here we clearly see that there has indeed been a transfer of charge to all the Si-Si bond regions, it is consistent with the analysis by Derlet et al. [69]. A dominant feature of Figure 9a,b is the concentration of charge between the Si1-S3/Al-Si2-Mg1-Si1 nearest neighbors, and to a lesser extent, between the Si3 and Mg1 nearest neighbors, indicating that covalency plays a role in this system, which was also reported in previous research [70]. Meanwhile, the charge distribution looks like a "charge loop", which can lead to the formation of an "Si network". Strong covalent bonds network can significantly increase the structural stability of β phases. Such a charge transfer to the bonding regions originates from the core regions of both atoms on the Mg and Si sites, in addition to the homogeneous interstitial region between the Mg atoms. The depletion of charge from the Mg3 sites indicates that for this system both metallicity and covalency are present in the bonding. Moreover, the charge transfer density between the Si3 and Si2 sites is slightly decreased, and the charge transfer density between the Mg2 and Si2 sites is increased, indicating the bonds between atoms on Mg2 and Si2 sites are covalent. As shown in Figure 9, the charge ionization of all Mg3 sites is strong, and when the Cu atom on the Si3 site charge ionization becomes stronger, it means both Mg and Cu valence electron are delocalized. The difference is that Mg uniformly provides charges to the surroundings to form a metallic environment [69], while the charges of Cu atoms are delocalized toward Si atoms in unit cells, forming a directional covalent like bond. Si-Si bond regions, it is consistent with the analysis by Derlet et al. [69]. A dominant feature of Figure 9a,b is the concentration of charge between the Si1-S3/Al-Si2-Mg1-Si1 nearest neighbors, and to a lesser extent, between the Si3 and Mg1 nearest neighbors, indicating that covalency plays a role in this system, which was also reported in previous research [70]. Meanwhile, the charge distribution looks like a "charge loop", which can lead to the formation of an "Si network". Strong covalent bonds network can significantly increase the structural stability of β″ phases. Such a charge transfer to the bonding regions originates from the core regions of both atoms on the Mg and Si sites, in addition to the homogeneous interstitial region between the Mg atoms. The depletion of charge from the Mg3 sites indicates that for this system both metallicity and covalency are present in the bonding. Moreover, the charge transfer density between the Si3 and Si2 sites is slightly decreased, and the charge transfer density between the Mg2 and Si2 sites is increased, indicating the bonds between atoms on Mg2 and Si2 sites are covalent. As shown in Figure  9, the charge ionization of all Mg3 sites is strong, and when the Cu atom on the Si3 site charge ionization becomes stronger, it means both Mg and Cu valence electron are delocalized. The difference is that Mg uniformly provides charges to the surroundings to form a metallic environment [69], while the charges of Cu atoms are delocalized toward Si atoms in unit cells, forming a directional covalent like bond. According to the analysis of the thermodynamic results of replacement energies and formation enthalpies in Section 3.1 "Structure stability", the stoichiometry of Cu-containing β″ phases in the precipitation sequence and the sequence of Cu atoms substituting sites in the β″ phases can be inferred. For the stable phases determined from the thermodynamics, the elastic properties of β″ phases with and without Cu were calculated in Section 3.2 "Elastic properties", further supporting the proposed stoichiometry. Besides, the "Phonon spectra" study in Section 3.3 shows that they are also dynamically stable. In summary, the proposed compositions Mg4Al3−xCuxSi4 (1 ≤ x ≤ 3) are reasonable, which is consistent with the results observed in the experiment [36]. In the Section 3.4 "Electronic structure", the origin for the stability of the Cu-containing β″ phases is analyzed from the perspective of electron interaction.

Conclusions
(1) The calculation of the formation enthalpies of 33 Cu-containing β" phases shows that the replacement order of Cu atoms in β" phases can be summarized as one Si3/Al site → two Si3/Al sites → two Si3/Al sites and one Mg1 site. According to the analysis of the thermodynamic results of replacement energies and formation enthalpies in Section 3.1 "Structure stability", the stoichiometry of Cu-containing β phases in the precipitation sequence and the sequence of Cu atoms substituting sites in the β phases can be inferred. For the stable phases determined from the thermodynamics, the elastic properties of β phases with and without Cu were calculated in Section 3.2 "Elastic properties", further supporting the proposed stoichiometry. Besides, the "Phonon spectra" study in Section 3.3 shows that they are also dynamically stable. In summary, the proposed compositions Mg 4 Al 3−x Cu x Si 4 (1 ≤ x ≤ 3) are reasonable, which is consistent with the results observed in the experiment [36]. In the Section 3.4 "Electronic structure", the origin for the stability of the Cu-containing β phases is analyzed from the perspective of electron interaction.

Conclusions
(1) The calculation of the formation enthalpies of 33 Cu-containing β phases shows that the replacement order of Cu atoms in β phases can be summarized as one Si3/Al site → two Si3/Al sites → two Si3/Al sites and one Mg1 site. (2) The Cu atoms strongly favor occupying one of each pair of Si3/Al sites and the most stable Cu-containing β phases were expected to have a stoichiometry of Mg 4+x Al 2−x CuSi 4 (x = 0, 1). In addition, taking into account the change of Cu content in β phases, the stoichiometry of Mg 4 Al 3−x Cu x Si 4 (1 ≤ x ≤ 3) may precipitate.
(3) The calculated mechanical properties show that all calculated β phases are mechanically stable. The incorporation of Cu atoms improves the values of bulk modulus (B), shear modulus (G), and Young's modulus (E) of β , respectively, and all β phases calculated show ductile behavior. Furthermore, the calculation of the phonon spectra shows that Mg 4+x Al 2−x CuSi 4 (x = 0, 1) are dynamically stable. (4) The electronic structure results shows that the Cu atom will join the Si network, and the bond between the Si atom and the Cu atom has the covalent property. The incorporation of Cu atom increases the electron interaction between the Mg2 and the Si3 atom, which may be one of the reasons why the incorporation of Cu atom increases the stability of the β phase structure.