
==== Front
Sci Adv
Sci Adv
sciadv
advances
Science Advances
2375-2548
American Association for the Advancement of Science

ado9697
10.1126/sciadv.ado9697
Research Article
Physical and Materials Sciences
SciAdv r-articles
Physical Sciences
Physical Sciences
Exceptional hardness in multiprincipal element alloys via hierarchical oxygen heterogeneities
Enhanced hardness by hierarchical heterogeneities
https://orcid.org/0000-0002-9929-4445
Beaudry David C. Conceptualization Data curation Formal analysis Investigation Methodology Supervision Validation Visualization Writing - original draft Writing - review & editing 1
https://orcid.org/0000-0001-6425-4331
Waters Michael J. Conceptualization Formal analysis Methodology Software Validation Visualization Writing - original draft Writing - review & editing 2
https://orcid.org/0000-0002-7350-8441
Valentino Gianna M. Formal analysis Investigation Visualization Writing - original draft Writing - review & editing 3
https://orcid.org/0009-0000-9185-3479
Foley Daniel L. Formal analysis Investigation Visualization 1
https://orcid.org/0000-0001-7808-7036
Anber Elaf Investigation 1
https://orcid.org/0000-0003-0448-3735
Rakita Yevgeny Formal analysis Visualization 1 4
Brandenburg Charlie J. Investigation 5
https://orcid.org/0000-0001-7786-397X
Couzinié Jean-Philippe Resources 6
https://orcid.org/0000-0002-8635-5069
Perrière Loïc Resources 6
https://orcid.org/0000-0001-6620-9390
Aoki Toshihiro Investigation 7
https://orcid.org/0000-0002-8388-7254
Knipling Keith E. Formal analysis Investigation Resources Visualization 8
https://orcid.org/0000-0002-6904-8171
Callahan Patrick G. Formal analysis Investigation Visualization 8
https://orcid.org/0000-0001-8863-4363
Redemann Benjamin W.Y. Formal analysis Investigation 1 9 10
https://orcid.org/0000-0002-8493-4630
McQueen Tyrel M. Data curation Formal analysis Funding acquisition Resources Supervision 1 9 10
https://orcid.org/0000-0001-5540-7084
Opila Elizabeth J. Conceptualization 5
https://orcid.org/0000-0003-0508-2175
Rondinelli James M. Conceptualization Formal analysis Funding acquisition Methodology Project administration Resources Software Supervision Validation Writing - original draft Writing - review & editing 2
https://orcid.org/0000-0001-5349-1411
Taheri Mitra L. Conceptualization Funding acquisition Methodology Project administration Resources Supervision Validation Writing - original draft Writing - review & editing 1 *
1 Department of Materials Science & Engineering, Johns Hopkins University, Baltimore, MD, USA.
2 Department of Materials Science & Engineering, Northwestern University, Evanston, IL, USA.
3 Department of Materials Science and Engineering, University of Maryland, College Park, MD, USA.
4 Department of Applied Physics and Applied Mathematics, Columbia University, New York, NY, USA.
5 Department of Materials Science & Engineering, University of Virginia, Charlottesville, VA, USA.
6 Univ Paris Est Creteil, CNRS, ICMPE, UMR 7182, 2 rue Henri Dunant, 94320 Thiais, France.
7 Irvine Materials Research Institute (IMRI), University of California, Irvine, Irvine, CA, USA.
8 Materials Science and Technology Division, U.S. Naval Research Laboratory, Washington, DC, USA.
9 Department of Chemistry, Johns Hopkins University, Baltimore, MD, USA.
10 William H. Miller III Department of Physics and Astronomy, Institute for Quantum Matter, Johns Hopkins University, Baltimore, MD, USA.
* Corresponding author. Email: mtaheri4@jhu.edu
20 9 2024
20 9 2024
10 38 eado969729 2 2024
15 8 2024
Copyright © 2024 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).
2024
The Authors
https://creativecommons.org/licenses/by-nc/4.0/ This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

Refractory multiprincipal element alloys (RMPEAs) are potential successors to incumbent high-temperature structural alloys, although efforts to improve oxidation resistance with large additions of passivating elements have led to embrittlement. RMPEAs containing group IV and V elements have a balance of properties including moderate ductility, low density, and the necessary formability. We find that oxidation of group IV-V RMPEAs induces hierarchical heterogeneities, ranging from nanoscale interstitial complexes to tertiary phases. This microstructural hierarchy considerably enhances hardness without indentation cracking, with values ranging between 12.1 and 22.6 GPa from the oxide-adjacent metal to the surface oxides, a 3.7 to 6.8× increase over the interstitial-free alloy. Our fundamental understanding of the oxygen influence on phase formation informs future alloy design to enhance oxidation resistance and obtain exceptional hardness while preserving plasticity.

Oxygen promotes heterogeneities in multiprincipal element alloys that bridge length scales to substantially enhance hardness.

http://dx.doi.org/10.13039/100000001 National Science Foundation DMR-1922234 http://dx.doi.org/10.13039/100000006 Office of Naval Research N00014-20-1-2368 http://dx.doi.org/10.13039/100000006 Office of Naval Research N00014-20-1-2368 http://dx.doi.org/10.13039/100000006 Office of Naval Research N00014-20-1-2368 http://dx.doi.org/10.13039/100000006 Office of Naval Research N00014-20-1-2368 http://dx.doi.org/10.13039/100000006 Office of Naval Research N00014-20-1-2368 http://dx.doi.org/10.13039/100000006 Office of Naval Research N00014-20-1-2368 http://dx.doi.org/10.13039/100000006 Office of Naval Research N00014-20-1-2368 http://dx.doi.org/10.13039/100023581 National Science Foundation Graduate Research Fellowship Program DGE2139757 http://dx.doi.org/10.13039/100023617 Naval Research Enterprise Internship Program http://dx.doi.org/10.13039/100000006 Office of Naval Research N00014-20-1-2368 http://dx.doi.org/10.13039/100000015 U.S. Department of Energy DE-AC02-05CH11231 http://dx.doi.org/10.13039/100000015 U.S. Department of Energy DE-AC02-06CH11357 http://dx.doi.org/10.13039/100009917 U.S. Naval Research Laboratory http://dx.doi.org/10.13039/100009917 U.S. Naval Research Laboratory http://dx.doi.org/10.13039/100000006 Office of Naval Research N00014-20-1-2368 http://dx.doi.org/10.13039/100000006 Office of Naval Research N00014-20-1-2368 http://dx.doi.org/10.13039/100000015 U.S. Department of Energy DE-AC02-05CH11231 http://dx.doi.org/10.13039/100000015 U.S. Department of Energy DE-AC02-06CH11357 http://dx.doi.org/10.13039/100000006 Office of Naval Research N00014-20-1-2368
==== Body
pmcINTRODUCTION

Refractory multiprincipal element alloys (RMPEAs), also called refractory high-entropy alloys or refractory complex concentrated alloys, consist of a combination of high melting point elements (Hf, Mo, Nb, Ta, Ti, V, W, and Zr) in near equimolar proportions, and promise strength and stability at high temperatures beyond that of incumbent dilute alloys (Ni- and Co-based superalloys and Nb-based alloys) (1–7). While the major motivation for developing this class of alloys is to enable hotter jet engines with cleaner emissions and higher efficiency, there has been increased interest in these alloys for nuclear fission, space travel, orthopedic implants, and other room temperature structural applications (1, 4, 8–11).

Arguably, the greatest limitation for RMPEAs is harmful internal oxidation without protective surface passivation at both transient and maximum operating temperatures. Understanding and predicting the phase evolution during oxidation in RMPEAs, especially those with group IV and V elements, have proven difficult due to the ambiguity of how elements with vastly different oxygen solubilities [i.e., Nb, 2.9 atomic % (at %) O; α-Ti, 33 at % O; α-Zr, 30 at % O at 1050°C (12–15)] and diffusivities behave when combined in a nondilute, heterogeneous chemical landscape. Most RMPEA oxidation studies have focused on the long-term phase formation without providing an understanding of initial phase evolution that affects long-term behavior. Attempts to alleviate the poor oxidation resistance have relied on adding conventional passivating elements (Al, Cr, and Si) in near-equimolar concentrations, which reduces ductility by the formation of deleterious intermetallic phases while also lowering the melting point and potential service temperatures (16–20). Alternative pathways are needed to control oxide formation and improve oxidation resistance while preserving ductility in RMPEAs.

