# Stability-superconductivity map for compressed Na-intercalated graphite

Shashi B. Mishra  Department of Physics, and Astronomy, Binghamton University-
SUNY, Binghamton, New York 13902, USA  Edan T. Marcial  Department of Physics,
and Astronomy, Binghamton University-SUNY, Binghamton, New York 13902, USA
Suryakanti Debata  Department of Physics, and Astronomy, Binghamton
University-SUNY, Binghamton, New York 13902, USA  Aleksey N. Kolmogorov
Department of Physics, and Astronomy, Binghamton University-SUNY, Binghamton,
New York 13902, USA  Elena R. Margine  Department of Physics, and Astronomy,
Binghamton University-SUNY, Binghamton, New York 13902, USA

###### Abstract

A recent ab initio investigation of Na-C binary compounds under moderate
pressures has uncovered a possible stable NaC4 superconductor with an
estimated critical temperature up to 41 K. We revisit this promising binary
system by performing a more focused exploration of Na-intercalated graphite
configurations, assessing the sensitivity of their thermodynamic stability to
density functional approximations at different $(T,P)$ conditions, and
examining their superconducting properties with the anisotropic Migdal-
Eliashberg formalism. The combinatorial screening of possible Na arrangements
reveals additional stable stoichiometries, i.e., Na3C10, NaC8, NaC10, and
NaC12, that redefine the previously proposed convex hulls for pressures up to
10 GPa. The evaluation of formation enthalpies with different van der Waals
functionals indicates that the proposed compounds might not be
thermodynamically stable at zero temperature but some of them could stabilize
due to the vibrational entropy or form via cold compression if graphite is
used as a starting material. Our more rigorous modeling of the electron-phonon
coupling in NaC4 confirms the material’s potential for high-temperature
superconductivity, with a critical temperature reaching 48 K at 10 GPa, and
reveals a well-defined two-gap structure unusual for an electron-doped
compound. By tracking the position of the intercalant nearly free electron
states with respect to the Fermi level in viable Na-C compounds, we map out
the range of pressures and compositions needed for strong electron-phonon
coupling and identify Na3C10 as an equally promising superconductor.

##  I Introduction

Carbon-based materials offer a rich playground in the quest for
thermodynamically stable superconductors at ambient and low pressure [1, 2, 3,
4, 5, 6, 7]. Graphite intercalation compounds (GICs) [8], in particular, have
been known to superconduct since the mid 1960s when superconductivity below 1
K was discovered in graphite intercalated with alkali metals (AC8 with A = K,
Rb, and Cs) [9]. Since then, considerable efforts have been made to find
superconductivity in GICs at elevated temperatures, and the search has been
focused on increasing the charge transfer from the intercalant atoms to the
host graphene sheets by varying the stoichiometry and the intercalant type. To
date, YbC6 and CaC6 with critical temperature ($T_{\rm c}$) values of 6.5 and
11.5 K stand as the highest temperature superconductors in this class under
ambient conditions [10, 11]. Superconductivity has also been induced or
enhanced by applying external pressure. For instance, LiC2 and NaC2 become
superconductors with $T_{\rm c}=1.9$ K at 3.3 GPa [12] and $T_{\rm c}=5.0$ K
at 3.5 GPa [13], while in CaC6 a maximum value of 15.1 K is reached at 7.5 GPa
[14].

Recently, an ab initio evolutionary structure search in the Na-C system has
identified a new GIC at the 1:4 composition [15]. The NaC4 phase with the oS20
Pearson symbol and Cmcm symmetry (Fig. 1(c)) was predicted to stabilize above
8.9 GPa and have a maximum $T_{\rm c}$ of 41.2 K at 5 GPa. While this work has
shown that the electron-phonon (e-ph) coupling is strong enough to yield a
$T_{\rm c}$ above the record 39 K value in MgB2 based on the Allen-Dynes
modified McMillan formula [16, 17], the nature of the superconducting gap has
not been investigated.

Our present study reexamines the stability and superconductivity of Na-C
candidate materials under pressures up to 10 GPa. Further exploration of the
vast configuration space was motivated by the well-known difficulty of
locating thermodynamically stable compounds even with the most advanced
structure optimization methods [18, 19, 20, 21]. For instance, our
unconstrained evolutionary searches yielded a number of new complex crystal
structures that have either been confirmed experimentally [22, 23, 24] or
inspired the creation of more stable derivatives based on the identified
favorable motifs [25, 26, 27, 28, 29]. In the latest investigation of layered
metal borocarbides, we generated viable quaternary phases via systematic
screening of possible metal decorations within ideal or modified honeycomb
frameworks [28]. Application of this strategy to Na-C in the present work
resulted in a significant revision of the previously proposed ground states.
The updated set at 10 GPa contains large-sized oP104-Na3C10, oP18-NaC8,
oS44-NaC10, and hP13-NaC12 phases that make mP12-NaC2 metastable and oS20-NaC4
marginally stable.

Figure 1: Select high-pressure Na-C phases found to be viable in previous [15]
(a,c) and present (b,d,e,f,g) ab initio studies. The top and side views show
the corresponding distributions of Na (large yellow spheres) intercalated
between C honeycomb layers (small black spheres).

Importantly, we demonstrate that the global stability of these van der Waals
solids is very sensitive to the treatment of the dispersive interactions.
While the Na-C are thermodynamically stable in the optB88-vdW approximation of
the non-local density functional [30], as was shown in Ref. [15], they have
positive formation enthalpies in our optB86b-vdW [31] calculations.
Fortunately, the inclusion of vibrational entropy brings them closer to
stability and makes Na3C10 a true ground state at room temperature. Moreover,
we argue that the feasibility of obtaining the binary layered intercalated
compounds should be evaluated with respect to graphite rather than diamond,
provided that the synthesis starts with the former material and the $(T,P)$
conditions are kept sufficiently low to avoid the $sp^{2}$ to $sp^{3}$
transformation. When referenced to the layered carbon polymorph, the Na-
intercalated variants have negative formation Gibbs free energies.

Finally, we probe the proposed materials for superconducting properties. We
analyze the position of the intercalant nearly free electron states across
different phases and pressures to establish conditions needed for strong e-ph
coupling. According to the constructed composition-pressure map, Na3C10
compound with a complex structure has promise to be a high-$T_{\rm c}$
superconductor. We further study oS20-NaC4 at 10 GPa using the fully
anisotropic Migdal-Eliashberg (aME) formalism [32, 33, 34]. The aME method has
been demonstrated to yield significantly higher estimates of the critical
temperature, by a factor of two to three, in related layered conventional
superconductors, i.e., metal borides and borocarbides [35, 36, 37, 32, 38, 39,
40]. We calculate the momentum dependence of the e-ph coupling and
superconducting energy gap on the Fermi surface. Our aME calculations yield a
$T_{\rm c}$ of 48 K, which is slightly larger than the estimate of 39 K
obtained in the previous work using an isotropic treatment of the e-ph
coupling [15].

##  II Methods

We examined the stability of the Na-C phases with VASP [41] using projector
augmented wave potentials [42]. By default, the non-local van der Waals
interactions were treated with the optB86b-vdW functional [31], but also
checked with alternative optB88-vdW [30] and r2SCAN+rVV10 [43, 44]
functionals. Dense Monckhorst-Pack $k$-point meshes with $\Delta k\sim
2\pi\times 0.025$ Å-1 and a plane-wave cutoff of 500 eV ensured numerical
convergence to typically within 1 meV/atom. To systematically screen possible
decorations of interlayer sites, we constructed supercells with up to 42 atoms
by expanding the stoichiometric MgB2 prototype and sequentially removing metal
atoms. For detection and elimination of equivalent configurations, we relied
on our radial distribution function fingerprint [24]. We calculated
thermodynamic corrections due to the vibrational entropy using the finite
displacement method implemented in PHONOPY [45], employing supercells sized
between 52 and 128 atoms and applying 0.1 Å displacements within the harmonic
approximation. Full structural information for select phases is provided in
the form of CIF files in Supplemental Material [46].

For calculating properties related to superconductivity, we employed the
Quantum ESPRESSO (QE) package [47] with the optB86b-vdW functional [31, 48,
49] and optimized norm-conserving Vanderbilt pseudopotentials (ONCVPSP) [50]
from the Pseudo Dojo library [51] generated with the relativistic Perdew-
Burke-Ernzerhof parametrization [52]. We used a plane-wave cutoff of 100 Ry, a
Methfessel-Paxton smearing [53] value of 0.02 Ry, and $\Gamma$-centered
k-grids of $12\times 12\times 8$ for oS52-Na3C10 and NaC12, $6\times 6\times
4$ for oP104-Na3C10, $16\times 16\times 16$ for NaC4 and NaC8, and $12\times
12\times 12$ for NaC10 to describe the electronic structure. The lattice
parameters and atomic positions were relaxed until the total enthalpy 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 [54] on q-meshes of $3\times 3\times 2$ for oS52-Na3C10, $4\times
4\times 4$ for NaC4 and NaC8, $3\times 3\times 3$ for NaC10, and $4\times
4\times 4$ for NaC12.

