Emergent B2 chemical orderings in the AlTiVNb and AlTiCrMo refractory high-entropy superalloys studied via first-principles theory and atomistic modelling

In a substitutional alloy with a fixed underlying crystal lattice, with the set of lattice positions denoted by , an arrangement of atoms (a ‘configuration’) can be uniquely specified by a set of site occupancies, , where if site i is occupied by an atom of chemical species α, and otherwise. Each lattice site is constrained to be occupied by one (and only one) atom, which leads to the condition

It is then natural to consider the average value of these site occupancies, which leads to the definition of the site-wise concentrations,

where denotes an average taken with respect to the relevant thermodynamic ensemble. Note that, by construction, we have that . These site-wise concentrations represent long-range order parameters classifying potential chemically ordered phases.

In the limit of high-temperature, where the free energy landscape is dominated by entropic contributions and the alloy is maximally disordered, the site-wise concentrations become spatially homogeneous, i.e. any atom can occupy any lattice site with probability equal to its average concentration in the alloy. This is equivalent to the statement that

where c_α is the overall (total) concentration of species α.

The ensemble average of the internal energy of this disordered phase can be described ab initio via application of the CPA [54–57] within the Korringa–Kohn–Rostoker (KKR) formulation [58–60] of DFT [61, 62]. The KKR method uses multiple scattering theory (MST) to construct the single-particle Green’s function for the Kohn–Sham equations, while the CPA seeks to construct an effective medium of electronic scatterers whose overall scattering properties approximate those of the disordered solid solution [54]. We emphasise that the CPA assumes a fixed underlying crystal lattice, with no direct incorporation of local lattice distortions, as are known to be present in many HEAs [63]. Despite this approximation, we note here that the CPA has been established as a powerful and efficient technique for describing the electronic structure and physical properties of HEAs. It has previously been shown to accurately capture details of their electronic structure [64], the nature of their heavily-smeared Fermi surfaces [12], aspects of their magnetism [65, 66], and also a range of transport [67–69] and structural [70–72] properties.

2.2.1. Describing chemical fluctuations: concentration waves

The KKR-CPA is not limited to studies of systems where all lattice sites have the same set of partial occupancies. It can also be applied to structures with multiple sublattices, each with different associated partial lattice site occupancies. For a particular (partially) ordered structure and a given set of (partial) sublattice occupancies, it is therefore possible to compute the total DFT energy of the system, . However, when assessing which chemically ordered structures are most energetically favourable, direct evaluation of the energy of all potential chemically ordered (and segregated) structures is both laborious and computationally expensive, requiring evaluation of a large number of (partially) ordered supercells. An alternative approach, which we use in this work is to consider the change of the internal energy of the system induced by an infinitesimal chemical perturbation applied to the disordered solid solution [50–53].

These perturbations can be expressed concisely using the language of concentration waves. We express an inhomogeneous set of site-wise concentrations (equation (2)) as a sum of a homogeneous term—the average concentrations of equation (3)—plus an inhomogeneous perturbation, writing

where denotes a spatially inhomogeneous chemical fluctuation. Given the translational symmetry of the underlying crystal lattice, it is convenient to express these chemical fluctuations in reciprocal space,

where denotes a static concentration wave with wavevector k and chemical polarisation Δc_α. (We note that, for most crystallographically ordered structures, the sum over k typically runs over a handful of high-symmetry k-vectors.)

This multicomponent concentration wave formalism provides a concise and elegant way of describing a range of chemically ordered structures, and has its roots in pioneering work on concentration waves in binary alloys by Khachaturyan [73] and Győrffy and Stocks [74]. Some one-dimensional examples of concentration waves and associated chemical polarisations are provided in figure 2. In three dimensions, when the underlying lattice is bcc, the wavevector(s) associated with a B2 chemical ordering shown in figure 1 is and equivalent, while the associated chemical polarisation, then describes which chemical species are moving to which sublattice.

2.2.2. Energetic cost of chemical fluctuations: the theory for multicomponent alloys

To quantify the energetic cost of such inhomogeneous chemical fluctuations, and the temperatures at which they may emerge, we begin by writing down an approximate expression for the grand potential, or Landau free energy, Ω. In general, this takes the form

where U is the internal energy, T the temperature, S the entropy, µ the chemical potential(s), and N the particle number(s). For the description of the alloy considered in this paper, the free energy is approximated via

In the above expression, the first term represents the average internal energy as obtained within the CPA, the second term represents the so-called entropy of mixing, and the third term represents the contribution to the free energy made by the sitewise chemical potentials, ν_iα, which, in principle, can vary for each chemical species and lattice site. Proceeding, we expand this free energy about our reference state—the homogeneous solid solution—in terms of the applied inhomogeneous fluctuation . This constitutes a Landau-type series expansion of the free energy, and is written

