Predicting Atomic Arrangement of Solute Clusters in Dilute Mg Alloys

Abstract

Theoretical prediction of the heterogeneous atomic structures in multicomponent alloys is one of the most challenging issues in materials science. Here we present a first-principles demonstration of this by constructing an on-lattice effective multibody potential model to describe the energetics of hexagonal close-packed Mg crystals containing stacking faults (SFs) with Al and Gd substitutions. Remarkably, we showed that our intuitive model can describe the segregation of solute atoms to SFs and reproduce the characteristic short-range chemical order in dilute Mg–Gd and Mg–Al–Gd systems, including long-period stacking ordered structures, in a manner consistent with recent scanning transmission electron microscopy measurements.

Keywords:

• Stacking Faults
• Solute-Atom Concentration
• Short-Range Order
• Long-Period Stacking Ordered Structures
• Ab Initio-Parameterized Potential

Alloying is important in the production of metallic materials for various scientific applications; it strengthens and enhances metal properties by adding other elements via processing at appropriate compositions and temperatures.[1,2] Due to their light weight, Mg-based alloys containing rare earth (RE) metals,[3] especially Mg–M–RE alloys (M: main group or transition metal) with long-period stacking ordered (LPSO) phases,[4,5] have recently been receiving increased attention because of their outstanding mechanical performances with high yield strength and ductility, which are believed to result from the formation of characteristic microstructures with solute-rich zones; in particular, LPSO phases are known to comprise stacking faults (SFs) at regular intervals where the solute atoms are enriched and even ordered due to their dissolving and binding nature in the matrix. However, because of the lack of experimental observations and theoretical models, details on the fundamental mechanism of solute-atom concentration and clustering at the atomic level have not yet been established.

To properly assess the role of solute-rich zones in the resultant strengthening behavior, identifying and clarifying the controlling factors for the formation and stabilization of solute-atom arrangements in the matrix system are necessary. Predicting the structure at the atomic level, especially in multicomponent systems, remains a fundamental challenge in condensed matter science.[6] State-of-the-art technologies for electron microscopy and imaging—for example, recent advances in high-angle annular dark-field scanning transmission electron microscopy (HAADF-STEM) techniques—have improved our fundamental understanding of the variation in composition and crystal phase containing an ordered arrangement of solute atoms in Mg-based alloys.[5,7–10] However, due to the limited sensitivity and spatial resolution, not all the features of ordering of solute atoms are expected to be readily discernible.

We propose an effective multibody potential model that defines lattice systems of dilute Mg–M–RE ternary alloys while considering the change in stacking of the close-packed planes to provide a theoretical framework for finding energetically favorable configurations from a bottom-up point of view. As an example, the Mg–Al–Gd system was parameterized based on first-principles density functional theory (DFT) calculations in an unambiguous manner, i.e. without any curve fitting. We applied this to atomistic Monte Carlo modeling to predict the segregation and clustering of solute atoms in the system and confirmed its validity by comparison with experimental results.[8,9]

To examine the solid solution properties of Mg–Al–Gd alloy crystals containing SFs, we performed total-energy minimization calculations using the Vienna ab-initio simulation package [11] with the projector-augmented-wave (PAW) construction for pseudopotentials.[12,13] PAW potentials are constructed so that the 3s electrons for Mg, 3s and 3p electrons for Al, and 5p, 5d, and 6s electrons for Gd are treated as valence electrons, and the remaining electrons are kept frozen. The exchange-correlation functional was parameterized in terms of Perdew and Wang's generalized gradient approximation.[14] A plane-wave basis set with 300.5 eV energy cutoff was adopted. We employed (where a₀=3.19 Å and c₀=5.18 Å are the lattice constants of the relaxed hexagonal Mg) orthorhombic supercells of Mg containing 96 atoms with different stacking sequences—hexagonal close-packed (hcp), face-centered cubic (fcc), single (type I₁) and double (type I₂) intrinsic SFs, etc.—in which 1–3 Mg atoms are replaced by Al and/or Gd atoms (see Supplementary Materials section). We applied a Monkhorst–Pack k-mesh [15] to carry out Brillouin-zone integration using the Methfessel–Paxton smearing method [16] (width: 0.2 eV) for ionic relaxations while using the tetrahedron method with Blöchl correction [17] to determine the total energies (E). The E values used as a reference were for relaxed states: the atomic positions and cell parameters are relaxed until their residual forces converge to less than 3 meV/Å.

