Theoretical prediction of high melting temperature for a Mo–Ru–Ta–W HCP multiprincipal element alloy

While rhenium is an ideal material for rapid thermal cycling applications under high temperatures, such as rocket engine nozzles, its high cost limits its widespread use and prompts an exploration of viable cost-effective substitutes. In prior work, we identified a promising pool of candidate substitute alloys consisting of Mo, Ru, Ta, and W. In this work we demonstrate, based on density functional theory melting temperature calculations, that one of the candidates, Mo0.292Ru0.555Ta0.031W0.122, exhibits a high melting temperature (around 2626 K), thus supporting its use in high-temperature applications.


INTRODUCTION
Among refractory metals, rhenium demonstrates a unique combination of desirable thermodynamic and mechanical properties: a high melting point (3448 K), excellent high-temperature strength, a good ablation and creep resistance, and no ductile-tobrittle transition 1,2 . The concurrence of low-temperature ductility and high-temperature strength makes Re-based alloys competitive candidates for structural materials with applications such as rocket engine nozzles, which involve rapid thermal cycling under high temperatures. However, rhenium is a rare element and its high cost limits its widespread use 3,4 . This situation has prompted the exploration of practical and efficient replacement for rhenium in rocket engine nozzle applications.
Previously, we showed that a straightforward design principle that trades off average valence electron count and cost considerations is useful in recognizing a favorable group of possible substitutes, the Mo-Ru-Ta-W quaternary alloys 5 , based on a thermodynamic model that integrates electronic structure calculations with the Calphad framework 6 , as well as through experimental synthesis and structural characterization of samples chosen in a favorable composition range. The computational thermodynamic model enabled the identification of alloy compositions on a so-called Re-equivalent plane that also preserve rhenium's hcp structure, and thus its desirable mechanical properties 7 . In other words, based on a rigid-band model, promising alloys have a composition on a Re-equivalent plane on which the average number of valence electrons per atom equals that of Re. The accuracy of this simple picture was assessed through high-throughput ab initio thermodynamic calculations, and with comparisons to established binary phase diagrams sections, and further corroborated by experimental synthesis and structural characterization demonstrating single-phase multiprincipal-element hcp solid-solution samples selected from a promising composition range. This region includes the Mo-Ru-Ta-W quaternary alloys having an hcp structure at compositions near Mo 0.3 Ru 0.512 Ta 0.064 W 0.125 and near Mo 0.292 Ru 0.555 Ta 0.031 W 0.122 , which straddle the Re-equivalent plane. Although Ru is also somewhat expensive, including it allows us to add more of the inexpensive elements such as Mo and W, thus still resulting in a net cost reduction. These specific compositions were also synthesized experimentally and their mechanical properties investigated 5,8 . We also showed that the proposed Ru-based alloys exhibit better thermal compatibility with Ir, which serves as an oxidation-resistant coating. 9 . However, it remains to be demonstrated that this proposed alloy exhibits a sufficiently high melting point for its intended application. An experimental determination of the melting point is complicated by the high temperatures involved. These alloys may have the potential to replace higher cost elements, like Re, in relatively lowertemperature thruster applications where cost is a major driver but performance, such as ablation and erosion, is adequate.
In this work, we employ density functional theory (DFT) [10][11][12] to calculate the melting temperature of Mo 0.292 Ru 0.555 Ta 0.031 W 0.122 . For accuracy assessment and improvement, we also compute the melting points of its constituting elements. We use an efficient small-cell coexistence method 13 and its implementation in the SLUSCHI code 14 , based on DFT molecular dynamics (MD), as shown in Fig. 1. This highly cost-effective and robust method allows us to perform expensive melting point calculations on multiple materials directly from DFT.
As shown in Table 1, our calculations suggest the need to use accurate pseudo-potentials where semicore elections s and p states were treated as valence states. With only valence elections, the DFT melting point errors were 210, 170, 330, and 200 K for Mo, Ru, Ta, and W, respectively. When semicore p (including s for W) electrons were further relaxed, the errors were significantly reduced to 80, 60, 100, and 170 K, respectively. The high accuracy of the latter served as strong indication that our DFT melting point calculation model will also accurately generate melting properties of Mo 0.292 Ru 0.555 Ta 0.031 W 0.122 , a promising Re-substitute alloy.
The hcp solid structure of Mo 0.292 Ru 0.555 Ta 0.031 W 0.122 was represented by a special quasi-random structure (SQS) 15  A recent study suggests that, due to underbinding, GGA may systematically underestimate melting point and thus GGA value may serve as a lower boundary of melting point 17 . We here make a melting temperature correction on Mo 0.292 Ru 0.555 Ta 0.031 W 0.122 , based on the known melting temperature errors of its alloying elements in pure form, where x i is the composition of each alloying element, and ΔT i is the corresponding DFT melting point error. The correction is 81 K, and the melting temperature of Mo 0.292 Ru 0.555 Ta 0.031 W 0.122 is thus estimated to be 2626 K. This is (slightly) lower than the weighted average of the elemental melting points (2850 K), thus suggesting a eutectic topology. While this melting temperature of 2626 K is significantly lower than that of rhenium (3459 K), this is not expected to be a limiting factor in applications where rhenium is used in conjunction with an iridium coating to improve oxidation resistance. Iridium's melting point is only 2719 K and thus becomes the more relevant limiting factor. Our proposed alloy's melting point would thus only lower the maximum operating temperature by about 100 K. Note that a typical rhenium/iridium combustion chamber operates at temperatures from 1925 up to 2200°C (2473 K) 18 , a temperature range our alloy should still be able to withstand without melting.
We should point out that the small-size coexistence method is inherently capable of finding the correct solid crystal structure, because the solid-liquid interface facilitates the nucleation of the structure with the lowest free energy even if we initially assumed the wrong solid. In other words, the method itself serves as a corroboration to verify the correct solid structure. In our simulation, we confirm that Mo 0.292 Ru 0.555 Ta 0.031 W 0.122 remains in hcp structure, thus further corroborating our previous computational and experimental findings 5 .
In addition to Mo-Ru-Ta-W alloys, we also show that Ru-Re-Ta and Ru-Re-W alloys also have favorable high melting temperatures, as summarized in Table 1, while also remaining hcp structures. These results indicate that attempting to increase the alloy's melting point by replacing the element with the lowest melting point (Mo) by Re does appear to be an effective strategy. It should be emphasized that the melting point reported here formally corresponds to the temperature, where the liquid and solid-free energies cross at a given fixed composition. This is the case because, over the timescale of our simulations, the species do not have the time to segrate from one phase to another, which can be readily verified by inspecting the final composition profile of the runs which terminate in a fully solid state. The real system may melt incongruently (with the solid and liquid having different compositions and the melting transformation taking place over a range of temperatures.) Our reported temperature lies in between the liquidus and the solidus, by construction.
In this system, the liquidus and solidus are likely very close to each other and thus close to our reported temperature. This is supported by the facts that (i) the alloy's main constituent (Ru) has a similar melting point (2550 K) as our composition of interest; (ii) additional calculations to estimate the melting at a composition near the middle point between Mo 0.292 Ru 0.555 Ta 0.031 W 0.122 (m.p. 2542 K) and pure Ru suggest a melting point approximately 2500 K. Therefore, the liquidus and solidus are indeed likely to both be very flat near the composition range of interest.
To summarize, we demonstrate, based on DFT calculations, that a potentially promising rhenium substitute alloy, Mo 0.292 Ru 0.555 -Ta 0.031 W 0.122 , exhibits a very high melting temperature of 2626 K, thus supporting its high-temperature applications. This prediction uses the known errors in the calculated melting temperature of the constituting elements to minimize any systematic bias. The calculations also corroborate our previous findings from both computation and experiment that Mo 0.292 Ru 0.555 Ta 0.031 W 0.122 has an hcp structure.

