|  Prospect of high-temperature superconductivity in layered metal
borocarbides†
---|---
|  Charlsey R. Tomassetti,a‡ Gyanu P. Kafle,a‡ Edan T. Marcial,a Elena R.
Margine,a and Aleksey N. Kolmogorov∗a
|  Delithiation of the known layered LiBC compound was predicted to induce
conventional superconductivity at liquid nitrogen temperatures but extensive
experimental work over the past two decades has detected no signs of the
expected superconducting transition. Using a combination of first-principles
stability analysis and anisotropic Migdal-Eliashberg formalism, we have
investigated possible LixBC morphologies and established what particular
transformations of the planar honeycomb BC layers are detrimental to the
material’s superconductivity. We propose that LixBC reintercalation with
select alkali and alkaline earth metals could lead to synthesis of otherwise
inaccessible metastable LixMyBC superconductors with critical temperatures
($T_{\rm c}$) up to 73 K. The large-scale exploration of metal borocarbides
has revealed that NaBC and Li1/2NayBC layered phases are likely true ground
states at low temperatures. The findings indicate that this compositional
space may host overlooked synthesizable compounds with potential to break the
$T_{\rm c}$ record for conventional superconductivity in ambient-pressure
materials.

##

††footnotetext: a Department of Physics, Applied Physics, and Astronomy,
Binghamton University-SUNY, Binghamton, New York 13902, USA††footnotetext: ∗
E-mail<EMAIL_ADDRESS>† Electronic Supplementary
Information (ESI) available: Additional figures and tables pertaining to
stability, structural, electronic, vibrational and superconducting properties
for relevant materials, x-ray diffraction comparisons with experiment, and
computational details. See DOI: 00.0000/00000000.††footnotetext: ‡ These
authors contributed equally to this work.

##  1 Introduction

Observations of superconducting transitions in ambient-pressure MgB2 1 and
compressed hydrides 2, 3, 4 have shown that high-critical temperature ($T_{\rm
c}$) conventional superconductivity is attainable in compounds with markedly
different chemistry and morphology. The shared feature underpinning their
exceptional superconducting properties is the presence of doped covalent
bonds. The partially filled electronic states in these materials generate a
sizable density of states (DOS) at the Fermi level while the hard vibrational
modes of strongly bonded light elements are able to withstand pronounced
softening from large electron-phonon (e-ph) coupling. These traits have been
widely used to guide the development of new superconductors but,
unfortunately, hole-doping of covalent bonds is detrimental to thermodynamic
stability and the issue has been circumvented successfully in only a handful
of cases. In synthesized H3S and LaH10 superhydrides with the respective
$T_{\rm c}$ of 203 K 2 and 250 K 3, formation of 3D metallized covalent
frameworks at ‘unnatural’ stoichiometries is forced via application of extreme
pressures. In the 55 K Q-carbon-based superconductor, the heavy B doping (27
at%) of the predominantly $sp^{3}$ amorphous C network is accomplished via
quenching of laser-heated samples 5. The quasi-2D MgB2 with a $T_{\rm c}$ of
39 K 1 possesses naturally hole-doped B layers with a significant amount of
key B-$\sigma$ states at the Fermi level because of the step-wise dependence
of DOS on the number of carriers in 2D systems. The proposed closely related
LiB with a similarly large electronic DOS and softened in-plane B phonon modes
6 forms just above 20 GPa and remains metastable at ambient conditions 7 but
the prediction of its superconductivity up to 34 K 8 has yet to be tested
experimentally.

A prime ambient-pressure candidate is a delithiated form of the bulk LiBC
material comprised of honeycomb BC layers. It was expected to supplant the
record-holding MgB2, with predicted $T_{\rm c}$ as high as 100 K 9. However,
after two decades, this material has never been observed to superconduct
despite extensive experimental effort 10, 11, 12, 13, 14, 15. The exact
reasons for why LixBC ($1>x$) has failed to live up to its full potential are
still not clear. Fogg et al. 14 sought to characterize the delithiation
process and observed expulsion of B with the eventual breakdown of the
quasi-2D BC network for $x$ values below $\sim 0.45$. Their density functional
theory (DFT) simulations revealed the increasing prevalence of a specific
defect in the BC lattice upon the reduction of Li content. A different end-
product configuration was later proposed by Kalkan and Ozdas 15. They inferred
that, much like in graphite intercalation compounds (GICs), Li extraction from
LiBC leads to the formation of Daumas-Hérold (DH)-type domains 16.

As the first step toward modeling the delithiation process, our recent ab
initio study focused on establishing the $(T,P_{\textrm{Li}})$ thermodynamic
conditions necessary to destabilize the starting LiBC material 17. The
constructed phase diagram proved to be in excellent agreement with the
synthesis conditions in successful delithiation experiments 11, 14.
Importantly, we considered all previously reported Li-B-C phases and
demonstrated that representative LixBC derivatives with $1>x\geq 0.5$ are only
metastable under ambient temperature and pressure, which indicates that the
covalent BC layers are strong enough to kinetically protect the delithiation
pathway 17.

Fig. 1:  Compositions of observed (grey points) and some of the proposed (red
points) compounds in the Li-M-B-C chemical space. The red (blue) triangles
highlight the location of possible holed-doped layered LixMyBC superconductors
partially reintercalated with alkali (alkaline earth) metals.

In this work, we address two questions pertaining to the superconducting
potential of hole-doped layered borocarbides. Firstly, we seek to attain a
better understanding of the LixBC morphologies produced during experiments and
use this knowledge to explain why no sample has ever been observed to
superconduct. Our systematic screening of LixBC configurations uncovers a
number of alternative motifs that become favored over the signature honeycomb
BC framework upon the Li extraction. The following examination of the
identified candidates’ properties within the anisotropic Migdal-Eliashberg
(aME) formalism reveals that ordered LixBC configurations should be
outstanding superconductors with $T_{\rm c}$ over 30 K but certain defects are
indeed detrimental to the e-ph coupling. Secondly, we explore the possibility
of obtaining synthesizable superconductors in the Li-M-BC composition subspace
shown in Fig. 1. We demonstrate that reintercalation of the metastable LixBC
with different metals (M = Na, Mg, K, or Ca) may be possible thermodynamically
provided that the material can withstand the process in its kinetically
constrained layered form. Moreover, the DFT results indicate that phases near
the Li1/2Na1/2BC composition are tantalizingly close to the convex hull and
could form from the elements. Our aME calculations show that ordered LixMyBC
configurations would have comparably high critical temperatures.

The expected accessibility of these (meta)stable layered phases via
reintercalation, ion exchange, electrochemical 18, 19, 20, 21, or direct
reactions makes the proposed materials class conceptually different from the
large set of recently considered superconductors with shared chemistry in
terms of synthesis conditions (e.g., at high pressures) and/or morphology
(e.g., in 2D or clathrate forms): Li4B5C2 22, Li3B4C2 23, 2D LiB2C2 24,
Li2xBC3 25, MxBC2 26, MgB3C3 27, 2D Mg2B4C2 28, X-Y-B2 29, Li-B-C 30, and
X-Y-B6C6 31. While the feasibility of obtaining our identified quaternary
phases depends on a number of kinetic and thermodynamic factors that require
further study, the findings dramatically expand the search space for
high-$T_{\rm c}$ ambient-pressure superconductors.

Fig. 2:  Schematic potential energy profile for Li1/2BC configurations. The
most favorable oI10 structure with ordered honeycomb BC layers is used as a
reference and its select transformations are illustrated with representative
structures. The heights of barriers separating local/global minima are not to
scale. The approaches used to explore the LixBC configurational space are
listed at the bottom. Our superconductivity calculations indicate that Li1/2BC
polymorphs with short- or long-range ordering of honeycomb BC layers and Li
vacancies should be high-$T_{\rm c}$ superconductors.

##  2 Methods

Stability analysis of Li-M-B-C phases was conducted with VASP 32 using
projector augmented wave potentials 33 and a 500 eV plane-wave cutoff. In
order to account for the dispersive interactions important in layered
materials 34, 35, 7, 36, we relied on the optB86b-vdW functional 37 but also
used Perdew-Burke-Ernzerhof (PBE) parametrization 38, optB88-vdW 39, and
r2SCAN+rVV10 40, 36 functionals to check the sensitivity of the results to the
DFT approximations. All structures were evaluated with dense ($\Delta k\sim
2\pi\times 0.025$ Å-1) Monkhorst–Pack ${\bf k}$-meshes 41.

Global structure optimizations were carried with an evolutionary algorithm
implemented in the MAISE package 42. In fixed-composition runs, randomly
initialized 16-member populations with up to 21 atoms per unit cell were
evolved for up to 200 generations using standard mutation and crossover
operations 42. For the systematic screening of possible decorations of
interlayer sites in metal borocarbides we started with stoichiometric MBC
supercells with up to 72 atoms and sequentially substituted or/and removed
metals keeping only unique configurations. Detection of similar structures for
elimination of duplicates in evolutionary optimizations and combinatorial
explorations was done with our fingerprint based on radial distribution
functions (RDFs) 43, 42. The thermodynamic corrections due to vibrational
entropy were evaluated within the finite displacement methods implemented in
PHONOPY 44. We used supercells between 56 and 256 atoms, with the average of
152 atoms over 80 phases, and applied 0.1 Å displacements within the harmonic
approximation. According to our previous quasi-harmonic approximation results
for related Li-B-C and Na-Sn materials 17, 45 and present tests for NaBC and
Li1/2Na1/2BC (Figs. S15-S16), volume expansion has an insignificant effect on
the formation free energies in the considered temperature range.

Detailed information about structure, stability, and $T_{\rm c}$ for key Li-M-
B-C phases discussed in the manuscript can be found in Table 1. To avoid
ambiguity when referring to different phases at the same composition, we use
Pearson symbols and space groups when needed. In addition, we indicate
morphology of considered LixBC variants with superscripts using abbreviations
explained in Fig. 2 and below. Full structural information for relevant DFT-
optimized Li-M-B-C phases is provided as CIF files.

