Main

In a strongly correlated system, the normal state of a superconductor is subject to diverse instabilities. This can result in unconventional superconductivity (SC) breaking a symmetry on top of global gauge symmetry1. Unconventional SCs play a key part in understanding correlated electron physics and can enable fault-tolerant quantum computations9. Despite its importance, it remains one of the biggest open questions in condensed matter physics owing to the intricate interplay between various broken symmetries and their effects on pairing mechanisms. In bulk materials10,11,12,13,14,15,16, studying doping-dependent behaviour relies on comparing different crystals, which can be complicated by varied defect levels. By contrast, two-dimensional materials allow for in situ tunable doping by electrostatic gating, which enables comparing hundreds of samples effectively in a reliable fashion. Especially, unconventional SC has been actively explored in moiré superlattices17,18,19,20,21,22,23 because of flexibility in interlayer twist angles. Nevertheless, the requirement of a precise twist angle and its spatial variation24 have posed challenges to reproducing phase diagrams from different devices25.

Rhombohedral N-layer graphene (RNG) has recently emerged as an ideal platform to study unconventional SC. It features a simple chemical ingredient and homogeneous crystalline structure, facilitating the search of unconventional SCs vulnerable to impurities3. Moreover, it hosts a flat low-energy band dispersion (E ∝ kN) that can be further fine-tuned by gate voltages26. These features foster enhanced electron correlation effects and various isospin-symmetry-broken states4,5, which can lead to SCs mediated by mechanisms other than phonons27,28,29. Although SCs in rhombohedral bilayer graphene (R2G)6,30,31 and trilayer graphene (R3G)32,33,34 have been studied, unconventional SCs in thicker RNG remain less explored. The latter, however, offers several advantages over its thinner counterparts. First, RNGs with N > 3 provide richer isospin-symmetry-broken states, such as the layer-antiferromagnet4,5 and valley-polarized half-metal35, in which exotic SCs can emerge. Second, isospin-symmetry-broken parent states are induced by lower gate electric fields for N > 3, reducing the risk of gate leakage. Third, RNGs with N > 3 have shown a rich family of quantum anomalous Hall (QAH) states8,36,37,38 that can lead to non-Abelian quasiparticles when combined with SC7,39,40.

Here, we report transport measurements on hole-doped rhombohedral tetralayer (R4G) and pentalayer graphene (R5G) without proximitized spin–orbit coupling (SOC) and moiré effects. We observed multiple superconducting states, highlighting several unconventional ones in R5G that are enhanced or induced by magnetic fields: (1) SC2 is strengthened by an in-plane magnetic field B∥; (2) SC4 is induced by a high B∥; (3) SC3 is enhanced by a small out-of-plane magnetic field B⊥ and is robust against B∥. These observations imply their unconventional nature. From a separate R4G/WSe2 device, we observed several new SCs induced by proximitized SOC effects. These states are located at relatively low displacement fields, and their parent states preserve the long mean free path of bare graphene—both facts will facilitate unconventional SC through interface engineering. The main findings in this work were reproduced in two R4G and five R5G devices and in two independent fridges from two different labs.

SC in rhombohedral multilayer graphene

Figure 1a,b shows the Rxx maps of R5G and R4G at the mixing chamber temperature T of 7 mK, highlighting several regions with vanishing resistance. We label each state as SC1–SC3. Figure 1c,d shows Rxx linecuts at varied temperatures. Rxx drops to zero as T decreases in a range of charge densities. The insets show dVxx/dI as a function of source–drain current I and B⊥. Peaks in dVxx/dI are visible at zero B⊥ and get suppressed by nonzero B⊥. Fermiology analyses (Methods) demonstrate that SC1 is developed at the phase boundary between a full-metal with an annular Fermi surface and a partially isospin-polarized phase (Fig. 1e,f), identically to SC1 in R3G32. Temperature-dependent Rxx, dVxx/dI and fermiology analyses for SC2 are shown in Extended Data Figs. 1 and 2.

Fig. 1: SC in hole-doped bare RNG.
Fig. 1: SC in hole-doped bare RNG.
Full size image

a,b, Rxx as a function of n and D/ε0 for devices R5G-1 (a) and R4G-1 (b), respectively, at zero magnetic field. The inset of a shows a magnified map corresponding to the black box. Several regions with vanishing Rxx are observed, including the SC1 akin to that observed in R3G and three new SC states labelled as SC2 and SC3. c,d, Rxx as a function of n at varied temperatures and fixed D/ε0 for SC1 in devices R5G-1 (c) and R4G-1 (d), respectively. The inset shows differential resistance as a function of d.c. current and out-of-plane magnetic field in R5G (n = −3.315 × 1012 cm−2, D/ε0 = 0.20 V nm−1) and R4G (n = −2.58 × 1012 cm−2, D/ε0 = 0.28 V nm−1). Vanishing Rxx, nonlinear dVxx/dI and the modulation of Ic by B⊥ suggest that these are superconducting states. e,f, Fast Fourier transform of Rxx(1/B⊥) as a function of n and the normalized frequency fν along the dashed lines in a and b, respectively. The relation f1 − f2 = 1/4 holds for SC1 in both R5G and R4G (grey range in the bottom panel), suggesting an annular-shaped full-metal as the parent state. The partially isospin-polarized phase (PIP) appears next to SC1. FFT, fast Fourier transform.

The decrease in Rxx, nonlinear dVxx/dI and response to B⊥ indicate that SC1 and SC2 in R5G and R4G are superconductors. Their phenomenology resembles SCs in moiré graphene17,18,19,20,21, twisted WSe2 (refs. 22,23) and thinner rhombohedral graphene6,30,31,32,33,34. However, the coherence length ξ and mean free path l of these SCs and their normal state (Methods) put themselves in the clean limit (ξ/l ≪ 1), in that ξ ≈ 200 nm, l ≈ 1.6 μm and ξ ≈ 300 nm, l ≈ 1 μm for SC1 in R5G and R4G, respectively. This is contrary to twisted graphene and WSe2 superconductors, in which the normal-state resistance is higher and ξ/l > 1. Rhombohedral graphene in the clean limit provides an ideal playground for exploring unconventional SCs, which are discussed in the following sections.

SC2 enhanced and SC4 induced by B ∥

Having SCs established in R5G and R4G, we investigate R5G under B∥. Figure 2a shows Rxx maps at B∥ = 0 and 8 T. We observe two differences comparing these maps: (1) SC2 is expanded in the n–D space at B∥ = 8 T; (2) a new low-resistance state indicated as SC4 is induced by high B∥. SC4 is well separated from SC2 and SC3 in the map, suggesting it is distinct from the superconductors existing at B∥ = 0 T. SC3 also survives at B∥ = 8 T and is eventually connected to SC2 at B∥ = 8.5 T in the other two devices (Fig. 3d and Extended Data Fig. 4g).

Fig. 2: SC2 and SC4 in R5G, respectively, enhanced and induced by B∥.
Fig. 2: SC2 and SC4 in R5G, respectively, enhanced and induced by B∥.
Full size image

a, Rxx maps at B∥ = 0 T (top) and 8 T (bottom). Whereas SC3 survives, SC2 is enlarged, and SC4 is induced by B∥. For comparison, the boundary of SC2 at B∥ = 0 T is drawn as white dashed curves. b, Rxx as functions of T at the (n, D) point (indicated as circles in a) with the minimal resistance at each field. The transition temperature is increased by applying B∥. The inset shows differential resistance as a function of I and B⊥ for SC2 at B∥ = 6.5 T. Fraunhofer-like modulations are observed. c, Fast Fourier transform of Rxx(1/B⊥) as a function of n and fν. Around SC2, f2 and f3 have merged and the relation f1 − 2f3 ≈ 1/2 holds, suggesting its normal state is a half-metal with an annular-shaped Fermi surface. d, Rxy as a function of B⊥ in the forward and backward sweeping directions in the normal state of SC2 (n = −1.07 × 1012 cm−2, D/ε0 = 58 mV nm−1, B∥ = 0 T, T = 160 mK). The absence of the anomalous Hall effect indicates zero valley polarization. e, Rxx maps magnifying into SC4 at B∥ = 8.0 T (top) and 8.835 T (bottom). The resistance of SC4 decreases as applying higher B∥. f, Rxx as a function of n at varied temperatures and D/ε0 = 140.1 mV nm−1 for SC4 at B∥ = 8.835 T. A clear decrease in Rxx is observed as T decreases. The inset shows differential resistance as a function of I and B⊥ for SC4 at B∥ = 8.835 T and the n value marked by the dashed arrow. Fraunhofer-like modulations are observed. All the data were taken from device R5G-2. FFT, fast Fourier transform.

Fig. 3: SC3 in R5G enhanced by B⊥.
Fig. 3: SC3 in R5G enhanced by B⊥.
Full size image

a, Rxx as a function of n at D/ε0 = 49.61 mV nm−1 and varied temperatures. The resistance vanishes as T decreases. b, Fast Fourier transform of Rxx(1/B⊥) as a function of n and fν, suggesting SC3 is developed from a half-metal with an annular-shaped Fermi surface. c, Rxy as a function of B⊥ in the forward and backward scanning directions in the normal state of SC3 (T = 160 mK). The absence of the anomalous Hall effect indicates zero valley polarization. d, Rxx map at B⊥ = 0 mT and B∥ = 8.5 mT. SC3 features ultrahigh PVR of >34. e, Rxx maps at B⊥ = 0 mT (top) and 1.8 mT (bottom). The SC3 region is enlarged by B⊥. f, Differential resistance as a function of I and B⊥ for n = −1.135 × 1012 cm−2 and D/ε0 = 57 mV nm−1. Ic is increased by small B⊥. g, Differential resistance as a function of I and T for n = −1.135 × 1012 cm−2 and D/ε0 = 57 mV nm−1 under B⊥ = 0 mT (top) and 1.8 mT (bottom). Both Ic and transition temperature increased with the small B⊥. Data in e–g show a rare enhancement of SC by B⊥. h, Description of a theoretical model. When B⊥ much larger than B∥ is applied, the spin vectors are steered along the z-axis, resulting in the valley-contrasting energy ±λ⟨σz⟩/2. This induces the valley imbalance and enhanced valley DOS, which can strengthen intravalley pairing. i, Critical currents (extracted from peaks in dVxx/dI from Extended Data Fig. 9f–k) as functions of B⊥ for B∥ tuned around a few mT. The enhancement of SC3 by B⊥ gets weaker under B∥ of only a few mT. The data in b–g are from device R5G-1 and the data in a,i are from device R5G-2. FFT, fast Fourier transform.

We now focus on the comparison of SC2 at B∥ = 0 and 8 T. First, the area of SC2 in the n–D space increased at B∥ = 8 T. This can be seen in the contour of SC2 at B∥ = 0 T, indicated as a white dashed curve, being smaller than SC2 at B∥ = 8 T. Second, Rxx dropped to lower values across the SC2 at B∥ = 8 T. Finally, Rxx starts to drop from the normal state at higher temperatures under B∥ = 8 T (Fig. 2b). Here, we compared the evolution of Rxx under varied temperatures at the (n, D) point featuring the minimal resistance at each field, because SC2 shifts to a different location in the n–D space by B∥. The (n, D) points where the data in Fig. 2b were taken are denoted by yellow and green circles in Fig. 2a.

The expansion of the phase in the n–D space, with the increase in the transition temperature and critical current Ic (Extended Data Fig. 1e) suggests that SC2 is enhanced by B∥. We performed fermiology analyses to understand its isospin structure. Figure 2c shows the fast Fourier transform of Rxx(1/B⊥) as a function of n and fν (Methods). Two branches f2 and f3 obey the relation f1 − (f2 + f3) ≈ 1/2, and they are merged around the emergence of SC2. This suggests the normal state of SC2 is an annular half-metal. Combining this with the absence of anomalous Hall effect in its normal state (Fig. 2d), we conclude that SC2 arises from a valley-unpolarized annular half-metal parent state.

Contrary to SC2 existing at zero magnetic field, SC4 is observed only at high B∥. Figure 2e shows the Rxx map around SC4 at B∥ = 8.0 T and 8.835 T, featuring a stripe-shaped region with low resistance. The resistance decreases as applying higher B∥, which was limited by the magnet used. Figure 2f shows Rxx linecuts at varied temperatures with Vxx–I characteristics. Nonlinear dVxx/dI is observed at the base temperature, and the reduction of Rxx disappears as temperature increases. Peaks in dVxx/dI disappear at B⊥ ≈ 6 mT.

The suppression of Rxx, nonlinear dVxx/dI and responses to B⊥ for SC4 are consistent with SC phenomenology. Although the fragility of SC4 complicates detailed analyses, we may take 70 mK at which Rxx starts decreasing as a conservative estimate of the transition temperature. This gives the Pauli-limit violation ratio (PVR) of about 68 at the highest field available, when compared with a zero-field Bardeen–Cooper–Schrieffer superconductor with the same transition temperature.

SC3 enhanced by B ⊥

Figure 3a shows Rxx for SC3, featuring vanishing resistance at the base temperature and a sharp rise as T increases. Together with nonlinear dVxx/dI (Fig. 3f), this suggests that SC3 is a superconductor.

We perform fermiology analyses to understand its Fermi surface structure. Figure 3b shows the fast Fourier transform of Rxx(1/B⊥) as a function of n and fν. The branch f1 spans across the whole carrier density, and the low-frequency feature is broken into two segments, f2 and f3. Although the low-frequency branch in the density range of SC3 is unclear, f2 and f3 at its two ends follow relations f1 − f2 ≈ 1/2 and f1 − f3 ≈ 1/2, suggesting a half-metal with an annular Fermi surface. Figure 3c shows Rxy as B⊥ is swept back and forth in the normal state of SC3, showing no anomalous Hall effect. SC3 is still observable at B∥ = 8.5 T (Fig. 3d).

These observations suggest that SC3 is developed from a valley-unpolarized half-metal with annular Fermi surface. Especially, SC3 exceeds the Pauli limit (about 0.24 T) by a factor of >34 (see Extended Data Fig. 7j for the estimation of the Berezinskii–Kosterlitz–Thouless (BKT) transition temperature TBKT). This ultrahigh PVR supports the unconventional nature of SC3.

Figure 3e shows Rxx maps at B⊥ = 0 and 1.8 mT. Unlike SC2 and most other SCs, SC3 is expanded by B⊥. Figure 3f shows dVxx/dI as a function of I and B⊥. The critical current is doubled by B⊥ = 1.8 mT before decreasing at higher fields. TBKT also increases from about 47 mK to 76 mK at this (n, D) point, as shown in comparison between maps at the two field values (Fig. 3g). Similar enhancements of Ic and TBKT are observed across the SC3 region (Extended Data Fig. 7a–f) and disappear under nonzero B∥ (Extended Data Fig. 8a,b).

The enhancement of SC by small B⊥ adds more exoticness to SC3. As explained in the Methods, mechanisms based on spin or valley magnetism, or natural Josephson junctions with the phase difference of π (that is, 0–π junction) cannot explain this observation. We sketch one possible scenario that consistently explains all of our observations (see Methods for details). Under zero magnetic fields, the spin-polarized and valley-unpolarized normal state possesses the spin vectors canted with the angle φ ≈ 90° (ref. 33) (Fig. 3h), because Hund’s coupling EH (about 2 meV) is much stronger than the intrinsic SOC λ of graphene41,42 (approximately 50 μeV). For nonzero B⊥ with B∥ = 0, the spin vectors can be steered along the z-axis even when the spin Zeeman energy (about 0.2 μeV for 1.8 mT) is much smaller than λ, because of the interplay with Hund’s coupling. The nonzero spin component ⟨σz⟩ leads to the valley-contrasting energy ±λ⟨σz⟩/2 and the valley imbalance in the hole occupation, as described in Fig. 3h. The detailed model calculations are presented in the Methods and Extended Data Fig. 9. This will result in an increase in the density of states (DOS) for one valley and a decrease for the other valley. Assuming intravalley pairing dominates SC3, the increased DOS for one valley can strengthen the superconductivity.

Transport measurements under carefully adjusted B⊥ and B∥ provide further evidence for the proposed scenario, although more experiments are required to fully test it. The spin vectors aligned to the z-axis by B⊥ will tend to be steered back to their easy plane (φ = 90°) when B∥ comparable to B⊥ is applied. As B⊥ of approximately 1.8 mT maximizes the strength of SC3, B∥ around the same value will weaken the enhancement. This hypothesis is confirmed in Fig. 3i, in which the maximum of Ic at B⊥ ≈ 1.8 mT is suppressed by B∥ of merely a few mT. Moreover, when stronger SOC is imposed, the spin vectors will be pinned more robustly and accordingly the same B⊥ will induce a smaller spin polarization and valley imbalance. This suggests the enhancement of SC3 is suppressed when strong SOC is proximitized to R5G. This prediction is consistent with the absence of the increase in Ic by B⊥ for SC3 in R5G/WS2 (Extended Data Fig. 9o).

New SCs induced by SOC

Finally, we examine the impact of proximitized SOC on SC in rhombohedral graphene. Figure 4a shows the Rxx map of device R4G/WSe2. Compared with bare R4G (Fig. 1b), we observe additional states with vanishing Rxx labelled as SC3–7. Figure 4b–g shows dVxx/dI as a function of I and B⊥ for these states. For each SC, dVxx/dI vanishes at small I and exhibits peaks at zero B⊥. These features are suppressed as B⊥ increases. Extended Data Fig. 10 shows temperature-dependent Rxx and dVxx/dI, which are aligned with SC phenomenology.

Fig. 4: Proliferation of SC in R4G by proximitized SOC effects.
Fig. 4: Proliferation of SC in R4G by proximitized SOC effects.
Full size image

a, Rxx map of a R4G/WSe2 device. Five new states (SC3–7) emerge while SC2 is suppressed. The names of SCs follow the convention of Fig. 1b. b–g, Differential resistance as a function of I and B⊥ for SC3 (n = −2.356 × 1012 cm−2, D/ε0 = 0.18 V nm−1), SC4 (n = −1.482 × 1012 cm−2, D/ε0 = 0.072 V nm−1), SC5 (n = −1.242 × 1012 cm−2, D/ε0 = 0.072 V nm−1), SC6 (n = −1.00 × 1012 cm−2, D/ε0 = 0.026 V nm−1), SC7 (n = −0.38 × 1012 cm−2, D/ε0 = − 0.24 V nm−1) and SC1 (n = −3.25 × 1012 cm−2, D/ε0 = − 0.42 V nm−1). Vanishing Rxx, nonlinear dVxx/dI (Extended Data Fig. 10) and their suppression by B⊥ suggest SC3–7 are superconductors. h, Square resistance as a function of n in the normal state of the SC1 of bare R4G (top, D/ε0 = 0.42 V nm−1), R4G/WSe2 at negative D (middle, D/ε0 = − 0.42 V nm−1) and R4G/WSe2 at positive D (bottom, D/ε0 = 0.469 V nm−1). Little change in Rs implies the preserved high quality of graphene with WSe2 at proximity.

These observations indicate SC3–7 are superconductors. Phenomenologically, the introduction of SOC proliferates regions with different isospin symmetries and renders more SCs at their phase boundaries. Simultaneously, SC1 is weakened when the holes are polarized to the WSe2-proximal side of R4G, evidenced by decreases in Tc, Ic and Bc,⊥ (Extended Data Fig. 10). Together with SC2 missing in R4G/WSe2, the diverse impacts of SOC on SC in R4G are similar to observations in R3G34.

We note that engineering SC states by WSe2 proximity keeps the high quality of bare graphene. Figure 4h compares normal-state resistance for SC1 in bare R4G, the WSe2-distant and WSe2-proximal side of R4G/WSe2 at comparable D fields. Although the strength of SC1 varies, all the normal-state square resistance is about 25 Ω/□ for the hole density where SC1 emerges. This suggests that the mean free path and the cleanliness of graphene are mostly unaffected by the WSe2 layer.

Discussions and outlook

In summary, we observed lots of superconductors in R4G and R5G, featuring three unconventional ones (SC2–4 in R5G) with unusual responses to magnetic fields. Especially, their remarkably high PVRs put them among superconductors with the highest PVR value reported2,43. These new SCs show several distinct characteristics compared with SCs in the literature. First, in-plane critical field exceeding the Pauli limit has been reported in transition-metal dichalcogenides (TMDs)43,44,45,46 and SOC-proximitized graphene31,37. There, strong Ising SOC pins the electron spins to the out-of-plane direction and makes the SC robust against B∥. However, these Ising SCs are killed by high B∥ before being fully spin-polarized. By contrast, SC2–4 are robust against B∥ with ultrahigh PVR and arise from half-metals. As the intrinsic SOC of graphene is much weaker than Hund’s coupling33,41,42 favouring the spin polarization of half-metals, this suggests SC2–4 are spin-polarized differently from Ising SC. Second, notable PVR has been reported in twisted trilayer graphene21 and R3G32 without proximitized SOC. However, their demonstrated PVRs are smaller than the lower bound of SC3 in R5G that is not even enhanced by B∥. Third, field-enhanced and field-induced SCs have been observed in uranium-based compounds47,48, organic superconductors49,50 and EuxSn1−xMo6S8 (ref. 51). Spin fluctuations or the Jaccarino–Peter effect52 has been proposed to explain these phenomena. However, bare graphene lacks magnetic elements and the aforementioned mechanisms cannot be directly applied. Considering the alteration of adjacent states under a high in-plane field, we propose that change in the normal states by B∥ may enhance SC2 and induce SC4 (Methods). This might also be the case for SC in R2G6 appearing only under nonzero B∥. Further efforts are necessary to fully understand these field-enhanced and field-induced SCs.

Rhombohedral multilayer graphene also provides a platform for further engineering SCs. A promising direction is proximations with QAH states to engineer non-Abelian quasiparticles7,39,40 in a lateral junction. QAH states with a wide range of Chern numbers and fractional QAH states have been discovered in rhombohedral graphene neighboured by TMDs8 or with moiré effects36,37,38. The coexistence of SCs and topological states alleviates complexities in combining them from different materials. Moreover, many superconductors in our work can be reached by low gate electric fields, compared with other graphene SCs6,30,31,53,54. This reduces the risk of gate leakage when accessing these SCs, which will facilitate realizing non-Abelian anyons and topological quantum computations.

Methods

Device fabrications

Rhombohedral graphene and hexagonal boron nitride (hBN) flakes were prepared onto SiO2/Si substrates by mechanical exfoliation. The rhombohedral domains of graphene were identified and confirmed with an infrared camera55, near-field infrared nanoscopy56 and Raman spectroscopy57, and subsequently isolated by cutting with a femtosecond laser. Van der Waals heterostructures were made following a dry transfer procedure. We picked up the top hBN, graphite, middle hBN and graphene using poly(bisphenol A carbonate) film on polydimethylsiloxane, and landed it on a prepared bottom stack consisting of an hBN and graphite bottom gate. The device was then etched into a multi-terminal structure using electron-beam lithography and reactive ion etching. Cr/Au was thermally evaporated for electrical connections to the source, drain and gate electrodes.

Electrical transport measurements

The devices were measured in either a Bluefors LD250 dilution refrigerator at MIT or a Leiden MNK126-700 dilution refrigerator at the University of Basel. In the setup at MIT, the d.c. and a.c. currents were generated by Keysight 33210A function generator with an a.c. frequency of 17.77 Hz. Stanford Research Systems SR830 lock-in amplifiers and Basel Precision Instruments (BASPI) SP1004 voltage preamplifiers were used to measure the longitudinal and Hall resistance Rxx and Rxy. Keithley 2400 source meters were used to apply top- and bottom-gate voltages, VTG and VBG. Twisted phosphorus bronze wires or thermocoax cables thermalized at each stage of the fridge went through two thermal meanders (one on the mixing chamber plate and the other on the cold finger), ceramic radiofrequency filters (DigiKey), four-stage RC filters, and silver epoxy filters (MFT25 from BASPI and a custom-made one) to cool down electrons. In the setup at the University of Basel, MFLI Zurich Instrument lock-in amplifiers were used. The same lock-in amplifier generated the a.c. current, whereas the d.c. current was provided by a BASPI DAC SP927. Moreover, the drain–source current was pre-amplified to a voltage signal by an SP983c I to V converter. Thermocoax cables thermalized at each stage of the fridge were used with two sets of silver epoxy filter58 (one on the mixing chamber plate and the other right next to the sample holder). Keithley 2400 or BASPI DAC SP927 was used to apply the gate voltages. The applied VTG and VBG were swept to adjust carrier density n = (CTVTG + CBVBG)/e and displacement field D = (CTVTG − CBVBG)/2, where CT and CB are top-gate and bottom-gate capacitance per unit area obtained from the Landau fan diagram. The longitudinal resistance Rxx is defined as four-terminal differential resistance dVxx/dI at zero d.c. current.

Shubnikov–de Haas oscillations and fermiology analyses

Shubnikov–de Haas oscillations were obtained by measuring Rxx as a function of B⊥ at the base temperature. For fermiology analyses, Rxx was Fourier-transformed as a function of 1/B⊥. The Fourier-transformed data are a function of the frequency f1/B. Then, f1/B is normalized by the frequency that corresponds to the full carrier density as fν = f1/B/(|n|h/e).

Extraction of coherence length and mean free path

