Defect-induced magnetism and Yu-Shiba-Rusinov states in twisted bilayer graphene

TL;DR

This study demonstrates how atomic defects, such as hydrogen adatoms and vacancies, induce local magnetic moments and Yu-Shiba-Rusinov states in twisted bilayer graphene, affecting its superconducting properties.

Contribution

It provides first-principles and tight-binding evidence that defects create magnetic moments and in-gap states in twisted bilayer graphene, revealing new defect-superconductivity interactions.

Findings

Hydrogen adatoms induce triple-point crossings.

Vacancies cause self-doping and symmetry breaking.

Defects lead to local magnetism and Yu-Shiba-Rusinov states.

Abstract

Atomic defects have a significant impact in the low-energy properties of graphene systems. By means of first-principles calculations and tight-binding models we provide evidence that chemical impurities modify both the normal and the superconducting states of twisted bilayer graphene. A single hydrogen atom attached to the bilayer surface yields a triple-point crossing, whereas self-doping and three-fold symmetry-breaking are created by a vacant site. Both types of defects lead to time-reversal symmetry-breaking and the creation of local magnetic moments. Hydrogen-induced magnetism is found to exist also at the doping levels where superconductivity appears in magic angle graphene superlattices. As a result, the coexistence of superconducting order and defect-induced magnetism yields in-gap Yu-Shiba-Rusinov excitations in magic angle twisted bilayer graphene.

Figures

Equations

$$H_0=\sum _{⟨ij⟩}t{c}_i^{†}{c}_j+\sum _{ij}{t}_{\perp }(r_i,r_j)c_i^{†}{c}_j$$

$$H_U=\sum _{i}U{c}_{i,\uparrow }^{†}{c}_{i,\uparrow }{c}_{i,\downarrow }^{†}{c}_{i,\downarrow }=\sum _{i}U{n}_{i,\uparrow }{n}_{i,\downarrow }$$

$$H_U^{MF}=\sum _{i}U[⟨n_{i,\uparrow }⟩n_{i,\downarrow }+⟨n_{i,\downarrow }⟩n_{i,\uparrow }-⟨n_{i,\uparrow }⟩⟨n_{i,\downarrow }⟩]$$

$$H=H_0+H_U^{MF}+H_{SC}$$

Full text

I Introduction

The tunability of grapheneCastro Neto et al. (2009) and two-dimensional materials has provided an outstanding solid-state platform to explore emergent physical phenomena. Yet, on top its standalone interest, graphene provides a powerful building block to create superstructures due to the weak van der Waals forces between layers.Geim and Grigorieva (2013) Those very same weak van der Waals forces allow to deposit graphene layers on top of each other with relative angles, in stark contrast with conventional bulk crystals. Twisted bilayer graphene (tBLG)Lopes dos Santos et al. (2007) is one the simplest structures that can be built in that fashion, and has opened the door towards realizing states of matter inaccessible in monolayer graphene. Rickhaus et al. (2018); Cao et al. (2018a, b); Sharpe et al. (2019); Cao et al. (2019); Jiang et al. (2019) This additional flexibility stems from the emergent moiré pattern that arises due to the twist between the layers, giving rise to an emergent electronic structure that can be controlled by the twist angle.Suárez Morell et al. (2010); Bistritzer and MacDonald (2011); San-Jose and Prada (2013); Ramires and Lado (2018)

A paradigmatic example of new physics associated with the emergent moiré band structure is the appearance of superconductivity in magic angle graphene ($\alpha\approx 1^{\circ}$) superlattices,Cao et al. (2018a) that exploit a chemical-free field effect doping of graphene without altering its structural integrity. Moreover, experimental observations reported strong correlated behavior,Cao et al. (2018b) anomalous Hall effect,Sharpe et al. (2019) strange metal behaviorCao et al. (2019), and rotational symmetry-breakingJiang et al. (2019) close to that regime. Extensive theoretical efforts are being directed towards deriving a faithful low-energy model,Kang and Vafek (2018); Koshino et al. (2018); Po et al. (2018); Guinea and Walet (2018) studying the possible electronic instabilities,Isobe et al. (2018); Xu and Balents (2018) and providing means of tailoring the electronic properties of such state.Chittari et al. (2018); Yankowitz et al. (2019a); Carr et al. (2018a) However, the influence of atomic defects in this structures has not been addressed in detail, and their potential impact in the low-energy properties has remained rather unexplored.Brihuega and Yndurain (2018); Ramires and Lado (2019)

