1. Introduction
Ultra high temperature ceramics, including the borides, carbides, and nitrides of the early transition metals, exhibit the highest melting temperatures of any known compounds, making them the primary candidates for applications exposed to extreme thermal conditions [
1], [
2]. The most severe of those conditions, however, exceeds the melting temperature of every oxide these ceramics can form, so refractoriness stops discriminating between candidates and the volatility of a molten scale takes its place. Two representative applications define the present study. A hypersonic leading edge is exposed to dissociated air at a small radius of curvature, where the stagnation point heat flux scales inversely with the leading edge radius and drives the material toward its thermal limits. A rocket nozzle throat is exposed to high temperature, high pressure combustion gases at the sonic plane, where the convective heat flux reaches its maximum throughout the engine. Both applications require a material that retains its structural integrity and dimensional stability under extreme thermomechanical loading.
High entropy carbides have substantially expanded the compositional design space for ultra high temperature ceramics. Sarker et al. introduced an entropy forming ability descriptor and demonstrated that nine of twelve five metal carbides crystallize as single phase rock salt structure when sintered near 2473 K, with carbon vacancies on the anion sublattice providing an additional contribution to the configurational entropy [
3]. Subsequent studies examined phase stability and mechanical behavior across transition metal carbides, reporting hardness values that significantly exceed those predicted by the rule of mixtures [
4], [
5]. The principal advantage of these materials is their compositional flexibility. Incorporating five or more metallic species on a single sublattice allows fracture toughness, sinterability, and thermal transport to be tailored simultaneously.
Most of these design strategies have focused on maximizing bulk single phase stability, although this is not the property that determines service lifetime. A ceramic that remains thermodynamically stable as a single phase at 3000 K can still fail through surface degradation, where the formation and evolution of the oxide scale determines the oxidation resistance. Zirconium carbide exhibits measurable oxidation at 1570 to 1770 K within minutes [
2] while the transition metal diborides form a protective borosilicate glass layer. This protection layer starts deteriorating above approximately 1370 K, causing the oxidation kinetics to transition from parabolic to paralinear and ultimately to linear behavior as the protective scale is progressively lost [
6], [
7], [
8]. Carbides lack silicon and cannot form a protective borosilicate glass. The oxidation behavior of the quinary carbide (HfZrTiTaNb)C and its diboride has been studied directly, confirming that oxide scale chemistry, rather than bulk phase stability, governs oxidation performance and service lifetime [
9]. In this class of materials, bulk phase stability is therefore well established, and the outstanding challenge lies in understanding and controlling the surface.
A second limitation is in the way these materials are typically designed. Many computational studies begin with density functional theory, investing significant computational effort in refining bulk stability margins before inexpensive screening methods have reduced the number of candidate compositions. Here, the order is reversed. The approach follows the process, structure, property, performance framework proposed by Olson [
10] and treats kinetics as the practical limit on thermodynamic stability. A phase may be thermodynamically stable yet remain unsuitable because its oxide scale volatilizes too rapidly. The present work develops a screening methodology that minimizes computationally intensive simulation: a closed form stability criterion reduces the search space, equilibrium thermochemistry identifies the oxide scale and its volatile species, and a dual regime recession model converts volatility into recession rates. The methodology is demonstrated for the (Hf,Zr,Ti,Ta,Nb)C system and selected lower component subsets. The superior oxidation resistance of hafnia and zirconia relative to tantala and niobia is already well established. The present work translates that qualitative ordering into quantitative recession rates under representative rocket nozzle throat conditions. It also quantifies the effect of hafnium content and derives the niobium oxide partial pressure from the JANAF thermochemical functions. The methodology is readily applicable to other oxidation limited ceramic systems for which vaporization data are available for the corresponding oxide scales.
2. Computational methodology
2.1. Configurational entropy stability screen
Each carbide is treated as a single rocksalt phase with mixing confined to the metal sublattice. The molar Gibbs energy of mixing per formula unit is written as a regular solution,
$G = 0.5xLx + RTΣ_{i}x_{i}lnx_{i}$
where x is the vector of metal fractions and L - the binary interaction parameters. Each interaction parameter is obtained from the mixing enthalpy of the corresponding equiatomic binary:
$L_{ij} = 4dH_{mix}(0.5)$
The ten metal pair enthalpies are calculated in the present work by density functional theory. All calculations employed Quantum ESPRESSO 7.3.1 with the Perdew Burke Ernzerhof exchange correlation functional and ONCV SG15 pseudopotentials. The computational procedure consisted of a variable cell relaxation followed by a single point energy calculation. The reference monocarbides were modeled in the two atom primitive rock salt unit cell, while each equiatomic binary carbide was modeled in an eight atom supercell with the two metal species ordered on the metal sublattice. The sensitivity of these enthalpies to the electronic structure method and to the ordered arrangement is examined in Section S7; neither changes the stability verdicts. The density functional calculations were performed on the Sol supercomputer at Arizona State University. No interaction parameter is taken from a prior thermodynamic assessment. At 0 K, the Gibbs free energy of mixing is equal to the enthalpy of mixing, allowing all interaction parameters to be determined directly without additional simulations. Single phase stability is evaluated from the eigenvalues of the Gibbs free energy Hessian over the composition simplex. A positive definite Hessian indicates a stable solution, and the temperature at which the smallest eigenvalue becomes zero defines the equiatomic spinodal temperature. The diboride is represented using a single stabilizing interaction parameter fitted to the reported equiatomic mixing enthalpy [
11]. The compositions considered include the equiatomic quinary carbide synthesized by Pak et al. [
12], three lower component carbide subsets, the binary HfC ZrC system, the equiatomic quinary diboride, and ZrB
2 with SiC as a benchmark material.
2.2. Oxide scale and volatile pressure
Each stable ceramic is paired with its service environment, and the oxide scale is identified from the thermodynamic stability and refractoriness of the candidate oxides. The recession rate is governed by the equilibrium partial pressure of the volatile metal bearing species above the oxide scale. General Gibbs energy minimization with Reaktoro [
13] and the NASA CEA database [
14] was attempted first and abandoned: that database carries neither the condensed oxide phases nor the volatile oxide species required for hafnium, zirconium, tantalum or niobium, so the refractory scales cannot be treated that way. The volatile species partial pressures reported in this work are instead obtained from third law thermodynamic calculations using tabulated Gibbs free energy functions together with the standard enthalpy of reaction at 298 K ( Table S1). Each calculation is performed for the vaporization reaction specific to the oxide of interest, rather than assuming congruent sublimation. This distinction is important because none of these refractory oxides sublimates congruently according to its stoichiometric formula. Hafnia vaporizes through the formation of HfO(g) and O(g). The standard enthalpy of formation for condensed HfO
2 is taken from combustion calorimetry measurements [
15], while the thermodynamic functions for gaseous HfO are taken from a combined quantum chemical and Knudsen effusion study [
16].
Zirconia vaporizes to both ZrO
2(g) and ZrO(g), while tantala undergoes incongruent vaporization from a melt to TaO
2(g) and oxygen [
17], [
18]. Niobium is treated using the oxide phase stable under the relevant service conditions. NbO
2 is used for the reducing nozzle throat environment, where the volatile species partial pressure is calculated from the JANAF thermodynamic functions using the third law method [
19], [
20] ( Tables S2a and S2b). Nb
2O
5 is used for the oxidizing hypersonic environment [
21]. This assignment is set by the oxygen potential rather than temperature, and the throat conditions fall near the NbO
2-Nb
2O
5 boundary, where the carrier pressure is insensitive to the oxide phase (Section S16). For oxides that evolve oxygen during vaporization the gas phase composition is solved subject to mass balance. The vaporization channels are therefore coupled and cannot be treated independently. All constituent oxides share a common gas phase and oxygen potential, and the metal bearing pressure is determined by the coupled equilibrium rather than by a linear combination of the pure oxide pressures. This distinction is important because zirconia is more volatile than hafnia, and the oxygen released by zirconia suppresses the HfO(g) vaporization channel, an effect that is not captured by a linear summation. The scale composition is treated using two limiting cases. The first assumes that the scale composition reflects the bulk metal ratio and produces the higher predicted recession rate. The second represents a steady state in which the vaporization rate of each metal equals its supply rate through recession, resulting in preferential enrichment of the least volatile component Unless otherwise stated, results are reported using the bulk composition limit, while both limiting cases are presented in Section S13. Hafnia and zirconia are treated as a continuous near ideal solid solution. The activities of the pentoxides are varied over two orders of magnitude to assess the sensitivity of the predicted recession rates. Water vapor introduces a second volatilization channel above the condensed oxide, MO
2(c) + n H
2O(g) = MO
x(OH)
y(g), which is second order in water for the tetrahydroxides. This channel is therefore a minor contribution of the water vapor pressures at one bar, as commonly considered in the environmental barrier literature, but becomes dominant at the substantially higher water vapor pressures characteristic of a hydrogen engine throat. Thermochemical functions for hafnium and zirconium hydroxides are taken from coupled cluster calculations [
22] (Section S15). Two service environments are considered: dry air at 2200 K for the hypersonic leading edge and the combustion gas of the selected engine at the isentropic nozzle throat conditions defined in Section 2.5. For the latter, the equilibrium gas composition and oxygen potential are calculated using Cantera [
23] (Section S16).
2.3. Evaporation coefficients
Free surface, or Langmuir, vaporization proceeds at a rate lower than that predicted from the equilibrium vapor pressure by a factor known as the evaporation coefficient. For refractory oxides, this coefficient is generally small and exhibits considerable scatter, with measured free surface vaporization rates commonly one to four orders of magnitude lower than the equilibrium limit. This observation provides the basis for the bounding approach adopted in the present work. An evaporation coefficient of unity defines a rigorous upper bound on the recession rate, while the actual recession rate must lie below this limit. Experimentally measured values include 0.002 to 0.006 for silica, approximately 0.11 for the MgO(100) surface, and about 0.5 for silicon monoxide [
24] ( Table S4). No accepted evaporation coefficient is available for hafnia or zirconia and a range of 0.01 to 0.1 is adopted as a working estimate and evaluated alongside the limiting case of unity.
2.4. Dual regime recession
The equilibrium partial pressure of a volatile species does not directly determine the recession rate. The relation between the two depends on the rate limiting transport mechanism, which differs between the two service environments. At the hypersonic leading edge, material is removed by free molecular vaporization, and the molar flux of species i is described by the Hertz Knudsen Langmuir equation: $J_{i} = α_{i} p_{i} \sqrt{2πM_{i}RT}$, where αi is the evaporation coefficient, pi the equilibrium partial pressure of species i, Mi its molar mass, R the gas constant and T the surface temperature. The Knudsen number over representative hypersonic flight conditions does not reach the free molecular threshold (Table S9), so this treatment is reported as a limiting case rather than as a general leading edge prediction. At the nozzle throat, the surrounding gas is sufficiently dense that transport through the boundary layer limits mass loss. The molar flux is $J_{i} = h_{m,i} p_{i}/(RT)$ where hm,i is the mass transfer coefficient and the concentration of the volatile species in the free stream is assumed to be negligible. Under these conditions, the evaporation coefficient does not influence the recession rate. The recession velocity is obtained from the total metal bearing flux according to $v=(Σ_{i}v_{i}J_{i})V_{m}$, where vi is the number of metal atoms per volatile species and Vm is the molar volume of the carbide per mole of metal atoms. Carbon is removed primarily as carbon monoxide and does not contribute to geometric recession.
2.5. Reference engines and throat transport
The nozzle throat mass transfer coefficient is calculated from the engine operating conditions rather than treated as an adjustable parameter. The flow is assumed to be sonic at the throat, allowing the local temperature and pressure to be obtained from the chamber conditions using the isentropic flow relations. The binary diffusivity of each volatile oxide species in the combustion gas and the gas viscosity are calculated using Chapman Enskog kinetic theory with the Neufeld collision integrals [
25], [
26]. The Sherwood number is determined from the Bartz turbulent convection correlation together with the Chilton Colburn heat and mass transfer analogy as S
h = 0.026 Re0.8 Sc1/3, evaluated at the nozzle throat using the throat diameter as the characteristic length [
27]. The mass transfer coefficient is then calculated from h
m = Sh D /L where (D) is the binary diffusivity and (L) is the throat diameter. Two rocket engines are considered to span a range of combustion environments: the RS-25, an oxygen hydrogen staged combustion engine with a reducing, steam rich throat, and the F-1, an oxygen kerosene engine with a carbon rich throat. The chamber pressure, chamber temperature, and throat diameter are taken from published engine specifications [
28], [
29]. The throat transport parameters derived for each engine are collected in Table S3.
2.6. Surface temperature
The throat gas temperature is different than the surface temperature, and the predicted recession rate depends exponentially on the latter. The gas side heat flux is $h_{g}(T_{aw} -T_{w})$, where hg is obtained from the Bartz form of the turbulent correlation and Taw is the recovery temperature. Because the gas temperature exceeds the surface temperature at the throat, heat is removed from the surface only by conduction through scale and liner and by the ablation products (Section S11). At the RS-25 wall temperature the resulting heat flux is 96 MW m-2 consistent with the published range.
3. Results
3.1. Stability screen
Thermodynamic stability does not eliminate any composition in the candidate set. The equiatomic quinary carbide has an equiatomic spinodal temperature of 2240 K, below the reported single phase synthesis temperature of 2473 K [
3], [
4], consistent with the experimental observation of a single phase material. A global two phase common tangent construction applied to the same Gibbs free energy predicts a miscibility gap boundary at 2253 K, only 13 K above the mean field spinodal. The local stability criterion provides an accurate estimate of the miscibility gap, with the equiatomic quinary composition lying close to the critical point. The lowest energy instability corresponds to partitioning between titanium and zirconium, which have the largest positive binary mixing enthalpy in the system, 17.3 kJ/mol. The titanium hafnium interaction is the second largest with 12.0 kJ/mol. This proximity to the critical point leaves only a limited stability margin. A compositional deviation of approximately 2.8 at.% from the equiatomic composition, toward simultaneous enrichment in titanium and zirconium, raises the spinodal above the synthesis temperature. Removing titanium substantially increases the stability margin. The equiatomic spinodal decreases by more than 1800 K to 381 K for the Hf Zr Ta Nb subset and falls further for the lower component systems. Several of these compositions exhibit negative binary mixing enthalpies for every interaction pair over the entire temperature range. Titanium removal therefore not only lowers the spinodal temperature but also eliminates the sensitivity to modest compositional deviations. The quinary diboride is stabilized in all binary pairs and exhibits no miscibility gap. No compositional spinodal instability is predicted for any candidate within the adopted rocksalt regular solution model (
Table 1). This test does not exclude long range ordering, ordered carbon vacancy phases, or decomposition into structures outside the assumed configurational space; it therefore constrains only the tendency for compositional unmixing and does not establish complete phase stability. All candidate materials satisfy the stability criterion.
3.2. Oxide scale and volatility in two environments
The key result is that one composition performs best in both environments, although the protection mechanism and overall viability is different. At the nozzle throat, oxidation resistance is controlled by oxide volatility and the presence of water vapor. At the throat gas temperature of 3248 K, the melting points of all candidate oxides, including hafnia and zirconia, are exceeded therefore the recession is determined by oxide vaporization. This is not, by itself, detrimental. Transition metal diborides possess an oxidation resistance from a liquid borosilicate scale [
6], [
7], [
8], and the recession model treats the oxide scale as a condensed phase regardless of whether it is solid or liquid. The melting point instead determines the extent to which a molten scale can resist removal by gas shear, a mechanical process that lies outside the thermochemical model developed in Section 2.4. The volatility ranking is unambiguous. Oxides richer in tantalum and niobium exhibit the highest vapor pressures and the greatest recession rates. Water vapor further accelerates recession by converting boria and silica into volatile hydroxide species, while the low oxygen partial pressure promotes active oxidation of silicon carbide through the formation of gaseous SiO. As a result, borosilicate forming ceramics perform poorly under nozzle throat conditions. Increasing the hafnium and zirconium content provides the greatest reduction in recession rate, and the candidate materials rank accordingly. Candidates are ranked by a surface protection index, SPI =-log
10(v/v
crit), where v is the predicted recession in the environment concerned and v
crit is the criterion defined in Section 3.5. Positive values indicate that the criterion is satisfied and each unit corresponds to one order of magnitude in recession rate, making the ranking directly reproducible from the stated equations and inputs. The index is calculated using the dissociation channel, which is available consistently for all metals. Because hydroxide thermochemical functions are unavailable for niobium, the hydroxide channel is reported separately and is not included in the ranking.
Fig.1 presents the resulting recession ranking, while
Fig.2 compares the oxide melting temperatures with the two service temperature ranges.
At the hypersonic leading edge, the ranking changes. In dry air at 2200 K, silica remains protective, allowing the ZrB2 + SiC system to outperform the diborides in the nozzle throat and rank near the middle of the candidate set. This improvement reflects the stability of the borosilicate scale under dry, oxygen rich conditions at lower temperatures. The boron only diboride, which cannot form a borosilicate scale, remains the least resistant because of continued boria vaporization. The carbides retain the same ordering, with hafnium and zirconium rich compositions providing the greatest oxidation resistance. Five of the seven candidates clear the adopted criterion here and none does at the throat, so the environment sets viability rather than the choice of composition.
3.3. Absolute recession from measured oxide volatility
With the equilibrium partial pressures established, the recession rate can be calculated. Each oxide follows the vaporization reaction of its equilibrium vapor species ( Section 2.2). Hafnia vaporizes by dissociation to HfO(g) and O(g). No gaseous HfO
2 has been detected, either by mass spectrometry at a sensitivity limit of approximately 2% of the HfO signal [
30] or by Knudsen effusion mass spectrometry over a hafnia containing ternary oxide [
31]. Zirconia produces both ZrO
2(g) and ZrO(g), while tantala undergoes incongruent vaporization to TaO
2(g) and oxygen. No gaseous tantalum pentoxide species is reported in the thermodynamic literature. Niobium is represented by NbO
2 under nozzle throat conditions and by Nb
2O
5 in air, consistent with the oxygen potential of each environment discussed in Section 2.2. For niobium the third law pressures agree with an independent extrapolation of experiment to within 6%. For the (Hf,Zr)C oxide scale, the equilibrium pressure at the hypersonic leading edge is 5.2 × 10
-4 Pa, corresponding to a free molecular recession rate of 0.22 μm/h at the bounding coefficient. Under nozzle throat conditions, the same oxide scale reaches 132 Pa and a recession rate of 265 μm/h in the RS-25 engine. Replacing zirconium with tantalum increases the pressure to 933 Pa and the recession rate to 1880 μm/h. Replacing it with niobium gives 1580 Pa and 3190 μm/h, while the four-component oxide scale containing both tantalum and niobium reaches 1820 Pa and 3670 μm/h. The F-1 engine produces the same ordering at approximately half the recession rate. The resulting ranking is C2_HfZr > C3_HfZrTa > C3_HfZrNb > C4_noTi reflecting the increasing volatility of the added oxides. Hafnium provides the largest reduction in recession rate. At the RS-25 nozzle throat, a pure zirconia scale recedes at 438 μm/h, the coupled equiatomic HfO
2-ZrO
2 scale at 265 μm/h, and a pure hafnia scale at 50 μm/h ( Table S6).
Fig. 3 summarizes the recession rates for the four highest ranked compositions, and
Table 2 and Table S5 list the corresponding values.
3.4. Consistency check against measured ablation
No directly comparable recession measurements exist for the hafnium rich (Hf,Zr)C composition recommended here under clean liquid rocket engine throat conditions. The comparison presented below is a consistency check rather than a validation. The objective is to determine whether the predicted recession falls where the governing mechanisms would place it relative to experiments conducted under different conditions. The model calculates only the thermochemical contribution to recession from the vaporization of volatile metal oxides in a clean combustion gas.
Table 3 compares these predictions with measured recession rates in progressively more complex environments. The atmospheric oxyacetylene tests probe the same thermochemical mechanism as the present model. In these experiments, recession increases with the content of volatile oxide formers, with the boride composite containing both boria and silica exhibiting the highest recession [
32], followed by the carbide containing silica [
33]. These experiments do not include the recommended composition and were conducted several hundred kelvin below nozzle throat temperatures and as a result provide a comparison of the proposed mechanism rather than a direct comparison of recession rates. The solid rocket motor measurements are approximately three orders of magnitude higher. This difference is expected because aluminized solid propellants introduce two additional recession mechanisms that are not included in the present model: enhanced convective mass transfer at chamber pressures of 50 to 100 bar and mechanical erosion by alumina particles in the exhaust [
34], [
35] (Section S17). These measurements therefore bound the model rather than validate it, and the predicted rates should be read as the thermochemical contribution to recession under the stated assumptions rather than as a complete service life prediction. High entropy carbide and diboride nozzle throat materials are now being evaluated in solid rocket motor tests, providing a foundation for future comparisons [
36].
3.5. Design map and sensitivity
Fig. 4 maps the required hafnium fraction against throat temperature and allowable recession. The required hafnium fraction depends on the allowable recession rate, which is not specified for liquid rocket engine throats in the available literature. Two independent estimates provide a reasonable range. A 1% increase in throat area corresponds to approximately a 1% reduction in chamber pressure at constant mass flow, giving an allowable recession rate of 87 μm/h over the certified 7.5 h service life of the RS-25. A NASA recession model of a silicon carbide thrust chamber liner in oxygen hydrogen combustion identified a recession rate of approximately 200 μm/h at 1973 K as negligible relative to the expected component lifetime [
37]. We adopt 100 μm/h, between the two. At that limit the equiatomic (Hf,Zr)C composition does not satisfy the RS-25 nozzle throat requirement and a hafnium fraction of 0.91 is required to reduce the predicted recession rate. At the cooler F-1 throat temperature, a hafnium fraction of 0.50 is sufficient, placing the equiatomic composition at the threshold. At the hypersonic leading edge, the predicted recession rates remain orders of magnitude below the adopted limit for all compositions considered. The hafnium requirement is therefore imposed by the nozzle throat rather than the leading edge.
The required hafnium fraction depends on the adopted limit, a design criterion rather than a material property, and on every other input the calculation carries. Propagating surface temperature, reaction enthalpies, transport and the criterion together gives medians of 0.79 for the RS-25 and 0.17 for the F-1, with fifth to ninety fifth percentile ranges spanning most of the join (Section S12). The ordering is robust where the values are not. Limits below the pure hafnia rate cannot be met by composition at all, which is the reusable engine case and motivates the graded architecture of Section 4.
The sensitivity to temperature has important design implications. The RS-25 and F-1 throat temperatures differ by only 105 K, approximately 3% on an absolute scale. Because the equilibrium vapor pressure varies exponentially with temperature, raising only the throat temperature from the F-1 to the RS-25 value at unchanged transport increases the required hafnium fraction from 0.50 to 0.87. The remaining rise to 0.91 follows from the higher mass transfer coefficient of the RS-25, which its chamber pressure, throat diameter and gas composition place 18% above that of the F-1. Temperature therefore supplies 0.37 of the total shift of 0.41, corresponding to more than 40% of the metal sublattice. A composition qualified on a kerosene engine is therefore not qualified on a hydrogen engine of comparable thrust.
Fig. 5 propagates the uncertainty reported for each reaction enthalpy. The uncertainty is applied as an additive shift in enthalpy (kJ/mol), consistent with the reported thermochemical data. A fractional perturbation is not used because the reaction enthalpies exceed 1000 kJ/mol, making such scaling physically unrealistic. Because the reaction enthalpy appears in the exponential term of the equilibrium constant, the resulting uncertainty is asymmetric and differs among the oxides. For the hafnia zirconia recession, the uncertainty corresponds to a factor of 1.3, giving an RS-25 recession rate between 238 and 297 μm/h about the nominal value of 265 μm/h. The niobium bearing recessions exhibit the largest uncertainty, a factor of 2.4, due to the larger uncertainty in the thermochemistry of gaseous NbO
2. Replacing vendor database values with evaluated thermochemical data reduced this uncertainty by approximately a factor of three. Despite these uncertainties, the compositional ranking is unchanged. The hafnium zirconium carbide remains the least volatile composition, the niobium bearing scales remain the most volatile, and the separation between them exceeds an order of magnitude across the entire uncertainty range. The design conclusions on hafnium enrichment and niobium additions are therefore insensitive to the thermodynamic uncertainty.
3.6. Effect of non-ideal scale mixing
The assumption that each oxide activity equals its cation sublattice fraction is the principal simplification remaining in the volatility model. Its influence is therefore evaluated directly. Hafnia and zirconia form a continuous near ideal solid solution, so the predicted recession of the (Hf,Zr)C scale is largely insensitive to the activity model. The coupled oxygen potential described in Section 2.2 is retained and continues to distinguish the result from a linear combination of the pure oxide vapor pressures. The uncertainty is due to the pentoxides, which may dissolve into the refractory oxide matrix with activities below unity thus lowering their volatility, or exsolve as a separate phase with activities approaching unity.
Fig. 6 varies the activity coefficients of the tantalum and niobium oxides over two orders of magnitude, taking each in the oxidation state stable in the environment concerned, Ta
2O
5 and NbO
2 at the nozzle throat. No candidate reaches the adopted service criterion at any activity in the range, although the ranking itself does not survive the whole sweep. Below an activity coefficient of approximately 0.1 the tantalum and niobium bearing scales fall below the pentoxide free binary: at 0.01 they recede at 194 and 207μm/h against 265 μm/h for the binary. The crossing is a dilution effect rather than a protection mechanism, since a pentoxide dissolved at low activity also dilutes the hafnia and zirconia that carry the volatility, and it is strongest for the composition holding the most pentoxide. As the pentoxide contribution diminishes, each composition approaches the recession rate set by its hafnia zirconia matrix, between approximately 130 and 180 μm/h depending on the HfO
2/ZrO
2 ratio. These results show that nonideal oxide mixing cannot qualify tantalum or niobium containing compositions for the hotter nozzle throat. The recession floor is determined by the hafnia zirconia matrix, and reducing that floor requires increasing the hafnium content rather than changing the oxide activity model.
3.7. Reproducing the experimental design rule at extreme temperature
The proposed screening pathway also reproduces a design rule established experimentally. Oxidation resistance in high entropy carbides at extreme temperatures has been extensively investigated, with the current state of the art reaching a recession rate of 2.7 μm/s at 3873 K for a (Hf,Ta,Zr,W)C carbide developed by Wen et al. through extensive compositional and microstructural optimization [
38]. The oxidation resistance of that material is provided by a hafnia and zirconia rich oxide backbone, the same oxide selection identified by the present thermodynamic descriptor.
Fig. 7 shows the equilibrium partial pressures of the volatile metal species at 3873 K for the four constituent metals. Tungsten forms no condensed oxide above 3000 K and exists as either liquid metal or as gaseous oxide with high volatility. The stable phase is determined primarily by the oxygen potential and not the temperature. Below approximately 5 × 10?? bar as imposed by a reducing torch environment tungsten is the least volatile metal present; a rocket throat operates at oxygen potentials several orders of magnitude higher. The predicted ranking matches the experimental design rule: hafnium and zirconium form the refractory, low volatility backbone; tantalum and, in particular, niobium increases oxide volatility. The descriptor therefore identifies the same compositional trend before synthesis and points toward hafnium and zirconium rich compositions.
The scope of this result should be stated clearly. The model predicts the thermochemical contribution to oxidation resistance and identifies the preferred metal selection. The exceptional performance reported by Wen et al. [
38] also depends on microstructural mechanisms, including a self healing molten oxide containing high melting tungsten particles, which are outside the scope of the present composition based model. Tungsten is a notable example because it is retained as metallic tungsten under reducing conditions rather than being removed as a volatile oxide. The thermodynamic screening approach and microstructural optimization are therefore complementary.
4. Discussion
The principal design outcome is that the service environment determines whether a composition is viable, rather than simply which composition is preferred. A material optimized for a hypersonic leading edge has oxidation resistance from a protective borosilicate glass that forms in dry, oxygen rich air. Under nozzle throat conditions, water vapor and low oxygen partial pressure destabilize this protection due to boria volatilization and active oxidation of silicon containing phases. A material optimized for the nozzle throat instead relies on the formation of refractory mixed hafnia-zirconia oxide layer. Such compositions also perform well at the leading edge, although they do not benefit from the self-healing borosilicate layer.
The predicted recession rates identify the role of each alloying element. Tantalum and niobium, which provide the configurational entropy that stabilizes high entropy carbides, also form the most volatile oxides and substantially increase the predicted recession rate. Hafnium has the opposite effect. The required hafnium content is best expressed as a design threshold rather than as a compositional preference, and is determined within the single-phase stability field established in Section 3.1.
Oxidation resistance increases monotonically as tantalum, niobium, titanium, boron, and silicon are removed from the surface oxide and the refractory hafnia and zirconia fraction increases. The same elements, however, contribute to fracture toughness, sinterability, and the configurational entropy that stabilizes high entropy carbides. Improving oxidation resistance and maximizing configurational entropy therefore require different compositions. The solution is structural rather than compositional: a high entropy carbide substrate providing mechanical performance and processability, combined with a hafnium rich surface or coating forming a refractory HfO2-ZrO2 oxide layer. A second route is already implemented in flight systems. Fuel rich injection along the wall is commonly used to reduce heat flux and also lowers the local water vapor pressure. Because recession scales approximately with the square of water vapor pressure, the same approach provides both thermal and chemical protection.
The proposed methodology reduces the computational effort required for materials screening. Density functional theory is used only to calculate the ten binary mixing enthalpies needed for the regular solution model. Oxide vapor pressures are obtained from measured thermochemical data, and recession rates are calculated using kinetic theory together with established mass transfer correlations. This sequence reserves the expensive calculations for what cannot be measured.
Several limitations define the scope of the present study. The vapor pressures of the hafnium, zirconium, and tantalum oxides are based on primary calorimetric measurements and evaluated thermodynamic compilations. The niobium pressure under nozzle throat conditions is calculated from the JANAF thermodynamic functions using the third law method [
19], [
20], consistent with the approach of Kamegashira et al., and agrees with an independent extrapolation of their measured vapor pressures to within 6%. Under oxidizing conditions, niobium is treated as Nb
2O
5, with vaporization data taken from Matsui and Naito [
21] and extrapolated only from 1988 to 2200 K. The scale composition is now solved from the steady-state flux balance rather than assumed equal to the bulk, which removes the principal free parameter of the earlier treatment. The enthalpy uncertainties propagated in Section 3.5 preserve the compositional ranking. The standard entropy of monoclinic hafnia illustrates the remaining uncertainty in the thermodynamic data. Reported values range from 59.33 J mol?1 K?1 from adiabatic calorimetry [
39] to 56.15 J mol
-1 K
-1 from low temperature calorimetry [
40], with an independent first principles assessment giving an intermediate value of 57.7 J mol
-1 K
-1 [
41]. Using the lowest reported entropy increases the predicted hafnia vapor pressure by approximately 20% and reduces the hafnia zirconia volatility contrast from about nine to seven, but changes the required hafnium fraction by no more than 0.01. The design thresholds are therefore insensitive to the available thermodynamic data. The nozzle throat calculations assume vaporization consistent with the stoichiometry of the condensed oxide. Lower oxide channels that are not represented would raise the equilibrium vapor pressures, so their omission makes the predicted recession rates an underestimate rather than a conservative bound. The consistency check in Section 3.4 compares the model with published measurements on related ceramics rather than direct measurements on the recommended composition. Shear driven drainage of a molten scale scales as the square of its thickness and becomes comparable to the thermochemical rate at micrometre thicknesses (Section S10), so the predicted rates are a lower bound on total recession. A recession measurement for a hafnium rich carbide in a clean liquid rocket engine environment would provide a direct test of the model. The transport calculations use estimated Lennard Jones parameters, and the sensitivity of the predicted rates to the mass transfer coefficient is reported in Section S14. None of these limitations requires additional density functional theory calculations to address.
The configurational entropy is treated as ideal although the same pair interactions that determine the mixing enthalpy also produce short range order and reduce it. The recursive entropy method of Hew et al. [
42] was applied to quantify this reduction resulting in less than % for the quinary carbide and negligible value above 2000 K. More details on the analysis are listed in Section S8. The stability conclusions in Section 3.1 are unaffected.
The two reference engines represent two points on a continuous design requirement. The horizontal axis of
Fig. 4 shows the throat gas temperature across the range relevant to liquid propellant engines. The contours shift only slightly between the mass transfer coefficients of the two engines, indicating that the chamber pressure, throat diameter, and gas composition together have smaller effect than the temperature.
Three distinct regimes follow. Below approximately 3060 K, every composition along the HfC-ZrC join satisfies the 100 μm/h recession criterion and increasing the hafnium content provides no additional benefit Between approximately 3060 and 3340 K, the required hafnium fraction increases rapidly with temperature, making the optimum composition strongly engine specific. Above approximately 3340 K even a pure hafnia scale exceeds the recession limit, so no composition on the join can meet the design criterion. Thus, improved performance must come from active cooling, protective coatings, or the graded architecture proposed here. The RS-25 throat lies only 92 K below this upper limit. Consequently, carbide surfaces operating in hydrogen fueled engines of this class are close to the maximum performance attainable through composition alone, indicating that oxide volatility imposes a more stringent design constraint on this material family than bulk phase stability.
5. Conclusions
Single phase stability does not determine the service life of high entropy carbides in rocket applications. Service life is controlled by the chemistry and volatility of the oxide layer. The thermodynamic and kinetic screening methodology developed here leads to three principal conclusions. First, all candidate carbides remain single phase during synthesis, and the hafnium and zirconium rich compositions selected by the oxidation analysis remain stable throughout the service temperature range. Material selection is therefore controlled by oxidation rather than bulk phase stability. Second, the preferred composition depends on the service environment. Second, the hafnium rich and zirconium rich carbide performs best in both service environments, while the environment determines whether any composition is viable. Five of seven candidates satisfy the adopted criterion at the leading edge and none at the throat. Under throat conditions, recession is determined by wall temperature and water vapor pressure rather than by the underlying carbide composition, which motivates the combined use of a graded architecture and a fuel rich wall film.
Third, the analysis provides quantitative thermochemical recession rates and the corresponding compositional requirements from an adopted service criterion rather than from an intrinsic material limit. The methodology limits density functional theory to the binary interaction enthalpies and otherwise relies on measured and evaluated thermochemical data. Direct recession measurements on a hafnium rich carbide in a well characterized combustion gas remain necessary for experimental validation. Since the dominant volatilization channel is determined by the gas at the wall, propellant chemistry becomes a materials design variable, and the present framework provides a basis for exploring that dependence.
Declaration of Generative AI and AI-assisted technologies in the writing process
During the preparation of this work the authors used Anthropic's Claude, through the Claude Code assistant, in an agentic but not autonomous capacity, as a limited facilitator that helps run the workflow and assists with coding, language editing, and grammar checking rather than as a research tool. The thermodynamic analysis and all verification and editorial checks are performed by proprietary in-house Odinzen software, a multi-program architecture independent of and not native to Claude that the authors designed and directed; the model runs the workflow while this author-built system, and not Claude, performs the research. The authors reviewed and approved all output and take full responsibility for this article. All figures, including the graphical abstract, are reproducible visualizations generated by this software directly from the underlying data, and no generative-AI image tool was used.
Declaration of Competing Interest
M.E. Bustamante and G. Bustamante are founders of Odinzen LLC, the company that develops the thermodynamic database reported in this work, and hold a financial interest in it. This is disclosed as a competing interest. K. Lilova declares no competing interest. The authors declare no other competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
The authors acknowledge Research Computing at Arizona State University for access to the Sol supercomputer [
43], which provided the high-performance computing resources used for the density functional theory calculations reported here. This research received no external funding.
Data and code availability
The screening, vaporization, transport, and recession programs that generate the figures and tables, together with the ten density functional pair enthalpies computed on the Sol supercomputer, are openly available at https://github.com/odinzen/rockets-public. The remaining thermodynamic data are from the cited public sources.
PDF