For calculating properties related to superconductivity, we employed the
Quantum ESPRESSO package 46 with the optB86b-vdW functional 37 and norm-
conserving pseudopotentials from the Pseudo Dojo library 47 generated with the
relativistic PBE parametrization 38. A plane-wave cutoff value of 100 Ry, a
Methfessel-Paxton smearing 48 value of 0.02 Ry, and $\Gamma$-centered
Monkhorst-Pack 41 k-meshes were used to describe the electronic structure. The
lattice parameters and atomic positions were relaxed until the total energy
was converged within $10^{-6}$ Ry and the maximum force on each atom was less
than $10^{-4}$ Ry/Å. The dynamical matrices and the linear variation of the
self-consistent potential were calculated within density-functional
perturbation theory 49 on irreducible sets of regular q-meshes. The optimized
lattice parameters for the investigated phases and the k\- and q-meshes used
are reported in Tables S1 and S4 50.

The EPW code 51, 52, 53, 54 was used to investigate the e-ph interactions and
superconducting properties. The Wannier interpolation 55, 56, 57 was performed
on uniform $\Gamma$-centered k-grids (see Table S4 50) using the Wannier90
code 55, 56 in library mode. We used three $p$ orbitals for every C atom as
projections for the maximally localized Wannier functions to accurately
describe the electronic structure of all the compounds under investigation
(see Fig. S22 50). The anisotropic Migdal-Eliashberg equations were solved on
fine uniform ${\bf k}$\- and ${\bf q}$-point grids (see Table S4 50), with an
energy window of $\pm 0.2$ eV around the Fermi level and a Gaussian smearing
value of 50 meV. The Matsubara frequency cutoff was set to 1.5 eV and we
employed a uniform sampling scheme when solving the isotropic or anisotropic
Migdal-Eliashberg equations. The sensitivity of $T_{\rm c}$ estimates to the
$\mu^{*}$ parameter was assessed by performing calculations not only with the
standard 0.10 value but also with 0.20.

##  3 Results and Discussion

Table 1:  Thermodynamic stability calculated with VASP 32 and PHONOPY 44 and superconducting critical temperature calculated with Quantum ESPRESSO 46 and EPW 52, 54 in select LixBC and LixMyBC phases. The full set of observed phases used in the construction of the convex hulls is given in Fig. S12 50. The phases defining the thermodynamic stability of LixBC and LixMyBC phases are given as a superscript number at each temperature, and the corresponding numbered reference phases are listed as follows: (1) hP6-LiBC, (2) oF32-LiB6C, (3) mS30-B4C, (4) hP4-C, (5) hP6-NaBC, (6) oP14-NaB5C, (7) cF4-Na, (8) cF4-K, (9) oF72-KC8, (10) oS80-MgB2C2, (11) oI60-MgB12C2, (12) tI20-CaB2C2, (13) cP7-CaB6, and (14) cF4-Ca. The critical temperatures were obtained by solving the isotropic Migdal-Eliashberg (iME) and anisotropic Migdal-Eliashberg (aME) with the EPW code for $\mu^{*}=0.10$ and $\mu^{*}=0.20$. Composition | Space group - | Distance to convex hull [meV/atom] |  $T_{\rm c}$ [K]
---|---|---|---
| Pearson symbol | 0 K | 0 K | 300 K | 600 K | 900 K | iME | iME | aME | aME
|  |  | (+ZPE) |  |  |  | 0.10 | 0.20 | 0.10 | 0.20
Li1/3BC3D |  $Amm2$-oS14 |  155.71,2,4 |  |  |  |  |  |  |  |
Li1/2BC3D |  $C/m$-mS20 |  118.21,2,4 |  111.41,3,4 |  111.01,3,4 |  109.31,3,4 |  107.01,3,4 |  |  |  |
Li1/2BCSRO |  $Pnnm$-oP10 |  210.61,2,4 |  195.11,3,4 |  191.51,3,4 |  182.81,3,4 |  172.61,3,4 | 17.4 | 9.8 | 42 | 32
Li1/2BCSRO |  $Imm2$-oI10 |  209.91,2,4 |  195.01,3,4 |  191.81,3,4 |  184.11,3,4 |  174.91,3,4 | 33.7 | 23.0 | 49 | 38
Li1/2BCRH |  $P\bar{3}m1$-hP5 |  244.61,2,4 |  229.51,3,4 |  227.51,3,4 |  221.71,3,4 |  214.41,3,4 | 49.8 | 36.7 | 84 | 71
Li1/2BCSRO |  $Pmma$-oP10 |  220.91,2,4 |  |  |  |  | 33.9 | 23.1 | 82 | 70
Li1/2BCSRO |  $P\bar{3}m1$-hP15 |  225.51,2,4 |  209.11,3,4 |  205.51,3,4 |  196.71,3,4 |  186.21,3,4 | 53.3 | 40.2 | 72 | 60
Li1/2BCSRO |  $P312$-hP15 |  219.51,2,4 |  |  |  |  |  |  |  |
Li1/2BCIB |  $P2/m$-mP20 |  169.91,2,4 |  159.61,3,4 |  157.41,3,4 |  152.01,3,4 |  145.61,3,4 | 0.5 |  |  |
Li1/2BCIB |  $Pmm2$-oP20 |  186.91,2,4 |  |  |  |  | 1.6 |  |  |
Li1/2BCBR |  $P/m$-mP15 |  205.81,2,4 |  197.41,3,4 |  194.11,3,4 |  186.81,3,4 |  178.91,3,4 | 0.0 |  |  |
Li1/2BCBR |  $P/m$-mP30 |  231.41,2,4 |  |  |  |  |  |  |  |
Li1/2BCBR |  $P1$-aP15 |  208.11,2,4 |  |  |  |  |  |  |  |
Li5/8BCSRO |  $C2/m$-mS42 |  148.31,2,4 |  |  |  |  | 40.5 | 28.4 | 59 | 48
Li2/3BCSRO |  $P6_{3}/mcm$-hP16 |  118.01,2,4 |  |  |  |  | 47.0 | 34.0 | 83 | 71
Li3/4BCSRO |  $Pmma$-oP22 |  85.91,2,4 |  78.91,3,4 |  77.81,3,4 |  74.71,3,4 |  71.01,3,4 | 42.6 | 29.6 | 82 | 72
Li5/6BCSRO |  $P\bar{3}1m$-hP17 |  55.61,2,4 |  |  |  |  | 28.6 | 17.3 | 59 | 49
Li1/2Na1/2BC |  $P\bar{3}m1$-hP6 |  3.51,5 |  3.51,5 |  3.61,5 |  11.41,4,6,7 |  22.21,4,6,7 |  |  |  |
Li1/2Na5/12BC |  $P2/m$-mP35 |  6.91,4,5,6 |  3.31,4,5,6 |  2.21,4,5,6 |  6.31,4,6,7 |  12.41,4,6,7 |  |  |  |
Li1/2Na3/8BC |  $P2/m$-mP23 |  9.61,4,5,6 |  5.01,4,5,6 |  3.51,4,5,6 |  6.01,4,6,7 |  10.21,4,6,7 | 28.5 | 17.1 | 75 | 67
Li1/2Na1/3BC |  $P\bar{3}m1$-hP17 |  20.81,4,5,6 |  16.01,4,5,6 |  15.21,4,5,6 |  17.91,4,6,7 |  21.91,4,6,7 | 35.5 | 23.1 | 77 | 63
Li1/2Na1/3BC |  $P2/c$-mP34 |  15.41,4,5,6 |  10.81,4,5,6 |  9.01,4,5,6 |  10.51,4,6,7 |  13.31,4,6,7 |  |  |  |
Li1/2Na1/4BC |  $P2/m$-mP11 |  56.71,4,5,6 |  |  |  |  | 32.9 | 21.6 | 74 | 62
Li1/2K1/4BC |  $P2/c$-mP22 |  35.61,2,8,9 |  |  |  |  |  |  |  |
Li1/2K1/6BC |  $P\bar{3}1m$-hP16 |  79.51,2,8,9 |  |  |  |  | 43.7 | 31.5 | 84 | 73
Li1/2Mg1/6BC |  $P\bar{3}1m$-hP16 |  61.41,4,10,11 |  |  |  |  | 33.7 | 21.7 | 62 | 49
Li1/2Ca1/6BC |  $P\bar{3}m1$-hP16 |  57.71,4,12,13 |  |  |  |  | 34.8 | 23.1 | 58 | 45
Li1/2Ca1/3BC |  $P\bar{3}m1$-hP17 |  74.81,12,14 |  |  |  |  | 2.4 |  |  |
Li1/2Ca1/2BC |  $P\bar{3}m1$-hP6 |  120.01,12,14 |  |  |  |  | 1.8 |  |  |

###  3.1 Background of LixBC

After Rosner et al. theorized that hole-doped lithium borocarbides could
exhibit superconductivity up to 100 K 9, several groups have synthesized
delithiated LiBC 10, 12, 13, 11, 14, 15 and concluded from resistivity
measurements that the material shows no sign of the expected transition down
to 2 K. The studies found that it is difficult to control the exact Li content
removed due to the dynamic and nonequilibrium nature of the delithiation
reaction. The typical route for obtaining LixBC involved high temperature Li
deintercalation 12, 11, 14, 15. At the higher 1500∘C end of the investigated
temperature range, experiments produced two groups of phases, one with Li
contents $x>0.8$ and the other with $0.6>x$ 12, 14, 15. The use of low-
pressure setups allowed Zhao et al. 11 to induce deintercalation at a lower
600-800∘C temperature and produce three hole-doped samples of nominal
compositions Li0.80BC, Li0.77BC, and Li0.63BC. The use of the oxidizing agent
NOBF4 in acetonitrile made it possible to extract Li at near-ambient 95∘C but
Raman spectra indicated that samples with $x=0.5$ and 0.65 consisted of a
disordered graphite-like BC material alongside unreacted LiBC 13.

The study by Fogg et al. 14 took a particular focus on the delithiation limit
of LixBC. After observing the expulsion of B-rich byproducts from their Li-
poor ($0.45>x$) samples, they conjectured that a type of defect consisting of
a B-C atom swap could happen at higher compositions ($x>0.45$). The increasing
occurrence of this defect with decreasing Li content was suggested to
kinetically ease the degeneration of Li0.45>xBC into non-layered C-rich and
B-rich phases. Their DFT calculations agreed that such a transformation would
become energetically favorable at Li content $x=0.5$, with recent work 17
corroborating that a low-symmetry structure with swapped B and C atoms becomes
slightly favored over the standard morphology with perfectly alternating B-C
sites. Bond rotation is a commonly observed transformation in covalent
honeycomb systems, such as the dislocation dynamics of the 5-7 Stone-Wales
defect pairs in graphene probed with low-voltage high-resolution transmission
electron microscopy 58 or the appearance of B-B bonds in LixBC, indicative of
B-C swaps, detected with spectroscopic methods in the discharged LiBC anode
material tested recently for Li-ion battery applications 59.

