Classifying superconductivity in Moiré graphene superlattices

Several research groups have reported on the observation of superconductivity in bilayer graphene structures where single atomic layers of graphene are stacked and then twisted at angles θ forming Moiré superlattices. The characterization of the superconducting state in these 2D materials is an ongoing task. Here we investigate the pairing symmetry of bilayer graphene Moiré superlattices twisted at θ = 1.05°, 1.10° and 1.16° for carrier doping states varied in the range of n = (0.5 − 1.5) · 1012 cm−2 (where superconductivity can be realized) by analyzing the temperature dependence of the upper critical field Bc2(T) and the self-field critical current Jc(sf,T) within currently available models – all of which start from phonon-mediated BCS theory – for single- and two-band s−, d−, p− and d + id-wave gap symmetries. Extracted superconducting parameters show that only s-wave and a specific kind of p-wave symmetries are likely to be dominant in bilayer graphene Moiré superlattices. More experimental data is required to distinguish between the s- and remaining p-wave symmetries as well as the suspected two-band superconductivity in these 2D superlattices.


Introduction
For an isotropic, spherical Fermi surface the density of states at the Fermi level is given by: (  ) = 8 ℎ 3 • ( * ) 2 •   (1) where h is the Planck constant, m * is the effective mass of the charge carriers, and vF is the Fermi velocity.Due to their large effective mass of  * ~200 •   (where me is the electron mass) heavy fermion superconductors [1] possess a robust superconducting state, characterized by high values of the upper critical field Bc2 [2] (significantly above the paramagnetic Pauli limit Bp) despite a low Fermi velocity,   ~5 • 10 3 /, in these compounds [2].
On the other hand, Eq. 1 prohibits a superconducting state in materials possessing a charge carrier effective mass of  * < 0.1 •   .The is because a large value of D(EF) requires relativistic values for the Fermi velocity,   ≳ 10 8 /, which has not yet been observed in any material.In single layer graphene (SLG) and other materials with Dirac cone Fermi surfaces, D(EF) is given by: which shows that D(EF) is inversely proportional to   2 .Therefore, the prerequisite to convert SLG and other planar honeycomb lattices [3] into intrinsic superconductors is a reduction of   .
Lopes dos Santos et al [4] were the first to understand that   can be suppressed in bilayer graphene by rotating the SLG sheets relative to each other by small twist angles .Detailed tight-binding model calculations performed by Bistritzer and MacDonald [5] showed that at the Dirac points   goes to zero when bilayer graphene is rotated by small angles, creating a Moiré superlattice.These graphene 2D structures are now called magic-angle twisted bilayer graphene (MATBG).
The discovery of intrinsic superconductivity in few-layer stanene (one of the closest analogues of graphene) which possesses a Fermi velocity of   = 4.5 • 10 4   ⁄ by Liao et al [6] gave clear experimental proof that the suppression of vF in Dirac cone materials makes it possible to convert these materials into intrinsic superconductors.Several months after Liao's work [6], intrinsic superconductivity in MATBG with   ~2 • 10 4   ⁄ and  * ~0.2 •   was reported by Cao et al [7].There, simultaneous vF suppression and  * enhancement was achieved in the way considered by Lopes dos Santos et al [4] as well as Bistritzer and MacDonald [5], i.e. by rotation of 2 stacked SLG sheets to form a Moiré superlattice.For instance, Cao et al [7] showed that MATBG with  = 1.16° exhibits zero resistance at Tc = 0.14 K, whereas an angle  = 1.05° has zero resistance at Tc ~ 1.2 K.More recently, several research groups have discovered a superconducting state in MATBG [8][9][10][11], including studies of MATBG at high-pressures (up to P = 4 GPa), as well as in trilayer graphene/boron nitride Moiré superlattices [12].
There are still many important questions regarding the superconducting state of MATBG that remain unanswered.In particular, the pairing mechanism and symmetry, as well as the number of superconducting bands, is not clear.While Cao et al [7] propose strong electronelectron correlations as the mechanism of the superconductivity, a full analysis with respect to existing theories for phonon mediation is still required.When treated as a phonon mediated superconductor the possibility of multiple superconducting bands arises from the fact that the charge carriers in MATBG experience two different phonon modes, one from intralayer covalent bonds, and a second mode due to the interlayer coupling (details can be found in [13,14]).
Here, we analyse temperature dependent measurements of the upper critical field and selffield critical current in bilayer twisted SLG structures reported by Cao et al. [7] and Lu et al. [11].This analysis is done in the existing phenomenology developed for phonon-mediated