Here, we uncover the fundamental mechanisms driving the initial transformations in oxidation of group IV-V RMPEAs at elevated temperatures that inform hypotheses on techniques to advance alloy development for long-term oxidation resistance. Equimolar NbTiZr is selected as a model system to examine phase stability, oxide formation, and microstructure at elevated temperatures because of its low density of 6.63 g/cm3, high strength and good ductility, and the presence of its comprising elements in half of all RMPEA compositions studied (2–5). Using a combination of high-resolution characterization and first-principles–based calculations, we determine how complex hierarchical microstructures evolve temporally and spatially during isothermal oxidation of NbTiZr at 1050°C in atmospheric air. The microstructural evolution presented in this work, which commences with oxygen-induced isothermal spinodal decomposition, results in an exceptional gradient hardness enhancement of up to 6.8× over the interstitial-free alloy while preserving plasticity. We confirmed the extension of this hardening phenomenon to other alloys and kinetic regimes by producing a comparable hardness increase in equimolar HfNbTiZr through a similar hierarchical phase separation that proceeded by nucleation and growth. These alloys were selected because the equimolar compositions consisting of a mix of group IV and V elements enable investigation into how varied metal-nonmetal solubilities and affinities of constituents manifest in a chemically complex system. This approach provides for broad applicability of these findings to other nondilute systems with comparable disparities. Although the hierarchical microstructures in the model systems reported here do not alone provide exceptional oxidation resistance, this study provides fundamental insight on initial phase transformations and identifies opportunities to alter the heterogeneities through oxygen-induced phase separations to concentrate alloying additions without hindering the ductility of the alloy. This may be an attractive alternative to the typical approach of adding embrittling levels of passivating elements to achieve enhanced oxidation resistance with disregard for consequences on the mechanical properties.

RESULTS

Order and miscibility in NbTiZr

To explore the effects of oxygen on chemical short-range order (SRO), we established a baseline by first determining the influence of thermal annealing on SRO in the unoxidized single-phase (fig. S1 and table S1) alloy. We began with a synchrotron x-ray total scattering experiment followed by a pair distribution function (PDF) analysis of the “homogenized” (1250°C, ~0.75 homologous temperature, TH) and homogenized + aged (1050°C, ~0.6 TH) conditions, which represent the dilute limit of oxygen from processing impurities. The relative increase in Zr-Zr first nearest-neighbor coordination at the expense of X-X pairs (X = Nb,Ti) for the 1050°C condition (Fig. 1A) indicates a tendency for SRO chemical clustering in the overall long-range disordered body-centered cubic (BCC) structure. Since this measurement was at the limit of resolution, we confirmed the qualitative PDF results with theoretical Monte Carlo (MC) simulations, using cluster expansion models fit to ab initio results, from which the SRO Warren-Cowley (WC) parameters were extracted at a temperature range of 700° to 1250°C (Fig. 1B). Between 1250° and 1050°C, the simulation of an oxygen-free alloy agrees with the experimental PDF results, showing a minute increase in the tendency toward Zr-Zr clustering. The simulation was repeated with 2 at % O included (Fig. 1B) and gave a stronger increase in the magnitude of the WC parameters with decreased temperature, suggesting that SRO in the alloy is sensitive to the presence of oxygen.

Fig. 1. Thermal and interstitial effects on short-range order and miscibility.

(A) PDF analysis of the change in short-range order between 1250° and 1050°C, with the gray line denoting the difference in the curves, X = Nb,Ti. (B) MC-derived WC parameters that confirm minor SRO increase between 1250° and 1050°C for 0 at % O and a more pronounced increase at 2 at % O. (C) MC-derived WC parameters for varying oxygen concentrations at 1050°C that show a drastic increase in SRO with increasing oxygen content. (D) Concentration versus distance from the MC simulations for NbTiZr +2 at % O at 1050°C, qualitatively matching the observed oxygen complexes. (E) APT reconstruction with oxygen complexes shown by the purple surface designating an iso-concentration of 8.5 at % O. (F) Concentration versus distance proximity histogram from the APT 8.5 at % O iso-concentration surfaces averaged for all interfaces, with SD error bars displayed. Positive distance values correspond to the interior of the complexes. (G) Quaternary isotherm at 1050°C of binodal tielines (spheres connected by solid lines) derived from the equilibrium simulations and lower spinodal points (tetrahedra connected by dashed line) from the paraequilibrium simulations. The gray vertical line demonstrates the equimolar NbTiZr undergoing isothermal oxygen uptake and crossing the lower spinodal line.

To determine the influence of dissolved interstitial oxygen on SRO and phase stability, MC simulations were performed on simulated alloys containing 0 to 20 at % O at 1050°C (Fig. 1C). The WC parameters extracted from the simulated data predict that metal-metal SRO is negligible up to 2 at % O, beyond which it increases rapidly with primarily Zr-Zr and Nb-Nb clustering. The cation-oxygen WC parameters indicate a strong preference for Zr-O coordination and reduced preference for Ti-O and Nb-O coordination between 5 and 10 at % O. These simulated WC parameters are qualitatively confirmed by atom-probe tomography (APT) of the single-phase base metal adjacent to the oxygen diffusion zone (ODZ), which reveals subnanometer oxygen-rich domains (Fig. 1E) that are enriched to an average of 13.9 at % O. The bulk concentration in the region is 6.7 at % O. Normalizing concentrations for only metals revealed that the composition changes in the domains were +4.8 at % Zr, −2.1 at % Ti, and −2.7 at % Nb relative to the bulk concentration (Fig. 1F). The domain composition matches that predicted by the equilibrium MC simulations with bulk-averaged 2 at % O (Fig. 1D). Further paraequilibrium (21, 22) MC simulations, which are detailed below and in Materials and Methods, account for kinetic variations between experiment and equilibrium simulations and find the corresponding spinodal upper limits to be 8 to 9 at % O (Fig. 1G, fig. S2, and table S2), which agrees well with the experimentally observed solubility limit at 6.7% O. Electron diffraction of the oxygen-enriched base metal revealed a single-phase BCC structure with evident diffuse scattering (Fig. 2E) that suggests fine real-space features without a structural phase transformation. This is comparable to previous reports that in BCC-ω systems the source of diffuse scattering is short-range atomic displacements along the long-range displacive transformation directions (23, 24). The driving force for long-range ω-type displacements in this alloy class is discussed in the sections to follow, which along with these results suggest that chemical clustering caused by initial oxygen uptake may also induce local ω-type deviations from BCC structural symmetry in the nanoscale complexes. SRO and nanoscale complexes are often precursors to phase formation, and in this case precede spinodal decomposition.

Fig. 2. Temporal evolution of phase transformations at the base metal interface.

(A) SEM-BSE images of a cross-sectioned oxidized sample showing pockets of spinodal decomposition extending into the base metal as observed after 10 min of oxidation at 1050°C, with an inset (yellow indicating highest intensity and black indicating lowest) showing the spinodal decomposition consuming a remaining pocket of single-phase base metal. The sample surface lies out of view and to the right. These micrographs were taken from the 30-min sample in a remaining uncoarsened region identical to the dominant microstructure observed at 10 min. (B) SEM-BSE images showing the onset of coarsening at the base metal interface after 30 min of oxidation at 1050°C, and an inset [colorized as in (A)] of the coarsening lamellae consuming a pocket of spinodal decomposition. (C) SEM-BSE image showing full coarsening into lamellae adjacent to the base metal after 3 hours of oxidation at 1050°C. (D) STEM-EDS map from a pocket of spinodal decomposition with clear Nb-Zr segregation across a diffuse interface. (E) TEM selected area diffraction of the base metal (lattice constant of 3.421 Å) and the spinodal phases (lattice constants of 3.43 and 3.33 Å) along the [ 1¯33 ] zone axis. Yellow arrows indicate arcs of diffuse intensity in the base metal diffraction. Cropped and magnified images of the reflections show clear peak splitting in the spinodal phases. (F) APT reconstruction of the coarsened lamellae with a representative location denoted by the red outline shown in (C). (G) STEM-EDS composition line profile from arrow shown in (D) with an integrated width of 16 nm. (H) Concentration versus distance proximity histogram of a 30 at % Nb iso-concentration surface from the reconstruction in F with SD error bars displayed. (I) Concentration versus distance from the MC simulations for NbTiZr +10 at % O at 1050°C, matching well with the APT proximity histogram.