The ability of borocarbides to form graphitic honeycomb networks raises the
question of whether they exhibit properties similar to that of GICs 60 or
transition metal dichalcogenides (TMDs) 61, including stage formation with a
DH-type domain intercalation mechanism 16. Staging is the phenomenon in which
intercalants are placed between $n$ layers of unintercalated host material,
with said material classified as stage-$n$. At low stages, the domain size is
expected to be large 60, which does not preclude the possibility of mixed-
stage samples and staging disorder. Based on the observation of a 12∘ (001)
powder x-ray diffraction (PXRD) peak in Li0.56BC, Kalkan and Ozdas 15 surmised
the existence of stage-2 ordering, known to occur in GICs 60 and TMDs 61 but
not yet detected in LixBC, as other studies reported reflections only down to
15-20∘ 62, 13, 11. Kalkan and Ozdas further proposed a mixed-stage phase of
24% (wt %) stage-1 LiBC and 76% stage-2 Li0.43BC as the best fitting model to
the PXRD data. Moreover, their Reitveld analysis offers evidence of shifting
and buckling of the BC layers, similar to that observed in Mg1/2BC 63.

###  3.2 Stability of LixBC

The investigation of the LixBC subspace in our previous study of the full Li-
B-C ternary focused on finding (i) the lowest-energy partially delithiated
configurations with ordered honeycomb BC layers in the $3/4\geq x\geq 1/4$
range by screening LiBC supercells with up to 18 atoms and (ii) alternative
stable morphologies with $1/2\geq x$ by performing unconstrained evolutionary
searches 17. The present thermodynamic stability analysis has been extended to
identify a wider variety of metastable LixBC structures relevant in the
context of kinetics-defined delithiation experiments and examine their
favorability across the composition range. The exploration of this
configuration space involved a suite of complementary strategies detailed
below and in the Methods Section.

First, we carried out a broader set of global structure searches using the
evolutionary algorithm implemented in the MAISE package 42. We (re)examined
key stoichiometries around $x=1/2$, namely, 2/3, 5/8, 1/2, 3/8, 1/3, etc., to
establish when non-honeycomb motifs become favored in the deintercalation
process. Previously, we observed that Li depletion below $x=1/3$ triggers the
BC layers to fuse into fully-connected 3D frameworks producing, e.g.,
oS14-Li1/3BC3D 120 meV/atom below the LiBC$\leftrightarrow$BC tie-line 17. The
present evolutionary optimizations of 16-member populations with up to 21
atoms per unit cell for up to 200 generations, totaling 30,000 local
optimizations across all compositions, helped us refine these findings. We
have now detected emergence of more stable 3D BC frameworks at a higher
$x=1/2$ Li concentration. For instance, the new mS20-Li1/2BC3D phase with a
porous BC scaffold hosting double rows of Li in 1D channels (Fig. 2), was
found to be nearly 90 meV/atom lower in energy (Fig. 3) than any half-filled
layered variant considered previously 14, 17. Inspection of local minima
produced in the evolutionary searches revealed that bridging between two
layers through vacant Li sites is favorable enough to occur at $x>1/2$. The
large unit cell sizes needed to represent such phases with nontrivial
compositions, e.g., 5:8:8 or 7:8:8, made it difficult to ensure that the
candidates found ‘from scratch’ were true ground states.

To pinpoint the crossover composition between the ordered and altered layered
morphologies, we proceeded with a combinatorial screening of intercalant
decorations in known and identified BC frameworks. The procedure involved
generating different LiBC (super)cells and sequentially removing the metal
atoms while keeping only non-equivalent configurations. Duplicate structures
were identified using our structural fingerprint based on the radial
distribution function 43, 42. We focused on the $1\geq x\geq 1/2$ range in
LixBC structures with three competing BC patterns: the original honeycomb
layers displaying short-range order (SRO), interlayer bridging (IB) between C
atoms, and defective layers comprising B-C bond rotations (BR).

Starting with perfectly ordered honeycomb BC layers, we created various
supercells to examine both AA’ and AA stacking variations corresponding to the
BC-CB and BB-CC placements of the B and C atoms along the $c$ axis,
respectively. Overall, we generated 12 supercells with up to 54 atoms and
ultimately evaluated over 2,500 distinct structures. At the $x=1/2$
composition, the more extensive screening reproduced the previously identified
oI10-Li1/2BCSRO phase with buckled layers as the most stable one with this
morphology. At two other previously sampled $x>1/2$ stoichiometries, the new
best hP16-Li2/3BC and oP22-Li3/4BC represented only minor improvements in
stability, by less than 1 meV/atom. The sampling of additional $x>1/2$
fractions revealed no particularly favorable ordered configurations, with the
cluster of red points in Fig. 3 reaching a modest minimum of -13.6 meV/atom
below the LiBC$\leftrightarrow$Li1/2BC tie-line at $x=2/3$.

Interestingly, for some Li vacancy arrangements at the lowest considered
filling fraction of $x=1/2$, the interlayer bridging occurred naturally upon
simple local relaxation of the starting planar geometries. To investigate the
possibility of bridging at higher Li contents, we prepared supercells with up
to 48 atoms extended along one in-plane direction, removed two neighboring Li
atoms, brought two C atoms in adjacent layers half-way to within the covalent
distance across the vacant metal sites, and systematically probed possible
arrangements of remaining Li atoms. After local optimization of over 600
candidate phases, we established that the IB configurations indeed become
preferred over the SRO phases below $x=2/3$. The increasing predominance of
the motif, manifested in the relative energy drop from $-12$ meV/atom at
$x=2/3$ down to $-46$ meV/atom by $x=1/2$ in Fig. 3, precipitates the eventual
linkage of layers into fully connected 3D frameworks.

Fig. 3:  Relative energies of LixBC phases referenced to LiBC and the lowest-
energy oI10-Li1/2BCSRO. The notation of the examined morphologies, shown in
different colors, is explained in Figs. 2 and 4. The gray solid line
represents the boundary of the convex hull. The lowest-energy Li1/2BC
structures for each morphology type are hP5RH, oI10SRO, mP15BR, mP20IB, and
mS203D (see Table 1)

.

Fig. 4:  Stability and structural parameters of Li-intercalated materials. (a)
Relative energies of Li1/2BC (Li1/2C6) phases referenced to the lowest-energy
oI10-Li1/2BCSRO (Rüdorff-Hofmann stage-2 Li1/2C6) as a function of the domain
length in the Daumas-Hérold model in units of the $a$ lattice constant, with
the armchair (zigzag) sets corresponding to DH boundaries along (perpendicular
to) the covalent bonds. (b) The composition-dependent $c$ lattice constants in
LixBC models fully relaxed in our DFT calculations (black stars) and reported
in experimental studies (other symbols). (c) $c$ lattice constant in Li1/2BC
as a function of the domain length.

The sampling of a larger set of unit cells with rotated BC bonds improved on
the aP15-Li1/2BCBR phase found previously with the evolutionary search. The
new mP15BR configuration at this composition is slightly more stable by 2.3
and 4.1 meV/atom compared to the aP15BR and oI10SRO variants, respectively. As
was discussed previously 14, 17, the likelihood of this defect type is heavily
suppressed at compositions above $x=1/2$. Our screening of unit cells with
higher Li content revealed that the most stable structure with rotated bonds
at $x=2/3$ was disfavored by nearly 50 meV/atom (Fig. 3).

Finally, we performed a comparative analysis of the stability and structural
trends in the half-filled borocarbide (Li1/2BC) and graphite (Li1/2C6)
intercalation compounds. It should first be noted that stage-2 hP13-Li1/2C6 is
a thermodynamically stable phase in the Li-C binary, while stage-2 hP5-Li1/2BC
is nearly 35 meV/atom higher in energy than our most stable oI10-Li1/2BCSRO.
In the DH domain model, intercalant islands left over from the delithiation
process define boundaries between different staging domains. To simulate the
energetics of increasingly large islands, we adopted a strip model with edges
along either armchair or zigzag orientations of the honeycomb BC lattice. A
representative armchair unit cell of length $a=10.80$ Å is shown in Fig. 2.
The resulting relative energies in both compounds are plotted in Fig. 4(a) as
a function of the domain length. As the strip size increases, the energy
gradually converges towards a Rüdorff-Hofmann-type (RH) model 64, a perfectly
boundaryless stage-2 configuration. These findings illustrate a stark
difference between the energetics of domain formation in the Li-intercalated
graphite and borocarbide compounds and disagree with the proposed explanation
of the stage-2 configuration emerging in the latter around $x=1/2$ 15.

Further evidence against the DH/RH models in LixBC can be seen in the
material’s structural response to the Li island formation summarized in Fig.
4(b-c). The $c$ lattice constants extracted from XRD or neutron diffraction
measurements 10, 11, 12, 13, 14, 15 show little variation during delithiation
and remain between 7.0 Å and 7.1 Å down to $x\approx 0.6$. At lower Li
concentrations, the reported values in samples obtained via different
synthesis routes disperse over a significant range, from 7.14 Å at $x=0.43$ 15
to 7.44 Å at $x=0.56$ 12, 14 and up to 7.61 Å at $x\approx 0.4$ 14. Our
calculated $c$-axis lattice constants agree well with the experimental values
for the $1\geq x\gtrsim 0.6$ compositions, which indicates a good performance
of the chosen van der Waals (vdW) functional for the description of the
interlayer interactions in LixBC. At $x=1/2$, the DFT results demonstrate that
the $c$ lattice constant would remain around 7.1 Å if the products maintained
the SRO configurations but would decrease dramatically, down to nearly 6.4 Å,
if the compound indeed adopted the stage-2 morphology. The trend differs from
the well-established expansion of the interlayer distance in the Li-
intercalated graphite en route from LiC6 (measured 3.687 Å 65 and our
calculated 3.603 Å) to stage-2 Li1/2C6 (measured 7.0229 Å 65 and our
calculated 6.936 Å).