The binding energies of the solute to the SF and the solute–solute pair (E_b) and the formation energies for the target system (E_f) are defined as follows:

where the E^hcp and E^SF represent the total energies for the hcp and stacking-faulted systems, respectively, X and Y represent two different solute atoms in the Mg supercell, and N is the number of lattice sites considered so that . For the SFs, six kinds of stacking sequence of the Mg layers are involved: hhh, hhc, chc, hch, hcc, and ccc. Here, h and c denote hcp (such as B in ABA) and fcc (such as B in ABC) layers, respectively. Note that a negative value of E_b corresponds to a favorable binding under this definition, and E_f is with respect to the reference energy E^ref of pure Mg, X, and Y in its equilibrium structure (hcp for Mg and Gd and fcc for Al). When a solute–solute pair is considered, two solute atoms are put at the first-, second-, and third-nearest neighbor (1nn, 2nn, and 3nn) distances on the hcp lattice, respectively. Due to the non-ideal c/a ratio of hcp Mg, we distinguished the two non-equivalent 1nn pairs for the solute: one in the same basal (0001) plane (1nn) and the other one out of the basal plane (1nn′). It should be noted that the calculated difference in binding energies less than 0.01 eV is beyond the accuracy of the usual DFT calculations, and thus allowed to use only for the purpose of avoiding the cancellation of significant digits in the numerical parameterization of the interatomic potential as stated below.

The solute binding to the SFs in Mg can be quantitatively characterized by the values listed in Table 1, which correspond to the difference between the dissolution energies of the solute atom in the hcp layer and in the faulted layer. Al and Gd dissolved in fcc-Mg more preferably than in hcp-Mg; in particular, the Gd-binding to the SFs was remarkable. eV at the type I₁ (hċh stacking) SF plane and eV at the type I₂ (hċch stacking) SF plane. Moreover, the binding energies of Al and Gd distinctly depend on the difference in stacking sequence of Al/Gd-substituted layers. It is noteworthy that favorable binding of Al and Gd was clearly observed even at normal hcp layers neighboring SFs, e.g. at the middle layer of ‘hḣc’ stacking, but the strength of the solute binding is weaker at the outer hcp layers than at the inner fcc layers of I₂-type SFs, as shown in Table 1. This indicates that the energetics of the dissolution and binding of Al and Gd in Mg is affected by changes in the stacking sequence over at least five neighboring layers along the [0001] direction.

Table 1.  Binding energies (in eV) between a solute atom (Al or Gd) and SF in the Mg lattice.

Stacking sequence		c^(1)	Al	Gd
hcp	hḣh (BAḂAB)	0	0.000	0.000
	hḣc (BAḂAC)	3	−0.010	−0.048
	cḣc (ACḂCA)	6	−0.027	−0.080
	hċh (BAḂCB)	6	−0.014	−0.050
fcc	hċc (BAḂCA)	9	−0.031	−0.118
	cċc (CAḂCA)	12	−0.056	−0.125

The letters h and c denote the hcp and fcc layers, respectively; the upper dots represent Al- or Gd-substituted layers and the underlines represent SFs. c⁽¹⁾ is the 1nn cubicity number of the substituted atom in Equation (8).