The site-wise chemical potentials, ν_iα, of equation (7) act as Lagrange multipliers in the linear response theory, serving to conserve the overall concentration of each chemical species. However, the variation of this term of the free energy is understood to be irrelevant to the underlying physics, so its derivatives are dropped in the analysis which follows [50–53]. Taken in conjunction with the translational symmetry of the solid solution, the requirement that overall concentrations of each chemical species be conserved by chemical fluctuations leads to the result that the first-order term of equation (8) is exactly zero. To second-order, the change in free energy induced by the applied chemical perturbation is written

The first term in the square brackets, , matrix which is both diagonal and positive definite, and is associated with entropic contributions to the free energy. The second term, is the second concentration derivative of the average energy of the disordered solid solution as described by the KKR-CPA.

To evaluate the quantity , a ring of coupled equations involving various CPA-related quantities must be solved self-consistently, carefully incorporating the rearrangement of the electrons due to the applied chemical perturbation. The relevant set of coupled equations have been defined in [50], and their solutions were first examined and discussed in [51]. It should be emphasised that the present scheme is a rigorous generalisation of an earlier theory used to examine the phase stability of binary alloys in a similar manner [75–77].

In the perturbative analysis described in [50], is evaluated in reciprocal space. The change in free energy of equation (9), therefore, is re-expressed as:

We refer to the matrix in square brackets, , as the ‘chemical stability matrix’. It can be shown to be directly related to an estimate of the atomic short-range order (ASRO) [50]. Moreover, it can be considered as an analogy of a Hessian matrix of second derivatives evaluated at a stationary point of an alloy’s free energy landscape. To infer a disorder–order transition, we start in the limit of high temperature, where entropy dominates the free energy landscape and all eigenvalues of this matrix are positive. (In this limit, the system is stable to applied chemical perturbations.) Proceeding, we progressively lower the considered temperature, searching, for a temperature at which the lowest lying eigenvalue of this matrix passes through zero for some k-vector in the irreducible Brillouin zone. At this temperature, denoted T_ord, the solid solution is unstable to a chemical perturbation described by a concentration wave with wavevector , and we infer the presence of a chemical disorder–order transition. The eigenvector associated with the eigenvalue passing through zero is the chemical polarisation of the concentration wave to which the solid solution is unstable, Δc_α. It describes which chemical species are partitioning themselves onto which sublattice(s).

2.2.3. Pairwise form of the alloy internal energy: the Bragg–Williams Hamiltonian

The above concentration wave analysis and associated linear response theory provides information about the initial inferred chemical ordering in the system, but it is also possible to go further and map the derivatives of the internal energy of the disordered alloy, , to a real-space effective pair interaction (EPI) describing the internal energy of a configuration with discrete lattice site occupancies. When appropriate sampling technique(s) are applied to this model, the phase behaviour of the alloy can be studied in detail. Moreover, equilibrated, lattice-based configurations can subsequently be extracted for illustration and further study.

The real space model and associated atom–atom EPIs are for the Bragg–Williams Hamiltonian [78, 79], which can be thought of as a multicomponent Lenz–Ising Hamiltonian [80], taking the form

where denotes the EPI between an atom of chemical species α on lattice site i and an atom of chemical species on lattice site j. Assuming interactions are homogeneous and isotropic, we can express this as a sum of interactions over first-nearest neighbours, second-nearest neighbours, etc, writing to denote the interaction between species α and on coordination shell n. The effective pairwise interactions, are recovered from by means of a inverse Fourier transform. The mapping from reciprocal-space to real-space and fixing of the gauge degree of freedom on the is specified in earlier works [50–53].

Using the pairwise Hamiltonian of equation (11), with atom–atom EPIs recovered from the ab initio data, we can employ lattice-based Monte Carlo simulations to further explore the alloy configuration space. Such simulations facilitate validation of the concentration wave analysis outlined above, and also allow us to search for further phase transitions occurring below any initial chemical ordering temperature. In this work, these atomistic Monte Carlo simulations consist of two different Markov chain random walks, both making use of Metropolis–Kawasaki dynamics [81], to explore the configuration space of the AlTiVNb and AlTiCrMo RSAs. Kawasaki dynamics ensure conservation of the overall concentration of each chemical species in the alloy by only performing swaps of pairs of atoms in the simulation cell. The sampling methods employed are Metropolis–Hastings Monte Carlo [82] and Wang–Landau sampling [83]. We outline the details of these two sampling algorithms below.

2.3.1. Metropolis–Hastings algorithm

The Metropolis–Hastings algorithm allows for a system to follow a chain of states to an equilibrium ensemble in a indeterminate but finite time [82, 84]. Phase equilibria are obtained by performing swap moves with a probability which depends on the energy difference between the initial and subsequent states and the simulation temperature.

For relaxation models such as the ones in this paper the time-dependent behaviour obeys

where is the probability of the system being in a state n at a time t, m is the final state and is the transition rate . When the system is in equilibrium, , and we obtain the expression