Therefore, staging does not seem to be a likely scenario behind the appearance
of the key (001) XRD peak in LixBC samples obtained in Kalkan and Ozdas’
experiments 15, and other possibilities may need to be considered to explain
the evident doubling of the period along the $c$-axis. We have simulated XRD
patterns for a variety of identified competing LixBC configurations and
plotted them in Fig. S1 50. It is clear that structures with broken hexagonal
symmetry caused by layer buckling, bond rotation, or interlayer bridging
cannot fit the observed signature peaks. The closest match is our
hP15-Li1/2BC${}^{\textrm{SRO}}$ ($P312$) phase, 9.6 meV/atom above
oI10${}^{\textrm{SRO}}$, which is a $\sqrt{3}\times\sqrt{3}$ supercell with
1/3 and 2/3 ordered Li populations in adjacent galleries (see Supporting
Information 50). In fact, all considered $\sqrt{3}\times\sqrt{3}$ supercells
for $5/6\geq x\geq 1/2$ produce well-matching reflection patterns except for
the weak (100) peak around 22∘ that is not present in the collected data (the
peak at the nearby 21∘ angle was attributed to the LiB3 delithiation byproduct
15).

Overall, our thermodynamic stability analysis indicates that LixBC prefers to
retain perfectly ordered honeycomb layers down to about $x=2/3$ and then shows
propensity to forming interlayer bridges that ultimately link BC layers into
extended 3D frameworks once the Li content reduces to $x=1/2$. The periodic
structures with rotated bonds become only marginally preferred over the SRO
variants at $x=1/2$. In fact, the results in Table 1 show that the inclusion
of zero point energy alone reverses the order in favor of the latter. The
vibrational entropy contributions bring the full set of the considered
morphology types at $x=1/2$ closer to the convex hull boundary by 24-39
meV/atom at $T=900$ K but neither make them globally stable and nor cause any
other significant changes in the relative ranking. The configurational
entropy, estimated in our previous study to reach 22 meV/atom at 900 K 17, may
promote disordered decorations of Li sites but delithiation at lower
temperatures and/or subsequent annealing are expected to counteract the
disorder of LixBC structures.

Given the variety of experimental conditions and the significant challenges of
simulating large-scale kinetic processes, such as staging or solid state
transformations involving covalent rebonding, it is difficult to establish
what particular metastable phases are produced in specific experiments.
Nevertheless, the library of possible LixBC phases gives us a chance to probe
the response of superconducting properties to a variety of factors.

###  3.3 Superconductivity of LixBC

The primary purpose of the following analysis is to establish how possible
compositional and morphological transformations taking place in the
delithiation process affect the LixBC derivatives’ superconducting properties.
It has been widely discussed that the starting LiBC compound is isoelectronic
and isostructural to MgB2 9, 66 but, being a semiconductor, it has no
signature hole-doped $\sigma$ states at the Fermi level essential for the
high-$T_{\rm c}$ superconductivity. In the delithiated form, the material does
develop multiple Fermi surface (FS) sheets of quasi-2D $p$-$\sigma$ or 3D
$p$-$\pi$ character shown in Fig. S2 50. So far, LixBC superconducting
properties have been investigated only in a handful studies with at best the
isotropic Eliashberg calculations 67, 66. Considering the importance of
anisotropic effects in the formation of multiple superconducting gaps that
define the critical temperature in related layered conventional
superconductors, e.g., MgB2 68, 69, 8, GICs 70, 71, and Li-Mg-B 8, we have
relied on the aME formalism to reexamine the e-ph coupling in LixBC as well.

The $T_{\rm c}$ dependence on the Li concentration in LixBC has been
previously considered only within the jellium model 66. The study found a
consistent boost in the isotropic e-ph coupling strength ($\lambda$) with the
increasing hole concentration, triggered by a rise in the DOS at the Fermi
level, $N(E_{\textrm{F}}$), and softening of the logarithmic average
frequency. In the present work, we have performed calculations for a series of
the most stable SRO configurations identified across the $5/6\geq x\geq 1/2$
compositions. Their representation with ordered unit cells allows one to probe
effects beyond the rigid band approximation arising not only from a non-
uniform charge distribution but also from variable shifts of different $s$ and
$p$ states on B/C atoms around Li vacancies, which proved to be important for
establishing the LixB optimal stoichiometry 72 and evaluating the LiB DOS
response to doping 73. Table S1 50 shows a slightly upward trend for the
important $N_{\sigma}(E_{\textrm{F}}$) with the number of holes but the total
DOS at the Fermi level exhibits a non-monotonous behavior. The analysis of the
phonon spectra (Fig. S3 50) reveals that the frequencies of the predominate
large-coupling bond-stretching modes at $\Gamma$ harden from about 80 meV
($x=5/6$) to about 90 meV ($x=1/2$) as $x$ decreases, which coincides with the
steady increase of the e-ph coupling strength from $0.7$ to 1.0 (see Fig. S3
50). However, the aME results indicate that superconducting properties are
sensitive to both composition and structure, as the $T_{\rm c}$ in SRO phases
with $5/6\geq x>1/2$ scatter between 48 K and 72 K (at $\mu^{*}=$0.20, see
Fig. S4 50) and even more so for various polymorphs at $x=1/2$, as shown in
Fig. 6(a,b). Note that in addition to the standard $\mu^{*}$ value of 0.10 74,
we used a value of 0.20 for all considered compounds to assess the sensitivity
of $T_{\rm c}$ on this semi-empirical parameter. The latter provides more
conservative estimates, lower by about 10 K (20%), for MgB2-type
superconductors 75, 53, 8 (see Table 1).

Fig. 5:  Properties of five Li1/2BC layered configurations, with the degree of
BC layer distortion increasing from (a) to (e). Each row contains the
electronic band structure and DOS with orbital characters (two left panels),
the phononic DOS, the Eliashberg spectral function $\alpha^{2}F(\omega)$ with
the integrated e-ph coupling strength $\lambda(\omega)$ (right panel), and
structural models with the buckling size.

To appreciate the impact of different geometric factors on the material’s
superconducting properties, we focused on the single $x=1/2$ composition near
the lowest experimentally achievable limit of Li content at which the compound
retains the stoichiometric honeycomb BC layers. The most natural deviation
from the ideal planar morphology, and apparently the only one supported by
measurements, is the corrugation of the covalent layers deduced from Raman and
PXRD characterizations of LiBC/LixBC 76, 15. It has been pointed out that the
BC layers are ‘floppy’ and, just like in Mg1/2BC 63, tend to pucker around
metal vacancies and/or stacking faults 76, 15. The observed behavior is
consistent with previous 15, 17 and current DFT results showing that
structures with uneven BC layers are thermodynamically favorable. In addition
to the most stable oP10SRO ($Pnnm$) and oI10SRO configurations with high
degree of buckling, we considered hP5RH, oP10SRO ($Pmma$), and hP15SRO
($P\bar{3}m1$) models with low out-of-plane distortions (see Fig. 5 and Table
S1 50). We defined the buckling magnitude as the deviation from the mean for
the B and C atom vertical positions in each layer and found it to range from
0.01 Å in hP15SRO ($P\bar{3}m1$) to 0.16 Å in oP10SRO ($Pnnm$).

The hexagonal-to-orthorhombic symmetry breaking has a pronounced effect on the
electronic states shown in Fig. 5. In all three orthorhombic models, the
important $\sigma$ bands split at the $Z$ point by up to 1 eV, and in oI10SRO
the upper $\sigma$ pair additionally splits at $\Gamma$ due to the asymmetric
distortion pattern that shortens half of the distances between B-B and C-C
atoms in adjacent layers (Fig. 5). While the resultant band energy shifts have
a moderate impact on the $N_{\sigma}(E_{\textrm{F}})$ contribution because of
the quasi-2D nature of the $\sigma$ states, the markedly different hole-doping
levels lead to substantial variations in the size of the signature hourglass-
shaped Fermi surfaces along the $\Gamma-Z$ line. In particular, the nearly
filled bottom pair of the bands in oP10SRO ($Pnnm$) generates Fermi surfaces
with dramatically reduced radii (Fig. S5 50).

Fig. 6:  Dependence of critical temperature in LixBC phases on the (a)
composition, (b) degree of BC layer distortion, and (c) BC layer
transformation. Diamond (circle) symbols denote values obtained with the
anisotropic (isotropic) ME formalism. Panel (c) shows that interlayer bridging
and bond rotation effectively suppress superconductivity in Li1/2BC.

The layer distortions also reshape the isotropic Eliashberg spectral function
around the frequencies of the bond-stretching B-C modes known to couple
strongly with the $\sigma$ states. In the structures with low buckling (Fig.
5(a-c)), the profile is dominated by a sharp peak around 90 meV that tapers
off at higher frequencies. The large BC layer corrugation causes the peak to
split into two broader lower-intensity peaks centered around 90 meV and
110–120 meV (Fig. 5(d,e) and Fig. S6 50). The reduction of the integrated e-ph
coupling strength in the two sets of considered structures arises from the
weaker coupling in both regions below and above the 90-meV peak:
0.75+0.40=1.15, 0.75+0.45=1.20, and 0.55+0.44=0.99 in hP15SRO ($P\bar{3}m1$),
oP10SRO ($Pmma$), and hP5RH with nearly flat layers versus 0.55+0.28=0.83 and
0.55+0.19=0.74 in the significantly buckled oI10SRO and oP10SRO ($Pnnm$)
configurations.

By solving the aME equations 53, we obtained the energy distribution of the
superconducting gap $\Delta_{\bf k}$ as a function of temperature to find the
$T_{\rm c}$ of each structure, shown in Fig. S7 50. We observed a substantial
near-linear decrease from $\sim$ 71 K to $\sim$ 32 K (at $\mu^{*}=$0.20) in
the $T_{\rm c}$ values going from the least to most buckled configurations.
The anisotropic domes in the structures with essentially flat BC honeycomb
layers are well separated, which is consistent with the superconductivity two-
gap profiles found in the MgB2 and LiB materials with planar B frameworks 8.
In contrast, the lower energy gap in the corrugated structures has a
considerably larger spread and almost reaches the region of the higher energy
gap peaks. The findings indicate that BC layer distortions are indeed
detrimental to superconductivity but not sufficient to suppress it entirely,
as summarized in Fig. 6(b).

