# Effect of Ti addition on the structural, thermodynamic, and elastic properties of Tix(HfNbTaZr)(1-x)/4 alloys

Asif Iqbal Bhatti, Marwa Al-Houcine, Sylvain Queyreau, David Tingaud

► To cite this version:

Asif Iqbal Bhatti, Marwa Al-Houcine, Sylvain Queyreau, David Tingaud. Effect of Ti addition on the structural, thermodynamic, and elastic properties of Tix(HfNbTaZr)(1-x)/4 alloys. 2024. hal-03911775v2

HAL Id: hal-03911775

https://hal.science/hal-03911775v2

Preprint submitted on 21 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Effect of Ti addition on the structural, thermodynamic, and elastic properties of Tix(HfNbTaZr)(1-x)/4 alloys

Asif Iqbal Bhatti, Marwa Al-Houcine, Sylvain Queyreau, David Tingaud

► To cite this version:

Asif Iqbal Bhatti, Marwa Al-Houcine, Sylvain Queyreau, David Tingaud. Effect of Ti addition on the structural, thermodynamic, and elastic properties of Tix(HfNbTaZr)(1-x)/4 alloys. 2024. hal-03911775v2

HAL Id: hal-03911775

https://hal.science/hal-03911775v2

Preprint submitted on 21 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Effect of Ti addition on the structural,
thermodynamic, and elastic properties of
$Ti_x(HfNbTaZr)_{(1-x)/4}$ alloys

Asif Iqbal Bhatti¹, Marwa Al-Houcine¹, David Tingaud¹, and Sylvain

Queyreau*¹

¹Université Sorbonne Paris Nord, Laboratoire des Sciences des
Procédés et des Matériaux, LSPM, CNRS, UPR 3407, 99 avenue
Jean-Baptiste Clément, F-93430 Villetaneuse, France

June 21, 2024

**Abstract**

The structure, thermodynamic and elastic properties of $Ti_x(HfNbTaZr)_{(1-x)/4}$ from Refractory High Entropy multicomponent Alloys to pure titanium are investigated employing comprehensive MCSQS realizations of the disordered atomic structure and DFT calculations. We showed that to model the random structure in a limited supercell, it is beneficial to probe a large space of random configurations with different atomic arrangements. Including a larger number

<footnote>*Corresponding author: sylvain.queyreau@lspm.cnrs.fr</footnote>
Effect of Ti addition on the structural,
thermodynamic, and elastic properties of
$Ti_x(HfNbTaZr)_{(1-x)/4}$ alloys

Asif Iqbal Bhatti¹, Marwa Al-Houcine¹, David Tingaud¹, and Sylvain

Queyreau*¹

¹Université Sorbonne Paris Nord, Laboratoire des Sciences des
Procédés et des Matériaux, LSPM, CNRS, UPR 3407, 99 avenue
Jean-Baptiste Clément, F-93430 Villetaneuse, France

June 21, 2024

**Abstract**

The structure, thermodynamic and elastic properties of $Ti_x(HfNbTaZr)_{(1-x)/4}$ from Refractory High Entropy multicomponent Alloys to pure titanium are investigated employing comprehensive MCSQS realizations of the disordered atomic structure and DFT calculations. We showed that to model the random structure in a limited supercell, it is beneficial to probe a large space of random configurations with different atomic arrangements. Including a larger number

<footnote>*Corresponding author: sylvain.queyreau@lspm.cnrs.fr</footnote>
of nearest-neighbor pair shells is in practice more effective in finding configurations corresponding to lower formation energies than considering many-body terms. In addition, we demonstrate the existence of weak-to-medium short-range chemical order for both equimolar compositions. The study of chemical order requires hundreds of MCSQS realizations to become meaningful, and SRO results are rationalized in terms of crystallographic structure of favorable element pairs and binary phase diagrams. When Ti is added to Tiₓ(HfNbTaZr)_(1-x)/4 alloys, the mixing enthalpy of the bcc phase decreases but remains slightly positive for all x. For large Ti content, a phase transition in favor of an hcp structure is suggested in agreement with the predictions of the Bo-Md diagram and mechanical stability analysis. Around x ≈ 0.50, a dual phase is predicted. The impact of local atomic configuration variability on elastic properties and bcc/hcp phase competition is discussed. Our results are in good agreement with the limited existing literature on these alloys, and provide a better understanding of it. The Ti content in this class of alloys could be a practical means of selecting the phase structure and tailoring the elastic properties to specific applications.

**Keywords**— High entropy alloys, Random alloys, Special Quasi-Random structures (SQSs), Density Functional Theory (DFT), mixing entropy, Phase stability, Body-centered cubic (bcc), Metastable alloys.

# 1 Introduction

Multicomponent High Entropy Alloys (HEAs), a relatively new class of materials, have recently attracted a great deal of attention due to their remarkable mechanical properties. These HEAs are typically composed of four or more elements at nearly the same concentration [1–3]. Despite consisting of a complex multicomponent system, they lead to the
of nearest-neighbor pair shells is in practice more effective in finding configurations corresponding to lower formation energies than considering many-body terms. In addition, we demonstrate the existence of weak-to-medium short-range chemical order for both equimolar compositions. The study of chemical order requires hundreds of MCSQS realizations to become meaningful, and SRO results are rationalized in terms of crystallographic structure of favorable element pairs and binary phase diagrams. When Ti is added to Ti<sub>x</sub>(HfNbTaZr)<sub>(1-x)</sub> alloys, the mixing enthalpy of the bcc phase decreases but remains slightly positive for all x. For large Ti content, a phase transition in favor of an hcp structure is suggested in agreement with the predictions of the Bo-Md diagram and mechanical stability analysis. Around x ≈ 0.50, a dual phase is predicted. The impact of local atomic configuration variability on elastic properties and bcc/hcp phase competition is discussed. Our results are in good agreement with the limited existing literature on these alloys, and provide a better understanding of it. The Ti content in this class of alloys could be a practical means of selecting the phase structure and tailoring the elastic properties to specific applications.

**Keywords**— High entropy alloys, Random alloys, Special Quasi-Random structures (SQSs), Density Functional Theory (DFT), mixing entropy, Phase stability, Body-centered cubic (bcc), Metastable alloys.

# 1 Introduction

Multicomponent High Entropy Alloys (HEAs), a relatively new class of materials, have recently attracted a great deal of attention due to their remarkable mechanical properties. These HEAs are typically composed of four or more elements at nearly the same concentration [1–3]. Despite consisting of a complex multicomponent system, they lead to the
formation of a simple single-phase solid solution, such as face-centered cubic (fcc), body-
centered cubic (bcc), or hexagonal close-packed (hcp). Among these alloys, fcc *Cantor*
*type* HEAs have been widely studied, but few studies exist on bcc HEAs, and even fewer
in the case of bcc refractory HEAs (RHEAs) [4, 5], sometimes named "Senkov" alloys.
The high degree of chemical complexity created by mixing different elements with different
atomic sizes and chemical properties increases the material's resistance to deformation and
fracture. In addition, the high entropy of the alloy leads to a variety of lattice distortions
due to atomic configurations. This arrangement of atoms promotes strain hardening and
increases resistance to dislocation movements [6–10]. Despite abundant literature, many
aspects of phase stability, exact atomic structure, and their relationship to macroscopic
properties remain unclear.

Recent experiments on HfNbTaTi, HfMoTaTiZr, and HfNbTaTiZr have shown that non-toxic and non-allergenic elements such as Nb, Ta, and Zr exhibit high corrosion resistance [11]. The latter alloy family, HfNbTaTiZr, will be the focus of the present paper. Properties such as high ductility, corrosion resistance, and strength at both room temperature and high temperatures can be improved by choosing suitable refractory elements (a mixture of bcc and hcp elements).

The *d*-electron theory has become a classical tool for predicting the structural stabil-
ity of Ti-based alloys [12] but it has rarely been applied to high-entropy Ti-based alloys.
Lilenstein *et al.* have synthesized a novel high entropy alloy (HEA) with the composition
Ti₃₅Zr₂₇.₅Hf₂₇.₅Nb₅Ta₅ using the *d*-electron alloy design approach based on the Bo-Md
diagram [13]. This alloy exhibited a remarkable transformation-induced plasticity effect,
resulting in a high normalized work hardening rate while retaining ductility when compared
to the reference equimolar Ti₂₀Zr₂₀Hf₂₀Nb₂₀Ta₂₀ alloy. Microstructure analysis after de-
formation revealed the large presence of stress-induced orthorhombic martensite. This
work shows that the phase stability and mechanical properties of HEA can be tailored by
a chemical design approach based on Bo and Md parameters. However, many of the pre-
formation of a simple single-phase solid solution, such as face-centered cubic (fcc), body-
centered cubic (bcc), or hexagonal close-packed (hcp). Among these alloys, fcc *Cantor*
*type* HEAs have been widely studied, but few studies exist on bcc HEAs, and even fewer
in the case of bcc refractory HEAs (RHEAs) [4, 5], sometimes named "Senkov" alloys.
The high degree of chemical complexity created by mixing different elements with different
atomic sizes and chemical properties increases the material's resistance to deformation and
fracture. In addition, the high entropy of the alloy leads to a variety of lattice distortions
due to atomic configurations. This arrangement of atoms promotes strain hardening and
increases resistance to dislocation movements [6–10]. Despite abundant literature, many
aspects of phase stability, exact atomic structure, and their relationship to macroscopic
properties remain unclear.

Recent experiments on HfNbTaTi, HfMoTaTiZr, and HfNbTaTiZr have shown that non-toxic and non-allergenic elements such as Nb, Ta, and Zr exhibit high corrosion resistance [11]. The latter alloy family, HfNbTaTiZr, will be the focus of the present paper. Properties such as high ductility, corrosion resistance, and strength at both room temperature and high temperatures can be improved by choosing suitable refractory elements (a mixture of bcc and hcp elements).

The *d*-electron theory has become a classical tool for predicting the structural stabil-
ity of Ti-based alloys [12] but it has rarely been applied to high-entropy Ti-based alloys.
Lilenstein *et al.* have synthesized a novel high entropy alloy (HEA) with the composition
Ti₃₅Zr₂₇.₅Hf₂₇.₅Nb₅Ta₅ using the *d*-electron alloy design approach based on the Bo-Md
diagram [13]. This alloy exhibited a remarkable transformation-induced plasticity effect,
resulting in a high normalized work hardening rate while retaining ductility when compared
to the reference equimolar Ti₂₀Zr₂₀Hf₂₀Nb₂₀Ta₂₀ alloy. Microstructure analysis after de-
formation revealed the large presence of stress-induced orthorhombic martensite. This
work shows that the phase stability and mechanical properties of HEA can be tailored by
a chemical design approach based on Bo and Md parameters. However, many of the pre-
dictions of the Bo-Md diagram need to be compared with the actual HEA phases observed
in experiments or simulations.

Departing from the initial picture of the random distribution of atoms, it is now clear
that HEA alloys exhibit local chemical fluctuations and various degrees of chemical short-
range ordering (SRO), and assessing them is crucial to understanding the microstructure-
property relationship. However, characterizing and quantifying composition fluctuations
is still a challenging task at the atomic scale, and no direct experimental observation
was available until very recently [14–16]. Yin *et al.* were certainly the first to show that
Monte Carlo simulations can thus provide an alternative route to study SRO in association
with DFT calculations [17], MC was then associated with cluster expansion model [18] or
molecular dynamics simulations [19], but there are still only a few studies of this type. An
increasing degree of SRO was found to lower the mixing energy of the HEA phase in the
corresponding alloys, while some alloys seem more prone to SRO than others. In the qua-
ternary bcc HfNbTiZr alloy, Hf-rich and Ti-rich clusters were experimentally highlighted
using atomic probe tomography and HR-HAADF-SEM observations [14, 15], in agreement
with the affinity of some of the Ti- and Hf-bond types observed in simulations [20, 21]. The
electronic density structure analysis performed on the same alloy [21] highlighted the pecu-
liar role of Ti-Zr bonds with a reduction of d electrons at the Fermi level, which stabilizes
the bcc phase.

Another fascinating aspect of HEA is to fully understand the origin of their large
strength, and a number of models have been successful in estimating the yield stress of
fcc and bcc alloys [22–24]. In deformed HfNbTaTiZr, TEM observations [25, 26] revealed
heterogeneous dislocation microstructures, with dense slip bands along {110} and {112}
planes and many debris loops. $\frac{a_0}{2}\langle 111 \rangle$ screw dislocations are mostly rectilinear, and the
large apparent activation volumes are consistent with the Peierls mechanism and large
friction with the lattice. Once again, the variation of the local chemical environment and
SRO existence seem key in rationalizing the dislocation behavior. EXAFS analysis of
dictions of the Bo-Md diagram need to be compared with the actual HEA phases observed
in experiments or simulations.

Departing from the initial picture of the random distribution of atoms, it is now clear
that HEA alloys exhibit local chemical fluctuations and various degrees of chemical short-
range ordering (SRO), and assessing them is crucial to understanding the microstructure-
property relationship. However, characterizing and quantifying composition fluctuations
is still a challenging task at the atomic scale, and no direct experimental observation
was available until very recently [14–16]. Yin *et al.* were certainly the first to show that
Monte Carlo simulations can thus provide an alternative route to study SRO in association
with DFT calculations [17], MC was then associated with cluster expansion model [18] or
molecular dynamics simulations [19], but there are still only a few studies of this type. An
increasing degree of SRO was found to lower the mixing energy of the HEA phase in the
corresponding alloys, while some alloys seem more prone to SRO than others. In the qua-
ternary bcc HfNbTiZr alloy, Hf-rich and Ti-rich clusters were experimentally highlighted
using atomic probe tomography and HR-HAADF-SEM observations [14, 15], in agreement
with the affinity of some of the Ti- and Hf-bond types observed in simulations [20, 21]. The
electronic density structure analysis performed on the same alloy [21] highlighted the pecu-
liar role of Ti-Zr bonds with a reduction of d electrons at the Fermi level, which stabilizes
the bcc phase.

Another fascinating aspect of HEA is to fully understand the origin of their large
strength, and a number of models have been successful in estimating the yield stress of
fcc and bcc alloys [22–24]. In deformed HfNbTaTiZr, TEM observations [25, 26] revealed
heterogeneous dislocation microstructures, with dense slip bands along {110} and {112}
planes and many debris loops. $\frac{a_0}{2}\langle 111 \rangle$ screw dislocations are mostly rectilinear, and the
large apparent activation volumes are consistent with the Peierls mechanism and large
friction with the lattice. Once again, the variation of the local chemical environment and
SRO existence seem key in rationalizing the dislocation behavior. EXAFS analysis of
pure binary and ternary alloys in the Ti-Ta-Hf-Nb-Zr system evidenced local chemical and
lattice distortions [27, 28]. MD investigations [29] of the dislocation structure in the Nb-
Ti-Zr system showed core structure variations along the dislocation line depending on the
local environment. A collision of cross kinks nucleated on different slip planes yields debris
(pairs of interstitial vacancies) in the wake of the dislocation gliding. Yin and collab. [23]
quantified the Peierls barrier for screw dislocations in the bcc MoNbTaW alloy by ab-initio.
The dislocation core energy roughly follows a Gaussian distribution depending on the local
environment of the dislocation, whose average and variance increase with the degree of
SRO. As a consequence, two types of Peierls potentials were identified, and the Peierls
barrier was found to increase with the degree of SRO. Finally, in-situ TEM observations
revealed [15] that local chemical fluctuations in the quaternary bcc HfNbTiZr alloy lead
to dislocation pinning, stimulating dislocation multiplication and cross-slip activity.

Ti content is expected to strongly impact the stability of the different possible phases
and local chemical fluctuations observed in the Hf-Nb-Ta-Ti-Zr system. In this paper, we
systematically investigate the effects of Ti content on the atomic structure, phase stability,
local chemical ordering, and elastic constants using Monte-Carlo Special Quasi random
(MCSQS) and ab-initio calculations. We employ DFT calculations as the derivation of re-
liable semi-empirical potentials for quinary alloys is still a challenging task [30]. To better
approximate the ideal random structure in a small DFT supercell, we considered multiple
atomic structures resulting from different MC SQS calculations, for each of which the en-
ergy was calculated. Interestingly, the average mixing energy and Global SRO parameter
converge rapidly to constant values with the number of different atomic configurations,
while the SRO parameters for each pair type require many more SQS realizations to con-
verge toward the ideal random structure where no pairs are favored. The comprehensive
set of data can also be further analyzed to identify favorable, neutral, and unfavorable
bond types corresponding to extremes in the mixing energy data. This solution consti-
tutes an alternative to a more complete but more demanding MC-DFT hybrid approach.
pure binary and ternary alloys in the Ti-Ta-Hf-Nb-Zr system evidenced local chemical and
lattice distortions [27, 28]. MD investigations [29] of the dislocation structure in the Nb-
Ti-Zr system showed core structure variations along the dislocation line depending on the
local environment. A collision of cross kinks nucleated on different slip planes yields debris
(pairs of interstitial vacancies) in the wake of the dislocation gliding. Yin and collab. [23]
quantified the Peierls barrier for screw dislocations in the bcc MoNbTaW alloy by ab-initio.
The dislocation core energy roughly follows a Gaussian distribution depending on the local
environment of the dislocation, whose average and variance increase with the degree of
SRO. As a consequence, two types of Peierls potentials were identified, and the Peierls
barrier was found to increase with the degree of SRO. Finally, in-situ TEM observations
revealed [15] that local chemical fluctuations in the quaternary bcc HfNbTiZr alloy lead
to dislocation pinning, stimulating dislocation multiplication and cross-slip activity.

Ti content is expected to strongly impact the stability of the different possible phases
and local chemical fluctuations observed in the Hf-Nb-Ta-Ti-Zr system. In this paper, we
systematically investigate the effects of Ti content on the atomic structure, phase stability,
local chemical ordering, and elastic constants using Monte-Carlo Special Quasi random
(MCSQS) and ab-initio calculations. We employ DFT calculations as the derivation of re-
liable semi-empirical potentials for quinary alloys is still a challenging task [30]. To better
approximate the ideal random structure in a small DFT supercell, we considered multiple
atomic structures resulting from different MC SQS calculations, for each of which the en-
ergy was calculated. Interestingly, the average mixing energy and Global SRO parameter
converge rapidly to constant values with the number of different atomic configurations,
while the SRO parameters for each pair type require many more SQS realizations to con-
verge toward the ideal random structure where no pairs are favored. The comprehensive
set of data can also be further analyzed to identify favorable, neutral, and unfavorable
bond types corresponding to extremes in the mixing energy data. This solution consti-
tutes an alternative to a more complete but more demanding MC-DFT hybrid approach.
In agreement with the Bo-Md diagram, tuning the Ti content in HfNbTaTiZr alloys has the potential to alter the alloy's stable phase. The knowledge gained from these investigations can be used to design high-entropy bcc alloys with properties tailored for specific applications. The investigation of the screw dislocation structure and property in correlation to its local chemical environment is left for a forthcoming publication. The remainder of the paper is organized as follows: Sec. 2 details the simulation conditions and a sensitivity study of the MCSQS parameters to construct the atomistic configurations of the alloys. Sec. 3 provides results and discussion regarding the SRO existing in quaternary and quinary equimolar alloys, the convergence of some average alloy properties with the number of MCSQS realizations, the relative stability of the different phases, the elastic properties as a function of Ti content, and possible correlations with local atomic arrangements. The last section summarizes the conclusions.

## 2 Methodology and Computational details

In the present article, we report the mixing energy, elastic properties, and structural parameters of a model of the Ti<sub>x</sub>(HfNbTaZr)<sub>(1-x)</sub> alloys calculated using Density Functional Theory (DFT) as implemented in Vienna Ab-initio Simulation Package (VASP) code [31]. The Projector Augmented Wave (PAW) pseudopotential was used to approximate the electron core energy [32] and the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation was used to approximate the exchange-correlation functional (XC) [33].

A lot of care has been paid to the construction of the disordered structure corresponding to the Ti<sub>x</sub>(HfNbTaZr)<sub>(1-x)</sub> alloys. The composition was varied from x = 0.0 to 0.81 by steps of about 0.1, while the rest of the composition is equally distributed among all the other elements Hf, Nb, Ta, and Zr. For all considered compositions, the supercell geometry was constrained so that the cubic bcc symmetry was preserved, as experiments suggest that most compositions correspond to a bcc single phase. Additional calculations using an hcp,
In agreement with the Bo-Md diagram, tuning the Ti content in HfNbTaTiZr alloys has the potential to alter the alloy's stable phase. The knowledge gained from these investigations can be used to design high-entropy bcc alloys with properties tailored for specific applications. The investigation of the screw dislocation structure and property in correlation to its local chemical environment is left for a forthcoming publication. The remainder of the paper is organized as follows: Sec. 2 details the simulation conditions and a sensitivity study of the MCSQS parameters to construct the atomistic configurations of the alloys. Sec. 3 provides results and discussion regarding the SRO existing in quaternary and quinary equimolar alloys, the convergence of some average alloy properties with the number of MCSQS realizations, the relative stability of the different phases, the elastic properties as a function of Ti content, and possible correlations with local atomic arrangements. The last section summarizes the conclusions.

## 2 Methodology and Computational details

