Phase-dependent dissipation and supercurrent of a graphene-superconductor ring under microwave irradiation

A junction with two superconductors coupled by a normal metal hosts Andreev bound states whose energy spectrum is phase-dependent and exhibits a minigap, resulting in a periodic supercurrent. Phase-dependent dissipation also appears at finite frequency due to relaxation of Andreev bound states. While dissipation and supercurrent versus phase have previously been measured near thermal equilibrium, their behavior in nonequilibrium is still elusive. By measuring the ac susceptibility of a graphene-superconductor junction under microwave irradiation, we find supercurrent response deviates from adiabatic ac Josephson effect as irradiation frequency is larger than relaxation rate. Notably, when irradiation frequency further increases above the minigap, the dissipation is enhanced at phase 0 where the minigap is largest and dissipation is minimum in equilibrium. We argue that this is evidence of the nonequilibrium distribution function which allows additional level transitions on the same side of the minigap. These results reveal that phase-dependent dissipation is more sensitive than supercurrent to microwave irradiation, and suggest a new method to investigate photon-assisted physics in proximitized superconducting system.

A junction with two superconductors coupled by a normal metal hosts Andreev bound states (ABSs), which shuttle the Cooper pairs from one superconducting bank to the other and whose energy spectrum depends on their phase difference (ϕ). The phase-dependent Andreev levels give rise to supercurrent I s periodic in ϕ, and the measurement of I s (ϕ), or the current-phase relation (CPR), has previously revealed, for example, the singlet/doublet transition in carbon nanotube quantum dots [1,2] and the helical edge states in topological materials [3,4]. If ϕ acquires a time-dependent ac component (δϕ), generated for example by an ac magnetic field in a ring geometry, finite-time relaxation of ABSs towards equilibrium causes delay in the current response and consequently a counter-intuitive dissipation appears [5]. Such dissipation involves two important mechanisms: One is the relaxation of thermally excited ABSs via inelastic scattering; the other is inter-level transitions induced by microwave photons [6][7][8]. In the weakly driven regime where δϕ π, the dissipation shed further light on the properties of the Andreev levels, revealing, for example, protected level crossing in topological junctions [9,10]. Strong driving power can significantly modify the distribution function from thermal equilibrium, activating additional level transitions [11,12]. However, this nonequilibrium state has so far mainly been investigated in dissipationless supercurrent response, finding enhanced critical current [11,12] or modified CPR [13,14]. Here we present the evolution of both the CPR and dissipation under microwave irradiation with different frequency and power, extracted simultaneously from the ac magnetic susceptibility of a phase-biased graphene-superconductor ring in long diffusive regime. The Andreev spectrum is characterized by a minigap whose size relative to the irradiation frequency plays an important role in the junction's response [15,16]. Graphene is chosen since it has lower density of states while maintaining similar meanfree path compared to conventional normal metals [17]. This gives access to similar minigap with a reduced supercurrent and the screening effect, enabling accurate measurement over a complete phase period [6].
The ac susceptibility is measured by coupling a graphene/superconductor ring to a superconducting resonator. Fig. 1(a) shows the scanning electron microscopy (SEM) images of the device. The left image shows a section of the resonator (the meander lines) made by ebeam lithography and sputtered molybdenum-rhenium (MoRe) on the undoped silicon substrate. The boronnitride/graphene/boron-nitride (BN/G/BN) stack is fabricated using exfoliated flakes and is then connected to MoRe via side contacts [18]. The details of the fabrication is given in [19]. The right image is a zoom of the junction region. The junction width is W = 5 µm and the length L = 950 nm. The Ti/Au top-gate covers 330 nm of the total graphene length. In Fig.1(b), the resonator is designed to have C r = 180 pF and L r = 40 nH. The loop at the end of the superconducting lines provides a coupling inductance L c = 355 pH. To maintain sufficiently high quality factor Q, the resonator is coupled to the rf coaxial cables on the dilution refrigerator through a coupling capacitor C cpl = 5.6 pF. At T = 12 mK, the resonance frequency ν r = 60 MHz and Q ∼ 300. The dc flux Φ (or the dc phase ϕ = ϕ 1 -ϕ 2 = 2πΦ/Φ 0 ) is set by the dc magnetic field through the area defined by the ring, and is added to the ac flux δΦ inducing the ac current δi ac . The ac susceptibility is thus defined as χ = δi ac /δΦ [6]. The rf power is heavily attenuated so that δΦ Φ 0 and χ does not depend on the rf power. On resonance, the real (dissipationless) and imaginary (dissipative) part χ and χ are linked to ν r and Q via χ = −2(L r /L 2 c )(δν r /ν r ) and χ = (L r /L 2 c )δ (1/Q) [6]. As Φ is swept, δν r and δ(1/Q) are simultaneously measured using a phase-locked feedback loop which maintains the resonator on resonance [6,20,21]. At sufficiently low frequency 60 MHz where χ /χ 1 (justified in [19]), χ (Φ) ≈ ∂I s /∂Φ = (2π/Φ 0 )(∂I s /∂ϕ) where I s is the supercurrent [6]. Thus, integrating the measured χ (Φ) yields I s (ϕ) of the junction.
We first explore the junction without microwave irradiation. Fig. 1(c) displays the CPRs I s (Φ) at T = 12 mK, 0.9 K and 1.5 K [χ (Φ) in inset]. From the CPRs we extract the critical current I c as well as its first three Fourier coefficients I c,n (n = 1,2,3), where I s (ϕ) ≈ I c,1 sin(ϕ) + I c,2 sin(2ϕ) + I c,3 sin(3ϕ). I c (T ) and |I c,n | are plotted in Fig.1(d). Note I c,2 is negative and |I c,2 | = −I c,2 . At T = 12 mK, |I c,2 | and |I c,3 | have non-negligible values, consistent with the skewed CPR in Fig. 1(c) while at T = 1.5 K only |I c,1 | is dominant, meaning the CPR is sinusoidal. Assuming long diffusive junction model at low temperature, we can fit I c (T ) by where Thouless energy E T h and normal resistance R are two fitting parameters [22]. The fitting yields R = 5.6 kΩ and E T h = 34 µeV (equivalently 400 mK or 8 GHz). The superconducting gap of MoRe is [22], the mean-free path is l e ≈ 100 nm. The minigap E g = 2 × 3.1E T h = 206 µeV (or 46 GHz) [23] is higher than the irradiation frequency accessible in later experiment. However, the above estimation assumes perfect interface and strongly diffusive transport. By numeric simulation [19], we show that E g = 2 × 3.1E T h overestimates the minigap in junction with imperfect contacts and weak disorder, which is plausible in BN-encapsulated graphene sample. The actual E g may thus fall within the energy range probed in the experiment. In [19], we also discuss the possibility of describing the data without irradiation with a ballistic model and find similar levels of agreement . However, the data with irradiation agree better with the diffusive model. Fig.1(e) shows |I c,n | versus V g at T = 12 mK. At V g = 2 V, the supercurrent is almost saturated at 44 nA, smaller than in previous studies [24][25][26] possibly because of insufficiently filtered radiation from environment and the higher normal resistance of graphene in the un-gated regions. Meanwhile, |I c,1 | ∼ 5|I c,2 | at V g = 2 V, similar to a uniformly gated sample [24], meaning the transparency between the MoRe and graphene is still enough to preserve I c,2 . For negative V g , |I c,1 | is much reduced and |I c,2 |, |I c,3 | become negligible, consistent with the formation of pn junctions [24,26]. Throughout the paper we focus on the high electron-doping regime with V g = 2 V for the highest susceptibility signal.
After exploring the junction without irradiation and finding the CPRs agree with previous experiments [24,26,27], we now present the data with irradiation. In Fig. 2(a, b), the CPRs are measured at two irradiation frequencies ν = 2 GHz and ν = 19 GHz, respectively. The power noted in the figure is converted to the normalized power s = αeV pk /hν (V pk is the peak voltage at the source and α includes the attenuation factor from source to junction). In both figures, a sign reversal of the supercurrent occurs as the irradiation power increases. However, for ν = 2 GHz, a strong second harmonic (halved periodicity) is observed during the reversal while almost no phase-dependence is seen for ν = 19 GHz in a similar situation. In Figs. 2(c, d) we plot the power dependence of the first two Fourier coefficients of the CPRs for ν between 2 GHz and 39 GHz. At low irradiation frequency, the junction shows the adiabatic ac Josephson effect [28] where ABSs follow instantaneously the oscillating δϕ and I c,n shows the Bessel function dependence. The experimental α is hard to calibrate accurately, and for each I c,1 (s), α is chosen such that its first zero coincides with s = 1.2 (the first zero of the Bessel function). The same α is then used for I c,2 (s) of the respective frequency. The 2 GHz data agrees well with the Bessel function, indicating the microwave drive is adiabatic [29]. It also shows that the electronic temperature is not significantly heated by irradiation. For higher frequencies, the agreement is less satisfactory. In particular, I c,2 at s = 1.2 decreases for higher ν, which is consistent with the disappearance of halved periodicity in the CPR for ν = 19 GHz. Using time-dependent Usadel equations incorporating a finite inelastic scattering rate γ [29], we calculate I c,1 (s) and I c,2 (s) at low temperature and ν lower and higher than γ. The results are plotted in Figs. 2(e, f). At ν = 0.2γ, I c,1 (s) and I c,2 (s) follow the Bessel function as expected, while for ν > γ they show deviation in qualitative agreement to Figs. 2(c,d). Comparing the theory (1.2γ and 2.5γ curves) with the experiment (5 GHz and 19 GHz curves), the experimental γ can be estimated as between 5/1.2 GHz = 4 GHz and 19/2.5 GHz = 8 GHz. Thus, hγ < ∼ E T h , reasonable for SNS junctions [6]. The finite ν r /γ gives rise to a nonzero dissipation χ [6]. Fig. 3(a) shows χ (Φ) taken simultaneously as the CPRs in Fig. 1(c) without irradiation. χ (Φ) peaks at 0.5Φ 0 and its height decreases with temperature. This higher dissipation is a result of the minigap closing which allows more excitation-relaxation events between Andreev states [6]. At low temperature, δχ /δχ ∼ 20 [here δχ = χ(0.5Φ 0 ) − χ(0)] is indeed similar to ν r /γ estimated independently from CPRs, as discussed in [19]. shows a gradual diminution in δχ for smaller power. This is compared with the CPR, as illustrated in Fig. 3(d): δχ (s = 0.2) is decreased by 80% of the un-irradiated δχ (s = 0) whereas the CPR is almost unaffected, demonstrating the much higher power sensitivity of χ than CPR. At high irradiation power, the dissipation response is drastically different for ν = 19 GHz: Instead of the flat χ (Φ) in Fig. 3(c), an additional lobe appears in Fig. 3(d) around flux 0 and the peak at 0.5Φ 0 becomes pronounced again, which turns χ (Φ) into a function quasi-periodic in Φ 0 /2. To our knowledge, this emergence of an enhanced dissipation at flux 0 under high irradiation frequency and power has never been reported and requires an explanation. While χ (Φ) with nonzero ν r can be calculated by similar Usadel equation approach [8] to that in Fig. 2 for the CPRs in the dc limit, it is not straightforward to include both the low-frequency ac flux and the high-frequency irradiation. In order to explain Fig. 3(c), we thus turn to the Kubo formalism, which has well described the χ (Φ) of diffusive junctions in previous experiments without irradiation [6,30,31]. In our junction, the rf frequency used for susceptibility measurement ν r is much smaller than the inelastic scattering rate γ, thus χ can be approximated by [19]: where J nm = (eh/m e i) n|∇|m is the off-diagonal element of the current operator, n is the n-th Andreev level. The effect of the microwave irradiation is included phenomenologically in the distribution function: (2) where ν is the irradiation frequency and σ stands for the irradiation power.
With the irradiation, the change in distribution function from the equilibrium one δf ( n ) = f n − f F D shown in Fig. 4(a) qualitatively reproduces the results based on the time-dependent Green's function approach near ϕ = 0 [29,32,33]: The Andreev state occupation is depleted (enhanced) within hν below (above) the minigap by irradiation photons. The total χ can thus be rewritten as χ = (1 − σ)χ 0 + σ/2[χ +hν + χ −hν ], where χ 0 and χ ±hν are calculated with Eq.1 using f F D (E) and f F D (E ± hν) respectively. We numerically calculate these terms in a diffusive SNS junction using the tight-binding method [31,34]. See details in [19]. The Andreev spectrum near E = 0 is plotted in Fig. 4(a) showing a minigap. Fig. 4(b) displays two examples |J 1,−1 | 2 and |J 23,44 | 2 whose corresponding (f n − f m )/( n − m ) are nonzero. They represent "interband" transition [red arrow] and "intraband" transition [blue arrow], respectively. As σ goes from 0 to 0.8, χ (ϕ) evolves from a 2π-periodic function to a quasi π-periodic function, which is the key feature in Fig. 3(c). By comparing Figs. 4(c) and (d), the rising lobe at phase 0 is attributed to the nonequilibrium terms χ ±hν and can be understood as follows: At low temperature, the equilibrium f F D is close to a step function, thus only interband terms |J n,−n | 2 contribute to Eq.(6). Since |J n,−n | 2 generally peaks at π (e.g. |J 1,−1 | 2 ) [31], χ 0 (π) is also high. Meanwhile, under irradiation with hν > E g , the distribution function is modified and states near the minigap are partially occupied, enabling more intraband transitions. The intraband terms |J nm | 2 , which in general have higher magnitudes around phase 0 (see |J 23,44 | 2 ), produces an enhanced χ ±hν (0). The same simulation with hν < E g only shows the 2π-periodic χ . The flat χ (Φ) at small power [s = 0.2 in Fig. 3(c)] is not captured by our simple model. However, this may be due to the detailed form of the distribution function deviating from Eq.2.
In conclusion, we measured the ac susceptibility of a phase-biased graphene/superconductor junction in response to the irradiation of microwave photons. For low irradiation frequency, the power dependence of the CPR Fourier coefficients follows the adiabatic response, while the phase dependance of the dissipation almost disappears at small irradiation power, demonstrating its higher power sensitivity than supercurrent. At higher irradiation frequency, the CPR responses deviate from adiabaticity and agree with the quasiclassical theory including finite relaxation rate. More remarkable effects manifest in dissipation as the irradiation frequency further increases above the minigap: The dissipation is enhanced around phase 0 which results in a π-periodic phase dependence. We argue that this new phenomenon happens as the system is driven strongly out of equilibrium. The phase-dependent intraband transitions, which are highly suppressed in equilibrium at low temperature, thus become possible with a nonequilibrium distribution function induced by microwave irradiation. The measurement exemplified here may lead to a more detailed understanding of other nonequilibrium physics in proximitized superconducting system [35][36][37]. It also offers insight to the development of novel superconducting quantum devices operating in microwave field [38][39][40].
Supplementary Materials: Supercurrent and dynamic dissipation of a phase-biased graphene-superconducting junction under microwave irradiation