One can anticipate the bond rotation and interlayer bridging transformations,
favored around $x=1/2$, to have more drastic impacts on the material’s
electronic and vibrational properties. We chose four representative lowest
energy structures to examine the two defect types: oP20-Li1/2BCIB with AA-
stacking, mP20-Li1/2BCIB with AA’-stacking, mP15-Li1/2BCBR, and
mP30-Li1/2BCBR. In the bridged configurations, one can still identify hole-
doped $\sigma$ states but their contribution to the DOS at the Fermi level is
reduced significantly (Fig. S8(a,b) and Table S1 50). The Eliashberg spectral
function shapes have certain similarities with those found for the most
distorted SRO structures at this composition, with the major peak around 90
meV either split or shifted up to about 120 meV, but the overall $\lambda$
values are reduced by more than a factor of two down to 0.2–0.3 (Fig. S9(a,b)
50). Given the small e-ph coupling strengths, we did not attempt aME-level
calculations of the $T_{\rm c}$ and found the iME values to be under 1 K.

In mP15-Li1/2BCBR, the bond rotation causes the $\sigma$ states to drop below
the Fermi level by almost 1 eV leaving only the $\pi$-type Fermi surfaces
(Fig. S8(c)). The observation is consistent with the insightful explanation of
the electronic DOS response to B-C swap offered by Fogg et al. 14. We further
investigated the bond rotation impact on the electronic structure by
constructing a larger mP30BR unit cell with diluted concentration of the
defect. The top of the $\sigma$ manifold remained at $-1$ eV but the structure
developed a $\sigma$ band above the Fermi level slightly dispersed between 0.5
eV and 1 eV (Fig. S8(d) 50). Notably, these structures with direct C-C bonds
have harder phonon modes reaching 170 meV, which may be used to detect the
defects. With no $\sigma$ Fermi surfaces present in either mP15BR or mP30BR,
the materials have even lower $\lambda\sim 0.2$, respectively, and negligible
$T_{\rm c}$ values (Fig. S9(c,d) 50). We conclude that the appearance of such
defects in delithiated structures can indeed explain the reported lack of
detectable superconductivity in these intriguing ternary compounds.

###  3.4 Stability of LixMyBC

Fig. 7:  Structural and stability properties of double-metal borocarbides. (a)
Dependence of the interlayer distance on the intercalant ion size 77 at the
3:1 Li:M ratio. (b) Energy difference between configurations with mixed and
segregated distributions of extra metals. (c,d) Relative energies of
Li1/2M5y/[2(1-y)]BC phases referenced to oI10-Li1/2BC and pure M. The stable
LiBC is marked with a solid circle. Metastable phases defining the local
convex hull in the kinetically-protected subspace of layered borocarbides are
shown as half-filled points connected with solid lines. The global convex
hulls are displayed with dashed lines.

The rationale for adding another metal into the metastable delithiated phases
is to help recover the ordered layered morphology of the starting LiBC
material by suppressing or healing possible defects in the BC honeycomb layers
while keeping the bonds hole-doped. Relying on the fact that the covalent
framework in LixBC withstands temperatures exceeding 1,500∘C, we constrained
the search for viable LixMyBC phases to the honeycomb layered structures that
may be kinetically accessible through reintercalation. We focused on the Na,
K, Mg, and Ca alkali/alkaline earth metals, closest to Li in terms of size and
valence, to investigate whether the insertion of the metals is a
thermodynamically downhill reaction. To efficiently explore the considerably
expanded configuration space with $x\geq y$ and $1\geq x+y\geq 1/2$, we first
generated all possible unique decorations of metal sites for 6 unit cells with
up to 48 atoms. Our examination of metal distribution stability revealed that
only Mg could co-exist with Li in the same galleries, with the larger Na, K,
and Ca showing a clear trend to segregate into M-only and Li-only intercalant
layers (Fig. 7(b)). Based on this observation, we were able to probe much
larger unit cells with up to 72 atoms for specific compositions, such as
Li1/2NayBC with $1/2>y$.

Figure 7(c,d) summarizes our DFT results on the stability of quaternary
compounds relative to the starting lowest-energy
oI10-Li1/2BC${}^{\textrm{SRO}}$ and pure M. Naturally, the majority of the
considered double-metal borocarbides have the most pronounced stabilization,
by $\sim 0.2$–0.3 eV/atom, at compositions that restore the 8-electron count
(1:1 for monovalent and 2:1 for divalent metals). Compounds with the largest K
are an exception, achieving a comparable relative energy of $-0.16$ eV/atom at
a lower total metal content, Li1/2K1/4BC. Importantly, all four quaternary
systems feature phases, denoted with half-filled symbols, that form a convex
hull with respect to the endpoints and include both electron-doped (e.g.,
Li1/2Mg1/2BC and Li1/2Ca1/2BC) and hole-doped (e.g., Li1/2Na1/4BC and
Li1/2Ca1/6BC) derivatives (similar profiles can be seen for other fixed Li
concentrations up to $x=3/4$ in Fig. S10 50).

One can expect the formation of these locally stable phases to depend strongly
on the morphology of the starting Li1/2BC. If each gallery in the ternary
compound is half-filled, as our DFT calculations indicate, then intercalation
of the extra metal would result in mixed-metal structures suboptimal for all
but Li-Mg-BC quaternaries. However, even these heterogeneous configurations
are well below the oI10-Li1/2BC-M tie-line (e.g., for Li1/2M1/6BC in Fig.
7(b,c,d)) meaning that there would be a thermodynamic force driving the larger
metals into the ternary matrix. The ensuing expansion of the interlayer
distances would likely increase the migration of the Li ions 78 and could lead
to an eventual metal segregation through the sample edges. If Li1/2BC does
have the starting stage-2 morphology, the larger metals would need to overcome
a higher barrier to enter the empty galleries. Fig. 7(a) shows that the
interlayer spacing correlates well with the ion size and would have to
increase considerably, e.g., nearly double from 2.8 Å in hP5-Li1/2BCRH to 5.3
Å in Li1/2K1/6BC, to accommodate the extra alkali or alkaline earth metals.

The potential advantage of such quaternaries over LixBC can be seen from the
following analysis of the defect stability in various compounds. To understand
the importance of the size and electron count factors, we first simulated the
aforementioned rotated bond and interlayer bridge defects in the ternary
analogs with half-filled metal sites. Fig. S11(c) 50 illustrates that both
defective structures become unstable by at least 100 meV/atom in the
electronically neutral M1/2BC compounds (with M = Mg or Ca) and get
significantly disfavored in the hole-doped M1/2BC (with M = K). We next
evaluated the relative stability of the defective structures with similar
contents but different types of the intercalant metals, e.g., Li2/3BC versus
Li1/2M1/6BC. The overall trend, particularly prominent for bond rotations, is
a destabilization of the non-honeycomb BC morphologies in the quaternary
compounds with the larger and more electron-rich metals. It should be kept in
mind that the imperfections appearing in the covalent layers during high-$T$
delithiation are likely protected by high kinetic barriers and might not
easily heal upon insertion of additional metals but, being less stable in the
quaternary compounds, they could be more responsive to annealing.

Fig. 8:  Distance to the convex hull as a function of temperature for proposed
Na and Li-Na borocarbides. Details on the hP6-Li1/2Na1/2BC,
mP35-Li1/2Na5/12BC, mP23-Li1/2Na3/8BC, and mP34-Li1/2Na1/3BC phases can be
found in Table 1. Different lines correspond to the lowest-free energy
combinations among all observed ground states (circled numbers) and proposed
phases (bracketed letters), with the full set of previously synthesized Li-Na-
B-C phases listed in Fig. S12 50. NaBC and Li1/2NayBC ($y=$ 1/2, 5/12, and
3/8) appear below the known convex hull at low temperatures, which is
illustrated with dashed green lines. In particular, NaBC is thermodynamically
stable up to 400 K, being below a facet of the convex hull defined by the
observed C, Na, and NaB5C phases, while hP6-Li1/2Na1/2BC is below the C, Na,
LiBC, and NaB5C mixture up to 300 K (the free energy of Na evaluated for the
solid state even above the melting point of 371 K represents an upper bound).
Although the three quaternary phases in panels (b-d) become metastable upon
inclusion of the predicted NaBC in the optB86b calculations, mP23-Li1/2Na3/8BC
is globally stable according to the optB88 results (Fig. S14 50). Fig. 9:
Distance to the global convex hull for LixMyBC phases with honeycomb
morphology at $T=0$ K.

Finally, we discuss the global stability of the identified quaternary phases
by constructing the full convex hulls based on DFT energies of the Li-B-C
materials listed in ref. 17 and all reported relevant compounds with the four
additional metals 79, 80, 81, 82. For phases (near) stable at $T=0$ K, we also
calculated free energies by taking into account the vibrational entropy and
constructed convex hulls at elevated temperatures.

In the Li-K-B-C system with the large mismatch of the alkali metal sizes, no
quaternary phases with honeycomb BC layers appeared within 35 meV/atom to the
boundary of the convex hull (see Figs. 7(c) and S12(b) 50). In the Li-Ca-B-C
chemical space, the double-metal layered borocarbides are much more
competitive because the large divalent Ca ions do not need to occupy every
vacant metal site to make the electron-deficient BC network whole.
Unfortunately, this trait has already been observed in the known Ca1/2BC
compound with a covalent framework comprised of 4-8 polygons 83. The ternary
phase makes the quaternary borocarbides metastable by at least 30 meV/atom
across the investigated composition range (Figs. 7(d) and S12(d)) and could be
the end product of the Li1/2BC intercalation with Ca if the barrier to
rebonding between the two related layered morphologies is not sufficiently
high. The Li-Mg-BC metal borocarbides are technically even closer to
stability, with some phases being only 8 meV/atom above the boundary of the
convex hull (Fig. 7(d) and S12(c)). Regrettably, the existence of the Mg1/2BC
phase with the honeycomb BC layers 63 means that the phase separation of Li-Mg
borocarbides into the corresponding ternaries would be achieved by a simple
diffusion of metal ions within galleries, which greatly reduces the chances of
obtaining stable ordered quaternary compounds. We are aware of only a few
attempts to synthesize mixed-metal borocarbides. Bharathi et al. 10 observed
the formation of an electron-doped Li0.8Mg0.2BC, while Renker et al. 84
obtained Li0.95Na0.05BC. Mori et al. 85 reported the successful synthesis of
Mg0.25Li0.40BC, although no superconductivity was observed in their samples
down to 1.8 K. Their intended doping of Mg1/2BC with B or Li (under high
pressure) resulted only in large amounts of neighboring phases MgB2 and LiBC.