During oxidation, the single-phase base metal initially undergoes spinodal decomposition followed by coarsening of the spinodal structure into aligned lamellae (Fig. 2). The supercell from our MC equilibrium model (Fig. 2I) shows that the introduction of oxygen to equimolar NbTiZr results in increased Nb-(Zr,O) segregation as the oxygen concentration approaches 10 at %. This was matched closely by scanning transmission electron microscopy–energy-dispersive x-ray spectroscopy (STEM-EDS) of a diffuse interface in the spinodal region (Fig. 2, D and G) and APT results of the lamellar region in the 3-hour microstructure (Fig. 2, F and H). We corroborate these results further by performing paraequilibrium (21, 22) MC simulations. The spinodal points, which delineate regions where the alloy becomes thermodynamically unstable with the free energy curve concave down, were found from the computed second derivative of free energy. The spinodal points plotted along the binodal tie lines in the quaternary 1050°C isotherm (Fig. 1G and fig. S2) confirm that the uptake of oxygen isothermally pushes the solid solution alloy into a thermodynamically unstable region where the predicted threshold is ~10 at % O (table S2). This contrasts with most reports of spinodal decomposition in which the variable that shifts the alloy to the unstable regime is a temperature change during quenching rather than composition change during oxidation. TEM diffraction of the spinodal phases (Fig. 2E) shows a decrease in diffuse scattering relative to the base metal. This corroborates that the (Zr,O) nanoclusters were the source of diffuse scattering in the base metal before spinodal decomposition when they separated into a disordered BCC phase with a distinct lattice parameter and diffraction peaks.

The induced miscibility gap and spinodal decomposition could be manipulated across temperature ranges and into varied composition space based on the findings reported here and established thermodynamic principles. The positive mixing enthalpy (ΔHmix) for Nb-Zr [4 kJ/mol (25)] is the primary condition for production of the oxygen-induced high-temperature miscibility gap, which is manifested in the overall alloy ΔHmix increase with oxygen in the paraequilibrium simulations (fig. S2). A more positive ΔHmix typically increases the temperature of the critical point (Tc) of a miscibility gap (26, 27). For exceptions to this trend, such as Nb-Zr [ΔHmix = 4 kJ/mol (25) and Tc = 977°C (5)] and Ta-Zr [ΔHmix = 3 kJ/mol (25) and Tc = 1783°C (5)], a large electronegativity difference (Δχ) may offer a better indicator that the solid solution is stable (2, 3, 11, 28). For identification of cation additions that raise Tc of the interstitial-free miscibility gap, thereby affecting the oxygen-induced gap, we suggest a selection method of initial elemental screening by empirical parameters [ΔHmix, Δχ, and Ω (29)] followed by CALPHAD for confirmation. This method can also predict which alloying additions of targeted atomic radii would segregate to each phase, which allows for control of the coherent spinodal through tuning of the resulting misfit strain and elastic moduli (30). While these cationic interactions are essential to invoking a high-temperature miscibility gap with interstitial incorporation, they must be viewed as the initial conditions, which are further modified by the cation-anion interactions.

The secondary condition to invoke an oxygen-induced high-temperature miscibility gap is varied cation-oxygen interaction parameters. Oxygen solubility, enthalpies of solution, electron affinities, and atomic radii all have important contributions to the effect of oxygen in these systems (31) and should be considered when modeling predicted behavior. However, we propose a more rapid screening parameter of chemical potentials of oxygen extracted from the phases on the oxygen-lean portion of the cation-oxygen convex hulls. These values for the Nb (−4.31 eV), Ti (−6.53 eV for α and −7.33 eV for β), and Zr (−5.92 eV for α and −6.24 for β) convex hulls with oxygen (32) show substantial variation between the main segregating species of Nb and Zr. The slightly larger chemical potential for Ti is likely offset by the intermediate Ti-Nb ΔHmix of 2 kJ/mol (25) allowing for solution of these two species in the lower oxygen phase with less influence of oxygen interactions. A more high-fidelity predictive approach is available in replication of our paraequilibrium MC simulations (Fig. 1G and fig. S2). One candidate alloying element for raising Tc through oxygen interactions is scandium, with binary mixing enthalpies (Sc-Nb 18 kJ/mol and Sc-Zr 4 kJ/mol) (25) and an oxygen chemical potential (−6.35 eV for α and −6.53 for β) (32) that suggest that it will segregate to the Zr phase. The large magnitudes of Sc-Nb repulsive forces and Sc-O attractive forces could increase the temperature at which the solid solution becomes unstable, elevating the Tc of the oxygen-induced miscibility gap. These guidelines follow the same thermodynamic principles that can apply to other alloys and other anions such as nitrogen, carbon, and boron. It should be noted that the chemical potentials of nitrogen follow similar disparities for group IV-V elements (Nb, −2.88 eV; Ti, −4.26 eV; and Zr, −3.9 eV) (32). However, the experimentally measured compositions and the simulations with only metals and oxygen prove that the influence of cation-nitrogen interactions is negligible in this particular composition space when oxygen is present in high fluxes. This can be attributed to roughly 50% larger oxygen chemical potentials than that of nitrogen, which highlights the need to consider competing nonmetal interactions when expanding this design strategy. While reparameterization of the cluster model used in the MC simulations for new elements is resource intensive, it would be prudent for future investigation.

Binodal phase separation and spinodal decomposition have not been observed in NbTiZr in previous studies due to kinetic limitations at the lower temperatures within the oxygen-free spinodal and binodal regions (2, 3, 11, 33). We find that spinodal decomposition during oxidation of NbTiZr at 1050°C occurs for the following reasons: (i) The cations have positive binary mixing enthalpies and varied oxygen affinities; (ii) interstitial oxygen alters thermodynamics to produce a high-temperature miscibility gap; (iii) oxygen clusters into nanoscale complexes with Zr, inducing a higher amplitude of compositional fluctuation and a larger driving force for spinodal decomposition; and (iv) the high-temperature miscibility gap overcomes the diffusional kinetic hindrance to spinodal decomposition reported in the O-free alloy. The interstitial complexes and the following spinodal decomposition lay the kinetic and thermodynamic foundation for the hierarchy of phase evolutions with further oxygen uptake that dictate mechanical and oxidation properties.

Oxygen-induced exceptional hardness

The temporal series of isothermal oxidation of NbTiZr (Fig. 2) provides insight into both hardness and phase evolution while uncovering a mechanism for controlling the surface-directed gradient in hardness. The following five distinct microstructural transformations occurred during oxidation that give rise to a hierarchical microstructure: Transformation 1, formation of interstitial complexes in the base alloy; Transformation 2, rapid spinodal decomposition of the single-phase base alloy; Transformation 3, separation of a third phase upon increased oxygen content; Transformation 4, conversion from metallic phases with dissolved oxygen to suboxides and monoxides; and Transformation 5, coarsening of the spinodal structure into lamellae, which includes the displacive BCC-to-ω transformation and formation of an internal TixOy scale. While further transformations will occur to complete the oxidation, our focus is the initial stages in these subsurface regions to determine previously unresolved phase transformations in the early stages of oxidation of group IV-V RMPEAs.

The nanoindentation hardness (Fig. 3) of the single-phase base metal with interstitial complexes had averages ranging between 10.1 ± 0.1 GPa 70 μm away from the spinodal interface and 12.1 ± 0.2 GPa at the interface, with 6.7 at % O per APT. These measured values display a 3.7× increase in hardness compared to that of the interstitial-free NbTiZr (2). After 10 min, the ODZ had clear spinodal decomposition (Transformation 2) near the base metal interface (Fig. 2) and a fine microstructure of oxides (Transformations 3 and 4) that lost bi-continuity near the surface. The spinodal microstructure displayed a hardness of 18.6 ± 0.6 GPa [a 5.6× increase over interstitial-free NbTiZr (2)], while the oxides just beneath the surface averaged 22.6 ± 0.4 GPa [a 6.8× increase over interstitial-free NbTiZr (2)] with a maximum measurement of 23.1 GPa. After 30 min, the spinodal structure at the base metal interface mostly coarsened into a narrow band of lamellae (Transformation 5) with a few locations still remaining uncoarsened. This coarsening consumed the entire base metal interface and extended further toward the surface to consume the inner subsurface monoxides after 3 hours of oxidation. The hardness decreased after 3 hours to 7.6 to 12.4 GPa due to the coarsening of the spinodal structure into lamellae. However, the elastic modulus remained relatively unchanged and ranged between 103 and 153 GPa across the layers (Fig. 3C).

Fig. 3. Evolution of mechanical properties with microstructural variations across oxidation times.