In the present article, we report the mixing energy, elastic properties, and structural parameters of a model of the Ti<sub>x</sub>(HfNbTaZr)<sub>(1-x)</sub> alloys calculated using Density Functional Theory (DFT) as implemented in Vienna Ab-initio Simulation Package (VASP) code [31]. The Projector Augmented Wave (PAW) pseudopotential was used to approximate the electron core energy [32] and the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation was used to approximate the exchange-correlation functional (XC) [33].

A lot of care has been paid to the construction of the disordered structure corresponding to the Ti<sub>x</sub>(HfNbTaZr)<sub>(1-x)</sub> alloys. The composition was varied from x = 0.0 to 0.81 by steps of about 0.1, while the rest of the composition is equally distributed among all the other elements Hf, Nb, Ta, and Zr. For all considered compositions, the supercell geometry was constrained so that the cubic bcc symmetry was preserved, as experiments suggest that most compositions correspond to a bcc single phase. Additional calculations using an hcp,
or orthorhombic (ortho) lattice, have also been conducted for high Ti content. Mimicking
a random structure in a supercell of sufficiently limited dimensions to be tractable by DFT
simulations raises an additional problem due to short-range spatial correlations among
atoms induced by the Periodic Boundary Conditions (PBC). To address this problem,
several techniques have been proposed in the literature based on the use of order parameters
to evaluate the disordered nature of the atomic configuration and minimize the spatial
correlations as much as possible [34–38]. Here, we use the Monte-Carlo Special Quasi-
random Structure (MCSQS) technique implemented in the Alloy Theoretic Automated
Toolkit (ATAT) code [39,40] to generate and analyze a large set of different random structure
realizations. We can thus perform a statistical study of the MCSQS parameters on the
resulting structure of different atomic configurations.

The supercell size was chosen to be a necessary tradeoff between physical relevance and
the numerical cost of the calculations. The supercell geometries depend on the investigated
composition and are as follows: For *x* = 0.11, 0.33, 0.40, 0.63 and 0.70, we generated a
supercell of 54 atoms based on 3 × 3 × 3 cubic bcc unit cells; for compositions corresponding
to *x* = 0.00, 0.50, and 0.81, a supercell of 128 atoms based on 4 × 4 × 4 cubic bcc unit
cells was used; and for *x* = 0.20, a supercell of 125 atoms based on 5 × 5 × 5 primitive bcc
unit cells was used. In the case of hcp structures, we employed supercells of 128 atoms for
*x* = 0.50, 0.81 (4 × 4 × 4 hexagonal unit cells), and 54 atoms for *x* = 0.63, 0.70 (3 × 3 × 3
hexagonal unit cells).

In all the reported simulations, the system was allowed to fully relax (to accommodate
the lattice distortions typical of HEA alloys), including atomic positions, supercell shapes,
and volumes (in an anisotropic fashion), using a conjugate gradient algorithm. Brillouin
zone integrations were performed with a Γ-centered k-point grid scheme. The plane-wave
cut-off energy was set to 600 eV. The first order Methfessel-Paxton method [41] with a
smearing parameter of 0.2 eV was used. The system was assumed to have converged
when the force on each atom was less than 1 meV/Å, and the stress tolerance was below
or orthorhombic (ortho) lattice, have also been conducted for high Ti content. Mimicking
a random structure in a supercell of sufficiently limited dimensions to be tractable by DFT
simulations raises an additional problem due to short-range spatial correlations among
atoms induced by the Periodic Boundary Conditions (PBC). To address this problem,
several techniques have been proposed in the literature based on the use of order parameters
to evaluate the disordered nature of the atomic configuration and minimize the spatial
correlations as much as possible [34–38]. Here, we use the Monte-Carlo Special Quasi-
random Structure (MCSQS) technique implemented in the Alloy Theoretic Automated
Toolkit (ATAT) code [39,40] to generate and analyze a large set of different random structure
realizations. We can thus perform a statistical study of the MCSQS parameters on the
resulting structure of different atomic configurations.

The supercell size was chosen to be a necessary tradeoff between physical relevance and
the numerical cost of the calculations. The supercell geometries depend on the investigated
composition and are as follows: For *x* = 0.11, 0.33, 0.40, 0.63 and 0.70, we generated a
supercell of 54 atoms based on 3 × 3 × 3 cubic bcc unit cells; for compositions corresponding
to *x* = 0.00, 0.50, and 0.81, a supercell of 128 atoms based on 4 × 4 × 4 cubic bcc unit
cells was used; and for *x* = 0.20, a supercell of 125 atoms based on 5 × 5 × 5 primitive bcc
unit cells was used. In the case of hcp structures, we employed supercells of 128 atoms for
*x* = 0.50, 0.81 (4 × 4 × 4 hexagonal unit cells), and 54 atoms for *x* = 0.63, 0.70 (3 × 3 × 3
hexagonal unit cells).

In all the reported simulations, the system was allowed to fully relax (to accommodate
the lattice distortions typical of HEA alloys), including atomic positions, supercell shapes,
and volumes (in an anisotropic fashion), using a conjugate gradient algorithm. Brillouin
zone integrations were performed with a Γ-centered k-point grid scheme. The plane-wave
cut-off energy was set to 600 eV. The first order Methfessel-Paxton method [41] with a
smearing parameter of 0.2 eV was used. The system was assumed to have converged
when the force on each atom was less than 1 meV/Å, and the stress tolerance was below
0.005 GPa. After relaxation, the mixing energy was evaluated using a single-point energy calculation with a dense k-point relative to the convergence of the mixing energy. A sensitivity study showed that a tolerance of 1 meV/atom, a 2 × 2 × 2 k-point grid, and a plane-wave cut-off energy of 600 eV were sufficient for the calculation of elastic and structural properties.

# 3 Results and Discussion

We employed the MCSQS approach to mimic the random structure of HEA in a DFT supercell. While MCSQS is now a well-established technique within the HEA community, its results can be impacted by a few parameters, in particular, the number of neighbor shells and the order of the many-body terms that are included to evaluate the disorder in the considered configuration. We thus start with a sensitivity study of the structure generated as a function of these parameters.

In the present work, we have varied the pairwise terms (hereafter called *p*) between the NN (Nearest Neighbor) shells, where *p* ranges from 2 to 10 NN. We have also varied the three-body (*t*) and four-body (*q*) terms (where *t* and *q* range from 0 to 1 NN).

<figure><img src="image_2.png" /><figcaption>Figure 1: Examples of supercells [(a)-(d)] obtained from MCSQS for the alloy Ti<sub>x</sub>(HfNbTaZr)<sub>(1-x)</sub>. Figs. (a) and (c) show a 3 × 3 × 3 supercell for *x* = 0.11, 0.33, 0.40, 0.60, and 0.70 containing 54 atoms with bcc and hcp structures, respectively. Fig. (b) and (d) show a 4 × 4 × 4 supercell for *x* = 0.00, 0.50, and 0.80, representing 128 atoms with a bcc and hcp structure, respectively.</figcaption></figure>

For each Ti content *x* and MCSQS parameters *p*, *t*, and *q* (hereafter denoted by the value triplet (*p*, *t*, *q*)), Fig. 1 shows some of the atomic configurations resulting from the MCSQS calculations obtained for different Ti content *x*. Our workflow consists of
0.005 GPa. After relaxation, the mixing energy was evaluated using a single-point energy calculation with a dense k-point relative to the convergence of the mixing energy. A sensitivity study showed that a tolerance of 1 meV/atom, a 2 × 2 × 2 k-point grid, and a plane-wave cut-off energy of 600 eV were sufficient for the calculation of elastic and structural properties.

# 3 Results and Discussion

We employed the MCSQS approach to mimic the random structure of HEA in a DFT supercell. While MCSQS is now a well-established technique within the HEA community, its results can be impacted by a few parameters, in particular, the number of neighbor shells and the order of the many-body terms that are included to evaluate the disorder in the considered configuration. We thus start with a sensitivity study of the structure generated as a function of these parameters.

In the present work, we have varied the pairwise terms (hereafter called *p*) between the NN (Nearest Neighbor) shells, where *p* ranges from 2 to 10 NN. We have also varied the three-body (*t*) and four-body (*q*) terms (where *t* and *q* range from 0 to 1 NN).

<figure><img src="image_2.png" /><figcaption>Figure 1: Examples of supercells [(a)-(d)] obtained from MCSQS for the alloy Ti<sub>x</sub>(HfNbTaZr)<sub>(1-x)</sub>. Figs. (a) and (c) show a 3 × 3 × 3 supercell for *x* = 0.11, 0.33, 0.40, 0.60, and 0.70 containing 54 atoms with bcc and hcp structures, respectively. Fig. (b) and (d) show a 4 × 4 × 4 supercell for *x* = 0.00, 0.50, and 0.80, representing 128 atoms with a bcc and hcp structure, respectively.</figcaption></figure>

For each Ti content *x* and MCSQS parameters *p*, *t*, and *q* (hereafter denoted by the value triplet (*p*, *t*, *q*)), Fig. 1 shows some of the atomic configurations resulting from the MCSQS calculations obtained for different Ti content *x*. Our workflow consists of
first generating 20 SQS for (p, 0, 0) only, then with (p, t, 0), and finally (p, t, q). In each case, the MCSQS calculation was initiated with different random seed values, leading to different final atomic configurations. The configuration was assumed to be final when the local minimum of the objective function was reached during the 120-hour runtime. MCSQS was terminated earlier when no further SQS structure was generated for more than 12 hours.

According to classical theory, the Gibbs free energy of the system is given by:

$$ \Delta G_x = \Delta H_x - (T\Delta S_{mix})_x - (T\Delta S_{vib})_x \quad (1) $$

In Eq. 1, $\Delta H_x$ is the mixing energy, and $\Delta S_{mix}$ and $\Delta S_{vib}$ are the mixing and vibrational entropy terms, respectively. The mixing energy, $\Delta H_x$, given in Eq. 2 is computed from the difference of energy of the SQS supercell ($E_{SQS}$) and the energy $E_i$ associated with the most stable phase of each element, i (bcc for Nb, and Ta and hcp for Hf, Ti, and Zr):

$$ \Delta H_x = (E_{SQS})_x - \left( \sum_i c_i E_i \right)_x \quad (2) $$

where $c_i$ is the atomic fraction of an alloy element, i. For alloys, the mixing entropic term is as follows:

$$ (\Delta S_{mix})_x = -k_B \sum_i c_i \ln(c_i) \quad (3) $$

where $k_B$ is the Boltzmann constant. It can be noted that Eq. 3 applies to perfect random solid solutions, which might not be the case for HEA alloys, where some level of chemical order may exist (see later). Eq. 3 is nevertheless commonly used within the HEA community to estimate the mixing entropy. The last term of Eq. 1, corresponding to the vibrational entropy, can be calculated from DFT using the so-called frozen phonon approach [42,43]. However, these types of calculations are computationally demanding, and
first generating 20 SQS for (p, 0, 0) only, then with (p, t, 0), and finally (p, t, q). In each case, the MCSQS calculation was initiated with different random seed values, leading to different final atomic configurations. The configuration was assumed to be final when the local minimum of the objective function was reached during the 120-hour runtime. MCSQS was terminated earlier when no further SQS structure was generated for more than 12 hours.

According to classical theory, the Gibbs free energy of the system is given by:

$$ \Delta G_x = \Delta H_x - (T\Delta S_{mix})_x - (T\Delta S_{vib})_x \quad (1) $$

In Eq. 1, $\Delta H_x$ is the mixing energy, and $\Delta S_{mix}$ and $\Delta S_{vib}$ are the mixing and vibrational entropy terms, respectively. The mixing energy, $\Delta H_x$, given in Eq. 2 is computed from the difference of energy of the SQS supercell ($E_{SQS}$) and the energy $E_i$ associated with the most stable phase of each element, i (bcc for Nb, and Ta and hcp for Hf, Ti, and Zr):

$$ \Delta H_x = (E_{SQS})_x - \left( \sum_i c_i E_i \right)_x \quad (2) $$

where $c_i$ is the atomic fraction of an alloy element, i. For alloys, the mixing entropic term is as follows:

$$ (\Delta S_{mix})_x = -k_B \sum_i c_i \ln(c_i) \quad (3) $$

where $k_B$ is the Boltzmann constant. It can be noted that Eq. 3 applies to perfect random solid solutions, which might not be the case for HEA alloys, where some level of chemical order may exist (see later). Eq. 3 is nevertheless commonly used within the HEA community to estimate the mixing entropy. The last term of Eq. 1, corresponding to the vibrational entropy, can be calculated from DFT using the so-called frozen phonon approach [42,43]. However, these types of calculations are computationally demanding, and
vibrational entropy is commonly neglected when compared to mixing entropy [1, 37, 44, 45].

This last term is thus disregarded in the present study.

## 3.1 Sensitivity study on MCSQS parameters

First, we conducted MCSQS calculations while neglecting the many-body terms, i.e., we included only the pairwise *p* NN shells, as these terms are expected to have the largest impact on the resulting SQS configuration. We then determined the effect of the number of *p* shells on the mixing energy. To evaluate this parameter, we have considered the *x* = 0.33 composition. Tab. 1 shows the average mixing energy and average lattice parameter as a function of *p* shells. The values given in Tab. 1 refer to the lowest mixing energy obtained over a relevant number of different atomic configurations. We will see later that 20 different SQS realizations are sufficient to assess the ideal solid solution mixing energy for the alloys and supercell sizes considered here. Interestingly, the mixing energy per atom decreases only weakly as the number of *p* NN shells increases up to 8 NN. For *p* = 10 NN, it surprisingly yielded a larger mixing energy per atom by about 7 meV/atom. In the latter case, it is probable that modeling SQS supercells with 10 pairs of shells within a 120-h time limit is not enough, or it might be that considering a larger number of shells within a finite supercell constrained by the given composition and size of the supercell overlaps due to periodic boundary conditions (PBC). As a result of this sensitivity study, we chose *p* = 8 NN as a compromise between MCSQS capacity to produce low energy configurations and the numerical cost to run the MCSQS.

In another set of calculations, we have investigated the possible effect of higher order many-body terms: ternary (t) and quaternary (q) interactions within 1st NN shells from MCSQS. For this, many-body terms were set to *t* = 1 NN, and *q* = 0 or 1 NN, while the pairwise terms were set to 6 and 8 NN shells. Fig. 2 shows the impact of including these higher order many-body terms on the mixing energy of $Ti_x(HfNbTaZr)_{(1-x)/4}$ for *x* ranging from 0.00 to 0.50. The figure reports the average (avg) mixing energy value, along with
vibrational entropy is commonly neglected when compared to mixing entropy [1, 37, 44, 45].

This last term is thus disregarded in the present study.

## 3.1 Sensitivity study on MCSQS parameters

First, we conducted MCSQS calculations while neglecting the many-body terms, i.e., we included only the pairwise *p* NN shells, as these terms are expected to have the largest impact on the resulting SQS configuration. We then determined the effect of the number of *p* shells on the mixing energy. To evaluate this parameter, we have considered the *x* = 0.33 composition. Tab. 1 shows the average mixing energy and average lattice parameter as a function of *p* shells. The values given in Tab. 1 refer to the lowest mixing energy obtained over a relevant number of different atomic configurations. We will see later that 20 different SQS realizations are sufficient to assess the ideal solid solution mixing energy for the alloys and supercell sizes considered here. Interestingly, the mixing energy per atom decreases only weakly as the number of *p* NN shells increases up to 8 NN. For *p* = 10 NN, it surprisingly yielded a larger mixing energy per atom by about 7 meV/atom. In the latter case, it is probable that modeling SQS supercells with 10 pairs of shells within a 120-h time limit is not enough, or it might be that considering a larger number of shells within a finite supercell constrained by the given composition and size of the supercell overlaps due to periodic boundary conditions (PBC). As a result of this sensitivity study, we chose *p* = 8 NN as a compromise between MCSQS capacity to produce low energy configurations and the numerical cost to run the MCSQS.

In another set of calculations, we have investigated the possible effect of higher order many-body terms: ternary (t) and quaternary (q) interactions within 1st NN shells from MCSQS. For this, many-body terms were set to *t* = 1 NN, and *q* = 0 or 1 NN, while the pairwise terms were set to 6 and 8 NN shells. Fig. 2 shows the impact of including these higher order many-body terms on the mixing energy of $Ti_x(HfNbTaZr)_{(1-x)/4}$ for *x* ranging from 0.00 to 0.50. The figure reports the average (avg) mixing energy value, along with
<table><caption>Table 1: Calculated average lattice parameters a', and mixing energy, ΔH, for x = 0.33 using a 54-atom supercell over 20 different initial SQS seed values (i) for each value of p along with the average (avg) mixing energy value and standard deviation (σ). The number of nearest neighbor (NN) shells is given by the following convention: p NN where p is the number of nearest neighbor shells for the pair.</caption><thead><tr><th>#Shells<br>(p,0,0) NN</th><th>min{ΔH} meV/atom</th><th>avg ± σ</th><th>a' Å</th></tr></thead><tbody><tr><td>2</td><td>84.73</td><td>88.5±3.2</td><td>3.382</td></tr><tr><td>4</td><td>82.37</td><td>89.0±3.9</td><td>3.384</td></tr><tr><td>6</td><td>81.27</td><td>89.7±3.6</td><td>3.383</td></tr><tr><td>8</td><td>80.37</td><td>86.8±3.8</td><td>3.383</td></tr><tr><td>10</td><td>87.30</td><td>90.5±2.6</td><td>3.382</td></tr></tbody></table>

the standard deviation (σ) and minimum and maximum energy values. Overall, average mixing energy decreases as Ti content x increases. Further, the average mixing energy over 20 seed values is weakly impacted when three and four-body terms are included. The standard deviation is relatively large for some of the composition and simulation conditions. This could be the result of the impact of favorable and unfavorable pair types (see below); it could also be due to the lower mechanical stability of the bcc phase in the presence of the hcp prone elements.

For a given Ti content x, the average mixing energy (shown in a blue dotted color line) decreases only weakly when including t and q terms. A similar conclusion was also reached by Tian et al. for CoCrFeMnNi, Ta-W, and CoCrNi alloys. Using the Similar Atomic Environment (SAE) approach, these authors showed that including many-body clusters did not provide much improvement in the cross-validation error between SAE and SQS for equimolar composition, implying that pair NN terms for atom pairs are relatively enough to model a random structure [46].

As we will see later, for a Ti content close to x ≈ 0.5 (see Fig. 8), a phase transition from bcc to hcp is expected. Therefore, we generated another 20 SQS structures with orthorhombic and hcp symmetry, respectively. The orthorhombic structure generated with (8, 1, 0) exhibits a slightly lower formation energy, 7 meV/atom less than the other configurations shown in Fig. 2. It shows that for low-symmetry structures, the (t, q) might play
<table><caption>Table 1: Calculated average lattice parameters a', and mixing energy, ΔH, for x = 0.33 using a 54-atom supercell over 20 different initial SQS seed values (i) for each value of p along with the average (avg) mixing energy value and standard deviation (σ). The number of nearest neighbor (NN) shells is given by the following convention: p NN where p is the number of nearest neighbor shells for the pair.</caption><thead><tr><th rowspan="2">#Shells<br>(p,0,0) NN</th><th>min{ΔH}</th><th>avg ± σ</th><th>a'</th></tr><tr><th colspan="2">meV/atom</th><th>Å</th></tr></thead><tbody><tr><td>2</td><td>84.73</td><td>88.5±3.2</td><td>3.382</td></tr><tr><td>4</td><td>82.37</td><td>89.0±3.9</td><td>3.384</td></tr><tr><td>6</td><td>81.27</td><td>89.7±3.6</td><td>3.383</td></tr><tr><td>8</td><td>80.37</td><td>86.8±3.8</td><td>3.383</td></tr><tr><td>10</td><td>87.30</td><td>90.5±2.6</td><td>3.382</td></tr></tbody></table>

the standard deviation (σ) and minimum and maximum energy values. Overall, average mixing energy decreases as Ti content x increases. Further, the average mixing energy over 20 seed values is weakly impacted when three and four-body terms are included. The standard deviation is relatively large for some of the composition and simulation conditions. This could be the result of the impact of favorable and unfavorable pair types (see below); it could also be due to the lower mechanical stability of the bcc phase in the presence of the hcp prone elements.

For a given Ti content x, the average mixing energy (shown in a blue dotted color line) decreases only weakly when including t and q terms. A similar conclusion was also reached by Tian et al. for CoCrFeMnNi, Ta-W, and CoCrNi alloys. Using the Similar Atomic Environment (SAE) approach, these authors showed that including many-body clusters did not provide much improvement in the cross-validation error between SAE and SQS for equimolar composition, implying that pair NN terms for atom pairs are relatively enough to model a random structure [46].

As we will see later, for a Ti content close to x ≈ 0.5 (see Fig. 8), a phase transition from bcc to hcp is expected. Therefore, we generated another 20 SQS structures with orthorhombic and hcp symmetry, respectively. The orthorhombic structure generated with (8, 1, 0) exhibits a slightly lower formation energy, 7 meV/atom less than the other configurations shown in Fig. 2. It shows that for low-symmetry structures, the (t, q) might play
<figure><img src="image_3.png" /><figcaption>Figure 2: Statistical plot showing the average mixing energy, minimal and maximal values, along with the standard deviation of the calculated mixing energy for Ti$_x$(HfNbTaZr)$_{(1-x)/4}$ alloys for various MCSQS runs. The number of shells included in the MCSQS calculations is given by the following convention: (p, t, q), where p, t, and q are the pair, three-body, and four-body terms, respectively. Maximum values are represented by black-filled circles, and minimum values by green circles. The average mixing energy for each composition is shown in a blue dotted color line, and the standard deviation is shown in a gray rectangle. For each $x = 0.00$ and $0.20$, 650 SQS calculations were performed, respectively.</figcaption></figure>

