Ab initio aided design of novel quaternary, quinary and senary high-entropy borocarbides

High-entropy materials have attracted considerable interest due to their unique, improved properties and large configurational entropy. Out of these, high-entropy ceramics (HECs) are of particular interest since the independent solubility of cations and anions results in increased configurational entropy. However, most HEC research considers only a single element occupying the anion sublattice, which limits the maximum attainable configurational entropy. Here, we expand our previous work on high-entropy borocarbides where both boron and carbon occupy the anion sublattice. By applying an ab initio based screening procedure, we identify six elements Li, Ti, V, Zr, Nb and Hf suitable for forming high-entropy borocarbides. With these elements, we propose six novel HEC compositions, and by computing their entropy forming ability, we identify that three are likely to form single-phase during synthesis. Material properties and lattice distortions for all proposed compositions are studied using density functional theory calculations with special quasirandom structures. The directional lattice distortions, a concept we introduce in this work, show that lattice distortions have an elemental and directional preference for certain HEC compositions. We also show that the novel inclusion of Li improves the mechanical properties of the proposed HECs, the details of which are studied using the electron localization function.


Introduction
High-entropy ceramics (HECs) are synthesized from a mixture of four or more multicomponent ceramic compounds following the concept of configurational entropy stabilization of a single-phase [1,2]. The crystal structure of HECs consists of a cation sublattice and anion sublattice [3], which offers independent solubility of cations and anions in one crystal structure, enhancing the configurational entropy [4]. The first reported HECs were high-entropy nitrides [5] and high-entropy oxides [6], followed by highentropy borides [7] and high-entropy carbides [8,9]. HECs have attracted attention in recent years due to their excellent chemical and mechanical properties such as good oxidation resistance [10], enhanced hardness [11][12][13] and superior thermal properties [14][15][16].
Computational studies have contributed toward the discovery and understanding of HECs through the screening of promising compositions and prediction of their properties, such as the development of the entropy forming ability (EFA) descriptor by Sarker et al. which can be used to predict whether a high-entropy composition is likely to form singlephase during synthesis [9]. The EFA descriptor is obtained via ab initio density functional theory (DFT) calculations, and was used by both Sarker et al. [9] and Harrington et al. [17] to explore the experimental synthesizability of many different quinary high-entropy carbides.
Recently, Kaufmann et al. [18] proposed a powerful machine learning (ML) method which uses both chemical attributes and thermodynamic data obtained via calculation of phase diagrams (CAL-PHAD) to predict EFA values for high-entropy carbide compositions. Their ML method show good accuracy compared to DFT calculations and was used to predict the synthesizability of 70 new high-entropy carbides containing chromium. CALPHAD, which is a semi-empirical calculation method, has itself been used to predict phases of high-entropy alloys [19][20][21][22] and in some cases even HECs [23]. However, CAL-PHAD relies on thermodynamic databases and requires expansion to include more possible compositions for border use.
In addition to its use for predicting HEC candidates, DFT has also been used to study their mechanical, electronic, and thermodynamic properties. Wang et al. [24] studied the high-entropy boride (TiZrHfTaM)B 2 system where M=Nb, Mo or Cr using DFT and showed that the highest strength and hardness was found for M=Nb. As a result of its stronger interatomic interactions, it also increases the thermal conductivity compared to the other compositions. Furthermore, Yang et al. used DFT to study the structural, mechanical and electronic properties of high-entropy (TiZrNbHfTa)C under pressure [25], and followed this up with a detailed study of the thermodynamic properties of the high-entropy (TiZrNbHfTa)X system [26], where X=C or N. They found that X=N increases the stability at high temperatures due to the larger entropy contributions, while X=C results in higher hardness due to the stronger chemical bonds. Apart from the quinary high-entropy carbides mentioned above, Jiang et al. [27] used DFT to study the quaternary high-entropy (ZrHfTaM)C system, where M=Ti or Nb. They found that M=Nb results in improved mechanical properties but also increased brittleness. A few studies even used classical molecular dynamics simulations to model the thermal conductivity of high-entropy oxides [28][29][30], where it was found that mass disorder alone cannot explain the low thermal conductivity observed in high-entropy oxides. Instead, the disorder in the interatomic forces must be considered as well. With some exceptions [31,32], most experimental and certainly most theoretical works have focused on HECs with only one element occupying the anion sublattice. However, with the inclusion of two or more elements occupying the anion sublattice, the configurational entropy significantly increases. In this paper, we expand on our previous works on the high-entropy borocarbide system [33][34][35] for which both boron and carbon atoms occupy the anion sublattice.
Using DFT, we screen elements from group 4 to 14 of the periodic table to find suitable candidates which can form high-entropy borocarbides. Screening was performed using three different systems based on the hexagonal LiBC unit cell (P6 3 /mmc). This crystal structure was chosen based on our previous experimental results where the synthesized high-entropy borocarbides show a similar structure [33][34][35]. These three systems represent different compositions in the anion sublattice, an equal ratio of boron and carbon, Fig. 1b, and the two extreme cases of only boron, Fig. 1a, or only carbon, Fig. 1c. It should be noted that the MC 2 structure is purely hypothetical and is used J Mater Sci (2022) 57:422-443 only to represent the case where metal atoms, M, are coordinated with a large amount of carbon anions. The unit cell of each system, MB 2 , MBC and MC 2 shown in Fig. 1 contains two formula units (FUs) each. All third to sixth row elements from group 4 to 14 was included in the screening process along with Li (as a reference). A total of 36 elements were screened where for each element, M, and all cations in the target structures ( Fig. 1) were replaced, resulting in 108 unique structures.
After relaxation of each structure following the procedure given in ''Screening of HEC candidate elements''. Four selections criteria were used to determine possible HEC candidate elements.