Nanoindentation plots of average elastic modulus and hardness for each x distance are displayed in the corresponding 2D heat maps. The SD for each x distance is shown as a shaded error region in the plots. Vertical lines represent microstructural transitions. Nanoindentation was performed on representative microstructures observed after (A) 10 min of oxidation at 1050°C for NbTiZr, collected from an identical uncoarsened region in the 30 min sample; (B) 30 min of oxidation at 1050°C for NbTiZr, collected from a coarsened region that does not contain the sample surface to better display the coarsening transition; (C) 3 hours of oxidation at 1050°C for NbTiZr, with black pixels indicating omitted data per the outlier screening described in Materials and Methods; (D) 30 min of oxidation at 1050°C for HfNbTiZr with a maximum value of 21.4 GPa. (A) and (C) extend almost exactly to the surface, while (B) and (D) are limited to a subsurface range focusing on the coarsening front.

In addition to the 6.8× increase in hardness presented in Fig. 3, high-magnification micrographs revealed no evidence of cracking (fig. S3) in any of the indented layers, which indicates moderate toughness. This is further evidenced by the plasticity observed via shear banding in the pileup region adjacent to the base metal indents (fig. S3). The indents in the oxides near the surface also show no cracking in or around the residual indent impression. Instead, homogeneous deformation is observed inside the indent impressions despite the fact that clear microstructural heterogeneities are present (fig. S3). In between is the spinodal region, marking the transition between microscale plastic instabilities (i.e., shear banding) and nanoscale strain accommodations (i.e., introduction of a high density of interfaces to delocalize plasticity). These observations provide a qualitative assessment of toughness and suggest that the hierarchical phase evolution in Transformations 1 to 4 results in nanoscale heterogeneities that accommodate strain and mitigate cracking in what would normally be brittle phases.

The 3.7× increase in hardness of the base metal is directly attributed to the oxygen nanodomains in Transformation 1 (Fig. 1). Nitrogen-rich domains were also found in this region, with an average size on the order of a nanometer (fig. S4). The bulk average nitrogen concentration was 0.1 at %, while the nitrogen domains contained 4.6 at % and the normalized cation levels varied by less than 1 at % from the bulk average. For comparison, the base metal hardness is more than 6× higher than standard refractory alloys, like Nb C-103 with a hardness of 1.5 GPa at ~1.5 at % oxygen (34). While our observed increase in hardness is consistent with interstitial bulk-dopant effects reported in other RMPEAs with lower levels of interstitial doping (8, 35–37), the hardness values measured in this study exceed the previously reported values for group IV-V alloys with oxygen or nitrogen additions (37). Our findings suggest that the exceptional hardness and absence of cracking are a result of a synergistic effect from the simultaneous presence of oxygen and nitrogen leading to hierarchical microstructure. In addition to the base metal, our measured NbTiZr oxide hardness of 22.6 ± 0.4 GPa has a 1.3 to 2.3× increase over common oxides and composites such as monoclinic ZrO2 [9.8 to 13 GPa (38)], cemented tungsten carbide [11.7 GPa (39)], and sapphire [17.4 GPa (40)]. We confirmed the broader applicability of this microstructural evolution and hardening behavior by oxidation and subsequent hardness measurements for equimolar HfNbTiZr (Fig. 3D and fig. S5). The hierarchical phase separation and maximum measured hardness of 21.4 GPa confims the broader applicability to compositions beyond equimolar NbTiZr.

Effect of oxygen on phase formation

The genesis of the interstitial complexes and spinodal decomposition were resolved by initial observations and models, but the mechanism behind the coarsening required further characterization and modeling. The Nb-rich, O-lean lamella exhibits the disordered BCC structure and the (Zr,O)-rich lamella is the hexagonal ω phase with some ordering (Fig. 4B and fig. S6), which is well studied in group IV dilute alloys (23, 24, 41, 42). Geometric phase analysis (GPA) performed on a high-angle annular dark-field (HAADF)–STEM image of the BCC-ω interface (Fig. 4, C and D) reveals a misfit dislocation spacing ranging from 10 to 14 planes, corresponding to a semi-coherent misfit strain that varied between 7.1 and 10%. This diffusional ω phase formation proceeds by collapse of BCC {111} planes and is favorable when β-stabilizers diffuse away from a location in group IV β alloys (41, 42), as occurs here with Nb due to the miscibility gap. The BCC-ω transformation and the resulting reduction in coherency cause the spinodal structure to coarsen into lamellae to reduce the total interfacial energy. The thermodynamic driving force for this planar collapse is oxygen concentration, as shown by our density functional theory (DFT)–calculated formation energy difference versus oxygen concentration for Zr2TiO (Fig. 4E). Despite this simplified model of the observed ω-phase composition in region II (45.0Zr-29.7Ti-21.8O at %), these calculations highlight how drastically oxygen enhances the BCC-to-ω collapse and ultimately causes the coarsening of the spinodal structure. The discovery that oxygen flux and the ω transformation promotes the coarsening reaction allows for control of the rate and degree of coarsening by selective doping to alter the interfacial strain. The reduction of interfacial misfit, and thereby the driving force for the coarsening reaction, could be achieved through the aforementioned methods for predicting to which phase a cation will segregate in combination with choosing a dopant with a target ionic radius. This could be accomplished through dopant segregation to either the BCC or the ω-phase, which would effectively control the desired gradient of mechanical properties and the density of semi-coherent interfaces through which oxygen and cations diffuse.

Fig. 4. Driving force for ω transformation and resultant coarsening.

(A) SEM-BSE, with yellow indicating highest intensity and black the lowest, showing the coarsening front in the 30-min sample and a representative area of phases captured in (B). The emergence of a (Ti,O)-rich scale at this front is evident. (B) STEM-HAADF of the two-phase lamellar interface from a location in the 3-hour sample equal distance from the base metal as the location shown in (A), with corresponding fast Fourier transforms of the Nb-rich BCC and Zr-rich ω phases. (C) Colorized inverse fast Fourier transform of a larger view field of (B) that shows misfit dislocation spacing. Masked reflections are (11¯0)β and (1¯101)ω. (D) GPA relative εxx strain map of the BCC-ω interface showing highly localized strain fields at the misfit dislocations. (E) Plot of the formation energy difference per formula unit of ω minus BCC, or the driving force for planar collapse, for varying oxygen molar fractions of Zr2TiOx.

The hierarchical microstructure after 3 hours of oxidation (Fig. 5) has clearly defined regions resulting from the transformation events that succeeded the spinodal decomposition. All three oxidation exposures contained a region near the spinodal/lamellae interface with the base metal where the Ti-rich phase emerged. After 3 hours, this third Ti phase presents as a ~10 nm sublamella that matches the hexagonal ω-type structure (Fig. 5C) and has a composition of 58.7Ti-24.2O-15.4H-0.8Nb-0.4Zr at %, which is the highest hydrogen content of any phase measured. APT measurements of H content should be viewed as relative between phases since absolute measurements are inflated by analysis chamber H2 condensation. While H has been measured in a suboxide and N has been measured in the base metal, the MC simulations and experimental results do not suggest that they drive the hierarchical phase separations as oxygen provides the necessary driving force alone. However, the effects of H and N on the oxide polymorph stabilities should not be discounted in this composition space with a large number of metastable structures. APT analysis of the Zr sublamella from region III in the 3-hour sample reveals a composition of 31.5Zr-12.6Ti-47.3O-6.0H-1.5N-0.1Nb at %. These results confirm that further oxygen enrichment of the initial spinodal decomposition caused a second phase separation between Zr and Ti, which is consistent with our MC simulations (fig. S7) that show the starting region II ω composition separates into Zr- and Ti-rich phases as oxygen levels increase. STEM-HAADF and bright-field (BF) analysis of the nanotwinned Zr monoxide in region III determined that this could be an orthorhombic suboxide of (Zr,Ti)O2 Srilankite (Pbcn) (43), polar ZrO2 (Pca21) (44), OI-ZrO2 (Pbca) (45), or a monoclinic suboxide of Baddeleyite ZrO2 (P21/c) (45) (Fig. 5C and fig. S8). These nanoscale heterogeneities in a complex energetic landscape have important implications for both dislocation motion and the diffusivity of oxygen into the bulk (36, 46).

Fig. 5. Hierarchical structure after 3 hours of oxidation displayed across length scales.