The binding energies between various solute pairs () at the 1nn, 2nn, and 3nn positions were also evaluated in the hcp-Mg matrix (Table 2). For the Al–Al pairs in hcp-Mg, all E_b values at the 1nn, 2nn, and 3nn distances were negligibly small; thus, the Al–Al binding is not a significant factor in the solute concentration. However, interactions between Al–Gd pairs at the 1nn positions are attractive, i.e. eV, where the local strain around the pair can be reduced by combining the atoms of different elements with small and large atomic radii. Gd–Gd binding was found to be repulsive at the 1nn distance but considerably attractive at the 2nn and 3nn ones, i.e. and −0.10 eV, respectively. This indicates that combinations of 1nn Al–Gd pairs and 2nn and/or 3nn Gd–Gd pairs are possible candidate units for the formation of local ordered clusters in Mg–Al–Gd systems.

Table 2.  Binding energies (in eV) for the 1nn, 2nn, and 3nn solute–solute pairs in the hcp-Mg matrix.

1nn	0.021	−0.136	0.022
1nn′	0.018	−0.148	0.066
2nn	−0.024	0.059	−0.172
3nn	−0.014	0.020	−0.104

Note: Negative energy corresponds to an attractive interaction.

To confirm the origin of attractive interactions between the solute atoms, charge density distributions around the 1nn Al–Gd and 2nn/3nn Gd–Gd pairs in the hcp-Mg system were analyzed (Figure 1). The large accumulation of electrons around the Al and/or Gd atoms and their joined shape suggest the presence of chemical interactions that lead to such strong binding energies. It is obvious that the Al–Gd and Gd–Gd pairs in Mg show binding features that may lead to the cross-linked distribution of charge density between the solute atoms. This implies that the binding nature of these solute pairs contributes to steady growth of short-range ordering and clustering among the solute atoms.

Figure 1. Charge density distributions around the (a) 2nn Gd–Gd, (b) 3nn Gd–Gd, and (c) 1nn Al–Gd pairs in the hcp-Mg matrix. White, gray, and yellow spheres represent Mg, Al, and Gd atoms, respectively.

To describe such characteristics concerning the energetics and configurations of the solute atoms, a multibody interatomic potential model for the Mg–M–RE system was constructed that considers the dependence of the strength of pair interaction with the local environment (i.e. bond order effect) in a simplified manner. We describe the potential of the system U as an approximation of the formation energy E_f in the form

where is the coordinates of atom i in the system; and U^2body, U^3body, and U^mbody are two-body, three-body, and multibody interaction energies defined in lattice systems, respectively. The main contribution to the total energy of the system (i.e. bonding and cohesion) comes from the sum of two-body energies in the form

where r_ij is the distance between atoms i and j and α and β are the element types of atoms i and j (, Al, or Gd). Here, is decomposed into its component energy terms:

where r^(l) is the lth nearest neighbor distance in the ideal hcp lattice of Mg; i.e. r⁽¹⁾, r⁽²⁾, and r⁽³⁾ are fixed as 3.2, 4.5, and 5.2 Å,respectively. is the two-body interaction energy between atoms i and j at distance r^(l), and δ is Dirac's delta function.

In the present model, U^3body and U^mbody are introduced as correction terms to the two-body additive approximation for the formation energy. The U^3body term has the following form:

where e^(m) is an additive three-body energy for the 1nn triplets of solute atoms with a bond angle of θ^(m) (see Supplementary Materials section), which represents the energy correction for the bonding state with a rather high concentration of solute atoms. Here, θ_jik is the bond angle between 1nn bonds ij and ik, α is the element type of solute atom i, and β is the element type of solute atoms j and k ( or Gd).

The U^mbody term has the following form:

where and are the multibody energies describing corrections to the formation energy with respect to the local environment of atom i. The indexes α and λ represent the element type and coordination state, respectively, of atom i ( or fcc). As in Equation (8), U^mbody is defined as a function of two kinds of coordination numbers that are introduced to identify the local geometrical and chemical environment of the system. is the number of fcc atoms in the 1nn shell of the atom i, which represents a ‘cubicity’ related to the strength of the fcc-like character in the 1nn local environment. The difference in stacking sequence in the atomic region can be distinguished by using a number from 0 (pure hcp) to 12 (pure fcc) (Table 1). Note that the number depends only on the geometry of the neighboring atoms and not on the element types. Further, represents the number of atoms of the same element in the 2nn shell of the atom i (ranging from 0 to 6), which is used to describe a ‘saturation effect’ in energy gain as the number of 2nn pairs increases. This concept is essential to describing the difference between formation energies of the D0₁₉ and β′ phases of Mg–Gd binary systems [7,8] and to capture the bonding nature of Gd atoms in the Mg matrix (see Supplementary Materials section).

Eventually, all potential parameters described above were obtained from relations concerning the binding and formation energies of the system listed in Tables 1–3. For more details on the potential parameterization, please refer to the Supplementary Materials section.

Table 3.  Reproduction and prediction of formation energies (in eV) from the Mg–M–RE interaction potential developed in this study compared to those from our ab initio DFT calculations.

	This work
System	U^2body	U^3body	U^mbody	Total	DFT (Ef)
Mg96 (hcp: ABABAB)	1.238	0.0	−1.238	0.0	0.0
Mg96 (I1-type SF: ABABAC)	1.238	0.0	−0.825	0.413	0.383
Mg96 (I2-type SF: ABACBC)	1.238	0.0	−0.413	0.825	0.823
Mg96 (fcc: ABCABC)	1.238	0.0	0.0	1.238	1.238
(hcp)	1.250	0.0	−1.182	0.069	0.069
(I1-type SF: ȦBABAC)	1.250	0.0	−0.783	0.467	0.437
(I1-type SF: AḂABAC)	1.250	0.0	−0.779	0.472	0.442
(I2-type SF: ȦBACBC)	1.250	0.0	−0.388	0.863	0.861
(I2-type SF: AḂACBC)	1.250	0.0	−0.384	0.866	0.865
(fcc)	1.250	0.0	0.0	1.250	1.250
(hcp)	1.068	0.0	−1.113	−0.045	−0.045
(I1-type SF: ȦBABAC)	1.068	0.0	−0.751	0.317	0.287
(I1-type SF: AḂABAC)	1.068	0.0	−0.749	0.319	0.289
(I2-type SF: ȦBACBC)	1.068	0.0	−0.406	0.662	0.660
(I2-type SF: AḂACBC)	1.068	0.0	−0.368	0.700	0.699
(fcc)	1.068	0.0	0.0	1.068	1.068
(hcp, 1nn)	1.283	0.0	−1.126	0.157	0.159
(hcp, 1nn′)	1.283	0.0	−1.126	0.157	0.155
(hcp, 2nn)	1.240	0.0	−1.126	0.114	0.114
(hcp, 3nn)	1.250	0.0	−1.126	0.124	0.124
(hcp, 1nn)	0.942	0.0	−0.988	−0.046	−0.068
(hcp, 1nn′)	0.942	0.0	−0.988	−0.046	−0.024
(hcp, 2nn)	0.726	0.0	−0.988	−0.262	−0.262
(hcp, 3nn)	0.795	0.0	−0.988	−0.194	−0.194
(hcp, 1nn)	0.939	0.0	−1.057	−0.118	−0.112
(hcp, 1nn′)	0.939	0.0	−1.057	−0.118	−0.125
(hcp, 2nn)	1.140	0.0	−1.057	0.083	0.083
(hcp, 3nn)	1.101	0.0	−1.057	0.043	0.043
[Figure S2(a)]	0.672	0.013	−0.933	−0.248	−0.312
[Figure S2(b)]	0.627	0.039	−0.933	−0.267	−0.269
[Figure S2(c)]	0.627	−0.085	−0.933	−0.390	−0.390
[Figure S2(d)]	0.672	0.013	−0.933	−0.248	−0.184
[Figure S2(e)]	0.455	0.057	−0.933	−0.420	−0.420
[Figure S2(f)]	0.627	0.039	−0.933	−0.267	−0.265
[Figure S2(g)]	0.524	−0.026	−0.933	−0.435	−0.435
[Figure S2(a)]	0.830	−0.012	−1.002	−0.184	−0.176
[Figure S2(b)]	0.810	0.021	−1.002	−0.171	−0.155
[Figure S2(c)]	0.810	0.037	−1.002	−0.155	−0.155
[Figure S2(d)]	0.830	−0.012	−1.002	−0.184	−0.192
[Figure S2(e)]	0.786	−0.002	−1.002	−0.217	−0.217
[Figure S2(f)]	0.810	0.021	−1.002	−0.171	−0.187
[Figure S2(g)]	0.796	−0.004	−1.002	−0.209	−0.209
(D019)	−0.738	0.0	0.403	−0.335	−0.335
(base-centered orthorhombic)	−0.514	0.0	0.139	−0.375	−0.375
(monoclinic, 18R-type LPSO)^a	−2.449	0.030	−0.155	−2.574^b	−2.541
(monoclinic, 18R-type LPSO)^c	−3.085	0.030	0.272	−2.783^b	−2.663