Upper critical field analysis: single superconducting band models
We analyse the temperature dependent perpendicular upper critical field Bc2,(T) measured by Cao et al [7,Figure 3e] where the magnetic field is applied in the perpendicular direction to the MATBG plane.First we present the standard literature analysis for Bc2(T) data by fitting the data to the 5 available models for Bc2,(T) and extracting Tc as well as the ξ(0).
However, we extend this analysis by noting that the first 5 models presented have an implicit assumption of the behaviour of the penetration depth λ(T).By dropping this assumption, we can express Bc2,(T) as a function of the band gap Δ(T) and its symmetry.
To begin, we fit the Bc2,(T) data to the Gorter-Casimir (GC) two fluid model [15,16] which is still in wide use [17,18]: where  0 = Another widely used model for fitting Bc2(T) data is the Werthamer-Helfand-Hohenberg (WHH) model [19,20]: a fit of the data to this model is also shown in Fig. 1(a).Baumgartner et al [21] proposed a simple and accurate analytical expression that matches the shape of the WHH model: This model will be designated as B-WHH and its fit of the data are shown in Fig. 1(b), where both values of Bc2,(T = 0 K) (i.e., original WHH and B-WHH) agree as expected.
Jones-Hulm-Chandrasekhar (JHC) proposed three models [22], two of which are often used to fit Bc2(T) data [23][24][25].The first is a combination of Ginzburg-Landau (GL) theory with the Gorter-Casimir expression for the temperature dependence of λ(T): which we shall refer to as the JHC model, a fit using this model is shown in Fig. 1(c).This model gives a little higher value for Bc2,(T = 0 K) and Tc in comparison with the other models, as well as the lowest value of ξab (0).The second model proposed by JHC is based on an expression from Gor'kov for Bc2(T) [26]: where Bc(T) is the thermodynamic critical field, and (0) is the ground state London penetration depth.Jones et al [22] again use the Gorter-Casimir form of λ(T) to produce a model of the following form: ) .
Eq. 7 will be referred to as the Gor'kov model.
However, rather than assume the behavior of the penetration depth, eq.6 can be used to build a model based directly on the BCS expression [16,27] for λ(T) as a function of the superconducting gap: where (T) is the temperature-dependent superconducting gap.Eq. 8 captures the effect of the gap symmetry on Bc2(T).An expression for (T) as a function of the wavefunction symmetry is given by Gross-Alltag et al [28].For now, we test s-wave symmetry where the gap takes the form: where ΔC/C is the relative jump in electronic specific heat at Tc, and  = 2/3.We can then use the standard GL expression: •   (), (10) to restate Gor'kov's expression, eq.6, as: • (1.77 − 0. ).
Then, by considering another GL expression: we can combine eqs.(8)(9)(10)(11) and eq. 12 into a single expression for the temperature dependent upper critical field: Eq. 13 is a function parameterized by four fundamental superconducting parameters, ab(0), (0), C/C, and Tc.It is critical to note that the role each parameter plays in the equation is well founded in BCS and GL theory.In particular ab(0) determines the absolute value of Bc2,(0), independently of the other parameters.Furthermore, Tc is tightly constrained to define a reduced temperature scale t = T/Tc, over which Gor'kov's expression, eq.6, is valid.The only new feature of eq. 13 is the introduction of the gap symmetry expression for λ(T), eq. 8.This introduction, and that of the Gross-Alltag et.al. expression for the gap (T) itself, only introduces degrees of freedom, (0) and C/C, whose behavior were phenomenologically assumed by the literature models.The advantage of this method is that the two new parameters can be deduced and explicitly checked against the bounds of the theory's validity.
A fit of the MATBG Bc2(T) data to Eq. 13 is shown in Fig. 1(e) where the superconducting gap ratio was found to be: and the relative jump in electronic specific heat: These deduced parameters, within their uncertainties, are remarkably close to the BCS weak-coupling limits of 43.This is the first quantitative evidence that intrinsic superconductivity in MATBG can be understood in the existing phenomenology of s-wave electron-phonon mediated superconductivity and a new phenomenology may not need to be developed.Overall, the deduced values for ab(0) and Tc using the six models are (see for details Supplementary Table S1):   (0) = 61.4± 1.7    = 0.500 ± 0.006 .