Results and discussion
As described above, selection criterion 1 concerns the formation energy, E f , which is shown in Fig. 2a for all systems. From the figure, it is clear that this criterion is rather relaxed and is fulfilled by many metals, Li, Al, Ti, V, Cr, Mn, Fe, Zr, Nb, Mo, Tc, Hf, Ta and W, several of which are known to form single-phase high-entropy ceramics [2]. For these metals, the MB 2 system has lower formation energy compared to the MBC and MC 2 systems, except for one metal, Li, where the formation energy is the lowest for the MBC system. Criterion 2 concerns the charge transfer, q Ã , between the cation and anion sublattices which signifies the strength of the interlayer bonding (ionic binding) between the respective lattices. Only elements with large charge transfer ( À 0:5 e -) for the MB 2 , MBC and MC 2 systems pass, and as seen in Fig. 2b, this is a farliy strict criterion which only six metals meet, namely Li, Ti, V, Zr, Nb and Hf.
Atomic radius is a common measure of the size of an atom, although several methods exist for Figure 1 Crystal structures of the three systems used for screening of HEC candidate elements. Here, the large gray spheres represent cations, the small green spheres boron and the small brown spheres carbon atoms (anions). a the MB 2 system, b the MBC system and c the MC 2 system, M denotes the element screened, e.g., M = Li, Al, Si, Ti, V, etc. Note that each system contains two formula units. determining it; one way is to measure the distance from the center of the nucleus to the boundary of the surrounding electron density (charge density). This boundary is not well-defined since the charge density gradually decreases as one moves away from the nucleus. Furthermore, in condensed matter systems, the charge density of atoms usually overlap to some extent, making it difficult to define an atomic radius. Richard Bader suggested a method for defining atomic volumes in a condensed matter system based on the charge density [36]. The so-called Bader analysis gives a measurement of the size (volume) of an atom, and by assuming a spherical charge density, an atomic radius can be estimated.
Depending on the local environment of an atom (e.g., its bonding, charge transfer, etc.), the charge density surrounding it will change and thus the size of the atom will change. This effect is shown in Fig. 2c where the atomic radius, as calculated using Bader analysis, changes depending on the system. Comparing the atomic radius in the cation sublattice for the MB 2 and MBC systems show that for most elements, the atomic radius is similar between the two systems. Some exceptions here are Al, Zn and Ga  Figure 2 a Formation energy, E f , b charge transfer, q Ã , between the cation and anion sublattices and c atomic radius, R Ã , of the cations for all investigate elements and structures in total 108 systems. Blue dots represent the MB 2 system, red dots the MBC and yellow dots the MC 2 system. which show large difference in their atomic radius depending on the system (MB 2 or MBC). The opposite can be seen for the MC 2 system, which compared to the MB 2 and MBC systems shows large differences in atomic radius for most elements. Except for Li, Ti, V, Zr, Nb and Hf, which for all three systems have similar atomic radii, and these six metals are the only ones which meet criterion 3.
The final criterion (4) concerns the lattice parameters a, b and c and states that these should be similar for each system, MB 2 , MBC and MC 2 . Figure 3a shows the lattice parameters a, b for each element and system where a clear separation between each system can be seen. The MB 2 system has the largest lattice parameters (except for Zr and Hf) followed by the MBC and MC 2 systems. This seems counterintuitive since the atomic radius is the largest for the MC 2 system followed by MBC and MB 2 systems, as shown in Fig. 2c. However, this simply implies that the a, b lattice parameters are primarily determined by the bond lengths B-B (1.76 Å ), B-C (1.61 Å ) and C-C (1.47 Å ) in the anion sublattice. The values in parentheses are the average bond lengths in the anion sublattice for all elements screened. Hence, the lattice parameters a, b are not of interest when selecting candidate elements, instead focus is on lattice parameter c which is shown in Fig. 3b. By comparing the lattice parameters for the three systems MB 2 , MBC and MC 2 , four metals are found for which the lattice parameter, c, is similar in all three systems, and they are Li, Ti, Zr and Hf. A closer look reveals a stark resemblance between Fig. 3b and Fig. 2c where a strong correlation between atomic radius, R Ã , and lattice parameter, c, can be found, see Figure S1b in the supplementary information. The same strong correlation cannot be found for R Ã and the lattice parameters a, b; see Figure S1a in the supplementary information. Thus, one might consider a reduced set of criteria when selecting candidate elements, consisting of criteria 1 and 4 (only for lattice parameter c). This is easier from a computation standpoint since it only requires relaxation of the system. From the screening process described above, we identify four metals that passes the full set of selection criteria (1, 2, 3 and 4), Li, Ti, Zr and Hf, and two metals that pass three of the four selection criteria, V and Nb. The rest of the elements screened passesnone or only one of the selection criteria and are thus  of no further interest. With the six elements that meet at least three selection criteria (assuming an equal atomic ratio between the cations), it is possible to create one senary high-entropy borocarbide (Li 1/6 Ti 1/ 6 V 1/6 Zr 1/6 Nb 1/6 Hf 1/6 )BC, here labeled as 6-HEC abc (a=Li, b=V and c=Nb). Furthermore by using the four metals which pass the full set of selection criteria and adding either V or Nb, it is possible to create two quinary high-entropy borocarbides (Li 1/5 Ti 1/5 V 1/  Table 1 with an additional quaternary HEC labeled 4-HEC 0 which has a composition based on the experimentally synthesized high-entropy borocarbide by Zhang et al. [33][34][35]37].