The Li-Na-B-C compositional space, on the other hand, may host true ground
states that have not yet been observed. The structure search by Miao et al.
demonstrated that hP6 is the lowest-energy configuration at the NaBC
composition but its global stability in the Na-B-C ternary system was not
evaluated 82. The Materials Project 79, 86, 87, 88 presently places this phase
(mp-1238775) 38 meV/atom above the convex hull facet defined by C, Na, and B4C
but the estimate is unreliable for several reasons. First, the graphite entry
(mp-2516584) contains only the Perdew-Burke-Ernzerhof (PBE) 38 total energy
and a manual evaluation of the distance based on the reported values in this
approximation produces 27.9 meV/atom, close to 25.0 meV/atom in our PBE
calculations. Second, none of the presently used Materials Project functionals
accounts for the dispersive interactions particularly important in the
considered layered materials 34, 35, 7, 36. Finally, the database is missing
the known NaB5C phase experimentally determined to have a simple cP7 structure
81 but difficult to simulate because of a random distribution of C on the B
sites. Having sampled $\sim 200$ unique decorations of different unit cells
with up to 35 atoms, we found an oP14 model to have the lowest energy, 23-40
meV/atom below the randomized $P1$ models generated by Delacroix et al. 81 for
calculating NMR spectra. The phase proved to be important in the construction
of the Na-B-C convex hull, as the
$\frac{4}{15}$C+$\frac{4}{15}$Na+$\frac{7}{15}$NaB5C mixture is $-16$ meV/atom
below $\frac{1}{4}$C+$\frac{1}{3}$Na+$\frac{5}{12}$B4C in our optB86b
calculations at $T=0$ K. Yet, NaBC is below either reference set, by $-20$ and
$-36$ meV/atom, respectively. Considering the diversity of the bonding types
and structural motifs defining the convex hull, one can expect the results to
be sensitive to the DFT approximation. With the optB88 39 (r2SCAN+rVV10 40,
36) functional, we find that NaBC is located even deeper, at $-49$ and $-74$
meV/atom ($-59$ and $-89$ meV/atom), below the two respective facets. Fig.
8(a) shows that the vibrational entropy quickly destabilizes this ternary
phase, which may explain why it has not been produced so far.

The double-metal Li1/2Na1/2BC compound in the simple hP6 structure also
appears stable relative to C, Na, NaB5C, and LiBC at $T=0$ K but is less
favorable with respect to LiBC and NaBC at all temperatures (Figs. 8(b)-9). It
is evident that for this midpoint quasi-ternary phase the Li and Na size
mismatch induces an elastic energy penalty. We investigated the interplay
between the strain and the electron count by considering Na-deficient
compounds shown in Fig. 8(c-e). The nearly fully-filled Li1/2Na5/12BC and
Li1/2Na3/8BC (see Figs. 8 and 9) prefer to have missing Na rows along the
zigzag direction, while mP34-Li1/2Na1/3BC achieves the highest stability via
stretching the Na sublattice to skip one hexagon row along the armchair
direction. Fig. 8(c) illustrates a clear benefit of removing a small amount of
the larger metal, which leads to not only a slightly lower energy at $T=0$ K
in Li1/2Na5/12BC but also a lower free energy at elevated temperatures due to
the softening of Na phonon modes in the less-crowded hole-doped material (Fig.
S13 50). Therefore, there is likely a range of temperatures around a few
hundred K and $1/2>y$ where Li1/2NayBC become true ground states, as our
optB88 calculations indicate for $y=3/8$ (Fig. S14 50), and may form from the
elements.

It should be noted that in the important NaB5C compound the configuration
entropy due to the disordered occupation of the B sites by C can noticeably
lower the free energy at high temperatures. However, our comparison of the
adopted oP14 representation with a considerable deviation from the cubic shape
($b/a=0.936$) against a far more uniform but low-symmetry aP189 model 81
revealed a significant energy penalty of $\sim 30$ meV/atom for NaB5C to have
the randomized B5C covalent framework that apparently remains frozen at low
temperatures. In contrast, the proposed Na-deficient Li1/2NayBC phases may
also benefit from the configurational entropy due to Na disorder at high
temperatures but should be able to become ordered through ion migration upon
annealing. It is also important to acknowledge that identification of
materials thermodynamically stable with respect to known ones only indicates
that additional phases, not necessarily the predicted ones, should exist in
the considered chemical system. Indeed, no search algorithm can guarantee the
identification of the full set of true ground states, especially in vast
multielement configuration spaces.

###  3.5 Superconductivity of LixMyBC

In the final set of calculations, we examined whether the quaternary layered
borocarbides identified as viable (meta)stable LixMyBC materials retain the
desired electronic and vibrational features to be good superconductors. The
phases chosen for this analysis are the lowest-energy stage-2 hP5-Li1/2BCRH
derivatives with extra metals reintercalated into the expanded empty galleries
at different Li:M ratios (4:3, 3:2, 2:1, or 3:1 for group-I and 2:1 or 3:1 for
group-II metals). The structural and electronic information of the six hole-
doped ($n$ = $+1/3$, $+1/4$, $+1/6$, and $+1/8$ holes/u.c.) and two electron-
doped ($n$ = $-1/6$ and $-1/2$ electrons/u.c.) phases is given in Table S1 and
Figs. S17 and S19 50. The electron doping levels achieved in Li1/2CayBC for
$y=$ 1/3 and 1/2 leads to the population of non-coupling interstitial states
and generates 3D Fermi surfaces (see Fig. S17 50). Unsurprisingly, the
resulting $T_{\rm c}$ of $\sim 2$ K obtained in the iME calculations are much
lower than those predicted for the corresponding hole-doped LixBC candidates.

The band structure and DOS plots in Fig. S19 50 indicate that the hole-doped
quaternary phases retain the signature 2D $\sigma$ bands. Keeping the electron
count constant by considering hP17-Li1/2Na1/3BC, Li1/2Mg1/6BC, and
Li1/2Ca1/6BC compositions against the control case of Li5/6BC, we could probe
the response of the electronic structure to the presence of different metals.
The $\sigma$-state contributions to $N(E_{\textrm{F}})$ remain the same in the
Na and Ca quaternaries and increase slightly in the Mg counterpart, while the
$\pi$-state contributions show a larger variability of $\sim 30$% because of
the different interlayer expansions and metal ion contents. The scan across
the LixMyBC series with doping levels ranging from $+1/8$ to $+1/3$ holes/u.c.
reveals $N(E_{\textrm{F}})$ variation by $\sim 40$%, from 0.20 states/(eV
atom) in hP17-Li1/2Na1/3BC to 0.29 states/(eV atom) in Li1/2K1/6BC, with no
clear trend.

The similarities between the ternary and quaternary hole-doped borocarbides
extend to the calculated Eliashberg spectral functions shown in Fig. S20 50.
Since the Li phonon modes couple with electrons weakly based on our
$\lambda$-resolved analysis, the frequency drop of the heavier metal
vibrations has little effect on the total e-ph coupling strength. The dominant
peaks present in all materials correspond to the in-plane bond-stretching B-C
modes, with the highest intensity peak occurring in the same range of
frequencies as seen in the LixBC ternary. The aME results summarized in Fig.
S21 50 demonstrate that the quaternary phases preserve the two highly
anisotropic superconducting gaps. The predicted $T_{\rm c}$ range from 45 K in
Li1/2Ca1/6BC to 73 K in Li1/2K1/6BC (at $\mu^{*}=0.20$), comparable to the
highest values obtained for the ternary LixBC phases (see Fig. 10).

Compared to MgB2, the quaternary borocarbides have similar 0.23 states/(eV
atom) values of $N(E_{\rm F})$, as shown in Figs. S19, S20, and S23 50. The
overall higher e-ph coupling strength in the latter can be attributed to the
contribution of low-frequency BC modes between 40-50 meV. For instance, below
the largest peak in $\alpha^{2}F(\omega)$, the modes contribute an average 42%
to total $\lambda$ in Li1/2MyBC, while only about 24% in MgB2.

Fig. 10:  Superconducting properties of metal borocarbides calculated at the
aME level with $\mu^{*}=0.20$. (a) $T_{\rm c}$ for select LixBC (grey
diamonds) and Li1/2MyBC (other symbols) phases as a function of hole doping.
For comparison, the $T_{\rm c}$ value of 40 K found for MgB2 at this level of
theory is displayed with a dotted line. (b) Superconducting gap in
Li1/2Na3/8BC.

##  4 Summary

The presented first-principles reexamination of LixBC provides new insights
into the seemingly well-studied material once considered to be one of the most
promising high-temperature ambient-pressure conventional superconductors.
Despite the extensive synthesis and characterization work 10, 12, 13, 11, 14,
15, our analysis indicates that the precise structures of some LixBC
derivatives obtained via delithiation at high temperatures remain unresolved.
To identify possible forms of the metastable hole-doped phases, we employed a
combination of structure search and analysis strategies (Fig. 2). Our modeling
of domain configurations has indicated that the stage-2 model proposed to
explain PXRD data 15 is not only energetically disfavored but also
inconsistent with the expected interlayer spacing dependence on the Li content
(Fig. 4). The evolutionary searches have uncovered significantly more stable
non-planar BC morphologies at $x=1/2$ including configurations with fully-
connected 3D frameworks. The combinatorial screening of possible Li
compositions and decorations has helped establish that the alternative BC
motifs become favored sooner than previously suspected 14, 17, around $x=5/8$
(Fig. 3). Even though the reported PXRD patterns clearly indicate that LixBC
retains its hexagonal layered structure down to $x=1/2$, the rebonded patterns
may appear in the form of local defects. We have used aME calculations to
quantify the effect of structural transformations on the critical temperature.
Surprisingly, the lowest-energy ordered layered configurations have been found
to have high $T_{\rm c}$ across the full examined $5/6\geq x\geq 1/2$ range.
The natural layer buckling has a relatively moderate effect on the
superconducting properties, reducing the highest $T_{\rm c}$ of $\sim$ 70 K by
roughly a factor of two (Table 1). The interlayer bridging suppresses $T_{\rm
c}$ down to a few Kelvin, while the intralayer B-C bond rotation extinguishes
the phonon-mediated superconductivity completely (Fig. 6). The absence of
superconductivity in annealed samples 11 suggests that such defects cannot be
completely eliminated, a factor that diminishes the prospect of high-$T_{\rm
c}$ superconductivity in the delithiated ternary borocarbide.