Figure 2: Stability of Na-C phases at 10 GPa calculated with (a) the
optB88-vdW functional or (b) the optB86b-vdW functional. The global (local)
convex hulls are denoted with solid (dashed) lines. The formation enthalpies
are shown with and without zero-point energy (ZPE) in blue and gray,
respectively.

To investigate e-ph interactions and superconducting properties of NaC4
further, we used the EPW code [55, 56, 32, 33]. The electronic wavefunctions
required for the Wannier interpolation [57, 58] were obtained on a uniform
$\Gamma$-centered $8\times 8\times 8$ k-grid. We used ten atom-centered
orbitals to describe the electronic structure of $\rm NaC_{4}$, with one $s$
orbital for each Na atom and one $p_{z}$ orbital for each C atom. The
anisotropic Migdal-Eliashberg equations were solved on fine uniform $80\times
80\times 80$ k\- and $40\times 40\times 40$ q-point grids, with an energy
window of $\pm 0.2$ eV around the Fermi level and a Matsubara frequency cutoff
of 1.0 eV. When solving the isotropic ME equations, the Dirac deltas of
electrons and phonons were replaced by Gaussians of width 50 meV and 0.3 meV,
respectively.

##  III Results and Discussions

###  III.1 Structure and stability of Na-C compounds

The global evolutionary searches by Hao et al. [15] provide an important
baseline for viable morphologies in the Na-C binary up to 10 GPa. Namely, only
graphite configurations intercalated with either one or two layers of Na
appear to compete for thermodynamic stability across the composition range. We
reproduce the negative formation enthalpy values evaluated with the optB88-vdW
functional for the reported oS20-NaC4 ($-17.3$ meV/atom) and mP12-NaC2
($-15.9$ meV/atom) phases relative to bcc-Na and diamond at 10 GPa. At the
same time, we observe a high sensitivity of the results to phonon
contributions and DFT approximations due to the starkly different bonding
types displayed by the binary and reference phases. The inclusion of the zero-
point energy (ZPE) alone shifts the values to $-36.2$ meV/atom and $-28.4$
meV/atom, respectively, and the vibrational entropy further stabilizes the
intercalated compounds (see Fig. 2). The beneficial effect of temperature may
be consequential because the formation enthalpies calculated in the
optB86b-vdW approximation are positive, at 28.2 and 23.7 meV/atom, and the
formation Gibbs free energies turn negative only at about 600 K. Considering
the wide range of enthalpy values of graphite relative to diamond, 59, 88, and
116 meV/atom produced by the optB88-vdW, r2SCAN+rVV10, and optB86b-vdW
functionals, respectively, it is evident that the systematic errors are too
large to make a definitive conclusion about the global stability of the Na-C
intercalated phases.

Fortunately, the diamond reference is likely irrelevant if graphite is used as
the staring material in the synthesis of the proposed binary compounds and the
sample is not heated to high temperatures required to induce the $sp^{2}$ to
$sp^{3}$ rebonding. The well-established phase diagram for elemental C [59]
and the experimental data on the kinetics of the transformation [59] indicate
that graphite survives up to $15-30$ GPa at room temperature. The increase in
the interlayer spacing and the transfer of the negative charge to the C layers
caused by Na intercalation may affect the kinetic barriers but are not
expected to appreciably facilitate the rebonding in the targeted $(T,P)$
range. Therefore, we find it more appropriate to define the synthesizability
of the GICs with respect to bcc-Na and graphite rather than diamond. Figure 2
shows that the two proposed phases are indeed thermodynamically favorable with
respect to the C layered polymorph, and the relative enthalpies agree well in
the optB86b-vdW and optB88-vdW approximations due to the similarity of the
elemental and binary morphologies.

Our screening of larger-sized intercalated structures offers further insight
into what phases may form under pressure in this binary system. Starting with
14 different supercells of the MgB2 or mP12-NaC2 prototypes, we generated
$\sim 750$ unique NaxC2 configurations and fully optimized them at 10 GPa. The
considered phases had hexagonal (7%), orthorhombic (48%), monoclinic (38%), or
trigonal (7%) unit cells with up to 52 atoms at 35 compositions below 1:2.

At the 1:4 stoichiometry corresponding to the previously proposed oS20 ground
state, we find an alternative mP20 structure that is nearly degenerate in
Gibbs free energy with oS20 in both optB86b-vdW and optB88-vdW calculations
(see Figs. S1 and S2 in Supplemental Material [46]). The two phases have
fairly uniform metal decorations but the monoclinic variant features a
pronounced interlayer shift that positions Na atoms near hexagon centers in
one of the two sandwiching C layers (Figs. 1(c) and 1(d)). The structural
frustration with respect to sheer in these and other examined Na-C phases
results in soft phonon modes that are difficult to converge.

The dominant compound at 10 GPa actually occurs at an unusual 3:10 composition
in all considered vdW approximations (Figs. 2 and S3 in Supplemental Material
[46]). The best oS52 structure with the $Cmcm$ symmetry in the generated pool
had multiple imaginary modes across the Brillouin zone. The displacement of
the atoms along the only imaginary mode at $\Gamma$ point corresponding to an
interlayer shift lowered the symmetry to oS52 ($C222_{1}$) and the energy by
0.7 meV/atom but did not make the structure dynamically stable (see Fig. S4
[46]). We doubled the shortest lattice constant ($a$) and repeated the
procedure, further lowering the symmetry to oP104 ($P2_{1}2_{1}2_{1}$) and
energy by 0.5 meV/atom. Unfortunately, the phonon dispersion still displays
branches with imaginary frequencies (see Fig. S5 [46]) and additional
distortions leading to, i.e., mP104 ($P2_{1}$), produce negligible gains in
energy. The situation is not uncommon and cannot always be resolved [60, 25,
61]. We consider the oP104 structure an adequate model for the frustrated
compound and use it in the thermodynamic stability analysis (the presence of
slightly imaginary modes has an insignificant effect on the vibrational
entropy). Figure 1(b) shows that, just as in oS20-NaC4, the Na atoms are not
in perfect registry with the honeycomb lattice but a third of them do reside
above and below the C hexagon centers.

The oP104-Na3C10 destabilizes the proposed double-layer mP12-NaC2 phase at 10
GPa but leaves room for possible GICs at the low-Na composition end, such as
NaC8 and NaC10. Figures 1(e) and 1(f) illustrate that the most favored
configurations feature nearly flat C layers and Na intercalating every other
gallery, known as stage-two GICs. The lower Na concentration allows the metal
sublattice to find a better balance between even distribution and registry
with the honeycomb lattice. Consequently, oS44-NaC10 and oP18-NaC8 are found
to be dynamically stable (see Fig. S4 [46]).

Figure 3:  Stability and electronic properties of Na1-xCx compounds calculated
with the optB86b-vdW functional. (a) Convex hulls with respect to bcc Na and
graphite at 2-10 GPa. (b) Stability ranges of the stable Na-C phases. The
shaded area marks the $(x,P)$ phase space in which Na-C materials are not
expected to have electron-doped Na-$s$ states important for superconductivity.
(c) Density of states in oP104-Na3C10. The dashed red line showing the bottom
of the Na-$s$ states. (d) Position of the Na-$s$ band edge for the
intercalated Na-C compounds at 2$-$10 GPa. The dashed red line marks an
approximate composition at which the Na-$s$ band edge is expected to align
with $E_{\rm F}$ at 10 GPa.

The resulting local convex hull (relative to graphite) at 0 K and 10 GPa is
defined by the four phases at the 3:10, 1:4, 1:8, and 1:10 compositions in all
optB86b-vdW, optB88-vdW, and r2SCAN+rVV10 approximations (Figs. 2 and S3
[46]). The set remains unchanged at elevated temperatures once the vibrational
entropy is included (Figs. S1 and S2 in Supplemental Material [46]).
Interestingly, oP104-Na3C10 is the only intercalated phase that is globally
stable at room temperature in both DFT treatments (Fig. 2). To probe what GICs
could be stable in the $1-10$ GPa range, we reoptimized relevant phases at
select pressures and constructed corresponding local convex hulls at 0 K (Fig.
3). The first material to stabilize under compression is actually a stage-two
hP13-NaC12 GIC derived from the $\sqrt{3}\times\sqrt{3}\times 2$ supercell of
MgB2. mP12-NaC2 may indeed form just above 3 GPa (a related LiB compound with
double metal layers has been synthesized via cold compression above 23 GPa
[62]). The four proposed materials between the 3:10 and 1:10 ratios require
higher pressures to be viable (Fig. 3(b)).

###  III.2 Electronic and vibrational properties of Na-C phases

To screen the Na-C binary system for possible high-$T_{\rm c}$
superconductivity, we investigated the electronic and vibrational properties
in the proposed set of intercalated phases, oS52 (oP104)-Na3C10, oS20
(mP20)-NaC4, oP18-NaC8, oS44-NaC10, and hP13-NaC12, at pressures up to 10 GPa.
First, we examined the response of the intercalant Na-$s$ states to
composition and pressure to determine conditions under which the band crosses
the Fermi level. The presence of the partially occupied nearly free electron
states has been linked to superconductivity in several layered compounds
(i.e., CaC6 and YbC6 GICs [63, 64, 65], bilayer C6CaC6 [66, 67], monolayer
LiC6 [68, 69], and LiB [70, 38]). The strong interaction of this band with the
surrounding lattice has been rationalized as a response to both quantum
confinement and electrostatic effects [68, 71, 72].