DEVICE FABRICATION
The superconducting meander lines are defined by ebeam lithography using chemical semi-amplified positive e-beam resist (CSAR) with ∼ 200 nm thickness for high spatial resolution [1], on a high-resistivity SiO 2 /Si substrate. Polycrystalline molybdenum-rhenium (MoRe) alloy is deposited by sputtering in Ar gas with thickness ∼ 60 nm. The graphene and BN flakes are exfoliated and identified via optical contrast on a Si substrate with 290 nm thick SiO 2 . The flakes are picked up by PDMS polymer to form the BN/G/BN stack which is then placed at the end of the MoRe lines [2]. Another e-beam lithography defines the top-gate region which goes across the MoRe lines. The shorting between the top-gate and the MoRe line is prevented by the top BN whose thickness is ∼ 20 nm measured by tapping mode atomic force microscopy. The top-gate is formed by thermally evaporated titanium/gold (5 nm/50 nm). A final e-beam step defines the contacts between the MoRe lines and graphene. After reactive-ion etching through the top BN to expose the graphene edge, the sample is annealed in high vacuum at 80 o C for 1 hour before MoRe sputtering to form the 1d contact to graphene [3].

TIGHT BINDING MODEL FOR THE SNS RING: BALLISTIC VS DIFFUSIVE MODELS
In Bogoliubov-de Gennes formalism, the Hamiltonian of the SNS junction is written in Nambu space as: where H is the system without superconductivity and E F is the Fermi level. The complex pairing potential ∆ = ∆ 0 e iϕ incorporates the phase ϕ of the superconducting electrodes. ∆ 0 is nonzero only in the superconducting regions. Since the data in the experiment is taken with high n-type doping, we approximate graphene lattice by 2d square lattice in the model. Discretizing Eq.1 [4] thus leads to: where i represents the i-th site and N is the total site number, σ x,z are the Pauli matrices operating on Nambu space. i and t ij are onsite potential of the i-th site and hopping energy between i-th and j-th sites respectively. Considering only the nearest-neighbor coupling, t ij = −tδ i,j±1 (δ is kronecker function). The onsite potential i = 4t−E F +V i where V i is the total electrostatic potential in the normal region. We set the Fermi level E F = 4t so that the conduction band is half-filled. In the ballistic case V i is set as 0. In the diffusive case V i is modelled by a uniformly distributed random number within the range [−D, D] (D is the disorder strength). The gating effect can also be included in V i by offseting the potential in the gated region. The Andreev spectrum is calculated by diagonalizing the Hamiltonian for each phases between 0 and 2π and the supercurrent is calculated according to [5]: where n is the n-th Andreev spectrum and f n is the equilibrium Fermi-Dirac distribution function for n . Throughout the simulation we set t = 1 and a = 1. We also set k B = 1 and = 1. The total length of the N region is L = 30, and width W = 4L = 120 (total site number is 3600) except for the 1d ballistic case (W × L = 1 × 30).