(A) SEM-BSE of microstructure extending from the base metal to the surface. Region I retains the single-phase BCC base metal with an oxygen concentration of 6.7 at % and a high density of nitrogen and oxygen interstitial complexes, as seen in Fig. 1 and fig. S4. Region II exhibits a two-phase lamellar layer with the Nb-rich A2 phase and a (Zr,O)-rich ω phase (as shown in Fig. 4), with Ti relatively even between them. Region III shows a three-phase lamellar structure due to a secondary phase separation within the ω lamella of region II into ~10-nm sublamellae with Zr/Ti segregation (Transformation 3). Region IV contains a finely-dispersed microstructure with monoxide compositions. (B) Left to right: STEM-EELS of region II (color legend for all chemical mapping displayed), STEM-EDS of region III with STEM-FFT of Nb-rich BCC and representative location, STEM-HAADF of the region III-IV interface (coarsening front, marked by white line) with inset FFT of the ω-type TixOy internal scale. (C) Colorized STEM-HAADF (adjusted gamma, radial difference filter applied) from region III [location given by black outline in (B)] showing nanotwins in the Zr sublamella with Zr-rich FFT from location shown and Ti-ω FFT from a representative location. APT reconstruction of the Zr/Ti sublamellae in region III with representative location shown in the black outline in (B) and corresponding proximity histogram for a 28 at % Zr iso-concentration surface (standard deviation shaded). (D) APT reconstruction of the monoxides in region IV with representative location shown by the red outline in (B), and a corresponding 1D concentration profile (SD shaded) along the length of 8.1-nm-diameter cylinder that passes through all three phases as shown in the inset. Iso-concentration surfaces are denoted by the color key and are 3.5Nb, 30Ti, 7Zr (at %).

APT from region IV near the interface with region III shows three fine monoxides (Fig. 5D). While the alloy does not undergo selective oxidation, there is a hierarchical sequence of oxidation due to local oxygen chemical potential in the microstructural layers and saturation of the more favorable suboxide. Region II sees selective oxidation of the Zr-Ti phase; region III sees saturation of the separated Zr phase and the onset of selective oxygen increase in the separated Ti phase; and region IV sees saturation of the Ti phase and oxidation of the Nb-phase into a monoxide. The NbO monoxide is only reported to have a stable polymorph of the rocksalt structure (12), which is inherently defective with both cation and anion vacancies. An ω-type TixOy suboxide emerges as an internal scale parallel to the surface along the advancing lamellar coarsening front (Fig. 5B and fig. S9), which also has inherent vacancies (47). The presence of anion vacancies increases the oxygen diffusivities and allows for increased transport of oxygen into the base metal. It is well established that doping oxides with cations of a higher charge state reduces the vacancy concentration to maintain charge neutrality (48). We propose the potential for using the oxygen-induced miscibility gap to concentrate minor amounts of dopants to the phases of interest.

Molybdenum is a potential candidate for NbO doping due to a higher expected charge state and a lower metal-oxygen binding dissociation energy than Nb, which would increase the energy barrier to oxygen vacancy formation (48, 49). Molybdenum has been reported to segregate to the Nb-rich BCC phase when added to oxygen-free NbTiZr (3) which, along with a relatively small oxygen chemical potential (−3.23 eV) (32), suggests that it would segregate to NbO rather than the Ti- or Zr-rich monoxides. The trivalent elements Sc and Y are potential candidates for reducing oxygen vacancies in the TixOy scale due to the higher charge state, large oxygen chemical potentials of −6.35 eV for both elements (32), and positive mixing enthalpies with Nb (25). The reduction of oxygen diffusion across this scale offers a unique opportunity to provide internal “passivation” if properly doped through the oxygen-induced miscibility gap due to its surface-parallel orientation.

DISCUSSION

Findings on NbTiZr and HfNbTiZr illustrate that oxidized RMPEAs with a mix of group IV and V elements form a hierarchy of layered heterogeneities that results in tunable high hardness without sacrificing plasticity. This unique microstructure-property evolution, enabled by nanoscale interstitial complexes, likely extends to other MPEA composition spaces with moderate binary cation immiscibilities and comparable disparities in oxygen solubility and affinity. Group IV elements have several suboxides on the M-O convex hull with many metastable polymorphs lying slightly above the hull. These phases precede the most stable oxide, have lower oxygen chemical potentials, and must be formed first during increase of oxygen content. For multicomponent systems, all the suboxides compete against a number of phases that are influenced by interfacial strains and solutioned cations, leading to the stabilization of previously unreported suboxides as we found here in the Zr-rich lamellar suboxide (Fig. 5C and fig. S8). The coupling of experiments and modeling provided an unprecedented in-depth analysis and mechanistic determination in this work that demonstrates previously unobtainable insight into interstitial-driven alloy design. The mechanistic understanding of initial subsurface phase transformations during oxidation of group IV-V RMPEAs established in this study lays the foundation for alloy development with more targeted control of phase formation in long-term oxidation studies. This could be achieved by targeted segregation of minor alloying additions without substantially affecting the base alloy composition and properties, which would be facilitated by the hierarchical phase separations and underlying mechanisms uncovered in this work.

MATERIALS AND METHODS

Alloy processing

The nominally 33.33Nb-33.33Ti-33.33Zr (at %) alloy was cast in an arc-melting furnace under argon atmosphere using high purity elemental Nb, Ti, and Zr. A water-cooled copper hearth was used in which the mold shape was a cross section of a cylinder with hemispherical ends. The ingot was flipped and re-melted five times to ensure the elements were distributed homogeneously. The ingot was then encapsulated in fused quartz under a partial Ar atmosphere and given a homogenization heat treatment at 1250°C for 24 hours, followed by an ice water quench in which the quartz was broken immediately upon entering the water. After homogenization, an ingot section was encapsulated, treated at 1050°C for 24 hours, and then quenched to investigate structure and order by x-ray analysis. A 5-mm cube for the bulk oxidation experiment was sectioned from the homogenized ingot, polished to 1200-grit SiC sandpaper, and oxidized in atmospheric air at 1050°C for 3 hours. The sample was removed from the furnace and “quenched” on a large aluminum heat sink held at room temperature. Oxidation was repeated with an ice water quench following 3 hours of oxidation to qualitatively confirm that the microstructure and extent of coarsening was not dependent on quenching rate. This process was repeated for 10 and 30 min of oxidation with the ice water quench. All cubes were then cross-sectioned and polished to 0.05-μm colloidal silica for scanning electron microscope (SEM) analysis and preparation of both TEM and APT lift-out samples.

Nanoindentation

Instrumented nanoindentation and the Oliver-Pharr method (50) were used to measure the elastic modulus and hardness values across the ODZ. The NanoBlitz 3D method was used on a KLA iNano system, outfitted with a diamond Berkovich tip, to generate more than 2800 indents in the specified area with a 5-μm spacing. A load of 50 mN (corresponding to 415- to 820-nm indentation depths as shown in fig. S10) and a Poisson’s ratio of 0.3 were used for all measurements. x and y dimensions (and indent counts) of the 5-μm spacing indentation areas for the NbTiZr “10-min”, 30-min, and 3-hour, and HfNbTiZr microstructures were 400 × 100 μm2 (80 × 20), 300 × 60 μm2 (60 × 12), 800 × 100 μm2 (160 × 20), and 320 × 100 μm2 (64 × 20), respectively. The 10-min microstructure area was gathered from an uncoarsened region in the 30-min sample due to identical microstructure of the small fraction of remaining spinodal decomposition areas in the 30-min sample. Hardness and modulus values were averaged across the direction parallel to the layer interfaces and perpendicular to the surface, such that 1D hardness and modulus plots can be examined across interfacial boundaries. The NanoBlitz method was repeated for a 10-μm indent spacing to ensure that the measured values did not change (fig. S3), confirming that there were no overlapping interaction volumes. For the NbTiZr 3-hour sample, datapoints with a modulus or hardness that deviated more than 30% from the average of the given x value (x axis is perpendicular to the surface) were omitted (and colored black in contour plots) to better represent the progression of mechanical properties across the layers without convolution from unparallel interfaces, uneven topography, surface contamination, or other confounding variables.

Microscopy and microanalysis

X-ray total scattering measurements were taken from ~1-mm-thick sections of the ingot. Both samples were homogenized, and one of them was aged, as described above but for 18 hours rather than 24 hours. The measurements were performed at the NSLS-ii synchrotron facility at the 28-ID-1 beamline using a Perkin Elmer detector. Since a PDF is a projection of the interatomic atom-atom distance weighted by the scattering form-factor, it is a good probe for the average short-, medium- and long-range order. When focusing on low r-regions and comparing two related PDFs, one can extract changes in the local chemical environment. The measurements were done in a PDF mode for total scattering, where the “camera length” of the detector was 21.2 cm, and the beamline radiation was 0.1665 Å. Calibration was done using Ni powder. The total scattering patterns were then azimuthally integrated using “pyFAI” package (51) and reduced and Fourier transformed to a PDF using “pdfgetx3” (52). The structures of the unoxidized homogenized and aged samples were analyzed by x-ray diffraction at room temperature using a Panalytical Empyrean X-ray diffractometer. The samples were measured from 20° to 90° and 30° to 90° 2θ using a Cu radiation source (45 kV and 40 mA). The diffraction patterns were refined using TOPAS5 Rietveld refinement (53) software (fig. S1, tables S1 and S3, and equation S1).