Our electronic band structure and DOS calculations show that only two of the
considered compounds, Na3C10 and NaC4, have partially occupied Na-$s$ states
in the considered pressure range. Just as in LiB [70], compression from 1 GPa
to 10 GPa raises these nearly free electron states by about 1.5 eV.
Fortunately, the band edge remains at $-1.3$ eV and $-0.7$ eV below the Fermi
level in the two respective Na-C compounds at 10 GPa (Figs. 3(c) and 4). The
data in Figs. 3(b) and 3(d) and Fig. S4 in Supplemental Material [46]
illustrate that none of the carbides with the lower metal content feature the
Na-$s$ states at the Fermi level in the corresponding regions of thermodynamic
stability. For instance, this band would be partially occupied in oP18-NaC8
below 4 GPa but the phase requires over 8 GPa to break the convex hull. To
construct an approximate $x(P)$ boundary between Na1-xCx compounds with and
without populated Na-$s$ states, we estimated the critical concentrations $x$
by linear interpolations at every 2 GPa pressure (the dashed red line across
panels (b) and (d) of Fig. 3 illustrates the determination of the point at 10
GPa). The near-linear shift of the critical compositions with pressure found
for the phases with different symmetries suggests little sensitivity of the
Na-$s$ energies to the particular decoration of the metal sites as long as the
arrangements within galleries are fairly uniform (see Fig. 1). The resulting
shaded area in Fig. 3(b), extrapolated down to the ambient pressure, marks the
$x(P)$ region in which Na-C phases are not expected to have superconductivity
arising primarily from the strong coupling to the intercalant states.

To check the validity of this observation, we evaluated the e-ph coupling by
calculating the Eliashberg spectral function for this set of materials using
the QE package and estimated their $T_{\rm c}$ using the Allen-Dynes equation
[17] with $\mu^{*}=0.1$. We denote the e-ph coupling strength and the critical
temperature obtained with QE as $\lambda^{\rm QE}$ and $T_{\rm c}^{\rm QE}$,
respectively, to differentiate them from the EPW values discussed in the next
section. Unsurprisingly, the three phases in the gray region probed within
their respective stability regions oP18 at 10 GPa, oS44 at 10 GPa, and hP13 at
5 GPa, were found to have low $\lambda^{\rm QE}$ of 0.20, 0.18, and 0.18,
respectively, with no measurable $T_{\rm c}^{\rm QE}$ values.

For oS20-NaC4 at 10 GPa, our calculations produced $\lambda^{\rm QE}$ and
$T_{\rm c}^{\rm QE}$ values of 1.13 and 33.4 K, comparable to 0.99 and 38.8 K
found previously for this material with similar settings except for the use of
a different functional [15]. In addition, the competing mP20-NaC4 structure
has lower $\lambda^{\rm QE}=0.89$ and $T_{\rm c}^{\rm QE}=28.5$ K. Due to the
dynamical instability of the Na3C10 large-sized representative models, we
could obtain only crude estimates of the compound’s superconducting properties
by disregarding all imaginary frequencies in the calculation of
$\alpha^{2}F(\omega)$ and $\omega_{log}$. The relatively high $\lambda^{\rm
QE}$ ($T_{\rm c}^{\rm QE}$) values of 1.39 (37.7 K), and 1.50 (29.6 K)
obtained at 10 GPa for the alternative oS52 ($Cmcm$) and oS52 ($C222_{1}$)
structural models indicate that Na3C10 should be as good a superconductor as
oS20-NaC4. Hence, we believe that the conclusions drawn from our following in-
depth analysis of the superconducting properties in NaC4 should be
transferable to the neighboring Na3C10 compound.

###  III.3 Superconducting properties of NaC4

Figure 4: (a) Electronic band structure and (b) density of states (DOS) in
(states/(eV$\cdot$spin$\cdot$atom)) of NaC4 at 10 GPa with orbital
decomposition. The size of the symbols in (a) is proportional to the
contribution of each orbital character. In (b), the total DOS is represented
by the gray shaded region and the projected DOS corresponding to Na-$s$,
C-$p_{z}$, and C-$sp^{2}$ orbitals is shown with red, green, and blue lines,
respectively. (c) Band structure of C4 structure (solid black) and Na
structure (dashed red) with the same lattice parameters as NaC4. In (c), C4
and Na bands are shifted downward and upward by $\sim$ 2 eV to match the Fermi
level of NaC4. (d) Charge density difference calculated using
$\Delta\rho(r)=\rho_{\rm NaC_{4}}(r)-\rho_{\rm C_{4}}(r)-\rho_{\rm Na}(r)$ and
generated using VESTA [73]. The red and green lobes represent the charge
accumulation and depletion regions, respectively, with an isosurface value set
to $7~{}\times 10^{-2}\,e/$Å3.

With the oS20-NaC4 phase determined to stabilize above 6.5 GPa, we perform a
more detailed examination of its properties at 10 GPa. A noticeable feature in
the electronic structure of NaC4 in Figs. 4(a) and 4(b) is the presence of a
van Hove singularity (vHS) peak near the Fermi level. This gives rise to a
sizable DOS at $E_{\rm F}$, with $N(E_{\rm F})=0.10$
states/(eV$\cdot$spin$\cdot$atom), a value closer to that of 0.11
states/(eV$\cdot$spin$\cdot$atom) reported for CaC6 [74]. The orbital DOS
decomposition shows that nearly 85% of $N(E_{\rm F})$ comes from the anti-
bonding C-$\pi^{*}$ states, while the remaining 15% arises predominantly from
the Na-$s$ states. As in CaC6, there is no notable contribution of the
C-$\sigma$ states to DOS at $E_{\rm F}$.

More insight into the electronic structure can be obtained by analyzing the
dispersion of NaC4 with Na atoms removed (empty graphite-like C4 structure,
solid black lines) and the one with C atoms removed (empty Na metal structure,
red dashed lines) depicted in Fig. 4(c). Comparing Figs. 4(a) and 4(c), we can
interpret the band structure of NaC4 as a superposition of the C4 and Na
bands. To facilitate this comparison, the C and Na bands in Fig. 4(c) are
shifted downward and upward by $\sim$2 eV, respectively, to match the Fermi
level in NaC4. Effectively, the relative shift is equivalent to a charge
transfer from the Na atoms to the graphene layers. This can be clearly
visualized from the charge density difference plot in Fig. 4(d), where the red
and green lobes represent the charge accumulation and depletion regions around
the C and Na atoms, respectively. Integrating the atomic-projected DOS, we
extract the Löwdin charge for each atom. In NaC4, we find a charge transfer of
0.21 electrons/C atom, whereas in CaC6 each C atom only gains 0.11 electrons
[63].

Figure 5: (a) Fermi surface (FS) of NaC4 at 10 GPa in the first Brillouin zone
generated with FermiSurfer [75]. (b$-$d) Band decomposed FS corresponding to
the three bands crossing $E_{\rm F}$.

Figures 5(a)-5(d) illustrate the total and band-decomposed Fermi surface (FS).
Comparison with the band structure in Fig. 4(a) reveals that the three bands
crossing the Fermi level give rise to a multi-sheet FS structure. The FS
sheets corresponding to band 1 (shown in blue) and band 2 (shown in green)
exhibit 3D nested shapes and originate from the C-$p_{z}$ states, as already
pointed out in Ref. [15]. These nested regions surround two bean-shaped FS
sheets (shown in red) centered around $\Gamma$ of predominant Na-$s$ character
(band 3).

Figure 6: (a) Phonon band structure of NaC4 at 10 GPa with decomposition into
in-plane and out-of-plane Na and C vibrations, respectively. The size of the
symbols is proportional to the amount of vibration. (b) Total phonon density
of states (PhDOS) (states/(meV$\cdot$atom)) shown as the gray-shaded region
along with atomic-projected PhDOS. (c) The isotropic Eliashberg spectral
function $\alpha^{2}F(\omega)$ and the cumulative electron-phonon coupling
strength $\lambda(\omega)$.

In Figs. 6(a) and 6(b), we plot the phonon dispersion and the atomic-projected
phonon density of states (PhDOS). Similar to CaC6 and several other GICs [63,
76, 74], the phonon dispersion can be divided into three regions: a low-
frequency region composed of mostly Na-related vibrations mixed with out-of-
plane Cz vibrations ($<30$ meV), an intermediate frequency region dominated by
out-of-plane Cz and Naz vibrations ($30-80$ meV), and a high-frequency region
of in-plane Cxy vibrations ($>100$ meV).