Ballistic model
Assuming ballistic transport, the experimental superconducting coherence length ξ = v F /∆ where v F = 10 6 m/s in graphene and ∆ = 1.8k B T c for MoRe (T c ∼ 6 K). Thus ξ ≈ 700 nm ∼ 2/3L (The experimental L = 950 nm). For tight-binding model, since E F = 4 − 2 cos(k F ) is set as 4, v F = ∂E/∂k = 2 sin(k F ) = 2. Choosing ∆ = 0.1 E F , thus ξ = v F /∆ = 20 = 2/3L similar to the relative length scale in the experiment. The Andreev spectra and their corresponding current phase relations (CPRs) are computed for three cases (according to Eq. 3) : (I) 1d ballistic case; (II) 2d ballistic case; (III) 2d ballistic case with doping inhomogeneity. Fig. 1(a) reproduces the 1d ballistic Andreev spectrum with linear levels as previous analytical results [6][7][8]. The level spacing at ϕ = 0 is ∆E/t = 0.124, close to the analytical ∆E/t = π v F /(L + ξ) = 0.125 [6][7][8]. Similar to the main text, the Fourier coefficient amplitudes of CPRs |I c,n |(n = 1, 2, 3) and the critical current I c at each temperature are extracted. |I c,n |(T )(n = 1, 2, 3) is plotted in Fig. 1(b) and fit to the expression [7]: where T c,n = v F /2π(L + 2ξ) and I c0,n is the zero temperature magnitude. The fittings almost exactly agree with the numerics. Also T c,n from the fittings are in good agreement with the analytical T c,n = v F /2π(L + ξ) for all n [7].
In order to see whether diffusive/ballistic transport regimes can be differentiated by the temperature dependence measurement of CPRs, it is illustrative to see how well the ballistic simulations can be fit to the diffusive model. In Fig. 1(c), I c (T ) is fit by the same expression used in the main text [9]: where the Thouless energy E T h and normal resistance R are two fitting parameters. b is a constant 7.7438 set by the ratio ∆/E T h . From the fit, E T h = 3.94 × 10 −3 t, which only accounts for 60% of T c,n and is not relevant to the real energy scales in the 1d ballistic case. Also, since E T h = v F l e /2L 2 , the mean-free path l e = 3.55 ≈ 1/8L, clearly not accurate for the ballistic system.
Extending from the 1d ballistic case to the 2d case by increasing W from 1 to 120 (6ξ), the Andreev spectrum in Fig. 1(d) shows many more levels coming from the modes whose momenta are not perpendicular to the SN interface. Comparing with the 1d spectrum (the blue dashed lines), one finds that the 2d levels are denser around the 1d levels, similar to the quasiclassical calculation [10]. In Figs. 1(e,f), by fitting |I c,n |(T )(n = 1, 2, 3) and I c (T ) to the ballistic and diffusive expression Eq.4 and Eq.5 respectively, T c,n is still close to the theoretical value v F /2π(L+ξ). Thus the characteristic temperature scale in 2d ballistic junction can be estimated by the 1d ballistic model. Meanwhile, E T h resulting from the diffusive fit is again around 60% of T c,1 and l e = 3.55 ≈ 1/10L, meaning the diffusive model is inaccurate in describing the 2d ballistic junction.
Finally, by introducing doping inhomogeneity, the Andreev spectrum is plotted in Fig. 1(c). Compared with the homogeneous case, a "soft gap" (a gap with nonzero but reduced density of state) appears, whose edge is highlighted by the solid black line. Indeed, the doping inhomogeneity can act as elastic scatterer which leads to a diffusive-like Andreev spectrum. Similar features have also been reported in the local density of states measured by the tunneling spectrum in a graphene/superconductor junction [11]. The fitting of |I c,n |(T )(n = 1, 2, 3) by Eq. 4 yields T c,n again close to 1d ballisitc theoretical value. The diffusive model, on the other hand, fits to I c (T ) with better quality and produces E T h similar to T c,n from the ballistic model. Therefore in the inhomogeneous case, both ballistic and diffusive model yield similar characteristic energy scales and do not indicate which one agrees better with the true transport regime. Also, l e = 3.55 ≈ 1/6L close to the length of the gate potential barrier 1/3L.