which is known as detailed balance. If the dynamics satisfy this equation then the simulation is able to reach equilibrium. To calculate the acceptance rate (attempt rate multiplied by acceptance rate), the only quantity needed is the energy difference between the initial and proposed state, , resulting in the Metropolis acceptance probability

where T denotes the simulation temperature, and is the usual Boltzmann constant. In practice, for the alloy simulations conducted in this work, the Metropolis–Hastings algorithm allows us to obtain sample atomic configurations from an equilibrated ensemble at a range of temperatures which are suitable for visualisation and further study. At a given temperature, we initialise a supercell containing the correct overall concentration of each chemical species, where the initial atomic site occupancies are randomly generated. Two lattice sites (not necessarily nearest neighbours) are then selected at random and the energy difference obtained by swapping their atomic occupancies is computed. The swap is accepted according to the transition probability of equation (14). This process is repeated until the internal energy of the simulation has reached a stable equilibrium by monitoring how the energy fluctuates about a stable energy average across a set number of previous states. Phase equilibria obtained using this method, allows for visualisation of alloy ordering and contextualise the results of Wang–Landau sampling method, the details of which follow.

2.3.2. Wang–Landau sampling

The Wang–Landau sampling method [83, 84] is a flat energy histogram method that has a wide applicability due to its ability to obtain the density of states (DOS) in energy from which information on thermodynamic quantities at any temperature can be obtained. In this context, the simulation DOS can be defined via the classical partition function. The conventional definition of the partition function with discretely labelled configurations i is rewritten as

where E is the energy, is the Boltzmann factor, T is the temperature and g(E) is the DOS. For the case of E_j, the energy represents a discretised energy bin within the energy histogram which is treated as a single energy macrostate. Wang–Landau sampling obtains an estimate for g(E) across a chosen energy range during a Monte Carlo simulation. The DOS can begin with a simple estimate such as and is improved throughout the course of the simulation. Atoms are swapped using the same method as in the Metropolis–Hastings case according to the probability

where E₁ is the initial energy and E₂ is the energy of the proposed swap. After each proposed swap, the DOS is updated according to

where E is the energy of the resultant state, f_k is a modification factor initially (k = 0) greater than 1, and k is the current iteration index of the Wang–Landau sampling algorithm. A histogram of the energies visited is maintained, , as is a measure of the ‘flatness’, F of the histogram,

Once F is below a given tolerance, sampling is interrupted and f is reduced for the next sampling iteration, e.g. . Then the histogram entries are set to zero and the sampling continues until the flatness tolerance is achieved again. This process continues until the modification factor, f, is sufficiently close to unity, e.g. , and it is decided that g(E) is known to an acceptable tolerance.

Once the DOS for the simulation has been obtained, the energy probability distribution at a given temperature can be found by using

from which we can extract a variety of system properties. One such property is the specific heat capacity (SHC), C, as a function of temperature, T, recovered via [84]

where E is the energy of a simulation at a particular temperature, and the relevant ensemble averages are taken using the energy probability distribution from equation (19).

2.3.3. Warren–Cowley ASRO parameters

To quantify the temperature-dependent ASRO in our Monte Carlo simulations, we use the Warren–Cowley ASRO parameters [85–87] adapted to the multicomponent setting, defined as

where n refers to the nth coordination shell, is the conditional probability of an atom of type q neighbouring an atom of type p on shell n, and c_q is the total concentration of atom type q. When , p–q pairs are disfavoured on shell n and, when they are favoured. The value corresponds to the ideal, maximally disordered solid solution.

These ASRO parameters are evaluated across the sampled energy range (in practice, by evaluating the average ASRO for configurations drawn from each ‘bin’ of the energy histogram) and then reconfigured to be written as a function of temperature

where is the energy probability distribution at a given temperature.

In combination, the SHC and Warren–Cowley ASRO parameters facilitate a detailed description of emergent phase transitions in a multicomponent alloy. The SHC data allows us to accurately identify the temperature at which a phase transition occurs, while the ASRO parameters help to quantify the nature of the transition in terms of atom–atom correlations.

A fundamental quantity of solid state physics is the electrical resistivity of a material, ρ. At its most basic level, this quantity allows one to distinguish between metals, semiconductors, and insulators. In metals and alloys, the electrical resistivity is directly connected to the electronic mean free path of states at the Fermi level, . In a crystalline, metallic material in the limit K, Bloch states are eigenvalues of the electronic Hamiltonian and the electronic mean free path diverges, , resulting in a vanishing resistivity, . However, in a disordered, substitutional alloy, translational symmetry is violated, Bloch’s theorem does not apply, and the electronic mean free path is finite even at T = 0 K. This results in a finite residual resistivity, ρ₀, defined as the limiting value of the resistivity as , which can provide insight into the effects of chemical disorder on the electronic structure of a material [68].

