High-temperature phonon-mediated superconductivity in monolayer Mg2B4C2

A two-dimensional material – Mg2B4C2, belonging to the family of the conventional superconductor MgB2, is theoretically predicted to exhibit superconductivity with critical temperature Tc estimated in the 47–48 K range (predicted using the McMillian-Allen-Dynes formula) without any tuning of external parameters such as doping, strain, or substrate-induced effects. The origin of such a high intrinsic Tc is ascribed to the presence of strong electron-phonon coupling and large density of states at the Fermi level. This system is obtained after replacing the chemically active boron-boron surface layers in a MgB2 slab by chemically inactive boron-carbon layers. Hence, the surfaces of this material are inert. Our calculations confirm the stability of 2D Mg2B4C2. We also find that the key features of this material remain essentially unchanged when its thickness is increased by modestly increasing the number of inner MgB2 layers.

Among all the Bardeen-Cooper-Schrieffer (BCS) type conventional superconductors, MgB 2 stands out with a record T c of 39 K, the highest reported T c at zero-pressure [34][35][36] . Such a high-T c in MgB 2 stems from the strong electron-phonon (el-ph) coupling occurring primarily due to the in-plane stretching of B-B bonds (i.e., E 2g phonon modes), which strongly couple with the selfdoped charge carriers from magnesium to boron atoms 26,[35][36][37][38] . Remarkably, only two (E 2g ) out of a total of nine phonon modes contribute strongly to the total el-ph coupling in MgB 2 [35][36][37][38][39][40][41][42][43][44] . Once the fundamental mechanism of such a high-T c in bulk MgB 2 was understood, which by the way was a subject of intense research for over a decade period 26,[35][36][37][38][39][40][41][42][43][44][45][46][47][48][49] , researchers started proposing ways to augment T c through rational materials design approach 25,37,43,44,[50][51][52] . Pickett and co-workers proposed that one can, in principle, achieve a much higher T c (than 39 K) by designing a MgB 2 -like stable material that has a similar Fermi surface as in MgB 2 , and in which more than two phonon modes couple to the electronic states near the Fermi level, thereby, resulting in a sizable total el-ph coupling 25,43,44 . This idea has been employed for the rational design of bulk superconductors with a good success rate [53][54][55][56][57][58][59][60][61][62][63][64] . The high-pressure superconductivity observed at 250 K in lanthanum hydride is one such example 65-68 . Despite the large success with the bulk conventional superconductors, two-dimensional intrinsic superconductors having a high-T c remained elusive. Notably, various attempts have been made to realize superconductivity in the 2D analogues of bulk MgB 2 13,32,51,69-74 . On the one hand, Xu and Beckman proposed a quasi-2D MgB 2 nanosheet with inert surfaces, which turns out to be a semiconductor with a bandgap of 0.51 eV resulting from the quantum confinement effects 13 . On the other hand, Bekaert et al. reported that a considerably high-T c of 20 K can be realized in monolayer MgB 2 without surface passivation, i.e., if only such a material with a highly chemically reactive surface could be made 73,74 . In a recent study, Bekaert et al. theoretically demonstrated that a MgB 2 monolayer can be stabilized by adding hydrogen adatoms. Interestingly, they find that the hydrogenation process leads to a high-T c of 67 K, which can be further boosted to over 100 K by means of a biaxial strain on the hydrogenated MgB 2 monolayer 32 . While an experimental validation of the predicted T c in monolayer MgB 2 is still missing, the aforementioned theoretical works markedly enhance our understanding of superconductivity in 2D materials.
In this work, we present a MgB 2 -like 2D material -Mg 2 B 4 C 2 , having charge neutral inert surfaces, which is predicted to superconduct at a strikingly high-T c in the 47-48 K range (predicted using the McMillian-Allen-Dynes theory [75][76][77], which is among the highest T c yet reported for an intrinsic 2D material without any doping, strain or substrate-induced effects. The main advantageous feature in 2D Mg 2 B 4 C 2 is the fact that, unlike in bulk MgB 2 , more than two phonon modes strongly couple to the electronic states near the Fermi level, thus, resulting in a substantially larger el-ph coupling (λ = 1.40) in monolayer Mg 2 B 4 C 2 compared with the bulk MgB 2 (λ bulk = 0.73 38 , and 0.61 36 ). We note that the estimated λ in monolayer Mg 2 B 4 C 2 is comparable to the predicted λ = 1.46) in hydrogenated MgB 2 monolayer 32 . Moreover, our calculations reveal nontrivial topological electronic features in Mg 2 B 4 C 2 exhibiting Dirac cones and practically gapless Dirac nodal lines at the Fermi level near the corner points of the hexagonal Brillouin zone (BZ), which enhance the density of states (DOS) at the Fermi level by almost 30% compared to that of bulk MgB 2 , hence, positively contributing towards a higher T c .