Diffusive model
The diffusive junction can be modeled by setting a large nonzero disorder strength D = 2.5t [12]. The Andreev spectrum is as Fig. 2(a). The spectrum clearly develops a phase-dependent minigap (highlighted by the black line) which is maximum at ϕ = 0 and closes at ϕ = π, characteristic of the long diffusive junction. From the spectrum, the minigap E g = 16 × 10 −3 . According to the analytical model [5], E g = 2×3.1E T h , thus the Thouless energye E T h = 2.62 × 10 −3 . Since E T h = v F l e /2L 2 , the mean-free path l e = 2.35. Therefore L ∼ 15l e and the long diffusive junction model is indeed applicable. Similar to the ballistic case, |I c,n |(T )(n = 1, 2, 3) and I c (T ) are plotted in Figs. 2(b,c), respectively. The fittings to the ballistic and diffusive expressions Eq.4 and Eq.5 yield similar E T h and T c,n . Similar results are ob-tained by adding doping inhomogeneity to the diffusive system [ Figs. 2(d-f)]. Therefore, although we know a priori the system is in the diffusive regime, the fittings to the supercurrent dependence on temperature alone cannot distinguish which one of the two models, ballistic or diffusive, applies better. In the experiment, E T h /h = 8 GHz and thus the minigap E g = 2 × 3.1E T h = 46 GHz, much higher than the irradiation frequency (19 GHz) with which π-periodic χ (ϕ) is observed. However, the relation E g = 2×3.1E T h is valid only in deeply diffusive regime. Figs. 2(g-i) display the Andreev spectrum, |I c,n |(T )(n = 1, 2, 3), and I c |(T ) in a weakly diffusive junction (disorder strength D = 0.5). The spectrum exhibits a soft minigap with E g = 18 × 10 −3 t (gap edges marked by black solid lines). From I c (T ), the diffusive model fitting yields E T h = 4.72 × 10 −3 t. Thus the minigap estimated by 2×3.1E T h = 30×10 −3 t, almost twice larger than the true minigap known directly from the spectrum. Therefore, the junction in the experiment is likely to operate in the weakly diffusive regime where E T h obtained by I c (T ) fitting overestimate the true minigap, thus E g < 2×3.1E T h . Fig. 3(a) compares directly the normalized experimental CPR with that calculated by tight-binding models, which shows better agreement with the diffusive case. Figs. 3(b,c) plots the same data as in Fig. 1(d) in the main text, fitted by the ballistic model Eq.4 and the diffusive model Eq.5, respectively. Similar to the tightbinding calculation, E T h and T c,n have similar values. The difference between diffusive and ballistic regimes is more clearly differentiated in the supercurrent response to irradiation, calculated by the quasiclassical Green's function approach [13,14]. Figs. 4(a,b) reproduce the Figs. 2(e,f) in the main text by the time-dependent Usadel equations incorporating finite relaxation rate. I c,1 and I c,2 are the Fourier coefficients of the CPR, and s = eV pk /hν is the normalized irradiation power (V pk is the peak voltage of the irradiation). Similarly, the time-dependent Eilenberger equations with finite relaxation can be constructed for the ballistic system [15]. The results are plotted in Fig. 4(c,d). As shown by the tight-binding simulations, the characteristic temperature k B T s = v F /2π(L + ξ) plays a similar role in the ballistic regime to that of E T h in the diffusive regime. The same irradiation frequencies as in Figs. 4(a,b) are used. The calculation conditions listed in the caption are normalized against k B T s . In order to be compared with the diffusive case, the s-axis for I c,1 is rescaled so that all curves have zeros at s = 1.2. The same scaling factor is then used for each I c,2 of respective irradiation frequency. The ballistic results show qualitative difference to diffusive results and the experiments. The calculation assumes perfectly ballistic system and perfect contact transmission. In reality, small disorder and imperfect contact may lead to I c,1 (s) and I c,2 (s) more similar to the diffusive case. From the Kubo formalism, the total χ is [16]:

COMPARISON BETWEEN BALLISTIC AND DIFFUSIVE MODELS WITH EXPERIMENT
where the three terms are called "Josephson term", "diagonal term", and "non-diagonal term" respectively. χ J is purely in-phase with the ac flux modulation and is proportional to the phase derivative of the supercurrent I s . χ D has both real and imaginary part. The finite imaginary part χ D = Im χ D indicates that there is a out-of-phase response thus a dissipation. This diagonal dissipation is caused by relaxation of thermally excited ABSs via inelastic scattering and has distinct phase dependence in which χ D is 0 at ϕ = 0 and π and has two peaks symetric to ϕ = π [17]. This contribution is important at high temperature (k B T ∼ E g ) and is maximum when the frequency of the ac flux ν r is comparable to the inelastic scattering rate γ. In the experiment k B T E g and ν r γ, this term is negligible. Indeed in the measured χ (Φ) without irradiation, no phase dependence expected from χ D is seen. The third term χ N D , resulting from inter-level transitions induced by microwave photons, also has real and imaginary part, and χ N D = Im χ N D corresponds to dissipation. χ N D is the same as defined in the main text and is the dominant contribution to the total χ . We note that total dissipationless χ which is directly proportional to the measured resonance frequency shift δν r is equal to χ J +χ N D , where χ N D is the real part of χ N D . Defining δχ = χ(π) − χ(0), in the limit k B T ∼ hν r hγ E T h , it is shown that δχ N D δχ N D ∼ (hν r /γ)δχ δχ [17] and thus the total χ can be approximated by χ J . From the main text, indeed we see δχ ≈ 20δχ , justifying the integration of the total χ (ϕ) to obtain I s (ϕ). In summary, in the experimental condition, the measured susceptibility can be approximated as χ ≈ χ J and χ ≈ χ N D .
For the tight-binding calculation of χ N D , the system is constructed using a square lattice with W × L = 50 × 30 sites and strong disorder D = 1.5∆. In order to include all terms with significant (f n − f m )/( n − m ) while improving computational efficiency, we cut off the energy at ∼ 0.8∆ (larger than the microwave irradiation energy hν = 0.5∆). The matrix element of the current operator J is computed by [12]: J nm = n|( /i)∇|m = i j,j n|x j , y j x j , y j |∇|x j , y j x j , y j |m = i j Ψ e * n (x j , y j )[Ψ e m (x j + 1, y j ) − Ψ e m (x j , y j )] where Ψ e n (x j , y j ) and Ψ h n (x j , y j ) correspond respectively to the electron and hole components of the wavefunction at site j.

LOW POWER IRRADIATION DATA
In the main text, the lowest irradiation power for both ν = 2 GHz and 19 GHz flattens χ (Φ). Fig. 5 displays another dataset taken at T = 12 mK and V g = 1 V starting from smaller irradiation power. Fig. 5(a) is taken at ν = 10 GHz, and there is a gradual decrease in χ (Φ 0 /2) − χ (0) as P increases and at high power χ (Φ). Fig. 5(b) is taken at ν = 19 GHz. There is also a gradual decrease in χ (Φ 0 /2) − χ (0) as P increases, while at high power χ (Φ) exhibits the halved periodic features same as the main text Fig. 3(c).