a role in modeling the un-equimolar composition. For the hcp structure, we see a similar trend regarding the evolution of the mixing energy as a function of the choice of (p, t, q). The lowest formation energies for the three structures considered at $x = 0.5$ are very close to each other. This will be discussed in more detail later.

The addition of three and four body terms up to 1 NN shell appears to have only a second-order effect in the present systems. The inclusion of a large number of NN shells for pair interactions has an important impact on the ability of the SQS to yield lower energy configurations. This sensitivity study of over 700 different calculations demonstrates
<figure><img src="image_3.png" /><figcaption>Figure 2: Statistical plot showing the average mixing energy, minimal and maximal values, along with the standard deviation of the calculated mixing energy for Ti$_x$(HfNbTaZr)$_{(1-x)/4}$ alloys for various MCSQS runs. The number of shells included in the MCSQS calculations is given by the following convention: (p, t, q), where p, t, and q are the pair, three-body, and four-body terms, respectively. Maximum values are represented by black-filled circles, and minimum values by green circles. The average mixing energy for each composition is shown in a blue dotted color line, and the standard deviation is shown in a gray rectangle. For each $x = 0.00$ and $0.20$, 650 SQS calculations were performed, respectively.</figcaption></figure>

a role in modeling the un-equimolar composition. For the hcp structure, we see a similar trend regarding the evolution of the mixing energy as a function of the choice of (p, t, q). The lowest formation energies for the three structures considered at $x = 0.5$ are very close to each other. This will be discussed in more detail later.

The addition of three and four body terms up to 1 NN shell appears to have only a second-order effect in the present systems. The inclusion of a large number of NN shells for pair interactions has an important impact on the ability of the SQS to yield lower energy configurations. This sensitivity study of over 700 different calculations demonstrates
the importance of performing multiple MCSQS calculations with different seed values to
evaluate the statistical differences associated with the local atomic configurations.

## 3.2 Convergence of selected average properties

Next, we focus on the quaternary and quinary systems. The atomic configuration obtained
from the MC SQS technique constitutes only one of the realizations of the system, and the
alloy composition and simulation sizes considered here only represent a part of the many
atomic arrangements seen in a random solution. Thus, different MC SQS realizations
could lead to different configurations and thus provide a more complete picture of the
structure of the considered systems. We chose the two equimolar compositions, i.e., $x =
0.00$ and $0.20$, for which we generated a large set of 650 distinct SQS structures using
the optimized SQS parameters ($p = 8, t = 0, q = 0$ NN pairs). SQS calculations were
run to produce a random structure for maximizing the objective function using the above-
mentioned convergence criteria. After this, the relaxation of each structure was done using
DFT, and the mixing energy was evaluated for all different initial configurations. The data
can thus be analyzed to see how many SQS calculations are needed to provide a realistic
picture of the ideal random solid solution. For this, we will consider two aspects: i) the
average mixing energy, and ii) the SRO parameter (see below).

Fig. 3 displays the energy per atom $E_{at}^i$ and associated scatter obtained from each of our MC SQS runs "s", taken in the order they were created. The figure also shows the running average of the energy per atom: $\bar{E}_{at}(k) = \frac{\sum_{s=1}^{k} E_{at}^s}{k}$. Interestingly, this latter quantity converges rapidly to a roughly constant average of $\bar{E}_{at}(k)$ after about 20 MC SQS realizations for both equimolar alloys. This convergence of the average energy is obviously very empirical here, as sorting the raw energy data in a different order might change the running average evolution but not its final value. For the moment, we will simply retain the trends in the evolution of quantities and the number of MC SQS realizations required for average convergence.
the importance of performing multiple MCSQS calculations with different seed values to
evaluate the statistical differences associated with the local atomic configurations.

## 3.2 Convergence of selected average properties

Next, we focus on the quaternary and quinary systems. The atomic configuration obtained
from the MC SQS technique constitutes only one of the realizations of the system, and the
alloy composition and simulation sizes considered here only represent a part of the many
atomic arrangements seen in a random solution. Thus, different MC SQS realizations
could lead to different configurations and thus provide a more complete picture of the
structure of the considered systems. We chose the two equimolar compositions, i.e., $x =
0.00$ and $0.20$, for which we generated a large set of 650 distinct SQS structures using
the optimized SQS parameters ($p = 8, t = 0, q = 0$ NN pairs). SQS calculations were
run to produce a random structure for maximizing the objective function using the above-
mentioned convergence criteria. After this, the relaxation of each structure was done using
DFT, and the mixing energy was evaluated for all different initial configurations. The data
can thus be analyzed to see how many SQS calculations are needed to provide a realistic
picture of the ideal random solid solution. For this, we will consider two aspects: i) the
average mixing energy, and ii) the SRO parameter (see below).

Fig. 3 displays the energy per atom $E_{at}^i$ and associated scatter obtained from each of our MC SQS runs "s", taken in the order they were created. The figure also shows the running average of the energy per atom: $\bar{E}_{at}(k) = \frac{\sum_{s=1}^{k} E_{at}^s}{k}$. Interestingly, this latter quantity converges rapidly to a roughly constant average of $\bar{E}_{at}(k)$ after about 20 MC SQS realizations for both equimolar alloys. This convergence of the average energy is obviously very empirical here, as sorting the raw energy data in a different order might change the running average evolution but not its final value. For the moment, we will simply retain the trends in the evolution of quantities and the number of MC SQS realizations required for average convergence.
To quantify the SRO in the simulated atomic configurations, we use the SRO parameter $\alpha_{ij}$ described by Fontaine et al. [47],

$$ \alpha_{ij} = \frac{\frac{n_j}{m_A} - c_j}{\delta_{ij} - c_j} \qquad (4) $$

Here, $i$ is the reference atom. $n_j$ is the number of atoms of the non-reference type among the $c_j$ atoms in the $i$-th coordination shell. $m_A$ is the concentration of the non-reference atom. If $i = j$, $\delta_{ij} = 1$; otherwise, the Warren-Cowley SRO parameter is returned. The total SRO, $\alpha = \sum_{ij} |\alpha_{ij}|$, is calculated by summing all pairwise terms around the reference atom in the bcc structure for the first three nearest neighbor coordination shells (1-3 NN), respectively. We employed the tool developed in [48] and extended it to perform the analysis, including the first three nearest neighbor coordination shells (3 NN).

Fig. 4 shows the evolution of the total SRO parameter $\bar{\alpha}(k) = \sum_{s=1}^{k} \alpha^s/k/N_{ij}$ as a function of the considered SQS realizations included in the running average, where $N_{ij}$ is the number of bond types. First, it can be noted that the SRO values remain rather small. The total SRO parameter (black curve) converges rapidly to a somewhat constant value when considering more than 20 different MC SQS realizations. The final global SRO value is not perfectly equal to zero as expected for the ideal random solid solution, but is negligibly close with a value of -0.01. Interestingly, the results of the SRO analyses on the 1NN and 2NN shells taken individually showed similar trends, but running average convergence was achieved for a smaller number of MC SQS realizations.

The figure also shows the running averages of the SRO parameters $\alpha_{ij}(k)$ corresponding to each bond type. Quite expected, these quantities require a larger number of MC SQS (between 20 and 100) before converging toward the ideal random solid solution ($\alpha_{ij} \to 0$). Conversely, the empirical convergence analysis carried out here highlights the need to perform and consider several MC SQS realizations, at least for the alloys and simulation conditions used here. If a similar variety of configurations were sought in a single simulation
To quantify the SRO in the simulated atomic configurations, we use the SRO parameter $\alpha_{ij}$ described by Fontaine et al. [47],

$$ \alpha_{ij} = \frac{\frac{n_j}{m_A} - c_j}{\delta_{ij} - c_j} \qquad (4) $$

Here, $i$ is the reference atom. $n_j$ is the number of atoms of the non-reference type among the $c_j$ atoms in the $i$-th coordination shell. $m_A$ is the concentration of the non-reference atom. If $i = j$, $\delta_{ij} = 1$; otherwise, the Warren-Cowley SRO parameter is returned. The total SRO, $\alpha = \sum_{ij} |\alpha_{ij}|$, is calculated by summing all pairwise terms around the reference atom in the bcc structure for the first three nearest neighbor coordination shells (1-3 NN), respectively. We employed the tool developed in [48] and extended it to perform the analysis, including the first three nearest neighbor coordination shells (3 NN).

Fig. 4 shows the evolution of the total SRO parameter $\bar{\alpha}(k) = \sum_{s=1}^{k} \alpha^s/k/N_{ij}$ as a function of the considered SQS realizations included in the running average, where $N_{ij}$ is the number of bond types. First, it can be noted that the SRO values remain rather small. The total SRO parameter (black curve) converges rapidly to a somewhat constant value when considering more than 20 different MC SQS realizations. The final global SRO value is not perfectly equal to zero as expected for the ideal random solid solution, but is negligibly close with a value of -0.01. Interestingly, the results of the SRO analyses on the 1NN and 2NN shells taken individually showed similar trends, but running average convergence was achieved for a smaller number of MC SQS realizations.

The figure also shows the running averages of the SRO parameters $\alpha_{ij}(k)$ corresponding to each bond type. Quite expected, these quantities require a larger number of MC SQS (between 20 and 100) before converging toward the ideal random solid solution ($\alpha_{ij} \to 0$). Conversely, the empirical convergence analysis carried out here highlights the need to perform and consider several MC SQS realizations, at least for the alloys and simulation conditions used here. If a similar variety of configurations were sought in a single simulation
box, this would represent sizes one or two orders of magnitude larger than what has been considered in the present study, and such simulation sizes are currently mostly out of reach for DFT calculations.

## 3.3 Analysis of the chemical short-range order of the bcc phase

The results of the previous section raise the question of whether the statistical fluctuation could be understood among all SQS structures that exhibit a local CSRO, i.e., whether the favorable pairings between elements are preferred, leading to a configuration with the lowest mixing energy. MC and DFT calculations have been associated recently to investigate the short-range order (SRO) in two bcc alloys [20, 22], MoNbTaW and HfNbTiZr, respectively. The development of SRO was shown to decrease the mixing energy of the solid solution phase. Here, we show that the large number of different SQS simulations allows us to sample the energy of a disordered atomic configuration in a less computationally demanding way. We continue focusing on the two equimolar alloys and employing the 650 SQS structures obtained from the optimized SQS parameters ($p = 8, t = 0, q = 0$ NN pairs). In Fig. 2, for $x = 0.00$ and 0.20, the low and high mixing energies of the configurations are shown, along with the standard deviation and average values.

To assess the SRO, we sort the different atomic structures with respect to the simulated energy configuration. Following the lines of previous sections, we checked the significance of this sorting by analyzing the average of all 650 different configurations. Thus, the first set of data corresponds to the average of all the data and is shown as green curves in Figs. 5, 6 and 7. As desired, the SRO parameter is close to zero for all pair types. A small deviation on the SRO parameter of 0.01 is found on this average, and this gives an estimate of the precision we can expect from our analysis. The fact that the SRO parameter is very close to zero means that the average of all 650 SQS configurations is a good approximation of the ideal random solid solution. This neutral curve, representing an ideal solid solution,
box, this would represent sizes one or two orders of magnitude larger than what has been considered in the present study, and such simulation sizes are currently mostly out of reach for DFT calculations.

## 3.3 Analysis of the chemical short-range order of the bcc phase

The results of the previous section raise the question of whether the statistical fluctuation could be understood among all SQS structures that exhibit a local CSRO, i.e., whether the favorable pairings between elements are preferred, leading to a configuration with the lowest mixing energy. MC and DFT calculations have been associated recently to investigate the short-range order (SRO) in two bcc alloys [20, 22], MoNbTaW and HfNbTiZr, respectively. The development of SRO was shown to decrease the mixing energy of the solid solution phase. Here, we show that the large number of different SQS simulations allows us to sample the energy of a disordered atomic configuration in a less computationally demanding way. We continue focusing on the two equimolar alloys and employing the 650 SQS structures obtained from the optimized SQS parameters ($p = 8, t = 0, q = 0$ NN pairs). In Fig. 2, for $x = 0.00$ and 0.20, the low and high mixing energies of the configurations are shown, along with the standard deviation and average values.

To assess the SRO, we sort the different atomic structures with respect to the simulated energy configuration. Following the lines of previous sections, we checked the significance of this sorting by analyzing the average of all 650 different configurations. Thus, the first set of data corresponds to the average of all the data and is shown as green curves in Figs. 5, 6 and 7. As desired, the SRO parameter is close to zero for all pair types. A small deviation on the SRO parameter of 0.01 is found on this average, and this gives an estimate of the precision we can expect from our analysis. The fact that the SRO parameter is very close to zero means that the average of all 650 SQS configurations is a good approximation of the ideal random solid solution. This neutral curve, representing an ideal solid solution,
<figure><img src="image_4.png" /><figcaption>Figure 3: The calculated total energy per atom for every 650 distinct SQS structures is shown in blue, while the running average, $\bar{E}_{at}(k) = \sum_{s=1}^{k} E_{at}^{s}/k$ as a function of the growing number $k$ of SQS realizations is shown with black line for (a) $x = 0.00$ and in (b) $x = 0.20$.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: The calculated total energy per atom for every 650 distinct SQS structures is shown in blue, while the running average, $\bar{E}_{at}(k) = \sum_{s=1}^{k} E_{at}^{s}/k$ as a function of the growing number $k$ of SQS realizations is shown with black line for (a) $x = 0.00$ and in (b) $x = 0.20$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Running averages of the pairwise $\alpha_{ij}$ and Total SRO parameters $\bar{\alpha}(k) = \sum_{s=1}^{k} \alpha^s/k/N_{ij}$ as a function of the growing number $k$ of SQS realizations included in the average, for (a) $x = 0.00$ and in (b) $x = 0.20$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Running averages of the pairwise $\alpha_{ij}$ and Total SRO parameters $\bar{\alpha}(k) = \sum_{s=1}^{k} \alpha^s/k/N_{ij}$ as a function of the growing number $k$ of SQS realizations included in the average, for (a) $x = 0.00$ and in (b) $x = 0.20$.</figcaption></figure>
will act as a reference level line for the next steps. Two bins were then considered: the configurations with the lowest energy are grouped and called the low energy (LE) set, while the configurations with the highest energy are called the high energy (HE) set. To determine the SRO for the LE SQS set, we averaged the SRO values of each pair type for a group of "LE" cases with the lowest energy and used the same strategy for the HE SQS set. We varied the size of these bins from 6 to 65 configurations (1 – 10% of the total data), and a larger bin size smooths the results by reducing the obtained averaged SRO parameter values. However, the sign of the SRO parameters, indicating whether a bond type is more or less present than in the ideal random structure, remained largely unchanged when the bin size was varied.

## SRO in the quaternary HfNbTaZr system

Fig. 5 shows the SRO analysis in a radar plot for $x = 0.00$, and the interpretation is admittedly not straightforward. This analysis was performed for the first 1-3 NN shells, where most of the chemical ordering can be observed. An SRO parameter with a negative sign indicates a pair type that is more abundant in the atomic structure under consideration than in the ideal random configuration. When comparing two different SQS realizations, an increase in mixing energy can be explained in two ways: a lower representation of favorable pair types and/or an increase in the contribution of unfavorable atomic pairs. Comparing the results for the LE and HE sets, we see that the SRO data for the two sets agree qualitatively well, i.e., favorable bonds are found in greater numbers in the LE sets ($α_{ij} < 0$) and unfavorable bond types in smaller numbers ($a_{ij} > 0$). However, the trend is less clear when the SRO parameter is close to zero, indicating a neutral influence of the corresponding pair.

Turning now to the specific interactions in quaternary alloys, among the 10 different bond types, the following pairs appear to be favorable (in decreasing order): Nb-Ta and Hf-Zr; the same atom pairs (e.g., Hf-Hf) show a slight interaction ($α_{ij} ≤ 0$); finally,
will act as a reference level line for the next steps. Two bins were then considered: the configurations with the lowest energy are grouped and called the low energy (LE) set, while the configurations with the highest energy are called the high energy (HE) set. To determine the SRO for the LE SQS set, we averaged the SRO values of each pair type for a group of "LE" cases with the lowest energy and used the same strategy for the HE SQS set. We varied the size of these bins from 6 to 65 configurations (1 – 10% of the total data), and a larger bin size smooths the results by reducing the obtained averaged SRO parameter values. However, the sign of the SRO parameters, indicating whether a bond type is more or less present than in the ideal random structure, remained largely unchanged when the bin size was varied.

## SRO in the quaternary HfNbTaZr system

Fig. 5 shows the SRO analysis in a radar plot for $x = 0.00$, and the interpretation is admittedly not straightforward. This analysis was performed for the first 1-3 NN shells, where most of the chemical ordering can be observed. An SRO parameter with a negative sign indicates a pair type that is more abundant in the atomic structure under consideration than in the ideal random configuration. When comparing two different SQS realizations, an increase in mixing energy can be explained in two ways: a lower representation of favorable pair types and/or an increase in the contribution of unfavorable atomic pairs. Comparing the results for the LE and HE sets, we see that the SRO data for the two sets agree qualitatively well, i.e., favorable bonds are found in greater numbers in the LE sets ($\alpha_{ij} < 0$) and unfavorable bond types in smaller numbers ($a_{ij} > 0$). However, the trend is less clear when the SRO parameter is close to zero, indicating a neutral influence of the corresponding pair.

Turning now to the specific interactions in quaternary alloys, among the 10 different bond types, the following pairs appear to be favorable (in decreasing order): Nb-Ta and Hf-Zr; the same atom pairs (e.g., Hf-Hf) show a slight interaction ($\alpha_{ij} \le 0$); finally,
the last pair types seem to be unfavorable (in decreasing order): Hf-Nb, Ta-Zr, Nb-Zr, and finally Hf-Ta. Interestingly, the classification of a bond "ij" can also be determined in a completely different way, namely by a linear fit of the configurational energy as a function of the parameter $a_{ij}$ calculated for all 650 SQS realizations, where favorable (or unfavorable) bonds are associated with a negative slope (positive) of the linear fit. This further strengthens confidence in the proposed SRO analysis.

Previous results strikingly indicate a simple correlation between the favorable nature of the bond types and the crystallographic reference structure of the two elements involved in the pair type. This is shown in Fig. 6(a), where an SRO parameter is defined as the sum of SRO values corresponding to pairs of elements depending on their structure: bcc-bcc, hcp-hcp, or bcc-hcp elements (at 0 K). While this analysis is an oversimplification of the more complex picture in Fig. 5, elements with identical reference crystallographic structures typically seem to lead to favorable bonds associated with a decrease in mixing energy, while those with different bcc-hcp crystallographic structures are unfavorable in this system. It would be interesting to nuance this rather simple sorting of SRO analysis by the fact that the elements here have a similar radius and chemical properties precisely chosen to promote the formation of a stable HEA, namely a unique solid solution. If we had bcc elements other than the elements in the quaternary alloys, the findings would be much more complex. What is more, it can also be noticed that the binary phase diagrams of the elements involved in the favorable pairs lead to a unique solid solution, regardless of composition, while the elements of the unfavorable pairs give rise to a mixture of binary solid solutions. Fig. 5(b) shows the SRO parameter $a_{ij}$ when performing the analysis of an increasing number of shells from one to three NN shells for the LE SQS realizations. Interestingly, the first NN shell is associated with rather weak SRO parameters for all bond types, so the chemical SRO is mostly present in the 2nd and 3rd NN shells, which have a larger number of atoms.
the last pair types seem to be unfavorable (in decreasing order): Hf-Nb, Ta-Zr, Nb-Zr, and finally Hf-Ta. Interestingly, the classification of a bond "ij" can also be determined in a completely different way, namely by a linear fit of the configurational energy as a function of the parameter $a_{ij}$ calculated for all 650 SQS realizations, where favorable (or unfavorable) bonds are associated with a negative slope (positive) of the linear fit. This further strengthens confidence in the proposed SRO analysis.

Previous results strikingly indicate a simple correlation between the favorable nature of the bond types and the crystallographic reference structure of the two elements involved in the pair type. This is shown in Fig. 6(a), where an SRO parameter is defined as the sum of SRO values corresponding to pairs of elements depending on their structure: bcc-bcc, hcp-hcp, or bcc-hcp elements (at 0 K). While this analysis is an oversimplification of the more complex picture in Fig. 5, elements with identical reference crystallographic structures typically seem to lead to favorable bonds associated with a decrease in mixing energy, while those with different bcc-hcp crystallographic structures are unfavorable in this system. It would be interesting to nuance this rather simple sorting of SRO analysis by the fact that the elements here have a similar radius and chemical properties precisely chosen to promote the formation of a stable HEA, namely a unique solid solution. If we had bcc elements other than the elements in the quaternary alloys, the findings would be much more complex. What is more, it can also be noticed that the binary phase diagrams of the elements involved in the favorable pairs lead to a unique solid solution, regardless of composition, while the elements of the unfavorable pairs give rise to a mixture of binary solid solutions. Fig. 5(b) shows the SRO parameter $a_{ij}$ when performing the analysis of an increasing number of shells from one to three NN shells for the LE SQS realizations. Interestingly, the first NN shell is associated with rather weak SRO parameters for all bond types, so the chemical SRO is mostly present in the 2nd and 3rd NN shells, which have a larger number of atoms.
SRO in HfNbTaTiZr system and comparison between the two systems

