Topological Superfluid Responses of Superconducting Dirac Semimetals

We demonstrate that topological constraints do not only dictate the geometric part of the superfluid stiffness, but can also govern the total superfluid stiffness. By introducing a general adiabatic approach for superfluid responses, we showcase such a possibility by proving that the stiffness of a superconducting Dirac cone in two dimensions (2D) is proportional to its topological charge. By relying on the emergent Lorentz invariance of Dirac electrons, we unify the superfluid stiffness and quantum capacitance in these systems. Based on this connection, we further predict a topological origin for the quantum capacitance of a Josephson junction where 2D massless Dirac electrons are sandwiched between two conventional superconductors. We show that the topological responses persist upon effecting strain, are resilient against weak disorder, and can be experimentally controlled via a Zeeman field. Remarkably, the nonuniversal topological quantization of the two superfluid responses, yet implies the universal topological quantization of the admittance modulus of the superconducting Dirac system in units of conductance. The quantum admittance effect arises when embedding the superconducting Dirac system in an ac electrical circuit with a frequency tuned at the absorption edge. These findings are in principle experimentally observable in graphene-superconductor hybrids.

We demonstrate that topological constraints do not only dictate the geometric part of the superfluid stiffness, but can also govern the total superfluid stiffness.By introducing a general adiabatic approach for superfluid responses, we showcase such a possibility by proving that the stiffness of a superconducting Dirac cone in two dimensions (2D) is proportional to its topological charge.By relying on the emergent Lorentz invariance of Dirac electrons, we unify the superfluid stiffness and quantum capacitance in these systems.Based on this connection, we further predict a topological origin for the quantum capacitance of a Josephson junction where 2D massless Dirac electrons are sandwiched between two conventional superconductors.We show that the topological responses persist upon effecting strain, are resilient against weak disorder, and can be experimentally controlled via a Zeeman field.Remarkably, the nonuniversal topological quantization of the two superfluid responses, yet implies the universal topological quantization of the admittance modulus of the superconducting Dirac system in units of conductance.The quantum admittance effect arises when embedding the superconducting Dirac system in an ac electrical circuit with a frequency tuned at the absorption edge.These findings are in principle experimentally observable in graphene-superconductor hybrids.
The exploration of the interplay between topology and quantum geometry through superfluid responses of timereversal invariant superconductors (SCs) [1] is currently in the spotlight [2][3][4][5][6][7].Such a pursuit got substantially boosted after experiments in magic angle twisted bilayer graphene (MATBG) provided evidences for a nonvanishing superfluid stiffness despite the almost flat energy dispersion which governs transport [8,9].This quite unusual result was subsequently understood in terms of lower bounds set by topological invariants dictating the MATBG band structure [10][11][12][13][14][15].These bounds impose constraints on the so-called geometric contribution to the stiffness, which becomes relevant when superconductivity is harbored by flat bands.However, no topological constraints have been so far predicted for the total superfluid stiffness.The potential discovery of systems whose total stiffness is equal or proportional to a topological invariant promises to deepen our understanding of topological platforms and set the stage for novel applications.
In this Letter, we show that topological constraints dictate the total superfluid stiffness of superconducting Dirac semimetals (SDSs).Notably, a quantized total superfluid stiffness has already been theoretically predicted for superconducting graphene [16][17][18][19], which constitutes a prototypical SDS in its Dirac regime [20].However, the quantization itself, along with its topological origin and resulting implications have not yet been discussed.After Refs.[16][17][18][19], the superfluid stiffness of a superconducting Dirac cone (D cone ) reads at charge neutrality as: and becomes quantized in units of ∆ ≥ 0, which is the intrinsic or proximity-induced conventional pairing gap felt by graphene [21][22][23][24][25][26][27].In the above, we set the reduced Planck constant ℏ and the electric charge unit e to unity.This quantization is rather puzzling as it cannot be explained by topological bounds of the type predicted for MATBG since, here, the dispersions are strongly non flat, and D cone receives equal contributions from both interband (geometric) and intraband (conventional) parts [4].
Here, we resolve this conundrum by proving that πD cone /∆ is a topological invariant.This result is a consequence of nontrivial topology for a single Dirac cone in the presence of a phase which twists its mass, similar to the Jackiw-Rossi model [28].As we further explain in our companion work [29], the quantization is also a manifestation of one-dimensional chiral anomaly [30,31].These properties, along with the arising Lorentz invariance, also render the quantum capacitance c Q of a Josephson junction sandwiching a Dirac cone topological and equal to: where υ D is the group velocity of the Dirac dispersion.arXiv:2211.07676v4[cond-mat.supr-con]13 Aug 2024 Both effects enjoy topological robustness as long as the chiral symmetry of the Dirac Hamiltonian is preserved.To back this claim, we show that Eqs. ( 1) and ( 2) persist upon adding nonuniform strain which leads to an energy spectrum that consists of relativistic pseudo-Landau levels (pLLs) [32][33][34][35][36][37][38][39].We show that for the pLL flat bands, D and c Q are purely of a geometric origin, and are solely carried by the zeroth pseudo-Landau levels (0pLLs).Hence, we reveal that the topological quantization dictating the total superfluid stiffness of the SDS, further imposes the quantization for the geometric part of the stiffness when the Dirac bands flatten due to strain.
We propose to experimentally detect these SDS effects in superconducting (un)strained graphene hybrids, such as the ones shown in Fig. 1.The quantized contribution of the Dirac electrons can be disentangled by applying a magnetic field which couples to graphene electrons only through a Zeeman energy scale B. The Zeeman energy sets the occupancy of the energy levels of the SC.Hence, by manipulating its strength |B|, one is in a position to controllably add or subtract the contribution of the Dirac part of the band structure.Most importantly, we bring forward that, at charge neutrality, the smoking gun signature of the topological effects predicted here is the observation of the universal topological quantization of the admittance modulus Y mod of the SDS, when it is embedded in an ac electrical circuit.Specifically, when the ac frequency ω is tuned at the absorption edge 2∆, Y mod equals the number of Dirac cones in units of conductance.
We commence our analysis by unifying the superfluid stiffness tensor elements D ij and the Josephson quantum capacitance (JQC) c Q in SDSs.The former quantities define the coefficients which link the charge current J (r) to the spatially dependent gauge invariant vector potential A(r) + ∇ϕ(r)/2, where A(r) denotes the electromagnetic vector potential and ϕ(r) is the superconducting phase.Therefore, we have the defining relations: where i, j = x, y for a two-dimensional system.The JQC can be defined in an analogous fashion, through the excess charge density ρ c (t) induced in a Josephson junction due to the presence of the gauge invariant voltage For a two-terminal Josephson junction as in Fig. 1(b), V (t) and ϕ(t) define the voltage and phase biases imposed across the junction.Therefore, we define c Q as the response coefficient which satisfies the relation: The apparent similarities between Eqs. (3) and (4) imply that the charge and current responses can be unified by here introducing the relativistic three-current given as J µ = D µν A ν − ∂ ν ϕ/2 with µ, ν = 0, 1, 2, metric tensor diag{1, −1, −1}, and D 00 = −c Q .Hence, the superfluid stiffness and the JQC are proportional in systems with an emergent Lorentz invariance.Therefore, for Dirac electrons with a "speed of light" υ D , Lorentz invariance leads to the constraint D xx,yy = υ 2 D c Q , which indeed holds for Eq.(2).Having settled the connection between D and c Q in SDSs, we now proceed by examining how nontrivial topology further imposes their nonuniversal quantization.
Our starting point is the Hamiltonian for a generic two-dimensional time-reversal invariant SC: where τ 1,2,3 denote Pauli matrices which are defined in Nambu space.The latter is spanned by electrons with spin up and momentum p, and, their time-reversed hole partners with spin down and momentum −p.From Eqs. ( 3) and ( 4) we infer that D µν can be identified with the coefficients relating the uniform curent J µ induced by a nonzero spatiotemporally-uniform phase gradient ∂ µ ϕ.In particular, to obtain expressions for the superfluid stiffness tensor elements D ij , we consider small deviations of the superconducting phase away from the value ϕ = 0, i.e., ϕ(r) ≈ (∂ x ϕ)x + (∂ y ϕ)y, with ∂ x,y ϕ being constants.Subsequently, we phase twist the pairing term according to ∆τ 1 → ∆τ 1 e −iϕ(r)τ3 and determine the spatially-uniform current component J i within linear response to the phase gradient ∂ j ϕ.In the same spirit, D 00 is found using linear response theory to a linearlyvarying time-dependent phase of the form ϕ(t) ≈ (∂ t ϕ)t.
By adopting this alternative approach, in Ref. 29 we show that D µν can be expressed as suitable response coefficients of the respective adiabatic Hamiltonian: Ĥ(p, ϕ) = ĥ(p)τ 3 + ∆τ 1 e −iϕτ3 (6) where the superconducting phase ϕ is now viewed as an additional synthetic momentum.Using the above adiabatic Hamiltonian, Ref. 29 provides concrete expressions for the elements D µν .Specifically, for D ij we find [29]: where "Tr" denotes trace over all internal degrees of freedom.For compactness, we employed the shorthand notation ´dP ≡ ´BZ corresponding to the frequency in imaginary time.The momenta are generally defined in a two-dimensional Brillouin zone (BZ), since the SC is considered to be a crystalline material.When evaluating the superfluid stiffness carried by a Dirac cone, we replace the BZ by R 2 [29].
In Eq. ( 7), we also introduced the group velocity υ(p) = ∂ p ĥ(p) of the normal phase Hamiltonian ĥ(p), along with the adiabatic matrix Green function: We note once again that the Zeeman energy B sets the Bogoliubov-Fermi level.The last ingredient of the adiabatic approach is the matrix Berry curvature function: where we suppressed the arguments of the Green functions for notational convenience.In spite of the fact that our adiabatic approach is fully-equivalent to the standard method of evaluating the superfluid stiffness, it is unique in uncovering the possible underlying topological properties of the system, since it is already expressed in the terms of a curvature function in a synthetic space [29].
We now proceed by recovering the result of Ref. 16 shown in Eq. ( 1), and proving that πD cone /∆ is a topological invariant.At charge neutrality, the chemical potential µ is set to be zero and the normal phase Hamiltonian for a given graphene valley λ = ±1 reads as [20]: where σ 1,2,3 denote Pauli matrices acting in the sublattice space spanned by the two interpenetrating triangular lattices of graphene [20].Note that, for graphene, electrons couple to holes of different valleys [40].Nonetheless, by a suitable choice of basis, each valley of graphene is described by a Hamiltonian of the form shown in Eq. ( 5).
The adiabatic Hamiltonian obtained for the λ = +1 valley of the graphene model in Eq. ( 10) takes the form: (11) Notably, the operator Π = σ 3 τ 3 establishes a chiral symmetry Ĥcone (p, ϕ), Π = 0. Its presence guarantees that there exists a basis in which Π = 1 σ τ 3 and the adiabatic Hamiltonian becomes block-off diagonal according to Ĥcone (p, ϕ) = Â(p, ϕ)τ + + Â † (p, ϕ)τ − [41].Here, we introduced the off-diagonal matrices τ ± = (τ 1 ± iτ 2 )/2, along with the upper off-diagonal Hamiltonian block: The topological properties of the adiabatic Hamiltonian are encoded in the topological invariant w 3 ∈ Z [41].This coincides with the winding number of Â(p, ϕ).We now write w 3 = ´2π 0 dϕ w 3 (ϕ)/2π, where we defined the winding number densities w 3 (ϕ) = ´dp w 3 (p, ϕ)/2π and: Notably, in the case of Dirac systems, w 3 (p, ϕ) and w 3 (ϕ) are independent of ϕ.Consequently, besides w 3 , also w 3 (ϕ) is quantized.In particular, for a phase twist of 2π we have w 3 (ϕ) = w 3 .The behaviors of w 3 (p, ϕ) and w 3 (ϕ) become relevant here since, for B = 0, we find [29]: where E(p) = (υ D p) 2 + ∆ 2 and we made use of the fact that D xx,yy = D cone and D xy,yx = 0.The above result highlights that the outcome for the superfluid stiffness is determined by the topological properties of the adiabatic Dirac Hamiltonian.In fact, we are now in a position to prove that πD cone /∆ is a topological invariant itself.For this purpose, we note that for the evaluation of w 3 (p, ϕ) we can linearize Â with respect to ϕ according to Â(p, ϕ) ≃ g(p) • σ + i∆ϕ1 σ .Here, we defined the vector g(p) = − υ D p x , −υ D p y , ∆ , where |g(p)| = E(p).Under these assumptions, we find: where we took into account that ϕ can be set to zero after ∂ ϕ is carried out.In the above, we introduced The antisymmetry of the integrand under the exchange p x ↔ p y further allows to obtain the expression: with the unit vector ĝ(p) = g(p)/|g(p)|.When the momentum space is compact, the above integral yields an integer, that coincides with the 1st Chern number of the negative eigenstate of Â0 (p).However, due to the Dirac nature of the adiabatic Hamiltonian, the integral yields 1/2 [42].Hence, the quantization of the superfluid stiffness directly probes the presence of a Weyl point at the origin of the synthetic coordinate space (p x , p y , ∆), since it is exactly there where Â0 (p) and |g(p)| become zero.
The origin of the nontrivial topology can be traced back to the properties of the Hamiltonian in the absence of a phase bias: This Hamiltonian features two chiral symmetries effected by the operators Π = σ 3 τ 3 and Π′ = τ 2 , which lead to a unitary symmetry with the generating operator Ô = σ 3 τ 1 .Employing the unitary transformation Π+τ 1 / √ 2 block diagonalizes the unitary symmetry operator according to Ôτ = τ 1 σ and the Hamiltonian into the blocks Ĥτ (p) = τ g(p) • σ.Since the occupied bands of the two massive Dirac Hamiltonians yield opposite fractional 1st Chern numbers with values τ /2, we conclude that the superfluid stiffness is proportional to the difference of these two 1st Chern numbers.The latter can be viewed as a fractional "spin" [43] or dipole [29] 1st Chern number.
It is important to examine the robustness of the quantized value of the superfluid stiffness against external perturbations.We first investigate situations which preserve chiral symmetry, as for instance the introduction of disorder in the pairing gap and the application of strain.The impact of weak and uncorrelated spatial disorder in the pairing gap is analyzed in our supplementary file [44]  π for B = 0. Across B/∆ = 1, the total superfluid stiffness drops by approximately 2∆/π as shown in (d).All energies are expressed in units of t, i.e., the nearest neighbor hopping in the absence of strain.Finally, for the numerics we have replaced the delta function entering in Eq. ( 28) by a Lorentzian with width Γ = 0.001.
within the first-order Born approximation [45].We find that the quantization persists, but in the disordered case ∆ is replaced by its spatially averaged value.
We now consider the presence of strain which varies linearly in space.This scenario is particularly interesting, since the low-energy description of strained graphene solely consists of pLLs which possess a perfectly flat dispersion, therefore allowing us to make connections with topological constraints associated with quantum geometry.We adopt a strain profile which conserves the p y momentum and yields the Hamiltonian [20,34,38,44]: where Here, ℓ B denotes the pseudomagnetic length.Each valley supports a single 0pLL which, for the given choice of strain profile, is an eigenstate of σ 3 with eigenvalue −1.The remaining spectrum of ĥB λ (p y ) consists of two families of non-zero-energy pLLs, with eigenenergies ε σ,n (p y ) = σε n (p y ) = σω B √ n for n ≥ 1.Each pLL sees a degeneracy per area which is equal to 1/2πℓ 2 B for a single valley and a single spin projection [20].Since the energy dispersions are flat, it is more convenient to deduce the contribution of each pLL to the stiffness using the standard expression for the interband part D inter [5].In the following we restrict to the λ = +1 valley and denote the corresponding pLL eigenvectors as |u α (p y )⟩, where α = (σ, n) compactly labels the two quantum numbers.As we show in Ref. 44, D inter can be expressed as a sum of contributions defined per level α: where W is the width of the sample.Here, the occupation of each level, where Θ is the Heaviside step function.The operator Mα (p y ) is defined as follows: Equation ( 18) holds as long as ĥB λ (p y )+ε α (q) for all pLLs, which is exactly the case here as we explain in Ref. 44.
By employing Eq. ( 18) at charge neutrality and B = 0, we confirm that now only the 0pLLs contribute to the stiffness, which is still given by Eq. (1).However, the stiffness is now purely geometric.To show this, we use that M0pLL (p y ) = 2 1 − u 0pLL (p y ) u 0pLL (p y ) .This allows us to express the stiffness in terms of the quantum metric of the 0pLL, i.e., g 0pLL (p The result of Eq. ( 1) is recovered by evaluating the integral ´dp y /W by accounting for the degeneracy of the pLLs.
The arising robustness of the superfluid stiffness for weak strains is indeed expected as long as the Dirac points in (p, ϕ) space of the adiabatic Hamiltonian persist.This holds when the pLL degeneracy is much smaller than unity, since then, the Dirac point is not removed but only "moves" in synthetic space.As a matter of fact, our adiabatic approach is still valid for weak strains and we can re-evaluate Eq. ( 14) for p y → p y + λx/ℓ 2 B and x → ϕ/∂ x ϕ.We find that corrections to Eq. ( 1) become negligible as long as ω B ≪ ∆ and ℓ B is sufficiently smaller than the sample's width W for a strip geometry.
So far we have restricted to an ideal (un)strained Dirac cone.However, in realistic experiments the nonconical part in graphene's band structure also contributes and spoils the universal behavior.Nonetheless, it is still possible to isolate the quantized part by introducing an inplane Zeeman field, which can selectively modify the occupation of the various states.In Fig. 2 we present numerical results for the superfluid stiffness of strained armchair graphene nanoribbons (GNRs) at charge neutrality and varying B. Notably, when B crosses the energy bands, the arising Bogoliubov-Fermi points lead to an additional contribution to the intraband part of the superfluid stiffness, with the latter being obtained using the formula [44]: The results of Fig. 2 verify that fields |B| > ∆ enable the experimental detection of the quantized contribution of the 0pLLs.In Fig. 2(d) we confirm the quantized jump of the total stiffness across |B| = ∆, where any deviation from the expected quantization is only due to numerical finite-size effects.Our analysis reveals that, armchair GNRs are better-suited for observing the quantized jump compared to zigzag GNRs.This is because zigzag GNRs harbor additional edge flat bands which appear even without strain [20] and spoil the quantization [44].
Another aspect that remains to be addressed concerns the impact of detuning away from charge neutrality due to a nonzero chemical potential µ which leads to . Now, all pLLs contribute to the stiffness.Specifically, the contribution of the 0pLLs now becomes: Besides a renormalized ∆, an additional factor emerges which diverges for |µ| = ω B /2. Thus, the stiffness can be strongly enhanced by tuning the system to this resonance.In fact, such resonances appear for all pLLs [44].
After investigating the superfluid stiffness, we now evaluate the quantum capacitance [46][47][48][49][50][51] for a Josephson junction fabricated by contacting an ideal strained monolayer graphene hybrid [52][53][54][55][56][57][58][59][60] to a conventional SC, as shown in Fig. 1(b).Using the model in Eq. ( 17) for strained graphene, we find that the two-valley JQC at charge neutrality and B = 0 takes the exact form [44]: with the Hurwitz zeta function ζ.Thus, at charge neutrality, weak strains This result marks the topological regime and is in accordance with the value set by the emergent Lorentz invariance.In contrast, for ω B ≫ ∆ only the 0pLLs contribute with c Q /c 0 → (ω B /2∆) 2 ≫ 1.For ∆ = 0.2 meV and υ D = 10 6 m/s, we find c 0 ≃ 5 nF/cm 2 .Therefore, tailoring the classical capacitance of the junction c Cℓ , so that c Cℓ ≫ c Q , in principle enables the detection of the underlying properties of the Dirac points.This may be possible by observing Coulomb-blockade-induced charging effects, which can be controlled by means of strain [44] and Zeeman field [29] engineering.
We conclude by discussing broader implications of the topological quantization of D and c Q demonstrated here.First of all, the invariance of the two response coefficients for weak strains implies that flat band SCs which are dictated by a solely-geometric quantized stiffness, can be adiabatically mapped to SDSs and the arising quantization linked to the quantum metric can be understood via topological constraints on the total SDS superfluid stiffness.Even more, based on the observation that moiré twisting can be effectively viewed as the application of strain [61], we infer that the nonstandard topological properties proposed here for strained SDSs can be also applicable to their moiré-twisted counterparts.
In this spirit, we expect that the geometric superfluid stiffness of MATBG [10][11][12][13][14][15] can be possibly attributed to the topological superfluid stiffness of untwisted bilayer graphene.Although we leave the verification of the above conjecture for a future work [62], we here stress that such a scenario is indeed plausible because the total superfluid stiffness of a number of s uncoupled graphene layers, or more general of a number of s uncoupled Dirac cones, satisfies the quantization law D (s) = |s|D cone [29].Therefore, the here-proposed connection between Dirac cones and flat band systems opens the door to predicting and linking distinct quantum materials and devices which yet share the same topological superfluid responses.
Aside from the above concept, the nonuniversal topological quantization for D (s) further implies the universal quantization of the admittance modulus Y mod of the SDS for ω = 2∆, i.e., at the absorption edge.This relies on the fact that for 0 <ω ≤ 2∆ and at zero temperature, a fully-gapped nondisordered SC behaves as a perfect inductor with Y mod = D/ω [63].Thus, in our case we find: where the quantum of conductance e 2 /h appeared after restoring e and ℏ.From the above, we confirm that a quantum admittance effect (QAE) emerges for ℏω = 2∆.Currently, however, it is challenging to experimentally observe the above QAE in graphene-SC hybrids, since the chemical potential can be only reduced down to ∼ 1 meV [64].This sets existing platforms in the antipodal regime |µ| ≫ ∆, thus, raising the question of how does the stiffness and the admittance behave in this limit.Re-evaluating the stiffness for µ ̸ = 0 yields that a higher-order Dirac cone with vorticity s carries a stiffness D (s) (|µ| ≫ ∆) ≈ |s||µ|/2π [44].In fact, for such a higher-order SDS, the stiffness is proportional to |s| for all µ [29].Since a gate potential V g modifies µ according to µ + V g , we propose to experimentally measure dD (s) /dV g ≈ |s|/2π which, while it is not a topological invariant, it is still approximately universally quantized.

Supplemental Material
Topological Superfluid Responses of Superconducting Dirac Semimetals Jun-Ang Wang It is important to investigate the robustness of the quantization of the superfluid stiffness for a single Dirac cone in the presence of disorder in the pairing gap.This is crucial, since the pairing gap sets the superfluid stiffness quantum.For this purpose, we consider that the pairing gap becomes modified as ∆ → ∆[1 + f (r)], where the function f (r) is considered to give rise to spatially uncorrelated disorder and can be equivalently expressed using the Fourier decomposition f (r) = ´dq (2π) 2 e iq•r f (q).In the remainder, we focus on weak disorder strengths by restricting to |f (r)| ≪ 1.Under this assumption, it is eligible to employ the Born approximation.We thus infer the self-energy Σ(ϵ, p) of the single particle Green function Ĝ(ϵ, p) due to disorder.This Green function is defined for the Hamiltonian in Eq. ( 5) of the main text through the formula Ĝ−1 (ϵ, p) = iϵ + B − Ĥ(p).Note that here it is also eligible to restrict ourselves to the study of the time-ordered Green function, instead of the retarded one, because we assume that the disorder strength is sufficiently weak to guarantee that possible level broadenings are much smaller than the pairing gap ∆.The validity of such an approach in the present case can be inferred from the analysis of prior studies on Dirac superconductors, e.g., see Ref. 45.According to this work, the Born approximation without selfconsistency is sufficient for examining the effects of disorder in this regime.Hence, under this condition and within the nonselfconsistent Born approximation (BA), we find the expression for the self-energy ΣBA (ϵ, p) in the absence of a Zeeman field (B = 0): Since in the present work we are mainly interested in the effects of disorder near the Dirac point at p = 0, we can set ΣBA (ϵ, p) ≈ ΣBA (ϵ, p = 0).By further assuming slowly varying disorder, we have f (q) ≈ f (0).These assumptions allow us to make progress with analytical manipulations.We find that terms ∝ τ 3 do not survive the angular integration in momentum space and we thus end up with the expression: where we introduced the superconducting coherence length ξ sc = υ D /∆ and an energy cutoff Λ.Since any topological properties are stemming from the vicinity of ϵ ∼ 0, we can approximate √ ∆ 2 + ϵ 2 ≈ ∆.In addition, in the remainder we consider the case Λ ≫ ∆ which further allows us to proceed with the approximation √ Λ 2 + ∆ 2 + ϵ 2 ≈ Λ.By defining the disorder strength γ = |f (0)| 2 ln Λ/∆ / ξ 2 sc (2π) 3 ≥ 0 we end up with the formula ΣBA (ϵ, p) ≈ −γ iϵ + ∆τ 1 and the Green function becomes: where we introduced the renormalization factor Z −1 = 1 + γ, along with the renormalized quantities υD = Zυ D and ∆ = (2Z − 1)∆.Note that the ∆ corresponds to the value of the pairing gap that one identifies using spectroscopic methods.In the situations of interest, ∆ can be inferred by means of spatially averaging the observed gap edge in the density of states of superconducting graphene using scanning tunneling microscopy.We now proceed with examining the impact of this type of spatial disorder on the superfluid stiffness.Firstly, we note that the paramagnetic current vertex, which is given by the term ∂ p Ĝ−1 BA (ϵ, p)τ 3 , 1 remains identical to the one for the clean case.Therefore, the paramagnetic current operator in the disordered case Ĵ (p) BA (p) continues to be given by Ĵ (p) BA (p) = Ĵ (p) clean (p) ≡ − υ(p) even in the presence of the pairing gap disorder.This additionally implies that the superconducting phase simply modifies the pairing gap in this disordered case by means of the shift ∆ → ∆e iϕ(r)τ3 .We are now in a position to define the adiabatic Green function in the absence of a Zeeman field according to Ĝ−1 BA (ϵ, p, ϕ) = Z −1 iϵ − Z ĥ(p)τ 3 − ∆τ 1 e −iϕτ3 .Next in line is to infer the modification on Fpi,ϕ (ϵ, p, ϕ) in the presence of this kind of disorder.Following once again the same steps as in Ref. 29, yields that FBA;pi,ϕ (ϵ, p, ϕ) = Z Fclean;pi,ϕ (ϵ, p, ϕ).This analysis implies that the superfluid stiffness can be now expressed as: with υ(p) = ∂ p ĥ(p), where we introduced ĥ(p) = Z ĥ(p), and ˆ F pj ϕ (ϵ, p, ϕ) is obtained by employing the adiabatic Green function ˆ G(ϵ, p, ϕ) = iϵ − ĥ(p)τ 3 − ∆τ 1 e −iϕτ3 −1 .Notably, the latter does not include the renormalization factor Z and is expressed solely in terms of the renormalized Hamiltonian quantities.The reason why we can obtain these simplified relations is because FBA;pi,ϕ (ϵ, p, ϕ) contains an equal number of Green functions and their inverses, so that the renormalization factors stemming from the Green functions cancel out.The remaining Z factor appearing in FBA;pi,ϕ (ϵ, p, ϕ) is then absorbed by redefining the current vertex.
In conclusion, we find that in the presence of the disorder in the pairing gap, the superfluid stiffness is given by the expressions obtained for the clean case but with renormalized Dirac velocity and pairing gap.Therefore, in the case of superconducting graphene, disorder in the pairing gap leads to a modified stiffness per Dirac cone which now becomes ∆/π.Hence, the quantization persists, but the pairing gap is given by the disorder averaged pairing gap.As we mentioned earlier, ∆ is actually the gap that one obtains experimentally using spectroscopic methods.We thus conclude that the clean and disordered cases do not essentially differ, as long as one defines the quantum of superfluid stiffness in terms of the experimentally observed pairing gap.

B. Strained Graphene Hybrids -Analytics
In the case of strained graphene at neutrality (i.e.µ = 0) we have contributions from two valleys (λ = ±1) and the respective Hamiltonians read as: where ω B = √ 2υ D /ℓ B .We remark that in the above Hamiltonian and in the remainder of this file, we have made for notational convenience the replacements x → x ⊥ and (p x , p y ) → (p ⊥ , q).We see that the Hamiltonian of one valley maps to the Hamiltonian of the other by taking q → −q.Therefore, in the following, we restrict for convenience to the valley with λ = 1 and multiply by a factor of 2 to account for both valleys in the final results.

C. Superfluid Stiffness: Definitions and Calculations
This section relies on the analysis of our accompanying work in Ref. 29, after considering the quasi-one-dimensional case, for a system in a strip geometry with width W .
Finally, we note that the contribution of a zero-energy flat band defined in the entire momentum space reads as: where P 0 = Θ(∆ − |B|) and g 0 (q) = ∂ q u 0 (q) 1 − P0 (q) ∂ q u 0 (q) defines the quantum metric tensor element of the flat band and P0 (q) the respective projector onto it.

Superfluid Stiffness at Charge Neutrality
Having obtained the general expression for D and the matrix elements between states of strained graphene within the Dirac-cone approximation, we proceed with the evaluation of the superfluid stiffness.Given that µ = 0 and by further considering B = 0, we have the following: where we accounted for both valleys by multiplying by a factor of 2 and suitably replaced the integral over q by means of considering the degeneracy of each Landau level.The energy scales appearing above read ε n,σ = σε n = σω B √ n and E n,σ = E n = ω 2 B n + ∆ 2 .Using the above expression, we can evaluate the contribution of each band to the superfluid stiffness.We first evaluate the contribution of the two zeroth pseudo-Landau levels: We continue with n = 1: We now complete the calculation with the evaluation of n ≥ 2.

Superfluid Stiffness away from Charge Neutrality
At this point it is important to evaluate the above contribution when a nonzero chemical potential and a nonzero Zeeman field are introduced.In our companion paper [29] we have examined the case of superconducting unstrained graphene in the Dirac regime.Specifically, we showed that for a higher-order Dirac point with vorticity s, the superfluid stiffness is given as D (s) (B, µ) = |s|D cone (B, µ), i.e., is |s| times larger than the respective stiffness of a single Dirac cone which carries vorticity of a single unit.The superfluid stiffness of a single Dirac cone for B = 0 and µ ̸ = 0 had been previously calculated in Ref. 17.In the same work it was shown that for B = 0 and |µ| ≫ ∆, one finds D cone (B = 0, |µ| ≫ ∆) = |µ|/2π.We now proceed with the case of strained graphene.We find that for a nonzero chemical potential µ all pLLs have now nonzero contributions to the superfluid stiffness.For the 0pLLs we obtain: Besides a renormalized ∆, an additional factor emerges which diverges for |µ| = ω B /2. Thus, the stiffness can be strongly enhanced by tuning the system to this resonance.In fact, such resonances appear for all pLLs, since we find: From the above, we observe that there exist resonances which are obtained by suitably tuning the chemical potential.Therefore, one can enhance the supercurrent drastically by approaching these.These resonances relate to the average energies of two pseudo-Landau levels in the absence of the chemical potential.Hence, this is an interesting feature that could be probed experimentally.Note also that as µ → 0 we have D n=1,σ → 0 for all values of B. Let us complete this analysis by considering the levels with n ≥ 2.
From the above we observe once again resonances for successive pseudo-Landau levels, since the resonance condition reads µ = ε n±1,σ ′ + ε n,σ /2.Even more, taking the limit µ → 0 yields that the contribution of the n ≥ 1 levels is exactly zero in this limit.In Fig. 3 we show numerical results.∆/ωB = 0.1, where only the contributions of the n = 0, 1, and 2 levels are considered.In the limit of zero chemical potential and away from the resonance, the contribution of the n ≥ 1 levels tends to zero.In the two limits, the superfluid stiffness undergoes a divergence corresponding to resonances.

D. Strained Graphene Hybrids -Numerics
We have performed additional calculations for the superfluid stiffness of both armchair and zigzag strained graphene nanoribbons.For all the numerics (including those of the main text) we have replaced the delta function appearing in Eq. ( 28) by a Lorentzian with width Γ = 0.001.
In Fig. 4 we present results for the same parameters as the ones presented in Fig. 2 of the main text, but for a zigzag strained graphene nanoribbon.Note that, in the zigzag case, we adopted a convention such that the pseudo-vector potential has only a nonzero y component, which is generated via considering tensile strain and increases linearly in the x direction.The latter also leads to an anisotropic Dirac velocity.As a result, the pLLs become now weakly dispersive.This is in contrast to the armchair case, where the strain does not affect the isotropic character of the Dirac velocity, thus rendering the pLLs entirely flat.We find that the quantization at B = 0 is already spoiled due to contributions originating from additional edge flat bands that appear for zigzag terminations.The extra edge modes appear even without the presence of strain.Moreover, the drop of the superfluid stiffness across |B| = ∆ is far from being quantized, thus implying that armchair graphene nanoribbons are more suitable for observing the topological quantization phenomena in the superfluid stiffness.
We have also carried out studies for further values of the parameters ∆ and ω B .These are presented in Fig. 5. From Note that all energies are expressed in units of t, i.e., the nearest neighbor hopping in the absence of strain.these results, we first find that increasing ∆/ω B spoils the quantization of the geometric contribution to the superfluid stiffness.In spite of this fact, we find that the drop of the superfluid stiffness across |B| = ∆ remains approximately quantized for armchair nanoribbons.In contrast, the extra edge flat bands appearing for zigzag nanoribbons contribute with a significant amount to the stiffness and introduce substantial deviations from the quantized value.

E. Josephson Quantum Capacitance for Strained Graphene Hybrids
To alternatively observe the underlying topological properties of bulk superconducting graphene in the hybrids under discussion, we consider that the superconducting gap picks up a time-dependent phase ϕ(t).This is possible by contacting the graphene hybrid of interest to a conventional superconductor through a dielectric which suppressess the Josephson coupling of the junction.ϕ corresponds to the phase difference appearing in the junction, while the junction sees a voltage bias V .Note that the gauge invariant electrostatic potential becomes V → V − ∂ t ϕ/2 (where we set e = ℏ = 1).The above coupling naturally leads to the generation of an excess charge density for a nonzero ∂ t ϕ.When ∂ t ϕ is constant, the above phenomenon leads to an additional contribution to the capacitive energy per area of the Josephson junction, which reads as Here, c JJ denotes the total capacitance per area which includes the so-called classical (c Cℓ ) and quantum (c Q ) parts.For the general class of quasi-one-dimensional systems of interest, the quantum capacitance c Q is given by the respective charge susceptibility and takes the form: Tr τ 3 Ĝ(ϵ, q)τ 3 Ĝ(ϵ, q) Tr τ 3 i(ϵ − iB) + ĥ(q)τ 3 + ∆τ 1 (ϵ − iB) 2 + ĥ2 (q) + ∆ 2 τ 3 i(ϵ − iB) + ĥ(q)τ 3 + ∆τ 1 where n labels the eigenstates of ĥ(q).Note that the above was derived under the assumption of an insulating system, i.e., the Bogoliubov-Fermi level set by B does not cross any band.For the case of hybrids harbouring the quantum capacitance Q stemming from the Dirac cones of the two by: 2 2πℓ 2 (n,σ) We focus on the B = µ = 0 case.The sum can be carried out exactly and yields: ζ is zeta function.expression following limiting in the weak strain regime ω B ≪ ∆ we find c µ=B=0 Q = 2 • ∆/πυ 2 D .In the antipodal limit ω B ≫ ∆, we obtain:

FIG. 2 .
FIG. 2. (a)Schematic illustration of the model adopted for an armchair graphene nanoribbon (GNR).We consider nonuniform strain which renders the x component of the pseudo-vector potential nonzero and increasing with increasing y, as indicated by the thickness of the green and purple lines.(b) Electronic band structure of (a) as a function of the conserved wave number kx.First, we numerically obtain the spectrum for an armchair graphene GNR of width W = 601 × √ 3a0, where a0 is the carbon-carbon distance, and ωB ≃ 0.045.Subsequently, we add a conventional superconducting gap ∆ = 0.02.pLLs emerge in the interval−W/2ℓ 2 B ≤ kx ≤ +W/2ℓ 2 B ofthe first Brillouin zone where a0 is set to unity for convenience.(c) and (d) show numerical results for the interband and total superfluid stiffness components as functions of an inplane Zeeman field B for the same values discussed in (b).The dashed vertical red lines indicate the energies for the pLLs in the presence of the nonzero ∆.From Fig. (c) we verify that Dinter ≈ 2∆π for B = 0. Across B/∆ = 1, the total superfluid stiffness drops by approximately 2∆/π as shown in (d).All energies are expressed in units of t, i.e., the nearest neighbor hopping in the absence of strain.Finally, for the numerics we have replaced the delta function entering in Eq. (28) by a Lorentzian with width Γ = 0.001.

FIG. 3 .
FIG.3.Dependence of the superfluid stiffness with the chemical potential for the two limits of ∆, (a) ∆/ωB = 10 and (b) ∆/ωB = 0.1, where only the contributions of the n = 0, 1, and 2 levels are considered.In the limit of zero chemical potential and away from the resonance, the contribution of the n ≥ 1 levels tends to zero.In the two limits, the superfluid stiffness undergoes a divergence corresponding to resonances.

FIG. 4 .
FIG. 4. (a)Schematic illustration of the model adopted for a zigzag graphene nanoribbon (GNR) under the influence of nonuniform strain.The strain is applied in a uniaxial manner so that the y component of the pseudo-vector potential is nonzero and increasing with increasing x, as indicated by the thickness of the green lines.(b) Electronic band structure of (a) as a function of the conserved wave number ky.We first numerically obtain the spectrum for an armchair GNR of width W = 601 × 3a0 and ωB ≃ 0.045.Subsequently we add a conventional superconducting gap ∆ = 0.02.Besides pLLs, additional edge flat bands emerge which are typical for zigzag GNRs and connect the two bulk Dirac cones.(c) and (d) correspond to the numerical results for the interband and total superfluid stiffness as a function of an applied inplane Zeeman field B for the same values discussed in (b).In (c) we also show the contribution of the additional edge modes, which spoils the expected quantization of the interband term for B = 0.The dashed vertical red lines indicate the energies for the pLLs in the presence of the nonzero ∆.Due to these extra edge modes the quantization is more difficult to see in zigzag compared to armchair GNRs.Note that all energies are expressed in units of t, i.e., the nearest neighbor hopping in the absence of strain.

5 FIG. 5 .
FIG. 5. Numerical results for the interband and total superfluid stiffness as a function of an applied inplane Zeeman field B for armchair [(a)-(d)] and zigzag [(e)-(h)]GNRs for the same strain strength ωB ≃ 0.045.The panels (a), (b), (e) and (f) are obtained for ∆ = ωB, while the panels (c), (d), (g) and (h) are obtained for ∆ = 2ωB.Note that all energies are expressed in units of t, i.e., the nearest neighbor hopping in the absence of strain.
1,2 , Mohamed Assili 1 , and Panagiotis Kotetes 1 1 CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China 2 School of Physical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China A. Effects of Disorder on the Quantization of the Superfluid Stiffness