A possible solution examined in this study is to repopulate vacant metal sites
with other metals via soft reactions constrained by the layered honeycomb
morphology. We have determined that the defects detrimental for
superconductivity have higher energies of formation if galleries are filled
with larger ions and, hence, the double-metal borocarbides might be more
amenable to annealing than the starting LixBC material. We have shown that the
intake of Na, K, Mg, or Ca is thermodynamically favorable along the kinetics-
restricted pathway and could lead to formation of various hole-doped phases,
such as Li1/2K1/6BC or Li1/2Ca1/6BC (Fig. 7), with estimated $T_{\rm c}$ up to
73 K (Fig. 10). Borocarbides based on Na appear to be the most promising. Our
DFT calculations indicate that Li1/2NayBC with $y$ just below 1/2 could be
true ground states and, hence, synthesizable from the elements (Fig. 8). In
contrast to LixBC, these phases are naturally hole-doped and could be long-
sought-after ambient-pressure materials with high-$T_{\rm c}$ conventional
superconductivity. We have also demonstrated an unexpected thermodynamic
stability of the ternary NaBC analog of LiBC at low temperatures. If the
semiconducting compound does form, it would be interesting to investigate
whether it could be desodiated and turned into a high-$T_{\rm c}$
superconductor. The findings indicate that layered metal borocarbides comprise
a materials class uniquely suited to host viable ambient-pressure conventional
superconductors rivaling the MgB2 archetype.

## Author Contributions

C. R. Tomassetti and G. P. Kafle contributed equally to this work.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

The authors acknowledge support from the National Science Foundation (NSF)
(Awards No. DMR-2320073 and DMR-2132586). This study used the Frontera
supercomputer at the Texas Advanced Computing Center through the Leadership
Resource Allocation (LRAC) award DMR22004. Frontera is made possible by NSF
award OAC-1818253 89. This work also used the Expanse system at the San Diego
Supercomputer Center through allocation TG-DMR180071 from the Advanced
Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS)
program 90, which is supported by NSF grants #2138259, #2138286, #2138307,
#2137603, and #2138296. The authors thank Igor Mazin for helpful discussions.