Some of the previous observations made for the quaternary alloy can be repeated here,
although the details may differ. Fig. 7(a) shows that the average SRO parameters for
all 15 bond types are close to zero when considering all 650 SQS realizations, consistent
with the ideal random solid solution. The LE and HE data point sets again complement
each other well, with some exceptions that we will discuss later. The following bond types
appear to be uniquely favorable in this new system ($\alpha_{ij} < 0$ for the LE set): Hf-Zr, Hf-
Ti, and Nb-Ta; the same element pairs "i-i" are typically neutral ($\alpha_{ij} = 0$) or weakly
favored, while the following bonds appear unfavorable ($\alpha_{ij} > 0$): Hf-Ta, Nb-Ti, and Ta-Zr.
Once more, the binary phase diagram of Hf-Ti, the additional favorable pair compared to
the quaternary, shows one solid solution, for all alloy compositions. The results for some
bond pairs, such as Ta-Ti and Ti-Ti, are inconclusive with our methodology. More SQS
realizations are certainly needed to analyze this complex quinary system.

The simple sorting of bond types as a function of the crystallographic structure of each
element is shown in Fig. 6(b). As with the quaternary alloy, the sum of the SRO parameters
is favorable for pairs of elements with the same bcc-bcc or hcp-hcp structures, although
the SRO parameter for hcp-hcp has a lower value here. The SRO parameter for bcc-hcp
element bonds is unfavorable and has a much larger value than in the quaternary system.
Finally, the SRO parameters for all 15 pair types are analyzed as a function of NN shells
and shown in Fig. 7(b). In the quinary system, the bond types are polarized even in the 1st
NN shell, where the SRO value closely follows the sorting by the crystallographic nature
of the elements involved in the bond (the agreement is better than when considering 1st
+ 2nd NN shells). The SRO data for other shells typically evolves uniformly from this 1st
NN shell to the following 2nd and 3rd NN shells. Thus, the chemical order in the quinary
system appears to be more complex and more spread than in the quaternary systems, as
the 1st, 2nd, and 3rd NN shells are involved.

The total SROs for $x = 0.00$ and $0.20$ are given in Tab. 2 and are calculated from the
SRO in HfNbTaTiZr system and comparison between the two systems

Some of the previous observations made for the quaternary alloy can be repeated here,
although the details may differ. Fig. 7(a) shows that the average SRO parameters for
all 15 bond types are close to zero when considering all 650 SQS realizations, consistent
with the ideal random solid solution. The LE and HE data point sets again complement
each other well, with some exceptions that we will discuss later. The following bond types
appear to be uniquely favorable in this new system ($\alpha_{ij} < 0$ for the LE set): Hf-Zr, Hf-
Ti, and Nb-Ta; the same element pairs "i-i" are typically neutral ($\alpha_{ij} = 0$) or weakly
favored, while the following bonds appear unfavorable ($\alpha_{ij} > 0$): Hf-Ta, Nb-Ti, and Ta-Zr.
Once more, the binary phase diagram of Hf-Ti, the additional favorable pair compared to
the quaternary, shows one solid solution, for all alloy compositions. The results for some
bond pairs, such as Ta-Ti and Ti-Ti, are inconclusive with our methodology. More SQS
realizations are certainly needed to analyze this complex quinary system.

The simple sorting of bond types as a function of the crystallographic structure of each
element is shown in Fig. 6(b). As with the quaternary alloy, the sum of the SRO parameters
is favorable for pairs of elements with the same bcc-bcc or hcp-hcp structures, although
the SRO parameter for hcp-hcp has a lower value here. The SRO parameter for bcc-hcp
element bonds is unfavorable and has a much larger value than in the quaternary system.
Finally, the SRO parameters for all 15 pair types are analyzed as a function of NN shells
and shown in Fig. 7(b). In the quinary system, the bond types are polarized even in the 1st
NN shell, where the SRO value closely follows the sorting by the crystallographic nature
of the elements involved in the bond (the agreement is better than when considering 1st
+ 2nd NN shells). The SRO data for other shells typically evolves uniformly from this 1st
NN shell to the following 2nd and 3rd NN shells. Thus, the chemical order in the quinary
system appears to be more complex and more spread than in the quaternary systems, as
the 1st, 2nd, and 3rd NN shells are involved.

The total SROs for $x = 0.00$ and $0.20$ are given in Tab. 2 and are calculated from the
<figure><img src="image_6.png" /><figcaption>Figure 5: Radar plot showing the average values for $x = 0$: low energy (LE) set averaged over fifteen data points, and high energy (HE) set averaged over fifteen data points, for (a) 1st + 2nd + 3rd (1-3) NN data points and (b) low energy (LE) set SQS realizations for the 1st (1) NN, 1st + 2nd (1-2) NN, and 1st + 2nd + 3rd (1-3) NN shells, respectively. The LE and HE SQS sets are an average of the first and last fifteen data points, respectively. "AVG" here refers to the average over all SQS realizations. The red dashed line represents the ideal random solid solution.</figcaption></figure>

sum of the first two and the sum of the first three NN shells. In our analysis, we find that for $x = 0.00$, we have the lowest SRO value for the lowest energy configuration compared to the highest energy configuration. For $x = 0.20$, the SRO values are higher for the configurations with the lowest energy, indicating the existence of SRO. When comparing the quaternary and quinary alloys, the SRO in the first system is mostly localized on the 2nd NN shell, and only 2 bond types (Nb-Ta and Hf-Zr) are favorable, 4 bond types are unfavorable, and the same element bonds are mostly neutral. In the quinary system, 4 bond types are favorable, 6 types are unfavorable, and 5 types are slightly less than zero, while the chemical ordering covers at least all three NN shells. Since, in a finite-sized supercell, the bond type contributions are related through the imposed overall alloy composition, the quaternary system seems more constrained, and this could explain why the total SRO is smaller in that system. In the study performed by Yin *et al*. using the Monte Carlo approach (MC), the author observed a large Warren-Cowley SRO value (in the range of 1.75 – 2.00) for bcc MoNbTaW RHEA for the lowest energy configuration by considering
<figure><img src="image_6.png" /><figcaption>Figure 5: Radar plot showing the average values for $x = 0$: low energy (LE) set averaged over fifteen data points, and high energy (HE) set averaged over fifteen data points, for (a) 1st + 2nd + 3rd (1-3) NN data points and (b) low energy (LE) set SQS realizations for the 1st (1) NN, 1st + 2nd (1-2) NN, and 1st + 2nd + 3rd (1-3) NN shells, respectively. The LE and HE SQS sets are an average of the first and last fifteen data points, respectively. "AVG" here refers to the average over all SQS realizations. The red dashed line represents the ideal random solid solution.</figcaption></figure>

sum of the first two and the sum of the first three NN shells. In our analysis, we find that for $x = 0.00$, we have the lowest SRO value for the lowest energy configuration compared to the highest energy configuration. For $x = 0.20$, the SRO values are higher for the configurations with the lowest energy, indicating the existence of SRO. When comparing the quaternary and quinary alloys, the SRO in the first system is mostly localized on the 2nd NN shell, and only 2 bond types (Nb-Ta and Hf-Zr) are favorable, 4 bond types are unfavorable, and the same element bonds are mostly neutral. In the quinary system, 4 bond types are favorable, 6 types are unfavorable, and 5 types are slightly less than zero, while the chemical ordering covers at least all three NN shells. Since, in a finite-sized supercell, the bond type contributions are related through the imposed overall alloy composition, the quaternary system seems more constrained, and this could explain why the total SRO is smaller in that system. In the study performed by Yin *et al*. using the Monte Carlo approach (MC), the author observed a large Warren-Cowley SRO value (in the range of 1.75 – 2.00) for bcc MoNbTaW RHEA for the lowest energy configuration by considering
<figure><img src="image_7.png" /><figcaption>Figure 6: Radar plot showing the average, the low energy (LE) set and the high energy (HE) set averaged over fifteen data points for the case 1st + 2nd + 3rd (1-3) NN data points for (a) $x = 0.00$ and (b) $x = 0.20$, respectively. The LE and HE SQS sets are averaged over fifteen data points when the data are regrouped into "simple pair categories", "bcc-bcc", "bcc-hcp" and "hcp-hcp" based on the stable phase for the individual atoms of the pair. "AVG" here refers to the average over all SQS realizations. The red dashed line represents the ideal random solid solution.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Idem as Fig. 5, but for $x = 0.20$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Radar plot showing the average, the low energy (LE) set and the high energy (HE) set averaged over fifteen data points for the case 1st + 2nd + 3rd (1-3) NN data points for (a) $x = 0.00$ and (b) $x = 0.20$, respectively. The LE and HE SQS sets are averaged over fifteen data points when the data are regrouped into "simple pair categories", "bcc-bcc", "bcc-hcp" and "hcp-hcp" based on the stable phase for the individual atoms of the pair. "AVG" here refers to the average over all SQS realizations. The red dashed line represents the ideal random solid solution.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Idem as Fig. 5, but for $x = 0.20$.</figcaption></figure>
the first two NN shells [49]. However, it should be noted that the quaternary system they studied is arguably simpler as it consists of only bcc elements.

<table><caption>Table 2: Total SRO of the lowest SQS configuration (LSQS) and highest SQS configuration (HSQS) for the combined 1st+2nd nearest neighbors (represented as 2 NN) and 1st+2nd+3rd nearest neighbors (represented as 3 NN) shells, respectively. Total SRO is determined by using SRO = $\sum_{ij} |\alpha_{ij}|$. #atoms indicates the number of atoms in a SQS supercell.</caption><thead><tr><th rowspan="2">Ti<sub>x</sub>(HfNbTaZr)<sub>(1-x)/4</sub></th><th rowspan="2">#atoms</th><th colspan="2">LSQS</th><th colspan="2">HSQS</th></tr><tr><th>2NN</th><th>3NN</th><th>2NN</th><th>3NN</th></tr></thead><tbody><tr><td>0.0</td><td>128</td><td>0.57</td><td>0.51</td><td>0.88</td><td>0.81</td></tr><tr><td>0.2</td><td>125</td><td>1.55</td><td>1.61</td><td>1.46</td><td>1.48</td></tr></tbody></table>

Our results regarding the chemical SRO are mostly in qualitative agreement with the limited existing literature, and in particular with another very recent MC investigation on the same quinary alloys [18]. The MC relied on a cluster expansion model fitted on DFT calculations and conducted at different finite temperatures. While the values of the WC SRO parameter are quantitatively different (and larger) for the various pair types, the authors concluded that Hf-Ti and Zr-Ti showed the strongest chemical preference, followed by Hf-Nb, Zr-Nb, Nb-Ta, and finally Ta-Ta. The SRO was analyzed at small (1st NN only) and larger length scales (1st to 4th NN), and the SRO appeared weaker when analyzed at the latter scale. The strong affinity among Ti-Zr, Ti-Hf and Nb-Hf elements was also observed in the similar HfNbTiZr bcc alloy from both simulation [20, 21] and experimental investigations [14, 15]. In experiments, Hf-rich and Ti-rich clusters were identified from HR-HAADF-SEM and APT measurements [14, 15]. The enthalpy of mixing for Ti-Zr and Hf-Nb atomic pairs close to zero or negative was suggested as an explanation [18, 21]. The electronic density structure analysis in [21] highlighted the role of Ti-Zr SRO in particular, with a reduction of *d* electrons at the Fermi level, which was found to stabilize the bcc phase.

Similarly to our findings, the quaternary HfNbTaZr alloy seemed to be less prone to local chemical ordering in [18] when compared to other similar Ti-based HEA alloys,
the first two NN shells [49]. However, it should be noted that the quaternary system they studied is arguably simpler as it consists of only bcc elements.

<table><caption>Table 2: Total SRO of the lowest SQS configuration (LSQS) and highest SQS configuration (HSQS) for the combined 1st+2nd nearest neighbors (represented as 2 NN) and 1st+2nd+3rd nearest neighbors (represented as 3 NN) shells, respectively. Total SRO is determined by using SRO = $\sum_{ij} |\alpha_{ij}|$. #atoms indicates the number of atoms in a SQS supercell.</caption><thead><tr><th rowspan="2">Ti<sub>x</sub>(HfNbTaZr)<sub>(1-x)/4</sub></th><th rowspan="2">#atoms</th><th colspan="2">LSQS</th><th colspan="2">HSQS</th></tr><tr><th>2NN</th><th>3NN</th><th>2NN</th><th>3NN</th></tr></thead><tbody><tr><td>x</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0.0</td><td>128</td><td>0.57</td><td>0.51</td><td>0.88</td><td>0.81</td></tr><tr><td>0.2</td><td>125</td><td>1.55</td><td>1.61</td><td>1.46</td><td>1.48</td></tr></tbody></table>

Our results regarding the chemical SRO are mostly in qualitative agreement with the limited existing literature, and in particular with another very recent MC investigation on the same quinary alloys [18]. The MC relied on a cluster expansion model fitted on DFT calculations and conducted at different finite temperatures. While the values of the WC SRO parameter are quantitatively different (and larger) for the various pair types, the authors concluded that Hf-Ti and Zr-Ti showed the strongest chemical preference, followed by Hf-Nb, Zr-Nb, Nb-Ta, and finally Ta-Ta. The SRO was analyzed at small (1st NN only) and larger length scales (1st to 4th NN), and the SRO appeared weaker when analyzed at the latter scale. The strong affinity among Ti-Zr, Ti-Hf and Nb-Hf elements was also observed in the similar HfNbTiZr bcc alloy from both simulation [20, 21] and experimental investigations [14, 15]. In experiments, Hf-rich and Ti-rich clusters were identified from HR-HAADF-SEM and APT measurements [14, 15]. The enthalpy of mixing for Ti-Zr and Hf-Nb atomic pairs close to zero or negative was suggested as an explanation [18, 21]. The electronic density structure analysis in [21] highlighted the role of Ti-Zr SRO in particular, with a reduction of *d* electrons at the Fermi level, which was found to stabilize the bcc phase.

Similarly to our findings, the quaternary HfNbTaZr alloy seemed to be less prone to local chemical ordering in [18] when compared to other similar Ti-based HEA alloys,
suggesting a particular role of Ti atoms in the SRO development. Our SRO analysis in Fig. 5(a) matches very well with the results in [18] on the same alloy (for data from 1st to 4th NN).

## 3.4 Possible phases according to the Bo-Md diagram

To predict the phase stability of our alloys, we reprise the *d*-electrons theory initially proposed to predict the phase stability of Ti-based alloys. Morinaga *et al.* originally proposed the theory to classify Ti-based alloys into α or α + β or γ phase as a function of solely two parameters: Bo and Md [50]. Bo is the bond strength between the Ti and the alloying elements, and Md is the metal's *d*-orbital energy level, which corresponds to the element's metallic radius and electronegativity [50]. Based on this, Hadi *et al.* extended the single β-phase field to the domain from $Bo \le 2.84$ to $Bo \le 2.96$ [12].

For Ti-based alloys, Bo and Md parameters are given as composition averages: $Bo = \sum c_i Bo_i$ and $Md = \sum c_i Md_i$. The $c_i$ is the concentration of the alloying elements. Tab. 3 lists the composition averages for the alloys considered here. The reference Bo and Md values for each element in bcc Ti are taken from Hadi *et al.* [12]. We will assume that this analysis is still valid for our alloys, in which Ti is not the main constituent in some of our nuances [13].

<table><caption>Table 3: The Bo and Md for the $Ti_x(HfNbTaZr)_{(1-x)/4}$ alloys are calculated using the reference values of Bo and Md from Hadi *et al.* [12].</caption><thead><tr><th>x</th><th>0.0</th><th>0.12</th><th>0.20</th><th>0.33</th><th>0.40</th><th>0.50</th><th>0.63</th><th>0.70</th><th>0.81</th></tr></thead><tbody><tr><td>Bo</td><td>3.110</td><td>3.074</td><td>3.046</td><td>3.003</td><td>2.979</td><td>2.95</td><td>2.908</td><td>2.884</td><td>2.850</td></tr><tr><td>Md</td><td>2.715</td><td>2.685</td><td>2.662</td><td>2.626</td><td>2.607</td><td>2.582</td><td>2.548</td><td>2.528</td><td>2.499</td></tr></tbody></table>

In the Bo-Md diagram shown in Fig. 8, the domain of stability associated with the various possible phases of Ti based alloys are drawn along with some examples from the literature [12, 13]. The solid line in red is extrapolated for $Ti_x(HfNbTaZr)_{(1-x)/4}$ alloys with the composition studied in the present work to pure Ti system. For compositions
suggesting a particular role of Ti atoms in the SRO development. Our SRO analysis in Fig. 5(a) matches very well with the results in [18] on the same alloy (for data from 1st to 4th NN).

## 3.4 Possible phases according to the Bo-Md diagram

To predict the phase stability of our alloys, we reprise the *d*-electrons theory initially proposed to predict the phase stability of Ti-based alloys. Morinaga *et al.* originally proposed the theory to classify Ti-based alloys into α or α + β or γ phase as a function of solely two parameters: Bo and Md [50]. Bo is the bond strength between the Ti and the alloying elements, and Md is the metal's *d*-orbital energy level, which corresponds to the element's metallic radius and electronegativity [50]. Based on this, Hadi *et al.* extended the single β-phase field to the domain from $Bo \le 2.84$ to $Bo \le 2.96$ [12].

For Ti-based alloys, Bo and Md parameters are given as composition averages: $Bo = \sum c_i Bo_i$ and $Md = \sum c_i Md_i$. The $c_i$ is the concentration of the alloying elements. Tab. 3 lists the composition averages for the alloys considered here. The reference Bo and Md values for each element in bcc Ti are taken from Hadi *et al.* [12]. We will assume that this analysis is still valid for our alloys, in which Ti is not the main constituent in some of our nuances [13].

<table><caption>Table 3: The Bo and Md for the $Ti_x(HfNbTaZr)_{(1-x)/4}$ alloys are calculated using the reference values of Bo and Md from Hadi *et al.* [12].</caption><thead><tr><th>x</th><th>0.0</th><th>0.12</th><th>0.20</th><th>0.33</th><th>0.40</th><th>0.50</th><th>0.63</th><th>0.70</th><th>0.81</th></tr></thead><tbody><tr><td>Bo</td><td>3.110</td><td>3.074</td><td>3.046</td><td>3.003</td><td>2.979</td><td>2.95</td><td>2.908</td><td>2.884</td><td>2.850</td></tr><tr><td>Md</td><td>2.715</td><td>2.685</td><td>2.662</td><td>2.626</td><td>2.607</td><td>2.582</td><td>2.548</td><td>2.528</td><td>2.499</td></tr></tbody></table>

In the Bo-Md diagram shown in Fig. 8, the domain of stability associated with the various possible phases of Ti based alloys are drawn along with some examples from the literature [12, 13]. The solid line in red is extrapolated for $Ti_x(HfNbTaZr)_{(1-x)/4}$ alloys with the composition studied in the present work to pure Ti system. For compositions
well above the dashed line in the Bo-Md diagram, a bcc structure is expected, while for compositions well below the same line, structures tend to become hcp. For $x = 0.0$ (the quaternary alloy without Ti) and $x = 0.2$ (the equimolar case), experimental studies exist and provide clear reference data in the shape of unambiguous identification of a single HEA phase of bcc structure [30,51], corresponding to well-defined X-ray diffraction peaks for the bcc structure. The Bo-Md diagram predictions are thus in agreement with these studies on HEA alloys. From this, it becomes clear that relevant MCSQS and DFT simulations must be provided not only for the bcc phase but also for the hcp phase at larger Ti content $x$. The transition from a bcc structure to an hcp structure is expected at a Ti content between 0.5 and 0.6.

<figure><img src="image_9.png" /><figcaption>Figure 8: Extrapolation of the Bo-Md diagram for $Ti_x$ compositions. Considered alloys are represented by a solid red line. The dotted black line shows the separation between beta and alpha phases, the inset figure is taken from reference [12] and compiles data for various conventional Ti alloys for illustrative purposes.</figcaption></figure>
well above the dashed line in the Bo-Md diagram, a bcc structure is expected, while for compositions well below the same line, structures tend to become hcp. For $x = 0.0$ (the quaternary alloy without Ti) and $x = 0.2$ (the equimolar case), experimental studies exist and provide clear reference data in the shape of unambiguous identification of a single HEA phase of bcc structure [30,51], corresponding to well-defined X-ray diffraction peaks for the bcc structure. The Bo-Md diagram predictions are thus in agreement with these studies on HEA alloys. From this, it becomes clear that relevant MCSQS and DFT simulations must be provided not only for the bcc phase but also for the hcp phase at larger Ti content $x$. The transition from a bcc structure to an hcp structure is expected at a Ti content between 0.5 and 0.6.