Figure 7: (a) Distribution of the electron-phonon coupling strength $\lambda_{\bm{k}}$. (b) Momentum-resolved electron-phonon coupling $\lambda_{\bm{k}}$ on the Fermi surface. Figure 8: (a) Isotropic superconducting gap of NaC4 at 10 GPa as a function of temperature according to the FSR (gray), FBW (blue) and FBW+$\mu$ (red) implementations of the Migdal-Eliashberg equations [33] for two different sets of Coulomb pseudopotential $\mu^{*}=0.1$ and 0.2. (b) and (c) Energy distribution of the anisotropic superconducting gap $\Delta_{\bm{k}}$ as a function of temperature for FSR (gray), FBW (blue) and FBW+$\mu$ (red) approaches for $\mu^{*}=0.1$ and 0.2, respectively. The black dashed lines in (b) and (c) represent a guide for the eye highlighting the two superconducting gap distributions. $N(E_{\rm F})$ |  $\omega_{log}$ |  $\lambda$ |  $\mu^{*}$ |  Semi-empirical $T_{\rm c}$ (K) |  Isotropic ME $T_{\rm c}$ (K) |  Anisotropic ME $T_{\rm c}$ (K)
---|---|---|---|---|---|---
(states/(eV$\cdot$spin$\cdot$atom)) | (meV) |  |  |  AD |  ML |  FSR |  FBW |  FBW+$\rm\mu$ |  FSR |  FBW |  FBW+$\rm\mu$
0.10 | 33.85 | 1.19 | 0.1 | 34.9 | 43.3 | 43.5 | 45.5 | 45.5 | 47.0 | 48.2 | 48.2
0.2 | 22.0 | 26.7 | 31.2 | 33.0 | 33.0 | 35.0 | 36.0 | 36.0
Table 1: Summary of superconducting properties of NaC4 at 10 GPa. The table
lists the DOS at the Fermi level $N(E_{\rm F})$, the logarithmic phonon
frequency $\omega_{\rm log}$, the e-ph coupling strength $\lambda$, the
Coulomb parameter $\mu^{*}$, and the superconducting critical temperature
$T_{\rm c}$. The $T_{\rm c}$ values are calculated using different
methodologies: the semi-empirical Allen-Dynes (AD) formula, the machine-
learned (ML) model, the isotropic and anisotropic Fermi surface restricted
(FSR) approximation, and the isotropic and anisotropic full bandwidth
implementation without (FBW) and with (FBW+$\mu$) variation of the chemical
potential $\mu$.

To investigate the e-ph interaction, we compute the isotropic Eliashberg
spectral function $\alpha^{2}F(\omega)$ and the cumulative e-ph coupling
strength $\lambda(\omega)$ shown in Fig. 6(c). Consistent with Ref. [15], we
find that low- and intermediate-frequency phonons play a key role in achieving
a high e-ph coupling in NaC4, and together they account for approximately 83%
of the total $\lambda=1.19$. Moreover, the relative contributions of these two
sets of modes are very similar, with the e-ph coupling strengths associated
primarily with Na and out-of-plane Cz vibrations making up 41% and 42% of
$\lambda$, respectively. The corresponding breakdown of the total e-ph
coupling strength into the three regions is illustrated by the horizontal
dashed lines in Fig. 6(c). Notably, the high-frequency in-plane Cxy modes
exhibit a significantly smaller contribution of only about 17%. This behavior
is similar to other GICs (i.e., CaC6, BaC6, SrC6 [63, 76, 77, 67]) and Ca-
intercalated bilayer graphene [67], where the largest contribution to
$\lambda$ has been associated with the modes below 80 meV.

The anisotropy in the e-ph coupling is quantified by evaluating the momentum-
resolved $\lambda_{\bm{k}}$ as defined in Ref. [32]. Figure 7(a) clearly
displays two peaks in the distribution of $\lambda_{\bm{k}}$, with a sharp
peak centered around 1 and a small peak centered around 3. To further analyze
the origin of these peaks, the variation of $\lambda_{\bm{k}}$ on the Fermi
surface is shown in Fig. 7(b). According to the color map, the first sharp
peak originates from the nested FS regions, while the second smaller peak
arises from the $\Gamma$-centered FS regions. In particular, by comparing with
the band structure and the band decomposed FS plots in Figs. 4(a) and 5, it
becomes evident that the first peak arises from the coupling of the phonon
modes with the C-$\pi^{*}$ states (band 1 and band 2), while the second peak
is mostly due to the intercalant Na-$s$ states (band 3) with some contribution
from the C-$\pi^{*}$ states (band 2).

In order to determine the superconducting critical temperature, we solve the
isotropic and anisotropic Migdal-Eliashberg equations [32, 33]. Since the
electronic DOS strongly peaks close to $E_{\rm F}$, we examine the
superconducting gap using both the Fermi surface restricted (FSR) and the full
bandwidth (FBW) approaches [33, 34]. In FSR, the DOS around $E_{\rm F}$ is
assumed constant, whereas in FBW the full energy dependence of the DOS is
considered, allowing the inclusion of e-ph scattering processes away from
$E_{\rm F}$. A recent study has demonstrated the importance of the FBW
treatment in the high-pressure H3S and the low-pressure BaSiH8
superconductors, which exhibit distinct DOS features that deviate
significantly from the constant-DOS assumption made in the FSR approximation
[34]. For a detailed description of the FSR and FBW approaches and their
implementation in the EPW code, please see Refs. [32, 56, 33, 34].

Figure 8(a) shows the isotropic superconducting gap $\Delta$ as a function of
temperature, calculated with the FSR, FBW, and FBW+$\mu$ methods using a
Coulomb parameter $\mu^{*}$ = 0.1 and 0.2. In FBW+$\mu$ approach, we update
the chemical potential self-consistently as discussed in Ref. [34]. Our
calculations yield $T_{\rm c}$ values of 43.5 (31.2) K and 45.5 (33.0) K, and
a zero-temperature superconducting gap of 7.62 (5.36) meV and 8.24 (5.70) meV
for FSR and FBW with $\mu^{*}$ = 0.1 (0.2), respectively. Furthermore, the
superconducting gap versus temperature curves for FBW and FBW+$\mu$ are
basically identical, thus leading to the same $T_{\rm c}$ values. This can be
understood as the Fermi level in NaC4 lies slightly away from the vHS peak, as
shown in Fig. 4(b), and therefore the DOS with and without updating $\mu$ are
very similar. For comparison, the commonly employed Allen-Dynes (AD) modified
McMillan formula [16, 17] and the machine-learned (ML) SISSO model [78] give
$T_{\rm c}$ values of 34.9 (22.0) K and 43.3 (26.7) K for $\mu^{*}$ = 0.1
(0.2), respectively, with the latter producing results that match very well
our isotropic predictions.

Given the anisotropy found in $\lambda_{\bm{k}}$, we calculate the anisotropic
superconducting gap $\Delta_{\bm{k}}$ as a function of temperature shown in
Figs. 8(b) and 8(c). We find two distinct gaps with mean values
$\Delta_{1}=7.25$ (5.00) meV and $\Delta_{2}=10.42$ (7.69) meV for $\mu^{*}$ =
0.1 (0.2) in the zero-temperature limit using the FSR approach. For
comparison, CaC6 displays a continuum anisotropic gap without separation into
distinct gaps [67, 77, 74]. The $\Delta_{1}$ gap is highly anisotropic and
varies over a wide energy range of 1.75 meV. In contrast, the $\Delta_{2}$ gap
is less pronounced (i.e., fewer energy states contribute to this gap) and has
a much narrower distribution of only 0.75 meV. The superconducting $T_{\rm c}$
is estimated to be 47.0 (35.0) K for $\mu^{*}$ = 0.1 (0.2), about 8 % (12 %)
higher than the isotropic FSR value. Application of the FBW and FBW+$\mu$
treatment only shifts slightly the gap distribution to higher energies while
retaining the overall shape at each temperature. The resulting $T_{\rm c}$
using the two FBW approaches is estimated to be 48.2 (36.0) K with $\mu^{*}$ =
0.1 (0.2). The $T_{\rm c}$ values obtained with different approaches are
summarized in Table 1.

Figure 9: (a) Momentum-resolved superconducting gap $\Delta_{\bm{k}}$ on the
Fermi surface and (b) normalized quasiparticle DOS for NaC4 at 12 K and 10 GPa
using the FSR approach with $\mu^{*}=0.1$. The data points in (b) correspond
to electronic states within $\pm$ 0.2 eV from the Fermi energy.

Figure 9(a) shows the momentum-resolved superconducting gap $\Delta_{\bm{k}}$
on the FS at 12 K within the FSR approximation using $\mu^{*}=0.1$. Comparing
this with Fig. 7(b), the regions of the FS with the lower gap
$\Delta_{1}=7.35$ meV can be associated with the C-$\pi^{*}$ states (band 1
and band 2). In contrast, the upper gap $\Delta_{2}=10.60$ meV primarily
arises from the Na-$s$ states (band 3) with some contribution from the
C-$\pi^{*}$ states (band 2). Similarly, in CaC6, the smallest and highest
values of the smeared multi-gap structure have been associated with the
C-$\pi^{*}$ and Ca-$s,d_{z^{2}}$ states [77, 74]. The two-gap structure of
NaC4 is also evident from the corresponding quasiparticle DOS plotted in Fig.
9(b).