Structural properties and thermodynamic stability
Using density functional theory (DFT) calculations with special quasirandom structures [38,39] (SQS), all HEC compositions in Table 1 were studied according to the procedure given in ''Construction of special quasirandom structures'', ''Relaxation of SQS and calculation of lattice distortions and material properties'' and ''Calculation of mechanical and thermal properties''. The structural properties for all HECs obtained after relaxation of both lattice parameters and atomic positions (second stage) using DFT are shown in Table 2. Here, lattice parameter a of the candidate HECs, 4-HEC a , 4-HEC b , 4-HEC c , 5-HEC ab , 5-HEC ac and 6-HEC abc , is slightly larger compared to that of the experimentally based composition, 4-HEC 0 . This indicates that lattice parameter a is predominantly determined by the covalent bonding in the anion sublattice. For lattice parameter c, the opposite is found, where all candidate HECs have a significantly smaller lattice parameter compared to that of the 4-HEC 0 .
The thermodynamic stability of a material can be estimated by calculating the formation energy, E f , and energy above hull, E h , respectively. Formation energy is calculated at 0 K and 0 atm. and can be considered as an approximation of the formation enthalpy at ambient conditions. As shown in Table 3, the formation energy for all candidate HECs is negative indicating that energetically favorable structures are formed. Energy above hull, E h , is a measurement of the thermodynamic stability of a material where E h ¼ 0 eV corresponds to the most stable material with respect to decomposition at the same, fixed chemical composition [40][41][42]. As expected, due to their high-entropy nature, all HECs have a positive energy above hull resulting in decomposition at high temperature into a mixture of borocarbides, borides, carbides and graphite as shown in Table 3. Interestingly, the HECs containing Hf ? Ta (4-HEC 0 ) and Hf ? Nb (4-HEC c , 5-HEC ac and 6-HEC abc ) are predicted to have decomposition products consisting of refractory bimetallic carbides HfTaC 2 and HfNbC 2, respectively.  Table 2 Structural properties for all HECs obtained after full relaxation of the SQS using DFT. Here, a, b and c are the lattice parameter per FU, V the volume per FU and q is the density The entropy forming ability (EFA) descriptor developed by Sarker et al. [9] is a measurement of the experimental synthesizability for a proposed HEC composition. EFA is defined as the inverse of the spread (standard deviation) in enthalpy (formation energy) per atom, including degeneracies, for all unique configurations of a target crystal structure with specified partial occupancies. The larger the spread in formation energy, the more energy is required to access the different configurations. Hence, the more difficult it is to introduce configurational disorder into the system, resulting in a lower EFA. High EFA values thus indicate a high probability of forming high-entropy single-phase, while low values indicate a low probability of forming high-entropy single-phase during synthesis. In their paper, Sarker et al. manage to achieve synthesis of single-phase HECs with EFA values as low as 50 eV -1 , where compositions with lower EFA values 45, 38 and 37 eV -1 formed multi-phase HECs.
In Table 3, the lowest EFA value is found for the 4-HEC 0 , with a value of 13.2 eV -1 , and thus it is predicted to have a very low probability of forming high-entropy single-phase during synthesis. Indeed, experimental synthesis using the same elements (Ti, Mo, Hf, Ta, B and C) shows the formation of multiphase HECs [34]. However, by adding SiC to the synthesis, the formation of a HEC composite was achieved [33] consisting of one SiC phase and one hexagonal high-entropy single-phase containing Ti, Mo, Hf, Ta, B and C. Note that although this highentropy single-phase has the same crystal structure as the 4-HEC 0 , theoretical results (see supplementary information to ref. [33]) suggests that it likely contains significantly more boron than carbon, an approximate ratio of 7 to 1.
The EFA for the candidate HECs listed in Table 3 predicts a high probability of forming high-entropy single-phase during synthesis for the 4-HEC c and 4-HEC a compositions, while the EFA for the 4-HEC b composition predicts a low probability. The EFA value for the 4-HEC c , 167 eV -1 , is higher than the highest previously reported value by Sarker et al. [9] 125 eV -1 for (Mo 1/5 Nb 1/5 Ta 1/5 V 1/5 W 1/5 )C. Although the EFA for the quinary and senary HECs was not determined, see ''Relaxation of SQS and calculation of lattice distortions and material properties'', the correlation between lattice distortions and EFA found by Sarker et al. can be used to roughly estimate their synthesizability. From the lattice distortions in Table 4, the EFA for the 5-HEC ac is estimated to be similar to that of the 4-HEC a , and thus it is expected to have a high probability of forming high-entropy single-phase during synthesis. The 5-HEC ab and 6-HEC abc have very similar lattice distortions which lie between that of the 4-HEC a and 4-HEC b , and thus the 5-HEC ab and 6-HEC abc are expected to have a medium probability of forming high-entropy singlephase during synthesis.