<figure><img src="image_9.png" /><figcaption>Figure 8: Extrapolation of the Bo-Md diagram for $Ti_x$ compositions. Considered alloys are represented by a solid red line. The dotted black line shows the separation between beta and alpha phases, the inset figure is taken from reference [12] and compiles data for various conventional Ti alloys for illustrative purposes.</figcaption></figure>
## 3.5 Structure and stability of the bcc phase as a function of the Ti content

The rest of the article focuses mainly on the properties obtained for the *best* MCSQS obtained among 20+ realizations and associated with the lowest energy (LE) for each composition and structure (bcc, hcp and ortho). These LE configurations are assumed to be more relevant to the atomic arrangements that might be obtained in real alloys or after the final homogenization annealing step. Where significant, the difference between LE results and those obtained with the Highest Energy (HE) configurations or average (AVG) sets will also be shown. Since these sets correspond to rather different local arrangements and SRO, as shown in Figs. 5 and 7 for equimolar alloys, this comparison can highlight the impact of local atomic arrangements and the contribution of favorable/unfavorable bonds (this latter correlation will be mostly qualitative, as more statistics would be required for the non-equimolar alloys). Finally, the AVG set approximates the ideal random solid solution, as shown in Fig. 3.

In this section, we present DFT results on the effects of Ti content on the structure and relative stability of the possible phases.

Fig. 9 shows the evolution of the lattice parameter averaged over the three dimensions of the supercell ($a' = (a+b+c)/3$) for LE configurations as a function of the Ti content of the alloys in the bcc structure. More details concerning the structural properties of the phases possible for the alloys considered can be found in the supplementary material S1. The cubic $a'$ decreases monotonically from 3.436 Å to 3.252 Å with increasing Ti content $x$ from 0.00 to 1.00. Our DFT calculations reproduce well the experimental lattice parameter for equimolar alloys $x = 0.00$ and $x = 0.20$ [30, 51], giving confidence in our methodology based on several sets of MCSQS and DFT simulations. Fig. 9(b) displays the average lattice parameters ($a' = (a+b)/2$, $c/a'$) when the hcp structure is now considered for high Ti contents ($x \ge 0.5$). The parameter $c$ remains practically constant, while $a'$ shows a
## 3.5 Structure and stability of the bcc phase as a function of the Ti content

The rest of the article focuses mainly on the properties obtained for the *best* MCSQS obtained among 20+ realizations and associated with the lowest energy (LE) for each composition and structure (bcc, hcp and ortho). These LE configurations are assumed to be more relevant to the atomic arrangements that might be obtained in real alloys or after the final homogenization annealing step. Where significant, the difference between LE results and those obtained with the Highest Energy (HE) configurations or average (AVG) sets will also be shown. Since these sets correspond to rather different local arrangements and SRO, as shown in Figs. 5 and 7 for equimolar alloys, this comparison can highlight the impact of local atomic arrangements and the contribution of favorable/unfavorable bonds (this latter correlation will be mostly qualitative, as more statistics would be required for the non-equimolar alloys). Finally, the AVG set approximates the ideal random solid solution, as shown in Fig. 3.

In this section, we present DFT results on the effects of Ti content on the structure and relative stability of the possible phases.

Fig. 9 shows the evolution of the lattice parameter averaged over the three dimensions of the supercell ($a' = (a+b+c)/3$) for LE configurations as a function of the Ti content of the alloys in the bcc structure. More details concerning the structural properties of the phases possible for the alloys considered can be found in the supplementary material S1. The cubic $a'$ decreases monotonically from 3.436 Å to 3.252 Å with increasing Ti content $x$ from 0.00 to 1.00. Our DFT calculations reproduce well the experimental lattice parameter for equimolar alloys $x = 0.00$ and $x = 0.20$ [30, 51], giving confidence in our methodology based on several sets of MCSQS and DFT simulations. Fig. 9(b) displays the average lattice parameters ($a' = (a+b)/2$, $c/a'$) when the hcp structure is now considered for high Ti contents ($x \ge 0.5$). The parameter $c$ remains practically constant, while $a'$ shows a
surprising non-monotonic evolution with Ti content. Finally, these parameters gradually converge towards those of pure Ti hcp for the highest *x*, with *a* = 2.94 Å and *c* = 4.64 Å.

Interestingly, the detail of the local atomic arrangement of the various elements in the supercell has little impact on the average lattice parameters. Anisotropic lattice parameters or geometries are slightly influenced by the composition of the atomic elements in the considered lattice directions (cf. Supp. Mat. and Tab. S1). In contrast, the average lattice parameters are virtually identical for the LE, AVG and HE configurations, suggesting average convergence of the lattice distortions typical of HEA systems, and thus these average parameters depend solely on alloy composition.

<figure><img src="image_10.png" /></figure>

Figure 9: Calculated average lattice parameters $a'$ as a function of Ti content $x$ are reported for the conventional unit cell. (a) displays the average lattice parameters for the bcc phase for the LE sets, while (b) presents the average lattice parameters for the hcp phase for the LE set. The red dots represent the experimental lattice parameters for the given composition [30, 51]. The tabulated data can be found in the supplementary section (see Tab. S1).

Fig. 10(a) shows the computed mixing energies of the possible phases of Ti content $x$ for the LE sets. Mixing energies are calculated using Eq. 2 by taking the energy of the most stable phase of each alloying element as a reference (namely, bcc for Nb, hcp for Ti, and so on). Interestingly, the mixing energies of all phases of the alloy compositions are slightly positive for all Ti content $x$. The corresponding values are, however small, below 100 meV/atom. This will be discussed further below.

The mixing enthalpy of the bcc phase decreases monotonically with the Ti content $x$. For $x = 0.20$, which corresponds to the equimolar case, the mixing energy of 84.4
surprising non-monotonic evolution with Ti content. Finally, these parameters gradually converge towards those of pure Ti hcp for the highest *x*, with *a* = 2.94 Å and *c* = 4.64 Å.

Interestingly, the detail of the local atomic arrangement of the various elements in the supercell has little impact on the average lattice parameters. Anisotropic lattice parameters or geometries are slightly influenced by the composition of the atomic elements in the considered lattice directions (cf. Supp. Mat. and Tab. S1). In contrast, the average lattice parameters are virtually identical for the LE, AVG and HE configurations, suggesting average convergence of the lattice distortions typical of HEA systems, and thus these average parameters depend solely on alloy composition.

<figure><img src="image_10.png" /></figure>

Figure 9: Calculated average lattice parameters $a'$ as a function of Ti content $x$ are reported for the conventional unit cell. (a) displays the average lattice parameters for the bcc phase for the LE sets, while (b) presents the average lattice parameters for the hcp phase for the LE set. The red dots represent the experimental lattice parameters for the given composition [30, 51]. The tabulated data can be found in the supplementary section (see Tab. S1).

Fig. 10(a) shows the computed mixing energies of the possible phases of Ti content $x$ for the LE sets. Mixing energies are calculated using Eq. 2 by taking the energy of the most stable phase of each alloying element as a reference (namely, bcc for Nb, hcp for Ti, and so on). Interestingly, the mixing energies of all phases of the alloy compositions are slightly positive for all Ti content $x$. The corresponding values are, however small, below 100 meV/atom. This will be discussed further below.

The mixing enthalpy of the bcc phase decreases monotonically with the Ti content $x$. For $x = 0.20$, which corresponds to the equimolar case, the mixing energy of 84.4
meV/atom determined in the present study is similar to the DFT calculations of Gao *et al.* [52], where a value of 86.5 meV/atom was reported. Fig. 10(a) also displays the mixing energies of pure Ti in the bcc and hcp phases. The mixing energy of Ti$_x$(HfNbTaZr)$_{(1-x)}$/4 alloys taken in the hcp structure smoothly decreases towards that of pure Ti in the hcp structure.

<figure><img src="image_11.png" /><figcaption>Figure 10: Calculated mixing energy (ΔH) and Gibbs free energy (ΔG) for Ti$_x$(HfNbTaZr)$_{(1-x)}$/4 alloys using the reference state as stable bcc and hcp phases for bulk systems at 0 K. In (a), mixing enthalpies and free energies at 300 K and 1000 K are given for the LE set only as a function of Ti content. Fig. (b) shows the impact of the local atomic arrangement onto the mixing enthalpy, through results corresponding to LE, AVG and HE sets. The inset highlights the extended overlap of mixing energies of the three possible phases.</figcaption></figure>

If we compare the relative stability of the possible phases according to the LE sets only for now at $x > 0.5$ in Fig. 10(a), mixing enthalpy for the hcp phase becomes lower than that of the bcc phase and decreases faster as Ti content $x$ is increased. This suggests that alloys with a large Ti content may exhibit an hcp structure or dual bcc and hcp phases. The mixing energy difference between the two structures is, however, very low (i.e., 2 meV/atom at the intersection, i.e., $x = 0.50$ and 5 meV/atom for $x = 0.63$). The mixing energy of the orthorhombic phase at $x = 0.50$ is slightly smaller than that of the bcc and hcp phases and could constitute a third possible structure. These phases could well be found in the synthesized unstrained samples or reveal themselves by phase transformation upon straining, as in [13].
meV/atom determined in the present study is similar to the DFT calculations of Gao *et al.* [52], where a value of 86.5 meV/atom was reported. Fig. 10(a) also displays the mixing energies of pure Ti in the bcc and hcp phases. The mixing energy of Ti$_x$(HfNbTaZr)$_{(1-x)}$/4 alloys taken in the hcp structure smoothly decreases towards that of pure Ti in the hcp structure.

<figure><img src="image_11.png" /><figcaption>Figure 10: Calculated mixing energy (ΔH) and Gibbs free energy (ΔG) for Ti$_x$(HfNbTaZr)$_{(1-x)}$/4 alloys using the reference state as stable bcc and hcp phases for bulk systems at 0 K. In (a), mixing enthalpies and free energies at 300 K and 1000 K are given for the LE set only as a function of Ti content. Fig. (b) shows the impact of the local atomic arrangement onto the mixing enthalpy, through results corresponding to LE, AVG and HE sets. The inset highlights the extended overlap of mixing energies of the three possible phases.</figcaption></figure>

If we compare the relative stability of the possible phases according to the LE sets only for now at $x > 0.5$ in Fig. 10(a), mixing enthalpy for the hcp phase becomes lower than that of the bcc phase and decreases faster as Ti content $x$ is increased. This suggests that alloys with a large Ti content may exhibit an hcp structure or dual bcc and hcp phases. The mixing energy difference between the two structures is, however, very low (i.e., 2 meV/atom at the intersection, i.e., $x = 0.50$ and 5 meV/atom for $x = 0.63$). The mixing energy of the orthorhombic phase at $x = 0.50$ is slightly smaller than that of the bcc and hcp phases and could constitute a third possible structure. These phases could well be found in the synthesized unstrained samples or reveal themselves by phase transformation upon straining, as in [13].
This picture of relative phase stability becomes more nuanced when fluctuations in the local atomic arrangements are taken into account. Fig. 10(b) also shows the enthalpies of mixing for the LE, AVG and HE sets for the different alloy compositions. Assuming that for a given overall alloy composition, the atomic microstructure of a large system consists of an assortment of various atomic configurations in the domain bounded by the LE and HE sets, the intersection between the bcc and hcp phases is not unique and now extends to a compositional range between 0.4 and 0.6 (see figure inset). This means that different phases can coexist in correlation with local chemical fluctuations in the alloy. Along these lines, it can be assumed that an as-cast alloy may contain a greater variety of atomic arrangements, including HE types, while a relaxed microstructure obtained after annealing may reinforce the local chemical order with favorable bonds associated with lower energies. The inset in the figure shows that the Ti content at which the bcc/hcp curves intersect shifts to the left when HE and LE sets are considered. It is therefore possible for a given sample to undergo a phase transition between the as-cast and annealed states if the alloy composition lies between these two points of intersection.

Such a structural evolution has been observed in a related Ti₃₅Zr₂₇.₅Hf₂₇.₅Nb₅Ta₅ alloy, where the as-cast alloy exhibited a metastable orthorombic structure, which transforms into an hcp structure as soon as residual stresses are removed during sample preparation [53]. The same alloy exhibits a bcc structure after homogenization processing [13]. These effects may be stronger for alloys that are more prone to chemical ordering, as in [23].

In the context of high-entropy alloys, the mixing entropy term has been thought to stabilize the HEA phase by lowering the Gibbs free energy. The consensus is now more nuanced, as it has been shown that high mixing entropy is neither a sufficient nor a necessary condition for the formation of a single-phase solid solution [54–56]. Fig. 10(a) shows that the Gibbs free energy calculated at RT accounting for the mixing entropy term as given in Eq. 3 is still not enough to reach negative values.

The fact that the mixing energy remains positive, albeit slightly, while a stable solid
This picture of relative phase stability becomes more nuanced when fluctuations in the local atomic arrangements are taken into account. Fig. 10(b) also shows the enthalpies of mixing for the LE, AVG and HE sets for the different alloy compositions. Assuming that for a given overall alloy composition, the atomic microstructure of a large system consists of an assortment of various atomic configurations in the domain bounded by the LE and HE sets, the intersection between the bcc and hcp phases is not unique and now extends to a compositional range between 0.4 and 0.6 (see figure inset). This means that different phases can coexist in correlation with local chemical fluctuations in the alloy. Along these lines, it can be assumed that an as-cast alloy may contain a greater variety of atomic arrangements, including HE types, while a relaxed microstructure obtained after annealing may reinforce the local chemical order with favorable bonds associated with lower energies. The inset in the figure shows that the Ti content at which the bcc/hcp curves intersect shifts to the left when HE and LE sets are considered. It is therefore possible for a given sample to undergo a phase transition between the as-cast and annealed states if the alloy composition lies between these two points of intersection.

Such a structural evolution has been observed in a related Ti₃₅Zr₂₇.₅Hf₂₇.₅Nb₅Ta₅ alloy, where the as-cast alloy exhibited a metastable orthorombic structure, which transforms into an hcp structure as soon as residual stresses are removed during sample preparation [53]. The same alloy exhibits a bcc structure after homogenization processing [13]. These effects may be stronger for alloys that are more prone to chemical ordering, as in [23].

In the context of high-entropy alloys, the mixing entropy term has been thought to stabilize the HEA phase by lowering the Gibbs free energy. The consensus is now more nuanced, as it has been shown that high mixing entropy is neither a sufficient nor a necessary condition for the formation of a single-phase solid solution [54–56]. Fig. 10(a) shows that the Gibbs free energy calculated at RT accounting for the mixing entropy term as given in Eq. 3 is still not enough to reach negative values.

The fact that the mixing energy remains positive, albeit slightly, while a stable solid
solution phase is observed for at least some of the alloy compositions in the experiments can be explained in several ways. First, the solid solution mixing energy is strongly impacted by the structure reference chosen for the elements, most of which are stable in hcp at 0 K while the others are stable in bcc structure. These HEA are, however, typically elaborated at very high temperatures > 2000 K, and all these elements exhibit a stable bcc structure at high temperature. The calculated Gibbs free energy nevertheless becomes negative at a temperature of about 1000 K, well below the elaboration temperature, for all considered alloy compositions (see Fig. 10(b)). Ti, Zr, and Hf are still stable in hcp structure at that temperature [57]. If we define the bcc structure as the reference one for all these elements, the mixing enthalpy will become clearly negative for all alloy compositions.

Finally, when considering high temperatures, vibrational entropy $S_{vib}$ could also play a stabilizing role in the solid solution phase. Interestingly, of all the considered compositions, the quinary equimolar alloy ($x = 0.20$) exhibits the lowest mixing energy, mostly due to a greater contribution from entropy.

## 3.6 Elastic constants for the bcc phase

In this final section, we evaluate the elastic constants of the bcc phase as a function of the Ti content. The energy-strain approach is a commonly used method to determine the elastic constants (ECs) at 0 K. In this study, the strain tensor was considered based on the symmetry of the relaxed supercell, as described in [58]. To accurately determine the unique elastic constants, a strain range of +/-2% was applied to the equilibrium structure, resulting in 11 deformed structures. These structures were then relaxed to obtain the atomic positions. To focus specifically on the bcc elastic properties, calculations were limited to alloys with a Ti content $x < 0.5$, where the bcc phase is expected to be the most stable.

Fig. 11(a) and Tab. 4 show the evolution of the ECs for the LE configurations measured for our alloy Ti compositions. $C_{11}$ decreases slowly at first and then rapidly with the
solution phase is observed for at least some of the alloy compositions in the experiments can
be explained in several ways. First, the solid solution mixing energy is strongly impacted
by the structure reference chosen for the elements, most of which are stable in hcp at 0 K
while the others are stable in bcc structure. These HEA are, however, typically elaborated
at very high temperatures > 2000 K, and all these elements exhibit a stable bcc structure
at high temperature. The calculated Gibbs free energy nevertheless becomes negative at
a temperature of about 1000 K, well below the elaboration temperature, for all considered
alloy compositions (see Fig. 10(b)). Ti, Zr, and Hf are still stable in hcp structure at that
temperature [57]. If we define the bcc structure as the reference one for all these elements,
the mixing enthalpy will become clearly negative for all alloy compositions.

Finally, when considering high temperatures, vibrational entropy $S_{vib}$ could also play a stabilizing role in the solid solution phase. Interestingly, of all the considered compositions, the quinary equimolar alloy ($x = 0.20$) exhibits the lowest mixing energy, mostly due to a greater contribution from entropy.

## 3.6 Elastic constants for the bcc phase

In this final section, we evaluate the elastic constants of the bcc phase as a function of the Ti content. The energy-strain approach is a commonly used method to determine the elastic constants (ECs) at 0 K. In this study, the strain tensor was considered based on the symmetry of the relaxed supercell, as described in [58]. To accurately determine the unique elastic constants, a strain range of +/-2% was applied to the equilibrium structure, resulting in 11 deformed structures. These structures were then relaxed to obtain the atomic positions. To focus specifically on the bcc elastic properties, calculations were limited to alloys with a Ti content $x < 0.5$, where the bcc phase is expected to be the most stable.

Fig. 11(a) and Tab. 4 show the evolution of the ECs for the LE configurations measured for our alloy Ti compositions. $C_{11}$ decreases slowly at first and then rapidly with the
increase of Ti content, while $C_{12}$ decreases at first and $C_{44}$ remains mostly unchanged. Born stability conditions on elastic constants for cubic systems also provide an alternative way of assessing the structural stability of alloys. We therefore checked that the elastic constants for $x$ from 0 to 0.4 satisfy the following conditions, i.e., $C_{11} > C_{12}$, $C_{11}+2C_{12} > 0$, and $C_{44} > 0$. These conditions are, however, not fulfilled by our calculated ECs for large Ti content $x > 0.7$ (when extrapolating the simulation results past $x > 0.4$). For these large Ti contents, some of the DFT calculations did not converge as the atomic configurations spontaneously transformed into an hcp lattice, confirming this instability and suggesting a correlation with the local atomic arrangement (see below). This loss of mechanical stability of the bcc phase is thus in agreement with the Bo-Md diagram prediction and our results regarding the mixing energies of the bcc and hcp phases in the previous section. The bcc phase of pure Ti becomes mechanically stable at high temperature, as shown by the experimental evaluation of ECs in [59].

A quantitative comparison can be made with ECs calculated in other DFT investigations [60, 61] or obtained from experiments [4, 59] for the two equimolar alloys and pure bcc Ti, where an overall good agreement is found. For the quaternary alloy, our results agree well with the DFT results from [61] with the exception of $C_{12}$, which is larger in [61] but still within an acceptable 15% difference margin typical of EC evaluations by DFT. For the quinary alloy, a nice agreement is also found with experimental data from [4]; the most significant discrepancy is found for $C_{44}$ (see below). A comparison can also be made with another simulation concerning equimolar quinary alloys carried out using the Coherent Potential Approximation (CPA) [60]. It yielded $C_{11} = 160.2$ GPa, $C_{12} = 124.4$ GPa, and $C_{44} = 62.4$ GPa. While the first two constants correspond well to experimental constants, the last value seems particularly high compared to experimental and our own data. The coherent potential approximation assumes a perfect mixing of atoms within the alloy, neglecting any local atomic ordering or clustering that may occur.

The isotropic bulk modulus $B$, shear modulus $G$, and Young's modulus $E$ can be
increase of Ti content, while $C_{12}$ decreases at first and $C_{44}$ remains mostly unchanged.

Born stability conditions on elastic constants for cubic systems also provide an alternative way of assessing the structural stability of alloys. We therefore checked that the elastic constants for $x$ from 0 to 0.4 satisfy the following conditions, i.e., $C_{11} > C_{12}$, $C_{11}+2C_{12} > 0$, and $C_{44} > 0$. These conditions are, however, not fulfilled by our calculated ECs for large Ti content $x > 0.7$ (when extrapolating the simulation results past $x > 0.4$). For these large Ti contents, some of the DFT calculations did not converge as the atomic configurations spontaneously transformed into an hcp lattice, confirming this instability and suggesting a correlation with the local atomic arrangement (see below). This loss of mechanical stability of the bcc phase is thus in agreement with the Bo-Md diagram prediction and our results regarding the mixing energies of the bcc and hcp phases in the previous section. The bcc phase of pure Ti becomes mechanically stable at high temperature, as shown by the experimental evaluation of ECs in [59].

A quantitative comparison can be made with ECs calculated in other DFT investigations [60, 61] or obtained from experiments [4, 59] for the two equimolar alloys and pure bcc Ti, where an overall good agreement is found. For the quaternary alloy, our results agree well with the DFT results from [61] with the exception of $C_{12}$, which is larger in [61] but still within an acceptable 15% difference margin typical of EC evaluations by DFT. For the quinary alloy, a nice agreement is also found with experimental data from [4]; the most significant discrepancy is found for $C_{44}$ (see below). A comparison can also be made with another simulation concerning equimolar quinary alloys carried out using the Coherent Potential Approximation (CPA) [60]. It yielded $C_{11} = 160.2$ GPa, $C_{12} = 124.4$ GPa, and $C_{44} = 62.4$ GPa. While the first two constants correspond well to experimental constants, the last value seems particularly high compared to experimental and our own data. The coherent potential approximation assumes a perfect mixing of atoms within the alloy, neglecting any local atomic ordering or clustering that may occur.

The isotropic bulk modulus $B$, shear modulus $G$, and Young's modulus $E$ can be
<figure><img src="image_12.png" /><figcaption>Figure 11: Evolution of the elastic moduli for $Ti_x(HfNbTaZr)_{(1-x)/4}$ alloys in the bcc structure as a function of the Ti content $x$. (a) the elastic constants $C_{ij}$ (GPa) and (b) the isotropic averages. Comparison with data from the literature is shown, corresponding to DFT calculations for the quaternary alloy [61] (half solid symbols) and experimental data for the quinary alloy [4] (hollow symbols) along with pure bcc Ti at high temperature $T = 1273$ K [59] for illustration. Corresponding data can be found in supplementary Tab. S2.</figcaption></figure>

evaluated from previous elastic constants for the various Ti content $x$. The shear modulus was obtained from the average of the Voigt and Reuss expressions, as suggested by Hill [62]. Fig. 11(b) shows the evolution of the isotropic elastic constants as a function of the Ti content. $B$ decreases monotonously when increasing the Ti content, while $G$ and $E$ are roughly constant for $x < 0.5$. Finally, Pugh's ratio ($B/G$) indicates the brittle or ductile nature of the material's behavior [63]. Materials that present a large $B$ and a low $G$ typically tend to exhibit ductile behavior. Mechanical testing conducted on the $x = 0.2$ [4] already demonstrated the ductile behavior of this composition. Based on our calculations, Pugh's ratio $B/G$ is > 1.75 for all Ti content $x$ from 0 to 0.4, suggesting a ductile behavior for this entire range of compositions.

Finally, we end this section by assessing whether the local atomic arrangement affects the elastic response of the considered alloys. To this aim, we evaluated ECs for the LE, one AVG and the HE configurations while focusing on equimolar compositions. The corresponding results are also presented in Tab. 4. For a given alloy composition, we expected initially the ECs to be only very weakly influenced by the local atomic arrangement, since the ECs are already an average measure of the collective response of the atoms contained
<figure><img src="image_12.png" /><figcaption>Figure 11: Evolution of the elastic moduli for Ti<sub>x</sub>(HfNbTaZr)<sub>(1-x)</sub> alloys in the bcc structure as a function of the Ti content x. (a) the elastic constants C<sub>ij</sub> (GPa) and (b) the isotropic averages. Comparison with data from the literature is shown, corresponding to DFT calculations for the quaternary alloy [61] (half solid symbols) and experimental data for the quinary alloy [4] (hollow symbols) along with pure bcc Ti at high temperature T = 1273 K [59] for illustration. Corresponding data can be found in supplementary Tab. S2.</figcaption></figure>

evaluated from previous elastic constants for the various Ti content x. The shear modulus was obtained from the average of the Voigt and Reuss expressions, as suggested by Hill [62]. Fig. 11(b) shows the evolution of the isotropic elastic constants as a function of the Ti content. B decreases monotonously when increasing the Ti content, while G and E are roughly constant for x < 0.5. Finally, Pugh's ratio (B/G) indicates the brittle or ductile nature of the material's behavior [63]. Materials that present a large B and a low G typically tend to exhibit ductile behavior. Mechanical testing conducted on the x = 0.2 [4] already demonstrated the ductile behavior of this composition. Based on our calculations, Pugh's ratio B/G is > 1.75 for all Ti content x from 0 to 0.4, suggesting a ductile behavior for this entire range of compositions.

Finally, we end this section by assessing whether the local atomic arrangement affects the elastic response of the considered alloys. To this aim, we evaluated ECs for the LE, one AVG and the HE configurations while focusing on equimolar compositions. The corresponding results are also presented in Tab. 4. For a given alloy composition, we expected initially the ECs to be only very weakly influenced by the local atomic arrangement, since the ECs are already an average measure of the collective response of the atoms contained
<table><caption>Table 4: Calculated values for elastic constants $C_{ij}$ (GPa) with Ti content $x$ for the investigated alloys in the bcc structure. Data in parentheses are extracted from the literature and correspond to DFT calculations for $x = 0$ [61], experimental data for $x = 0.2$ [4] and pure bcc Ti at T = 1273 K [59]. Data for the LE, AVG and HE configurations are also provided for equimolar alloys to assess the impact of the local atomic arrangement on the elastic response of these systems. Zener Anisotropy, $A_z$, coefficient is also given.</caption><thead><tr><th rowspan="2">x</th><th colspan="3">0.0</th><th colspan="3">0.11</th><th colspan="3">0.20</th><th rowspan="2">0.33</th><th rowspan="2">0.40</th><th rowspan="2">1.0 ( [59])</th></tr><tr><th>LE</th><th>AVG</th><th>HE ( [61])</th><th>LE</th><th>AVG</th><th>HE ( [4])</th><th>LE</th><th>AVG</th><th>HE ( [4])</th></tr></thead><tbody><tr><td>C<sub>11</sub></td><td>162</td><td>161</td><td>158 (169)</td><td>155</td><td>156</td><td>151</td><td>145 (172±6)</td><td>154</td><td>147</td><td>70 (98)</td><td></td><td></td></tr><tr><td>C<sub>12</sub></td><td>109</td><td>107</td><td>109 (127)</td><td>108</td><td>105</td><td>103</td><td>105 (108±1.5)</td><td>100</td><td>100</td><td>123 (83)</td><td></td><td></td></tr><tr><td>C<sub>44</sub></td><td>41</td><td>39</td><td>36 (45)</td><td>46</td><td>39</td><td>38</td><td>32 (28±1.5)</td><td>45</td><td>41</td><td>39 (38)</td><td></td><td></td></tr><tr><td>A<sub>z</sub></td><td>1.5</td><td>1.4</td><td>1.5 (2.1)</td><td>1.9</td><td>1.5</td><td>1.6</td><td>1.6 (0.87)</td><td>1.7</td><td>1.9</td><td>- (5)</td><td></td><td></td></tr></tbody></table>

in the supercell subjected to deformation. Furthermore, the prediction of ECs from DFT simulations is known to be limited in accuracy and is strongly influenced by the choice of XC function. That said, the results in Tab. 4 show a consistent and slight decrease of ECs ($C_{44}$ in particular) when measured for the LE set compared to the AVG and HE sets. These differences between the EC calculated for the LE and HE may be explained in two ways. i) The deeper energy wells corresponding to LE configurations may be circumvented by steeper energy variations, resulting in stiffer ECs. ii) HE configurations contain relatively more of unfavorable bonds, in particular bcc-hcp ones (cf. Figs. 6 and 7), possibly making the system less mechanically stable and lowering the ECs. Interestingly, the ECs ($C_{44}$) obtained for the HE configurations are in better agreement with experimental determinations for the quinary alloy [4], which were evaluated on the as-cast state that may well contain more high-energy atomic arrangements than in the final annealed state.

