Cooper-Pairs Distribution Function of Untwisted-Misaligned Bilayer Graphene

1. Introduction

Twisted Bilayer Graphene (TBLG) exposed to a constant electric field undergoes a zero-resistance state at a temperature below 1.7K [1]. The TBLG has remarkable properties; superconductivity [1], Mott insulator behavior [2], suppressing group velocity at certain angles and its relationship with electronic correlation [3,4,5,6]. The relative twisted angle between layers, in 2D materials such as TBLG [7,8], has become an undeniable phenomenon and a new field. In this way, magic angle value provided a novel knob for tuning superconducting properties [1] and many-body effects [2,9,10]. By contrast, recently superconductivity was reported in non-twisted rhombohedral trilayer graphene, a system without moiré patterns [11]. This could lead one to think that a specific piece of a TBLG could originate the superconductor state. Nevertheless, to achieve this specific piece, from a technical point of view, it is easier to twist a set of layers than to achieve very tiny displacements between them.

These quite specific structures drive towards providing needed characteristics to the superconductivity become possible. In the end, all of them could have in common a fundamental characteristic: the Cooper pairs. To set superconductivity up, the electrons belonging to a Cooper pair and their interaction by phonon must keep up energy and momentum conservation laws. Also, this kind of interaction must have enough energy to overcome repulsive Coulomb interaction so that the superconductor phase becomes feasible. Each of these conditions is corresponded with a probability. If then a sum of overall electronic energies are performed over the product of these probabilities, it turns to end a simultaneous probability to get a Cooper pair as a function of lattice vibrational energy and $T_c$, namely $D_{cp}(\omega ,T_c)$ function [12]. This function is built from the well-established Eliashberg spectral function and phonon density of states, calculated by first principles (for more details see reference [12]). Cooper-pairs distribution function $D_{cp}(\omega ,T_c)$ promises to provide information about the superconductor state through the determination of the spectral regions of Cooper-pairs formation. In addition, from $D_{cp}(\omega ,T_c)$ function is possible to obtain the N${}_{cp}$ parameter, which is proportional to the total number of Cooper-Pairs formed at temperature $T_c$ [13,14]. $D_{cp}(\omega ,T_c)$ function has been used to evaluate some superconductor systems [12,13,14,15].

In this paper, we aim to contribute to the understanding of superconducting fundamentals of challenging structures such as TBLG. We study an untwisted-misaligned bilayer graphene (UMBLG) structure, which mimics a specific region of a TBLG lattice, using Cooper-pairs distribution function $D_{cp}(\omega ,T_c)$. Here, we report and analyze $D_{cp}(\omega ,T_c)$ function of UMBLG structures in the presence of a perpendicular external electric field, at temperature $T=0.5$ K. The findings show that superconductivity in these structures responds to slight relative displacement between graphene sheets and to external electric fields, as in magic angle TBLG.

2. Materials and Methods

In order to build the $D_{cp}(\omega ,T_c)$ function, calculations of electronic density states, vibrational density states, and the Eliashberg function of the concerned structure are required. To perform these ab initio calculations, we first relaxed the internal degrees of freedom and the lattice vectors of the crystal structure using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) quasi-Newton algorithm [16,17]. From these relaxed structure configurations, electronic and phonon band structures, electron (DOS) and phonon (PhDOS) densities of states, and the Eliashberg function ${\alpha }^{2}F\left(\omega \right)$ are calculated, for this study, in the presence of perpendicular electric fields; $E_0$, 3$E_0$, and 5$E_0$, with $E_0=9$ mV/Å. We used a kinetic energy cut-off of 70 Ry for the expansion of the wave function into plane waves and 280 Ry for the density. To integrate over the Brillouin zone, we used for the electronic integration a k-grid of 32 × 32 × 1 and the phononic integration a q-grid of $8\times 8\times 1$, according to the Monkhorst–Pack scheme [18]. We performed the calculations using the pseudopotential plane-wave (PW) method of Perdew et al. [19], the generalized gradient approximation (GGA), and a Troullier and Martins [20] norm-conserving pseudopotential. The cut-off and grids were chosen big enough to obtain a good precision in ${\alpha }^{2}F\left(\omega \right)$ calculated within the density-functional perturbation theory (DFPT) frame [21,22]. We used the Quantum Espresso code [23] for all these calculations.

2.1. Cooper-Pairs Distribution Function