MATBG ( = 1.16°)
Eq. 12 magnetic flux density (mT) Once a value of the ground state gap has been obtained another standard BCS expression [16,27] can be used to calculate vF for the sample M1 ( = 1.16°) based on the deduced values of ab(0) and (0): The value given in Eq. 16 is about two orders of magnitude lower than the Fermi velocity of pure SLG ( , ~ 10 6   ⁄ ) [29].
The Fermi temperature, TF, can be calculated in the usual way: where m * is the effective mass of charge carriers.In the original paper by Cao et al [7], m * was not measured for sample M1 ( = 1.16°).A simple assumption is that at the same doping, MATBG with different  should have the same effective mass, m*.By assuming this, and using Bc2(T) as an indication of the doping state we can assume that for sample M1: (we present full analysis for this value for sample M2 in the following section).
Corresponding to this m* value the Fermi temperature, TF, is: We plot this data point for MATBG sample M1 in an Uemura plot [30] (Fig. 2) where the transition temperature Tc is defined as the temperature at which order parameter phase coherence is established i.e. when the sample resistance becomes zero, R= 0  (rather than R= 0.5Rn).This definition is in accordance with all other data points (for other materials) and follows the general logic of the Uemura et al [30,31] definition of the transition temperature.
From R(T) data for sample M1 presented in [7,Figure. 1b] we determine the transition temperature to be Tc = 0.136 K.One therefore obtains: which fits within the boundary of: established by Uemura et al. [30,31] for all other unconventional superconductors.