The cross-sectioned and polished samples were imaged in a Thermo Fisher Scientific Helios G4 UC SEM. Regions of interest were targeted for TEM lift-outs using a Ga focused ion beam in the same Helios microscope. Beam energies of 30, 8, and 5 kV were used during the thinning process. Extended low-energy exposures followed at 2 and 1 kV to minimize the Ga damage layer. TEM analysis, including STEM and electron energy-loss spectroscopy (EELS), of region II was performed in a Thermo Fisher Scientific Spectra 300 kV fitted with an X-CFEG cold field-emission source, panther STEM detection system, GIF Continuum 1066 EELS detector, Super-X symmetrical quad-crystal EDS system, and fifth-order-corrected S-corr Cs probe corrector. STEM HAADF/BF imaging and EDS of regions III and IV were conducted by the double-aberration-corrected JEOL ARM 300CF, equipped with dual 100-mm2 silicon drift detectors operating at 300 kV. The convergence semi-angle for STEM imaging and EDS was 25.6 mrad and the collection angles for STEM HAADF and BF imaging were 54 to 220 and ~14.2 mrad, respectively. S/TEM analysis of the spinodal structure was performed on a single-aberration-corrected JEOL GrandARM2 at an accelerating voltage of 300 kV. Diffraction patterns were collected on a JEOL OneView detector with a camera length of 500 mm. STEM-EDS was collected with JEOL dual-EDS detectors. All EDS analysis was standardless and thereby not considered perfectly precise in the quantitative results displayed. More reliable quantitative compositions can be found in the APT results. High-resolution STEM images from region III were filtered with the radial difference filter produced by HREM Research Inc.

Further high-resolution analysis was done by APT, using a Cameca 4000X Si local electrode atom probe with a 355-nm ultraviolet pulsed laser. Specimens for APT were prepared using standard lift-out and milling procedures (54, 55) using a Ga ion beam in an FEI Nova 600, as well as in a Thermo Fisher Scientific Helios G4 UC. The conditions for the base metal (region I) run were a vacuum level of ~1 × 10−11 torr, 40 K specimen base temperature, 60 pJ nominal laser pulse energy, a pulse repetition rate of 250 kHz, and a detection rate of 0.01 ions per pulse (1%). The conditions for the two-phase lamellar (region II) run were a vacuum level of ~2 × 10−11 torr, 40 K specimen base temperature, 40 pJ nominal laser pulse energy, a pulse repetition rate of 250 kHz, and a detection rate of 0.01 ions per pulse (1%). The conditions for the run of the Zr and Ti sublamellae (region III) were a vacuum level of ~1.5 × 10−11 torr, 35 K specimen base temperature, 5 pJ nominal laser pulse energy, a pulse repetition rate of 250 kHz, and a detection rate of 0.01 ions per pulse (1%). The conditions for the outer oxide (region IV) run were a vacuum level of ~2 × 10−11 torr, 30 K specimen base temperature, 20 pJ nominal laser pulse energy, a pulse repetition rate of 250 kHz, and a detection rate of 0.01 ions per pulse (1%). Data reconstruction and analysis were performed using the Cameca Integrated Visualization and Analysis Software version 3.8.10. Reconstruction was done using default values of 1.65 for the image compression factor and 3.30 for the field factor. A combination of methods was used including voltage evolution with evaporation field of 26 V nm−1 (region I), shank evolution with evaporation field of 37 V nm−1 (region II), and tip profile (regions III and IV). Evaporation fields were determined from the most probable charge state of majority constituents, and the tip profiles were constructed using 100kX SEM images of the pre-analysis tips.

Theory and simulation

To understand the equilibrium phase evolution of NbTiZr alloys with dissolved oxygen, MC simulations were performed with a cluster model for NbTiZr atoms on a BCC lattice with a sublattice of vacant- and oxygen-filled octahedral interstitial sites. The cluster model was fit to the total energies of a number of unique orderings represented by small supercells, where the total energies were found after relaxing the structures using DFT calculations.

DFT parameters

All DFT calculations were performed using the Perdew-Burke-Ernzerhof (PBE) functional (56) with the Vienna Ab Initio Simulation Package version 5.4.4. A plane-wave cutoff of 600 eV was used with an augmentation grid cutoff of 2400 eV for 1:2 spatial resolution. The k-point grid was automatically chosen to sample the Brillouin zone at a minimum density of 20,000 k-points per reciprocal Å3. Structural relaxations were performed with simultaneous cell relaxation, while only time-reversal symmetry was applied to the k-point grid. The self-consistency loop for the electronic solver was terminated once the change in total energy was less than 10−8 eV per atom between iterations to ensure accurate forces. Structural relaxations were terminated when none of the forces exceeded 5 meV/Å.

Cluster model fitting

Cluster model/expansions and structure enumerations were performed with icet (57) version 1.3. The training dataset was a structural enumeration (58) of symmetrically inequivalent supercells of the prototype structure, which was the BCC conventional cell with two metal sites and six octahedral interstitial sites. Because of the large number of permutations, the oxygen occupation was limited to 20% of interstitial sites (~43 at %) and only ~15% of the thousands of 2× supercells were included. Structures were created with a BCC lattice parameter of 3.397 Å, which corresponds to the average of the equilibrium BCC lattice parameters of Nb, Ti, and Zr as found in our DFT calculations. Of the 1205 structures sampled, 800 could be mapped back to the original structure after relaxation. The large distortions that prevent automatic mapping to the prototype lattice are driven by strong and symmetry breaking interactions with oxygen interstitials and by the dynamic instability of the elemental Ti and Zr BCC phases. Regardless, all structures on the convex hull were mappable to the ideal lattice, except for several oxygen-free structures where the aforementioned dynamic instability of Ti and Zr did not permit stable BCC structures. With a range of structures with widely varying levels of distortion from the prototype lattice, a set of hyperparameter-optimal cluster models was created, where each one was fit to progressively larger subsets of the structures as determined by increasing thresholds of displacement from ideal prototype lattice. By using the optimal cluster model for each distortion threshold, we remove the ambiguity introduced by the freedom of choice of hyperparameters. At increments of 0.05, 0.1, 0.2, 0.3, 0.5, and 0.8 Å of permitted average atomic displacement from ideal, the optimal cluster model was found by brute-force looping over the pair and triplet range cutoffs in increments of 0.2 Å up to a maximum of 8 Å. The triplet cutoff range was restricted to be less than the pair cutoff range. Cluster models were fit using automatically tuned LASSO regularization and 16-fold k-means cross-validation. Cluster model formation energies were referenced to the pure alloys and the spin-polarized O2 molecule’s formation energy. The optimal cluster model with a threshold 0.2 Å of average atomic displacement was chosen as a balance between an acceptable cross validation score of 13.6 meV/atom and the inclusion of 492 structures in the training set. The pair and triple cutoff ranges are 5.0 and 4.2 Å, respectively, with 161 parameters. Effective cluster interaction (ECI) data are available in table S4 and displayed in fig. S11.

MC simulations

MC simulations were performed with mchammer software included with icet. The variance-constrained, semi-grand canonical ensemble (59) was used with concentration deviation penalties, κ, of 300/atom and 16,000/atom for the metal lattice and the oxygen interstitial site sublattice, respectively. The metal lattice and the oxygen interstitial site sublattice were chosen for MC steps with probabilities of 10 and 90% to ensure good convergence with respect to the higher energy interactions of the oxygen sublattice. The simulation of the equilibrium phase separation of the lamella was simulated with a series of MC simulations of varying oxygen content using a 30 × 6 × 6 unit cell structure with periodic boundary conditions.

Paraequilibrium/equilibrium calculations using 12 × 12 × 12 unit cell structures were performed in composition intervals along the tie lines found connecting the equilibrium phases in the previous lamella simulations. For each equilibrium simulation, 8× paraequilibrium (21, 22) simulations were used to ensure ergodic sampling of the frozen lattice. Each simulation was initialized with a different, oxygen-free, equimolar structures, which were equilibrated at the temperature of interest. Thermodynamic integrations were performed along each tieline to obtain the free-energy change of segregation, i.e., the difference between single-phase paraequilibrium and two-phase equilibrium structures. Fifth-order univariate splines were fit to the free-energy change of segregation as function of composition along the tie line to and analytically differentiated to obtain smooth second derivatives of free energy as third-order splines, a procedure that compensates for the noisy nature of MC data. The spinodal points where the free energy becomes concave down, i.e., thermodynamically unstable, were found by numerically finding the roots of the second derivative splines of free energy.