The coherence length ξ of a superconductor and the mean free path l of its normal state are extracted as ξ \(=\sqrt{\frac{h/(2e)}{2\pi {B}_{c,\perp }}}\) and \(l\cong \frac{h}{{e}^{2}}\frac{1}{4{k}_{F}{R}_{s}}\), respectively. Here, Bc,⊥ is the out-of-plane critical field, \({R}_{s}={R}_{n}W/L\) is the square resistance of normal-state resistance Rn, \({k}_{{\rm{F}}}\approx \sqrt{\pi n}\) is the Fermi wavevector, L and W are the length and width of the sample, h is Planck’s constant and e is the elementary charge32.

Measurements under high B ∥

Devices R5G-1, R5G-2 and R5G-3 were measured in the fridges at MIT and the University of Basel, and SC1-3 appearing at zero magnetic fields were observed in both of the fridges. The fridge at the University of Basel is equipped with a Cryomagnetics two-axis magnet providing the magnetic fields both in-plane and perpendicular to the sample. Here, the quench boundaries limited the highest available fields. As the perpendicular component of the in-plane field due to a sample misalignment could be compensated with the perpendicular magnet using the sample as a magnetometer, we obtained the data for R5G under B∥ from the fridge at the University of Basel. The data for R4G/WSe2 under B∥ (Supplementary Fig. 1) were obtained from the fridge at MIT using a three-axis superconducting magnet (American Magnetics).

Fraunhofer-like interference patterns and nonzero residual resistance

The residual resistance relative to the normal-state resistance in some of the superconducting states, such as SC2 and SC4 in R5G, is relatively large compared with the others. The ratios of residual resistance to normal-state resistance for superconductors reported in this work are presented in Supplementary Table 1. This nonzero residual resistance has been observed in several graphene SCs30,31,32,59. Fraunhofer-like interference patterns that arise from phase-coherent transport can serve as additional evidence for SC. There can be non-SC islands in a μm-size device because of sample inhomogeneity, even when (n, D) are tuned to a superconducting phase. Phase-coherent charge transport is affected by the Aharonov–Bohm phase, and thus the interference pattern appears. We found these interference patterns and periodic modulations of SC by B⊥ in most of the superconducting states reported here: SC1 (Extended Data Fig. 6c) in R4G; SC1 (Extended Data Figs. 3d and 5f), SC2 (Fig. 2b, Extended Data Figs. 3e and 5g,h), SC3 (Extended Data Figs. 5i and 7i) and SC4 (Fig. 2f and Extended Data Fig. 3f) in R5G; SC1 (Fig. 4g), SC3 (Fig. 4b), SC4 (Fig. 4c) and SC6 (Fig. 4e) in R4G/WSe2. We note that a limited number of crystalline rhombohedral graphene devices32,33 have demonstrated Fraunhofer-like patterns in their superconducting states. This is in contrast to moiré graphene devices17,18,19,20,60,61 in which Fraunhofer-like patterns could be observed more universally and were more pronounced. Possible reasons are (1) less inhomogeneity in a crystalline graphene system and (2) a lower out-of-plane critical field that kills SC even before the interference pattern is developed. For example, SC1 at negative D in R4G/WSe2 shows the first side-peak at B⊥ ≈ 1.1 mT (Fig. 4g). However, SC5 (Fig. 4d) and SC1 at positive D (Extended Data Fig. 10s) are already killed at this value of B⊥ and the interference patterns do not appear.

Ruling out possible mechanisms for the enhancement of SC3 by B ⊥

In summary, SC3 is enhanced by a small out-of-plane magnetic field when B∥ is zero. Similar values of B∥ do not enhance the SC3. The enhancement by B⊥ also disappears in the presence of nonzero B∥. Here, we discuss several scenarios that cannot fully explain these observations. First, B⊥ can couple to orbital ferromagnetic moments and lower the energy of the SC state compared with its competing states. This is the case for chiral SC developed at high displacement fields53. However, unlike chiral SC, the absence of magnetic hysteresis in the normal (Fig. 3c) and superconducting (Extended Data Fig. 8c) states of SC3 suggests the vanishing orbital ferromagnetism. Second, even without a spontaneous valley polarization, the valley-contrasted orbital magnetization can be generated by finite B⊥ through the valley Zeeman effect and induce the valley polarization. The DOS in one of the valleys may increase, and thus SC3 is enhanced. This is, however, contradictory to the disappearance of the enhancement under nonzero B∥, because it is unlikely that the in-plane field affects the valleys. Third, spin Zeeman energy can shift the Fermi level to a van Hove singularity whose enhanced DOS can develop stronger SC. However, the absence of the increase in Ic under B∥ (Extended Data Fig. 8f) rules out this scenario. Finally, for a percolative superconductor, natural junctions formed by superconducting and normal-state islands in a sample can lead to an increase in Ic under tiny B⊥. This occurs when a π phase difference across the junction, which is called a 0–π junction, is compensated by the B⊥-induced Aharonov–Bohm phase. This phenomenon has been observed in engineered Josephson junctions62,63,64,65. The anomalous SC phenomenology under B⊥ observed in some twisted bilayer graphene devices18,60 and twisted MoTe2 (tMoTe2) nanojunctions66 has also been understood as coming from a 0–π junction. For the tMoTe2 nanojunctions, it was proposed that both spin-singlet and spin-triplet condensates exist across the Pd7MoTe2–tMoTe2–Pd7MoTe2 junction (where Pd7MoTe2 is a superconductor without anomalous behaviour), because tMoTe2 breaks the inversion symmetry. Provided these two condensates have a phase difference of π, a 0–π junction is formed across the nanojunction and the SC transport phenomena may be enhanced by applying B⊥. Regarding SC3 in R5G, if there are 0–π junctions in the devices due to sample inhomogeneity, the optimal B⊥ value resulting in the strongest enhancement will notably depend on (n, D). This is because the normal-state area connecting two superconducting regions in the junction will expand as the SC gets weaker. However, the value of B⊥ showing the strongest SC3 is independent of (n, D) in our data (Extended Data Fig. 7a–f).