In this paper, we study the impact of atomic defects in tBLG with angles ranging from large angles $\alpha\approx 5^{\circ}$ to the magic angle $\alpha\approx 1^{\circ}$. First-principles methods and effective real-space tight-binding models are used to describe two paradigmatic defects in graphene systems, namely carbon vacanciesUgeda et al. (2010); Yazyev (2010); Palacios et al. (2008); Pereira et al. (2008); Yazyev and Helm (2007) and chemisorbed H atomsBrihuega and Yndurain (2018); Gonzalez-Herrero et al. (2016); García-Martínez et al. (2017); Lopez-Bezanilla et al. (2009); Ramires and Lado (2019) as shown in Fig. 1. In particular, the chemical modification introduced by mono-hydrogenation creates flat bands at the Fermi energy, whose unbalanced electronic occupation yields local magnetic moments. In the superconducting state of magic-angle tBLG, such magnetic moments are shown to in-gap Yu-Shiba-Rusinov excitations, dramatically impacting the spectral properties of the superconducting state.

The manuscript is organized as follows. In Sec. II the emergence of magnetism as a result of both C vacancies and mono-hydrogenation is analyzed with first-principles calculations. Sec. III shows that chemisorbed H can be treated with a low-energy model that reproduces accurately the first-principles results. In Sec. IV we show that the interplay between impurity induced magnetism and superconductivity gives rise to in-gap Yu-Shiba-Rusinov states. Finally, Sec. V summarizes our conclusions.

II Defect-induced states from first principles

II.1 General considerations

We first study the effect of structural and chemical modification of tBLG with first-principle methods. Moiré lattices are created by twisting a graphene layer over the other by a commensurate angle, so that a unit cell with coincidence site lattice points is createdCampanera et al. (2007); Lopes dos Santos et al. (2007); Dai et al. (2016). In the chosen tBLG geometry, one layer of AA-bilayer graphene is rotated an angle $\theta=5.09^{\circ}$ with respect to the other layer. Internal coordinates of pristine and modified tBLGs separated initially by the van der Waals distance of 3.3Å were fully relaxed. The momentum mismatch between the layers yields two linearly dispersive bands, which are degenerate in the $\Gamma-K-M$ high-symmetry path of the BZ.Suárez Morell et al. (2010); Lopes dos Santos et al. (2007); Bistritzer and MacDonald (2011); Latil and Henrard (2006); Latil et al. (2007)

First-principles calculations provide a powerful insight on effect of structural relaxation on the electronic structure of defective systems, together with the properties of the remaining unpaired electron. Chemisorption of a single H atom Brihuega and Yndurain (2018); Gonzalez-Herrero et al. (2016); García-Martínez et al. (2017); Lopez-Bezanilla et al. (2009); Ramires and Lado (2019) is a simple way of creating an imbalance between the number of A and B sites in graphene, with bare modification of the bipartite nature of graphene latticeGargiulo et al. (2014); Katoch et al. (2018); Soriano et al. (2011); Gmitra et al. (2013); Ulybyshev and Katsnelson (2015); Kochan et al. (2014); Wong et al. (2018); Farjam et al. (2011); Matis et al. (2013); Vergés and de Andres (2010); Gargiulo et al. (2014); Lei et al. (2008). The imbalance between the number of sites in the two sublattices has an important consequence that follows from Lieb’s theoremLieb (1989): a bipartite lattice with more A sites than B sites will show a number of zero modes proportional to such difference. In a similar fashion, a vacant C site breaks the lattice symmetry by removing one atom, but creating addtional distortionsYazyev (2008); Zhang et al. (2016); Lehtinen et al. (2004); Sofo et al. (2012); Nakhmedov et al. (2019); Miranda et al. (2016); Wang and Pantelides (2012). In an oversimplified picture, both approaches are equivalent to removing a $\pi$ electron. However, the sp3 rehybridization of an atom or the creation of dangling $\sigma$-bond have completely different consequences on the localization an dispersion of the electronic state. This feature can be clearly seen in the charge redistribution of Fig. 2, where it is shown that the vacancy triggers an atomic reconstruction that breaks $C_{3}$ symmetry (Fig. 2a), whereas the hydrogen ad-atom creates a redistribution that conserves the original three-fold rotation (Fig. 2b).