# 4 Conclusions

In this work, we have studied the structure, phase stability, and elastic behavior of $Ti_x(HfNbTaZr)_{(1-x)/4}$ alloys as a function of the Ti content $x$. The atomistic structure, corresponding to the disordered structure of the HEA solid solution, was approximated using a large number of different MCSQS realizations, while the atomic relaxation and
<table><caption>Table 4: Calculated values for elastic constants $C_{ij}$ (GPa) with Ti content $x$ for the investigated alloys in the bcc structure. Data in parentheses are extracted from the literature and correspond to DFT calculations for $x = 0$ [61], experimental data for $x = 0.2$ [4] and pure bcc Ti at T = 1273 K [59]. Data for the LE, AVG and HE configurations are also provided for equimolar alloys to assess the impact of the local atomic arrangement on the elastic response of these systems. Zener Anisotropy, $A_z$, coefficient is also given.</caption><thead><tr><th rowspan="2">x</th><th colspan="3">0.0</th><th colspan="3">0.11</th><th colspan="3">0.20</th><th rowspan="2">0.33</th><th rowspan="2">0.40</th><th rowspan="2">1.0 ( [59])</th></tr><tr><th>LE</th><th>AVG</th><th>HE ( [61])</th><th>LE</th><th>AVG</th><th>HE ( [4])</th><th>LE</th><th>AVG</th><th>HE ( [4])</th></tr></thead><tbody><tr><td>C<sub>11</sub></td><td>162</td><td>161</td><td>158 (169)</td><td>155</td><td>156</td><td>151</td><td>145 (172±6)</td><td>154</td><td>147</td><td>70 (98)</td><td></td><td></td></tr><tr><td>C<sub>12</sub></td><td>109</td><td>107</td><td>109 (127)</td><td>108</td><td>105</td><td>103</td><td>105 (108±1.5)</td><td>100</td><td>100</td><td>123 (83)</td><td></td><td></td></tr><tr><td>C<sub>44</sub></td><td>41</td><td>39</td><td>36 (45)</td><td>46</td><td>39</td><td>38</td><td>32 (28±1.5)</td><td>45</td><td>41</td><td>39 (38)</td><td></td><td></td></tr><tr><td>A<sub>z</sub></td><td>1.5</td><td>1.4</td><td>1.5 (2.1)</td><td>1.9</td><td>1.5</td><td>1.6</td><td>1.6 (0.87)</td><td>1.7</td><td>1.9</td><td>- (5)</td><td></td><td></td></tr></tbody></table>

in the supercell subjected to deformation. Furthermore, the prediction of ECs from DFT simulations is known to be limited in accuracy and is strongly influenced by the choice of XC function. That said, the results in Tab. 4 show a consistent and slight decrease of ECs ($C_{44}$ in particular) when measured for the LE set compared to the AVG and HE sets. These differences between the EC calculated for the LE and HE may be explained in two ways. i) The deeper energy wells corresponding to LE configurations may be circumvented by steeper energy variations, resulting in stiffer ECs. ii) HE configurations contain relatively more of unfavorable bonds, in particular bcc-hcp ones (cf. Figs. 6 and 7), possibly making the system less mechanically stable and lowering the ECs. Interestingly, the ECs ($C_{44}$) obtained for the HE configurations are in better agreement with experimental determinations for the quinary alloy [4], which were evaluated on the as-cast state that may well contain more high-energy atomic arrangements than in the final annealed state.

# 4 Conclusions

In this work, we have studied the structure, phase stability, and elastic behavior of $Ti_x(HfNbTaZr)_{(1-x)/4}$ alloys as a function of the Ti content $x$. The atomistic structure, corresponding to the disordered structure of the HEA solid solution, was approximated using a large number of different MCSQS realizations, while the atomic relaxation and
energy calculations were performed using DFT calculations for each MCSQS configuration.

The following results were obtained:

* We conducted a systematic sensitivity study on the number of NNs and on the order of the many-body terms employed to approximate the disordered atomic configuration in the MCSQS. In agreement with other studies, we found that the number of NN pairs has a stronger effect than higher bond orders to obtain SQS configurations corresponding to lower mixing energies.

* We show that the ideal solid solution is better approximated by considering several MCSQS realizations of the atomic disordered structure. The average mixing energy and global Short Range Order (SRO) parameter converge, however, faster than the SRO parameters for each bond type. The averages over multiple MCSQS configurations will converge faster toward the ideal solid solution representation for simpler alloy compositions and/or bigger supercells.

* The mixing energy of the alloys in the bcc structure is slightly positive; the value is, however, smaller than the absolute value of the bulk energy of individual elements. The bcc HEA phase observed experimentally for the equimolar quinary alloy could thus be metastable or stabilized thanks to SRO and/or vibrational entropy.

* We have proposed an original analysis of the short-range chemical ordering in the HfNbTaZr and Ti(HfNbTaZr) systems by making use of the large number of MCSQS realizations and DFT evaluations of the corresponding energies. Our SRO analysis strategy can be employed to quickly separate very favorable from unfavorable atomic pairs and is particularly efficient for simpler systems (i.e., quaternary alloys).
energy calculations were performed using DFT calculations for each MCSQS configuration.

The following results were obtained:

* We conducted a systematic sensitivity study on the number of NNs and on the order of the many-body terms employed to approximate the disordered atomic configuration in the MCSQS. In agreement with other studies, we found that the number of NN pairs has a stronger effect than higher bond orders to obtain SQS configurations corresponding to lower mixing energies.

* We show that the ideal solid solution is better approximated by considering several MCSQS realizations of the atomic disordered structure. The average mixing energy and global Short Range Order (SRO) parameter converge, however, faster than the SRO parameters for each bond type. The averages over multiple MCSQS configurations will converge faster toward the ideal solid solution representation for simpler alloy compositions and/or bigger supercells.

* The mixing energy of the alloys in the bcc structure is slightly positive; the value is, however, smaller than the absolute value of the bulk energy of individual elements. The bcc HEA phase observed experimentally for the equimolar quinary alloy could thus be metastable or stabilized thanks to SRO and/or vibrational entropy.

* We have proposed an original analysis of the short-range chemical ordering in the HfNbTaZr and Ti(HfNbTaZr) systems by making use of the large number of MCSQS realizations and DFT evaluations of the corresponding energies. Our SRO analysis strategy can be employed to quickly separate very favorable from unfavorable atomic pairs and is particularly efficient for simpler systems (i.e., quaternary alloys).
* We found that Ti-Zr, Ti-Hf, and Nb-Ta are particularly favored at short and longer distances (1 and 3 NN), probably leading to local fluctuations of the composition in the alloy. This could explain the waviness of dislocations and debris observed in experiments. Our SRO analysis is also mostly in agreement with the few existing studies relying on a more advanced MC method.

* To simplify the complex SRO analysis, elements with identical crystallographic structures are typically favorable, while elements with different crystallographic structures are unfavorable. The SRO in HfNbTaZr is mainly located on the 2nd NN shell, while in HfNbTaTiZr all the first three NN shells are involved. Few bond pairs are favorable in the quaternary system, and this could explain the weak global SRO in this system. Short-range ordering appears to be present in both systems, although it is stronger in the quinary system.

* The content of Ti has various effects on the structure or properties of the $Ti_x(HfNbTaZr)_{(1-x)/4}$ system. The average bcc lattice parameter and the mixing energies decrease monotonically with the Ti composition. Ti content can be tailored to adapt the elastic properties to specific applications ($C_{11}$ and $B$) by a few tens of GPa. The Ti content also controls the structure of the HEA phase, with a bcc phase stable, predicted as ductile, for $x < 0.5$. For Ti content larger than 0.5, the hcp phase becomes the most stable. Close to the quinary equimolar composition, a dual bcc-hcp phase or deformation induced phase transformation may be possible. A third possible phase in the shape of an orthorhombic structure is discarded when the mixing entropy term is considered. These results are consistent with the so-called Bo-Md diagram and mechanical stability analysis.
* We found that Ti-Zr, Ti-Hf, and Nb-Ta are particularly favored at short and longer distances (1 and 3 NN), probably leading to local fluctuations of the composition in the alloy. This could explain the waviness of dislocations and debris observed in experiments. Our SRO analysis is also mostly in agreement with the few existing studies relying on a more advanced MC method.

* To simplify the complex SRO analysis, elements with identical crystallographic structures are typically favorable, while elements with different crystallographic structures are unfavorable. The SRO in HfNbTaZr is mainly located on the 2nd NN shell, while in HfNbTaTiZr all the first three NN shells are involved. Few bond pairs are favorable in the quaternary system, and this could explain the weak global SRO in this system. Short-range ordering appears to be present in both systems, although it is stronger in the quinary system.

* The content of Ti has various effects on the structure or properties of the $Ti_x(HfNbTaZr)_{(1-x)/4}$ system. The average bcc lattice parameter and the mixing energies decrease monotonically with the Ti composition. Ti content can be tailored to adapt the elastic properties to specific applications ($C_{11}$ and $B$) by a few tens of GPa. The Ti content also controls the structure of the HEA phase, with a bcc phase stable, predicted as ductile, for $x < 0.5$. For Ti content larger than 0.5, the hcp phase becomes the most stable. Close to the quinary equimolar composition, a dual bcc-hcp phase or deformation induced phase transformation may be possible. A third possible phase in the shape of an orthorhombic structure is discarded when the mixing entropy term is considered. These results are consistent with the so-called Bo-Md diagram and mechanical stability analysis.
* We have assessed the impact of variability of local atomic arrangements on the structure, phase stability and elastic behaviour of these alloys. Local atomic configuration may affect the relative stability and explain some phase observed in experiments, while structural and elastic properties are mostly impacted by Ti content.

* The present work is qualitatively consistent with trends found in HEA systems and quantitatively consistent with the limited literature on this system. The average bcc lattice parameter and elastic constants agree with the experimental values obtained for the equimolar alloy.

# 5 Acknowledgements

We are grateful to the University Sorbonne, Paris Nord, Labex SEAM, ANR (Agence Nationale de la Recherche), and CGI (Commissariat à l'Investissement d'Avenir) for financing this work through Labex SEAM (Science and Engineering for Advanced Materials and devices), ANR-10-LABX-0096 and ANR-18-IDEX-0001. Further, we would like to thank Ing. Nicolas Greneche for HPC support, MAGI, Univ. Sorbonne Paris Nord, and HPC resources. We would like to thank Céline Varvenne, Guy Dirras and Pierre-Antoine Geslin for their helpful suggestions and discussions that improved the manuscript.
* We have assessed the impact of variability of local atomic arrangements on the structure, phase stability and elastic behaviour of these alloys. Local atomic configuration may affect the relative stability and explain some phase observed in experiments, while structural and elastic properties are mostly impacted by Ti content.

* The present work is qualitatively consistent with trends found in HEA systems and quantitatively consistent with the limited literature on this system. The average bcc lattice parameter and elastic constants agree with the experimental values obtained for the equimolar alloy.

# 5 Acknowledgements

We are grateful to the University Sorbonne, Paris Nord, Labex SEAM, ANR (Agence Nationale de la Recherche), and CGI (Commissariat à l'Investissement d'Avenir) for financing this work through Labex SEAM (Science and Engineering for Advanced Materials and devices), ANR-10-LABX-0096 and ANR-18-IDEX-0001. Further, we would like to thank Ing. Nicolas Greneche for HPC support, MAGI, Univ. Sorbonne Paris Nord, and HPC resources. We would like to thank Céline Varvenne, Guy Dirras and Pierre-Antoine Geslin for their helpful suggestions and discussions that improved the manuscript.
References

[1] M C Gao, C S Carney, N. Doğan, P D Jablonksi, J A Hawk, and D E Alman. Design of refractory high-entropy alloys. *JOM*, 67:2653–2669, 2015.

[2] D. B. Miracle and O. N. Senkov. A critical review of high entropy alloys and related concepts, 1 2017.

[3] Yuji Ikeda, Blazej Grabowski, and Fritz Körmann. Ab initio phase stabilities and mechanical properties of multicomponent alloys: A comprehensive review for high entropy alloys and compositionally complex alloys, 1 2019.

[4] G. Dirras, L. Lilensten, P. Djemia, M. Laurent-Brocq, D. Tingaud, J. P. Couzinié, L. Perrière, T. Chauveau, and I. Guillot. Elastic and plastic properties of as-cast equimolar tihfzrtanb high-entropy alloy. *Materials Science and Engineering A*, 654:30–38, 1 2016.

[5] O. N. Senkov, S. Rao, K. J. Chaput, and C. Woodward. Compositional effect on microstructure and properties of nbtizr-based complex concentrated alloys. *Acta Materialia*, 151:201–215, 6 2018.

[6] Michael C. Gao, Changning Niu, Chao Jiang, and Douglas L. Irving. Applications of special quasi-random structures to high-entropy alloys. *High-Entropy Alloys: Fundamentals and Applications*, pages 333–368, 1 2016.

[7] Akira Takeuchi, Kenji Amiya, Takeshi Wada, Kunio Yubuta, Wei Zhang, and Akihiro Makino. Entropies in alloy design for high-entropy and bulk glassy alloys. *Entropy*, 15:3810–3821, 2013.

[8] H. W. Yao, J. W. Qiao, J. A. Hawk, H. F. Zhou, M. W. Chen, and M. C. Gao. Mechanical properties of refractory high-entropy alloys: Experiments and modeling. *Journal of Alloys and Compounds*, 696:1139–1150, 2017.
References

[1] M C Gao, C S Carney, N. Doğan, P D Jablonksi, J A Hawk, and D E Alman. Design of refractory high-entropy alloys. *JOM*, 67:2653–2669, 2015.

[2] D. B. Miracle and O. N. Senkov. A critical review of high entropy alloys and related concepts, 1 2017.

[3] Yuji Ikeda, Blazej Grabowski, and Fritz Körmann. Ab initio phase stabilities and mechanical properties of multicomponent alloys: A comprehensive review for high entropy alloys and compositionally complex alloys, 1 2019.

[4] G. Dirras, L. Lilensten, P. Djemia, M. Laurent-Brocq, D. Tingaud, J. P. Couzinié, L. Perrière, T. Chauveau, and I. Guillot. Elastic and plastic properties of as-cast equimolar tihfzrtanb high-entropy alloy. *Materials Science and Engineering A*, 654:30–38, 1 2016.

[5] O. N. Senkov, S. Rao, K. J. Chaput, and C. Woodward. Compositional effect on microstructure and properties of nbtizr-based complex concentrated alloys. *Acta Materialia*, 151:201–215, 6 2018.

[6] Michael C. Gao, Changning Niu, Chao Jiang, and Douglas L. Irving. Applications of special quasi-random structures to high-entropy alloys. *High-Entropy Alloys: Fundamentals and Applications*, pages 333–368, 1 2016.

[7] Akira Takeuchi, Kenji Amiya, Takeshi Wada, Kunio Yubuta, Wei Zhang, and Akihiro Makino. Entropies in alloy design for high-entropy and bulk glassy alloys. *Entropy*, 15:3810–3821, 2013.

[8] H. W. Yao, J. W. Qiao, J. A. Hawk, H. F. Zhou, M. W. Chen, and M. C. Gao. Mechanical properties of refractory high-entropy alloys: Experiments and modeling. *Journal of Alloys and Compounds*, 696:1139–1150, 2017.
[9] Anna M. Manzoni and Uwe Glatzel. High-entropy alloys: Balancing strength and ductility at room temperature. *Reference Module in Materials Science and Materials Engineering*, 1 2020.

[10] Yuan Yuan, Yuan Wu, Zhi Yang, Xue Liang, Zhifeng Lei, Hailong Huang, Hui Wang, Xiongjun Liu, Ke An, Wei Wu, and Zhaoping Lu. Formation, structure and properties of biocompatible tizrhfnbta high-entropy alloys. *Materials Research Letters*, 7:225-231, 6 2019.

[11] S. Gurel, M. B. Yagci, B. Bal, and D. Canadinc. Corrosion behavior of novel titanium-based high entropy alloys designed for medical implants. *Materials Chemistry and Physics*, 254:123377, 11 2020.

[12] Mohamed Abdel-Hady, Keita Hinoshita, and Masahiko Morinaga. General approach to phase stability and elastic properties of β-type ti-alloys using electronic parameters. *Scripta Materialia*, 55:477-480, 9 2006.

[13] Lola Lilensten, Jean Philippe Couzinié, Julie Bourgon, Loïc Perrière, Guy Dirras, Frédéric Prima, and Ivan Guillot. Design and tensile properties of a bcc ti-rich high-entropy alloy with transformation-induced plasticity. *Materials Research Letters*, 5:110-116, 3 2017.

[14] Zhifeng Lei, Xiongjun Liu, Yuan Wu, Hui Wang, Suihe Jiang, Shudao Wang, Xidong Hui, Yidong Wu, Baptiste Gault, Paraskevas Kontis, Dierk Raabe, Lin Gu, Qinghua Zhang, Houwen Chen, Hongtao Wang, Jiabin Liu, Ke An, Qiaoshi Zeng, Tai Gang Nieh, and Zhaoping Lu. Enhanced strength and ductility in a high-entropy alloy via ordered oxygen complexes. *Nature*, 563:546-550, 11 2018.

[15] Yeqiang Bu, Yuan Wu, Zhifeng Lei, Xiaoyuan Yuan, Honghui Wu, Xiaobin Feng, Jiabin Liu, Jun Ding, Yang Lu, Hongtao Wang, Zhaoping Lu, and Wei Yang. Local
[9] Anna M. Manzoni and Uwe Glatzel. High-entropy alloys: Balancing strength and ductility at room temperature. *Reference Module in Materials Science and Materials Engineering*, 1 2020.

[10] Yuan Yuan, Yuan Wu, Zhi Yang, Xue Liang, Zhifeng Lei, Hailong Huang, Hui Wang, Xiongjun Liu, Ke An, Wei Wu, and Zhaoping Lu. Formation, structure and properties of biocompatible tizrhfnbta high-entropy alloys. *Materials Research Letters*, 7:225-231, 6 2019.

[11] S. Gurel, M. B. Yagci, B. Bal, and D. Canadinc. Corrosion behavior of novel titanium-based high entropy alloys designed for medical implants. *Materials Chemistry and Physics*, 254:123377, 11 2020.

[12] Mohamed Abdel-Hady, Keita Hinoshita, and Masahiko Morinaga. General approach to phase stability and elastic properties of β-type ti-alloys using electronic parameters. *Scripta Materialia*, 55:477-480, 9 2006.

[13] Lola Lilensten, Jean Philippe Couzinié, Julie Bourgon, Loïc Perrière, Guy Dirras, Frédéric Prima, and Ivan Guillot. Design and tensile properties of a bcc ti-rich high-entropy alloy with transformation-induced plasticity. *Materials Research Letters*, 5:110-116, 3 2017.

[14] Zhifeng Lei, Xiongjun Liu, Yuan Wu, Hui Wang, Suihe Jiang, Shudao Wang, Xidong Hui, Yidong Wu, Baptiste Gault, Paraskevas Kontis, Dierk Raabe, Lin Gu, Qinghua Zhang, Houwen Chen, Hongtao Wang, Jiabin Liu, Ke An, Qiaoshi Zeng, Tai Gang Nieh, and Zhaoping Lu. Enhanced strength and ductility in a high-entropy alloy via ordered oxygen complexes. *Nature*, 563:546-550, 11 2018.

[15] Yeqiang Bu, Yuan Wu, Zhifeng Lei, Xiaoyuan Yuan, Honghui Wu, Xiaobin Feng, Jiabin Liu, Jun Ding, Yang Lu, Hongtao Wang, Zhaoping Lu, and Wei Yang. Local
chemical fluctuation mediated ductility in body-centered-cubic high-entropy alloys.
*Materials Today*, 46:28–34, 6 2021.

[16] Ruopeng Zhang, Shiteng Zhao, Jun Ding, Yan Chong, Tao Jia, Colin Ophus, Mark Asta, Robert O. Ritchie, and Andrew M. Minor. Short-range order and its impact on the crconi medium-entropy alloy. *Nature*, 581:283–287, 5 2020.

[17] Sheng Yin, Jun Ding, Mark Asta, and Robert O Ritchie. Ab initio modeling of the energy landscape for screw dislocations in body-centered cubic high-entropy alloys. *npj Computational Materials*, 6, 2020.

[18] Kaihui Xun, Bozhao Zhang, Qi Wang, Zhen Zhang, Jun Ding, and En Ma. Local chemical inhomogeneities in tizrnb-based refractory high-entropy alloys. *Journal of Materials Science & Technology*, 135:221–230, 2 2023.

[19] Shuai Chen, Zachary H. Aitken, Subrahmanyam Pattamatta, Zhaoxuan Wu, Zhi Gen Yu, David J. Srolovitz, Peter K. Liaw, and Yong Wei Zhang. Simultaneously enhancing the ultimate strength and ductility of high-entropy alloys via short-range ordering. *Nature Communications* 2021 12:1, 12:1–11, 8 2021.

[20] Bozhao Zhang, Jun Ding, and En Ma. Chemical short-range order in body-centered-cubic tizrhfnb high-entropy alloys. *Applied Physics Letters*, 119, 11 2021.

[21] S. D. Wang, X. J. Liu, Z. F. Lei, D. Y. Lin, F. G. Bian, C. M. Yang, M. Y. Jiao, Q. Du, H. Wang, Y. Wu, S. H. Jiang, and Z. P. Lu. Chemical short-range ordering and its strengthening effect in refractory high-entropy alloys. *Physical Review B*, 103:104107, 3 2021.

[22] Binglun ; Yin, Francesco ; Maresca, William Curtin, and Acta Materialia. Vanadium is an optimal element for strengthening in both fcc and bcc high-entropy alloys. 2020.

[23] Binglun Yin and William A. Curtin. First-principles-based prediction of yield strength in the rhirpdptnicu high-entropy alloy. *npj Computational Materials*, 5:14, 12 2019.
chemical fluctuation mediated ductility in body-centered-cubic high-entropy alloys.
*Materials Today*, 46:28–34, 6 2021.

[16] Ruopeng Zhang, Shiteng Zhao, Jun Ding, Yan Chong, Tao Jia, Colin Ophus, Mark Asta, Robert O. Ritchie, and Andrew M. Minor. Short-range order and its impact on the crconi medium-entropy alloy. *Nature*, 581:283–287, 5 2020.

[17] Sheng Yin, Jun Ding, Mark Asta, and Robert O Ritchie. Ab initio modeling of the energy landscape for screw dislocations in body-centered cubic high-entropy alloys. *npj Computational Materials*, 6, 2020.

[18] Kaihui Xun, Bozhao Zhang, Qi Wang, Zhen Zhang, Jun Ding, and En Ma. Local chemical inhomogeneities in tizrnb-based refractory high-entropy alloys. *Journal of Materials Science & Technology*, 135:221–230, 2 2023.

[19] Shuai Chen, Zachary H. Aitken, Subrahmanyam Pattamatta, Zhaoxuan Wu, Zhi Gen Yu, David J. Srolovitz, Peter K. Liaw, and Yong Wei Zhang. Simultaneously enhancing the ultimate strength and ductility of high-entropy alloys via short-range ordering. *Nature Communications* 2021 12:1, 12:1–11, 8 2021.

[20] Bozhao Zhang, Jun Ding, and En Ma. Chemical short-range order in body-centered-cubic tizrhfnb high-entropy alloys. *Applied Physics Letters*, 119, 11 2021.

[21] S. D. Wang, X. J. Liu, Z. F. Lei, D. Y. Lin, F. G. Bian, C. M. Yang, M. Y. Jiao, Q. Du, H. Wang, Y. Wu, S. H. Jiang, and Z. P. Lu. Chemical short-range ordering and its strengthening effect in refractory high-entropy alloys. *Physical Review B*, 103:104107, 3 2021.

[22] Binglun ; Yin, Francesco ; Maresca, William Curtin, and Acta Materialia. Vanadium is an optimal element for strengthening in both fcc and bcc high-entropy alloys. 2020.

[23] Binglun Yin and William A. Curtin. First-principles-based prediction of yield strength in the rhirpdptnicu high-entropy alloy. *npj Computational Materials*, 5:14, 12 2019.
[24] Céline Varvenne, Aitor Luque, and William A. Curtin. Theory of strengthening in fcc high entropy alloys. *Acta Materialia*, 118:164-176, 10 2016.

[25] L. Lilensten, J. P. Couzinié, L. Perrière, A. Hocini, C. Keller, G. Dirras, and I. Guillot. Study of a bcc multi-principal element alloy: Tensile and simple shear properties and underlying deformation mechanisms. *Acta Materialia*, 142:131-141, 1 2018.

[26] J. Ph Couzinié, L. Lilensten, Y. Champion, G. Dirras, L. Perrière, and I. Guillot. On the room temperature deformation mechanisms of a tizrhfnbta refractory high-entropy alloy. *Materials Science and Engineering: A*, 645:255-263, 10 2015.

[27] Lola Lilensten, Karine Provost, Loïc Perrière, Emiliano Fonda, Jean Philippe Couzinié, Fabienne Amann, Martin Radtke, Guy Dirras, and Ivan Guillot. Experimental investigation of the local environment and lattice distortion in refractory medium entropy alloys. *Scripta Materialia*, 211:114532, 4 2022.

[28] Yuan Yuan Tan, Tong Li, Yu Chen, Zhong Jun Chen, Ming Yao Su, Jing Zhang, Yu Gong, Tao Wu, Hai Ying Wang, and Lan Hong Dai. Uncovering heterogeneity of local lattice distortion in tizrhfnbta refractory high entropy alloy by sr-xrd and exafs. *Scripta Materialia*, 223:115079, 1 2023.

[29] S I Rao, B Akdim, E Antillon, C Woodward, T A Parthasarathy, and O N Senkov. Modeling solution hardening in bcc refractory complex concentrated alloys: Nbtizr, nb1.5tizr0.5 and nb0.5tizr1.5. *Acta Materialia*, 168:222-236, 2019.

[30] Xiusong Huang, Lehua Liu, Xianbao Duan, Weibing Liao, Jianjun Huang, Huibin Sun, and Chunyan Yu. Atomistic simulation of chemical short-range order in hfnbtazr high entropy alloy based on a newly-developed interatomic potential. *Materials & Design*, 202:109560, 4 2021.

[31] G. Kresse and J. Furthmüller. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. *Physical Review B*, 54:11169-11186, 10 1996.
[24] Céline Varvenne, Aitor Luque, and William A. Curtin. Theory of strengthening in fcc high entropy alloys. *Acta Materialia*, 118:164-176, 10 2016.

[25] L. Lilensten, J. P. Couzinié, L. Perrière, A. Hocini, C. Keller, G. Dirras, and I. Guillot. Study of a bcc multi-principal element alloy: Tensile and simple shear properties and underlying deformation mechanisms. *Acta Materialia*, 142:131-141, 1 2018.

[26] J. Ph Couzinié, L. Lilensten, Y. Champion, G. Dirras, L. Perrière, and I. Guillot. On the room temperature deformation mechanisms of a tizrhfnbta refractory high-entropy alloy. *Materials Science and Engineering: A*, 645:255-263, 10 2015.

[27] Lola Lilensten, Karine Provost, Loïc Perrière, Emiliano Fonda, Jean Philippe Couzinié, Fabienne Amann, Martin Radtke, Guy Dirras, and Ivan Guillot. Experimental investigation of the local environment and lattice distortion in refractory medium entropy alloys. *Scripta Materialia*, 211:114532, 4 2022.

[28] Yuan Yuan Tan, Tong Li, Yu Chen, Zhong Jun Chen, Ming Yao Su, Jing Zhang, Yu Gong, Tao Wu, Hai Ying Wang, and Lan Hong Dai. Uncovering heterogeneity of local lattice distortion in tizrhfnbta refractory high entropy alloy by sr-xrd and exafs. *Scripta Materialia*, 223:115079, 1 2023.

[29] S I Rao, B Akdim, E Antillon, C Woodward, T A Parthasarathy, and O N Senkov. Modeling solution hardening in bcc refractory complex concentrated alloys: Nbtizr, nb1.5tizr0.5 and nb0.5tizr1.5. *Acta Materialia*, 168:222-236, 2019.

[30] Xiusong Huang, Lehua Liu, Xianbao Duan, Weibing Liao, Jianjun Huang, Huibin Sun, and Chunyan Yu. Atomistic simulation of chemical short-range order in hfnbtazr high entropy alloy based on a newly-developed interatomic potential. *Materials & Design*, 202:109560, 4 2021.

[31] G. Kresse and J. Furthmüller. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. *Physical Review B*, 54:11169-11186, 10 1996.
[32] G. Kresse and D. Joubert. From ultrasoft pseudopotentials to the projector augmented-wave method. *Physical Review B*, 59:1758-1775, 1999.

[33] John P. Perdew, Kieron Burke, and Matthias Ernzerhof. Generalized gradient approximation made simple. *Physical Review Letters*, 77:3865-3868, 1996.

[34] Colin W. Glass, Artem R. Oganov, and Nikolaus Hansen. Uspex-evolutionary crystal structure prediction. *Computer Physics Communications*, 2006.

[35] Paul Soven. Coherent-potential model of substitutional disordered alloys. *Physical Review*, 156:809-813, 4 1967.

[36] Michael C. Gao, Pan Gao, Jeffrey A. Hawk, Lizhi Ouyang, David E. Alman, and Mike Widom. Computational modeling of high-entropy alloys: Structures, thermodynamics and elasticity, 10 2017.

[37] Hongquan Song, Fuyang Tian, Qing Miao Hu, Levente Vitos, Yandong Wang, Jiang Shen, and Nanxian Chen. Local lattice distortion in high-entropy alloys. *Physical Review Materials*, 1:23404, 2017.

[38] A. E. Carlsson and P. A. Fedders. Maximum-entropy method for electronic properties of alloys. *Physical Review B*, 34:3567-3571, 9 1986.

[39] A. van de Walle, P Tiwary, M. de Jong, D.L. Olmsted, M Asta, A Dick, D Shin, Y Wang, L.-Q. Chen, and Z.-K. Liu. Efficient stochastic generation of special quasirandom structures. *Calphad*, 42:13-18, 9 2013.

[40] Axel van de Walle. Multicomponent multisublattice alloys, nonconfigurational entropy and other additions to the alloy theoretic automated toolkit. *Calphad*, 33:266-278, 6 2009.

[41] M. Methfessel and A. T. Paxton. High-precision sampling for brillouin-zone integration in metals. *Physical Review B*, 40:3616-3621, 8 1989.
[32] G. Kresse and D. Joubert. From ultrasoft pseudopotentials to the projector augmented-wave method. *Physical Review B*, 59:1758-1775, 1999.

[33] John P. Perdew, Kieron Burke, and Matthias Ernzerhof. Generalized gradient approximation made simple. *Physical Review Letters*, 77:3865-3868, 1996.

[34] Colin W. Glass, Artem R. Oganov, and Nikolaus Hansen. Uspex-evolutionary crystal structure prediction. *Computer Physics Communications*, 2006.

[35] Paul Soven. Coherent-potential model of substitutional disordered alloys. *Physical Review*, 156:809-813, 4 1967.

[36] Michael C. Gao, Pan Gao, Jeffrey A. Hawk, Lizhi Ouyang, David E. Alman, and Mike Widom. Computational modeling of high-entropy alloys: Structures, thermodynamics and elasticity, 10 2017.

[37] Hongquan Song, Fuyang Tian, Qing Miao Hu, Levente Vitos, Yandong Wang, Jiang Shen, and Nanxian Chen. Local lattice distortion in high-entropy alloys. *Physical Review Materials*, 1:23404, 2017.

[38] A. E. Carlsson and P. A. Fedders. Maximum-entropy method for electronic properties of alloys. *Physical Review B*, 34:3567-3571, 9 1986.

[39] A. van de Walle, P Tiwary, M. de Jong, D.L. Olmsted, M Asta, A Dick, D Shin, Y Wang, L.-Q. Chen, and Z.-K. Liu. Efficient stochastic generation of special quasirandom structures. *Calphad*, 42:13-18, 9 2013.

[40] Axel van de Walle. Multicomponent multisublattice alloys, nonconfigurational entropy and other additions to the alloy theoretic automated toolkit. *Calphad*, 33:266-278, 6 2009.

[41] M. Methfessel and A. T. Paxton. High-precision sampling for brillouin-zone integration in metals. *Physical Review B*, 40:3616-3621, 8 1989.
[42] Dario Alfè. Phon: A program to calculate phonons using the small displacement method. *Computer Physics Communications*, 180:2622-2633, 12 2009.

[43] Atsushi Togo and Isao Tanaka. First principles phonon calculations in materials science. *Scripta Materialia*, 108:1-5, 11 2015.

[44] Zhenyu Du, Jie Zuo, Nanyun Bao, Mingli Yang, Gang Jiang, and Li Zhang. Effect of ta addition on the structural, thermodynamic and mechanical properties of cocrfeni high entropy alloys. *RSC Advances*, 9:16447-16454, 5 2019.

[45] Damian Sobieraj, Jan S. Wróbel, Tomasz Rygier, Krzysztof J. Kurzydłowski, Osman El Atwani, Arun Devaraj, Enrique Martinez Saez, and Duc Nguyen-Manh. Chemical short-range order in derivative cr-ta-ti-v-w high entropy alloys from the first-principles thermodynamic study. *Physical Chemistry Chemical Physics*, 22:23929-23951, 11 2020.

[46] Fuyang Tian, De Ye Lin, Xingyu Gao, Ya Fan Zhao, and Hai Feng Song. A structural modeling approach to solid solutions based on the similar atomic environment. *Journal of Chemical Physics*, 153, 2020.

[47] D. de Fontaine. The number of independent pair-correlation functions in multicomponent systems. *Journal of Applied Crystallography*, 4:15-19, 1971.

[48] Sarath Menon, Grisell Díaz Leines, and Jutta Rogal. pyscal: A python module for structural analysis of atomic environments. *Journal of Open Source Software*, 4:1824, 11 2019.

[49] Sheng Yin, Yunxing Zuo, Anas Abu-Odeh, Hui Zheng, Xiang Guo Li, Jun Ding, Shyue Ping Ong, Mark Asta, and Robert O. Ritchie. Atomistic simulations of dislocation mobility in refractory high-entropy alloys and the effect of chemical short-range order. *Nature Communications*, 12, 12 2021.
[42] Dario Alfè. Phon: A program to calculate phonons using the small displacement method. *Computer Physics Communications*, 180:2622-2633, 12 2009.

[43] Atsushi Togo and Isao Tanaka. First principles phonon calculations in materials science. *Scripta Materialia*, 108:1-5, 11 2015.

[44] Zhenyu Du, Jie Zuo, Nanyun Bao, Mingli Yang, Gang Jiang, and Li Zhang. Effect of ta addition on the structural, thermodynamic and mechanical properties of cocrfeni high entropy alloys. *RSC Advances*, 9:16447-16454, 5 2019.

