The material-design strategy of mixing five or more elements in a single lattice with the random-site occupancy has brought about ✩ Notice of copyright: This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy . The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paidup, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of the federally-sponsored research in accordance with the DOE Public Access Plan ( http://energy.gov/downloads/ doe-public-access-plan ).
the recent emergence of high entropy alloys (HEAs), a new class of metallic materials [1][2][3][4][5] . Possessing a simple crystal structure with a high degree of chemical disorder, the ideal HEA presents a single-phase solid solution with exceptional engineering properties [4,[6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21] . Multi-principal element HEAs have the potential of achieving a good balance of desired properties not attained in most conventional alloys containing one or two principal elements, including the high strength [22] , reasonable ductility [23] , high hardness [24] , corrosion [25,26] and fatigue resistance [18 , 19 , 27-30] , and thermal stability [4 , 31-33] . In addition, due to the high softening resistance at high temperatures, some HEAs are promising structural materials for elevated-temperature nuclear and aerospace applications [34 , 35] .
Despite the perception that single-phase face-centered-cubic (FCC) alloys exhibit good ductility at the expense of strength [4] , the Al 0.3 CoCrFeNi HEA, a high-temperature single-phase FCC structure, shows a decent combination of strength, ductility, and toughness [15 , 16 , 36-38] . The as-cast phase and microstructure of the Al x CoCrFeNi (0 ≤ x ≤ 2.0) system were experimentally investigated by Wang et al. [39] to elucidate the effects of Al on the alloy system. For Al x CoCrFeNi where 0 ≤ x ≤ 0.5, the system forms a single FCC phase. When 0.5 < x < 0.9, Al x CoCrFeNi transforms to a twophase mixture of FCC and BCC structures, and a BCC single phase where 0.9 ≤ x ≤ 2.0. Li et al. [7] studied the correlations between the microstructure and mechanical behavior of an Al 0.3 CoCrFeNi HEA fiber fabricated by hot rotary forging and hot drawing. Their results showed a remarkable tensile strength and ductility at room temperature, and its strength and ductility even increase at cryogenic temperatures.
In this study, we employ a combination of in-situ neutron diffraction (ND), first-principles calculations, and machine learning (ML) to investigate the elastic properties of the Al 0.3 CoCrFeNi HEA, including elastic moduli and anisotropy. A reduction in the stiffness of the HEA is revealed accompanying the severe lattice distortion in the HEA. The aim of the study is two-fold. First, while there is a consensus that first-principles predictions of elastic constants are reliable for ordered structures [40] , limited results are available to draw the same conclusion for HEAs. A review by Huang et al. [41] compared ab initio methods for studying the elastic parameters of HEAs. They found that effective medium methods, such as Coherent Potential Approximation (CPA) as well as supercell methods, are generally suitable, but the effective medium model is more powerful for magnetically-disordered HEAs. Using ab initio techniques, the study also presents the dopant-concentration dependence of the elastic anisotropy [41] . The present study offers a rare opportunity to assess the agreement between experimental and computational results for the elasticity of HEAs. The second aim of this study is to understand the nature of interatomic interactions and their effects on the elastic properties of materials. With the existence of exotic structural features, such as the severe lattice distortion induced by the large-size differences in the constituent elements of the Al 0.3 CoCrFeNi HEA, it is scientifically meaningful to investigate the structure-property relationships from the perspective of the elastic response.
First-principles methods based on the density functional theory (DFT) excel in predicting the elastic properties of pure metals and ordered alloys [42 , 43] . To utilize this tool for HEAs, there is a challenging issue with constructing relatively-small supercell models that describe disordered arrangements of multiple species. This problem can be solved with the special quasi-random structure (SQS) [44] , which builds a supercell that statistically mimics the most relevant, near-neighbor pair and multi-site correlation functions of the disordered solid solution [45] . In the literature, however, there are very limited DFT investigations of the elastic properties of HEAs with experimental validations [46] . ND is employed to characterize the structural evolution of the Al 0.3 CoCrFeNi HEA during mechanical deformation. Due to a high level of penetration, even in relatively-heavy elements, in-situ ND can be used to measure the internal strain evolution of the oriented grain families in polycrystalline materials [47] . ND has been employed to study the deformation mechanism of the HEAs at both room and high temperatures [48][49][50][51][52] . For example, Huang et al. [52] studied the lattice elasticity of an FCC-structured CoCrFeMnNi by in situ ND experiments, and the HEA shows the orientation-dependent lattice elasticity. In addition to providing the information on the grain-level stress and strain heterogeneity, this technique reveals the overall role of the elastic and plastic anisotropy [53] . By integrating experiments and computations, these results provide valuable assessment of the DFT prediction of elastic constants for HEAs. In this work, we also build a ML model using the gradient-boosted trees (GB-Trees) algorithm to aid the understanding of the elastic properties of the Al 0.3 CoCrFeNi HEA.
Alloy ingots with a nominal composition of Al 0.3 CoCrFeNi were prepared by vacuum-induction melting a mixture of high-grade metals with purities greater than 99.95 weight percent (wt%). The specimen was then hot-isostatic-pressed at 1204 °C and 103 MPa for 4 h, encapsulated in an argon triple-pumped quartz tube and homogenized at 1200 °C for 2 h, air cooled, aged at 700 °C for 500 h, and finally quenched. The hot-isostatic-pressing treatment reduces casting defects such as voids which is necessary for precise mechanical test results. The homogenization treatment resulted in a homogeneous atomic distribution which was confirmed using atom probe tomography (APT). The aging treatment at 700 °C was necessary to determine the phase stability at 700 °C and removed the polycrystalline texture that was present prior to aging.
In-situ ND measurements of the cylindrical dog-bone specimen of Al 0.3 CoCrFeNi HEA (50 mm in length and 6.35 mm in diameter) were taken under tension at room temperature. An MTS load-frame on the VULCAN Engineering Diffractometer was used at the Spallation Neutron Source (SNS), Oak Ridge National Laboratory (ORNL) 54 , 55] . The ND instrument uses the time-of-flight (TOF) measurement, which allows for the ND measurements with a diffraction pattern, covering a wide range of d-spacings without the rotation of samples or detectors. VULCAN is equipped with two detectors, designated as Banks 1 and 2 at ±90 °. These detectors record diffraction patterns associated with the lattice planes parallel to the axial and transverse directions, respectively. The incident neutron beam, with a 5 mm × 5 mm slit size, illuminates the sample. A constant load-control mode with a stepwise-loading sequence was used during the measurement of the diffraction patterns. The measurement time to obtain the ND data was 10 min at each stress level up to 140 MPa. At larger stress levels, the control mode was converted from the load to displacement-control mode. The collected data were analyzed by single-peak fitting, using the VULCAN Data Reduction and Interactive Visualization software (VDRIVE) program [56] .
In the ND experiment, the HEA was uniaxially loaded in the elastic regime, and each lattice-plane ( hkl ) strain, hkl , was simultaneously measured as a function of the applied stress in both loading and transverse directions using below equation [49 , 57]
where d hkl denotes the lattice spacing of the ( hkl )-plane during deformation, and d 0 hkl is the stress-free lattice spacing of the plane. Using the Kroner model [57 , 58] , a regression is fit on , where E hkl and ν hkl are the diffraction elastic moduli and the Poisson's ratios, respectively, through setting elastic constants as parameters satisfying below equations:
where B and G hkl are the bulk and shear moduli, respectively. The single-crystal elastic constants are determined, using the leastsquares fitting over the different hkl directions, minimizing the
value of a cost function shown in below equation:
where n is the number of ( hkl ) planes applied to fit the Kroner's model;
(
First-principles calculations were performed with the Vienna Ab-initio Simulation Package (VASP) [59 , 60] using the projector augmented wave (PAW) method [61] . The exchange-correlation energy was described with the generalized gradient approximation (GGA) in the Perdew-Becke-Ernzehof (PBE) parameterization [62] . A plane-wave cutoff of 700 eV and Monkhorst-Pack k-point grid of 5 × 5 × 5 was used for all calculations. Chemical disorder was modeled with SQS [44] . The generation of SQS was based on the Monte Carlo-simulated annealing with an objective function to find the closest match of correlation functions of a disordered state [63] . The convergence of elastic constants was tested with a series of SQS with different sizes. A 64-atom SQS (4 Al, 15 Cr, 15 Fe, 15 Co, and 15 Ni) was selected for analysis in the work.
The elastic tensor was calculated, using a computational workflow based on the stress-strain method described in Ref. [40] . Starting with a relaxed SQS of the HEA, a set of distorted structures were generated using 3 × 3 Green-Lagrange strain tensors of varying magnitudes at ± 0.5% and ± 1%. For each distorted structure, the 3 × 3 stress tensor is computed by DFT. The elastic tensor of the SQS is calculated from the relationship between the stress and strain tensors:
where S ij and E ij denote the stress and strain, respectively, at row, i , and column, j , in the stress tensor. C ij are the elements in the elastic tensor with subscripts using the Voigt-notation. Each element of the elastic tensor is calculated from a linear fit of the calculated stresses over the range of strain magnitudes. Since SQS does not conserve the point-group symmetry, a projection technique is employed to approximate the elastic tensor with cubic symmetry for the Al 0.3 CoCrFeNi HEA. The projection method offers a computationally-efficient approach to calculate elastic constants for random alloys and it is shown that the approximated elastic tensor is accurate and converges quickly with the size of the SQS [64] . Convergence of the elastic constants was confirmed for the 64-atom supercell. The projected elastic tensor elements for cubic lattices are calculated, using below equation [64] :
Using the elements of the projected elastic tensor, the elastic moduli of the HEA are calculated according to Table 1 .
Two ML models using the Gradient Boosting Trees (GB-Trees) algorithm are trained on 6826 ordered inorganic compounds from the Materials Project database [66 , 68] to predict the Voigt-Reuss-Hill (VRH) averages of bulk and shear moduli [58 , 65 , 66] . To prevent the higher moduli materials from overly affecting the models, the target data was log-normalized, which also reduces skew in the distribution. The elastic constants were calculated from stress and strain relations computed using DFT and described in detail in Ref. [40] . The database contains the results of DFT calculations performed with VASP; detailed settings used for VASP are described in Ref. [66 , 68] . Our models are compared against benchmark models published in Ref. [66] which uses the Gradient Boosting Machine Local Polynomial Regression (GBM-Locfit) framework trained on 1940 ordered inorganic structures from the Materials Project database [66 , 68] . The main difference is that the models in this work use decision trees as a base learner within the gradient boosting framework, whereas the benchmark models use local polynomial regressions within the gradient-boosting framework. Models were trained using Scikit-Learn [69] .
The GB-Trees models built in this work use features generated using structural and compositional features. Properties that were calculated for each compound, such as the density, or cohesive energy per atom, are denoted as structural features. Properties that relate to elements but not compounds, such as atomic radius or group number, were combined using weighted holder means to generate descriptors denoted as compositional features.
No assumption is made about the optimum way to average the elemental properties to generate the compositional features. Hence, a collection of weighted Holder means, with powers ranging from -4 to 4, were computed for each compositional feature. The arithmetic average corresponds to a power of 1, a harmonic average corresponds to a power of -1, a geometric average corresponds to the zeroth power, and the Euclidean average corresponds to a power of 2. The effect of the power is such that the average is shifted to smaller values for smaller powers, and towards larger values for larger powers [70] . An example where the type of averaging clearly makes a difference is in estimating the density of a compound from the density of constituent elements. If the constituent elements' densities are weighted by mass fraction, then the harmonic mean more accurately estimates the compound's density as
where the terms, μ p ( x ), is the Holder mean to the power of p of the property, x, w i are the atomic weights of the elements, n is the number of constituent elements, and x i are the properties. The Yeo-Johnson's power transformation was applied to each descriptor to make the training data distribution more gaussian [71] . The data transformation accounts for differences in units, and skew in distributions.
In total 67 features were computed for each compound. For ML models to be able to predict bulk and shear moduli accurately, the models must learn patterns relating feature values to the modulus values. The more complex the underlying relationship is the more data is needed to learn complex patterns. However, learning on a dataset with unrelated, redundant, or noisy features can decrease model performance because the relevant pattern was obfuscated with data from irrelevant features. Additionally, minimizing the number of features will improve the interpretability of the model itself. The GB-Trees model makes predictions by fitting shallow decision trees, called weak learners, sequentially where each weak learner makes a prediction and the next weak learner is fit to minimize the error residual of the previous weak learner [66] . Each weak learner has high bias and low variance and by using a weighted sum of all the weak learner predictions bias is reduced while keeping variance low and accurate predictions can be made [72] . There is a bias-variance trade-off with the number of features considered; if the model contains too many features, the variance will be high, and the model may be overfit to the training data causing it to underperform when making predictions on unseen data. On the other hand, if the model contains too few features the bias will be high, and the model may be underfit, which is also undesirable [73] .
A feature subset that optimizes the model performance is generated by using a multi-objective optimized genetic algorithm feature selection methodology [74 , 75] . The genetic algorithm procedure is as follows: First, the algorithm populates a random pool of feature subsets. This pool is denoted as a generation. Feature subsets are paired up and each pair makes two new feature subsets by randomly swapping features between the two subsets. Next, all the feature subsets are sorted using a method called non-dominated sorting, in which all the feature subsets that either have the best prediction score for a given number of features, or the smallest number of features for a prediction score are selected to be included in the next generation. This selection process continues until the original generation size is achieved and the process repeats until the average model prediction score of the population converges. This non-greedy method algorithm maximizes model prediction performance and minimizes model complexity simultaneously [74 , 75] . Multiple genetic algorithm searches were performed each with a population size of 100 and most searches converged within 50 generations. A smaller feature subset improves the model interpretability and generalizability. Two separate models for predicting bulk modulus and shear modulus were made. The resulting feature subsets with used to predict the bulk ( B ) and shear ( G ) moduli are: cohesive energy, E c , density, ρ, Holder mean of the group number of atoms, μ 4 ( g ), the Holder mean of atomic radii, μ 1 ( r ), for B , and μ 2 ( r ) for G , and the Holder mean of electronegativity, μ -4 (X ) , for B , and μ 4 ( X ) for G . Table 2 shows the generated feature subsets for each GB-Trees model selected by the genetic algorithm and feature importance values.
Fitting the GB-Trees model to the training data itself is a separate task from feature selection and is explained in the following description. The model training has hyperparameters that controls how the model fits to the data. In order to prevent overfitting, a nested cross-validation scheme was used to tune hyperparameters and fit the model in separate cross-validation loops. The model fitting is done in the outer loop with 5-folds. Within each fold of the outer loop, the hyperparameters of the model are tuned using 2-folds; this cross-validation step is the inner loop. This nested cross-validation procedure is needed to prevent an overestimation of model performance which would occur if the same set of data used to tune hyperparameters was used to fit and score the model. The nested cross-validation method is used to fit the models on 80% of the total dataset. The remaining 20% was used to evaluate the model performance on an untouched test set. This ensures that the dataset used for fitting the model is disjoint from the dataset that is used to score the fitted model performance. The models were trained to minimize the mean-squared-error (MSE), with a learning rate of 0.15 (a weighting applied to each weak learner), and a maximum tree depth of 4.
To measure the lattice-strain changes with respect to the stress for the Al 0.3 CoCrFeNi HEA system, in-situ ND experiments subjected to continuous tension are performed. The 0.2% yield stress was determined as 300 MPa for the alloy. Fig. 1 (a) and (b) presents typical ND patterns of Al 0.3 CoCrFeNi under tension detected by the Banks 1 and 2 detectors, respectively at room temperature. More detailed descriptions of the ND analysis of the same sample can also be found in Ref. [76] . In addition to showing the small presence of a second phase, Fig. 1 (a) and 1(b) indicates that the sample is texture free.
The lattice strain for each ( hkl ) plane is calculated, based on the shift of the peaks during loading. The single-crystal elastic constants of the FCC phase are calculated by fitting the Kroner model to the ND peaks belonging to the FCC reflections. The aging heat treatment is necessary for the removal of polycrystalline texture which was present in the homogenized sample which is important because the determination of the single-crystal elastic constants by fitting the Kroner model to the FCC peaks has an underlying assumption of a texture free polycrystalline sample. The aged sample at 700 °C has a minor amount of secondary BCC phase (7.0% area fraction) as shown in the SEM image in Fig. 2 . More detailed descriptions of the structural and compositional characterization of the same sample can also be found in Ref. [76] .
To quantitatively and efficiently predict the elastic properties for the Al 0.3 CoCrFeNi HEA from first principles, we used a 64atom FCC SQS. The medium-sized SQS mimics the near-neighbor pair and triplet correlation functions of the five-element random solid solution. The volume, shape, and atomic positions were optimized for the FCC Al 0.3 CoCrFeNi SQS. The elastic constants of the HEA were, then, calculated with the SQS using the strain-stress method with DFT. The convergence of the elastic constants was tested carefully with respect to the plain-wave cutoff energy and k-point density. The largest difference in elastic constants between the last two k-points grids tested was 4.59 GPa for C 11 , which is less than 2% of the final value of C 11 . Table 3 gives the computational elastic constants of the Al 0.3 CoCrFeNi HEA from DFT and experimental elastic constants obtained from the in-situ ND experiments. Elastic constants calculated from the Exact Muffin-Tin Orbitals Coherent Potential Approximation (EMTO-CPA) for the same HEA [77] are also listed for comparison.
Fig. 3 illustrates the comparison of calculated and experimental values for the elastic constants of the FCC Al 0.3 CoCrFeNi HEA. As shown in Table 3 and Fig. 3 , contrary to the CPA results, the elastic constants of the Al 0.3 CoCrFeNi HEA obtained from the SQS model agree very well with those from ND experiments. The difference between DFT predictions and ND measurement is within 5% for all elastic constants. When compared with DFT predictions, CPA overestimates all elastic constants to a much higher degree. The C 44 value from CPA is 37% higher than the ND measurement.
The small difference between the elastic constants of DFT and ND experiments can be understood from the different conditions under which these results are obtained. Given the fact that the elastic moduli are affected by temperature due to lattice-thermal vibrations [78] , the slight overestimation in the DFT results, shown in Table 3 , is expected because DFT calculations are performed at 0 K, whereas ND is performed at about room temperature, 293 K. The temperature effect is also related to the slight underestimation The excellent agreement between the single-crystal elastic constants estimated by the Kroner model and DFT calculations provides the convincing evidence that a medium-sized SQS can predict the elastic properties of HEAs accurately, even with a large degree of size differences in the constituent elements and a minor amount of a second phase. When considering the intrinsic uncertainties in DFT energetics and SQS choices, the first-principles approach combining the medium-sized SQS and symmetry tensor projection still offer an efficient and reliable route to predict the elastic properties of HEAs.
To determine the effects of lattice distortion on the elastic properties, a set of first-principles elastic-constant calculations comparing a volume-only optimized SQS with atoms on the ideal lattice sites was compared to a fully-optimized SQS structure which introduces the effect of lattice distortion. The volume-optimized SQS is found the have a lower lattice constant (3.551 Å) than the fullyoptimized SQS (3.555 Å). The volume-only optimized SQS supercell is shown to have the same degree of overestimation in elastic constants as the CPA results with as shown in Fig. 4 . The bulk and shear moduli are plotted in Fig. 4 with the fully-relaxed SQS, CPA, and ML predictions.
The mechanical stability of the Al 0.3 CoCrFeNi alloy can be assessed with the Born-Huang elastic stability criteria [79] , which states a mechanically stable structure satisfies, C 11 -C 12 > 0, C 11 + 2C 12 > 0, and C 44 > 0. According to the Born's stability condition [80] , as long as the eigenvalues remain positive, the structure remains stable. Thus, the linear-stability analysis of the 4 . All the eigenvalues of the elasticity tensors are positive, meaning that the HEA can retain its structures under small deformation. Another phenomenological model, Pugh's criterion [81] predicts whether a polycrystalline material will fail in a ductile or brittle manner using the ratio of bulk and shear moduli. All values for the B/G ratio in Table 3 are greater than 1.75, suggesting that the Al 0.3 CoCrFeNi HEA is likely to behave in a ductile manner [82][83][84] . Such prediction is in good agreement with previous studies showing the Al 0.3 CoCrFeNi HEA with the superior ductility [8] . Moreover, the Cauchy's pressure, C 12 -C 44 , is commonly related to the ductility behavior of materials, as it infers the character of the atomic bonding in metallic materials. Where metallic bonding corresponds to positive Cauchy pressure, and bonding with an angular character corresponds to negative Cauchy pressure [85] . Contrary to the CPA result, the Cauchy pressures obtained from the DFT calculations are positive, revealing that Al 0.3 CoCrFeNi shows a strong metallic bond and ductility behavior. This feature is also consistent with the experimental analysis [36 , 86] .
The lattice-strain response of the (111), ( 200), (311), (331), and (220) planes, with respect to the applied tensile stress, are plotted in Fig. 5 . Description of the lattice strain evolution of the same sample can also be found in Ref. [76] . Each line indicates the response of grains with a lattice direction corresponding to the line legend. The various slopes of the curves, within the elastic region, relates to the elastic anisotropy of different ( hkl ) grain families in the FCC phase [87] . For the (331) grain family, the plotted loading is in the linear-elastic region. The curve of the (200) lattice grains, however, exhibits a slight downward deflection at about 350 MPa, signifying an elastic-to-plastic transition, and load being shared from the yielding FCC grains. The curves for the (111) and (220) lattice grains present a slight upward deflection at about 40 0 MPa. The (20 0) and (111) lattice planes correspond to the extremes of the elastic stiffness in elastically-anisotropic FCC materials, as it can be seen for Al 0.3 CoCrFeNi in the figure, i.e., the (200) and (111) grain orientations possess the low and high directional strength-to-stiffness ratios, respectively [88] . As a result, the ( 200) and ( 111) orientations are the first and the last to yield, respectively. The (200) grains cannot take on more stress whereas the (111) grains take on more stress, which explains the behavior of the stress-lattice responses in Fig. 5 . Fig. 5 shows lattice strains in loading and transverse directions obtained from grains oriented such that the individual ( hkl ) plane normals are parallel and perpendicular, respectively, to the load axis. It can be observed that the lattice strain changes in the transverse direction are significantly smaller than those in the axial direction. The response of the individual peaks in the transverse direction is almost linear, with a slope opposite in sign to that obtained in the load direction, up to about 300 MPa, where the yield of preferentially-oriented grains starts. Some of the reflections, in the transverse direction, show more lattice strains than others, which is consistent with the differences in lattice strains in the loading direction, presented in Fig. 5 .
The description of spatial dependency of the elastic properties of the Al 0.3 CoCrFeNi HEA in Fig. 6 is derived from the DFTcalculated single-crystal elastic constants. In the three-dimensional (3D) representation, properties for an elastically-isotropic systems present a spherical shape, and any deviation from a spherical geometry suggests the presence of anisotropy [89] . Table 5 and Fig. 6 reflect vast anisotropic elastic properties of the Al 0.3 CoCrFeNi HEA system.
The extent of anisotropy and the variation of the elastic moduli with the directions of the HEA can be demonstrated by the degree of deformation of the sphere. The shear modulus and Poisson's ratio depend on both the direction of the applied stress and a perpendicular measurement direction [90] , which is represented in the figures by maximum and minimum surfaces. In certain directions, the Poisson's ratio is negative. Fig. 6 (b), (d), and (f) illustrates the 3D surfaces of the crystallographic-orientationdependent Young's and shear moduli, and Poisson's ratios for all possible crystallographic directions to show the weak and stiff directions in the crystal structure. In Fig. 6 (d), the transparent green outer surface and the solid magenta-yellow inner surface represent the positive maximum and minimum values of shear moduli, respectively. In addition, Fig. 6 (f) displays three surfaces; the outer transparent green, the middle transparent magenta-yellow, and the inner solid black surfaces represent the positive maximum, positive minimum, and negative values of Poisson's ratios, respectively. Fig. 6 (a), (c), and (e) represents the two-dimensional projections onto the XY, XZ , and YZ planes of Young's and shear moduli and Poisson's ratios. The magnitude of Young's modulus in a specific direction illustrates the strength of chemical bonds in that direction, suggesting that the Al 0.3 CoCrFeNi HEA is more demanding to be stretched along the body-diagonal directions than the axial ones. The macroscopic Young's modulus calculated from Table 3 has a value of 210.0 GPa. This modulus is close to the average of the maximum and minimum values in the [001] and [111] orientations, which is 225.45 GPa. This trend may be expected for a polycrystalline material without a crystallographic texture. The shear modulus of Al 0.3 CoCrFeNi shows an opposite trend with respect to the Young's modulus, illustrated in Fig. 6 (c) and(d). The shear modulus remarkably depends on the stress direction, and the Poisson's ratio has similar characteristics. The lowest and highest orientations for shear moduli are [111] (45.9 GPa) and [001] (135.5 GPa), respectively. The opposite trends of Young's and shear moduli in [001] and [111] contribute substantially to the large anisotropy of the cubic crystal, which is ascribed to the electronic structures and atomic arrangements. Lastly, the 3D contour surface and the projection representations of the Poisson's ratio for the Al 0.3 CoCrFeNi HEA crystal are shown in Fig. 6 (e) and (f), respectively. The Poisson's ratio of the crystal Al 0.3 CoCrFeNi is maximized, ν = 0.728, as the stretching along the diagonal directions, [100], induces the lateral contraction along axial directions, [010]. Al 0.3 CoCrFeNi also exhibits a negative Poisson's ratio that is minimized, ν = -0.155, in the < 110 > directions. The appearance of negative Poisson's ratios along certain crystallographic directions in the HEA is in line with various studies of materials with the large anisotropy [91][92][93] . Lethbridge et al. surveyed the experimental elastic constants of 472 materials and found a correlation between the extreme Poisson's ratios and elastic anisotropy, irrespective of the crystal symmetry [91] . Also, a study by Wang et al. determined that compounds with negative Poisson's ratios had metallic bonding as opposed to compounds that did not have negative Poisson's ratios and had bonds of more covalent nature [93] .
Using ML models to predict material properties has the advantage of being extremely fast, compared to computing the property values using DFT. In the present work, a GB-Trees algorithm is used to build a model to predict the bulk and shear moduli. The features of the model are subsequently analyzed to gain insights into the dependence of the predictions on feature values. The ML model in the current work is benchmarked against a model published in Ref. [66] , the GBM-Locfit model, which uses local-polynomial regression functions for base-learners, as opposed to using decision trees. The GB-Trees model has improved predictions for the bulk modulus, compared to the benchmark GBM-Locfit model. The experimentally measured and predicted VRH averages of the elastic bulk and shear moduli are listed in Table 5 . The VRH averages of elastic moduli are averages of the Voigt and Reuss average moduli for polycrystalline materials, which correspond to upper and lower bounds, respectively. The Al 0.3 CoCrFeNi bulk modulus prediction by the ML model created in the present work is improved, compared to the bulk modulus prediction by the GBM-Locfit model. The prediction of shear modulus is similar in both models. For the Al 0.3 CoCrFeNi HEA, the GB-Trees model predicted a bulk modulus of 161 GPa, which has a relative error of 6%, which is an improvement on the 208 GPa benchmark prediction with a relative error of 20%. For the Al 0.3 CoCrFeNi HEA the GB-Trees model predicted a shear modulus of 73 GPa, which has a relative error of 10%, and the benchmark predicted 85 GPa, which has a relative error of 5% as shown in Table 6 .
We note that the training datasets for the GB-Trees model and the benchmark models were different and that neither contains random alloy or HEA data. To make more equal comparisons between the two models learning curves are used to evaluate the dependence of the model learning performance based on the size of the training dataset. In Fig. 7 Considering that the models are trained with no HEA data in the dataset, the favorable model performance indicates that the model has good generalizability and can serve as a rapid screening tool to make predictions for compounds and alloys without the need to do further DFT elastic constant calculations. The bulk and shear moduli of compounds can be predicted with the GB-Trees models without having to perform DFT elastic-constant calculations. The models made in this work have relative errors of less than 5% for bulk and less than 10% for shear moduli and for over half of the testing dataset compounds. It may be possible to further reduce the relative errors of predictions for disordered alloys like HEA's by introducing alloys with disorder, and lattice distortion into the training dataset, and also, include descriptors for disorder, lattice distortion [66] .
The effect of interatomic interactions on mechanical properties of the Al 0.3 CoCrFeNi HEA is elucidated by examining nearest neighbor distances and magnetic moments of elemental species. Fig. 8 shows the distribution of the interatomic distances of the nearest-neighbor atoms for the FCC Al 0.3 CoCrFeNi HEA calculated from the optimized 64-atom SQS. The SQS does not contain any nearest-neighbor Al-Al pair, but it represents a population that allows probing the statistical distribution of the nearest-neighbor distances for other homoatomic and heteroatomic pairs in the HEA. The wide spread of interatomic distances, even between atoms of the same species, such as Cr, is rare in conventional alloys. We note that the interatomic distances of Cr-Cr pairs are related to the local environment of the Cr atoms, which also strongly affects the magnetic moments of Cr atoms. The magnetic frustration of Cr due to its antiferromagnetism has been observed in other HEA's such as CrMnFeCoNi HEA and its quaternary derivative of NiFe-CrCo [94 , 95] . Given a local environment with mostly-ferromagnetic nearest neighbors (such as Fe, Co, and Ni), another Cr atom as a nearest neighbor results in geometric frustration. The bond length of a Cr-Cr nearest neighbor pair is shown in Fig. 9 to be affected by its local environment, and the degree of magnetic frustration. This indicates that the magnetic property of constituent elements can also have a large impact on lattice distortion. When the Cr atoms are surrounded by atoms with a large positive magnetic moment, it can attain a large negative magnetic moment. However, when there are more Cr nearest neighbors, they become magnetically frustrated and tend to have magnetic moments closer to zero. The red circles around the Cr data points in Fig. 9 (a) and (b) represents those atoms with an Al atom in its nearest neighbor.
There are more Cr atoms with frustrated magnetic moments in the unrelaxed supercell, compared to the relaxed supercell.
Fig. 10 plots the magnetic-moment difference between the Cr nearest neighbors against the Cr nearest neighbors' interatomic distance for both the relaxed SQS supercell (red dots) and the unrelaxed SQS supercell (blue dots). The number of red and blue dots is the same, but the blue dots are heavily overlapped. In the relaxed supercell, there seems to be a relation that the more antiparallel Cr-Cr pairs have smaller interatomic distances, compared to the magnetically-frustrated Cr nearest neighbors. As previously mentioned, the magnetic moment of the Cr atom is heavily influenced by the local environment. In the case of the magneticallyfrustrated Cr atoms, the presence of an Al atom instead of another Fe, Co, or Ni atom will reduce the average magnetic moment of the nearest neighbors which may affect the bond lengths and lattice distortion as well.
This study compares the single-crystal elastic constants of an HEA with the composition of Al 0.3 CoCrFeNi, obtained by ND experiments, first-principles calculations, and ML model precdictions. Obtaining single crystal elastic constants requires sufficiently large single-crystals that are homogeneous and defect-free. This requirement makes the availability of the single-crystal elastic constant data for HEAs sparse. The integrated approach of this study provides the valuable experimental validation of elastic-constants values calculated using first-principles methods. Using the Kroner model [57 , 96] , the single-crystal elastic constants are estimated from diffraction elastic constants of a polycrystalline Al 0.3 CoCrFeNi HEA sample. Comparing the estimated single-crystal elastic constants with the results from different ab initio methods, it is found that the SQS supercell method with relaxed atomic positions has elastic-constant values that are within 5% of those from the ND experiments. The CPA model and the SQS model without the atomic position relaxation both overestimated elastic-constant values. Within CPA, the random alloys are treated by assuming the average occupations of lattice sites and do not consider lattice distortion [97] . The ideal-lattice assumption made in CPA is not supported by the experimental observations in which it is demonstrated that the interatomic distances of the atomic pairs in HEAs are essentially different [98] . In addition, CPA also does not specifically consider the charge transfer between atoms, which can lead to errors in describing the interatomic bonding. In the SQS model, the most relevant local states of the disordered phases are captured [45] . The lattice distortion results in large distributions in interatomic bond lengths for all the different atomic specie pairs. The largest difference in the atomic radii is between an Al-Cr pair with a difference of 0.15 angstrom [99] . It was expected that the elastic moduli of the HEA would be increased due to the strain the interatomic bonds are experiencing within the distorted lattice. However, the comparison between the volume-only optimized SQS supercell and the fully-optimized SQS supercell revealed a decrease in the elastic moduli. The lattice constant of the fully-relaxed SQS model was larger than that of the volume-optimized SQS model, suggesting that the severe lattice distortion in HEAs can reduce the average interatomic bond strength and reduce stiffness as in the case of the FCC Al 0.3 CoCrFeNi HEA. Some of the reduction in stiffness from the volume-optimized SQS and the fully-optimized SQS may be attributed to the internal relaxations of the atoms; as explained by a relaxation term in the elastic tensor calculation defined by Lutsko [100] . The effect of the relaxation term and the lattice distortion effect are coupled since the relaxation term is a direct result of internal heterogeneity i.e. atomic size differences; which gives support to the prospect of tuning elastic properties via lattice distortion engineering. In addition to atomic size differences, magnetic properties are another factor to consider with respect to lattice distortion effects. For example, the Cr nearest neighbor bond lengths are affected by the magnetic moments of other atomic species in the local environment. Toda-Caraballo et al. proposed a methodology that can compute the distribution of interatomic distances of HEAs with the input of unit-cell parameters and bulk moduli of elements involved [101] . Here, our study provides the concrete evidence to corroborate their argument on the linkage between the lattice distortion and elastic properties of HEAs. Lee et al. designed a single-phase BCC solid-solution phase refractory HEA with exceptional yield strength as well as ductility which was attributed to solid-solution hardening originating in turn from the lattice distortion in the refractory HEA [102] . The discrepancy in the elastic constant of the FCC Al 0.3 CoCrFeNi HEA between the SQS supercells (one with the atomic position relaxation and one without) suggests the importance of such lattice distortion in determining the elastic properties of HEAs.
The ML models are leveraged here to uncover patterns between features and the bulk and shear moduli learned by the models. Interpretable ML models are useful in gaining insight into the underlying physics that governs the observed behavior, and may serve as guides for adjusting compositions in designing HEAs with desirable characteristics. To quantify how much the model predictions are influenced due to a change in a given feature value, feature importances are computed. In this work, we use Shapley Additive Explanations (SHAP) values to measure feature importance [103 , 104] . The SHAP value of a feature represents is the average marginal contribution of a feature value. For example, given a modulus to predict for a single sample, the SHAP method takes the difference between the predicted modulus and the average of all the modulus predictions in the training dataset. Then the method determines how much each feature affects the prediction. The sum of all the contributions of each feature equals the difference between the modulus prediction of the current instance and the average modulus prediction of all instances. Figs. 11 and12 show scatterplots of the SHAP values for all 5 features and for the bulk and shear moduli-prediction models, respectively. Each dot in a scatterplot represents a single prediction instance. The y -axis shows the SHAP value assigned to that feature in that particular instance, and the corresponding feature value can be read on the x -axis. A SHAP value of 0 means that for that instance which had a given feature value, that feature did not contribute to a change in the model prediction from the average prediction. Positive or negative SHAP values on the other hand, mean that the given feature value made the model prediction go up or down, respectively. The scatter points are also mapped to a color bar corresponding to a second feature so that the interaction of pairs of features and their effect on model predictions can be observed. In both the bulk and shear modulus prediction models, an increase in moduli generally corresponds to a decrease in cohesive energy, the average group number of the constituent elements, or the average atomic radus of the constituent elements. Conversely, an increase in both bulk and shear modulus predictions are seen for an increase in density. The average electronegativity of the constituent elements, however, has a divergent behavior between the bulk and shear modulus predictions. The bulk modulus predictions increases with increasing average electronegativity, hitting a peak. With further increased average electronegativity, the bulk modulus decreases. The shear modulus predictions, however, decrease with increasing average electronegativity. It can be pointed out that the power of the Holder mean of electronegativity for the bulk modulus predictions is -4, whereas it is 4 for the shear modulus predictions. The effect of this is that the bulk modulus predictions are more dependent on the electronegativity of the most electronegative element, and the shear modulus predictions are more dependent on the electronegativity of the least electronegative element. These general trends can aid the development of new design rules similar to how empirical rules have been developed historically in the materials science and metallurgy i.e., Hume-Rothery rules. The ML work done in this study, illustrates how the ML-model selection and interpretation techniques can be generalized on a larger scale to understand underlying relationships between properties.
The feature subsets for the bulk and shear modulus prediction models are listed in Table 2 with the computed global SHAP values which are computed by averaging the absolute values of all the individual SHAP values from each of the training samples. For both bulk and shear moduli predictions, the group number and cohesive energy features are the largest contributors to the model output indicating that changes in these features cause the larges changes in the model prediction.
The predictions from the ML models can be combined with empirical models to predict materials properties that are difficult to predict from first-principles, for example, the hardness of the Al x CrCoFeNi HEA. Experimental studies on the effect of varying the Al content in the Al x CrCoFeNi HEA on its mechanical properties have found that for the as-cast samples with the FCC phase, Vicker's hardness increases with increasing the Al content from x = 0 to x = 0.5 [105] . Vicker's hardness can be predicted, using known shear and bulk moduli, and Eq. ( 12) published by Niu et al., which takes both shear and bulk moduli into account and has good agreement for both cubic and non-cubic materials [83 , 106 , 107] :
where H v is Vicker's hardness. The increase in the predicted bulk and shear moduli, and the predicted hardness is positively correlated with the experimentally-observed increase in hardness, as shown in Table 7 . The ML models for predicting bulk and shear moduli has its own errors, and using those predictions in the above model can propagate errors. However, this qualitative agreement in trends between the predicted and measured hardness values illustrates the ease and utility of using the ML model to study trends and feature dependence of model outputs either for design or screening purposes.
In conclusion, the work explored the elasticity of FCC Al 0.3 CoCrFeNi, an HEA with high strength and good ductility, with experimental and computational methods. Good agreement in elastic constants was observed from ND, first-principles calculations, and ML models. Anaysis of the computational results found that a severe lattice distortion presents in this HEA, leading to an reduction in its stiffness. A high degree of eastic anisotropy was revealed through the first-principles determination of the elastic constants, which also agreed with the observed lattice strain evolution during the in-situ ND tension experiments. It is demonstrated that ML models can cpature the correlations between material features and make accurate predictions of elastic constants of the HEA. This study on the elastic constants of the Al 0.3 CoCrFeNi HEA provides a much-needed experimental validation of computational tools and encourages more studies using this integrated approach to actively guide the search for new alloy design. The ML models can be interpreted either for developing new design rules or understanding relations between physical quantities of interest.