Conventional superconductivity can be explained by an attractive electron-electron interaction (Cooper pair) mediated by the lattice phonons. Cooper-pairs formation is induced by the simultaneous occurrence of specific physical conditions, which can be described as a product of probabilities in terms of electronic and vibrational states, and electron-phonon interaction of a system, that summed over all electronic states defines Cooper-pair distribution function $D_{cp}(\omega ,T_c)$, given by

$$D_{cp}\left(\omega ,T\right)={\int }_{E_F-{\omega }_c}^{E_F+{\omega }_c}{\int }_{E_F-{\omega }_c}^{E_F+{\omega }_c}{g}_{ep}^s(T,ϵ,\omega )\times g_{ep}^b(T,{ϵ}^{\prime }+\omega ,\omega )\times {\alpha }^2\left(\omega \right)dϵd{ϵ}^{\prime },$$

(1)

where $g_{ep}^s(T,ϵ,\omega )=g_e^o(T,ϵ)g_e^v(T,ϵ+\omega )g_p^n(T,ϵ)$ and $g_{ep}^b(T,ϵ,\omega )=g_e^o(T,ϵ)g_e^v(T,ϵ-\omega )g_p^{n+1}(T,ϵ)$ are the simultaneous probability at T related with the presence of electronic states occupied and available (near to the Fermi level), as well as optimal vibrational frequencies, which are candidates for the consolidation of a attractive interaction potential, that drives Cooper pair formation. Further, ${\alpha }^2\left(\omega \right)={\alpha }^{2}F\left(\omega \right)/N_{ph}\left(\omega \right)$ establishes the electron-phonon coupling probability [24,25], through the well-known Eliashberg spectral function and phonon density of states. In this sense, $D_{cp}(\omega ,T_c)$ function provides the spectral range where Cooper pairs could be formed.

Furthermore, from $D_{cp}(\omega ,T_c)$ function it is possible to estimate the total number of Cooper pairs formed at the temperature through a quantity proportional to it, the parameter $N_{cp}\left(T\right)$, given by

$$N_{cp}\left(T\right)={\int }_0^{{\omega }_c}{D}_{cp}\left(T,\omega \right)d\omega ,$$

(2)

Where $\omega$ is phonon energy and ${\omega }_c$ is a cutoff energy, so that to $\omega >{\omega }_c$$D_{cp}(\omega ,T_c)$ is negligible.

2.2. Structure

Figure 1 illustrates the AB-staked TBLG structure. The top view, in Figure 1a, shows two hexagonal layers, where one is rotated with respect to the other at a small angle ${\theta }_m$, and Figure 1b illustrated the side view showing the corresponding interlayer distance d. Figure 2 summarizes the geometric concepts used to determine the structural, electronic and vibrational properties of the AB-stacked TBLG structure. The rotation angle ${\theta }_m$ in a commensurate rotation condition ($\theta ={\theta }_m$) [26] is determined as following:

$$cos\left({\theta }_m\right)={\displaystyle \frac{3m^2+3m+1/2}{3m^2+3m+1}},$$

(3)

with $m=0,1,2,3,...$, allows long-wave-length Moiré pattern with ${\lambda }_m=a/(2sin({\theta }_m/2))$ [27].

With the purpose of establishing a possible UMBLG structure associated with a TBLG structure, we initially consider a supercell of an AB-stacked BLG rotated at a specific small angle $\theta$ (Figure 2b), in which were identified two equivalent regions where AB-staked TBLG configuration are situated particularly around a rotation axis (see Figure 2a). In the midpoint of these regions, corresponding to ${\lambda }_m/2$, we encounter a lattice arrangement akin to Misaligned Bilayer Graphene, which is observed as two AB-stacked layers with a parallel displacement between them, as illustrated in Figure 2c. We established that the value of this displacement $D_m$ is linked to the twisted angle ${\theta }_m$ of the TBLG structure associated. Moreover, this displaced structure can be easily built from an untwisted bilayer graphene (Figure 2d). Table 1 shows the geometrical characteristics of three configurations, and the respective displacement values $D_m$ associated to each twisted angle ${\theta }_m$.

The structure proposed in this study consists of four carbon atoms (two per-layer) whose positions and internal parameters generate two displaced-layers graphene layers, which reproduce the AB-stakes BLG structure in the zone ${\lambda }_m/2$. In this cell unit, a of 12 Å vacuum gap was used to eliminate interactions between periodic images perpendicular to the slab surface.