[45] Damian Sobieraj, Jan S. Wróbel, Tomasz Rygier, Krzysztof J. Kurzydłowski, Osman El Atwani, Arun Devaraj, Enrique Martinez Saez, and Duc Nguyen-Manh. Chemical short-range order in derivative cr-ta-ti-v-w high entropy alloys from the first-principles thermodynamic study. *Physical Chemistry Chemical Physics*, 22:23929-23951, 11 2020.

[46] Fuyang Tian, De Ye Lin, Xingyu Gao, Ya Fan Zhao, and Hai Feng Song. A structural modeling approach to solid solutions based on the similar atomic environment. *Journal of Chemical Physics*, 153, 2020.

[47] D. de Fontaine. The number of independent pair-correlation functions in multicomponent systems. *Journal of Applied Crystallography*, 4:15-19, 1971.

[48] Sarath Menon, Grisell Díaz Leines, and Jutta Rogal. pyscal: A python module for structural analysis of atomic environments. *Journal of Open Source Software*, 4:1824, 11 2019.

[49] Sheng Yin, Yunxing Zuo, Anas Abu-Odeh, Hui Zheng, Xiang Guo Li, Jun Ding, Shyue Ping Ong, Mark Asta, and Robert O. Ritchie. Atomistic simulations of dislocation mobility in refractory high-entropy alloys and the effect of chemical short-range order. *Nature Communications*, 12, 12 2021.
[50] M. Morinaga, N. Yukawa, T. Maya, K. Sone, and H. Adachi. Theoretical design of titanium alloys, 1988.

[51] O. N. Senkov, J. M. Scott, S. V. Senkova, D. B. Miracle, and C. F. Woodward. Microstructure and room temperature properties of a high-entropy tanbhfzrti alloy. *Journal of Alloys and Compounds*, 509:6043-6048, 5 2011.

[52] Michael C. Gao, Peter K. Liaw, Jien Wei Yeh, and Yong Zhang. *High-entropy alloys: Fundamentals and applications*. Springer International Publishing, 1 2016.

[53] L. Lilensten, J.P. Couzinié, L. Perrière, J. Bourgon, N. Emery, and I. Guillot. New structure in refractory high-entropy alloys. *Materials Letters*, 132:123-125, October 2014.

[54] F. Otto, Y. Yang, H. Bei, and E. P. George. Relative effects of enthalpy and entropy on the phase stability of equiatomic high-entropy alloys. *Acta Materialia*, 61:2628-2638, 4 2013.

[55] Mathilde Laurent-Brocq, Loïc Perrière, Rémy Pirès, and Yannick Champion. From high entropy alloys to diluted multi-component alloys: Range of existence of a solid-solution. *Materials and Design*, 103:84-89, 8 2016.

[56] I. A. Tomilin and S. D. Kaloshkin. 'high entropy alloys'-'semi-impossible' regular solid solutions? *Materials Science and Technology (United Kingdom)*, 31:1231-1234, 7 2015.

[57] J.D. Cox. *CODATA Key Values for Thermodynamics*. Hemisphere Publ. Co., NY, 1989.

[58] P. T. Jochym and K. Parlinski. Ab initio lattice dynamics and elastic constants of zrc. *European Physical Journal B*, 15:265-268, 5 2000.
[50] M. Morinaga, N. Yukawa, T. Maya, K. Sone, and H. Adachi. Theoretical design of titanium alloys, 1988.

[51] O. N. Senkov, J. M. Scott, S. V. Senkova, D. B. Miracle, and C. F. Woodward. Microstructure and room temperature properties of a high-entropy tanbhfzrti alloy. *Journal of Alloys and Compounds*, 509:6043-6048, 5 2011.

[52] Michael C. Gao, Peter K. Liaw, Jien Wei Yeh, and Yong Zhang. *High-entropy alloys: Fundamentals and applications*. Springer International Publishing, 1 2016.

[53] L. Lilensten, J.P. Couzinié, L. Perrière, J. Bourgon, N. Emery, and I. Guillot. New structure in refractory high-entropy alloys. *Materials Letters*, 132:123-125, October 2014.

[54] F. Otto, Y. Yang, H. Bei, and E. P. George. Relative effects of enthalpy and entropy on the phase stability of equiatomic high-entropy alloys. *Acta Materialia*, 61:2628-2638, 4 2013.

[55] Mathilde Laurent-Brocq, Loïc Perrière, Rémy Pirès, and Yannick Champion. From high entropy alloys to diluted multi-component alloys: Range of existence of a solid-solution. *Materials and Design*, 103:84-89, 8 2016.

[56] I. A. Tomilin and S. D. Kaloshkin. 'high entropy alloys'-'semi-impossible' regular solid solutions? *Materials Science and Technology (United Kingdom)*, 31:1231-1234, 7 2015.

[57] J.D. Cox. *CODATA Key Values for Thermodynamics*. Hemisphere Publ. Co., NY, 1989.

[58] P. T. Jochym and K. Parlinski. Ab initio lattice dynamics and elastic constants of zrc. *European Physical Journal B*, 15:265-268, 5 2000.
[59] Massel Ledbetter, Hirotsugu Ogi, Satoshi Kai, Sudook Kim, and Masahiko Hirao. Elastic constants of body-centered-cubic titanium monocrystals. *Journal of Applied Physics*, 95:4642-4644, 5 2004.

[60] E. Fazakas, V. Zadorozhnyy, L. K. Varga, A. Inoue, D. V. Louzguine-Luzgin, Fuyang Tian, and L. Vitos. Experimental and theoretical study of ti20zr20hf20nb20x20 (x = v or cr) refractory high-entropy alloys. *International Journal of Refractory Metals and Hard Materials*, 47:131-138, 11 2014.

[61] Wenqiang Feng, Yang Qi, and Shaoqing Wang. Effects of short-range order on the magnetic and mechanical properties of feconi(alsi)x high entropy alloys. *Metals*, 7:482, 11 2017.

[62] R. Hill. The elastic behaviour of a crystalline aggregate. *Proceedings of the Physical Society. Section A*, 65:349-354, 5 1952.

[63] S.F. Pugh. Xcii. relations between the elastic moduli and the plastic properties of polycrystalline pure metals. *The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science*, 45:823-843, 8 1954.
[59] Massel Ledbetter, Hirotsugu Ogi, Satoshi Kai, Sudook Kim, and Masahiko Hirao. Elastic constants of body-centered-cubic titanium monocrystals. *Journal of Applied Physics*, 95:4642-4644, 5 2004.

[60] E. Fazakas, V. Zadorozhnyy, L. K. Varga, A. Inoue, D. V. Louzguine-Luzgin, Fuyang Tian, and L. Vitos. Experimental and theoretical study of ti20zr20hf20nb20x20 (x = v or cr) refractory high-entropy alloys. *International Journal of Refractory Metals and Hard Materials*, 47:131-138, 11 2014.

[61] Wenqiang Feng, Yang Qi, and Shaoqing Wang. Effects of short-range order on the magnetic and mechanical properties of feconi(alsi)x high entropy alloys. *Metals*, 7:482, 11 2017.

[62] R. Hill. The elastic behaviour of a crystalline aggregate. *Proceedings of the Physical Society. Section A*, 65:349-354, 5 1952.

[63] S.F. Pugh. Xcii. relations between the elastic moduli and the plastic properties of polycrystalline pure metals. *The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science*, 45:823-843, 8 1954.
Supplementary Information

Tab. S1 shows the structural parameters for the configurations corresponding to the lowest mixing energy structures. The relaxed lattice parameters (a, b, c) deviate slightly from the perfect bcc lattice. They are slightly anisotropic (with an average difference of about 3%), and the angles (α, β, γ) between the basis vectors of the supercell slightly differ from the expected 90° (by typically less than 1°). These small distortions are unsurprising for HEA systems, as differences in the local chemical order are expected along the supercell axis, especially for the small dimensions dealt with in DFT.

The known stable experimental lattice parameters at RT for stable elements are $a_{Ti}$ = 2.95 Å, $c/a_{Ti}$ = 1.58 Å, $a_{Zr}$ = 3.23 Å, $c/a_{Zr}$ = 1.59 Å, $a_{Hf}$ = 3.196 Å, $c/a_{Hf}$ = 1.58 Å, $a_{Nb}$ = 3.30 Å, and $a_{Ta}$ = 3.31 Å. The lattice constants for a pure Ti (hcp) unit cell are a = b = 2.950 Å and c = 4.687 Å. For large Ti content, typically $x \ge 0.5$, the Bo-Md diagram suggests that hcp and/or orthorhombic (ortho) structures may be stable. We thus performed another set of MCSQS+DFT calculations for these crystallographic structures. We adopted the same methodology as discussed in Sec. 2. To assess the prediction of the Bo-Md diagram, we generated a set of MCSQS structures, assuming the hcp phase for x ranges from 0.5 to 0.81.

The bottom of Tab. S1 (hcp and ortho) shows the lattice geometry obtained for these simulations. The simulated lattices are again slightly distorted when compared to a perfect hcp lattice, certainly because the different element arrangements are probed along the supercell axis. Surprisingly, no clear trend was found regarding the evolution of the c/a ratio as the Ti content increased.
Supplementary Information

Tab. S1 shows the structural parameters for the configurations corresponding to the lowest mixing energy structures. The relaxed lattice parameters (a, b, c) deviate slightly from the perfect bcc lattice. They are slightly anisotropic (with an average difference of about 3%), and the angles (α, β, γ) between the basis vectors of the supercell slightly differ from the expected 90° (by typically less than 1°). These small distortions are unsurprising for HEA systems, as differences in the local chemical order are expected along the supercell axis, especially for the small dimensions dealt with in DFT.

The known stable experimental lattice parameters at RT for stable elements are $a_{Ti}$ = 2.95 Å, $c/a_{Ti}$ = 1.58 Å, $a_{Zr}$ = 3.23 Å, $c/a_{Zr}$ = 1.59 Å, $a_{Hf}$ = 3.196 Å, $c/a_{Hf}$ = 1.58 Å, $a_{Nb}$ = 3.30 Å, and $a_{Ta}$ = 3.31 Å. The lattice constants for a pure Ti (hcp) unit cell are $a = b = 2.950$ Å and $c = 4.687$ Å. For large Ti content, typically $x \ge 0.5$, the Bo-Md diagram suggests that hcp and/or orthorhombic (ortho) structures may be stable. We thus performed another set of MCSQS+DFT calculations for these crystallographic structures. We adopted the same methodology as discussed in Sec. 2. To assess the prediction of the Bo-Md diagram, we generated a set of MCSQS structures, assuming the hcp phase for $x$ ranges from 0.5 to 0.81.

The bottom of Tab. S1 (hcp and ortho) shows the lattice geometry obtained for these simulations. The simulated lattices are again slightly distorted when compared to a perfect hcp lattice, certainly because the different element arrangements are probed along the supercell axis. Surprisingly, no clear trend was found regarding the evolution of the $c/a$ ratio as the Ti content increased.
<table><caption>Table S1: Detailed analysis of the structure of the atomic configurations after relaxation for Ti<sub>x</sub>(HfNbTaZr)<sub>(1-x)</sub>/4 alloys as a function of the Ti content x and the choice of crystallographic structure. The average lattice parameters a' are reported for the conventional unit cell. The number of atoms is given for the SQS supercell generated from MCSQS. Known experimental values are reported from [30, 51] and are given in parentheses.</caption><thead><tr><th rowspan="2">Phases</th><th rowspan="2">x</th><th rowspan="2">#atoms</th><th colspan="3">Å</th><th rowspan="2">α°</th><th rowspan="2">β°</th><th rowspan="2">γ°</th><th colspan="2">Å</th></tr><tr><th>a</th><th>b</th><th>c</th><th>a'</th><th>a'</th></tr></thead><tbody><tr><td rowspan="10">bcc</td><td>0.00</td><td>128</td><td>3.454</td><td>3.417</td><td>3.438</td><td>89.6</td><td>90.0</td><td>90.0</td><td>3.436(3.446)</td><td></td></tr><tr><td>0.11</td><td>54</td><td>3.429</td><td>3.408</td><td>3.425</td><td>89.5</td><td>89.8</td><td>89.9</td><td>3.420</td><td></td></tr><tr><td>0.20</td><td>125</td><td>3.398</td><td>3.417</td><td>3.405</td><td>89.7</td><td>89.6</td><td>89.8</td><td>3.407(3.404)</td><td></td></tr><tr><td>0.33</td><td>54</td><td>3.392</td><td>3.372</td><td>3.381</td><td>89.7</td><td>90.1</td><td>89.9</td><td>3.381</td><td></td></tr><tr><td>0.40</td><td>54</td><td>3.371</td><td>3.336</td><td>3.403</td><td>89.5</td><td>89.9</td><td>89.6</td><td>3.370</td><td></td></tr><tr><td>0.50</td><td>128</td><td>3.379</td><td>3.321</td><td>3.361</td><td>89.9</td><td>90.8</td><td>90.0</td><td>3.353</td><td></td></tr><tr><td>0.63</td><td>54</td><td>3.313</td><td>3.345</td><td>3.339</td><td>90.6</td><td>88.6</td><td>90.7</td><td>3.332</td><td></td></tr><tr><td>0.70</td><td>54</td><td>3.300</td><td>3.305</td><td>3.350</td><td>89.2</td><td>90.4</td><td>90.5</td><td>3.318</td><td></td></tr><tr><td>0.81</td><td>128</td><td>3.073</td><td>3.421</td><td>3.418</td><td>85.8</td><td>89.9</td><td>90.1</td><td>3.306</td><td></td></tr><tr><td>1.00</td><td>2</td><td>3.252</td><td>3.252</td><td>3.252</td><td>90.0</td><td>90.0</td><td>90.0</td><td>3.252</td><td></td></tr><tr><td rowspan="5">hcp</td><td>0.50</td><td>32</td><td>3.430</td><td>4.684</td><td>4.707</td><td>86.9</td><td>88.2</td><td>88.0</td><td></td><td></td></tr><tr><td>0.50</td><td>128</td><td>2.933</td><td>2.942</td><td>4.680</td><td>90.1</td><td>89.7</td><td>110.9</td><td>c/a' = 1.593</td><td></td></tr><tr><td>0.63</td><td>54</td><td>3.318</td><td>2.961</td><td>4.543</td><td>89.9</td><td>90.0</td><td>123.8</td><td>c/a' = 1.447</td><td></td></tr><tr><td>0.70</td><td>54</td><td>2.953</td><td>3.285</td><td>4.533</td><td>89.7</td><td>89.6</td><td>123.7</td><td>c/a' = 1.453</td><td></td></tr><tr><td>0.81</td><td>128</td><td>2.934</td><td>2.997</td><td>4.687</td><td>90.0</td><td>89.9</td><td>119.3</td><td>c/a' = 1.580</td><td></td></tr></tbody></table>

<table><caption>Table S2: Elastic constants $C_{ij}$ (GPa) and Zener anisotropy coefficient $A_z$ calculated for the bcc SQS supercell. The elastic constants of bcc Ti at high temperature are also given for comparison. The experimental EC with isotropic Zener Coefficient is given in parentheses for the known composition.</caption><thead><tr><th>x</th><th>0.0 ([61])</th><th>0.11</th><th>0.20 ([4])</th><th>0.33</th><th>0.40</th><th>1.0 ([59])</th></tr></thead><tbody><tr><td>C<sub>11</sub></td><td>162 (169)</td><td>155</td><td>156 (172±6)</td><td>154</td><td>147</td><td>70(98)</td></tr><tr><td>C<sub>12</sub></td><td>109 (127)</td><td>108</td><td>105 (108±1.5)</td><td>100</td><td>100</td><td>123(83)</td></tr><tr><td>C<sub>44</sub></td><td>42 (45)</td><td>46</td><td>39 (28±1.5)</td><td>45</td><td>41</td><td>39(38)</td></tr><tr><td>A<sub>z</sub></td><td>1.5 (2.1)</td><td>1.9</td><td>1.5 (0.87)</td><td>1.7</td><td>1.9</td><td>- (5.1)</td></tr></tbody></table>
<table><caption>Table S1: Detailed analysis of the structure of the atomic configurations after relaxation for Ti<sub>x</sub>(HfNbTaZr)<sub>(1-x)</sub>/4 alloys as a function of the Ti content x and the choice of crystallographic structure. The average lattice parameters a' are reported for the conventional unit cell. The number of atoms is given for the SQS supercell generated from MCSQS. Known experimental values are reported from [30, 51] and are given in parentheses.</caption><thead><tr><th rowspan="2">Phases</th><th rowspan="2">x</th><th rowspan="2">#atoms</th><th colspan="3">Å</th><th rowspan="2">α°</th><th rowspan="2">β°</th><th rowspan="2">γ°</th><th colspan="2">Å</th></tr><tr><th>a</th><th>b</th><th>c</th><th>a'</th><th>a'</th></tr></thead><tbody><tr><td rowspan="10">bcc</td><td>0.00</td><td>128</td><td>3.454</td><td>3.417</td><td>3.438</td><td>89.6</td><td>90.0</td><td>90.0</td><td>3.436(3.446)</td><td></td></tr><tr><td>0.11</td><td>54</td><td>3.429</td><td>3.408</td><td>3.425</td><td>89.5</td><td>89.8</td><td>89.9</td><td>3.420</td><td></td></tr><tr><td>0.20</td><td>125</td><td>3.398</td><td>3.417</td><td>3.405</td><td>89.7</td><td>89.6</td><td>89.8</td><td>3.407(3.404)</td><td></td></tr><tr><td>0.33</td><td>54</td><td>3.392</td><td>3.372</td><td>3.381</td><td>89.7</td><td>90.1</td><td>89.9</td><td>3.381</td><td></td></tr><tr><td>0.40</td><td>54</td><td>3.371</td><td>3.336</td><td>3.403</td><td>89.5</td><td>89.9</td><td>89.6</td><td>3.370</td><td></td></tr><tr><td>0.50</td><td>128</td><td>3.379</td><td>3.321</td><td>3.361</td><td>89.9</td><td>90.8</td><td>90.0</td><td>3.353</td><td></td></tr><tr><td>0.63</td><td>54</td><td>3.313</td><td>3.345</td><td>3.339</td><td>90.6</td><td>88.6</td><td>90.7</td><td>3.332</td><td></td></tr><tr><td>0.70</td><td>54</td><td>3.300</td><td>3.305</td><td>3.350</td><td>89.2</td><td>90.4</td><td>90.5</td><td>3.318</td><td></td></tr><tr><td>0.81</td><td>128</td><td>3.073</td><td>3.421</td><td>3.418</td><td>85.8</td><td>89.9</td><td>90.1</td><td>3.306</td><td></td></tr><tr><td>1.00</td><td>2</td><td>3.252</td><td>3.252</td><td>3.252</td><td>90.0</td><td>90.0</td><td>90.0</td><td>3.252</td><td></td></tr><tr><td rowspan="5">hcp</td><td>0.50</td><td>32</td><td>3.430</td><td>4.684</td><td>4.707</td><td>86.9</td><td>88.2</td><td>88.0</td><td></td><td></td></tr><tr><td>0.50</td><td>128</td><td>2.933</td><td>2.942</td><td>4.680</td><td>90.1</td><td>89.7</td><td>110.9</td><td>c/a' = 1.593</td><td></td></tr><tr><td>0.63</td><td>54</td><td>3.318</td><td>2.961</td><td>4.543</td><td>89.9</td><td>90.0</td><td>123.8</td><td>c/a' = 1.447</td><td></td></tr><tr><td>0.70</td><td>54</td><td>2.953</td><td>3.285</td><td>4.533</td><td>89.7</td><td>89.6</td><td>123.7</td><td>c/a' = 1.453</td><td></td></tr><tr><td>0.81</td><td>128</td><td>2.934</td><td>2.997</td><td>4.687</td><td>90.0</td><td>89.9</td><td>119.3</td><td>c/a' = 1.580</td><td></td></tr></tbody></table>

<table><caption>Table S2: Elastic constants $C_{ij}$ (GPa) and Zener anisotropy coefficient $A_z$ calculated for the bcc SQS supercell. The elastic constants of bcc Ti at high temperature are also given for comparison. The experimental EC with isotropic Zener Coefficient is given in parentheses for the known composition.</caption><thead><tr><th>x</th><th>0.0 ([61])</th><th>0.11</th><th>0.20 ([4])</th><th>0.33</th><th>0.40</th><th>1.0 ([59])</th></tr></thead><tbody><tr><td>C<sub>11</sub></td><td>162 (169)</td><td>155</td><td>156 (172±6)</td><td>154</td><td>147</td><td>70(98)</td></tr><tr><td>C<sub>12</sub></td><td>109 (127)</td><td>108</td><td>105 (108±1.5)</td><td>100</td><td>100</td><td>123(83)</td></tr><tr><td>C<sub>44</sub></td><td>42 (45)</td><td>46</td><td>39 (28±1.5)</td><td>45</td><td>41</td><td>39(38)</td></tr><tr><td>A<sub>z</sub></td><td>1.5 (2.1)</td><td>1.9</td><td>1.5 (0.87)</td><td>1.7</td><td>1.9</td><td>- (5.1)</td></tr></tbody></table>