Due to the varying distance between C atoms across the parallel layers, the effect of a H atom depends on its location.Mapasha et al. (2012); Katoch et al. (2018); Yuan et al. (2010); Ferreira et al. (2011); Dahal et al. (2008) Similarly, a vacant site interacts differently with the opposite graphene layer depending on the stacking region where the C atom was removed from.Palacios and Ynduráin (2012); Andelkovic et al. (2018); Ulman and Narasimhan (2014) In the following, we address in detail the electronic reconstruction associated to each one of the defects introduced.

II.2 Carbon vacancy

Within the DFT formalism, a vacancy is more than a simple removal of a lattice site, and reconstruction of the whole tBLG structure is consideredUgeda et al. (2010). Atomic rearrangement is present in real systems, and thus removing a C atom out of one of the graphene sheet has a non-trivial effect on the distribution of electronic states in the low-energy band diagram. Fig. 3a displays a paramagnetic calculation of a fully relaxed tBLG structure with a vacant site. The electronic state at the Fermi level corresponds to the defect state, which removes the degeneracy of the bands that form the Dirac cones and creates an almost degenerate (4 meV gaped) triple-crossing at the $K$ point. The vacancy state creates a major disruption on one of the sets of linear bands, pushing them at higher and lower energies.

By introducing the spin degree of freedom, a magnetic ground state with an associated spin splitting of the defect band is found (Fig.3b). A magnetic moment of 1.5 $\mu_{B}$ is created as a result of the p-type doping character that unbalances the $\pi$-bands filling at the Fermi level. A vacancy in a tBLG polarizes the bilayer by distorting the dispersion of the $\pi$-bands and introducing a p-doping on the opposite layer that unbalances the electronic population of the valence bands. Furthermore, the defect creates two spin-dependent dispersionless bands at very different energies, (-0.5 and 1.35 eV for spin-$\uparrow$ and spin-$\downarrow$ respectively) as show in Fig.3b, and whose spatial extension is reduced to the immediate surrounding of the vacancy. Real-space representation of the charge density distribution in Fig.2a shows that the unpaired electron is located in the $\sigma$-bonds of the C atoms at the vacant site. The highly localized nature of the defect, together with the location of its electronic states in the band diagram is expected to yield a greatly enhanced signal in a scanning tunnel microscope experiment.Brihuega and Yndurain (2018); Ugeda et al. (2010)

Removing one electron out of the system barely changes the shape of the electronic bands, but empties completely the spin-$\uparrow$ band, leaving the other partially filled (see Fig.3c). A total magnetic moment of 0.9 $\mu_{B}$ is obtained, similarly to the case of monolayer grapheneUgeda et al. (2010). Removing an additional C atom near to the removed site creates a divacancy, that restores the non-magnetic state of the structure by effectively removing the unpaired electrons of the single-vacancy. As shown in Fig. 3d, a defect state similar to the paramagnetic single vacancy state is located in the vicinity of the Fermi level and the flat bands at low and high energies disappear.

II.3 Chemisorbed hydrogen atom

We now move on to consider the effect of a single H atom attached to one of the graphene layers. Relaxation of the entire structure results on the adsorbent C atom site pulled out of the graphene plane, adopting the pyramidal geometry characteristic of the sp3 hybridization. Moreover, the bond lengths between the anchoring carbon atom and its first neighboring atoms were elongated to 1.49 Å, and the C-H bonding distance is of 1.14 Å. Interestingly, the atomic reconstruction preserves the local $C_{3}$ symmetry (Fig. 2b), in strike comparison with the vacancy case. The atomic positions in the opposite layer are barely affected, and the band structure of the system shows a triple point crossing at the K point (Fig. 4a). In a spin polarized calculation, a spin splitting appears in the band structure, breaking time reversal symmetry as a consequence of the net magnetic moment triggered by the presence of the hydrogen zero mode (Fig. 4bc). In particular, the exchange field generates a shift of $\sim$0.25 meV up and down of the impurity band. Interestingly, for any attachment position, the linear dispersion typical of graphene is preserved, although the Dirac cones are shifted in energy and a meV gap removes the Dirac point.