3. Results and Discussion

Relaxation procedure leads to a lattice constant of $a=2.47$Å, which is in good agreement with other theoretical reports [22,28]. On the other hand, we find out the inter-layer distance $d=4.40$Å, Figure 1b. By contrast, experimental and theoretical reports are 3.35 Å and 3.31 Å [28,29] to AB-stake BLG. The discrepancy is due to the bilayer has been put isolated. To do so, we set the AB-stakes bilayer up into a cell with a vertical length of more than eleven times the lattice constant to dismiss the inter-bilayer interaction. We used the Quantum Espresso code [23] for all these calculations.

Figure 3 shows the Eliashberg function ${\alpha }^{2}F\left(\omega \right)$ and Phonon density of states (PhDOS) calculated for AB-stacked Bilayer Graphene proposed in presence of a perpendicular electric field proportional to $E_0=9.0$ mV/Å, and AB-stacked Untwisted-Misaligned Bilayer Graphene (associated with ${\theta }_{33}=0.{987}^{\circ }$ and ${\theta }_{32}=1.{018}^{\circ }$ angles of a TBLG) in presence of a perpendicular electric field of $5E_0$. These spectra show four characteristic peaks located approximately at 60, 80, 170 and 200 meV, which are consistent with the theoretical calculations of AB-stacked Bilayer Graphene (with 10 000 atoms) reported by Choi et al. [30,31] using ${\theta }_m=1.{08}^{\circ }$, and PhDOS of TBLG at ${\theta }_m={13}^{\circ }$ and ${21}^{\circ }$ (with 28 and 76 atoms, respectively) reported by Cocemasov et al. [32].

It is observed from Figure 3 that the ${\alpha }^{2}F\left(\omega \right)$ and PhDOS are nearly insensitive to the small displacements associated with twist-angle (${\theta }_m$) as well as to the presence of an electrical field. Besides, the insensitive of PhDOS (TBLG) to the twist-angle was observed by Cocemasov et al. [32]. Therefore, it is almost not possible to establish significant differences induced by displacements (associated with ${\theta }_m$) or $E_0$, more than those observed below 5 meV (see inset of Figure 3). At first sight, these slight differences do not seem to offer significant information about variations in the MUBLG properties. However, with these spectra, it is possible to obtain the Cooper-pairs distribution function $D_{cp}(\omega ,T_c)$, which allows evaluating the system in the superconducting state [12,14].

Figure 4 shows the $D_{cp}(\omega ,T_c)$ function of AB-stacked Untwisted-Misaligned Bilayer Graphene (associated with angles ${\theta }_m=0^{\circ },0.{987}^{\circ }$ and $1.{018}^{\circ }$) in presence of an electrical field. $D_{cp}(\omega ,T_c)$ function of MUBLG suggests that Cooper-pairs are possible in the low-frequency vibration region. By contrast, BLG with a perpendicular electric field of magnitude $E_0$ (Figure 4a), it is in the normal state (non-superconducting). It is inferred because it reports a $D_{cp}(\omega ,T_c)\approx 0$. Namely, the superconducting phase is theoretically observing if $D_{cp}(\omega ,T_c)$ function is not zero, which, in this case, is achieved when Bilayer Graphene is displaced (associated with a rotation small angle) in the presence of the electric field. This condition matches the experimental reports (superconductivity in magic-angle graphene superlattices [1]). It suggests the origin of superconductivity in this kind of structure could be a little misaligned piece. It is also evident from Figure 4 that the $D_{cp}(\omega ,T_c)$ function reveals a significant sensitivity to displacements associated with small variations of ${\theta }_m$ as well as to the presence of an electric field (note the scale differences between Figure 4 a) and b). This is possibly due to the correlational analysis of the electronic and phononic properties of the system that the $D_{cp}(\omega ,T_c)$ function achieves, which cannot be determined from an individual evaluation of such properties, as was evidenced with the results presented in Figure 3. However, the sensitivity of the $D_{cp}(\omega ,T_c)$ function is a criterion that requires more calculations to be validated, in addition to those already carried out (H${}_3$S [12], Nb [13] and D${}_3$S [14]).