Density functional theory melting point calculations
We use an efficient small-cell extension of the coexistence method 13 and its implementation in the SLUSCHI code 14 , based on DFT MD. This highly efficient method makes it possible to perform, directly from first principles, expensive melting point calculations. The method runs solid-liquid coexisting simulations on small-size systems, and the melting temperatures are rigorously inferred based on statistical analysis of the system's fluctuations, namely, the temperaturedependence of the probability that simulations terminate with the system in a fully liquid state. This probability can be calculated analytically as a function of unknown thermodynamic parameters (including the melting point) that can then be determined from a fit to the observed temperature-dependent probabilities 13 . The accuracy (typically with an error smaller than 100 K), robustness and efficiency of the method have been demonstrated in a range of materials 13,14,19-24 . In particular, the small-cell coexistence method 13 and the SLUSCHI code 14 was employed to computationally predict the material with the highest melting point, which was subsequently confirmed by independent experiments [25][26][27] . DFT calculations were performed by the Vienna Ab initio Simulation Package (VASP) 28 , with the projector-augmented-wave (PAW) 29 implementation and the generalized gradient approximation (GGA) for exchange-correlation energy, in the form known as Perdew, Burke, and Ernzerhof 30 . Since the simulations were performed at hightemperature conditions, we used accurate pseudo-potentials where the semicore s and p states were treated as valence states.

DATA AVAILABILITY
All data generated or analyzed during this study are included in this published article.
Received: 2 July 2020; Accepted: 28 November 2020; 1 PAW and PBE were employed. "v", "pv", "sv" and dates denote different pseudo-potentials. "v" stands for valence elections only, while "pv" means semicore elections p states were also treated as valence states. For W, even s semicore elections were relaxed. Supercell size was approximately 20 × 10 × 10 Å 3 in all calculations. Default energy cutoffs were set and Pulay stress was included. 2 Title, date and type do not uniquely characterize a pseudo-potential. Here we identify a pseudo-potential with a 32-bit CRC function of the POTCAR 31 . 3 Number of atoms in the supercell. 4 Total number of MD trajectories sampled. 5 Physical time spent. 6 A special Baldereschi mean-value point in the Brillouin zone, (1/4,1/4,1/4). 7 Automatic generation of k-points (see VASP manual).