Figure 10: (a) Total DOS (states/(eV$\cdot$spin$\cdot$atom)) of NaC4 in the
energy range of $-$0.2 to 0.1 eV. (b) Superconducting critical temperature
$T_{\rm c}$ as a function of Fermi level shift $\Delta\varepsilon_{\rm F}$
obtained with the AD formula, the ML model, and the isotropic and anisotropic
ME equations in the FSR and FBW approaches with a $\mu^{*}=0.1$.

Finally, we explore the effect of electron and hole doping on the $T_{\rm c}$
in the light of the vHS peak present in the DOS of NaC4 near the Fermi level
(see Fig. 10(a)). Previous theoretical studies on H3S have shown an
enhancement in $T_{\rm c}$ when the Fermi level approaches the maximum of the
vHS peak [79, 80, 34]. To investigate the dependence of $T_{\rm c}$ with
doping, we shift the Fermi level by $\Delta\varepsilon_{\rm F}$ in the $-$0.2
eV to 0.1 eV range in increments of 0.05 eV, corresponding to changes in the
electron number of $\sim$ 0.076. Figure 10(b) shows the $T_{\rm c}$ values
obtained by solving the isotropic and anisotropic ME equations within the FSR
and FBW approaches, and those obtained with the semi-empirical AD and ML
formulas. Overall, all treatments produce a similar trend for the variation of
$T_{\rm c}$ with doping. In particular, the $T_{\rm c}$ versus
$\Delta\varepsilon_{\rm F}$ curve obtained with the anisotropic approach
resembles the shape dependence of the DOS, the $T_{\rm c}$ reaching a maximum
value at $\Delta\varepsilon_{\rm F}=-0.1$ eV when the Fermi level lies in very
close proximity to the vHS peak. For comparison, the effect of doping on
$T_{\rm c}$ is less pronounced when employing the isotropic ME treatment or
the semi-empirical formulas, with the AD formula consistently underestimating
the $T_{\rm c}$ values by about 7 K.

##  IV Summary

We have explored the stability and superconductivity of Na-intercalated
graphite compounds under moderate pressures using ab initio methods. Through
systematic screening, we identified new stable stoichiometries, Na3C10, NaC8,
NaC10, and NaC12, which redefine the previously established convex hulls up to
10 GPa. While some compounds may lack stability at zero temperature, they
could be high-temperature ground states or form through cold compression of
graphite. Our analysis of the nearly free electron states important for
superconductivity in intercalated compounds narrows down the promising
candidates to two. One of them, NaC4 proposed previously by Hao et al. [15],
has been re-analyzed within the anisotropic Migdal-Eliashberg theory and
predicted to be a two-gap superconductor with a $T_{\rm c}$ of 48 K, 24% above
the prior isotropic Eliashberg estimate. The other compound, Na3C10 introduced
in the present study, appears to be prone to interlayer shifts and may develop
domain walls to become dynamically stable. According to our estimates, it
should have $T_{\rm c}$ at least comparable to that in NaC4 at 10 GPa. The ab
initio findings reveal the Na-C system may host viable compounds with
unexpectedly high $T_{\rm c}$ for electron-doped materials.

###### Acknowledgements.

S.B.M. and E.R.M. acknowledge support from the National Science Foundation
(NSF) Award No. DMR-2035518. E.T.M. and A.N.K. acknowledge support from the
NSF Award No. DMR-2320073. 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
[81]. This study 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 [82], which is
supported by NSF grants #2138259, #2138286, #2138307, #2137603, and #2138296.