The KKR-CPA naturally captures the broken translational symmetry of a substitutionally disordered alloy, and provides a means by which to evaluate the Green’s function, G, for such a system [60]. From this Green’s function, one can then use the linear response Kubo–Greenwood formula [88, 89] applied to the KKR-CPA effective medium [90] to evaluate the resistivity at both zero and finite temperatures, where the latter can also incorporate the effects of thermally induced ionic displacements within an ‘alloy analogy’ model [60]. The Kubo–Greenwood formula for the conductivity tensor, σ_µν, implemented in SPR-KKR [60, 91] is given by

where V is the simulation cell volume, is the retarded Green’s function at the Fermi level, and j_µ the current density operator with Cartesian coordinate indices µ and ν. Angled brackets, , indicate the CPA average over substitutional disorder. From the conductivity tensor, it is a simple matter to recover the resistivity of a material. In this study, we employ this approach to evaluate the influence of (partial) chemical order on the electronic transport properties of the examined alloys.

We begin our analysis by performing self-consistent DFT calculations to model the electronic structure of the disordered solid solutions. We use the all-electron SPR-KKR package [60, 91] to construct the self-consistent potentials of the KKR [58–60] formulation of DFT [61, 62], using the CPA [54–56] to average over chemical disorder. We perform scalar-relativistic calculations within the atomic sphere approximation (ASA) [92], employing an angular momentum cutoff of for basis set expansions, and a 32-point semi-circular contour in the complex plane to integrate over valence energies. Sampling for integrations over the Brillouin zone during the self-consistency cycle used a parameter of NKTAB = 5000, resulting in a dense k-point mesh over the first Brillouin zone. We use the local density approximation (LDA) and the exchange-correlation (XC) functional of Vosko, Wilk, and Nusair [93]. Both systems are treated as non-magnetic; for AlTiCrMo a self-consistent, spin-polarised, disordered local moment calculation [94] (to represent the paramagnetic state) was tested but the local moments collapsed during the self-consistency cycle, which is consistent with calculations for elemental Cr [95]. We therefore believe that non-magnetic calculations are representative of the electronic structure of these alloys at typical processing temperatures. Further details of the self-consistent calculations and the relevant data can be found in the open-access dataset associated with this work [96].

We optimise the lattice parameter and associated unit cell volume for the disordered bcc (A2) phase of both alloys by calculating the total DFT energy across a range of lattice parameters and fitting the results to a stabilised jellium equation of state [97] as implemented in the Atomic Simulation Environment [98]. This procedure also allows us to estimate the bulk modulus for both alloys. Our optimised lattice parameters and calculated bulk moduli are shown in table 1. Our optimised bcc cubic lattice parameters for the AlTiVNb and AlTiCrMo alloys of 3.113 Å and 3.025 Å, respectively, compare reasonably with the experimentally determined values of 3.186 Å [24] and 3.100 Å [26], though are slight underestimates. This underestimation is typically expected when using LDA exchange correlation functionals on materials containing transition metals while, conversely, generalised gradient approximation (GGA) functionals typically slightly overestimate lattice constants [99]. In addition, our DFT calculations do not account for the modest thermal expansion of the lattice which would be present at room temperature.

Proceeding, in figure 3, we plot the electronic DOS and Bloch spectral function (BSF) along high-symmetry lines of the Brillouin zone of the bcc lattice for both AlTiVNb and AlTiCrMo simulated in the chemically disordered A2 phase. The DOS and BSF are plotted for an energy range around . The BSF can be thought of as a k-resolved DOS [60]. For a pristine crystal, it reduces simply to the conventional bandstructure. However, in a system with broken translational symmetry, such as the substitutional alloys of this work, pristine bands are smeared out by the disorder. The degree of smearing of bands can be related to the average lifetime (or mean free path) of electronic states in a material, which leads naturally to calculations of quantities such as the residual resistivity of a material.

Zoom In Zoom Out Reset image size

Considering first the BSF data for the two alloys, it can be seen that the broken translational symmetry, associated with the substitutional disorder, heavily smears electronic states across the entire considered energy range. It can also be seen that the narrow d-bands associated with the transition metals are shifted down relative to the Fermi level, , for AlTiCrMo as compared to AlTiVNb, which is associated with the increased valence of Cr and Mo as compared to V and Nb.

When considering the species-resolved DOS, contrasts should be made between the sp states associated with Al (a ‘simple’ metal); those associated with Ti, V, and Cr (3d metals); and those associated with Nb and Mo (4d metals). In elemental Al, electronic states are nearly free-electron like in character, and the DOS can be expected to be approximately proportional to . For both of the alloys considered here, that the species-resolved DOS for Al diverges from this behaviour and displays localised peaks around the d-states of the transition metals is indicative of the formation of hybridised p–d bonding states between Al and the other elements present in the compositions. It should also be noted that the width of the d-bands associated with Ti, V, and Cr (3d transition metals) are narrower than those associated with Nb and Mo (4d transition metals). In previous studies, we have found that both bandwidth differences [52] and p–d hybridisation [100] can drive strong atomic ordering tendencies in multicomponent alloys such as the ones considered in this work.

To study the nature of ASRO in the solid solution, as well as to infer the temperature at which the solid solution becomes unstable to chemical fluctuations and a long-range crystallographically ordered structure emerges, we look to the S(2) theory for multicomponent alloys, as outlined in section 2.2. We use a computational implementation of this theory of which the details have been discussed extensively in earlier works [50–53]. This methodology has previously been applied with success to the Cantor alloy and its derivatives [51, 101], the refractory HEAs without Al present in the composition [52, 102], and to the Al_xCrFeCoNi system [100].

Shown in the top row of figure 4 are computed values of along high-symmetry lines of the irreducible Brillouin zone of the bcc lattice, while the bottom row shows the eigenvalues of the chemical stability matrix constructed from the ab initio data for both alloys along the same high-symmetry directions and evaluated at a temperature of T = 2600 K. The quantity can be thought of as the Fourier transform of an EPI between different chemical species in the alloy, indicating the strength and nature of various atom–atom correlations in the solid solution. Considering the data for (a) AlTiVNb and (b) AlTiCrMo, we see that there are strong interactions in both alloys for Al–Al, Al–Ti, Al-4d, Ti–Ti, and 4d–4d elemental pairs. (Note that Nb and Mo are the 4d elements present in the two compositions.) However, despite their presence in both compositions, the Al–Al, Al–Ti, and Ti–Ti data, although qualitatively similar, are substantially different in numerical value between the two alloys. This emphasises that extrapolation of atom–atom interactions from binary subsystems [103, 104] may not be a reliable approach for modelling the thermodynamics of HEAs.

Zoom In Zoom Out Reset image size

Proceeding, we now consider the eigenvalues of the chemical stability matrix, for (c) AlTiVNb and (d) AlTiCrMo evaluated at a temperature of T = 2600 K, where both matrices are strictly positive definite, i.e. all eigenvalues are greater than zero. For both alloys, the eigenvalue with lowest energy occurs at the H point, corresponding to and associated with a B2 chemical ordering as illustrated in figure 1. That, at a temperature of T = 2600 K, the lowest-lying eigenvalue for AlTiVNb is lower than that for AlTiCrMo indicates that the B2 chemical ordering temperature for AlTiVNb will be higher than that for AlTiCrMo. The local minimum eigenvalue at the P point (associated with B32 ordering tendencies [52]) for AlTiCrMo is perhaps suggestive of some weaker secondary ordering tendencies in this alloy.

As outlined in section 2.2, we then search for the temperature at which an eigenvalue passes through zero for some wavevector in the irreducible Brillouin zone, as this provides a mean-field estimate of the chemical ordering temperature. Table 2 gives the predicted ordering temperature, wavevector, and associated eigenvector describing the chemical orderings for both AlTiVNb and AlTiCrMo RSAs. For both alloys, the wavevector describing the ordering is , associated with a B2 chemical ordering. For AlTiVNb this ordering is predicted to occur at 2574 K, while for AlTiCrMo it is predicted to occur at 2107 K. The eigenvectors or ‘chemical polarisations’ associated with these orderings then provide insight as to which chemical species are preferentially sitting on which sublattice. For AlTiVNb, Al and Nb are predicted to move to one sublattice, while Ti and V move to the other, with Al and Ti having stronger site preferences than V and Nb. For AlTiCrMo by far the strongest site preference is for Ti, with Al and Mo moving to the other sublattice, and Cr remaining comparatively disordered, indicated by the small value of .

Table 2. Predicted chemical ordering temperatures (), associated concentration wavevectors (), and chemical polarisations (Δc_α) for the two RSAs considered in this work. For both alloys, the wavevector describing the chemical ordering is , which describes a B2 chemical ordering imposed on the bcc lattice. For AlTiVNb, the chemical polarisation of the concentration wave suggest that the B2 ordering will have one sublattice rich in Al and Nb, while the other will be rich in Ti and V, with Al and Ti having the strongest site preferences. For AlTiCrMo, Al and Mo are anticipated to move to one sublattice, Ti the other, with Cr remaining relatively disordered and spread across both sublattices.

Alloy	(K)	()	Structure
AlTiVNb	2574	(0,0,1)	B2	0.700	−0.571	−0.361	0.232
AlTiCrMo	2107	(0,0,1)	B2	0.521	−0.792	−0.044	0.315

The temperature and nature of these B2 chemical orderings are in agreement with existing literature, both experimental and computational. (We note that our computed transition temperatures in table 2 are likely to be overestimates as they are computed within a single-site or ‘mean-field’ theory.) Experimentally, [25] reports that AlTiVNb samples annealed for 24 h at 1200 ^∘C were composed entirely of a B2 phase, while samples annealed at 1000 ^∘C and 800 ^∘C contained a majority B2 phase with some Nb₂Al-like σ-phase precipitates. Computationally, Körmann et al [37] simulated the phase stability of AlTiVNb across a wide temperature range using DFT calculations, a machine-learned interatomic potential, and Monte Carlo simulations. They report a B2 chemical ordering occurring at approximately 1700 K, with Al and Ti expressing strong site preferences, and Nb and V expressing weaker preferences, in good agreement with our own predictions. For AlTiCrMo, [26] reports combined experimental results and thermodynamic calculations. The thermodynamic calculations estimate a B2 ordering temperature of just over 1100 ^∘C (≈ 1500 K), in alignment with our result, though these authors also report that it is challenging to identify the B2 ordered structure in their experimental samples via x-ray diffraction due to counterbalancing of sublattice occupancies and atomic form factors in their predicted B2 ordered structure.