For an H atom sitting in the AB stacking region, one of the electronic bands associated with the defect is almost fully populated with one electron, whereas a similar band with opposite spin quantum number is almost empty. The total magnetic moment is of 0.8 $\mu_{B}$. A similar behavior is observed for an H atom sitting on a AA stacked C atom, although the resulting magnetic moment decreases to 0.5 $\mu_{B}$. A small n-type doping is induced in the layer opposite to the one with the H atom. This last feature can be inferred from the band diagrams of Fig.4b and c, where the upper Dirac cone becomes slightly filled. This electron donor character demonstrates that both layers are partially hybridized through the defect state and a partial charge transfer occurs between layers. Finally, we consider the effect of a doping of a single electron per moiré supercell in this regime $\alpha\approx 5^{\circ}$. In that situation, we observe that the magnetic moment of the impurity state is completely quenched (Fig.4d), yielding a fully occupied impurity band. As it will shown in Sec. IV, the phenomenology at $\alpha\approx 1^{\circ}$ is dramatically different due to the appearance of the magic angle flat bands.

The previous phenomenology applies to twisted bilayers with large twisting angle $\alpha\approx 5^{\circ}$. In the following we will focus in smaller angles tBLG, that exhibit a substantially different behavior but for which DFT-based calculations are computationally demanding due to the large number of atoms in the unit cell. To circumvent such a limitation, the same phenomenology will be explored within a real-space tight-binding model, that allows to ultimately consider the effects of interactions in mono-hydrogenated magic angle superlattices with a more affordable computational cost and similar accuracy.

III Defect-induced states within a tight-binding model

We now move on to consider the effect of a single-H atom attached to the tBLG surface with an effective tight-binding model of the form

[TABLE]

where $\langle ij\rangle$ denotes sum over first neighbors in the same layer. The interlayer hopping $t_{\perp}(\vec{r})$ is modulated according to the structure of tBLG, which is taken as $t_{\perp}(\mathbb{r}_{i},\mathbb{r}_{j})=t_{\perp}\frac{(z_{i}-z_{j})^{2}}{|\mathbb{r}_{i}-\mathbb{r}_{j}|^{2}}e^{-\beta(|\mathbb{r}_{i}-\mathbb{r}_{j}|-d)},$ with $\mathbb{r}_{i}$ the position of site $i$ and $d$ the distance between layers. As a reference, the hopping parameters in graphene are $t\approx 3$ eV and $t_{\perp}\approx 0.4$ eV.Sboychakov et al. (2015) For computational convenience, in the tight-binding calculations an enhanced interlayer hopping is used, which allows us to recover the physics of small angle graphene superlattice with smaller units cells.Liu et al. (2015); Carr et al. (2018b)

This low-energy tight-binding model aims to simplify the effect of removing a single $p_{z}$ orbital by just eliminating a site out of the tBLG. As demonstrated in the DFT calculations above, this is equivalent to attaching a H atom to one of the bilayer surfaces. Furthermore, in the following we will be interested in the localization properties of the states in the moiré supercell. This can be easily computed within the tight-binding framework by means of the inverse participation ratio ($IPR$) of each Bloch wave function $\Psi_{n}$, that provides an estimation of the degree of localization of a state in a moiré unit cell: $IPR=\sum_{i}|\Psi_{n}(i)|^{4}$. For a fully extended state in the unit cell the $IPR$ is a small value of $1/N$, where N is the number of atoms, whereas for a localized state it yields a large finite value even for large unit cells.

The emergence of magnetism is explored introducing a local Hubbard interaction term of the form

[TABLE]

where $n_{i,\uparrow}=c^{\dagger}_{i,\uparrow}c_{i,\uparrow}$ and $n_{i,\downarrow}=c^{\dagger}_{i,\downarrow}c_{i,\downarrow}$ The previous interaction is solved at the mean-field level $H_{U}\rightarrow H_{U}^{MF}$ with

[TABLE]

We first address the case of a tBLG with a rotation angle of $\alpha\approx 3^{\circ}$, taking a single H atom per moiré unit cell (Fig. 5). Fig.5a shows that the paramagnetic band structure displays a flat band meeting the linearly dispersive bands, yielding a triple-point at the $K$ and $K^{\prime}$ points. According to the color-resolved diagrams of the $IPR$ (Fig. 5a), the nearly flat band at $E=0$ is the most localized state in the supercell, which corresponds to the zero mode that build around the impurity (Fig. 5c). Including electronic interactions in the model leads to a spin-dependent band splitting (Fig.5b), with the consequent emergence of magnetism (Fig.5d). The real-space representation of the states associated to the flat band in Fig.5c demonstrates that they are the same that become magnetic (Fig. 5d). This results exemplify that the phenomenology of mono-hydrogenated tBLG obtained from the DFT calculations is captured by the low-energy tight-binding model.