Lattice distortions
The fully relaxed SQS (second stage) show a variety of lattice distortions depending on the HEC composition as shown in Fig. 4. Here, the 4-HEC 0 , Fig. 4a, shows severe lattice distortions in both the cation and anion sublattices, while the candidate HECs, Fig. 4bg, show much smaller lattice distortions. Interestingly, the 4-HEC b , 5-HEC ab and 6-HEC abc , Fig. 4c, e, g, all of which contain V, show larger lattice distortions compared to the rest of the candidate HECs. For Table 3 Formation energy, E f , energy above hull, E h , entropy forming ability (EFA) and the most stable decomposition products for all HECs.
Here, E f , E h and the decomposition products are calculated using the fully relaxed SQS and given in eV per FU, for details see ''Relaxation of SQS and calculation of lattice distortions and material properties''. The EFA is calculated following the method described by Sarker et al. [9] and is given as the inverse of the energy per atom,  Table 4 is smaller for the candidate HECs compared to the 4-HEC 0 , in Table 4 Lattice distortions quantified for each HECs. Here, E r is the relaxation energy per FU and DV is the change in volume per FU during relaxation. d is the ALD for all atoms, d M the ALD for the cations and d B , d C the ALD for the B and C atoms, respectively. kÁd lÁDV is a measurement of the contribution to the relaxation energy for the ALD, d, and the volume change, DV  (2) is the energy gained by the HEC during relaxation (second stage as described in ''Relaxation of SQS and calculation of lattice distortions and material properties'') and is a combination of the change in volume, DV, Eq. (3) and the ALD, d, Eq. (7). If a linear relation is assumed, then the relaxation energy can be described as Fitting Eq. (1) against the data in Table 4 gives k ¼ À0:873 eV/Å and l ¼ 0:0846 eV/Å 3 , and by comparing the contribution to the relaxation energy for each term, see kÁd lÁDV in Table 4, it is clear that the relaxation energy is dominated by the lattice distortions. Plotting the ALD against the corrected relaxation energy, E r ¼ E r À l Á DV, as shown in Fig. 5a, shows a good agreement between the ALD and the relaxation energy. This relation can be described by d ¼ 1 k E r which can be approximated as d % 1 k E r for small changes in volume during relaxation. Here, a larger (more negative) relaxation energy corresponds to larger lattice distortions; this is as a result of the atoms moving from a non-equilibrium position to an equilibrium position lowering the total energy of the system.
To further study the lattice distortions, we calculate the directional lattice distortions (DLDs) for each atom in the HECs. As given in ''Relaxation of SQS and calculation of lattice distortions and material properties'', two DLDs are defined, one parallel, d k , and one perpendicular, d ? , to the c-axis, see Eq. (4) and (5). For the DLDs of each element in the HECs, a probability density function (PDF), P d x¼ k;? f g À Á , is estimated using the kernel density estimation method available in the seaborn [44] Python library (seaborn.kdeplot) which provides an estimation of the continuous distribution of the lattice distortions for each element based on the calculated DLDs for that element.
The PDFs for the 4-HEC 0 shown in Fig. 5b have a wide probability distribution (spread) with multiple overlapping peaks. Here, the spread in Pðd k Þ is larger compared to that of Pðd ? Þ, which shows that the lattice distortions are more severe along the c-axis compared to in the ab-plane. Looking at the Pðd k Þ for each element in Fig. 5b shows that the spread in Pðd k Þ is the largest for Mo compared to the other elements which indicates that the lattice distortions are most severe for Mo. An inspection of the relaxed SQS in Fig. 4a shows the same even if the large lattice distortions for the 4-HEC 0 makes this less visual apparent.
The DLDs for the candidate HECs shown in Fig. 5c-h have a much narrower spread compared to those for the 4-HEC 0 , and interestingly, the spread in Pðd k Þ is similar to that for Pðd ? Þ which proves that for the candidate HECs, the lattice distortions are small and directionally uniform, in agreement with Fig. 4. However, some exceptions to this can be found where the Pðd k Þ for V in all candidate HECs and for Ti in the 4-HEC c has a larger spread compared to Pðd ? Þ. This shows that for these two cations, especially V, the lattice distortions are larger compared to the rest of the cations in the candidate HECs and has a preferential distortion direction along the c-axis.