Though the eigenvectors of table 2 describe initial (infinitesimal) change in partial lattice site occupancies, they can be used to infer partially ordered phases by allowing the size of the fluctuation to ‘grow’ until one sublattice contains (at least) one atomic species whose concentration reaches zero. This condition identifies the largest permitted chemical fluctuation consistent with that concentration wave’s chemical polarisation. We define an order parameter η to quantify the degree of (partial) ordering, where η = 0 corresponds to the disordered solid solution, and η = 1 corresponds to the largest permitted chemical fluctuation consistent with the predicted chemical polarisation. For both alloys, the predicted partial occupancies of the two non-equivalent lattice sites as a function of η are provided in the supplementary material [105].

To assess the impact of these chemical orderings on the mechanical properties of the alloys, we take the obtained partial lattice site occupancies for the case η = 1 and compute revised lattice parameters and bulk moduli. These self-consistent calculations were again performed using SPR-KKR [60, 91], with the same computational parameters as for calculations on the disordered phases, the results of which are provided in table 3. When these data are compared to the calculations performed on the disordered alloys (table 1), it can be seen that the predicted chemical ordering for AlTiVNb results in a small lattice contraction and consequent increase in bulk modulus, while for AlTiCrMo the ordering results in a small lattice expansion and consequent decrease in the bulk modulus. However, these changes are minor, and we conclude that the predicted chemical orderings do not significantly impact the alloys’ mechanical properties.

However, as a result of the B2 chemical orderings, the DFT internal energy of the alloys are found to be lowered by 49.6 meV atom⁻¹ and 37.2 meV atom⁻¹ for AlTiVNb and AlTiCrMo respectively. When we consider the electronic DOS for both alloys in their predicted B2 chemically ordered structures, shown in figure 5, we see how this chemical ordering impacts the alloys’ electronic structure. For both systems, there is a modest reduction in the DOS at the Fermi level, , and the peak in the DOS below is also shifted to lower energies. Some sharper features at low energies are also seen to emerge for both alloys, suggestive of the formation of more localised bonding states. Finally, for AlTiVNb, a new peak in the DOS above can be seen to emerge.

Zoom In Zoom Out Reset image size

From the ab initio reciprocal space data, we Fourier transform and fit atom–atom EPIs for both of the considered RSAs. These EPIs are for the Bragg–Williams Hamiltonian, equation (11). We assume that the interactions are homogeneous (translationally invariant) and isotropic, and write to denote the EPI between chemical species α and chemical species at nth nearest neighbour distance.

We plot as a function of coordination number, n, for a fit to the first ten coordination shells of the bcc lattice for AlTiVNb and AlTiCrMo in figure 6. For both alloys, it can be seen that interactions are dominated by the first few coordination shells and tail off quickly with increasing distance. It should be noted that, despite the fact that Al and Ti are present in both of the considered compositions, Al–Al, Al–Ti, and Ti–Ti interactions are different between the two systems. This confirms our earlier assertion that extrapolating atom–atom interactions from data for binary alloys [103, 104] (the so-called ‘pseudobinary’ approximation) may be unreliable in the multicomponent setting.

Zoom In Zoom Out Reset image size

For the Monte Carlo simulations which are to follow, we choose to truncate our fitted interactions and perform a fit limited to the first six coordination shells of the bcc lattice, corresponding to real-space radial ‘cutoffs’ of 6.226 Å and 6.050 Å for AlTiVNb and AlTiCrMo respectively. Such cutoffs are supported by the data shown in figure 6 and are also consistent with typical radial cutoffs chosen for machine-learned interatomic potentials [39]. These truncated EPIs are explicitly tabulated in the supplementary material [105] and are also provided in the open-access dataset associated with this study [96].