We now move on to study the fate of the impurity bands as the magic angle regime is approached. In particular, at the magic angle regime, an additional set of bands is expected to appear, namely the flat AA bandsSuárez Morell et al. (2010). In this situation, it is expected that the impurity band will heavily hybridize with the AA bands. To dive into this regime, a single H atom in a tBLG with twisting angle of $\alpha\approx 1.3^{\circ}$ (Fig. 6) is first considered. As shown in Fig. 6a), c), and e), a set of additional bands with highly reduced bandwidth appear close to the charge neutrality point. Again, a triple point is observed at the $K$ and $K^{\prime}$ points in the band structure, for an impurity deposited on the AA region (Fig. 6ab), on the AB region (Fig. 6cd) and in the AB/BA boundary (Fig. 6ef). Interestingly, two additional bands appear close to charge neutrality point, having a substantially strong $IPR$. Those new bands appear due to the interplay of AA band with the impurity bands, that drifts spectral weight from the impurity state slightly above and below the charge neutrality point. Importantly, the small bandwidth of the AA bands suggests that the system may accommodate extra electrons in the AA modes, instead of filling the impurity modes. This last feature becomes important in the next section when dealing with the interplay between magnetism and magic angle superconductivity.

IV Defect-induced Yu-Shiba-Rusinov states in twisted bilayer graphene

As the twist between the two layers approaches a rotation angle of $\alpha\approx 1^{\circ}$, the density of states of pristine twisted bilayer graphene is largely enhanced, triggering an electronic instability when the system is doped with nearly 2 holes per unit cell.Cao et al. (2018a); Lu et al. (2019); Yankowitz et al. (2019b) In the following, we will explore the interplay between the vacancy induced magnetism presented above, and the superconducting state found in magic angle superlattices.Cao et al. (2018a); Lu et al. (2019); Yankowitz et al. (2019b) To demonstrate that a magnetic defect has an impact on the tBLG superconducting propertiesCao et al. (2018a); Yankowitz et al. (2019b); Lu et al. (2019), the low-energy spectrum of a magnetic symmetry-broken tBLG is now studied by taking the Hamiltonian

[TABLE]

where $H_{SC}=\sum_{i}\Delta_{i}[c_{i,\uparrow}c_{i,\downarrow}+c^{\dagger}_{i,\downarrow}c^{\dagger}_{i,\uparrow}]$. For the sake of concreteness, we take $\Delta_{i}=\Delta$, realizing a spatially uniform s-wave superconducting singlet pairing.Wu et al. (2018); Choi and Choi (2018); Polshyn et al. (2019); Lian et al. (2018) We point out that additional superconducting symmetries have been proposed for magic angle superlattices besides the one consider above.Xu and Balents (2018); Isobe et al. (2018)

We first consider the case without superconducting order, namely $\Delta=0$. The non-interacting band structure for an angle $\alpha\approx 1^{\circ}$ tBLG mono-hydrogenated in the AB region is shown in Fig.7a. The triple-point survives in this small angle regime and the flat band is strongly hybridized with the nearly flat Dirac cones, located in the AA region. In addition, new bands appear close to charge neutrality carry spectral weight of the impurity state. The fate of the impurity-induced magnetic state triggered by electron interactions is considered next. When introducing interactions and doping with two holes per moiré unit cell, the impurity induced magnetic moment survives (Fig.7b) and e)). This strikingly compares with the large angle scenario, in which doping with a single electron per moiré unit cell already destroys the impurity-induced magnetism.

We now unveil the effect of superconductivity in the previous setup. For that sake, a superconducting term is included on top of the normal state Hamiltonian. In a system with time-reversal symmetry and uniform pairing, the BdG eigenvalues $\epsilon_{k}$ are given by $\epsilon_{k}=\pm\sqrt{E^{2}_{k}+\Delta^{2}}$, where $E_{k}$ are the Bloch eigenvalues in the absence of superconductivity. Namely, the BdG band structure of a nonmagnetic system shows a gap of $2\Delta$ in the spectra. This is observed in Fig.7c, where the superconducting field was included on top of the paramagnetic non-interacting band structure.