## Notes and references

  * Nagamatsu _et al._ 2001 J. Nagamatsu, N. Nakagawa, T. Muranaka, Y. Zenitani and J. Akimitsu, _Nature_ , 2001, 410, 63.
  * Drozdov _et al._ 2015 A. P. Drozdov, M. I. Eremets, I. A. Troyan, V. Ksenofontov and S. I. Shylin, _Nature_ , 2015, 525, 73–76.
  * Drozdov _et al._ 2019 A. P. Drozdov, P. P. Kong, V. S. Minkov, S. P. Besedin, M. A. Kuzovnikov, S. Mozaffari, L. Balicas, F. F. Balakirev, D. E. Graf, V. B. Prakapenka, E. Greenberg, D. A. Knyazev, M. Tkacz and M. I. Eremets, _Nature_ , 2019, 569, 528–531.
  * Bernstein _et al._ 2015 N. Bernstein, C. S. Hellberg, M. D. Johannes, I. I. Mazin and M. J. Mehl, _Phys. Rev. B_ , 2015, 91, 060511.
  * Bhaumik _et al._ 2017 A. Bhaumik, R. Sachan, S. Gupta and J. Narayan, _ACS Nano_ , 2017, 11, 11915–11922.
  * Kolmogorov and Curtarolo 2006 A. N. Kolmogorov and S. Curtarolo, _Phys. Rev. B_ , 2006, 73, 180501.
  * Kolmogorov _et al._ 2015 A. N. Kolmogorov, S. Hajinazar, C. Angyal, V. L. Kuznetsov and A. P. Jephcoat, _Phys. Rev. B_ , 2015, 92, 144110.
  * Kafle _et al._ 2022 G. P. Kafle, C. R. Tomassetti, I. I. Mazin, A. N. Kolmogorov and E. R. Margine, _Phys. Rev. Mater._ , 2022, 6, 084801.
  * Rosner _et al._ 2002 H. Rosner, A. Kitaigorodsky and W. E. Pickett, _Phys. Rev. Lett._ , 2002, 88, 127001.
  * Bharathi _et al._ 2002 A. Bharathi, S. J. Balaselvi, M. Premila, T. Sairam, G. Reddy, C. Sundar and Y. Hariharan, _Solid State Commun._ , 2002, 124, 423–428.
  * Zhao _et al._ 2003 L. Zhao, P. Klavins and K. Liu, _J. Appl. Phys._ , 2003, 93, 8653–8655.
  * Fogg _et al._ 2003 A. M. Fogg, J. B. Claridge, G. R. Darling and M. J. Rosseinsky, _Chem. Commun._ , 2003, 3, 1348–1349.
  * Fogg _et al._ 2003 A. M. Fogg, P. R. Chalker, J. B. Claridge, G. R. Darling and M. J. Rosseinsky, _Phys. Rev. B_ , 2003, 67, 245106.
  * Fogg _et al._ 2006 A. M. Fogg, J. Meldrum, G. R. Darling, J. B. Claridge and M. J. Rosseinsky, _J. Am. Chem. Soc._ , 2006, 128, 10043.
  * Kalkan and Ozdas 2019 B. Kalkan and E. Ozdas, _ACS Appl. Mater. Interfaces_ , 2019, 11, 4111–4122.
  * Daumas and Hérold 1969 N. Daumas and A. Hérold, _CR Acad. Sci. Ser. C_ , 1969, 268, 373–375.
  * Kharabadze _et al._ 2023 S. Kharabadze, M. Meyers, C. R. Tomassetti, E. R. Margine, I. I. Mazin and A. N. Kolmogorov, _Phys. Chem. Chem. Phys._ , 2023, 25, 7344–7353.
  * Isaev _et al._ 1997 Y. Isaev, N. Lenenko, L. Gumileva, A. Buyanovskaya, Y. Novikov and E. Stumpp, _Carbon_ , 1997, 35, 563–566.
  * Abramchuk _et al._ 2017 M. Abramchuk, C. Ozsoy-Keskinbora, J. W. Krizan, K. R. Metz, D. C. Bell and F. Tafti, _J. Am. Chem. Soc._ , 2017, 139, 15371–15376.
  * Uppuluri _et al._ 2018 R. Uppuluri, A. S. Gupta, A. S. Rosas and T. E. Mallouk, _Chem. Soc. Rev._ , 2018, 47, 2401–2430.
  * Bahrami _et al._ 2019 F. Bahrami, W. Lafargue-Dit-Hauret, O. I. Lebedev, R. Movshovich, H.-Y. Yang, D. Broido, X. Rocquefelte and F. Tafti, _Phys. Rev. Lett._ , 2019, 123, 237203.
  * Bazhirov _et al._ 2014 T. Bazhirov, Y. Sakai, S. Saito and M. L. Cohen, _Phys. Rev. B_ , 2014, 89, 045136.
  * Gao _et al._ 2015 M. Gao, Z.-Y. Lu and T. Xiang, _Phys. Rev. B_ , 2015, 91, 045132.
  * Gao _et al._ 2020 M. Gao, X.-W. Yan, Z.-Y. Lu and T. Xiang, _Phys. Rev. B_ , 2020, 101, 094501.
  * Quan and Pickett 2020 Y. Quan and W. E. Pickett, _Phys. Rev. B_ , 2020, 102, 144504.
  * Hayami and Tanaka 2020 W. Hayami and T. Tanaka, _AIP Advances_ , 2020, 10, 065213.
  * Pham and Nguyen 2023 T.-T. Pham and D.-L. Nguyen, _Phys. Rev. B_ , 2023, 107, 134502.
  * Singh _et al._ 2022 S. Singh, A. H. Romero, J. D. Mella, V. Eremeev, E. Muñoz, A. N. Alexandrova, K. M. Rabe, D. Vanderbilt and F. Muñoz, _npj Quantum Mater._ , 2022, 7, 37.
  * Wang _et al._ 2022 R. Wang, Y. Sun, F. Zhang, F. Zheng, Y. Fang, S. Wu, H. Dong, C.-Z. Wang, V. Antropov and K.-M. Ho, _Inorg. Chem._ , 2022, 61, 18154–18161.
  * Zheng _et al._ 2023 F. Zheng, Y. Sun, R. Wang, Y. Fang, F. Zhang, S. Wu, C.-Z. Wang, V. Antropov and K.-M. Ho, _Phys. Rev. B_ , 2023, 107, 014508.
  * Geng _et al._ 2023 N. Geng, K. P. Hilleke, L. Zhu, X. Wang, T. A. Strobel and E. Zurek, _J. Am. Chem. Soc._ , 2023, 145, 1696–1706.
  * Kresse and Furthmüller 1996 G. Kresse and J. Furthmüller, _Phys. Rev. B_ , 1996, 54, 11169.
  * Blöchl 1994 P. E. Blöchl, _Phys. Rev. B_ , 1994, 50, 17953.
  * Kolmogorov and Crespi 2005 A. N. Kolmogorov and V. H. Crespi, _Phys. Rev. B_ , 2005, 71, 235415.
  * Lebègue _et al._ 2010 S. Lebègue, J. Harl, T. Gould, J. G. Ángyán, G. Kresse and J. F. Dobson, _Phys. Rev. Lett._ , 2010, 105, 196401.
  * Ning _et al._ 2022 J. Ning, M. Kothakonda, J. W. Furness, A. D. Kaplan, S. Ehlert, J. G. Brandenburg, J. P. Perdew and J. Sun, _Phys. Rev. B_ , 2022, 106, 075422.
  * Klimeš _et al._ 2011 J. Klimeš, D. R. Bowler and A. Michaelides, _Phys. Rev. B_ , 2011, 83, 195131.
  * Perdew _et al._ 1996 J. P. Perdew, K. Burke and M. Ernzerhof, _Phys. Rev. Lett._ , 1996, 77, 3865.
  * Klimeš _et al._ 2010 J. Klimeš, D. R. Bowler and A. Michaelides, _J. Phys. Condens. Matter_ , 2010, 22, 022201.
  * Furness _et al._ 2020 J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew and J. Sun, _J. Phys. Chem. Lett._ , 2020, 11, 8208–8215.
  * Monkhorst and Pack 1976 H. J. Monkhorst and J. D. Pack, _Phys. Rev. B_ , 1976, 13, 5188.
  * Hajinazar _et al._ 2021 S. Hajinazar, A. Thorn, E. D. Sandoval, S. Kharabadze and A. N. Kolmogorov, _Comput. Phys. Commun._ , 2021, 259, 107679.
  * Kolmogorov _et al._ 2010 A. N. Kolmogorov, S. Shah, E. R. Margine, A. F. Bialon, T. Hammerschmidt and R. Drautz, _Phys. Rev. Lett._ , 2010, 105, 217003.
  * Togo and Tanaka 2015 A. Togo and I. Tanaka, _Scr. Mater._ , 2015, 108, 1–5.
  * Thorn _et al._ 2023 A. Thorn, D. Gochitashvili, S. Kharabadze and A. N. Kolmogorov, _Phys. Chem. Chem. Phys._ , 2023, 25, 22415–22436.
  * Giannozzi _et al._ 2017 P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu and S. Baroni, _J. Phys. Condens. Matter_ , 2017, 29, 465901.
  * van Setten _et al._ 2018 M. J. van Setten, M. Giantomassi, E. Bousquet, M. J. Verstraete, D. R. Hamann, X. Gonze and G.-M. Rignanese, _Comput. Phys. Commun._ , 2018, 226, 39.
  * Methfessel and Paxton 1989 M. Methfessel and A. T. Paxton, _Phys. Rev. B_ , 1989, 40, 3616.
  * Baroni _et al._ 2001 S. Baroni, S. De Gironcoli, A. Dal Corso and P. Giannozzi, _Rev. Mod. Phys._ , 2001, 73, 515.
  * 50 _See Supporting Information for Figs. S1-S19 and Table S1-S4. A description for each figure or table is given in the body of the paper, wherever Ref. [49] is cited._
  * Giustino _et al._ 2007 F. Giustino, M. L. Cohen and S. G. Louie, _Phys. Rev. B_ , 2007, 76, 165108.
  * Poncé _et al._ 2016 S. Poncé, E. R. Margine, C. Verdi and F. Giustino, _Comput. Phys. Commun._ , 2016, 209, 116.
  * Margine and Giustino 2013 E. R. Margine and F. Giustino, _Phys. Rev. B_ , 2013, 87, 024505.
  * Lee _et al._ 2023 H. Lee, S. Poncé, K. Bushick, S. Hajinazar, J. Lafuente-Bartolome, J. Leveillee, C. Lian, J.-M. Lihm, F. Macheda, H. Mori, H. Paudyal, W. H. Sio, S. Tiwar, M. Zacharias, X. Zhang, N. Bonini, E. Kioupakis, E. R. Margine and F. Giustino, _npj Comput. Mater._ , 2023, 9, 156.
  * Marzari _et al._ 2012 N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza and D. Vanderbilt, _Rev. Mod. Phys._ , 2012, 84, 1419.
  * Pizzi _et al._ 2020 G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi and J. R. Yates, _J. Phys. Condens. Matter_ , 2020, 32, 165902.
  * 57 A. Marrazzo, S. Beck, E. R. Margine, N. Marzari, A. A. Mostofi, J. Qiao, I. Souza, S. S. Tsirkin, J. R. Yates and G. Pizzi.
  * Warner _et al._ 2012 J. H. Warner, E. R. Margine, M. Mukai, A. W. Robertson, F. Giustino and A. I. Kirkland, _Science_ , 2012, 337, 209–212.
  * Peng _et al._ 2023 H. Peng, L. Chen, Y. Wang, B. Cao, D. Li and Y. Chen, _ACS Appl. Mater. Interfaces_ , 2023, 15, 25567–25574.
  * Dresselhaus and Dresselhaus 2002 M. S. Dresselhaus and G. Dresselhaus, _Adv. Phys._ , 2002, 51, 1–186.
  * Friend and Yoffe 1987 R. Friend and A. Yoffe, _Adv. Phys._ , 1987, 36, 1–94.
  * Nakamori and ichi Orimo 2003 Y. Nakamori and S. ichi Orimo, _J. Alloys Compd._ , 2003, 370, L7–L9.
  * Wörle and Nesper 1994 M. Wörle and R. Nesper, _J. Alloys Compd._ , 1994, 216, 75–83.
  * Rüdorff and Hofmann 1938 W. Rüdorff and U. Hofmann, _Z. fur Anorg. Allg. Chem._ , 1938, 238, 1–50.
  * Vadlamani _et al._ 2014 B. Vadlamani, K. An, M. Jagannathan and K. S. R. Chandran, _J. Electrochem. Soc._ , 2014, 161, A1731–A1741.
  * Dewhurst _et al._ 2003 J. K. Dewhurst, S. Sharma, C. Ambrosch-Draxl and B. Johansson, _Phys. Rev. B_ , 2003, 68, 020504.
  * An _et al._ 2002 J. M. An, S. Y. Savrasov, H. Rosner and W. E. Pickett, _Phys. Rev. B_ , 2002, 66, 220502.
  * Choi _et al._ 2002 H. J. Choi, D. Roundy, H. Sun, M. L. Cohen and S. G. Louie, _Nature_ , 2002, 418, 758–760.
  * Floris _et al._ 2007 A. Floris, A. Sanna, M. Lüders, G. Profeta, N. N. Lathiotakis, M. A. L. Marques, C. Franchini, E. K. U. Gross, A. Continenza and S. Massidda, _Physica C: Supercond._ , 2007, 456, 45.
  * Sanna _et al._ 2012 A. Sanna, S. Pittalis, J. K. Dewhurst, M. Monni, S. Sharma, G. Ummarino, S. Massidda and E. K. U. Gross, _Phys. Rev. B_ , 2012, 85, 184514.
  * Margine _et al._ 2016 E. R. Margine, H. Lambert and F. Giustino, _Sci. Rep._ , 2016, 6, 21414.
  * Kolmogorov and Curtarolo 2006 A. N. Kolmogorov and S. Curtarolo, _Phys. Rev. B_ , 2006, 74, 224507.
  * Calandra _et al._ 2007 M. Calandra, A. N. Kolmogorov and S. Curtarolo, _Phys. Rev. B_ , 2007, 75, 144506.
  * Pellegrini _et al._ 2022 C. Pellegrini, R. Heid and A. Sanna, _JPhys Materials_ , 2022, 5, 024007.
  * Golubov _et al._ 2002 A. A. Golubov, J. Kortus, O. V. Dolgov, O. Jepsen, Y. Kong, O. K. Andersen, B. J. Gibson, K. Ahn and R. K. Kremer, _J. Phys. Condens. Matter_ , 2002, 14, 1353.
  * Hlinka _et al._ 2003 J. Hlinka, V. Železný, I. Gregora, J. Pokorný, A. M. Fogg, J. B. Claridge, G. R. Darling and M. J. Rosseinsky, _Phys. Rev. B_ , 2003, 68, 220510.
  * Shannon 1976 R. D. Shannon, _Acta Cryst. A_ , 1976, 32, 751–767.
  * Kang _et al._ 2006 K. Kang, Y. S. Meng, J. Bréger, C. P. Grey and G. Ceder, _Science_ , 2006, 311, 977–980.
  * Jain _et al._ 2013 A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder and K. A. Persson, _APL Materials_ , 2013, 1, 011002.
  * Emery _et al._ 2005 N. Emery, C. Hérold, M. d’Astuto, V. Garcia, C. Bellin, J. F. Marêché, P. Lagrange and G. Loupias, _Phys. Rev. Lett._ , 2005, 95, 087003.
  * Delacroix _et al._ 2021 S. Delacroix, F. Igoa, Y. Song, Y. L. Godec, C. Coelho-Diogo, C. Gervais, G. Rousse and D. Portehault, _Inorg. Chem._ , 2021, 60, 4252–4260.
  * Miao _et al._ 2016 R. Miao, G. Huang and J. Yang, _Solid State Commun._ , 2016, 233, 30–34.
  * Albert and Schmitt 1999 B. Albert and K. Schmitt, _Inorg. Chem._ , 1999, 38, 6159–6163.
  * Renker _et al._ 2004 B. Renker, H. Schober, P. Adelmann, P. Schweiss, K.-P. Bohnen and R. Heid, _Physical Review B_ , 2004, 69, 052506.
  * Mori and Takayama-Muromachi 2004 T. Mori and E. Takayama-Muromachi, _Current Applied Physics_ , 2004, 4, 276–279.
  * Jain _et al._ 2011 A. Jain, G. Hautier, S. P. Ong, C. J. Moore, C. C. Fischer, K. A. Persson and G. Ceder, _Phys. Rev. B_ , 2011, 84, 045115.
  * Wang _et al._ 2021 A. Wang, R. Kingsbury, M. McDermott, M. Horton, A. Jain, S. P. Ong, S. Dwaraknath and K. A. Persson, _Scientific Reports_ , 2021, 11, 15496.
  * Aykol _et al._ 2018 M. Aykol, S. S. Dwaraknath, W. Sun and K. A. Persson, _Science Advances_ , 2018, 4, eaaq0148.
  * Stanzione _et al._ 2020 D. Stanzione, J. West, R. T. Evans, T. Minyard, O. Ghattas and D. K. Panda, in _Frontera: The Evolution of Leadership Computing at the National Science Foundation_ , ACM, New York, NY, USA, 2020, pp. 106–111.
  * Boerner _et al._ 2023 T. J. Boerner, S. Deems, T. R. Furlani, S. L. Knuth and J. Towns, Practice and Experience in Advanced Research Computing, New York, NY, USA, 2023, p. 173–176.