Mechanical and thermal properties
The elastic tensor for all HEC compositions in Table 1 was calculated using DFT with SQS, see ''Calculation of mechanical and thermal properties'' for details. Using the elastic tensors, available in the supplementary information, the mechanical and thermal properties can be predicted. Here, the bulk, shear and Young's modulus are predicted using the upper and lower bounds for polycrystalline materials (see supplementary information for details) are given in Table 5, where the highest stiffness, as measured by the Young's modulus, is found for the 4-HEC a and for the rest of the HECs the stiffness is significantly lower 5-HEC ac J 5-HEC ab J 6-HEC abc [ 4-HEC c J 4-HEC b [ 4-HEC 0 . Pugh's ratio, k, and Poisson's ratio, t, can indicate whether a material is ductile or brittle, the common values listed in the literature [45,46] for the transition between brittle to ductile are k  Table 5, all HECs except 4-HEC a are predicted to be ductile. The universal elastic anisotropy [47], A U , is a measurement of the directional dependency of the mechanical properties, and as shown in Table 5, A U is large for all HECs except the 4-HEC a which suggest that for most HECs, the mechanical properties are highly directionally dependent. Additional anisotropic indices A 1 , A 2 and A 3 for hexagonal systems [48], see Table S3 in the supporting information, further supports the high directional dependency of the mechanical properties. To further study this, the directional Young's modulus was calculated following the method described in ''Calculation of mechanical and thermal properties''. Figure 6 shows the Young's modulus in the ð001Þ, ð010Þ and ð100Þ planes for all HECs, from which it is evident that the stiffness along the crystal axes, a, b and c for all candidate HECs are similar and around 100 GPa higher in each axis compared to the 4-HEC 0 . All quinary and senary HECs along with the 4-HEC a show a very uniform Young's modulus in the ð001Þ plane but for the rest of the quaternary HECs the stiffness is less uniform in this plane.
Where the candidate HECs differ significantly from each other is in the off-axis stiffness, especially in the ð010Þ and ð100Þ planes. For all quaternary HECs, except the 4-HEC a , the stiffness in the ð010Þ and ð100Þ planes reach a minimum at 45°(135°e quivalent) from the c-axis. Because of their layered structure, a low stiffness when load is applied at 45-degrees from the c-axis indicates slippage between the cation and anion sublattices in the HECs. This would also correspond to a low shear modulus, which is indeed the case for the 4-HEC 0 , 4-HEC b and 4-HEC c but not for the 4-HEC a . The same phenomenon can be seen in the ð100Þ plane for the quinary and senary HECs. Note that the higher stiffness at 45°in the ð010Þ plane as compared to the same angle in the ð100Þ plane for the quinary HECs is likely an artifact of the finite size SQS used.
The Grü neisen parameter [49], c, is a measure of the anharmonic effects in a solid, such as thermal expansion and temperature dependence of phonon frequencies as well as linewidths. As shown in Table 5, the 4-HEC a has a relatively small value of c which indicates low anharmonic effects, while the 5-HEC ab , 5-HEC ac and 6-HEC abc have larger but mutually similar values. The remaining quaternary 4-HEC b and 4-HEC c as well as the 4-HEC 0 have higher anharmonic effects as shown by the Grü neisen parameter.
Theoretical prediction of the hardness of materials is a challenging problem since hardness is affected by many factors which requires extensive calculations to account for. Still, there are many semi-empirical models fitted to experimental data and theoretically calculated elastic properties, which can be used to calculate hardness [50]. But since the accuracy of these when applied to high-entropy materials is largely unexplored, we have opted to investigate the hardness using the Debye temperature. The Debye temperature is the temperature of the highest possible normal mode vibration in a crystal and can be related to the hardness of a material [51][52][53][54]. From Table 5, it is clear that the 4-HEC a has the highest Debye temperature followed by the 5-HEC ab J 5- Table 5 Mechanical and thermal properties for all HECs calculated using DFT. Here, B, G and Y denote the bulk, shear and Young's modulus measured in GPa. k is the Pugh's ratio, m the Poisson's ratio, A U the universal elastic anisotropy and c the Grüneisen parameter. H is the Debye temperature and the minimal thermal conductivity is estimated using the two methods of Clarke, j Clarke min , and Cahill, j Cahill min . Y=j min is the Braun figure of merit with j min ¼ j Clarke Y=j min (GPa m K W -1 ) HEC ac J 6-HEC abc [ 4-HEC c J 4-HEC b [ 4-HEC 0 . Thus, the 4-HEC a is predicted to have the highest hardness followed by the quinary and senary HECs, while the rest of the quaternaries (4-HEC c , 4-HEC b and 4-HEC 0 ) are predicted to be softer. The minimum thermal conductivity of a material, in the high-temperature limit, can be estimated using two different methods derived by Clarke [55] and Cahill [56], see supplementary information for details, and are listed in Table 5 for all HECs. Here, the Cahill method predicts a higher minimum thermal conductivity compared to the Clarke method, and the 4-HEC 0 is found to have the lowest minimum thermal conductivity using both methods. The minimum thermal conductivity for the candidate HECs is within 0.7 W m -1 K -1 of each other, where the lowest minimum thermal conductivity is found for the 4-HEC b and the highest for the 4-HEC a . For materials used as thermal barriers, it is desirable to have very low thermal conductivity, which generally comes at the expense of material stiffness. High-entropy materials on the other hand can achieve both high stiffness and low thermal conductivity due to large phonon scattering caused by disorder in the interatomic force constants [57]. Braun et al. used the ratio of the Young's modulus to the thermal conductivity as a figure of merit to measure the thermal barrier performance of entropystabilized oxides [29]. Table 5 shows the Braun figure of merit for all HECs which here is calculated as Y=j min with j min ¼ j Clarke min þ j Cahill min À Á =2. Note that j min is an estimate of the thermal conductivity in the hightemperature limit, and since stiffness is reduced at high temperatures [58], Y=j min can be seen as an upper limit for the figure of merit at high temperatures. From

Electron localization function
The interactions between the atoms inside of a material, whether they are covalent, metallic, ionic, etc., can be characterized by the electron localization function [61][62][63] (ELF) which is a relative measurement of electron localization and ranges from 0 to 1.
High ELF values ( [ 0:7) correspond to localized electrons in core or boning regions (covalent bonds) or lone pairs. ELF values between 0.7 and 0.2 correspond to an electron localization like that of the electron-gas and indicate metallic bonds. While ELF values larger than 0.2 is consistent with chemical bonding (shared-electron interactions), values lower than 0.2 are seen for physical binding (unsharedelectron interactions) such as ionic binding or dispersion (van der Waals binding).
Analysis of the ELF for the HECs in the (001) plane of their anion sublattice (left most panels in Fig. 7) shows high ELF values ( [ 0:8) for all HECs with attractor basins located between the B and C atoms. This proves the existence of covalent bonding between the atoms in the anion sublattice. Figure 7 clearly shows that the size of the attractor basins varies within the same 2D slice which indicates that the strength of the covalent B-C bonds varies throughout the material for all HEC compositions. Analysis of the ELF in the (001) plane of the cation sublattice (center panels in Fig. 7) shows ELF values ranging from around 0.2 to 0.5 with large attractor basins extending throughout the sublattice. This is indicative of metallic bonding and can be seen for all elements and all HEC compositions, with one exception. For the 4-HEC a , 5-HEC ab , 5-HEC ac and 6-HEC abc (Fig. 7b, e-g) all of which contain Li, the ELF around the Li atoms has a very low value (\0:05), see the dotted hexagon in Fig. 7e. Such low values correspond to highly delocalized electrons which is indicative of physical binding.
Analysis of the ELF in the (110) plane, which slices through both the cation and anion sublattice, can provide information on the interlayer bonding of the HECs (right most panels in Fig. 7). Here, physical binding, with a minimum between the cations and anions is present (dotted triangle in Fig. 7c), which is predominantly ionic in nature due to the charge transfer, see Table S1 in the supplementary information. In addition, the ELF has large attractor basins with values ranging from around 0.3 to 0.6 which connects the anion sublattices together, see the dotted rectangle in Fig. 7d. This implies that the HECs have metallic-like interlayer bonding between the anion sublattices in conjunction with physical binding (ionic binding) which occurs between the cation and the anion sublattices. The ELF in the (110) plane also shows the effect of Li on the interlayer bonding, where Li completely or partially disturbs the metallic-like interlayer bonding, depending on its neighboring cations. This effect is particularly evident for the 4-HEC a , see the dotted rectangle in Fig. 7b, where Li causes a discontinuity, i.e., ELF values lower than 0.2 in the surrounding attractor basins that connects the anion sublattices together.
Thus, the presence of Li in the HECs causes large local disturbances in the otherwise uniform ELF representative of metallic bonding. This is a result of Li creating a region around it with significantly higher delocalized electrons compared to those surrounding other metal atoms. During sliding of the sublattices in HEC, this region of delocalized electrons deforms the metallic-like interlayer bonding between the anion sublattices. This deformation is larger around the Li atoms compared to the other cations, which results in a significant increase in total energy during sliding for the HECs containing Li. Sliding is therefore more energetically unfavorable for the HECs containing Li leading to an increase in their shear modulus. We thus argue that Li acts as an anchoring mechanism that reduces the slippage between the cation and anion sublattices. This will improve the mechanical properties of the HEC, and as shown in Table 5, the shear modulus is indeed higher for all HECs containing Li (4-HEC a , 5-HEC ab , 5-HEC ac and 6-HEC abc ) compared to those without (4-HEC 0 , 4-HEC b and 4-HEC c ). More extensive work on Li's effect on the shear modulus in MBC systems will be presented elsewhere.
Our detailed study of lattice distortions shows that for the 4-HEC a , 4-HEC c and 5-HEC ac , the directional lattice distortions (DLDs), a concept introduced in this work, are uniform with respect to element and direction except for Ti in the 4-HEC c , for which the DLDs are larger parallel to the c-axis. This is however offset by the fact that the DLDs perpendicular to the c-axis are much smaller compared to the other HECs, resulting in the 4-HEC c having the smallest average lattice distortion (ALD) of all suggested HEC compositions. For the remaining 4-HEC b , 5-HEC ab and 6-HEC abc , the ALDs are larger, and from analysis of their DLDs, we find that they dominate parallel to the c-axis and that V contributes most to the increase in the ALDs.
Calculation of mechanical and thermal properties shows that the 4-HEC a has a Young's modulus of 426 GPa, which is significantly higher than the rest of the proposed HEC compositions. The calculated Debye temperature predicts that the 4-HEC a has the highest hardness of the suggested compositions. Furthermore, the low thermal conductivity of the 4-HEC a (2.47 W m -1 K -1 ) in combination with its high mechanical properties results in a Braun figure of merit of 179 GPa m K W -1 , making the 4-HEC a an excellent candidate for use as thermal barrier coating.
Analysis of the electron localization function (ELF) shows that the atoms in the anion sublattice are covalently bound, while the atoms in the cation sublattice are metallically bound, and between the cation and anion sublattices, physical binding occurs which is predominantly ionic in nature. Interestingly, we find a previously unreported metallic-like interlayer bonding which connects the anion sublattices together. Furthermore, the inclusion of Li in the HECs causes large local changes in the otherwise uniform ELF representative of metallic bonding. These local changes increase the sliding resistance between the cation and anion sublattices resulting in an increased shear modulus.
Based on the results of our theoretical study, we recommend that future work focus on the synthesis and characterization of the 4-HEC a , (Li 1/4 Ti 1/4 Zr 1/ 4 Hf 1/4 )BC, due to its desirable mechanical and thermal properties. The 4-HEC c , (Ti 1/4 Zr 1/4 Nb 1/4 Hf 1/ 4 )BC and the 5-HEC ac (Li 1/5 Ti 1/5 Zr 1/5 Nb 1/5 Hf 1/5 )BC are also of interest because of their high probability of forming single-phase during synthesis and their low lattice distortions.
The plane wave energy cutoff, 520 eV, k-point density length parameter, R k ¼ 30 Å , and choice of pseudopotential followed that of the Materials Project [70][71][72] in order for the calculations to be comparable with those available in their database. For all calculations, the electronic self-consistent loop was converged to 10 À5 eV using no symmetry constraints (ISYM = 0) and high precision (PREC = Accurate). To reduce the noise in forces, an additional support grid for the evaluation of the augmentation charges was used (ADDGRID = .TRUE), and since the calculations include 3d elements, aspherical contributions to the gradient corrections were also included (LASPH = .TRUE).
The Methfessel-Paxton scheme [73] (ISMEAR = 1) for partial occupancies was used and the smearing width was set to 0.05 eV.

Screening of HEC candidate elements
As given in the ''Introduction'' section, screening was performed using three systems representing different compositions in the anion sublattice (Fig. 1). For each structure, the lattice parameters were relaxed using VASP (ISIF = 6) until the change in total energy was smaller than 10 -4 eV. The same settings as described above was used with the addition of spin polarization (ISPIN = 2) since some systems are magnetic. After relaxation, the formation energy, =FU, for each system was calculated by subtracting from its total energy, E tot , the energy per atom, E i , of each element's most stable unary phase available through the Materials Project [72] (database version 2019.05) using the Python Materials Genomics [74,75] (pymatgen) Python library. Bader analysis [76] was performed to determine the atomic radius, R Ã , and charge transfer, q Ã , for the cations. To generate the required input for Bader analysis, single-point calculations (NSW = 0) were performed on the optimized systems with LCHARG = .TRUE and LAECHG = .TRUE.

Construction of special quasirandom structures
The substitutional disorder present in high-entropy systems such as HECs makes DFT modeling of these systems challenging. When using periodic boundary conditions, as is the case for VASP, special periodic structures are needed whose correlation functions for the first nearest neighbor shells are as close as possible to those of a disordered system. Such periodic structures called special quasirandom structures [38,39] (SQS) are supercells built from the unit cell of the target HEC. For sufficiently large SQS, periodicity errors only exist between distant neighbors, and thus the local substitutional disorder is captured. However, such ideal SQS are only obtainable for equiatomic systems with few components such as binary and possible ternary random alloys. For HECs, ideal SQS are not obtainable at supercell sizes suitable for DFT calculations (hundreds of atoms). Instead, several SQS must be generated for each HEC composition where a good SQS (for a particular composition) can be selected based on criteria such as mechanical stability. Another more rigorous approach for determining a good SQS is to check the variance in a chosen property (total energy, bandgap, mechanical properties, etc.) for all possible permutations of the same SQS (same composition). However, when the number of components and the size of the SQS increase, this method becomes very computationally demanding.
Here, the Monte Carlo based method (mcsqs) for generating SQS available through the Alloy Theoretic Automated Toolkit [77] was used with the hexagonal LiBC unit cell (P6 3 /mmc) shown in Fig. 1b as the input structure. With increasing substitutional disorder, much larger SQS are needed [78,79] and thus to reduce size of the as generated SQS, substitutional disorder was only considered in the cation sublattice. Specifically, the alternating B-C structure for anion sublattice shown in Fig. 1b was kept for all SQS. It would of course be beneficial to include substitutional disorder in the anion sublattice since it would increase the entropy of mixing, but this is left for future studies. The SQS size used varies based on the number of components of the HEC and was chosen as follows: a 4 9 4 9 2 supercell (64 FU, 192 atoms) for quaternary, a 5 9 5 9 2 supercell (100 FU, 300 atoms) for quinary and a 4 9 4 9 3 supercell (96 FU, 288 atoms) for senary HECs. For each HEC composition, five SQS were generated using mcsqs considering only pair correlation functions within a cutoff of 5.5 Å (mcsqs -2 = 5.5). All five SQS for each HEC composition were then fully relaxed and the mechanical properties were calculated using DFT as given in ''Calculation of mechanical and thermal properties''. The mechanically unstable SQS, as determined by eigenvalues of the elastic tensor, were discarded. And of the remaining ones, one was chosen as representative based on the mechanical properties (highest Young's modulus).

Relaxation of SQS and calculation of lattice distortions and material properties
The generated SQS were all relaxed with VASP using the settings described in the beginning of this section; here, spin polarization was not used (ISPIN = 1) since all HEC compositions are non-magnetic. Relaxation was performed using the conjugate gradient method (IBRION = 2) in a two-stage manner. In the first stage, only the lattice parameters were relaxed until the maximum stress tensor component [80,81] was smaller than 10 À2 eV/Å 3 resulting in a relaxed structure with internal forces but close to zero external stress and zero lattice distortions (no distortion of the lattice sites). In the second stage, the SQS were further relaxed allowing for relaxation of both lattice parameters and atomic positions resulting in a structure with close to zero external stress and internal forces but with lattice distortions. Using the structures and their energy obtained from the two relaxation steps, the relaxation energy, E r , and relaxation volume, DV, were calculated as where E, V is the energy, volume per FU of the fully relaxed SQS (second stage) and E 0 , V 0 is the energy, volume per FU of the SQS allowing for relaxation of lattice parameters only (first stage). From the two-stage relaxation method described above, the concept of directional lattice distortions (DLDs) is proposed. These are the displacements from the ideal lattice sites u 0 , v 0 and w 0 (first stage) to the final lattice sites u, v and w (second stage) calculated as In Eq. (4), a, b and c are the lattice vectors of the fully relaxed SQS (second stage). Due to the layer structure of the HECs, the DLDs are defined as parallel, d k , and perpendicular, d ? , to the c-axis.
The total lattice distortion for each atom is then defined as Finally, the lattice distortion for the whole SQS is defined as the average lattice distortion (ALD) for all atoms where N is the number of atoms in the SQS and d i is the total lattice distortion of the i th atom as defined by Eq. (6). The ALD for the boron, d B , and carbon, d C ; atoms or the cations, d M , can be calculated in the same manner as d.
Using the fully relaxed SQS (second stage), the electron localization function [82] (ELF) was calculated, and Bader analysis was performed to determine the atomic radius, R, and charge transfer, q, for each cation in the HEC. Furthermore, the formation energy, was calculated using the total energy, E tot , of the fully relaxed SQS (second stage) in combination with the energy per atom, E i , of each element's most stable unary phase available through the Materials Project (database version 2019.05). Similarly, the energy above hull, E h , was calculated considering all the possible phases containing any combination of Li, Ti, V, Zr, Nb, Mo, Hf, Ta, B or C available through the Materials Project (database version 2019.05). Calculations were performed using the pymatgen Python library for materials analysis which in addition to E f and E h also calculates the expected decomposition products based on the composition of the HEC.
For each HEC composition, the entropy forming ability (EFA) descriptor developed by Sarker et al. was determined following the method described in ref. [9] with one exception. The enthalpies, H i , for the unique supercells generated by the automatic flow partial occupation (AFLOW-POCC) algorithm [83] was here calculated using the Materials Project database instead of the AFLOW database [84]. Note that the definition of enthalpy, H i , used by Sarker et al. is the same as the definition of formation energy, E f , used here. From the AFLOW input files (PARTCAR files available in the supplementary information) created based on the LiBC unit cell (Fig. 1b), the AFLOW-POCC algorithm generates 335 configurations (24 atoms each) for the quaternary HECs, 4923 configurations (30 atoms each) for the quinary HECs and 280,721 configurations (36 atoms each) for the senary HEC. Due to the large number of configurations generated, the EFA was only calculated for the quaternary HECs.

Calculation of mechanical and thermal properties
Mechanical properties of a system can be derived from the 6 9 6 elasticity tensor, C, and the compliance tensor, S ¼ C À1 . Here, calculation of the elastic tensor follows the method described by Jong et al. [85], where the lattice parameters of the fully relaxed SQS (second stage) are deformed using six unique deformation modes where each deformation is applied individually to the fully relaxed SQS. For each deformation mode, four different deformation magnitudes were used, which resulted in 24 deformed structures per SQS. For each deformed structure, VASP was used to relax the atomic positions while keeping the lattice parameters fixed and to calculate the stress tensor. The elastic tensor is then obtained by a least-squares-fit of the applied strain and calculated stress tensors for all 24 deformed structures. Here, the pymatgen Python library was used to generate the deformed structures and fit the elastic tensor using the applied strain and the calculated stress tensors.
From the elasticity tensor, mechanical properties were determined including the bulk modulus, B, shear modulus, G, and Young's modulus, Y, along with the Pugh's ratio, k, Poisson's ratio, t, Grü neisen parameter [49], c, universal elastic anisotropy [47], A U and the additional anisotropic indices A 1 , A 2 and A 3 for hexagonal systems [48]. The directional dependency of the material properties was investigated by calculating the directional Young's modulus using the elastic tensor in combination with the ELATE [86] Python library. Note that for the calculation of the directional Young's modulus, ELATE requires a symmetric elasticity tensor, C Ã , which was estimated as C Ã ¼ C þ C T À Á =2. Thermal properties were also derived from the elasticity tensor including the Debye temperature [87], H and minimal thermal conductivity [55,56]

Data availability
Data that support this work is available from Dr. Daniel Hedman (daniel.hedman@ltu.se) upon reasonable request.

Code availability
DFT calculations were performed with the Vienna Ab initio Simulation Package (www.vasp.at). All other code that was developed for this work is available from Dr. Daniel Hedman (daniel.hedman@ltu.se) upon reasonable request.

Declarations
Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licen ses/by/4.0/.