WC parameters, αij, are calculated using the total number of nearest-neighbor pairs of atom types, i and j, to find the pair probability, Pij in regularly sampled snapshots from the equilibrated MC simulations. The WC parameters are calculated for each snapshot before averaging asαij=1−Pijcicj

where ci and cj are the concentrations of element types, i and j.

Acknowledgments

STEM data acquisition of region II was performed by A. Carlsson, Applications Scientist at Thermo Fisher Scientific. STEM imaging and EDS of regions III and IV were conducted in the UC Irvine Materials Research Institute. S/TEM of the base metal and spinodal regions was conducted in the Johns Hopkins University Materials Characterization and Processing Center. PDF data acquisition was performed by K. A. from NSLS-ii at the Brookhaven National Laboratory. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under contract no. DE-AC02-05CH11231 using NERSC award BES-ERCAP0023827.

Funding: This work was supported by the Multidisciplinary University Research Initiative (MURI), Office of Naval Research, Department of the Navy grant N00014-20-1-2368 (D.C.B., M.J.W., D.L.F., E.A., Y.R., C.J.B., B.W.Y.R., T.M.M., E.J.O., J.M.R., and M.L.T.); U.S. Naval Research Laboratory, Office of Naval Research, Department of the Navy (K.E.K. and P.G.C.); National Research Scientific Computing Center, Lawrence Berkeley National Laboratory, Department of Energy contract DE-AC02-05CH11231 (M.J.W. and J.M.R.); Center for Nanoscale Materials Cluster, Office of Basic Energy Sciences, Department of Energy contract DE-AC02-06CH11357 (M.J.W. and J.M.R.); Designing Materials to Revolutionize our Future Program, National Science Foundation grant DMR-1922234 (Y.R.); Graduate Research Fellowship Program, National Science Foundation grant DGE2139757 (D.C.B.); and Naval Research Enterprise Internship Program, Office of Naval Research, Department of the Navy (D.C.B.).

Author contributions: Conceptualization: D.C.B., M.J.W., E.J.O., J.M.R., and M.L.T. Methodology: D.C.B., M.J.W., J.M.R., and M.L.T. Software: M.J.W. and J.M.R. Validation: D.C.B., M.J.W., J.M.R., and M.L.T. Formal analysis: D.C.B., M.J.W., D.L.F., Y.R., K.E.K., P.G.C., G.M.V., B.W.Y.R., T.M.M., and J.M.R. Investigation: D.C.B., D.L.F., E.A., C.J.B., T.A., K.E.K., P.G.C., G.M.V., and B.W.Y.R. Resources: J.-P.C., L.P., K.E.K., T.M.M., J.M.R., and M.L.T. Writing: D.C.B., M.J.W., G.M.V., J.M.R., and M.L.T. Visualization: D.C.B., M.J.W., D.L.F., Y.R., K.E.K., G.M.V., and B.W.Y.R. Project administration: J.M.R. and M.L.T.

Competing interests: D.C.B., D.L.F., and M.L.T. disclose a provisional patent application filed with the USPTO on this work: Application Number 63/574,943 filed on 5 April 2024. The other authors declare that they have no competing interests.

Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials.

Supplementary Materials

This PDF file includes:

Figs. S1 to S11

Tables S1 to S4

Equation S1
==== Refs
REFERENCES AND NOTES