To validate the results of the above concentration wave analysis, and to further explore the configuration space below the initial chemical ordering temperatures, we perform lattice-based Monte Carlo simulations using the BraWl package [106] as outlined in section 2.3. Using the atom–atom EPIs obtained from the ab initio data and plotted in figure 6, we perform Monte Carlo simulations for both alloys using the Wang–Landau sampling algorithm [83]. The AlTiCrMo and AlTiVNb density of state histograms were obtained using a system consisting of bcc cubic unit cells for a total of 1024 atoms in an initially random configuration with periodic boundaries. The applied Wang–Landau tolerance was with a desired flatness of 80%. The AlTiCrMo simulation had 1024 uniform energy bins across total energy range of width 98 meV atom⁻¹. The AlTiVNb simulation had 1024 uniform energy bins across a total energy range of width 128 meV atom⁻¹. We do not consider total energies because the alloy theory evaluates derivatives of the alloy internal energy, and it is therefore not meaningful to consider total energies, only relative differences in energies between structures. The width (or, equivalently, number) of energy bins for the Wang–Landau simulations was chosen such that the algorithm was just over the verge of being able to easily bias the system out of energy bins. The energy range was chosen based on the minimum and maximum energies obtained through equilibrium Metropolis–Hastings Monte Carlo simulations at the highest and lowest temperatures of interest. Figure 7 shows plots of energy probability distribution histograms, Warren–Cowley ASRO parameters and SHC curves from our simulations, while figure 8 shows visualised configuration for AlTiCrMo and AlTiVNb obtained using Metropolis–Hastings Monte Carlo equilibration at selected temperatures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In AlTiVNb, with decreasing temperature, an initial peak in the SHC can be seen at approximately 2500 K, indicating a phase transition, with further transitions appearing below approximately 600 K. The initial ordering temperature of 2500 K is slightly reduced compared to the value of 2574 K predicted by the earlier concentration wave analysis, which is consistent with the fact that the concentration wave analysis provides a mean-field estimate of ordering temperatures. For the high temperature phase transition, the strongest ASRO is for Al–Al, Al–Ti, and Ti–Ti pairs. Al–Ti pairs are favoured on the first coordination shell and disfavoured on the second coordination shell, while Al–Al and Ti–Ti pairs are disfavoured on the first coordination shell and favoured on the second. This is indicative of a B2 ordering driven by Al and Ti, with Al moving to one simple cubic sublattice, Ti the other, and V and Nb expressing weaker site preferences. This analysis of the transition is supported by the visualised equilibrated configurations at 3000 K and 1200 K shown in figure 8. The peaks in the heat capacity observed below approximately 600 K are associated with additional atom–atom correlations, clustering tendencies, and eventual phase decomposition, as can be seen in the configuration equilibrated at low temperature in figure 8. These results suggest that heat treatment (annealing) of this alloy at temperatures between 1600 K and 1200 K should result in formation of a uniform B2-structured material.

These findings for AlTiVNb are consistent with [37], which reports a B2 chemical ordering occurring at approximately 1700 K, with Al and Ti expressing strong site preferences, and Nb and V expressing weaker preferences, in alignment with our own predictions. In addition, these authors report a secondary peak in their heat capacity data emerging around 600 K, consistent with our own findings. That our predicted B2 ordering temperature of 2300 K is higher than their prediction 1700 K we understand as originating in differences in lattice parameters and choices of XC functional. [37] accounts for thermal expansion and uses a lattice parameter of a = 3.23 Å, in combination with a GGA XC functional, whereas we use a DFT-optimised lattice parameter of a = 3.113 Å and an LDA XC functional. In previous work, we have found that a marginally contracted (expanded) lattice parameter results in an increased (decreased) chemical ordering temperature in our concentration wave analysis [53, 102], and it is also known that LDA XC functionals typically over-bind transition metals [99]. We understand these factors as the origin of the discrepancy between the two predicted ordering temperatures.

In AlTiCrMo, with decreasing temperature, a peak in the SHC can be seen at approximately 1800 K, indicating a phase transition, with a further transition occurring at approximately 400 K. The initial transition temperature of approximately 1800 K is reduced compared to the value of 2107 K predicted by the concentration wave analysis, which is again consistent with fact that the concentration wave analysis provides a mean-field estimate of chemical ordering temperatures. For the higher of the two ordering temperatures, the atom–atom correlations indicate that Al–Ti and Ti–Mo pairs are favoured on the first coordination shell, while Ti–Ti pairs are very heavily disfavoured. These pair preferences are consistent with the B2 ordering predicted by the earlier concentration wave analysis. Similarly to AlTiVNb, the transition occurring at lower simulation temperatures is associated with the emergence of additional sublattice atom–atom correlations and eventual phase decomposition. This is evidenced by the ASRO parameters on the second coordination shell showing that Al–Mo and Ti–Cr pairs are favoured, which is supported by the equilibrated low-temperature configuration shown in figure 8. These results suggest that heat treatment (annealing) of this alloy at temperatures between 1200 K and 800 K should result in formation of a uniform B2-structured material. That both AlTiVNb and AlTiCrMo eventually segregate into multiple competing phases at low temperatures emphasises that entropy plays an important role in stabilising these single-phase systems.