U^2body, U^3body, and U^mbody represent the two-body, three-body, and multibody interaction energies, respectively, broken down from the total formation energy.

^aTwo formula units per unit cell containing one cluster.

^bNot fitted but predicted.

^cFour formula units per unit cell containing two clusters.[9]

To prove the validity and effectiveness of the developed potential, thermal equilibrations of the Mg–Gd and Mg–Al–Gd systems with and without SFs were carried out by using Markov chain Monte Carlo (MCMC) simulations.[18] To accelerate the sampling process for such complex multicomponent systems, we employed the MCMC method without detailed balance.[19] For the equilibrium simulations of the Mg–Gd and Mg–Al–Gd supersaturated solid solutions, orthorhombic supercells of Mg containing 1344 atoms with 2H stacking were used. As an initial configuration, solute atoms were randomly distributed on the sites of a hcp lattice with a composition of Mg–5.4at%Gd (), Mg–12.5at%Gd (), Mg–4.0at%Al–5.4at%Gd (), and Mg–9.4at%{ }Al–12.5at%Gd () with an Al/Gd ratio of ¾. The simulations starting from multiple (typically five) initial configurations with different sets of random positions of solute atoms were conducted for 5×10⁶ steps at 500 K via successive swapping of neighboring pair elements under periodic boundary conditions.

For the Mg–Gd system (Figure 2), two kinds of characteristic local clusters were observed: linear zigzag patterns (along the direction) and localized hexagonal patterns. These corresponded to local structures of β′-Mg₇Gd (base-centered orthorhombic) [7] and (D0₁₉), respectively. Both are arranged from 2nn and 3nn Gd–Gd pairs, and the coexistence of the two kinds of short-range orders (corresponding to β″ phase) was experimentally observed by HAADF-STEM images of an Mg–5at%Gd alloy at early stages of aging at 200°C.[8] The linear zigzag patterns became dominant over the hexagonal ones as the aging time increased as a reflection of the lower formation energy of the β′ structure than that of the D0₁₉-type structure (see Supplementary Materials section). For the Mg–Al–Gd system (Figure 3), the linear zigzag patterns disappeared, and only rather isolated hexagonal patterns survived. Note that the Al enrichment was distinctly confined inside the hexagonal arrangements formed by the 2nn and 3nn Gd–Gd pairs. This indicates that the addition of Al stabilizes the local D0₁₉-type structure because the energy gain by the formation of 1nn Al–Gd pairs overcomes the energy loss caused by the change in Gd arrangements with high .