1 F. Liu, P. K. Liaw, Y. Zhang, Recent progress with BCC-structured high-entropy alloys. Metals 12 , 501 (2022).
2 O. N. Senkov, S. Rao, K. J. Chaput, C. Woodward, Compositional effect on microstructure and properties of NbTiZr-based complex concentrated alloys. Acta Mater. 151 , 201–215 (2018).
3 O. N. Senkov, C. Zhang, A. L. Pilchak, E. J. Payton, C. Woodward, F. Zhang, CALPHAD-aided development of quaternary multi-principal element refractory alloys based on NbTiZr. J. Alloys Compd. 783 , 729–742 (2019).
4 O. N. Senkov, D. B. Miracle, K. J. Chaput, J.-P. Couzinie, Development and exploration of refractory high entropy alloys—A review. J. Mater. Res. 33 , 3092–3128 (2018).
5 J. Zýka, J. Málek, J. Veselý, F. Lukáč, J. Čížek, J. Kuriplach, O. Melikhova, Microstructure and room temperature mechanical properties of different 3 and 4 element medium entropy alloys from HfNbTaTiZr system. Entropy 21 , 114 (2019).33266830
6 J. W. Yeh, S. K. Chen, S. J. Lin, J. Y. Gan, T. S. Chin, T. T. Shun, C. H. Tsau, S. Y. Chang, Nanostructured high-entropy alloys with multiple principal elements: Novel alloy design concepts and outcomes. Adv. Eng. Mater. 6 , 299–303 (2004).
7 B. Cantor, I. T. H. Chang, P. Knight, A. J. B. Vincent, Microstructural development in equiatomic multicomponent alloys. Mater. Sci. Eng. A 375-377 , 213–218 (2004).
8 Z. Lei, Y. Wu, J. He, X. Liu, H. Wang, S. Jiang, L. Gu, Q. Zhang, B. Gault, D. Raabe, Z. Lu, Snoek-type damping performance in strong and ductile high-entropy alloys. Sci. Adv. 6 , eaba7802 (2020).32596465
9 Y. Tian, W. Zhou, Q. Tan, M. Wu, S. Qiao, G. Zhu, A. Dong, D. Shu, B. Sun, A review of refractory high-entropy alloys. Trans. Nonferrous Met. Soc. Chin. 32 , 3487–3515 (2022).
10 M. Todai, T. Nagase, T. Hori, A. Matsugaki, A. Sekita, T. Nakano, Novel TiNbTaZrMo high-entropy alloys for metallic biomaterials. Scr. Mater. 129 , 65–68 (2017).
11 T. E. Whitfield, G. J. Wise, H. J. Stone, N. G. Jones, The influence of the Nb:Ta ratio on the microstructural evolution in refractory metal superalloy systems. Appl. Phys. Lett. 119 , 211901 (2021).
12 K. Naito, T. Matsui, Review on phase equilibria and defect structures in the niobium-oxygen system. Solid State Ion. 12 , 125–134 (1984).
13 J. L. Murray, H. A. Wriedt, The O−Ti (oxygen-titanium) system. J. Phase Equilibr. 8 , 148–165 (1987).
14 J. P. Abriata, J. Garcés, R. Versaci, The O−Zr (oxygen-zirconium) system. Bull. Alloy Phase Diagrams 7 , 116–124 (1986).
15 E. Gebhardt, R. Rothenbacher, Investigation of the niobium-oxygen system. II. Solution of oxygen in niobium and precipitation of oxide from supersaturated mixed crystals. Z. Metallk. 54 , 4128117 (1963).
16 B. Gorr, S. Schellert, F. Müller, H.-J. Christ, A. Kauffmann, M. Heilmaier, Current status of research on the oxidation behavior of refractory high entropy alloys. Adv. Eng. Mater. 23 , 2001047 (2021).
17 T. M. Butler, K. J. Chaput, J. R. Dietrich, O. N. Senkov, High temperature oxidation behaviors of equimolar NbTiZrV and NbTiZrCr refractory complex concentrated alloys (RCCAs). J. Alloys Compd. 729 , 1004–1019 (2017).
18 T. M. Butler, O. N. Senkov, T. I. Daboiku, M. A. Velez, H. E. Schroader, L. G. Ware, M. S. Titus, Oxidation behaviors of CrNb, CrNbTi, and CrNbTaTi concentrated refractory alloys. Intermetallics 140 , 107374 (2022).
19 K.-C. Lo, H. Murakami, U. Glatzel, J.-W. Yeh, S. Gorsse, A.-C. Yeh, Elemental effects on the oxidation of refractory compositionally complex alloys. Int. J. Refract. Met. Hard Mater. 108 , 105918 (2022).
20 C.-H. Chang, M. S. Titus, J.-W. Yeh, Oxidation behavior between 700 and 1300 °C of refractory TiZrNbHfTa high-entropy alloys containing aluminum. Adv. Eng. Mater. 20 , 1700948 (2018).
21 J. M. Rahm, J. Löfgren, P. Erhart, Quantitative predictions of thermodynamic hysteresis: Temperature-dependent character of the phase transition in Pd–H. Acta Mater. 227 , 117697 (2022).
22 J. M. Rahm, J. Löfgren, E. Fransson, P. Erhart, A tale of two phase diagrams: Interplay of ordering and hydrogen uptake in Pd–Au–H. Acta Mater. 211 , 116893 (2021).
23 J. M. Sanchez, D. De Fontaine, Anomalous diffusion in omega forming systems. Acta Metall. 26 , 1083–1095 (1978).
24 D. De Fontaine, N. E. Paton, J. C. Williams, The omega phase transformation in titanium alloys as an example of displacement controlled reactions. Acta Metall. 19 , 1153–1162 (1971).
25 A. Takeuchi, A. Inoue, Classification of bulk metallic glasses by atomic size difference, heat of mixing and period of constituent elements and its application to characterization of the main alloying element. Mater. Trans. 46 , 2817–2829 (2005).
26 R. F. Brebrick, Miscibility gap and critical point behavior for a cubic margules solution model. Calphad 2 , 17–34 (1978).
27 J. E. Morral, S.-L. Chen, High entropy alloys, miscibility gaps and the rose geometry. J. Phase Equilibr. Diffus. 38 , 319–331 (2017).
28 H. Huang, Y. Sun, P. Cao, Y. Wu, X. Liu, S. Jiang, H. Wang, Z. Lu, On cooling rates dependence of microstructure and mechanical properties of refractory high-entropy alloys HfTaTiZr and HfNbTiZr. Scr. Mater. 211 , 114506 (2022).
29 X. Yang, Y. Zhang, Prediction of high-entropy stabilized solid-solution in multi-component alloys. Mater. Chem. Phys. 132 , 233–238 (2012).
30 J. E. Morral, S. Chen, Stability of high entropy alloys to spinodal decomposition. J. Phase Equilibr. Diffus. 42 , 673–695 (2021).
31 C. H. Belcher, B. E. MacDonald, D. Apelian, E. J. Lavernia, The role of interstitial constituents in refractory complex concentrated alloys. Prog. Mater. Sci. 137 , 101140 (2023).
32 J. E. Saal, S. Kirklin, M. Aykol, B. Meredig, C. Wolverton, Materials design and discovery with high-throughput density functional theory: The open quantum materials database (OQMD). JOM 65 , 1501–1509 (2013).
33 K. C. H. Kumar, P. Wollants, L. Delaey, Thermodynamic assessment of the Ti-Zr system and calculation of the Nb-Ti-Zr phase diagram. 206, 121–127 (1994).
34 M. Sankar, R. G. Baligidad, D. V. V. Satyanarayana, A. A. Gokhale, Effect of internal oxidation on the microstructure and mechanical properties of C-103 alloy. Mater. Sci. Eng. A 574 , 104–112 (2013).
35 Z. Lei, X. Liu, Y. Wu, H. Wang, S. Jiang, S. Wang, X. Hui, Y. Wu, B. Gault, P. Kontis, D. Raabe, L. Gu, Q. Zhang, H. Chen, H. Wang, J. Liu, K. An, Q. Zeng, T.-G. Nieh, Z. Lu, Enhanced strength and ductility in a high-entropy alloy via ordered oxygen complexes. Nature 563 , 546–550 (2018).30429610
36 E. Ma, X. Wu, Tailoring heterogeneities in high-entropy alloys to promote strength–ductility synergy. Nat. Commun. 10 , 5623 (2019).31819051
37 Y. X. Ye, B. Ouyang, C. Z. Liu, G. J. Duscher, T. G. Nieh, Effect of interstitial oxygen and nitrogen on incipient plasticity of NbTiZrHf high-entropy alloys. Acta Mater. 199 , 413–424 (2020).
38 Y. Al-Khatatbeh, K. K. M. Lee, B. Kiefer, Phase relations and hardness trends of ZrO2 phases at high pressure. Phys. Rev. B Condens. Matter. Mater. Phys. 81 , 214102 (2010).
39 R. W. Armstrong, The hardness and strength properties of WC-Co composites. Materials 4 , 1287–1308 (2011).28824143
40 E. J. Haney, G. Subhash, Static and dynamic indentation response of basal and prism plane sapphire. J. Eur. Ceram. Soc. 31 , 1713–1721 (2011).
41 D. Fontaine, Simple models for the omega phase transformation. Metall. Trans. A. 19 , 169–175 (1988).
42 J. Ballor, T. Li, F. Prima, C. J. Boehlert, A. Devaraj, A review of the metastable omega phase in beta titanium alloys: The phase transformation mechanisms and its effect on mechanical properties. Int. Mater. Rev. 68 , 26–45 (2022).
43 U. Troitzsch, A. G. Christy, D. J. Ellis, The crystal structure of disordered (Zr,Ti)O2 solid solution including srilankite: Evolution towards tetragonal ZrO2 with increasing Zr. Phys. Chem. Miner. 32 , 504–514 (2005).
44 A. Kersch, M. Falkowski, New Low-Energy Crystal Structures in ZrO2 and HfO2. Phys. Status Solidi Rapid Res. Lett. 15 , 2100074 (2021).
45 Y. Al-Khatatbeh, K. K. M. Lee, From superhard to hard: A review of transition metal dioxides TiO2, ZrO2, and HfO2 hardness. J. Superhard Mater. 36 , 231–245 (2014).
46 Z. Wang, S. Huang, H. Lu, J. Liu, I. V. Alexandrov, K. Luo, J. Lu, The role of annealing heat treatment in high-temperature oxidation resistance of laser powder bed fused Ti6Al4V alloy subjected to massive laser shock peening treatment. Corros. Sci. 209 , 110732 (2022).
47 N. S. H. Gunda, B. Puchala, A. Van Der Ven, Resolving phase stability in the Ti-O binary with first-principles statistical mechanics methods. Phys. Rev. Mater. 2 , 033604 (2018).
48 H. Jiang, D. A. Stewart, Using dopants to tune oxygen vacancy formation in transition metal oxide resistive memory. ACS Appl. Mater. Interfaces 9 , 16296–16304 (2017).28436217
49 K. A. Moltved, K. P. Kepp, The chemical bond between transition metals and oxygen: Electronegativity, d-orbital effects, and oxophilicity as descriptors of metal–oxygen interactions. J. Phys. Chem. C 123 , 18432–18444 (2019).
50 W. C. Oliver, G. M. Pharr, An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments. J. Mater. Res. 7 , 1564–1583 (1992).
51 G. Ashiotis, A. Deschildre, Z. Nawaz, J. P. Wright, D. Karkoulis, F. E. Picca, J. Kieffer, The fast azimuthal integration Python library: pyFAI. J. Appl. Cryst. 48 , 510–519 (2015).25844080
52 P. Juhás, T. Davis, C. L. Farrow, S. J. L. Billinge, PDFgetX3: A rapid and highly automatable program for processing powder diffraction data into total scattering pair distribution functions. J. Appl. Cryst. 46 , 560–566 (2013).
53 L. B. McCusker, R. B. Von Dreele, D. E. Cox, D. Louër, P. Scardi, Rietveld refinement guidelines. J. Appl. Cryst. 32 , 36–50 (1999).
54 K. Thompson, D. Lawrence, D. J. Larson, J. D. Olson, T. F. Kelly, B. Gorman, In situ site-specific specimen preparation for atom probe tomography. Ultramicroscopy 107 , 131–139 (2007).16938398
55 M. Schaffer, B. Schaffer, Q. Ramasse, Sample preparation for atomic-resolution STEM at low voltages by FIB. Ultramicroscopy 114 , 62–71 (2012).22356790
56 J. P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple. Phys. Rev. Lett. 77 , 3865–3868 (1996).10062328
57 M. Ångqvist, W. A. Muñoz, J. M. Rahm, E. Fransson, C. Durniak, P. Rozyczko, T. H. Rod, P. Erhart, ICET—A Python library for constructing and sampling alloy cluster expansions. Adv. Theory. Simul. 2 , 1900015 (2019).
58 G. L. W. Hart, R. W. Forcade, Generating derivative structures from multilattices: Algorithm and application to hcp alloys. Phys. Rev. B 80 , 014120 (2009).
59 B. Sadigh, P. Erhart, Calculation of excess free energies of precipitates via direct thermodynamic integration across phase boundaries. Phys. Rev. B 86 , 134204 (2012).