As an example of how chemical orderings, such as the B2 orderings predicted in this work for the AlTiVNb and AlTiCrMo RSAs, can affect materials’ properties, we now consider calculations of the residual resistivity for both alloys and examine how this intrinsic physical quantity is affected by the ordering. For the B2 orderings predicted for AlTiVNb and AlTiCrMo by the concentration wave analysis in section 3.2, we calculate the residual resistivity, ρ₀ as a function of atomic long-range order parameter, η, where η = 0 corresponds to the disordered solid solution, and η = 1 corresponds to the largest permitted chemical fluctuation consistent with the predicted chemical polarisation, i.e. the maximal B2 ordering. (Explicit partial lattice site occupancies as a function of η are provided in the supplementary material [105].) For these calculations, we again use the SPR-KKR package [60, 91], with broadly the same settings as for our self-consistent calculations in section 3.1, but now with increased k-point sampling (NKTAB = 150 000) and angular momentum cutoff (), as is necessary to capture the subtle details of the electronic scattering states around the Fermi level, . For both alloys, these resistivity data are plotted in figure 9.

Zoom In Zoom Out Reset image size

For the case η = 0, i.e. the A2 (disordered bcc) solid solution, AlTiVNb has a calculated residual resistivity of 100.4 cm, smaller than the value of 121.2 cm calculated for AlTiCrMo. (Note that the conductivity tensor is diagonal due to isotropy of the alloys with , resulting in an isotropic resistivity.) We attribute this difference in calculated resistivity to A2 AlTiCrMo having a reduced electronic DOS and increased presence of heavily smeared states around the when compared to A2 AlTiVNb, as evidenced in figure 3. For the case η = 1, i.e. the predicted partially ordered B2 structures, both alloys have an increased residual resistivity compared to the disordered A2 phase. For AlTiVNb the predicted residual resistivity increases to 128.7 cm, while for AlTiCrMo it increases to 158.4 cm. Between η = 0 and η = 1, the resistivity varies smoothly, monotonically increasing with increasing atomic order parameter. Although counter-intuitive—chemical orderings usually represent a restoration of a degree of translational symmetry and a consequent reduction in resistivity—such behaviour is consistent with earlier calculations demonstrating that a degree of ASRO can increase the residual resistivity of some transition metal alloys [108]. This is known as the ‘komplex’ or K-state phenomenon, and was first observed in Ni–Cr alloys [109]. In part, we believe that the increase in residual resistivity for both alloys can be attributed to a reduction in the electronic DOS at for the B2 ordered structures compared to the disordered A2 phase, as seen in figure 5. However, this does not represent the complete picture, and further analysis of the electronic structure of the alloys is required.

Understanding the underlying physical mechanisms governing the resistivity increase with increasing η within the Kubo–Greenwood formalism remains challenging, as this calculation yields only the components of the conductivity tensor. To gain a more intuitive understanding, one may employ the semi-classical Boltzmann transport equation within the relaxation-time approximation [110, 111]. In this framework, discussed in detail in the supplementary material [105], the zero temperature electrical conductivity is expressed in terms of the DOS at the Fermi level, , as

where e is the charge on the electron, represents the electronic group velocity component in the ith direction, is the k-dependent electronic relaxation time, and denotes an average taken over k lying in the first Brillouin zone with the energy fixed at . In addition, conductivity is directly influenced by k-space smearing of the spectral function and the corresponding reduction in the electronic mean free path [12].

A site- and orbital-resolved analysis of the electronic DOS—presented in the supplementary material [105]—reveals that, with the increase in residual resistivity, there is a concurrent reduction in the number of delocalised sp states at . This is in alignment with a simple picture of the residual resistivity being proportional to the DOS at the Fermi level, [67]. In particular, a decomposition of the sp DOS by atomic sublattices shows that, in AlTiVNb, the relative decrease of the sp states at on the 1b sublattice is in perfect agreement with the increase in relative resistivity. For AlTiCrMo, a similar trend is observed, albeit sublattice-inverted. We understand this decrease in available sp states at the 1b (1a) sites for AlTiVNb (AlTiCrMo) as originating from Al moving to the 1a (1b) lattice site in the predicted B2 ordering. Since, for both alloys, both sublattices contribute to the total residual resistivity, the asymmetric decrease in sp DOS at E_F between sublattices 1a and 1b in AlTiVNb and AlTiCrMo, respectively, suggests associated changes to the electronic mean free path. An analysis of the BSF along high-symmetry directions and an k-space cut through the Fermi surface resolved into its site-contributions supports this interpretation. As shown in the supplementary material [105], for AlTiVNb the 1a lattice exhibits a slight k-space broadening (smaller change of intensity in BSF with ) compared to the 1b lattice. For AlTiCrMo, however, the results are less clear without a direct calculation of the Fermi surface average of the electronic mean free path [12].

Thus, we conclude that the increase in residual resistivity with increasing B2 chemical order parameter, η, in the alloys studied is primarily governed by two mechanisms: (i) a site-specific reduction in the electronic DOS associated with sp bands at the Fermi level, which is the dominant effect, and (ii) a broadening or change of the spectral function in the k space, which leads to a shortening of the electronic mean free path, but which makes a less pronounced contribution and is more challenging to quantify.