Proposed mechanism for the enhancement of SC3 by B ⊥ based on SOC

Here, we provide a physical picture compatible with all our observations on SC3 in R5G. The normal state of SC3 is a spin-polarized and valley-unpolarized half-metal. Owing to the intrinsic SOC strength λ much smaller than the Hund’s coupling EH in bare R5G, the spins from the two valleys are largely polarized in-plane with a small canting to the out-of-plane direction33,67. When a large-enough B⊥ is applied and the spins are polarized out of plane, the SOC projected to one layer of R5G (suggested by the fermiology analysis), \({H}_{\mathrm{SOC}}=-\frac{\lambda }{2}{\tau }_{z}{\sigma }_{z}\), where τz is the valley, and σz is the z-component of the spin, induces a valley-contrasting shift in the Fermi energy. Then, the DOS of one valley increases, whereas that of the other valley decreases. Assuming the dominant pairing mechanism is intravalley and the two valleys are mostly decoupled68 owing to the large momentum difference between them, one of the intravalley pairing channels can be enhanced by the increased DOS69. As the global transport of a superconductor is dominated by the stronger pairing channel, SC can be enhanced by B⊥. The SOC-induced shift in the Fermi energy is weaker when the spin is mostly aligned in-plane by B∥, and thus the enhancement by B⊥ is not observed.

We provide calculations supporting the aforementioned scenario. We consider the free energy density F (refs. 33,67,70):

$$\begin{array}{c}F=\frac{\kappa }{2}({n}_{+}^{2}+{n}_{-}^{2})-{J}_{{\rm{H}}}{{\bf{n}}}_{+}{\boldsymbol{\cdot }}{{\bf{n}}}_{{\boldsymbol{-}}}-\frac{\lambda }{2}({n}_{+}^{z}-{n}_{-}^{z})\\ \,-\,{h}_{\perp }({n}_{+}^{z}+{n}_{-}^{z})-{h}_{\Vert }({n}_{+}^{x}+{n}_{-}^{x})-{h}_{V}({n}_{+}-{n}_{-}).\end{array}$$

Here, n+ and n− describe the total spin polarization vectors whose magnitude is the hole density and direction corresponds to the spin direction for each valley (+ for K, − for K′). The first term features the largest energy scale from the band dispersion and long-range Coulomb interaction. The second term describes Hund’s coupling, where JH ≈ 2 meV/(1012 cm−2) (ref. 33). The third term accounts for the SOC of the system. For bare R5G, this term corresponds to the intrinsic Kane–Mele SOC projected to one layer of the graphene. We take λ = 50 μeV for bare R5G33,41,42,71. When the SOC is proximitized to graphene by TMD, λ is of the order of 1 meV depending on the relative angle between graphene and TMD31. The fourth and fifth terms are the spin Zeeman energy for B⊥ and B∥, respectively. Here, h⊥ = gSμBB⊥ and h∥ = gSμBB∥ with the spin g-factor gS = 2 and Bohr magneton μB. The final term accounts for the valley Zeeman effect, where hV = gVμBB⊥. The two remaining parameters, κ and gV, are affected by the band structure. We used a 10-band Slonczewski–Weiss–McClure-type tight-binding model and calculated the average energy per hole ⟨Eh⟩ and gV. We assumed a half-metal state with n = −1.1 × 1012 cm−2, D/ε0 = 50 mV nm−1, which corresponds to where the SC3 was observed. Two sets of tight-binding parameters taken from refs. 53,72 gave consistent results of ⟨Eh⟩ ≈ 1.4 meV and gV ≈ 6. The value of ⟨Eh⟩ is used to estimate κ as \(\frac{\kappa }{2}({n}_{0}^{2}+{n}_{0}^{2})=\frac{\kappa {n}^{2}}{4}=\langle {E}_{{\rm{h}}}\rangle |n|\) and thus \(\kappa =4\langle {E}_{{\rm{h}}}\rangle /|n|\) ≈ 5 meV/(1012 cm−2).

We can calculate the spin canting angle (Fig. 3h) and the valley imbalance of hole occupations by minimizing F as a function of n±. We write the spin polarization vectors as

$${{\bf{n}}}_{+}=({n}_{0}+\delta )(\sin {\theta }_{+},0,\cos {\theta }_{+}),\,{{\bf{n}}}_{-}=({n}_{0}-\delta )(\sin {\theta }_{-},0,\cos {\theta }_{-}).$$

Here, θ+ and θ− are the spin canting angles and δ is the valley imbalance. The total hole density is |n| = 2n0. For zero magnetic fields, the solution is \({\theta }_{+}=\pi -{\theta }_{-}={\cos }^{-1}\left(\frac{\lambda }{2{J}_{{\rm{H}}}|n|}\right)\) and δ = 0 (ref. 33). When external magnetic fields are applied, we cannot assume θ+ = π − θ− and δ = 0, but can numerically find the solution.

We first focus on the case when B∥ = 0 but B⊥ is nonzero. The calculated spin polarization along the z-axis, \({s}_{z}=({n}_{+}^{z}+{n}_{-}^{z})/|n|\), and the valley imbalance δ normalized by n0 are shown in Extended Data Fig. 9a,b. The spin polarization is saturated to unity at around B⊥ = 1.1 mT, regardless of the valley Zeeman effect. This field value is reasonably close to the optimal B⊥ (1.4–1.8 mT) found in experiments. The valley imbalance also rapidly increases up to the same field value and shows the cusp regardless of gV. This increase leads to the substantial change in the DOS of each valley and thus the enhancement of SC3 as explained previously. When the valley Zeeman term is turned on, the valley imbalance keeps increasing slowly above the optimal field, although sz is saturated to unity. Nevertheless, this increase in δ is much slower than that driven by the spin canting, and the orbital depairing effect will probably become comparable to or more important than the valley Zeeman effect in experiments.

We can expect that the spins are steered ineffectively when B∥ comparable to B⊥ is applied. This can be directly shown in the calculation, including the spin Zeeman effect from B∥, as shown in Extended Data Fig. 9c. The increase in the valley imbalance by B⊥ becomes much slower under B∥ of a few mT. This is manifested in the SC3 phenomenology as the B⊥-induced enhancement disappearing by a tiny B∥ of a few mT, which is observed in Fig. 3i and Extended Data Fig. 9f–k. We note that there has been no example showing that graphene SCs are affected by this tiny B∥, and thus this peculiar observation well supports the validity of our model.

The theory can also be checked by increasing the SOC strength through the proximation to a TMD. This is because the spin vectors are pinned more robustly when stronger SOC is imposed, as the Ising SOC serves as a valley-contrasting effective magnetic field. This is explicitly shown in Extended Data Fig. 9d, where sz was computed with λ = 0.05 and 1.0 meV. The spin polarization decreases substantially in the strong-SOC case. Accordingly, the valley imbalance arising from the valley-contrasting energy ±λsz/2 is suppressed. This is demonstrated in Extended Data Fig. 9e, where the valley imbalance at the optimal B⊥ is reduced by an order of magnitude (from 0.66% to 0.05%).

To test the scenario with varied SOC strength, we characterized device R5G/WS2 shown in Extended Data Fig. 9l. SC3 was weak (Extended Data Fig. 9n), but nonlinear Vxx–I relations and peaks in dVxx/dI were observed. Extended Data Fig. 9o shows the evolution of dVxx/dI under B⊥ for SC3, featuring no enhancement behaviour. This is consistent with the prediction from our theory as discussed above.

Effects of B ∥ on a superconducting state

The in-plane magnetic field affects a superconducting state in two ways. First, it affects the free-energy difference between the SC and normal state, ΔE, depending on the spin configuration of a Cooper pair. For a spin-singlet Cooper pair, ΔE decreases as B∥ is applied due to spin Zeeman energy. For a spin-triplet Cooper pair, ΔE is mostly unaffected by B∥. Second, B∥ modifies the details of a parent state. For instance, both superconducting and normal states move in the n–D space by high B∥ (Fig. 2a). The modified parent state can give rise to a larger ΔE and thus a more stable superconductivity. As noted, when a Cooper pair is spin-triplet, the superconducting state does not gain energy from the normal state by B∥. In a naïve picture, therefore, B∥ will not enhance SC by the first effect. This leads to the hypothesis that the second effect (that is, modification of normal state) may result in the enhancement of SC2 and induction of SC4 by B∥. Further efforts are needed to separate the two coexisting effects here and understand the mechanism of the unusual responses of SC2 and SC4 under B∥.

Ruling out other mechanisms for nonlinear dV xx/dI

Nonlinear Vxx–I characteristics, one of the manifestations of superconductivity, can result from other mechanisms. We discuss how to exclude such possible mechanisms other than SC. First, nonlinear dVxx/dI can occur by the Joule heating effect from the drain–source current73. In this case, dVxx/dI monotonically increases by increasing the d.c. current, because the Joule heating power increases with the d.c. current. Our dVxx/dI data show peaks at a certain threshold current, which is inconsistent with this scenario. Second, nonlinear dVxx/dI can arise from a Schwinger-like effect74,75,76,77. The analogy of the Schwinger effect in quantum electrodynamics can come from the displacement of the Fermi surface by a large drain–source current. However, this current-induced displacement does not disappear by a small B⊥, which is inconsistent with our observations where nonlinear dVxx/dI is killed by a small B⊥. The same logic has been applied for distinguishing two different origins of nonlinear dVxx/dI in moiré graphene superconductors76,77.

Recent reports on unusual SCs in rhombohedral graphene

During the preparation and after the submission of this work, there have been several recent reports on unusual SCs in rhombohedral graphene70,78,79,80,81,82,83.