Figure 2. Solute clusters distribution for the (a) Mg–5.4at%Gd and (b) Mg–12.5at%Gd alloy systems with 2H stacking projected on the (0001) plane. Only Gd atoms are displayed in yellow spheres. (A) Initial configuration with all atoms distributed randomly over the hcp lattice. (B) Clustering because of a 5×10⁶-step MCMC run at 500 K. Bars between the solute atoms represent the 2nn Gd–Gd pairs.

Figure 3. Solute clusters distribution for the (a) Mg–4.0at%Al–5.4at%Gd and (b) Mg–9.4at%Al–12.5at%Gd alloy systems with 2H stacking projected on the (0001) plane. Only solute atoms are displayed: Gd atoms in yellow and Al atoms in gray spheres. (A) Initial configuration with all atoms distributed randomly over the hcp lattice. (B) Clustering because of a 5×10⁶-step MCMC run at 500 K. Bars between the solute atoms represent the 1nn Al–Gd pairs.

We performed MCMC simulation of the Mg–Al–Gd system with an SF by employing a monoclinic supercell of Mg containing 144 atoms with 18R stacking, where an I₂ SF was introduced in the middle of six atomic layers (i.e. hhcchh). At the beginning of the simulation, solute atoms were randomly distributed on the sites of a lattice with a composition of Mg–8.3at%Al–11.1at%Gd (). Figure 4 and the supplemental movie (SupplMovie.mpg) show that solute atoms gathered near each other due to the formation of 1nn Al–Gd and 2nn Gd–Gd pairs and also gradually segregated into the SF layer during the configuration sampling. The formation and in-plane alignment of L1₂-type solute clusters (Al₆Gd₈×2) were clearly achieved after equilibration of the system, which is completely consistent with recent HAADF-STEM measurements for the Mg–Al–Gd LPSO phases.[9] The obtained L1₂-ordered configuration was not a fitting target in the development of our potential model but was predicted.

Figure 4. Evolution of potential energy for the Mg–Al–Gd alloy system with 18R stacking at 600 K. The inset shows snapshots of the simulated solute atom clusters during MCMC sampling.

Table 3 lists the formation energies of the various Mg–Al–Gd configurations obtained from the potential developed here. We found that this interaction potential provides reasonable agreement with the results from the DFT calculations. The U^mbody term (i.e. f⁽¹⁾ and f⁽²⁾) proposed here played a crucial role in reproducing the energetics of the various Mg–Al–Gd configurations. This suggests that this term pertinently extracts the essential features of the inherent multibody nature of the dilute Mg–M–RE system and that the complex and cooperative interactions in the system can be characterized and understood from the viewpoint of the cubicity and chemical occupancy in the local environment with a range of 2r⁽¹⁾=6.4 Å at most. In this study, the effects of long-range interactions (>6.4 Å) on the system were omitted and thus the energetic preference between the different stacking polytypes in the LPSO phases [9] was not distinguished. The analysis about the effect of such interactions on the formation of long-period stacking and chemical order will be the subject of the future work.

We confirmed that our approach can successfully describe the formation and arrangement of Al/Gd clusters in a system and provide clear-cut insights into the short-range chemical order in Mg-based alloys, including LPSO structures, by taking into account the identified multibody effects on the system. The considerable agreement between the ab initio potential-based predictions and HAADF-STEM images is encouraging with regard to predictive validity as a physics-based framework for alloy exploration and explanation; this will help improve our understanding of these multicomponent systems and provide further clarity on the organizing principles that govern local atomistic ordering in various alloy systems. We firmly believe that such efforts will open a new design space in the advanced alloy engineering.

Supplementary online material

A more detailed information on calculations is available at http://dx.doi.org/10.1080/21663831.2013.838705.

Acknowledgements

This study was supported by Grants-in-Aid for Scientific Research on Innovative Area, No. 23109004; Young Scientists (A), No. 24686072; Scientific Research (A), No. 23246025; and Challenging Exploratory Research, No. 25630013. We also appreciate the support from the Elements Strategy Initiative for Structural Materials, No. 22656030, and useful discussions with M. Fronzi.