Figure 2.
A plot of Tc versus TF obtained for materials representative of the various superconducting families.Data was taken from Uemura [28], Cao et al [7], as well as [24,25].
Unfortunately, the existing data set is insufficient to reveal irreconcilable differences between the models analyzed.While this is reasonable considering how low the temperature must be taken, it is also unfortunate that Bc2,(T) data is not available for sample M2; where the higher critical temperature would lead to a lower acheivable range of reduced temperature t= T / Tc.
It is possible to fit the Bc2,(T) data to a model with d-, p-and d+id-gap symmetry by substituting the relevant expressions for (T) and (T) in Eq. 13, however, the Bc2,(T) data for the sample M1 [7] was not measured to low enough temperatures (T / Tc < 0.33) to accurately determine the fit parameters.
To resolve this issue we look at self-field critical current density data Jc(sf,T) which was measured by Cao et al [7] for the sample M2 ( = 1.05°) down to T = 0.05 K.These measurements, taking into account that Tc ~ 1.2 K, were taken down to a reduced temperature of T/Tc = 0.04 << 0.33.Therefore, the raw Jc(sf,T) data provides a meaningful insight into the superconducting gap symmetry of MATBG, which we will explore in a later section.
Upper critical field analysis: two superconducting band model.As we already mentioned above, several authors have proposed a two-band superconducting state in MATBG which originates from two different phonon modes.One is due to intralayer covalent bonds, and the other is due to interlayer coupling (an extended reference list and discussion can be found in Refs.13,14).It is interesting to note that in our previous papers [17,[32][33][34] we found that many atomically thin superconductors, including SLG and topological insulators (with a proximity induced superconducting state) have a two-superconducting band state.
Looking closer at the upper critical field, Bc2(T), fit to the single band GC model in Fig. 3 (a), it can be seen that in the temperature range of T = 0.3 -0.4 K the experimental Bc2(T) data is well below the fitted curve.A fit of the data to a decoupled-two-band GC model (which we proposed in Refs.17): is shown in Fig. 3 (b).The deduced parameters, see Fig. 3 (b), agree well with the mutual parameter interdependence criteria, which remains low and varies from 0.69 to 0.97 (the definition of this parameter can be found elsewhere [17]).Giving strong evidence that MATBG is a two band superconductor.)  (23) where V0 is an instrumental offset, k is a linear term used to accommodate incomplete current transfer in short samples, n is the flux creep exponent, and Vc is the voltage criteria which was chosen to be 10 V.The superconducting state can be defined through the V(I) fit to Eq. 23 by a n > 1 criterion, where the normal state corresponds to n = 1 (Ohm's law).Physically speaking, this approach defines detection of the superconducting phase coherence by detection of flux vortex flow in the sample.This distinction is necessary as data in Cao et.al.
The resulting Jc(sf,T) = Ic(sf,T)/(4ab) (where 2a is the sample M2 width, 2a = 1.05 m) is shown in Fig. 4. In this figure the Jc(sf,T) data was fitted using different superconducting gap symmetries by the equation [45]: where 0 is the magnetic permeability of free space.Eq. 24 is valid for thin superconductors if 2b < ) [17,18,[45][46][47] and this should be the case for 1.0 nm thick MATBG.
Unfortunately, raw Bc2,(T) experimental data, from which to deduce ab(0) for sample M2, is unavailable.However, it can be seen that ab(0) for both sample M1 and sample M2 is practically identical if we consider values at optimal doping.This follows by inspection of Bc2,(T) in [7, Figure 3e] for sample M1, where the comparison of Bc2,(T) for both samples is displayed in [7, Figure 3a] and [7, Figure 3b].Thus, for the analysis of sample M2 by Eq. 23 we will use the value ofab(0) = 61.4nm deduced from sample M1.This assumption is also supported by a recent report from Lu et al. [11] (in their Extended Data Figure 3), who measured Bc2(T=16 mK) for a broad range of doping states for MATBG with  ~ 1°.We analyse the Lu et al. [11] data in the following section.
However, of more importance is whether such gap symmetries are supported by experimental data.Therefore, we fit the available data using an extended BCS model with different expressions for the gap symmetry and compare the deduced parameters with weakcoupling BCS limits.Based on this we can infer which gap symmetries can potentially explain the measured behaviour.Here, we fit the experimental Jc(sf,T) data for sample M2 to Eq. 24 by utilizing different expressions for the symmetry of the superconducting gap.For these expressions, we again draw from the approach proposed by Gross-Alltag et al [28] for s-, d-, and p-wave symmetries.For d+id symmetry we use an approach proposed by Pang et al [57].
Equations for (T) and (T) for s-wave symmetry have already been presented in Eqs.8,9 respectively, and in Fig. 4(a) a fit of the Jc(sf,T) data to single band s-wave model is shown.
All the deduced parameters are presented in Table I.
Next, we performed a fit to a d-wave gap symmetry model.Using a 2D cylindrical Fermi surface the equation for the London penetration depth is as follows: where the superconducting energy gap, (T,), is given by: where m(T) is the is the maximum amplitude of the k-dependent d-wave gap given by Eq. 9,  is the angle around the Fermi surface subtended at (, ) in the Brillouin zone (details can be found elsewhere [58,59]).In Eq. 9 the value of  = 7/5 [58].The fit to this model does not converge, as the value m(0) tends toward an infinitely large value.
Table 1.Deduced parameters for single-band models applied to MATBG sample M2 [7] doped at nn=-1.44 10 12 cm -2 where an effective mass of charge carriers m * /me=0.1637±0.0154 was used, which was obtained from the analysis presented in Fig. 5.For the case of d+id we fixed the ratio of C/C = 0.995 for both gaps.The ground-state coherence length was assumed to be (0) = 61.4± 1.7 nm.

London penetration depth (m)
A fit can still be obtained by using a low-T asymptote of the single-band d-wave model which is presented in Fig. S4 It can be seen (Fig. S4, Table I) that the deduced m(0) = 2.3 ± 0.5 meV is unacceptably large.These results indicate that phonon-mediated d-wave symmetry in MATBG is not supported by experimental data and should be omitted from further consideration.
Fitting a p-wave gap symmetry model is more complicated (compared with s-and d-wave) because in this case the gap function is given by [28]: where, (T) is the superconducting gap amplitude, k is the wave vector, and l is the gap axis.
The electromagnetic response depends on the mutual orientation of the vector potential and the gap axis.In an experiment this is given by the orientation of the crystallographic axes compared with the direction of the electric current.There are two different p-wave pairing states: "axial" where there are two point nodes, and "polar" where there is an equatorial line node.The shape of the London penetration depth, (T) for p-wave polar Al and axial Al cases are difficult to distinguish from the s-wave counterpart, and the p-wave axial Al case is difficult to distinguish from the dirty d-wave case [28,59].
Despite these difficulties there is still a possibility to make a distinction based on the values of the deduced fundamental superconducting parameters, in particular by considering the ratios of and C/C.These two ratios are given in Table 2 for p-wave and other gap symmetries [28].
The result for only one p-wave configuration (axial Al) is shown in Fig. 4(b) as this has the only meaningfully deduced parameters.The other three fits are presented in Figs.S5.The deduced parameters for all cases are presented in Table 1.
The gap equation for the case of d+id gap symmetry is considered elsewhere [48,49,53,56].In this paper we use the approach proposed by Pang et al [57]: where 1(0) and 2(0) are the two d-wave gap amplitudes,  = 7/5, and  is the angle around the Fermi surface subtended at (, ) in the Brillouin zone, and αBCS,d+id-wave is the double band gap ratio (details can be found elsewhere [56,57]).
Because the experimental Jc(sf,T) dataset was not dense and the two gap model has many parameters, we were forced to reduce the number of free-fitting parameters so as not to overfit the data.We therefore have chosen to fix: 1. C/C was set to the weak-coupling limit of BCS theory for d-wave superconductors, C/C = 0.959.
2. The double band gap ratio  BCS,+−wave was set to its weak-coupling limit which is 1.925.
The fit values obtained for this d+id model (Table 1) give the two values of the gap ratios: where 2•Δ1 (0)/kBTc far exceeds the weak coupling limit of 3.85 for d+id symmetry superconductor [48].This therefore suggests that d+id symmetry should also be excluded from further consideration.
The analysis of the self-field critical current density in MATBG within a single band model shows that: 1. MATBG has an equal chance to be moderately strong coupled s-wave or p-wave superconductor.
2. The ground-state London penetration depth is independent of gap symmetry: 3. The GL parameter (0) is: The Fermi velocity vab,F and the Fermi temperature TF for the sample M2 ( = 1.05°) based on the deduced values of ab(0) and (0) for all gap symmetries have been calculated using Eqs.16 and 17 respectively, values are presented in Table 1.
All the calculated TF values are in presented in Table III together with the ratio Tc/TF.This ratio varied in the range: with the lower limit corresponding to the p-wave fit and upper limit corresponding to the swave fit.It should be noted that both lower and upper limits of this ratio are within the range (Eq.21) established for all unconventional superconductors by Uemura et al [30,31].
This value is also in the same range as other unconventional superconductors [30,31].
where each band is described by Eq. 24.Because each band has four free-fitting parameters, i.e. (0), Δ(0), ΔC/C and Tc, we were forced to reduce the number of parameters, and as we found in [17], the most appropriate approach is to equalize ΔC/C for both bands, i.e.: We also assume that for both bands c = 35.6 (Eq.36).The result of fits to two-band s-wave and two-band p-wave axial Al models are shown in Figs.6(a) and 6(b), respectively.There is experimental evidence that at T ~ 0.5 K a new superconducting band opens.The deduced parameters for both fits are in Table 3.It can be seen from Fig. 6 and Table 3, that there is not a significant differences between transition temperatures, Tc1 and Tc2, of each model.nor is there a significant difference comparing the difference between the BCS ratios of

MATBG phase diagram.
Here we use Eq.23 to analyse data reported by Lu et al. [11] on the phase diagram of MATBG with  ~ 1° measured for a wide range of doping states.We note that Lu et al. [11] clearly observe at least six superconducting domes shown in [11,  By using a sample D1 width, 2a = 2 m, and thickness, 2b = 1 nm, we calculate Jc(sf, T=16 mK) for four doping states where Ic(sf, T=16 mK) data was reported.By using the deduced ab(T=16 mK) for the same doping state we numerically solved Eq. 23 and deduced ab(T=16mK), c(T=16mK),  ,, , and  ,,   (Fig. 7).

IV. Conclusions
In this paper we use existing BCS and GL models to analyze raw experimental data from MATBG reported by Cao et al. [7] and Lu et al. [11].Surprisingly enough, the results of our analysis show that MATBG has a very low charge carrier effective mass,  * = (0.164 ± 0.015) •   , and is located in the Uemura plot (Fig. 2), next to the heavy fermion superconductors [1], particularly to UBe13 which has  * ~200 •   [2].This places MATBG in the same band of TF/Tc values as all other unconventional superconductors as categorized by Uemura, contrary to [7, Figure 6].
Our analysis of the temperature dependent upper critical field and the self-field critical current experimental data show that three of four p-wave as well as d-, and d+id-wave symmetries should be excluded from further consideration as possible phonon-electron mediated pairing symmetries in MATBG.Furthermore, the analysis indicates that MATBG is a moderately strong coupled two-band superconductor with s-or p-wave symmetry which can be categorized using the established phenomenology of superconductivity.Because graphene has planar honeycomb lattice of sp 2 bonded carbon atoms, our findings of either s-or p-wave pairing symmetry have supporting background evidence.
Supplementary Table S2.Deduced Ic(sf,T) from fit of V(I) curves to Eq. 23.

Temperature
067 • 10 −15 Wb is the superconducting flux quantum; composed of Plank's constant h and the electron charge e is, and ab(T) is the coherence length in the a-b plane.The fit is shown in Fig. 1(a).

Figure 4 .
Figure 4.The self-field critical current density, Jc(sf,T), for sample M2 ( = 1.05°) from the work of Cao et al [7] and a fit of the data to three single-band models.For all models ab(0) = 61.4nm was used.(a) s-wave fit, the goodness of fit R = 0.9964; (b) p-wave axial Al fit, R = 0.9971; (c) d+id-wave fit, R = 0.9967.Deduced parameters are listed in TableI.

Figure 5 .
Figure 5. Raw experimental m * (n)/me data with a linear fit where the goodness of fit was R=0.911.The green spot indicates the extrapolated m * (n)/me at n = -1.44•10 12cm -2 .

Figure 6 . 6 "
Figure 6.The self-field critical current density, Jc(sf,T), for sample M2 ( = 1.05°) from the work of Cao et al [7] and a fit of the data to s-and p-wave two-band models.For both models and both bands c(0) = 35.6 (Eq.35) was used.(a) s-wave fit, the goodness of fit R = 0.9965; (b) p-wave axial Al fit, R = 0.9965.Deduced parameters are in TableIV.

Figure 7 .
Figure 7. Analysis of the superconducting phase diagram of MATBG with  ~ 1.1°.(a) the upper critical field, Bc2,(16 mK), and deduced ab(16 mK) by Eq. 1.(b) deduced ab(16 mK) and c for four doping states for which Ic(sf, 16mK) was reported by Lu et al. [11].(c) Cooper pairs surface density,  ,, , and the ratio of  ,,   for four doping states for which Ic(sf, An interesting result is that for all four different superconducting domes, the ratio of: significantly over the whole phase diagram.What is more surprizing is that two superconducting domes located near nn = -1.73 10 12 cm -2 and nn = 1.11 10 12 cm -2 have practically the same  ,,   ≅ 0.015, and more or less close values for ab(T=16mK) and c(T=16mK).
superconductors, to test if experimental data necessitates the development of a new phenomenology.We investigate s-, d-, p-and d+id pairing symmetry scenarios within single and two-band models.As a result, we found that phonon mediated d-wave, d+id-wave and three of the four p-wave pairing symmetries should be excluded from further consideration.More experimental data is required to confirm which of the remaining two ( i.e., s-wave or axial Al p-wave) pairing scenarios this material exhibits.As graphene has a planar honeycomb lattice of sp 2 bonded carbon atoms, it is unsurprising that phonon-mediated superconductivity would exhibit pairing symmetry of either s-or p-wave in MATBG.

Table 3 .
, which are different by an order of magnitude.This is a reasonable result given the interlayer charge carrier concentration is much lower in comparison with interlayer one, as the two SGL are only interacting via weak van der Waals forces.Deduced parameters for two-band models for MATBG sample M2 (Ref.7) doped at nn=-1.44 10 12 cm -2 where for both bands we assumed an effective mass of charge carriers, m * /me=0.1637±0.0154,and c = 35.6 (Eq.36).