## References

  * Jishi and Dresselhaus [1992] R. A. Jishi and M. S. Dresselhaus, Superconductivity in graphite intercalation compounds, Phys. Rev. B 45, 12465 (1992).
  * Ekimov  _et al._ [2004] E. Ekimov, V. Sidorov, E. Bauer, N. Mel’Nik, N. Curro, J. Thompson, and S. Stishov, Superconductivity in diamond, Nature 428, 542 (2004).
  * Okazaki  _et al._ [2015] H. Okazaki, T. Wakita, T. Muro, T. Nakamura, Y. Muraoka, T. Yokoya, S.-i. Kurihara, H. Kawarada, T. Oguchi, and Y. Takano, Signature of high Tc above 25 K in high quality superconducting diamond, Appl. Phys. Lett. 106, 052601 (2015).
  * Heguri  _et al._ [2015] S. Heguri, N. Kawade, T. Fujisawa, A. Yamaguchi, A. Sumiyama, K. Tanigaki, and M. Kobayashi, Superconductivity in the graphite intercalation compound BaC6, Phys. Rev. Lett. 114, 247201 (2015).
  * Zhou  _et al._ [2021] H. Zhou, T. Xie, T. Taniguchi, K. Watanabe, and A. F. Young, Superconductivity in rhombohedral trilayer graphene, Nature 598, 434 (2021).
  * Bhaumik  _et al._ [2017] A. Bhaumik, R. Sachan, S. Gupta, and J. Narayan, Discovery of high-temperature superconductivity (Tc = 55 K) in B-doped Q-Carbon, ACS Nano 11, 11915 (2017).
  * Emery  _et al._ [2008] N. Emery, C. Hérold, J.-F. Marêché, and P. Lagrange, Synthesis and superconducting properties of CaC6, Sci. Technol. Adv. Mater. 9, 044102 (2008).
  * Dresselhaus and Dresselhaus [1981] M. Dresselhaus and G. Dresselhaus, Intercalation compounds of graphite, Adv. Phys. 30, 139 (1981).
  * Hannay  _et al._ [1965] N. B. Hannay, T. H. Geballe, B. T. Matthias, K. Andres, P. Schmidt, and D. MacNair, Superconductivity in graphitic compounds, Phys. Rev. Lett. 14, 225 (1965).
  * Weller  _et al._ [2005] T. E. Weller, M. Ellerby, S. S. Saxena, R. P. Smith, and N. T. Skipper, Superconductivity in the intercalated graphite compounds C6Yb and C6Ca, Nat. Phys. 1, 39 (2005).
  * 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, Superconductivity of bulk CaC6, Phys. Rev. Lett. 95, 087003 (2005).
  * Belash  _et al._ [1989] I. Belash, A. Bronnikov, O. Zharikov, and A. Pal’nichenko, Superconductivity of graphite intercalation compound with lithium C2Li, Solid State Commun. 69, 921 (1989).
  * Belash  _et al._ [1987] I. Belash, A. Bronnikov, O. Zharikov, and A. Palnichenko, On the superconductivity of graphite intercalation compounds with sodium, Solid State Commun. 64, 1445 (1987).
  * Gauzzi  _et al._ [2007] A. Gauzzi, S. Takashima, N. Takeshita, C. Terakura, H. Takagi, N. Emery, C. Hérold, P. Lagrange, and G. Loupias, Enhancement of superconductivity and evidence of structural instability in intercalated graphite C6Ca under high pressure, Phys. Rev. Lett. 98, 067002 (2007).
  * Hao  _et al._ [2023] C.-M. Hao, X. Li, A. R. Oganov, J. Hou, S. Ding, Y. Ge, L. Wang, X. Dong, H.-T. Wang, G. Yang, X.-F. Zhou, and Y. Tian, Superconductivity in compounds of sodium-intercalated graphite, Phys. Rev. B 108, 214507 (2023).
  * McMillan [1968] W. L. McMillan, Transition temperature of strong-coupled superconductors, Phys. Rev. 167, 331 (1968).
  * Allen and Dynes [1975] P. B. Allen and R. C. Dynes, Transition temperature of strong-coupled superconductors reanalyzed, Phys. Rev. B 12, 905 (1975).
  * Oganov [2011] A. R. Oganov,  _Modern Methods of Crystal Structure Prediction_ (John Wiley & Sons, 2011).
  * Oganov  _et al._ [2019] A. R. Oganov, C. J. Pickard, Q. Zhu, and R. J. Needs, Structure prediction drives materials discovery, Nat Rev Mater 4, 331 (2019).
  * Wang and Ma [2014] Y. Wang and Y. Ma, Perspective: Crystal structure prediction at high pressures, The Journal of Chemical Physics 140, 040901 (2014).
  * Boeri and et al. [2022] L. Boeri and et al., The 2021 room-temperature superconductivity roadmap, J. Phys.: Condens. Matter 34, 183002 (2022).
  * Kolmogorov  _et al._ [2010] A. N. Kolmogorov, S. Shah, E. R. Margine, A. F. Bialon, T. Hammerschmidt, and R. Drautz, New superconducting and semiconducting Fe-B compounds predicted with an ab initio evolutionary search, Phys. Rev. Lett. 105, 217003 (2010).
  * Kolmogorov  _et al._ [2012] A. N. Kolmogorov, S. Shah, E. R. Margine, A. K. Kleppe, and A. P. Jephcoat, Pressure-driven evolution of the covalent network in CaB6, Phys. Rev. Lett. 109, 075501 (2012).
  * Hajinazar  _et al._ [2021] S. Hajinazar, A. Thorn, E. D. Sandoval, S. Kharabadze, and A. N. Kolmogorov, MAISE: Construction of neural network interatomic models and evolutionary structure optimization, Comput. Phys. Commun. 259, 107679 (2021).
  * Shah and Kolmogorov [2013] S. Shah and A. N. Kolmogorov, Stability and superconductivity of Ca-B phases at ambient and high pressure, Phys. Rev. B 88, 014107 (2013).
  * Kharabadze  _et al._ [2022] S. Kharabadze, A. Thorn, E. A. Koulakova, and A. N. Kolmogorov, Prediction of stable Li-Sn compounds: boosting ab initio searches with neural network potentials, npj Comput Mater 8, 1 (2022).
  * Singh  _et al._ [2022] J. Singh, A. Behatha, S. Kharabadze, A. N. Kolmogorov, G. Vaitheeswaran, and V. Kanchana, Prediction of ground state structures and robust Weyl fermionic states in MnRhP, J. Phys. Chem. C 126, 17328 (2022).
  * Tomassetti  _et al._ [2024a] C. R. Tomassetti, G. P. Kafle, E. T. Marcial, E. R. Margine, and A. N. Kolmogorov, Prospect of high-temperature superconductivity in layered metal borocarbides, J. Mater. Chem. C 12, 4870 (2024a).
  * Kharabadze  _et al._ [2023] S. Kharabadze, M. Meyers, C. R. Tomassetti, E. R. Margine, I. I. Mazin, and A. N. Kolmogorov, Thermodynamic stability of Li–B–C compounds from first principles, Phys. Chem. Chem. Phys. 25, 7344 (2023).
  * Klimeš  _et al._ [2009] J. Klimeš, D. R. Bowler, and A. Michaelides, Chemical accuracy for the van der waals density functional, J. Phys.: Condensed Matter 22, 022201 (2009).
  * Klimeš  _et al._ [2011] J. Klimeš, D. R. Bowler, and A. Michaelides, Van der waals density functionals applied to solids, Phys. Rev. B 83, 195131 (2011).
  * Margine and Giustino [2013] E. R. Margine and F. Giustino, Anisotropic Migdal-Eliashberg theory using Wannier functions, Phys. Rev. B 87, 024505 (2013).
  * 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. Tiwari, M. Zacharias, X. Zhang, N. Bonini, E. Kioupakis, E. R. Margine, and F. Giustino, Electron–phonon physics from first principles using the EPW code, npj Comput. Mater. 9, 2057 (2023).
  * Lucrezi  _et al._ [2024] R. Lucrezi, P. P. Ferreira, S. Hajinazar, H. Mori, H. Paudyal, E. R. Margine, and C. Heil, Full-bandwidth anisotropic migdal-eliashberg theory and its application to superhydrides, Commun. Phys. 7, 33 (2024).
  * Liu  _et al._ [2001] A. Y. Liu, I. I. Mazin, and J. Kortus, Beyond Eliashberg Superconductivity in MgB2: Anharmonicity, Two-Phonon Scattering, and Multiple Gaps, Phys. Rev. Lett. 87, 087005 (2001).
  * Choi  _et al._ [2002] H. J. Choi, D. Roundy, H. Sun, M. L. Cohen, and S. G. Louie, The origin of the anomalous superconducting properties of MgB2, Nature 418, 758 (2002).
  * Choi  _et al._ [2009] H. J. Choi, S. G. Louie, and M. L. Cohen, Prediction of superconducting properties of CaB2 using anisotropic eliashberg theory, Phys. Rev. B 80, 064503 (2009).
  * Kafle  _et al._ [2022] G. P. Kafle, C. R. Tomassetti, I. I. Mazin, A. N. Kolmogorov, and E. R. Margine, Ab initio study of Li-Mg-B superconductors, Phys. Rev. Mater. 6, 084801 (2022).
  * Wang  _et al._ [2023] J. Wang, M. Wang, X. Liu, M. Jiang, and L. Liu, Covalent bond inducing strong electron-phonon coupling superconductivity in MgB2-type transition metal diboride WB2, Phys. Rev. Mater. 7, 074804 (2023).
  * Tomassetti  _et al._ [2024b] C. R. Tomassetti, G. P. Kafle, E. T. Marcial, E. R. Margine, and A. N. Kolmogorov, Prospect of high-temperature superconductivity in layered metal borocarbides, J. Mater. Chem. C 12, 4870 (2024b).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994).
  * Furness  _et al._ [2020] J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, and J. Sun, Accurate and Numerically Efficient r2SCAN Meta-Generalized Gradient Approximation, J. Phys. Chem. Lett. 11, 8208 (2020).
  * 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, Workhorse minimally empirical dispersion-corrected density functional with tests for weakly bound systems: r2SCAN+rVV10, Phys. Rev. B 106, 075422 (2022).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scr. Mater. 108, 1 (2015).
  * [46] See Supplemental Material for detailed figures and structural information. Figures S1 and S2 illustrate the distance to the convex hull for selected Na-C phases using optB86b-vdW and optB88-vdW functionals. Figure S3 presents the stability analysis using r2SCAN+rVV10 functional. Figure S4 shows the electronic structure and phonon dispersion of selected stable structures calculated using Quantum ESPRESSO. Figure S5 shows the phonon dispersion of oP104 calculated using PHONOPY and VASP. CIF files for select Na-C phases are also included.
  * Giannozzi and et al. [2017] P. Giannozzi and et al., Advanced capabilities for materials modelling with Quantum ESPRESSO, J. Phys.: Condens. Matter 29, 465901 (2017).
  * Thonhauser  _et al._ [2015] T. Thonhauser, S. Zuluaga, C. A. Arter, K. Berland, E. Schröder, and P. Hyldgaard, Spin signature of nonlocal correlation binding in metal-organic frameworks, Phys. Rev. Lett. 115, 136402 (2015).
  * Berland  _et al._ [2015] K. Berland, V. R. Cooper, K. Lee, E. Schröder, T. Thonhauser, P. Hyldgaard, and B. I. Lundqvist, van der Waals forces in density functional theory: a review of the vdW-DF method, Rep. Prog. Phys. 78, 066501 (2015).
  * Hamann [2013] D. R. Hamann, Optimized norm-conserving vanderbilt pseudopotentials, Phys. Rev. B 88, 085117 (2013).
  * van Setten  _et al._ [2018] M. van Setten, M. Giantomassi, E. Bousquet, M. Verstraete, D. Hamann, X. Gonze, and G.-M. Rignanese, The PseudoDojo: Training and grading a 85 element optimized norm-conserving pseudopotential table, Comput. Phys. Commun. 226, 39 (2018).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Methfessel and Paxton [1989] M. Methfessel and A. T. Paxton, High-precision sampling for brillouin-zone integration in metals, Phys. Rev. B 40, 3616 (1989).
  * Baroni  _et al._ [2001] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Phonons and related crystal properties from density-functional perturbation theory, Rev. Mod. Phys. 73, 515 (2001).
  * Giustino  _et al._ [2007] F. Giustino, M. L. Cohen, and S. G. Louie, Electron-phonon interaction using Wannier functions, Phys. Rev. B 76, 165108 (2007).
  * Poncé  _et al._ [2016] S. Poncé, E. Margine, C. Verdi, and F. Giustino, EPW: Electron–phonon coupling, transport and superconducting properties using maximally localized Wannier functions, Comput. Phys. Commun. 209, 116 (2016).
  * Marzari  _et al._ [2012] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Maximally localized Wannier functions: Theory and applications, Rev. Mod. Phys. 84, 1419 (2012).
  * Pizzi and et al. [2020] G. Pizzi and et al., Wannier90 as a community code: new features and applications, J. Phys. Condens. Matter 32, 165902 (2020).
  * Sundqvist [2021] B. Sundqvist, Carbon under pressure, Phys. Rep. 909 (2021).
  * Kolmogorov and Curtarolo [2006] A. N. Kolmogorov and S. Curtarolo, Prediction of different crystal structure phases in metal borides: A lithium monoboride analog to MgB2, Phys. Rev. B 73, 180501(R) (2006).
  * Sandoval  _et al._ [2016] E. D. Sandoval, S. Hajinazar, and A. N. Kolmogorov, Stability of two-dimensional BN-Si structures, Phys. Rev. B 94, 094105 (2016).
  * Kolmogorov  _et al._ [2015] A. N. Kolmogorov, S. Hajinazar, C. Angyal, V. L. Kuznetsov, and A. P. Jephcoat, Synthesis of a predicted layered LiB via cold compression, Phys. Rev. B 92, 144110 (2015).
  * Calandra and Mauri [2005] M. Calandra and F. Mauri, Theoretical explanation of superconductivity in C6Ca, Phys. Rev. Lett. 95, 237002 (2005).
  * Csányi  _et al._ [2005] G. Csányi, P. B. Littlewood, A. H. Nevidomskyy, C. J. Pickard, and B. D. Simons, The role of the interlayer state in the electronic structure of superconducting graphite intercalated compounds, Nat. Phys. 1, 42 (2005).
  * Mazin [2005] I. I. Mazin, Intercalant-driven superconductivity in YbC6 and CaC6, Phys. Rev. Lett. 95, 227001 (2005).
  * Mazin and Balatsky [2010] I. Mazin and A. Balatsky, Superconductivity in Ca-intercalated bilayer graphene, Philos. Mag. Lett. 90, 731 (2010).
  * Margine  _et al._ [2016] E. Margine, H. Lambert, and F. Giustino, Electron-phonon interaction and pairing mechanism in superconducting Ca-intercalated bilayer graphene, Sci. Rep. 6, 21414 (2016).
  * Profeta  _et al._ [2012] G. Profeta, M. Calandra, and F. Mauri, Phonon-mediated superconductivity in graphene by lithium deposition, Nat. Phys. 8, 131 (2012).
  * Zheng and Margine [2016] J.-J. Zheng and E. R. Margine, First-principles calculations of the superconducting properties in Li-decorated monolayer graphene within the anisotropic migdal-eliashberg formalism, Phys. Rev. B 94, 064509 (2016).
  * Calandra  _et al._ [2007] M. Calandra, A. N. Kolmogorov, and S. Curtarolo, Search for high ${T}_{c}$ in layered structures: The case of LiB, Phys. Rev. B 75, 144506 (2007).
  * Boeri  _et al._ [2007] L. Boeri, G. B. Bachelet, M. Giantomassi, and O. K. Andersen, Electron-phonon interaction in graphite intercalation compounds, Phys. Rev. B 76, 064510 (2007).
  * Margine and Crespi [2006] E. R. Margine and V. H. Crespi, Universal behavior of nearly free electron states in carbon nanotubes, Phys. Rev. Lett. 96, 196803 (2006).
  * Momma and Izumi [2011] K. Momma and F. Izumi, VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data, J. Appl. Cryst. 44, 1272 (2011).
  * Sanna  _et al._ [2012] A. Sanna, S. Pittalis, J. K. Dewhurst, M. Monni, S. Sharma, G. Ummarino, S. Massidda, and E. K. U. Gross, Phononic self-energy effects and superconductivity in CaC6, Phys. Rev. B 85, 184514 (2012).
  * Kawamura [2019] M. Kawamura, FermiSurfer: Fermi-surface viewer providing multiple representation schemes, Comput. Phys. Commun. 239, 197 (2019).
  * Calandra and Mauri [2006] M. Calandra and F. Mauri, Possibility of superconductivity in graphite intercalated with alkaline earths investigated with density functional theory, Phys. Rev. B 74, 094507 (2006).
  * Sanna  _et al._ [2007] A. Sanna, G. Profeta, A. Floris, A. Marini, E. K. U. Gross, and S. Massidda, Anisotropic gap of superconducting CaC6: A first-principles density functional calculation, Phys. Rev. B 75, 020511(R) (2007).
  * Xie  _et al._ [2022] S. Xie, Y. Quan, A. Hire, B. Deng, J. DeStefano, I. Salinas, U. Shah, L. Fanfarillo, J. Lim, J. Kim, _et al._ , Machine learning of superconducting critical temperature from Eliashberg theory, npj Comput. Mater. 8, 14 (2022).
  * Sano  _et al._ [2016] W. Sano, T. Koretsune, T. Tadano, R. Akashi, and R. Arita, Effect of Van Hove singularities on high-${T}_{\mathrm{c}}$ superconductivity in ${\mathrm{H}}_{3}\mathrm{S}$, Phys. Rev. B 93, 094525 (2016).
  * Quan and Pickett [2016] Y. Quan and W. E. Pickett, Van Hove singularities and spectral smearing in high-temperature superconducting ${\mathrm{H}}_{3}\mathrm{S}$, Phys. Rev. B 93, 104526 (2016).
  * Stanzione  _et al._ [2020] D. Stanzione, J. West, R. T. Evans, T. Minyard, O. Ghattas, and D. K. Panda, Frontera: The Evolution of Leadership Computing at the National Science Foundation, in  _Practice and Experience in Advanced Research Computing_ (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, ACCESS: Advancing Innovation: NSF’s Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support, in  _Practice and Experience in Advanced Research Computing_, PEARC ’23 (Association for Computing Machinery, New York, NY, USA, 2023) p. 173–176.

Supplemental Material:
Stability-superconductivity map for compressed Na-intercalated graphite

FIG. S1: Relative Gibbs free energy at 10 GPa as a function of temperature for
select Na-C phases calculated with the optB86b-vdW functional. Each line shows
relative stability with respect to a combination of known (circled numbers) or
proposed (bracketed letters) phases. The dashed lines correspond to phase
combinations containing the diamond ground state at this pressure, which we
argue to be not relevant for Na intercalation of graphite via cold
compression. Since all solid lines measuring the relative stability with
respect to graphite are below zero (shown in green), all four proposed phases
are expected to be synthesizable in the considered temperature range. In
addition, oP104-Na3C10 is globally stable above about 200 K. FIG. S2: Relative
Gibbs free energy at 10 GPa as a function of temperature for select Na-C
phases calculated with the optB88-vdW functional (see Fig. S1 caption for
further details). All four phases appear synthesizable in the considered
temperature range if graphite is used as a starting material, with two of
them, oP104-Na3C10 and oS20-NaC4, being globally stable with respect to
diamond. FIG. S3: Stability of C-Na phases at 10 GPa and 0 K calculated with
the r2SCAN+rVV10 functional [44]. The global (local) convex hulls are denoted
with solid (dashed) lines. The set of four locally stable phases is consistent
with those obtained with the optB88-vdW and optB86b-vdW functionals (see Fig.
2 in the manuscript). FIG. S4: Electronic band structures, DOS, phonon
dispersion, and PhDOS calculated with Quantum ESPRESSO for (a) oS52
($Cmcm$)-Na3C10, (b) oS52 ($C222_{1}$)-Na3C10, (c) oS20-NaC4, (d) mP20-NaC4,
(e) oP18-NaC8, (f) oS44-NaC10, and (g) hP13-NaC12 at 10 GPa. FIG. S5: Phonon
dispersion for oP104-Na3C10 at 10 GPa calculated with the finite displacement
method using PHONOPY and VASP. The linear response calculations proved to be
particularly demanding in this case.

Listing 1: CIF file for mP12-NaC2 containing the VASP optimized lattice
constants and atomic positions at 5 GPa.

⬇

_symmetry_Int_Tables_number 11

_cell_length_a 4.9727911481359230

_cell_length_b 4.3033130471733356

_cell_length_c 7.2283164338578878

_cell_angle_alpha 90.0000000000000000

_cell_angle_beta 98.7971251687718848

_cell_angle_gamma 90.0000000000000000

_symmetry_Int_Tables_number 11

_chemical_formula_sum

’C Na ’

loop_

_atom_site_label

_atom_site_type_symbol

_atom_site_site_symmetry_multiplicity

_atom_site_occupancy wyckoff

_atom_site_fract_x

_atom_site_fract_y

_atom_site_fract_z

C1 C 4 f 0.3739409661740112 0.9170379731911351 0.9946890523514327

C2 C 4 f 0.8750958655346466 0.9166666666666666 0.0006343559857812

Na1 Na 2 e 0.5359695585021954 0.2500000000000000 0.6650636316635571

Na2 Na 2 e 0.0402706172980845 0.7500000000000000 0.6740356198602769

#End

Listing 2: CIF file for oS52-Na3C10 ($Cmcm$) containing the VASP optimized
lattice constants and atomic positions at 10 GPa.

⬇

_symmetry_Int_Tables_number 63

_cell_length_a 4.3345876014862990

_cell_length_b 12.4744252875898436

_cell_length_c 8.6445697968283355

_cell_angle_alpha 90.0000000000000000

_cell_angle_beta 90.0000000000000000

_cell_angle_gamma 90.0000000000000000

_symmetry_Int_Tables_number 63

_chemical_formula_sum

’C Na ’

loop_

_atom_site_label

_atom_site_type_symbol

_atom_site_site_symmetry_multiplicity

_atom_site_occupancy wyckoff

_atom_site_fract_x

_atom_site_fract_y

_atom_site_fract_z

C1 C 16 h 0.8333333333333334 0.2001572180448292 0.9922351866168100

C2 C 16 h 0.3333333333333333 0.1000940496660334 0.5120353617970680

C3 C 8 e 0.1669435501513448 0.0000000000000000 0.5000000000000000

Na1 Na 4 c 0.5000000000000000 0.2484765693898439 0.7500000000000000

Na2 Na 4 c 0.5000000000000000 0.4049998440550624 0.2500000000000000

Na3 Na 4 c 0.0000000000000000 0.4420250481179195 0.7500000000000000

#End

Listing 3: CIF file for oS52-Na3C10 ($C222_{1}$) containing the VASP optimized
lattice constants and atomic positions at 10 GPa.

⬇

_symmetry_Int_Tables_number 20

_cell_length_a 4.3283454520037719

_cell_length_b 12.4550961726904426

_cell_length_c 8.7019554570706106

_cell_angle_alpha 90.0000000000000000

_cell_angle_beta 90.0000000000000000

_cell_angle_gamma 90.0000000000000000

_symmetry_Int_Tables_number 20

_chemical_formula_sum

’C Na ’

loop_

_atom_site_label

_atom_site_type_symbol

_atom_site_site_symmetry_multiplicity

_atom_site_occupancy wyckoff

_atom_site_fract_x

_atom_site_fract_y

_atom_site_fract_z

C1 C 8 c 0.8565797045459753 0.2003656544749290 0.9906896905986997

C2 C 8 c 0.3098836960032756 0.1001369377800074 0.5094431237821340

C3 C 4 a 0.1436128152677867 0.0000000000000000 0.5000000000000000

C4 C 8 c 0.6902220277336751 0.3001567132545233 0.0045171073435595

C5 C 8 c 0.1433277497844330 0.4001932826437908 0.4831244691225791

C6 C 4 a 0.8096430435082667 0.0000000000000000 0.5000000000000000

Na1 Na 4 b 0.5000000000000000 0.0695855689837757 0.2500000000000000

Na2 Na 4 b 0.5000000000000000 0.2418185439696074 0.7500000000000000

Na3 Na 4 b 0.5000000000000000 0.4098049512768895 0.2500000000000000

#End

Listing 4: CIF file for oP104-Na3C10 containing the VASP optimized lattice
constants and atomic positions at 10 GPa.

⬇

_symmetry_Int_Tables_number 19

_cell_length_a 8.6571339789676092

_cell_length_b 8.7003343806980844

_cell_length_c 12.4566915428260288

_cell_angle_alpha 90.0000000000000000

_cell_angle_beta 90.0000000000000000

_cell_angle_gamma 90.0000000000000000

_symmetry_Int_Tables_number 19

_chemical_formula_sum

’C Na ’

loop_

_atom_site_label

_atom_site_type_symbol

_atom_site_site_symmetry_multiplicity

_atom_site_occupancy wyckoff

_atom_site_fract_x

_atom_site_fract_y

_atom_site_fract_z

C1 C 4 a 0.4507620388737408 0.2412722355502409 0.2002836109717944

C2 C 4 a 0.7159835056491347 0.7604690269915223 0.1001229603435971

C3 C 4 a 0.7990938581497673 0.7500000000000000 0.0000000000000000

C4 C 4 a 0.5339395625350393 0.2553336754710428 0.3000932659378034

C5 C 4 a 0.7993055002287193 0.7339035928609262 0.4001689503665201

C6 C 4 a 0.4660637764204765 0.7500000000000000 0.0000000000000000

C7 C 4 a 0.7992359572009847 0.7413680362824165 0.7997428117503628

C8 C 4 a 0.5339920062649526 0.2607247608217022 0.8999032511803474

C9 C 4 a 0.4508054785026872 0.2500000000000000 0.9999865914032778

C10 C 4 a 0.7160632668574110 0.7552653191550148 0.6999122605521108

C11 C 4 a 0.4507371346340321 0.2340265058040031 0.5998443164547779

C12 C 4 a 0.7838311343020977 0.2497670407480256 0.0000000000000000

C13 C 4 a 0.4507684383074065 0.2588296507732223 0.7997299759256646

C14 C 4 a 0.7159641189005654 0.7392555506996387 0.8999278505968743

C15 C 4 a 0.5339370177679471 0.2449503259371396 0.6999032756511808

C16 C 4 a 0.7993444489652789 0.7657196483495666 0.5998772927784869

C17 C 4 a 0.7992555146565010 0.7587831215649031 0.2003039368210861

C18 C 4 a 0.5339820054295076 0.2395603206925400 0.1000893926225570

C19 C 4 a 0.7160624746705367 0.7449077831784620 0.3000949858428325

C20 C 4 a 0.4507027849798250 0.2657543547107649 0.4001480943255583

Na1 Na 4 a 0.6252816249467941 0.0000000000000000 0.9281883876001086

Na2 Na 4 a 0.8744455038088575 0.5000000000000000 0.2592461429000107

Na3 Na 4 a 0.8754102156613045 0.9999773581896462 0.7412876496637536

Na4 Na 4 a 0.6245455364933862 0.5000000000000000 0.4102936122942250

Na5 Na 4 a 0.6256409368228546 0.0000000000000000 0.5901517763312539

Na6 Na 4 a 0.8750732446945194 0.5000000000000000 0.5696379329270868

#End

Listing 5: CIF file for oS20-NaC4 containing the VASP optimized lattice
constants and atomic positions at 10 GPa.

⬇

_symmetry_Int_Tables_number 63

_cell_length_a 8.4461896681214643

_cell_length_b 4.9787180730519980

_cell_length_c 4.3158958682925386

_cell_angle_alpha 90.0000000000000000

_cell_angle_beta 90.0000000000000000

_cell_angle_gamma 90.0000000000000000

_symmetry_Int_Tables_number 63

_chemical_formula_sum

’C Na ’

loop_

_atom_site_label

_atom_site_type_symbol

_atom_site_site_symmetry_multiplicity

_atom_site_occupancy wyckoff

_atom_site_fract_x

_atom_site_fract_y

_atom_site_fract_z

C1 C 16 h 0.2435236207386140 0.6250134563772747 0.0829306596834911

Na1 Na 4 c 0.0000000000000000 0.2884085523109192 0.2500000000000000

#End

Listing 6: CIF file for mP20-NaC4 containing the VASP optimized lattice
constants and atomic positions at 10 GPa.

⬇

_symmetry_Int_Tables_number 11

_cell_length_a 4.9763713072779092

_cell_length_b 4.2971403813973037

_cell_length_c 8.6995515646003394

_cell_angle_alpha 90.0000000000000000

_cell_angle_beta 99.7172128079009497

_cell_angle_gamma 90.0000000000000000

_symmetry_Int_Tables_number 11

_chemical_formula_sum

’C Na ’

loop_

_atom_site_label

_atom_site_type_symbol

_atom_site_site_symmetry_multiplicity

_atom_site_occupancy wyckoff

_atom_site_fract_x

_atom_site_fract_y

_atom_site_fract_z

C1 C 4 f 0.3748660758264677 0.9173775863282644 0.5006096736221578

C2 C 4 f 0.8723223735758482 0.9179787848993116 0.4887063160384884

C3 C 4 f 0.6277710174186169 0.4180513681935950 0.0114755861022431

C4 C 4 f 0.1251066518153792 0.4174037974332522 0.9995767940923272

Na1 Na 2 e 0.4840379642283494 0.2500000000000000 0.7490266449917584

Na2 Na 2 e 0.0180017815718756 0.7500000000000000 0.7504474160489637

#End

Listing 7: CIF file for oP18-NaC8 containing the VASP optimized lattice
constants and atomic positions at 10 GPa.

⬇

_symmetry_Int_Tables_number 51

_cell_length_a 4.2849416211107201

_cell_length_b 7.0408405569361925

_cell_length_c 4.9454256864695196

_cell_angle_alpha 90.0000000000000000

_cell_angle_beta 90.0000000000000000

_cell_angle_gamma 90.0000000000000000

_symmetry_Int_Tables_number 51

_chemical_formula_sum

’C Na ’

loop_

_atom_site_label

_atom_site_type_symbol

_atom_site_site_symmetry_multiplicity

_atom_site_occupancy wyckoff

_atom_site_fract_x

_atom_site_fract_y

_atom_site_fract_z

C1 C 8 l 0.5833333333333334 0.8017078173988708 0.1243533376894139

C2 C 8 l 0.9166666666666666 0.7962083500714385 0.6253601678460757

Na1 Na 2 f 0.2500000000000000 0.5000000000000000 0.6984680181220927

#End

Listing 8: CIF file for oS44-NaC10 containing the VASP optimized lattice
constants and atomic positions at 10 GPa.

⬇

_symmetry_Int_Tables_number 65

_cell_length_a 4.2701616893271561

_cell_length_b 12.3372454224847932

_cell_length_c 7.0227230974644392

_cell_angle_alpha 90.0000000000000000

_cell_angle_beta 90.0000000000000000

_cell_angle_gamma 90.0000000000000000

_symmetry_Int_Tables_number 65

_chemical_formula_sum

’C Na ’

loop_

_atom_site_label

_atom_site_type_symbol

_atom_site_site_symmetry_multiplicity

_atom_site_occupancy wyckoff

_atom_site_fract_x

_atom_site_fract_y

_atom_site_fract_z

C1 C 16 r 0.8333333333333334 0.2998818841363473 0.2986779735933709

C2 C 16 r 0.6660824722040530 0.3998340922640508 0.2940550137969069

C3 C 8 o 0.6672155408982263 0.0000000000000000 0.7077668535232888

Na1 Na 4 i 0.0000000000000000 0.3807971938647299 0.0000000000000000

#End

Listing 9: CIF file for hP13-NaC12 containing the VASP optimized lattice
constants and atomic positions at 5 GPa.

⬇

_symmetry_Int_Tables_number 191

_cell_length_a 4.2827286628374583

_cell_length_b 4.2827286628374583

_cell_length_c 7.2725399518830214

_cell_angle_alpha 90.0000000000000000

_cell_angle_beta 90.0000000000000000

_cell_angle_gamma 120.0000000000000142

_symmetry_Int_Tables_number 191

_chemical_formula_sum

’C Na ’

loop_

_atom_site_label

_atom_site_type_symbol

_atom_site_site_symmetry_multiplicity

_atom_site_occupancy wyckoff

_atom_site_fract_x

_atom_site_fract_y

_atom_site_fract_z

C1 C 12 n 0.6666666666666666 0.0000000000000000 0.2063351542498993

Na1 Na 1 b 0.0000000000000000 0.0000000000000000 0.5000000000000000

#End