RESULTS AND DISCUSSION Material design strategy
We start by describing our rationale for design of a stable MgB 2like 2D superconductor having inert surfaces. Generally, layered vdW materials can be exfoliated to produce their 2D analogues 78 . Although bulk MgB 2 has a layered structure, it is not a vdW material. Bulk MgB 2 crystallizes in space group P 6 /mmm (#191) containing alternating layers of Mg and B atoms stacked along the c ! lattice direction, as shown in Fig. 1(a) 42 . The bonding between the Mg and B atoms is purely ionic, which means that Mg atoms donate two electrons to B atoms, thereby making each Mg 2 + and each B 1 − . Since a B − is isoelectronic to a charge-neutral carbon atom, a B-B sheet is structurally analogous to a single layer graphene, but it has a different ordering of bands than those of graphene. A simple exfoliation of MgB 2 into a 2D slab with B (or Mg) termination would yield a highly reactive electron-rich (or hole-rich) surface layer that is chemically unstable.
We propose that one can passivate the charged surface layers in the MgB 2 slab by systematically substituting one boron by one carbon atom at the top and bottom surfaces of the slab. Figure 1 (b) shows the top and side views of a Mg 2 B 4 C 2 monolayer designed using the aforementioned strategy. Strikingly, we find that modestly repeating the intermediate Mg-B layers, i. e. , the layers sandwiched between the top and bottom surfaces (highlighted using dashed rectangle in Fig. 1(b)), thereby making thicker slabs of (MgB 2 ) n C 2 while remaining in the quasi-2D limit, n being the total number of Mg layers, retains the key features of the Mg 2 B 4 C 2 monolayer. The electronic bandstructures calculated up to n = 5 are shown in the Supplementary Fig. 1. This feature could be particularly useful in the experimental realization of 2D superconductivity in Mg 2 B 4 C 2 . We note that MgB 2 monolayer can also be passivated by an appropriate hydrogenation process 32 .
Mg 2 B 4 C 2 monolayer, shown in Fig. 1(b), belongs to the layer group p3m1 (#72) having DFT (PBE) optimized lattice parameters a = b = 2.87 Å. The absolute thickness between the top and bottom atomic layers is 7.14 Å, whereas, the interlayer spacing between the adjacent Mg and B-B, and Mg and B-C (C-B) layers is 1.8 Å, and~1.7 Å, respectively. We note that the inversion symmetry is preserved due to the inverted ordering of the top and bottom layers in the structure shown in Fig. 1(b). However, one could break the inversion symmetry by replicating the top and bottom layers, i.e., by making the top and bottom layers alike, either both as B-C or both as C-B. Our calculations suggest that the structure with inversion symmetry is energetically more favorable (5 meV/f.u.) than the structure with broken inversion symmetry; although both structures are dynamically, elastically, and mechanically stable since they exhibit all positive phonon frequencies, positive elastic constants, and satisfy the Born-Huang mechanical stability criteria (see Supplementary Table 2 and Supplementary Fig. 6). The only qualitative difference in the electronic properties of the structure with broken inversion symmetry is a small lifting of some band degeneracies at the K high symmetry point (see Supplementary Fig. 3). This effect is analogous to the application of a perpendicular electric field to a bilayer graphene 79 .
In this article, hereafter, we focus only on ground state structure of a monolayer Mg 2 B 4 C 2 with preserved inversion symmetry. We note, all other possible atomic configurations of this composition are higher in energy. Furthermore, our exfoliation energy calculations (see Supplementary Table 3) suggest that the reported monolayer Mg 2 B 4 C 2 belongs to the "easily exfoliable" category, as classified by Mounet et al. 80 .
Topological electronic properties of Mg 2 B 4 C 2 monolayer After describing the crystal structure and its stability, we now focus on the topological electronic properties of Mg 2 B 4 C 2 monolayer. We begin by summarizing the key features of the electronic structure of bulk MgB 2 42 from which Mg 2 B 4 C 2 monolayer is derived. As shown in Fig. 2(a), the Fermi surface of MgB 2 is composed of boron p orbitals, where p x,y orbitals hybridize with s orbitals to form strong covalent in-plane σ bonds at the zone center, while the unhybridized p z orbitals form relatively weak outof-plane π bonds at zone boundaries (Mg acts as electron donor). Due to such a distinct Fermi-surface geometry, two superconducting gaps exists in bulk MgB 2 : (i) the stronger σ gap of 7 meV, and (ii) the weaker π gap of~2-3 meV 35,36,39,46,[81][82][83][84][85] . Different symmetries of the σ and π bonds largely suppress the impurity scattering in MgB 2 39,41,42,85 . Since the basic structure and charge neutrality of MgB 2 is preserved in monolayer Mg 2 B 4 C 2 , the electronic band structure of monolayer Mg 2 B 4 C 2 qualitatively resembles with that of the bulk MgB 2 , as shown in Fig. 2(a, b), but with some additional features. For instance, there is a set of degenerate σ bands (σ outer ) present at Γ below the Fermi level arising from the p x,y orbitals of the outer boron-carbon layers. The other set of degenerate σ bands (σ inner ) at Γ that cross the Fermi level (also present in MgB 2 ) are formed by the p x,y orbitals of the inner boron-boron layer. These two sets of σ bands are almost parallel and split by~1.6 eV at Γ. Since the σ outer bands are completely occupied, they should, in principle, have no contribution in superconductivity, unless there is a large external field applied in a FET-like geometry 86 .
In addition to the set of σ bands at Γ, we notice the presence of Dirac-like band crossings at the K point, as well as along the highsymmetry directions near the K point of monolayer Mg 2 B 4 C 2 . Regardless of their topological nature, these band crossings at the Fermi level, highlighted using a dashed magenta box in Fig. 2(b), yield a large DOS at the Fermi level (almost 30% larger than in bulk MgB 2 ), which contributes substantially to the total el-ph coupling in the studied monolayer. We note that the Dirac-like crossing at K is also present in bulk MgB 2 , but it is situated wellabove the Fermi level 87 . The Dirac-like band crossings in Mg 2 B 4 C 2 monolayer are formed by highly dispersing p z orbitals of carbon and boron atoms (see Supplementary Fig. 4). Thus, the Fermi surface of Mg 2 B 4 C 2 monolayer, shown in Fig. 2(c), embodies three main features: (i) two-hole pockets at Γ (one circular and another that takes the shape of the BZ) composed of σ bonded boron p x,y orbitals, (ii) an electron pocket at M formed by boron p z orbitals, and (iii) intertwined electron and hole pockets at the K point and along K-M high-symmetry line, formed by π bonded carbon and boron p z orbitals. We note that all these pockets show very strong coupling to the phonon modes, and, as a result, they play the key role in governing superconductivity in Mg 2 B 4 C 2 monolayer, as we discuss later. Furthermore, the sharp and well-defined (almost flat) boundaries of the charge-carrier pockets at the Fermi surface set up the stage for the possible realization of Kohn-like divergencies 88 , and charge-density-wave ordering 89,90 in this 2D system, which is beyond the scope of present work and calls for a more comprehensive attention in the future.
By plotting the energy bandgap (E gap ) distribution in the vicinity of the K points, we discover presence of a triangular nodal line in the vicinity of each K point, as shown in Fig. 2(d). However, this is not a truly gapless nodal line since a small E gap (~5 meV) exists due to the subtle breaking of M z mirror symmetry. It is worth noting that the Dirac point at K is protected by the C 3v rotation, inversion, and time-reversal symmetries; a small gap opens at Dirac points when the inversion symmetry is broken by making the top and bottom B-C layers identical 91 . Although there are theoretical proposals suggesting the possibility of topological superconductivity in Dirac semimetals 92 , we think that the so-far studied models are quite simple, and this topic requires a more thorough examination before any exotic effects can be confidently claimed here.
In order to prove the nontrivial topological nature of Dirac points, we compute the Berry phase along a k-loop enclosing the gapless point at K, as marked using dashed lines in Fig. 2(d). Our calculations yield a nontrivial Berry phase of π for Dirac points at K. We note, this exercise could not be performed for the Dirac nodal line near K because enclosing the nodal line residing in the k x -k y plane would require a k-loop encircling along k z and k z is not defined for a 2D system. Nevertheless, the presence of time-reversal and spatialinversion symmetries of Mg 2 B 4 C 2 monolayer enables us to determine the Z 2 topological invariants using the Fu-Kane criterion 93 . The inversion parity eigenvalues of the electronic wavefunction of all 12 occupied bands at four time-reversal invariant momenta (TRIM) points are given in Table 1. The product of all parity eigenvalues (δ) at each TRIM is also listed in Table 1. We find that the Z 2 topological index is nontrivial due to δ = −1 at three TRIM points. Here, we note that bulk MgB 2 has a weak Z 2 topological index (0; 001) due to the band-inversions occurring at the Γ and A (0, 0, 0.5) high-symmetry points of 3D hexagonal BZ 87 . Robust topological surface states have recently been experimentally observed in bulk MgB 2 94 .
Since the nontrivial topology in 2D systems is often manifested in the gapless 1D edge states, we further confirm the nontrivial topological features of monolayer Mg 2 B 4 C 2 by computing the local density of states at (100) and (010) edges of 60 unit cell thick nano-ribbons. Topologically nontrivial 1D edge states connecting band-crossing points were obtained at both (100) and (010) edges, as shown in Fig. 2(e), thus, proving the nontrivial topology of the Mg 2 B 4 C 2 monolayer.

Electron-phonon coupling and superconductivity in Mg 2 B 4 C 2
We find that the roots of superconductivity in Mg 2 B 4 C 2 monolayer are the same as in bulk MgB 2 35-44 . However, the main advantageous factor in Mg 2 B 4 C 2 is that, in addition to the doubly degenerate E 2g modes that govern superconductivity in MgB 2 , numerous other phonon modes strongly couple to the electronic states near the Fermi level yielding a much larger overall el-ph coupling, and thus, resulting in a considerably higher T c .
The calculated phonon spectrum of Mg 2 B 4 C 2 monolayer, shown in Fig. 3(a), contains a total of 24 phonon modes (8 atoms/cell) Fig. 2 Electronic structure. Atomic orbitals projected electronic band structure of (a) bulk MgB 2 , and (b) monolayer Mg 2 B 4 C 2 calculated without spin-orbit coupling (SOC) along the high symmetry direction of BZ. Cyan, red, and blue colors represent the contribution from the s, p x,y , and p z orbitals, respectively. See Supplementary Fig. 4 for more details. c Calculated Fermi surface of monolayer Mg 2 B 4 C 2 . Light pink/ green, and grey colors depict hole/electron, and intertwined electron-hole pockets, respectively. d Energy bandgap (E gap ) plotted in color scale (eV units) in the vicinity of a K high-symmetry point. The dashed circle marks the k-loop along which Berry phase was computed. e The local electronic density of states of the (100) and (010) edge states spectrum. Red/White color denotes the states near the edge/interior of the 2D system. Topological nontrivial edge states are marked using arrows.
having the following mode symmetry at Γ: Here, A 1g and E g are Raman-active modes, whereas, A 2u and E u are infrared-active modes. In Fig. 3(c-f), we show the atomic vibration patterns for the four phonon modes, namely, three nondegenerate A 1g modes (indices 14, 18, and 19) and one degenerate E g mode (indices [16][17], which exhibit the dominant el-ph coupling. All these A 1g modes correspond to the out-of-plane vibrations of the Mg, inner B-B, and outer B-C layers, while the E g mode corresponds to the in-plane stretching of the inner B-B layer. The A 1g modes primarily modulate the el-ph coupling associated with the π bonded p z orbitals contributing to the electron and hole pockets located at the BZ boundaries. Whereas, the doubly degenerate E g mode couples with the σ bonded p x,y orbitals forming the hole pockets located at Γ. Here, it is worth noting that the higher frequency E g modes (indices [21][22] that correspond to the in-plane stretching of the outer B-C layers do not make a significant contribution to the overall el-ph in this system, which is as expected since these modes modulate the occupied σ outer bands located well-below the Fermi level at Γ [see Fig. 2(b)]. However, these modes may participate in superconductivity when the system is doped with p-type charge carriers 32 .
Since the electronic and vibrational band structures of inner B-B and outer B-C layers are essentially independent of each other, we predicate that the reported properties of the studied Mg 2 B 4 C 2 monolayer would be retained even when the number of the inner B-B layers are repeated (until a critical thickness), thus making the system thicker. This feature might greatly simplify the eventual realization of superconductivity in Mg 2 B 4 C 2 .
To quantify the superconducting properties of Mg 2 B 4 C 2 monolayer, we employ the McMillian-Allen-Dynes theory derived from the isotropic Migdal-Eliashberg formalism 75-77 which relies on the calculation of the el-ph coupling matrix elements within DFT. The calculated matrix elements correspond to the transition probabilities of different Kohn-Sham states induced by a change in the potential due to a small ionic displacement. Thus, these matrix elements provide the main ingredients to calculate the elph coupling strength and the Eliashberg spectral function α 2 F(ω) as a function of the phonon frequency ω. Since the physical process behind the phonon-mediated superconductivity is the exchange of a phonon between two electrons, a strong el-ph coupling is desired to achieve a high-T c in a BCS superconductor. Theoretical details of such calculations are explained in numerous other papers 82,95,96 .
In Fig. 3(a), we plot the calculated phonon linewidth λ(q, n) for each phonon mode n at each wave vector q using blue color. Note that the plotted phonon linewidth is scaled down by a factor of two to avoid large overlap with the neighboring phonon branches. The largest contribution to the total el-ph coupling strength comes from three nondegenerate A 1g modes and one doubly degenerate E g mode, as marked in Fig. 3(a). We note that the A 2u mode (index 15), marked using`×' in Fig. 3(a), does not contribute to the total el-ph coupling, although it appears buried in the large λ(q, n) overlap from the E g mode. Notably, in addition to the aforementioned A 1g and E g phonon modes, various other modes make relatively smaller contributions to the overall el-ph coupling strength, as revealed by the Eliashberg spectral function α 2 F(ω) plot shown in the right panel of Fig. 3(a).
In addition to the el-ph coupling, the net phonon linewidth λ(q, n) can have some contribution from the phonon-phonon (ph-ph) interactions owing to the phonon anharmonicity 36 . Therefore, we thoroughly investigate ph-ph interactions by computing ph-ph linewidth using the ab-initio molecular dynamics simulations. In this approach, we mapped the forces, obtained from the Fig. 3 Electron-phonon coupling. a Calculated phonon spectrum of Mg 2 B 4 C 2 monolayer with phonon linewidths λ(q, n) plotted using shaded blue color. To avoid large overlap of λ(q, n) with the phonon spectra, we have divided the intensity by a factor of two. The colored circles mark the three out-of-plane nondegenerate A 1g modes (indices 14, 18, and 19), and the yellow diamond marks one in-plane doubly degenerate E g mode at Γ. These modes exhibit dominant el-ph coupling. The atomic displacement patterns corresponding to these modes are shown in (c-f). The nondegenerate A 1u mode (index 15) marked using symbol`× ' does not contribute to the total el-ph coupling, although it appears to be buried in the large λ(q, n) of the E g mode. The numerals 14,15,16,17,18, and 19 denote the phonon mode index as counted from the lowest to the highest frequency modes (i.e.,1-3 for acoustic modes). Mg atoms are omitted in (c) for the sake of clarity. The Eliashberg spectral function α 2 F(ω) along with the el-ph coupling constant λ in plotted in the right panel of (a). b Estimated T c as a function of the μ * parameter.
finite-temperature molecular dynamics simulations, evaluated in a 3 × 3 × 1 supercell onto a model Hamiltonian describing the lattice dynamics. This temperature dependent effective potential (TDEP) technique 97,98 enabled us to calculate the third-order response from the effective renormalized interatomic force constants. Our calculations revealed that the ph-ph linewidths are an order of magnitude smaller than the el-ph linewidths. The maximum value of obtained ph-ph linewidth is~2 meV, which is much smaller compared to the el-ph linewidth values that are typically larger than~70 meV in the studied system. This result implies that, although the system inherits some anharmonic effects, we can safely discard the ph-ph contributions in the study of its superconducting properties.
Based on the BCS theory of superconductivity and above results, we estimate the critical temperature T c using the McMillian-Allen-Dynes formula 99-101 : where ω log is the logarithmic averaged phonon frequency, λ is the total el-ph coupling constant, and μ * is the effective screened Coulomb repulsion constant with a typical value ranging from 0.04 to 0.16 (see Table 2) 21,35,36,38,56 . We obtain λ by integrating the cumulative frequency-dependent el-ph coupling λ(ω) given by the following expression: We find a fairly large value of λ = 1.40, which is considerably larger than the one reported for bulk MgB 2 (λ bulk = 0.73 38 , and 0.61 36 ). We observe that the estimated T c does not vary drastically as a function of μ * , as shown in Fig. 3(b). This is consistent with an earlier work by Choi et al. 36 , which reported that the superconducting properties of MgB 2 are not very sensitive to the μ * parameter within the isotropic McMillian-Allen-Dynes formalism. We note that for bulk MgB 2 , μ * = 0.05 has been used to get the correct estimate of T c~4 0 K 38 . Therefore, using the McMillian-Allen-Dynes formula 99-101 , we estimate the T c of Mg 2 B 4 C 2 monolayer to be in the range 47-48 K without any doping or strain. Our results are consistent with a recent study 32 in which T c = 67 K and λ = 1.46 was predicted in a hydrogenated MgB 2 monolayer by solving the fully anisotropic Eliasberg equations. We argue that the predicted T c in Mg 2 B 4 C 2 monolayer can be further enhanced by biaxial strain 17,32 or by p-doping 32 . In passing, we would like to mention that the predicted T c could moderately vary if a fully anisotropic Migdal-Eliashberg theory 32,36,82,102 or SC-DFT [103][104][105][106] is employed. This is particularly important here because the applicability of the McMillian-Allen-Dynes formula becomes limited in the case of large el-ph coupling. In order to highlight the novelty of our results, in Table 2 we list the theoretical superconducting parameters along with the estimated T c for some reported 2D phonon-mediated superconductors (see Supplementary Table 1). The good agreement between the experimental data for LiC 6 102,107,108 , 2H-NbSe 2 7, [109][110][111] , and C 6 CaC 6 112,113 and the corresponding theoretical results obtained from the McMillian-Allen-Dynes theory validate the the predictive power of the McMillian-Allen-Dynes theory.

SUMMARY
In summary, we present a 2D material Mg 2 B 4 C 2 , similar to MgB 2 , but with inert surfaces obtained by the replacement of outer B-B layers by B-C layers. Our calculations suggest that this structure is dynamically, elastically and mechanically stable. It also features a nontrivial topological electronic band structure together with a large el-ph coupling (λ = 1.40), which is more than twice as large as that of the bulk MgB 2 and comparable to that of in a hydrogenated monolayer MgB 2 32 . Use of the standard McMillian-Allen-Dynes theory predicts the superconducting transition temperature T c to be in the range of 47-48 K without any doping or tuning of external parameters such as strain. To the best of our knowledge, this is among the highest predicted intrinsic T c in a conventional BCS-type 2D superconductor. In addition to the large el-ph coupling, the presence of sharp and well-defined flat boundaries of the charge-carrier pockets at the Fermi surface imply the possible realization of Kohn-like divergencies and Table 2. Listing of superconducting parameters required for the prediction of T c using the McMillian-Allen-Dynes formula for some reported 2D phonon-mediated superconductors (data for bulk MgB 2 are included for comparison). This table includes data of effective Coulomb screening parameter μ * , electronic DOS at the Fermi level N(E F ) (in states/spin/Ry/cell), logarithmic averaged phonon frequency ω log (in K), total electronphonon coupling constant λ, and estimated T c (in K). Experimental T c values are noted in the table (see the Supplementary Table 1  charge-density-wave ordering in this 2D system, which calls for a dedicated study in future.

DFT calculations
The electronic bands structure and phonon calculations were performed using density-functional theory (DFT) as implemented in the VASP [114][115][116][117] . The phonopy 118 and PyProcar 119 tools were used for the post-processing of data. The Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional 120 and PAW pseudo-potentials 121,122 were used. The employed k-point grid for self-consistent calculations was 30 × 30 × 1, and the cutoff for the kinetic energy of plane waves was set to 700 eV. A vacuum of thickness~30 Å was added to avoid the periodic interactions along the c-axis. Since the spin-orbit coupling (SOC) effects were found to be negligible in the studied system, SOC was not included in the reported calculations. The elastic and mechanical properties were analyzed using the MechElastic code 123,124 . The exfoliation energy was calculated using four different exchange-correlation approximations: the (PBE) GGA approximation 120 , the SCAN 125 meta-GGA, vdW-DF2 GGA functional 126 , and SCAN together with the rVV10 correlation functional (SCAN+rVV10) 127 .
The topological properties of Mg 2 B 4 C 2 were studied by fitting the DFT calculated bandstructure to a real space tight-binding Hamiltonian obtained using the maximally localized Wannier functions (MLWFs) approach 128,129 . The local density of states at (100) and (010) edges were calculated for 60 unit cells thick nano-ribbons using the WannierTools package 129 with vacuum added along the c-axis of the ribbon. For the electron-phonon coupling matrix elements calculations, we used the abinit package [130][131][132][133][134] . We employed norm conserving pseudopotentials (using the ONCVPSP scheme of Hamann 135 ), and a plane wave basis set up to kinetic energies of 35 Ha. Cell parameters were optimized by using the PBE exchange-correlation functional as in VASP calculations. We used a uniform grid of 18 × 18 × 1 for the ground state calculations, and a phonon grid of 9 × 9 × 1 for the phonon part. A total of 288 el-ph matrix elements were calculated. Calculations of the phonon interatomic force constants, and the el-ph coupling matrix elements performed in this work used the second-order perturbation theory 136,137 . The temperature dependent effective potential (TDEP) technique 97,98 was used to study the phonon-phonon interactions and phonon anharmonic effects.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding authors upon reasonable request.

CODE AVAILABILITY
The first-principles DFT calculations were performed using the privately-licensed VASP code, and the ABINIT code, which is available at https://www.abinit.org under the GNU General Public License.