{'fraction_non_alphanumeric': 0.02702777014022546, 'fraction_numerical': 0.01525982952983228, 'mean_word_length': 4.91590761223162, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Integrating Machine Learning with Mechanistic Models for Predicting the Yield Strength of High Entropy Alloys

Accelerating the design of materials with targeted properties is one of the key materials informatics tasks. The most common approach takes a data-driven motivation, where the underlying knowledge is incorporated in the form of domain-inspired input features. Machine learning (ML) models are then built to establish the input-output relationships. An alternative approach involves leveraging mechanistic models, where the domain knowledge is incorporated in a predefined functional form. These mechanistic models are meticulously formulated through observations to validate specific hypotheses, and incorporate elements of causality missing from data-driven ML approaches. In this work, we demonstrate a computational approach that integrates mechanistic models with phenomenological and ML models to rapidly predict the temperature-dependent yield strength of high entropy alloys (HEAs) that form in the single-phase face-centered cubic (FCC) structure. Our main contribution is in establishing a quantitative relationship between the HEA compositions and temperature-dependent elastic constants. This allows us to improve the treatment of elastic constant mismatch to the solid solution strengthening effect in the mechanistic model, which is important for reliable prediction of yield strength as a function of temperature in single-phase FCC-based HEAs. We accomplish this by combining Bayesian inference with ensemble ML methods. The outcome is a probability distribution of elastic constants which, when propagated through the mechanistic model, yields a prediction of temperature-dependent yield strength, along with the uncertainties. The predicted yield strength shows good agreement with published experimental data, giving us confidence in applying the developed approach for the rapid search of novel FCC-based HEAs with excellent yield strength at various temperatures.


I. INTRODUCTION
High entropy alloys (HEAs) that form in single-phase face centered cubic (FCC) structure show high ductility and fracture toughness. [1][2][3][4][5][6] However, unlike their body centered cubic (BCC) counterparts, the HEAs that form in the single-phase FCC structure have relatively low yield strength that restricts their application in high temperature environments. This shortcoming has triggered numerous experimental and computational efforts in the literature to develop strategies based on solid solution strengthening, grain size, and precipitation hardening that will improve the yield strength of FCC-based HEAs. 7,8 In strictly single-phase alloys, solid solution strengthening and grain size reduction are the two major mechanisms for hindering the dislocation motion and increasing the strength. Among the two, the solid solution strengthening effect is considered as one of the core properties of HEAs. Misfit between the atoms in the solid solution originating from the differences (or mismatch) due to atomic sizes, elastic constants, and electronic charge densities are argued as key variables in governing the solid solution strengthening effects in HEAs. [9][10][11] Fundamental computational approaches based on first principles calculations and molecular dynamics simulations have remained major drivers in the prediction of mechanical properties of HEAs. [12][13][14][15][16][17] However, these approaches are computationally expensive and cannot be used to rapidly navigate the vast HEA search space in an efficient manner. Inspired by theory and computations, there are several studies in the literature that have attempted to capture the essence of solid solution strengthening mechanism in compact mechanistic models. These mechanistic models are meticulously formulated through observations to validate specific hypotheses. Toda-Caraballo and Rivera-Díaz-del-Castillo developed a model for predicting the contribution of solid solution strengthening to the yield strength in HEAs. 18 In their work, the authors extended the Labusch's approach 19 to multicomponent alloys (including the HEAs) by developing a description of the variable unit cell parameter and atomic size misfit. Okamoto et al. uncovered a correlation between average atomic displacements and athermal yield strength (normalized by shear modulus) in three-, four-and five-component alloys. 20 Curtin and co-workers have developed mechanistic models that include the effects of defect-dislocation interaction and thermally-activated glide to predict the temperaturedependence of yield strength in HEAs. 21 This model is based on the effective-medium theory of strengthening for random alloys, which works under the assumption that each alloying element can be regarded as a "solute" in an effective "solvent" environment that represents the average alloy properties. 10 Oh et al. developed a model, where they argue that the solid solution strengthening is driven by fluctuations in the solute-dislocation interaction energy. 11 The authors correlate solid solution strengthening to atomic-level pressure, whose origin was linked to the charge transfer between the nearest neighbor atoms in the crystal structure.
In addition, there are also other mechanistic models discussed in the literature that have treated the solid solution strengthening effect in varying degrees of complexities. 8,22 As an alternative approach to using mechanistic models, many data-driven machine learning (ML) models have also been developed for predicting the yield strength of HEAs. [23][24][25][26][27][28][29] Unlike the mechanistic models, these ML models are not constrained to specific single-phase FCC and/or BCC HEAs. Thus, a vast majority of the published ML work combines all kinds of HEAs (single-phase, dual-phase, multi-phase, and mixed-phase) into one data set for establishing the quantitative composition-property relationships. While some authors have reported data on new experiments to validate the ML model predictions, others have used post hoc model explainability methods to uncover the hidden trends.
The foundation of this paper is the mechanistic model by Curtin and co-workers that was developed on the basis of the effective-medium theory. 21,30 This model takes as input the lattice constants, elastic constants, and misfit volume to predict the temperature-dependence of yield strength of single-phase FCC and BCC HEAs. All three input quantities exhibit temperature-dependence. However, it is not common in the literature to include their temperature dependence when predicting the yield strength. The objective of this work is to develop a novel computational pipeline that enables rapid prediction of temperaturedependent elastic constants of HEAs, along with the associated uncertainties. In the literature, several strategies exist to extract and predict the HEA elastic constants. [31][32][33] However, none of them has been shown to predict their temperature-dependence and quantify uncertainties.
One of the simplest and naïve approaches involve using the elastic constants of the pure elements and taking linear combination of weighted contribution of the alloying elements in the HEA alloy composition to represent HEA elastic constant. 8 While simple, any deviation from linearity will introduce error in the yield strength prediction. Therefore, uncertainty quantification is crucial to reliably interpret the results. It is also common to use first principles calculations to predict the elastic constants at 0 K 34 and use them to predict the yield strength of single-phase HEAs. More recently, surrogate ML models have also been built using the first principles calculation data to alleviate the computational cost. 31 While the surrogate models can account for uncertainty due to data sampling, they do not capture the temperature-dependence unless this data is explicitly included in the training data set.
Therefore, use of 0 K elastic constants to predict the high-temperature yield strength may not be optimal especially when the HEAs show strong temperature-dependence of elastic constants. For instance, experimental data on the shear modulus of CrMnFeCoNi HEA shows a decreasing trend as a function of temperature; it decreases from 85 GPa at 0 K to 68 GPa at 673 K. 21,[35][36][37] If we use the 0 K elastic constants to predict the high temperature yield strength, then the predictions will be severely over-estimated. In addition, if the prediction uncertainties are not quantified then the over-estimated prediction will likely mislead the HEA design. More rigorous treatment of temperature-dependent elastic constants exist within the first principles calculations based on the quasi-harmonic approximation 38 and molecular dynamics simulations 15 exists. However, these calculations are expensive.
Here, we start from the work of Varshni who proposed a phenomenological model for the temperature-dependence of elastic stiffness constants. 39 We extract the parameters of the phenomenological model using Bayes' theorem by calibrating the model with the available experimental data on temperature-dependent elastic constants for pure metals and multicomponent alloys. Ensemble support vector regression (eSVR) based ML models are then used to establish the relationship between the chemical compositions and posterior distributions of the Varshni model parameters. Finally, statistical sampling is used to propagate the predicted elastic constant distribution through the mechanistic model, which in turn predicts the temperature-dependent yield strength, along with the uncertainties, for known and unexplored FCC-based HEAs. Comparison of predicted yield strength with published experimental data shows good predictive performance, which gives us confidence in applying the developed approach to rapidly search for novel FCC-based HEAs with excellent yield strength at room and high temperatures. The developed predictive capability can be accessed via an interactive web application (https://adaptivedesign.shinyapps.io/AIRHEAD/) to guide the design novel HEAs with targeted properties.

A. Data sets
We consider three data sets for this work. Data set 1 is an experimental data set of temperature-dependent elastic constants for metals and alloys, whose crystal structures are known. We have surveyed the literature and compiled this data set for 14 materials: 12 of them are pure metals (Al, Cu, Au, Pb, Fe, Mg, Mo, Si, Ag, Ta, Th, and W), 39 one medium entropy alloy (CrCoNi 41 ) and one HEA (CrMnFeCoNi 42 ). Data set 1 calibrates the Varshni model parameters based on Bayes' theorem. Data set 2 contains experimentally measured room temperature yield strength data of as-cast FCC-based HEAs. We collected a total of 29 unique FCC-based HEAs, whose as-cast yield strength data are experimentally measured at room temperature. [43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60] Not all compositions have equiatomic concentrations of constituent elements. In addition, experimentally measured room temperature elastic constants data do not exist for the compositions in Data set 2. In Data set 2, we also did not remove the contribution from the grain size effect. Data set 3 contains experimentally measured temperature-dependent yield strength data for CrCoFeNi and CrCoFeMnNi alloys. We note that we do not use any of the compositions in Data sets 2 and 3 to calibrate the mechanistic model. Thus, the compiled yield strength data sets (Data sets 2 and 3) are solely used to validate the models and demonstrate the predictive capability of our approach. All three data sets are given in the Supplemental Material. 61

B. Varshni model
In the theory developed by Born and co-workers, it was shown that the temperaturedependence of the elastic constants originate from the variation of the lattice potential energy due to anharmonicity. 62 According to the theory, in the two extreme limits, the lattice contribution to the elastic constants should vary as T 4 at very low temperatures and as T at high temperatures. Varshni, in the year 1970, proposed a phenomenological model for explaining the variation of elastic stiffness constants with temperature. 39 The Varshni model can be written as, where C 0 ij is the elastic constant at 0 K, s ij and t ij are both fitting parameters, T is the temperature in Kelvin (K) and the subscript ij=11, 12 or 44 (since our focus is on cubic materials in this work). Empirically, for temperatures above Θ D 50 (where Θ D is the Debye temperature), Varshni showed that Equation 1 had the smallest average percentage difference with the experimental measurements. For each C ij component, we fit C 0 ij , s ij and t ij using the experimental data. While it is common to apply least squares method to extract the parameters of the Varshni model, we use Bayes' theorem for this task (discussed in Section II C) because it also allows us to quantify uncertainties in the extracted model parameters at the expense of a little extra computational cost.

C. Bayesian inference
In Bayesian inference, we use Bayes' theorem [63][64][65][66] to update the parameter distributions in the Varshni's phenomenological model given some experimental data. The mathematical expression of Bayes' theorem can be written as: where P (θ) encodes the prior probability distribution (our perceived knowledge about the model parameters) before the data is collected, P (X|θ) encodes the likelihood function that which can be calculated as follows: Since the likelihood function and prior distributions are generally known, θ can be estimated using numerical methods. In this work, we assumed Gaussian priors for the Varshni model parameters (C 0 ij , s ij and t ij ) with the following initial values: 20), and a noise ∼ N (0, 10). The mean values for the priors (C 0 lit , s lit and t lit ) are set to be equal to the least square fit from literature (when the data is available).
These would constitute an informed prior. When the least square fit is not available from the literature, the prior distribution was informed by the posteriors of those materials whose mean value are calculated. The standard deviation of priors for C 0 ij , s ij and t ij are chosen to be 50, 1, and 20, respectively. Bayesian inference with adaptive Metropolis-Hastings Markov Chain Monte Carlo (MCMC) algorithm is used to find the best match for C 0 ij , s ij and t ij from the experimental calibration data set. [67][68][69] The outcome from Bayesian inference is a posterior distribution. Our problem formulation results in producing nine posterior distributions. For each C 11 , C 12 and C 44 elastic constant, we have a set of {C 0 ij , s ij and t ij }, where ij = 11, 12, and 44. All Bayesian inference calculations are performed using the open source BayesianTools 70 package as implemented in the R-language.

D. Machine Learning
The goal of ML is to establish a quantitative relationship between the HEA composition and the C 0 ij , s ij and t ij parameters of the phenomenological Varshni model. We used the dataset (with 14 instances) that was developed for Bayesian calibration to accomplish the ML goal. The input descriptors used to represent each metal or an alloy composition in the data set are based on elemental properties. 71 Initially, we built 36 descriptors and removed redundant descriptors based on the Pearson correlation coefficient (PCC) analysis. Only those descriptors whose pairwise PCC was less than 0.4 was retained. Our PCC analysis down-selected the following three descriptors for ML model building: • "meanNdValence" is the weighted average of the number of d-electrons in the valence shell taken over all elements in the chemical composition, • "meanNValence" is the weighted average of the number of electrons in the valence shell taken over all elements in the chemical composition, and • "meanAtomicVolume" is the weighted average of the approximated atomic volume taken over all elements in the chemical composition.
We used an ensemble of support vector regression (eSVR) models with bootstrap resampling capability to learn the relationship between chemical composition and Varshni's model parameters. 72 Non-linear radial basis function (RBF) kernel was employed because of its improved generalization ability. 73,74 Each ensemble was initialized using a modified bootstrap resampling method, where a Monte Carlo sampler was used to take samples from the posterior distribution of Bayesian inference, rather than use the mean value. We explored a large number of bootstrap resamples (5,10,20,50,100,200) to define our ensemble. For each SVR model in the ensemble, a grid search was used to optimize the hyperparameters of the RBF kernel. We trained a total of nine independent eSVR models to predict C 0 ij , s ij , and t ij for each of the C 11 , C 12 and C 44 elastic constants. The optimized number of bootstrap resamples for each C 0 ij , s ij , and t ij parameter is listed in Supplementary Table S1. The final prediction from eSVR is also a probability distribution defined by the ensemble mean and standard deviation. We utilized the -SVR method as implemented in the e1071 R-package. 75

E. Mechanistic Yield Strength Model
The mechanistic model 76 for predicting the temperature-dependent yield strength of FCC-based HEAs can be written as, where the uniaxial yield stress (σ y ) is described as the product of shear stress τ y0 (T ) times M = 3.06, which is the Taylor factor for random FCC alloy. In Equation 4, k B is Boltzmann constant,ε 0 andε are reference strain rate and experimental strain rate, respectively. Thė where Γ refers to the dislocation line tension, and the line tension parameter α is set to 1/8 for FCC alloys; the δ-parameter includes information about effects of the misfit volume: For simplicity, we used Vegard's Law (i.e., linear combination of the atomic volume of the constituent elements) to approximate the misfit volume (∆V n ). We set the Burgers' vector of FCC lattice to be a 0 / √ 2, where a 0 is the lattice constant of the HEA. We approximated a 0 using Vegard's law, i.e., we used linear combination of the lattice constants of the constituent elements to approximate the a 0 of the HEA. Data for atomic volumes and lattice constants are taken from Materials Project database. 40 These quantities reflect the F m3m space group of FCC metals with some representing the ground state whereas the others in their hypothetical (or metastable) FCC structure. Thus, the misfit volume and Burgers' vector are treated as deterministic parameters, which we acknowledge is a gross approximation that can be improved. The quantity P in Equation 6 and Equation 7 is a function of elastic constants, which is informed from the combined Bayesian calibration and eSVR modeling pipeline shown in Figure 1. For FCC alloys, P can be written in terms of Voigt-averaged shear modulus (µ ave ) and Poisson's ratio (ν ave ): F. Integrated Modeling Workflow Figure 1 shows the integrated modeling workflow that combines the various approaches we have discussed thus far in the Methods section (Section II). Bayesian inference is used to calibrate the phenomenological Varshni model parameters from experimental data and quantify uncertainties in those parameters. The eSVR models are then used to establish the quantitative relationships between HEA compositions and the nine Bayesian inference calibrated model parameters (C 0 ij , s ij and t ij , where ij=11, 12 and 44). In total, nine eSVR models are built for predicting the three elastic constants: C 11 , C 12 and C 44 , along with their associated uncertainties. After establishing the predictive capability for each of the three elastic constants, the probability distribution is propagated to the mechanistic model to give the final probability distribution of temperature-dependent yield strength.

G. Global Sensitivity Analysis
To understand the relative importance of factors that controls the Varshni model and the mechanistic yield strength model, we performed the sensitivity analysis. We utilized the well-known Sobol's method, which is a variance-based method that quantifies how variance in input parameters contribute to the variance of the output. 78 We calculated the contributions due to the main effect (also known as the 1 st -order Sobol index), 2 nd -order effect, 3 rd -order effect and the total effect. 79 The expression for the 1 st -order Sobol index is given in Equation 11, where V i represents the conditional variance due to factor i and V (Y ) represents the unconditional variance of the output. This 1 st -order index shows the contribution of a single input factor i to the output variance. 80 Expression for the 2 nd -order index is given in Equation 12, where V ij indicates the variance due to the interaction of factor i with j. Thus, the 2 nd -order index can capture the variance contribution due to the interaction between two factors. Similarly, a 3 rd -order index is given in Equation 13, and it indicates the variance due to the interaction between factor i, j and k. The total effect index can be calculated using Equation 14, where E(Y |X ∼i ) is the expected value of Y when all factor, except X i , are fixed. The total effect can be interpreted as the total contribution of an input factor from both the 1 st -order effect as well as any higher-order interaction terms. 81 The Sobol sensitivity analysis is performed using the open source sensitivity package 82 as implemented in R-language.

III. RESULT AND DISCUSSION
We first use the Bayes' theorem to calibrate the Varshni model parameters (C 0 ij , s ij and t ij ) from available experimental data on temperature-dependent elastic constants. Table I shows a representative result of the Bayesian calibration procedure for the CrMnFeCoNi Cantor alloy, where we compare the outcome with the traditional least squares fit. Since the least square fits for the C 12 elastic constant is not given in the literature, we cannot compare it with the Bayesian calibration. The data shown in Table I reveals that the mean value from the Bayesian calibration fitting agrees closely with that of the least squares fitting procedure, which is expected. In addition to the mean value, the Bayesian calibration also provide the uncertainties associated with the mean value. It is important to note that the uncertainties are not readily accessible using the least squares fitting procedure. While the uncertainties for the C 0 ij and s ij parameters in the Varshni model are fairly small, we consistently find a large uncertainty for the t ij parameter. The source for relatively large uncertainty in t ij can be attributed to any of the following reasons: (1) the nature of Varshni model itself, models. From Figure 2, we infer that the eSVR predictions are relatively more accurate for the C 0 ij parameter (have high goodness-of-fit as informed by R 2 ), when compared to that of the s ij and t ij parameters. In addition, Figure 2 shows that the data is not sampled uniformly throughout the search space. Presence of non-uniform sampling indicates the eSVR models did not learn well in the regions of sparse data distribution. Thus, more data is needed in specific regions of the search space for improving the eSVR performance that opens the door for exploring the principles of adaptive materials design 83 to strategically collect data on temperature-dependent elastic constants.
The synergy between Bayesian calibration and eSVR models enable rapid prediction of the nine Varshni model parameters, along with the associated uncertainties. This then allows us to use the Equation 1 to rapidly predict the temperature-dependent C ij for any arbitrary HEA composition. Thus, for each HEA composition, we have a Gaussian distribution of C 11 , C 12 and C 44 elastic constants parameterized by mean and standard deviation. We use this distribution to randomly draw 1,000 samples of C 11 , C 12 and C 44 . We then predict the yield strength at 300 K for each of the HEA in Data set 2 using the mechanistic yield strength model (Equation 4). Since we have 1,000 C 11 , C 12 and C 44 samples for each HEA, we end up with 1,000 predictions of yield strength for each HEA composition. We then estimate the sample mean and standard deviation from these predictions. The outcome displaying the yield strength predictions (mean and standard deviation) at 300 K is shown in Figure 3. Most of the data points fall on the X=Y line, which is encouraging. Thus, for most of the as-cast HEA compositions that form in the single-phase FCC structure, the mechanistic model informed by the temperature-dependent elastic constants can capture the room temperature yield strength trend fairly well. Nonetheless, there are also outliers that appear in Figure 3, which we discuss next.
The prominent outlier in Figure 3  gives us an opportunity for validation. However, in the case of RhIrNiPdPtCu, only the room temperature yield strength data is experimentally reported. Our high temperature predictions for RhIrNiPdPtCu remain to be experimentally validated.
In Figure 4a, we compare the temperature-dependent yield strength prediction with the experimental results and predictions from Varvenne et al. 21 It is known that the mechanistic mechanistic model underestimates the yield strength at 77 K due to the theory neglecting small-scale fluctuations that contribute to the strengthening at zero-temperature and the fact that low temperature prediction is more sensitive to the line tension parameter, α. 21 However, the high temperature predictions show good agreement with the experimental data. Our yield strength prediction agrees well with the theoretical prediction of Varvenne et al 21 at both low and high temperatures. In Figure 4b, we show the predicted elastic constants as a function of temperature that was used as an input to predict the yield strength of the CrCoFeNi alloy.
All three elastic constants are predicted to show a negative slope with temperature. The C 11 and C 44 show a marginally higher temperature-dependence compared to C 12 . In the work of Varvenne et al, the authors used experimentally measured temperature-dependent shear modulus and Poisson's ratio for making the predictions. In Supplemental Table S2, we compare the experimentally measured temperature-dependent shear modulus for CrCoFeNi with that of the Bayesian-eSVR predictions from this work. We find excellent agreement between the experimentally measured and predicted values. In Figure 5(a), we compare the temperature-dependent yield strength prediction for the five-component CrCoFeMnNi alloy with the experimental data and the theoretical prediction of Varvenne et al. 21 Our predicted yield strength is marginally smaller than that of the predictions from Varvenne et al. 21 Similar to CrCoFeNi, the CrCoFeMnNi alloy also shows a decreasing slope for the elastic constants as a function of temperature (see Figure 5b). In Supplemental Table S2, we also compare the experimentally measured temperature-dependent shear modulus for CrCoFeMnNi with that of the Bayesian-eSVR predictions from this work.
Our predictions overestimate the shear modulus at all temperatures, which we attribute as one of the reasons for the marginal discrepancy in the yield strength prediction shown in Figure 5(a).
In our opinion, the DFT approach is more rigorous than the use of Vegard's law to estimate the misfit volume. In terms of the elastic constants, Yin and Curtin reported C 11 , C 12 , and C 44 as 289, 176 and 112 GPa, respectively. Our predicted temperature-dependent elastic constants is shown in Figure 6(b). Given the inherent limitations within the Varshni model (Equation 1) in accurately capturing the elastic constants trend at temperatures below Θ D /50 (see Section II B for additional information), we cannot reliably and quantitatively compare our C ij predictions with that of the 0 K DFT work from Yin and Curtin. Nonetheless, we observe qualitative agreement between our work and that of Yin and Curtin, where the C 11 > C 12 > C 44 trend is preserved.
Unlike the CrCoFeNi and CrCoFeMnNi alloys, the RhIrNiPdPtCu alloy has relatively large uncertainties in the predicted values of C 11 , C 12 and C 44 . This is because the elastic constants prediction for the RhIrNiPdPtCu alloy is an extrapolation into a new and previously unexplored composition space for our models. When compared to the CrCoFeNi and CrCoFeMnNi alloys, the RhIrNiPdPtCu alloy has the following similarities and differences in the predicted C 11 , C 12 and C 44 elastic constants: (1) The magnitude of the C 11 value is more or less the same, (2) The mean value of the C 12 is marginally higher in RhIrNiPdPtCu, (3) The mean value of the C 44 is marginally lower in RhIrNiPdPtCu. Thus, we conjecture that the origin of high yield strength in the RhIrNiPdPtCu alloy compared to the CrCoFeMnNi or CrCoFeNi alloys is likely not due to the elastic constant mismatch between the constituent atoms, but due to the atomic size (or volume) mismatch effect. The predicted low (below 300 K) and high temperature (above 300 K) yield strength data shown in Figure 6a and the temperature-dependent elastic constants data shown in Figure 6b are reported for the first time in this paper. These results await experimental validation from the community.
In order to understand the source for large uncertainty in yield strength prediction of the RhIrNiPdPtCu alloy (Figure 6a), we performed Sobol's global sensitivity analysis for the yield strength prediction at 300 K. The input parameters are the eSVR predicted C 11 , C 12 and C 44 elastic constants at 300 K. Although we assumed the C 11 , C 12 and C 44 elastic constants to be normally distributed, the Sobol analysis requires samples to be drawn from a uniform distribution U [0, 1] hypercube. We note that this simplification will not impact the generality of the findings because Sobol showed that any input space can be transformed onto unit hypercube. 79 A total number of 2 × 10 6 samples are drawn to perform the analysis.
The Sobol indices are shown in Figure 7. The total effect index is much greater than that of the main effect index for all three input parameters, which indicates that the larger variance (or error bar) observed in the yield strength predictions in Figure 6a originates from the higher-order Sobol indices (eg., interaction between C 11 , C 12 and C 44 ). The Sobol indices are also calculated for the 300 K yield strength prediction of the five-component CrCoFeMnNi and total effect (red diamond) is the S T i term in Equation 14. alloy. The result is shown in Supplementary Figure S1. We find that the Sobol indices for the CrCoFeMnNi alloy follow the same pattern as that of the RhIrNiPdPtCu alloy (shown in Figure 7). The C 11 , C 12 and C 44 elastic constants have the same degree of contribution to the variation in the yield strength predictions. Thus, the sensitivity analysis reveals that the large uncertainty in the 300 K yield strength prediction seen in the RhIrNiPdPtCu alloy compared to the CrCoFeMnNi alloy is mainly due to the large uncertainty from the eSVR predictions for the C 11 , C 12 and C 44 elastic constants.
We now use the Sobol global sensitivity analysis to understand the source of variance for the predicted C 11 , C 12 and C 44 elastic constants from the context of the Varshni model at 300 K. We use the RhIrNiPdPtCu alloy to discuss the sensitivity analysis. In Figure 8, we show the calculated 1 st -order, 2 nd -order and 3 rd -order Sobol index. We find that the variance of t ij parameter is the most dominant input factor contributing to the variance of C 11 , C 12 and C 44 elastic constants predictions, and the interaction between s ij and t ij is identified as the second dominant contributor. Thus, large uncertainties in the t ij parameters from the eSVR model are expected to contribute to the large variation in the C 11 , C 12 and C 44 predictions In an earlier work, 85 we built an interactive web application (https://adaptivedesign. shinyapps.io/AIRHEAD/), where interested readers can input an alloy composition (in the front end) and a trained ML model (in the back end) returns a prediction along with uncertainties about what phase the composition will likely form. We have now updated the web app with new sets of models discussed in the paper, where in addition to phase information the app will also output the predicted elastic constants (C 11 , C 12 and C 44 ) and yield strength at a temperature set by the user. Since the yield strength work discussed in this paper is specific to single-phase FCC HEAs and a particular mechanistic model, we have added a note in the app to inform the users about these restrictions. We believe that web apps provide an excellent avenue to share the research results and the developed predictive capability with the community. This also allows the community to objectively evaluate our approach with the literature data and create a pathway for future improvement through crowd-sourcing.

IV. SUMMARY
In this paper, we demonstrated a novel computational framework for rapidly predicting the yield strength, along with the associated uncertainties, of single-phase FCC-based HEAs with small experimental data using a mechanistic model. The key contribution is the development of a computational approach that shows the integration of a Bayesian calibrated phenomenological model with ensemble ML methods. The main outcome is the rapid prediction of temperature-dependent elastic constants of FCC-based HEAs, which is also equipped with the uncertainty quantification capability. By using these temperature-dependent elastic constants as input to the mechanistic model, we believe that we have improved the contribution from the elastic constant mismatch to the solid solution strengthening effect in single-phase FCCbased HEAs. The ultimate impact lies in more reliable prediction of temperature-dependent yield strength in the family of single-phase FCC-based HEAs. In principle, the developed approach can be extended to predict the yield strength of single-phase BCC HEAs. However, the validity of Vegard's law to capture the δ-parameter (and misfit volume) needs a more critical evaluation.
We now focus on some of the limitations of our approach, which also present future opportunities for improvement. Our temperature-dependent elastic constants predictions were based on a data set that contained 14 materials. As more temperature-dependent elastic constants data become available, the predictive capability of our model is expected to improve. In addition to elastic constants, the mechanistic model also needs the lattice constant and misfit volume as inputs. We expect both these quantities to depend on the coefficient of thermal expansion, and therefore their temperature-dependence should also be a factor that must be considered for improving the predictive capability. According to Equation 7 and Equation 9, the final expression of σ y includes the misfit volume raised to the power of 4/3, which suggests that the misfit volume is the highest order term in the analytical model. This also makes the misfit volume parameter (∆V n ) one of the most sensitive parameters among others. In this work, the misfit volume is treated as a linear combination of atomic volume from the Materials Project database, which we acknowledge is a gross approximation. Another factor that contributes to uncertainty originates from the choice of approximation used to calculate the shear modulus (µ ave ) and Poisson's ratio (ν ave ). In addition to Voigt average (the method that we use in this paper), there are also alternative methods (eg., Voigt-Reuss-Hill approximation 86 ) that may impact the various terms that appear in Equation 4. In this paper, we used the MCMC sampling method for uncertainty propagation because it is fairly straightforward to implement. However, it suffers from low convergence rate. Improving the sampling strategy is also a desired direction for future work. Thus, there are promising future opportunities that can build upon this work to develop a comprehensive Bayesian framework and predict the properties of HEAs. reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein.