electron-phonon coupling and the anisotropy of the superconducting gap is the first step towards answering these questions. However, the 2D nature of MATBG, the small moiré Brillouin zone, and the relatively low energy scales make the use of many standard techniques for investigating bulk materials, such as calorimetry, angle-resolved photoemission spectroscopy (ARPES) or neutron scattering, challenging or impossible. Similar concerns are relevant for all 2D superconductors, where in many cases the mechanism and character of superconductivity also remain undetermined [11][12][13][14] .

Superconducting mesoscopic devices have proven to be a useful characterization tool of the material they are built of [15][16][17] . In particular, Josephson junctions (JJs) have been used as a probe of electronic thermalization rates 17 , and the superfluid density, through characterizing the kinetic inductance in thin-film devices 16 . In the case of MATBG, superconducting devices have already proven instrumental for probing the charge of the Cooper pairs 18,19 , the long-range coherence of the superconducting condensate 19 and its orbital magnetic properties 20 .

Here, we use a gate-defined Josephson junction (JJ) in MATBG 18,20,21 to extract electron-phonon coupling, thermodynamic, and superfluid properties of MATBG across its phase diagram. Biasing the junction with a combination of DC and AC currents we probe the dynamics of both the electronic quasiparticles and the superfluid of MATBG. We show that the measured timescales governing the junction's transition between resistive and superconducting states are directly related to the microscopic properties of the material, such as electronic cooling power due to phonons, specific heat, and superfluid density. The gatetunability of the device allows us to probe these quantities across the density-tuned phase diagram of MATBG, for chemical potential both within and outside the flat bands. These measurements allow to constrain the electron-phonon coupling of MATBG, and the current bias dependence of the superfluid density is incompatible with isotropic pairing.

Our experimental setup and theoretical models can be applied to other two-dimensional materials, establishing a tool for the study of two-dimensional superconductors at cryogenic temperatures. They advance the state of the art for MATBG, directly probing the superfluid density, rather than spectroscopic 8 properties of the superconducting state. They also allow for the characterization of the thermal relaxation across the whole phase diagram of MATBG and at cryogenic temperatures. Previous studies either accessed only parts of the phase diagram 22 or were limited to higher temperatures 23 , where correlation effects are suppressed 1,2 .

Our device is a JJ electrostatically defined in MATBG, with a twist angle of 1.06 ∘ ± 0.04 ∘ , also studied in reference 18 (Fig. 1a). The global carrier density n, tuned by the back gate, is set to n = -1.73 × 10 -12 cm -2 , at which the bulk has its highest critical current, 250 nA (See SI). Two layers of top gates, separated by a layer of Al 2 O 3 , tune the local density in the central region, allowing us to fine-tune the details of the junction.

For each value of electron density in the central region (Fig. 1b) we analyze the current-voltage (I/V) characteristic. For densities in the central region close to n j = -2 we observe a gradual onset of resistance above a critical current value, consistent with bulk superconductivity (see also discussion of Fig. 1f below). For all other densities, we universally observe a hysteretic I/V trace with two characteristic voltage jumps ΔV, as shown in Fig. 1c. The two jumps correspond to switching from the superconducting to the resistive state (increasing current bias, blue line) and retrapping back (decreasing current bias, blue dashed). Together with Shapiro step measurements 18 this indicates the formation of a weak superconducting link between the left and right parts of the device, where the weak link region can switch between resistive and superconducting states. From the band structure of MATBG 2 (see inset of Fig. 1(d), for a schematic), the weak link region is expected to be metallic except for a narrow range of voltages placing the chemical potential into the gap between the flat and dispersive bands. Such assessment is consistent with the observation of a positive excess current, I ex , 24,25 in the resistive state of a large portion of the phase diagram (green curve in Fig. 1d). In analogy to conventional superconductors 17 , the dynamic response of such metallic weak links should give access to the dynamics of the electronic quasiparticles and the superconducting condensate in MATBG. We probe the dynamics of our weak links by adding a small AC current component to the DC current flowing through the junction. Sweeping the frequency across three orders of magnitude (0.1-100 MHz), we focus on the changes in the I/V characteristics, as shown in Fig. 1e. At low frequencies, the AC drive brings the two hysteresis branches closer together, which can be understood as follows. The abrupt character of switching and retrapping with DC bias suggests that the junction will undergo a change whenever the total current I DC + I RF (t) reaches the critical value for switching (I sw ) or retrapping (I re ). Consequently, one expects the switching to occur prematurely at I sw -I RF , and the retrapping to occur at a higher DC bias, I re + I RF , reducing the size of the hysteresis loop.

For increasing frequency, the effect of AC bias gradually disappears (Fig. 1e), with a different rate for switching and retrapping. This indicates that both processes, in fact, do not occur instantaneously and are characterized each by a certain rate, which we denote as Γ re and Γ sw . We note that under switching (retrapping) rate we mean the characteristic scale for the switching (retrapping) current dependence on AC bias frequency, with a precise definition of that scale given below. At highest frequencies, the AC drive effect is absent, indicating that neither switching nor retrapping processes are fast enough to occur over one AC drive period. The switching and retrapping rates that can be extracted from Fig. 1e reflect the properties of superconducting MATBG. We now turn to their physical interpretation.

We can first rule out switching and retrapping driven only by the dynamics of the superconducting phase difference across the junction, exemplified by, e.g., the Resistively and Capacitively Shunted Junction (RCSJ) model 24 . In that case, the characteristic frequency is fixed by the Josephson relation to 2eΔV/ℏ. For our weak links it is of the order of 10 GHz, several orders of magnitude larger than the frequencies used in our experiments. The RCSJ model also predicts the switching rate to be smaller than the retrapping one, inconsistent with experimental observations (see additional discussion in Supplementary Information). We therefore conclude that our experimental observations require a mechanism beyond the RCSJ model to explain the switching and retrapping charateristics.

Such an alternative mechanism, for both the retrapping and the hysteresis in metallic weak links is the heating of the electrons in the junction, followed by their thermalization 17,26 . In this case the retrapping branch at I < I sw is characterized by a higher temperature than the switching one due to the Joule heating in the resistive state (Fig. 2a, b). This overheating reduces the weak link critical current for the retrapping branch, leading to a hysteresis. Most importantly, retrapping back into the superconducting state requires the electronic temperature to equilibrate to base temperature, a process, depicted in Fig. 2c, that has been directly demonstrated in superconductor-normal metal-superconductor junctions 17 .

While there are several mechanisms for energy dissipation in graphene, at low temperatures the dominant one is the coupling between electrons and acoustic phonons. In particular, thermalization can occur via diffusion of hot electrons into the leads, emission of blackbody photons or interaction of electrons with acoustic phonons (as the optical ones are frozen out) 27 . The first mechanism is suppressed by the presence of a superconducting gap 28 in the leads in our case, while the second one has been estimated to be negligible in MATBG 29,30 . This suggests that the dominant heat loss mechanism is via coupling to phonons, in agreement with conventional SNS junctions 17,26 .

The above mechanism on its own, however, still implies that switching occurs with the Josephson rate 2eΔV/ℏ, which is inconsistent with our observations, as detailed above. To understand the switching dynamics in our devices we now turn to the case without a central gate voltage, i.e., where the sample is homogeneously superconducting at the optimal density. We observe a frequency-dependent IV characteristic (Fig. 1f), despite the absence of a weak link. Note that there is no hysteresis, ruling out overheating as its origin.

In addition to these observations, it has been shown that a supercurrent can flow in MATBG in narrow superconducting paths separated by normal regions 31 . The normal region thus forms a resistive shunt R bulk coupled in parallel to the superconducting regions (purple shaded path in Fig. 2a). At a non-zero frequency ω, the superfluid impedance is purely inductive due to the inertia of the Cooper pairs (blue shaded mechanism in Fig. 2d) and given by Z sc = jωL kin , with the kinetic inductance L kin / m * n s e 2 , where m * is the effective mass, e the electron charge, and n s is the superfluid density. At frequencies larger than R bulk L kin , the impedance of the superconducting branch becomes higher than the resistance of the normal bulk and the AC current flows through the non-superconducting regions (purple shaded mechanism in Fig. 2d). Intriguingly, L kin in MATBG is expected to be large 19 due to two unique properties: extremely low electron densities, and high effective mass 1 . This explains our observation of a rather low characteristic switching rate in Fig. 1f. The same mechanism applies for MATBG weak links -the kinetic inductance of bulk MATBG is then coupled in series to the junction (Fig. 2a).

Using the ideas outlined above, we construct a model to describe the non-equilibrium dynamics of the Josephson junction. Importantly, this model allows us to relate the observed switching and retrapping rates, Γ sw and Γ re , to the microscopic and thermodynamic properties of MATBG. The dynamics of the current-biased junction are described by: C

We highlight that this description has not been previously used to analyze either the MATBG 18,19,21 or conventional 17 Josephson junctions; in what follows below we show that it allows to describe the Josephson junction data in a self-contained way without invoking results of other measurements [29][30][31] . Equation ( 1) describes a Josephson junction with a phase difference φ, a temperature-dependent critical current I J (T), and a fixed excess current value I ex shunted by resistance R J (Fig. 2(a), dashed box). For results in the main text we assume R J ≪ R bulk , the general case is discussed in Supplementary Information. We note that the form of I J (T) has not been determined experimentally; we assume that it is a decreasing function of temperature with a single characteristic scale T J that can be estimated to be of the order 0.1 K based on the disappearance of interference in superconducting quantum interference devices (SQUIDs) 19 . In the main text, we focus on an empirical model

Á θð1 À T=T J Þ that correctly captures the highfrequency asymptotic behavior of the retrapping current; we provide a discussion of different models and their general properties in the Supplementary Information.

Equation ( 2) describes the evolution of the electronic temperature T with respect to the base temperature. The left-hand side represents the total power dissipated in the link, C el being the electronic heat capacity. On the right-hand side, the first term corresponds to Joule heating, while the second one is the electronic heat loss (G th ) attributed, as discussed above, to electron-phonon interactions. The processes relevant for the description of the Josephson effect occur at T ≈ T J (see Supplementary Information), such that the value of the thermal conductivity G th can be approximated by its value at T = T J . The final equation describes the shunting of the junction by the resistive quasiparticles of bulk MATBG (Fig. 2 (a,d)). The current I sc (t) is the full external current driven through the weak link.

Remarkably, we find that the model defined by Eq. ( 1)-( 3) captures all of the behaviors observed in the experiment. As an example, we consider a highly nonlinear regime where the RF amplitude is larger than the hysteresis I sw À I re . For a range of DC bias values the junction spends part of the AC period in the resistive regime and part of it being superconducting, resulting in a double step in voltage, as shown in Fig. 3a. Such voltage values are the average between the resistive and superconducting voltages weighted by the percentage of the time spent by the junction in each regime. Fig. 3d shows a simulated trace in the same regime, demonstrating remarkable agreement between the model and the experiment. As we increase the frequency of the current bias across the junction we recover the regular hysteresis (Fig. 3a, b, black line). The model captures the evolution of the I/V traces as the bias frequency increases, as is shown in Fig. 3e. Even finer details of the experimental data, discussed in the Supplementary Information are captured by the model. These comparisons confirm that our model accurately describes the dynamics of our junction.

To extract the retrapping and switching rates, Γ re and Γ sw , for a given density from the experimental data, we fit the evolution of the retrapping and switching currents as a function of bias frequency. An analysis of the data, discussed in the Supplementary Information, demonstrates that both currents asymptotically approach a constant high-frequency value as 1/ω. To fit the results at all frequencies, we use the following functions:

. That allows to characterize the corresponding rates (see Fig. 3(c), gray lines). The model described in Eqs.

(1),( 2),(3), reproduces correctly the asymptotic behavior of the switching current, while for the retrapping current the result depends on the particular form of I J (T) (see Supplementary Information). For a fixed density in the junction, we extract the switching and retrapping currents for all frequencies and fit the results. In the example shown in Fig. 3(c), for a density of -4.5 × 10 12 cm -2 , we extract Γ re = 0.52 MHz and Γ sw = 2.75 MHz. Therefore, the weak-link dynamics of our junction gives us access to the quasiparticle thermalization rate and kinetic inductance of MATBG (Fig. 2c, d).

We now provide a physical interpretation of the rates, Γ sw and Γ re that allows us to connect them to the properties of MATBG. We begin with the switching rate Γ sw . From Eq. ( 3) we identify the switching rate as Γ sw = R bulk /L kin ∝ n s (see additional discussion in Supplementary Information). Assuming that the resistance of normal regions R bulk does not strongly depend on T or bias strength, Γ sw -1 ∝ L kin , which allows to probe the superfluid stiffness of MATBG.

Before discussing the thermalization rate of the weak link, we note that for ω ≫ Γ sw the AC part of the current does not reach the junction at all: I sc ≈ I DC . Thus, for Γ re > Γ sw , the kinetic inductance would set the rate for both switching and retrapping. However, as shown in Fig. 4, we have Γ re strictly smaller than Γ sw for all densities (note the different y-axis in Fig. 4a, c), confirming that we can interpret the former as a thermalization rate.

The equation governing thermalization in the device in Eq. ( 2) contains two implicit frequency scales: γ G th C el and k

. Importantly, the hysteresis size for DC driving depends on their ratio γ/k (which is proportional to G th , but independent of C el ), while the retrapping rate Γ re depends on both, allowing in principle, to determine both scales.

These observations allow for a qualitative discussion of the results (Fig. 4(a,b)) across the MATBG phase diagram. The noticeable peaks in ΔI/I s , Fig. 4(b), occur near the band insulator (BI) and charge neutrality points (CNP) and indicate a suppressed thermal conductance G th . This can be attributed to a lower density of states near these points compared to other concentrations, expected from previous experiments 32,33 and theoretical analysis 34 . In contrast, Γ re (Fig. 4(a)), which also depends on C el (Fig. 4(a)) shows weaker features at these concentrations, indicating a simultaneous reduction of G th and C el , again consistent with a suppressed density of states. Remarkably, the minimum of Γ re occurs for densities within the dispersive band. A potential explanation for this behavior is an increased resistance R J to the mismatch between flat-band electrons outside the junction and dispersive ones within it. Deeper into the dispersive band, this effect can be offset by an increased G th .

The quantitative nature of Γ re depends on the particular form of I J (T); for the square-root model introduced above and γ/k < 1/2 we obtain an analytical result for the retrapping rate:

Furthermore, the ratio between DC retrapping current and switching is

(see Supplementary Information). We stress that the observed I r and I s are rather close to one another, which results in γ and k being effectively of the same order of magnitude. For larger values of γ the model predicts an 1/ω 2 dependence of the retrapping current under AC bias; therefore the model should not be applicable for

. However, in the absence of direct measurements of I J (T), we will use this model to estimate C el and G th .

We observe, across the whole density range, three sets of values of Γ re : 0.5 MHz, 1 MHz and 1.5 MHz, corresponding to the chemical potential of the link tuned to the dispersive band, lower flat band and upper flat band, respectively. The change in the hysteresis width, Fig. 4(b) is relatively smaller. Using the analytical formula given above, we can estimate for ΔI/I sw ≈ 0.5 that γ ≈ 0.8-2.3 MHz. Reference 22 , where laser-mediated heating of a MATBG sample allows for the extraction of the same quantity, reports a value of γ ≈ 2 MHz. The fact that such different methods for extracting the thermalization rate agree on the obtained value strengthens both of them as reliable characterization tools.

This result already provides an important insight into the lowtemperature behavior of electron-phonon coupling in MATBG when contrasted with those at higher temperatures. In particular, the cooling rate has been found to be of the order of hundreds of GHz above 5 K with a very weak temperature dependence 23 , attributed to effective moiré umklapp scattering 35 (that is related to folding of the acoustic Fig. 3 | Extraction of the switching and rertapping rates. a I/V traces of the junction at bias frequencies of 0.1 MHz (red) and 100 MHz (black) in the regime where the effective AC amplitude is higher than the hysteresis. ΔI re and ΔI sw highlight the change in retrapping and switching currents, respectively, between the two bias frequencies. b I/V traces of the junction at bias frequencies of 8 MHz (red) and 100 MHz (black) in the regime where the effective AC amplitude is higher than the hysteresis. The mismatch in retrapping current between the red and black curves is probably due to a charge jump (note it is of the order of a few pA). c Switching and retrapping currents as a function of AC bias frequency. d-f Numerical simulations of our device in the same regime as the data shown in (a-c). The grey dashed line in (c) is a fit to the functional forms provided in equations ( 1)-( 3).

phonons by the moiré lattice) explaining the linear-in-temperature resistivity 4,10,35 . The strong difference with our result at T ~TJ ≈ 100 mK suggests a suppression of the cooling rate much stronger than linearin-temperature. This result is consistent with electron-phonon scattering at 100 mK being in the Bloch-Gruneisen regime where umklapp scattering is suppressed 35 and resistivity from electron-phonon scattering should follow a stronger power-law dependence on temperature 36 . This excludes electron-phonon scattering as the origin of linear-in-temperature resistance at low temperatures 37 .

In the case of superconductivity, the most relevant quantity when discussing electron-phonon coupling is the dimensionless coupling constant, which we note here λ. The temperature relaxation rate at low temperatures is related to the strength of the coupling to acoustic phonons [38][39][40] . While this coupling does not take the contribution of optical phonons into account, it is expected to be of the same order of magnitude as the full coupling constant 38 . To obtain an estimate we use a Dirac electron model 4,39,40 , motivated by the theoretical 34 and experimental evidence 32,33 for their presence in MATBG bands, in particular in the vicinity of CNP (consistent with the peak in ΔI/I s discussed above). One

, where s is the acoustic phonon velocity. Using T el ~TJ ~0.1 K from the extinction temperature of SQUID oscillations 19 , s ≈ 20 km/sec (the value for singlelayer graphene 41 is expected to be close to that in MATBG 42 ), k F = ffiffiffiffiffiffi ffi πn p

for n ~1 × 10 -12 cm -2 and γ ~1 MHz we obtain λ ~10 -3 . Several comments are in order regarding this estimate of λ. We begin by stressing this estimate may not be directly comparable to other transport measurements (e.g., resistivity) as our estimates of λ stem from the electron-phonon cooling rate. While the relation between the way electron-phonon coupling enters the cooling rate, the resistivity, and the superconducting pairing is established in the Dirac model 4,39,40 , it is yet to be determined (and may be different) in the MATBG bands. Moreover, the potential inhomogeneity of the twist angle across the sample may reduce the average thermal relaxation rate, since regions that are away from the magic angle are expected to have lower density of states and thus slower thermal relaxation. Additionally, the above estimate is assuming the system is in the Bloch-Gruneisen regime T ≪ c s k F ; it has been however demonstrated 35,43 that the crossover to this regime may be quite different in MATBG than in single-layer graphene and, in particular, it occurs at much lower temperatures. In fact, the strong reduction of the thermal relaxation rate in our experiments with respect to the values at 5 K 23 may present a first demonstration of this crossover occurring in MATBG. We note that this estimate assumes coupling to lowest-energy acoustic phonons only, and does not address optical or any other phonons that are frozen out at the experimental temperatures. Finally, the Dirac model should be applicable only around certain fillings; in our case, a signature of Dirac physics is observed near the CNP in enhanced ΔI/I s indicating less efficient heat relaxation due to lower density of states near a Dirac point. We can further estimate G th and C el taking I J -I exc ~5 nA, ΔV ~20 μV from Fig. 3. The result is G th ~250 fW/K and C el ~5 × 10 19 J/K. From the junction area and n ~1 × 10 12 cm -2 one expects above 10 3 electrons, with the usual Sommerfeld expression

≪1, while in our case this implies k B T E F $ 1, that may be related to large residual entropy of interacting states of MATBG 44,45 . Both G th and C el are much higher than those expected in monolayer graphene 40 , consistent with strongly suppressed bandwidth and elec- in MATBG at T ≈ T J are of the same order as those predicted for monolayer graphene 40 .

We now discuss the switching rate Γ sw (Fig. 4c), related to the superfluid stiffness in the bulk of our MATBG device. Importantly, the AC measurements are still performed at a finite DC bias, thus, our measurements reveal the superfluid density at a finite current bias, n s (I DC ) ≈ n s (I sw ). Since n s (I sw ) is a decreasing function of current, the steep increase in Γ sw at the edges of the lower flat band is explained by the decreasing critical current of the junction (Fig. 1 (d)). On the contrary, the decrease of Γ sw for densities in the top flat and dispersive bands, is unexpected -at such low critical currents n s (I sw ) ≈ n s (0) should be density-independent and large. We suggest that this observation can be explained by the kinetic inductance of proximityinduced superconductivity in the junction region. Being very weak, the proximity-induced superfluid has an extremely large kinetic inductance that is in parallel to the smaller one from the bulk TBG, effectively shunting it. Furthermore, since Γ sw ∝ R bulk , changes in the resistance with concentration can also affect its magnitude. In particular, this provides a plausible explanation of the peak of Γ sw near CNP, where the normal-state resistance peaks.

Let us now return to the densities within the lower flat band, where Γ sw is related to the superfluid density of bulk MATBG. The dependence of Γ sw ∝ n s as a function of I DC , shown in Fig. 5, gives important information about the nature of the superconducting gap in MATBG. Current biasing a superconductor produces a Doppler shift 46,47 of the quasiparticle bands in a superconductor, see inset of Fig. 5 (ΔE = v s × ℏk, where v s is the superfluid velocity and k the quasiparticle momentum). For an isotropic superconductor, depicted in the inset of Fig. 5, this does not affect the quasiparticle occupations until a critical value of bias current is reached. As a result, n s (I) dependence is highly nonlinear with an abrupt drop close to the critical current 48 . For a highly anisotropic or nodal superconductor, across its nodal axis in real space, the quasiparticle band structure presents cones instead of a gap in density (Fig. 5, inset). A small shift originating from a finite bias current, leads to a finite generation of quasiparticle pairs, thus reducing the superfluid density before breaking down the superconducting condensate. As Fig. 5 shows, the relation between superfluid density and bias current is linear in the case of MATBG in the range I dc ∈ [0.6I c , 0.95I c ]. This result is inconsistent with the behavior expected of an isotropic superconducting gap, ruling in favor of a highly anisotropic or nodal pairing state in MATBG.

In conclusion, we have developed a method for characterizing electron dynamics in twisted bilayer graphene by combining electrostatic and radiofrequency current bias of an electrically-defined Josephson junction. This has resulted in a study of out-of-equilibrium dynamical properties of electrons at cryogenic temperatures in MATBG. Our results demonstrate the presence of two distinct characteristic timescales guiding the dynamics of electrons, attributed to the thermalization of electrons and the kinetic inductance of the superconducting condensate, respectively. A phenomenological model, capturing these processes, is found to describe the data well and allows to relate the measurement outcomes to the physical properties of the electrons in the material for a wide range of densities. In particular, we discuss the estimates for electron-phonon coupling (that may have bearing on theories of superconductivity [4][5][6][7] or strange metal 36,37 ), specific heat and superfluid density. The current bias dependence of the deduced superfluid stiffness points towards the superconducting gap of the material being anisotropic. The technique we developed in this work can be applied to a wide range of gate-tunable superconducting 2D materials, introducing a general way to access important thermodynamic quantities, such as specific heat and superfluid stiffness. In addition to being a valuable addition to experimental probes of 2D materials, we demonstrated a controllable driving of a correlated electronic system, opening the path to the realization of out of equilibrium states of electrons.

While writing this manuscript the authors became aware of two works where, by different experimental means, some of the quantities studied in this work are also probed 49,50 .

To fabricate the device, we begin by assembling a so called 'stack' of two dimensional materials. Such flakes are exfoliated mechanically using a micro manipulator, polydimethylsiloxane/polycarbonate stamps, a moving stage and an optical microscope. The stack consists of a bottom layer of few-layer graphite, which is used as a global backgate, 24 nm-thick hexagonal boron nitride (hBN), the magic angle graphene, and hBN again, this time with a thickness of 27 nm. The twisting and stacking are performed following the standard procedures of the field, with a cutting step of the graphene flake 51 . The graphene-cutting step is performed with a tungsten needle with a tip diameter of 2.5 μm. Once the stack is in place we move on to the lithography phase of the fabrication process. Gold contacts are defined by a combination of electron-beam lithography, reactive ion etching and electron-beam evaporation. We use chromium as adhesion layer for the gold (10/70 nm). Top gates are defined by electron-beam lithography and evaporation. The graphene is finally etched to define the mesa. A layer of aluminium oxide is deposited through atomic layer deposition. Finally, we define another layer of local gold gates with a chromium adhesion layer (10/110 nm). As a last remark, the device used for this work is the same one as the one presented in reference 18 .

In this work, there are however two modifications with respect to the setup of the aforementioned reference. The first one is that the AC bias is not sent to the central gate but to one of the leads, using a bias T to be able to send both AC and DC signals to the same contact. The other difference is that, because for the AC measurements we present in this study we need a higher degree of precision than for the ones presented in reference 18 , we must ensure that the AC amplitude reaching the device is neither frequency-dependent nor sample- resistance-dependent. The details of such procedure are given in the following section.

When applying a radiofrequency bias to our junction two different aspects must be taken into account. The first one is the evolution of the amplitude as a function of its frequency. Indeed, having a frequencydependent amplitude reaching our junction would make it impossible to disentangle such effects from the physical mechanisms taking place at the junction level. We performed simulations of our circuit using the software LTspice and obtain an evolution in amplitude of our radiofrequency bias reaching the junction of less than 5%, for frequencies ranging from 100 kHz to 100 MHz. We thus conclude that, in comparison to the experimental results, these effects can be neglected.

Because the junction changes its state during the acquisition of an I/V trace, we must also ensure that the change in resistance triggered by the RF biasing does not significantly affect the AC current flowing through the device. In order to achieve that, we place a 100 kΩ resistor in series between the RF feed line and the device. Like for the previous effect, this leads to a variation in AC amplitude reaching the device of the order of a few percents across the whole resistivity range of the device. We thus can also neglect this effect taking into account the precision of the claims made in our analysis. Supplementary Fig. 1 shows a schematics of the electronics set-up. We neglect in this schematics the capacitances from the MATBG to the gold or graphite gate electrodes because they are smaller than 30 fF 18 . They thus present, at the frequencies at which we bias our sample (up to 100 MHz), an impedance orders of magnitude higher than any other impedance in the set-up.