In the presence of magnetic moments, states inside the superconducting gap can appear due to the breaking of time reversal symmetry. In particular, Yu-Shiba-Rusinov statesBalatsky et al. (2006); Ruby et al. (2016, 2018, 2015); Wehling et al. (2008) are in-gap modes in the superconducting spectrum that appear as a result of the time-reversal symmetry breaking introduced by local magnetic moments. Adding a uniform superconducting pairing $\Delta$ to the self-consistent solution of Fig.7b, two in-gap bands are observed in the BdG spectrum, as shown in Fig.7d. Interestingly, the Yu-Shiba-Rusinov state (Fig. 7f) reflects the spatial structure of the magnetic state (Fig. 7e). Ruby et al. (2016, 2018, 2015) This results highlight that the impurity-induced magnetic moment of the tBLG can coexist with the superconducting state, generating in-gap Yu-Shiba-Rusinov states that could be observed by means of scanning tunnel microscope.Ruby et al. (2016, 2018, 2015)

Finally, it is interesting to note that the energy of this Yu-Shiba-Rusinov states cannot be easily estimated as in conventional metals hosting magnetic impurities.Balatsky et al. (2006) In those instances, the energy of the in-gap state generated in the superconducting states takes the formShiba (1968); Balatsky et al. (2006) $\epsilon=\Delta\frac{1-(JS\pi\rho/2)^{2}}{1+(JS\pi\rho/2)^{2}}$ where $J$ is the exchange coupling between the impurity and the conduction electrons, $S$ is the spin of the impurity, and $\rho$ is the density of states in the normal state. In mono-hydrogenated tBLG the density of states in the normal state is divergent due to the existence of nearly flat bands, and the onset of magnetism changes substantially the density of states and the Fermi energy, invalidating the previous procedure. This phenomenology is also found in hydrogenated. monolayer grapheneLado and Fernández-Rossier (2016) Nevertheless, the emergence of Yu-Shiba-Rusinov states is a robust feature, as shown in the exact solution presented above.

V Conclusions

To summarize, we have presented a detailed description of the impact of two types defects on the electronic and magnetic properties of twisted bilayer graphene. Localized states are found for both a vacancy and mono-hydrogenation of tBLG, which yield spatially localized magnetic states. Both types of defects remove effectively a pz electron out of the bilayer. On one hand, mono-hydrogenation defect states modify the low-energy band diagram at the charge neutrality point creating a triple-point. On the other hand, a carbon vacancy yields impurity bands at different energies, distorting the $\pi$ bands at the Fermi level. Interestingly, the presence of impurities leads to the emergence of spatially localized magnetic moments whose magnitude depends on the location of the defect.

In the magic angle regime, doping with two holes per moiré unit cell does not destroy the hydrogen induced magnetic moment, suggesting that such a defect-induced magnetism may coexist with the superconducting state. Including a superconducting term in the Hamiltonian, the interplay between defect-induced magnetism and superconductivity gives rise to in-gap Yu-Shiba-Rusinov states in the twisted bilayer. Our results demonstrate that atomic imperfections may have a substantial impact in the spectral properties of superconducting twisted graphene superlattices.

Acknowledgments

Los Alamos National Laboratory is managed by Triad National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy under Contract No. 89233218CNA000001. This work was supported by the U.S. DOE Office of Basic Energy Sciences Program E3B5. We thank I. Brihuega, E. Cortes-del Rio and A. Ramires for fruitful discussions. A.L.-B. acknowledges the computing resources provided on Bebop, the high-performance computing clusters operated by the Laboratory Computing Resource Center at Argonne National Laboratory. J. L. L. acknowledges financial support from the ETH Fellowship program.

Appendix

Appendix A First-principles calculations

The description of the coupling between the graphitic layers before and after the introduction of an impurity was conducted through self-consistent calculations with the SIESTA code within a localized orbital basis set scheme. Spin-polarized calculations were conducted using a double-$\zeta$ basis set, and the local density approximation (LDA) approach for the exchange-correlation functional was used. Atomic positions of systems formed by over 508 atoms were fully relaxed with a force tolerance of 0.01 eV/Å. The integration over the Brillouin zone (BZ) was performed using a Monkhorst sampling of 7 $\times$ 7 $\times$ 1 k-points. The radial extension of the orbitals had a finite range with a kinetic energy cutoff of 50 meV. A vertical separation of 25 Å in the simulation box prevents virtual periodic parallel layers from interacting.