$D_{cp}(\omega ,T_c)$ function of MUBLG shows a quasi-gaussian form, which is appreciably affected by ${\theta }_m$ and the presence of an electric field. It is observed that an increase in the intensity of the electric field leads to an increase in the area under the spectral line and moves it slightly towards higher values of frequencies, for both angles. However, the effect was more marked for ${\theta }_m=0.{987}^{\circ }$ (note the difference in the scale of the vertical axes in Figure 4. $D_{cp}(\omega ,T_c)$ is only observed at frequencies below $0.4$ meV (centered around $0.05$ meV), which seems consistent with the low critical temperature of TBLG ($Tc=0.5$ K [1]).

Furthermore, from $D_{cp}(\omega ,T_c)$ function it is possible to get an estimate of the total number of Cooper pairs formed at temperature T${}_c$ through a quantity proportional to it, $N_{cp}$ parameter [13,14]. The comparison of $N_{cp}$ parameters obtained from each $D_{cp}(\omega ,T_c)$ as a function of an electrical field, for both misaligned structures which are associated with two angles, is shown in Figure 5.

It is observed in Figure 5 that an electric field has a greater effect on $N_{cp}$ parameters obtained of UMTBLG associate with the angle $\theta =0.{987}^{\circ }$, compared to $\theta =1.{018}^{\circ }$. This could be because the electric field affects the flat band physics (dispersion and topological properties) [33,34], which are associated with the superconducting phase. Since a high $N_{cp}$ value implies a high $T_c$ [13,14], hence it can be suggested from Figure 4 that an increase in the intensity of the electric field could lead to an increase in $T_c$ in TBLG, for a specific angle of rotation. Besides, calculations to $D_{31}=0.122988$ nm (${\theta }_{31}=1.{05}^{\circ }$) were performed also. Outcomes showed the regarding $D_{cp}(\omega ,T_c)$ function much less than $D_{32}=0.123024$ nm (${\theta }_{32}=1.{018}^{\circ }$) in case, so it was not included in Figure 4 and Figure 5.

Finally, despite the consistent results shown up to now, some expression of the characteristics flat bands of TBLG [33] are not observed directly in the version of $D_{cp}(\omega ,T_c)$ used in this study. This requires further analysis.

4. Conclusions

In this study, we reported for the first time, to the best of our knowledge, the Cooper-pairs distribution function $D_{cp}(\omega ,T_c)$ of Untwisted-Misaligned Bilayer Graphene (UMBLG), under the influence of external electric fields. The UMBLG structure was proposed as a structural geometric approximation that mimics a characteristic region of a Twisted Bilayer Graphene (TBLG) structure. This UMBLG structure was modeled as two AB-stacked layers with a parallel displacement $D_m$ between them. We consider three displacements $D_{33}=0.12295$ nm, $D_{32}=0.12302$ nm, and $D_{31}=0.12299$ nm, which are linked with the rotation angles (near to magic-angle) ${\theta }_{33}=0.{987}^{\circ }$, ${\theta }_{32}=1.{018}^{\circ }$, and ${\theta }_{31}=1.{05}^{\circ }$, respectively, associated to Twisted Bilayer Graphene (TBG). $D_{cp}(\omega ,T_c)$ function of Bilayer Graphene in the presence of an electric field was also calculated, as a reference. As it was expected, the reference structure showed a negligible $D_{cp}\left(\omega {,}_c\right)$ function due to it is not a superconductor one. However, $D_{cp}\left(\omega {,}_c\right)$ function of UMBLG structures with $D_{33}$ and $D_{32}$ displacements showed characteristics indicative of superconducting behavior, which became more noticeable with the increase of the electric field. This behavior was greater for $D_{33}$(${\theta }_{33}=0.{987}^{\circ }$) than for $D_{32}$(${\theta }_{32}=1.{018}^{\circ }$). Furthermore, $D_{cp}\left(\omega {,}_c\right)$ function of these structures suggests that Cooper-pairs are possible in the low-frequency vibration region. In addition, the calculation of the $N_{cp}$ parameters allowed us to suggest that for specific displacements, an increase in $N_{cp}$ could be induced by applying and increasing the external perpendicular electric field. This implies an increase in the number of Cooper-pairs and therefore an increase in the superconducting critical temperature. The displacement $D_{31}$ in the UMBLG structure indicated that it does not exhibit superconductivity under this particular configuration.All these results demonstrated that the structural geometric approximation appears to be consistent with the presence of the superconducting state observed in TBLG at small angles, which allowed evaluation of the possible effect of the electric field.