Anton V Khvalyuk
email: anton.khvalyuk@lpmmc.cnrs.fr

Thibault Charpentier

Nicolas Roch

Benjamin Sacépé
email: benjamin.sacepe@neel.cnrs.fr

Mikhail Feigel'
email: mvfeigel@gmail.com

Man Near

Mikhail V Feigel'man

Near-power-law temperature dependence of the superfluid stiffness in strongly disordered superconductors

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Superconducting superinductors, proposed about decade ago [START_REF] Mooij | Phase-slip flux qubits[END_REF][START_REF] Mooij | Superconducting nanowires as quantum phase-slip junctions[END_REF][START_REF] Douçot | Physical implementation of protected qubits[END_REF][START_REF] Brooks | Protected gates for superconducting qubits[END_REF][START_REF] Groszkowski | Coherence properties of the 0-π qubit[END_REF][START_REF] Grunhaupt | Granular aluminium as a superconducting material for high-impedance quantum circuits[END_REF] as important elements of quantum circuits, now constitute an intensively developing sub-field in the physics of strongly disordered superconductors, as some selected examples [START_REF] Astafiev | Coherent quantum phase slip[END_REF][START_REF] Rotzinger | Aluminium-oxide wires for superconducting high kinetic inductance circuits[END_REF][START_REF] De Graaf | Charge quantum interference device[END_REF][START_REF] Nierce | High kinetic inductance nbn nanowire superconductors[END_REF][START_REF] Zhang | Microresonators fabricated from high-kinetic-inductance aluminum films[END_REF][START_REF] Shaikhaidarov | Quantized current steps due to the a.c. coherent quantum phase-slip effect[END_REF] demonstrate. Superconducting films used for the construction of superinductors must combine high kinetic inductance per square L K with low dissipation in the microwave frequency range. Large L K corresponds to small 2D superfluid stiffness Θ = (2ec/ℏ) 2 /L K , and the latter can be achieved close to the Superconducting-Insulator Transition (SIT). The condition of low losses naturally points to a family of superconducting materials in which the SIT occurs without closing the single-particle spectral gap [START_REF] Sacepe | Quantum breakdown of superconductivity in low-dimensional materials[END_REF]: Indium Oxide [1], Titanium Nitride [START_REF] Sacépé | Pseudogap in a thin film of a conventional superconductor[END_REF], Niobium Nitride [START_REF] Mondal | Enhancement of the finite-frequency superfluid response in the pseudogap regime of strongly disordered superconducting films[END_REF], and, possibly, granular Aluminum [START_REF] Levy-Bertrand | Electrodynamics of granular aluminum from superconductor to insulator: Observation of collective superconducting modes[END_REF]. Indeed, the absence of low-energy quasiparticles naturally decreases the absorption of microwave electromagnetic field. However, it does not guarantee the absence of other channels of dissipation.

Far away from the SIT, moderately dirty superconductors described by the semiclassical BCS-like theory have a sharp gap ∆ in the excitation spectrum, leading to exponentially low density of quasiparticles at low temperatures, i.e. n (T ) ∝ exp(-∆/T ). One then expects similar temperature dependence in all other physical quantities, e.g., Θ ∝ ´dE ν (E) E e -E/T ∼ T n (T ). However, such a connection between the DoS and the temperature dependence of Θ is not observed experimentally at strong disorder. In particular, it was found in Ref. [START_REF] Coumou | Electrodynamic response and local tunneling spectroscopy of strongly disordered superconducting tin films[END_REF] that the temperature dependence of the kinetic inductance per square L K ∝ 1/Θ in strongly disordered TiN films is much stronger than the one predicted within usual Mattis-Bardeen model [START_REF] Mattis | Theory of the anomalous skin effect in normal and superconducting metals[END_REF] for the single-electron DoS extracted by Scanning Tunneling Spectroscopy in the same experiment.

In the present paper we report even more striking behavior of the low-temperature suppression of superfluid stiffness, δΘ(T ) = Θ(0) -Θ(T ), in strongly disordered amorphous InO x , as deduced from the dispersion law of one-dimensional plasmon waves in a long superconducting stripe. Namely, we observe the power-law-like temperature dependence δΘ/Θ ∼ (T /T 0 )

b with b ∼ 1.6 -1.7 and T 0 ∼ (1 -7) × 10 1 K. This unusual dependence as well as the large magnitude of the effect confidently defy any semiclassical explanation based on mean thermal quasiparticle density.

Addressing δΘ by a semiclassical approach is additionally hindered by the fact that strongly disordered InO x is known [1,[START_REF] Dubouchet | Collective energy gap of preformed cooper pairs in disordered superconductors[END_REF] to posses a hard gap ∆ P ∼ 5 K in the singleparticle DoS even above the transition temperature T c ∼ 2 K.

arXiv:2311.15126v1 [cond-mat.supr-con] 25 Nov 2023

As a result, the electron-hole quasiparticle excitations can be safely neglected at T ≪ T c and thus cannot account for the experimentally observed suppression of the superfluid stiffness Θ with temperature reported in this paper.

To understand the latter, one instead should turn to the properties of the collective modes. While certain contribution to δΘ comes from the aforementioned 1D long-wavelength plasmonic excitations, a simple calculation presented in Section II below quickly demonstrates that the associated effect is too weak at low temperatures. One should thus address the behavior of short-range collective excitations. This latter issue was initially considered within the approximate analytical theory of Ref [START_REF] Feigel | Microwave properties of superconductors close to the superconductor-insulator transition[END_REF], where it was found that low-energy collective excitations are expected to exist in a broad range of disorder, not necessarily close to the SIT. Unfortunately, the approximation of the space-independent order parameter ∆ employed in Ref. [START_REF] Feigel | Microwave properties of superconductors close to the superconductor-insulator transition[END_REF] was later found to be to rather crude [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF], prompting the issue of short-range low-energy collective modes to be reconsidered.

In the present paper, we show that a proper account of strong inhomogeneity of the order parameter ∆(r) in a specific model of a pseudo-gaped superconductor allows one to describe near-power-law temperature dependence of the superfluid stiffness Θ, with the order of magnitude of the effect comparable to that in the experimental data. We find that the exponent b of the power law decreases with the increase of disorder, with b ∈ (1.5, 3) in a wide range of disorder parameters. We predict such behavior in a broad range of low temperatures, λ∆ 0 ≤ T ≤ ∆ 0 /2, where ∆ 0 is the typical energy scale of the order parameter, and λ ≪ 1 is the dimensionless Cooper coupling constant. The broad distribution of the order parameter P (∆), similar to the one found in Ref. [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF], plays a crucial role in our theoretical description. Our theoretical results are in a good qualitative agreement with the experimental data.

The paper is organized as follows: Section II presents the experimental results; Section III formulates the theoretical model; the results of the calculations (both numerical and analytical) are presented in Section IV. Qualitative comparison between experimental and theoretical results, discussions of our findings, and conclusions are present in Section V. The paper is supplemented by a number of Appendices containing various technical details of the analytical approach.

II. SUPPRESSION OF SUPERFLUID STIFFNESS IN STRONGLY DISORDERED AMORPHOUS INDIUM OXIDE RESONATORS A. Experimental results

The low-temperature evolution of superfluid stiffness in strongly disordered superconductors can be probed experimentally by studying the shift in resonance frequency of a microwave resonator due to the increase of kinetic inductance L K ∝ 1/Θ with temperature [START_REF] Coumou | Electrodynamic response and local tunneling spectroscopy of strongly disordered superconducting tin films[END_REF][START_REF] Grünhaupt | Loss mechanisms and quasiparticle dynamics in superconducting microwave resonators made of thin-film granular aluminum[END_REF], a method originally developed for the field of microwave kinetic inductance detectors (MKIDs) [START_REF] Day | A broadband superconducting detector suitable for use in large arrays[END_REF].

In this work we fabricated open-ended λ/2 microstrip res-onators made from five strongly disordered amorphous indium oxide thin films of different disorder and constant thickness d = 40 nm. The film disorder is characterized in situ at cryogenic temperatures using a co-deposited Hall bar and is shown to be increasingly strong, as evidenced by the significant reduction of critical temperature T c and enhancement of normalstate resistance above the superconducting transition. Details of sample disorder are summarized in Table I. Importantly the disorder range shown here is known to be characterized by the presence of a pseudogap and spatial inhomogeneities of the order parameter [1].

The microwave resonators are long (L = 3.5 mm) and narrow (w = 1 µm) indium oxide strips deposited on a silicon dielectric substrate under which a gold metallic layer acts as a ground plane. Through capacitive coupling to a microwave feedline, a collective motion of Cooper pairs in the resonator can be excited by an AC drive, giving rise to the propagation of plasmon waves with velocity v = 1/ √ lc, (where l and c are kinetic inductance and capacitance per unit length respectively) [START_REF] Kulik | Surface-charge oscillations in superconductors[END_REF][START_REF] Mooij | Propagating plasma mode in thin superconducting filaments[END_REF][START_REF] Camarota | Experimental evidence of one-dimensional plasma modes in superconducting thin wires[END_REF].

The open boundary conditions at each ends of the strip lead to Fabry-Pérot-like standing wave resonances with nearly linear dispersion relation f (k n ) ∼ (1/2π)vk n , where k n = nπ/L is the wave vector for mode n.

The superfluid stiffness of a given InO x resonator can be extracted from the velocity of plasmons since v ∝ 1/ √ l ∝ √ Θ. Upon increase of temperature the superfluid density decreases, leading to a decrease of resonance frequency. The evolution of Θ(T ) with temperature can therefore be extracted from the relative frequency shift δf (T )/f = [f (20 mK) -f (T )] /f (20 mK) ≈ (1/2) δΘ(T )/Θ defined with respect to the lowest achievable temperature. Further details on the experimental setup and samples can be found in [START_REF] Charpentier | Quantum circuits and the superconductorinsulator transition in a strongly disordered superconductor[END_REF].

The main experimental data are shown in Fig. 1a which displays the transmission as a function of temperature for sample TC007-3. The resonance exhibits a frequency shift upon increasing the temperature from 20 mK up to 0.8 K. The resulting relative frequency shift δf (T )/f is plotted in Fig. 1b together with that of the other four samples listed in Table I. We readily see that δf (T )/f is not exponentially suppressed at low temperatures, in stark contrast with usual BCS dirty superconductors in which the superfluid density is reduced as δΘ(T ) ∝ exp {-∆/T }, due to thermally activated quasiparticles. In log-log scale shown in Fig. 1c, these data exhibit a nearly power-law dependence, δf (T )/f ∝ (T /T 0 ) b , at the lowest temperature with an exponent b ∼ 1.7 independent of disorder. This non-BCS power-law dependence of the frequency shift is the central result of this work.

Inspecting Fig. 1c in the high temperature range, we see that deviations from this power law with a stronger frequency shift occur for temperatures above about 0.5 K, as shown by the upward curvature of the data in Fig. 1c. We conjecture that these high-temperature deviations most likely relate to the thermally activated quasiparticles and can be phenomenologically described by standard Mattis-Bardeen (MB) theory to account for thermal activation. We thus describe the entire temperature Table I. Critical temperature, normal-state resistance and zerotemperature superfluid stiffness of disordered indium oxide samples measured in this work. Power law exponent b and temperature scale T0 extracted from the temperature-induced frequency shift are also displayed. Where possible, the energy gap ∆1 corresponding to the best fit to Mattis-Bardeen theory at higher temperatures is presented.

Sample Tc (K) Rn (kΩ/□) Θ(0) (K) T0 (K) b ∆1 (K) TC002 3
δf (T ) f = T T 0 b + δf (T ) f MB , (1)
where (δf (T )/f ) MB is the MB contribution accounting for the high-T deviations. The resulting dashed lines in Fig. 1b fit well data with b = 1.6 -1.7 and ∆ 1 ∼ 2 T c , which is a reasonable approximation of the single particle gap for the moderately disordered samples [1]. For the samples TC007-3 and TC007-2 (in blue and light blue on Fig. 1b andc), there is no MB contribution in the measured temperature range. Interestingly, we found that the T 0 values resulting from the fits scale linearly with the low temperature superfluid stiffness Θ(0) that we extracted from the plasmon dispersion of the resonators. Fig. 2 shows T 0 as a function of Θ(0) together with a linear fit of slope 5.5. This particular dependence points to a phase-fluctuation origin of this power law suppression of the frequency shift and thus of the superfluid stiffness.

These experimental observations suggest the existence of low-energy excitations with energies E ≪ T c which cannot be captured by the semiclassical BCS-like theory. Indeed, for a given density of bosonic excitation states ν (E), the frequency shift is given by ℏδf (T ) ∼ ´dE ν (E) E/ e βE -1 , so for a power-law dependence of δf (T ) ∼ T b with b ∼ 2 one needs a finite density of states at E ∼ T , viz., ν

(E ∼ T ) ∼ E b-2 .
Otherwise, all thermal effects are exponentially suppressed.

B. Thermal excitation of one-dimensional plasma waves

One possible source of low-energy excitations is provided, in principle, by one-dimensional plasmon waves as we discuss below. We first describe the behavior of plasmon excitations at T = 0, and then calculate the effect of the latter on the superfluid stiffness at low T .

It is well known that the application of an electromagnetic drive on a superconductor in the dirty limit leads to a non-linear current-phase relation [START_REF] Maki | The behavior of superconducting thin films in the presence of magnetic fields and currents[END_REF] (see for instance, Eqs. [START_REF] Zhang | Microresonators fabricated from high-kinetic-inductance aluminum films[END_REF][START_REF] Shaikhaidarov | Quantized current steps due to the a.c. coherent quantum phase-slip effect[END_REF][START_REF] Sacepe | Quantum breakdown of superconductivity in low-dimensional materials[END_REF] in Ref. [START_REF] Maki | The behavior of superconducting thin films in the presence of magnetic fields and currents[END_REF], derived in the framework of the Gor'kov's equations). For a one-dimensional superconducting wire in the dirty limit and at low frequencies ω ≪ ∆ such non-linearity translates into the appearance of higher-order terms in the expansion of the 2D current density j 2D with respect to vector potential A:

j 2D = - ρ S d c A 1 -αξ 2 2e ℏc A 2 , (2)
where α = π 4 + 3 4π ≈ 1.02, ρ S is the superfluid density, d is the film thickness, and ξ is the dirty-limit superconducting coherence length. Eq. ( 2) is analogous to the current-phase relation for a Josephson junction I(φ) = I c sin φ ≈ I c φ -φ 3 /6 , where φ is the superconducting phase-difference across the junction and I c is the critical current.

The Hamiltonian describing long-wavelength plasmons along the wire can be split into two parts: H = H 0 + δH. H 0 is the Hamiltonian related to the linear part of the current-phase relation, and can be diagonalized in a basis of normal modes as H 0 = n ℏω n a † n a n where a † n and a n are bosonic creation and annihilation operators. The effect of nonlinearity in Eq. ( 2) translates into the perturbation δH = -α 4 Θξ 2 ´d2 r |∇φ| 4 , where ∇φ is the superconducting phase gradient.

The relevant part of the Hamiltonian then takes the form (see derivation in Appendix F):

H 0 + δH = n ℏω ′ n a † n a n - ℏ 2 n,m K nm a † n a n a † m a m , (3)
with ω ′ n = ω n -(K nn + m K nm ) /4 being the renormalized frequency, and the Kerr coefficients K nm for the strip geometry defined as

K nm = 3α 1 - 1 4 δ nm ξ 2 Lw ℏω n Θ(0) ω m , (4)
where α = π 4 + 3 4π ≈ 1.02 relates to the strength of the current-phase nonlinearity, L and w are the strip length and width, respectively. Eq. ( 3) is well-known in quantum optics and is used to describe the interaction of a given mode n with itself (through the self-Kerr coefficient K nn ), and with another mode m (via the cross-Kerr coefficient K nm ). The corresponding Kerr effect is seen as the reduction of a normal mode's frequency due to occupation of other modes:

ω n → ω ′ n - 1 2 m K nm N m , (5)
where N m = ⟨a † m a m ⟩ is the bosonic occupation number of a given normal mode.

A result similar to Eq. ( 3) for a chain of Josephson junctions can be found in Refs. [START_REF] Weißl | Kerr coefficients of plasma resonances in josephson junction chains[END_REF][START_REF] Yu | Kerr nonlinearity in a superconducting josephson metamaterial[END_REF], where authors find a good agreement between experimental and theoretical Kerr coefficients. Both Eq. ( 4) and the model of [START_REF] Yu | Kerr nonlinearity in a superconducting josephson metamaterial[END_REF] with short-range capacitive coupling give the exact same result if one identifies Θ(0) with the Josephson junction energy E J , and Lw/(2ξ 2 ) with the number of junctions in he chain N , and sets α = 1/6, corresponding to the coefficient of the cubic term in the expansion of sin φ.

We emphasize that the results above are applicable to a homogeneous superconducting strip made of a moderately dirty superconductor in the diffusive limit. The presence of a pseudogap and other non-trivial consequences of strong disorder are, therefore, completely neglected.

We now discuss how one-dimensional plasmons induce a frequency shift as a function of temperature. Upon increase of temperature, the thermal population of plasmonic modes is increased, following the Bose-Einstein distribution N m = e ℏβωm -1

-1 , where β = 1/k B T . As a result, the frequency of a given mode n is reduced due to the interaction with other thermally populated modes. Using Eqs. ( 5) and ( 4), the total frequency shift with temperature can be expressed as

δf (T ) f = 3 2 α ξ 2 Lw m≥1 ℏω m /Θ(0) e ℏβωm -1 . (6)
To calculate the sum over modes, the dispersion law ω (k) of the plasmonic modes with the wave number k is required, and the latter depends on the electrostatic properties of the system. While generally one expects logarithmic corrections to the linear dispersion law due to long-range Coulomb interaction [START_REF] Mooij | Propagating plasma mode in thin superconducting filaments[END_REF], in the present experimental setup the plasmonic modes with kh ≪ 1 are not affected due to the screening from the ground plain at distance h ∼ 300 µm, see Appendix F for details. For the lowest temperatures, one can thus approximate the plasmonic spectrum by a purely linear dispersion relation, corresponding to ω m = mω 1 , where ω 1 is the fundamental mode's angular frequency. For temperatures higher than ℏω 1 /k B , the sum in Eq. ( 6) then evaluates to

δf (T ) f ≈ 3 2 α ξ 2 Lw ℏω 1 Θ(0) π 2 6 ℏω 1 k B T -2 = T T K 2 , (7)
leading to a power-law frequency shift with the temperature scale T K given by

k B T K = 2 π √ α √ Lw ξ Θ(0) ℏω 1 . (8)
Since ω 1 ∝ Θ(0), T K should scale with the superfluid stiffness as T K ∝ [Θ(0)] 3/4 . Eq. ( 7) is applicable while

ℏω 1 ≪ k B T ≪ ℏω 1 L/h ∼ 10 ℏω 1 , which translates to 0.02 K T 0.
2 K for the parameters of the present experimental setup. At higher temperatures, the logarithmic corrections to the plasmonic spectrum should be taken into account, resulting in a weak (logarithmic in T ) enhancement of the effect. Nevertheless, Eq. ( 7) allows one to correctly estimate the magnitude of the frequency shift due to the plasmonic resonances.

In particular, the estimate Eq. ( 7) predicts a power-law frequency shift at low temperatures with an exponent b = 2 close to the experimental value of ∼ 1.6. However, the magnitude of the effect is much smaller than observed experimentally: using a reasonable estimate for the coherence length in disordered Indium Oxide ξ ≈ 5 nm [START_REF] Sacépé | High-field termination of a cooper-pair insulator[END_REF] and the experimental values of Θ(0) and ω 1 /(2π) ∼ 0.4 GHz, one obtains T K ∼ 3900 K for the lowest disorder and T K ∼ 1300 K for the highest disorder, both of which are two orders of magnitude larger than the values of T 0 observed experimentally (see table I and Fig. 2).

We now discuss theoretically how the features of strong disorder can lead to collective modes that suppress the superfluid stiffness at low temperatures.

III. MODEL AND THEORETICAL APPROACH

In the present section, we outline the theoretical approach that consistently describes the superfluid stiffness Θ in a strongly disordered superconductor. Subsections III A and III B present the Hamiltonian of a pseudogapped superconductor and the corresponding current operator; Subsection III C discusses the problem of macroscopic electromagnetic response of a disordered superconductors and presents a semiphenomenological connection between Θ and the statistics of the microscopic current response; Subsection III D then provides a way to calculate the latter in a particular disorder realization by means of a certain generalization of Belief Propagation, and Subsection III E describes the numerical procedure for calculating the statistics of the microscopic response. As a result, one obtains a controllable approach to calculate the temperature dependence Θ (T ) /Θ (0) for various disorder strengths.

A. Model Hamiltonian

As demonstrated both experimentally [1, 20] and theoretically [START_REF] Feigel'man | Fractal superconductivity near localization threshold[END_REF], the materials in question feature localized Cooper pairs even above the transition temperature, whereas the singleparticle spectrum exhibits a spectral gap ∆ P several times larger than the bulk transition temperature. Quasiparticle excitation are thus practically absent at low temperatures, and the low-energy physics is governed by hopping of Cooper pairs between localized single-particle states. This can be captured by the following pseudo-spin Hamiltonian:

H = - i 2ξ i S z i - ⟨ij⟩ 4J ij S x i S x j + S y i S y j . (9)
Here, i is the index of the single-particle state, S α i are the pseudo-spin operators derived from the fermionic operators as

2S z i = a † i↓ a i↓ + a † i↑ a i↑ -1, S + i = a † i↑ a † i↓ , S - i = a i↓ a i↑ .
Then, ξ i are random energies of single-particles states, with probability density at the Fermi level ξ = 0 expressed in terms of the true single-particle density of states ν 0 per spin projection as P (ξ = 0) = ν 0 /n, with n being the electron concentration, and the summation ⟨ij⟩ goes over all pairs of states that interact due to Cooper attraction with amplitude J ij . The latter is given by the corresponding matrix element of the form J ij = g ´dr |ψ i (r)| 2 |ψ j (r)| 2 , with g being an interaction constant, and ψ i (r) corresponding to single-particle wave function of state i. While the value of J ij generally vanishes for sites that are localized sufficiently far apart from each other (further than several localizations lengths ξ loc ), the randomness of the ψ i (r) renders the magnitudes of the matrix elements between spatially close states random. To simplify this situation, one adopts three approximations: i) all pairs of states are classified as either strongly interacting or not interacting at all, ii) each state interacts with a large constant number K + 1 of other states, such that 1 ≪ K + 1 ≪ N loc , with N loc = nV loc being the total number of other states available within the localization volume and iii) the interaction amplitude can be replaced by a constant value J ij = λn/2ν 0 K, where λ ≪ 1 is the dimensionless Cooper coupling constant. While these approximations might seem too crude, a detailed analysis shows [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF] that the ignored effects mostly amount to renormalization of physical quantities and inessential corrections, with one notable exception of approximation iii) discussed later. An extended discussion and derivation of this model can be found in [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF][START_REF] Feigel'man | Fractal superconductivity near localization threshold[END_REF], while Ref. [START_REF] Feigel'man | Superconductorinsulator transition and energy localization[END_REF] directly addresses the Superconductor-Insulator Transition in this model.

The proposed structure of interaction matrix elements J ij suggests the notion of interaction graph, in which each singleparticle state represents a vertex, while each pair of interacting states ⟨ij⟩ corresponds to an edge. In essence, Hamiltonian [START_REF] Rotzinger | Aluminium-oxide wires for superconducting high kinetic inductance circuits[END_REF] reduces the problem to hopping of Cooper pairs (equivalent to hard-core bosons) along this interaction graph, with each site representing a point in real 3D space with some approximate coordinate r i (e.g. the center of mass of the corresponding localized wave function ψ i ), and each edge ⟨ij⟩ describing the tunneling amplitude between the two "points" r i and r j . Crucially, this graph features locally tree-like topology, i.e. for a given site i the neighborhood of size m tree ∼ ln N loc / ln K ≫ 1 hops along the graph is unlikely to contain any loops, which is the direct consequence of the sparsity of the graph controlled by the parameter K/N loc ≪ 1 [START_REF] Bollobás | Random graphs. Number 73 in Camridge studies in advanced mathematics[END_REF]. As a result, the graph is indistinguishable from a portion of the Bethe lattice as long as any local quantity is concerned, allowing one to apply the rich palette of methods available for analysis of tree-like systems (as done e.g. in Ref. [START_REF] Feigel'man | Superconductorinsulator transition and energy localization[END_REF]). On the other hand, at distance r ∼ ξ loc √ m tree ≫ ξ loc from a given point r i the actual 3D structure of the model reveals itself via the presence of many loops of length ≥ 2m tree hops containing the chosen site i in the interaction graph. This connection between the structure in the real 3D space and the interaction graph allows one to calculate spatially resolved quantities, as demonstrated e.g. in Ref. [START_REF] Mirlin | Localization transition in the anderson model on the bethe lattice: Spontaneous symmetry breaking and correlation functions[END_REF], where a similar graph model was used in the framework of single-particle Anderson localization problem. The model ( 9) possesses a natural superconducting energy scale ∆ 0 ≈ e -1/λ n/ν 0 , as explained in [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF]. In what follows, all energy quantities in the problem are expressed in units of ∆ 0 . Along with this scale goes the dimensionless disorder strength κ, also introduced in [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF]:

κ = λn 2ν 0 ∆ 0 K = λe 1/λ 2K . (10)
We expect that all qualitative properties of the model are defined by physical parameters, such as temperature T and the dimensionless strength of disorder κ. In particular, the value of other microscopic parameters, such as λ and K, are not essential as long as the values of κ and T (in units of ∆ 0 ) are set.

B. Current operator

The only possible way to transfer charge in the system described by the Hamiltonian of Eq. ( 9) is hopping of the Cooper pairs between different states. This implies that the current density operator is described as

j (r, t) = 1 2 e I e D e (r) , (11)
with the sum going over all directed edges of the interaction graph, hence the factor 1/2. Here, I e is the operator of the current through a given directed edge e = i → j:

I e = 2eV e - 4e 2 c N e A e , (12)
where V e plays the role of the velocity operator

V i→j = J ij S x i S y j -S x j S y i , (13)
and the second term represents the diamagnetic contribution to the current:

N ij = J ij S x i S x j + S y i S y j , (14)
A e = ˆd3 r (D e (r) , A (r))

with A (r) being the vector potential. Both in j (r) and in the diamagnetic term, D e (r) is a short-range vector field that translates the graph topology to the real space, i.e. describes the distribution of the current density j (r) induced by the process of hopping of a Cooper pair from one site to the other. The exact value of D e (r) is expressed via the response of the interaction matrix elements J ij to external vector potential and it thus also inherits the randomness of the matrix elements themselves (see Subsection III A). Importantly for us, D e is antisymmetric w.r.t the edge direction, viz. D i→j = -D j→i , and it also obeys the following exact identity due to charge conservation in real space:

2divD i→j (r) = |ψ i (r)| 2 -|ψ j (r)| 2 , (16)
where ψ i (r) are single-particle wave functions.

C. Macroscopic superfluid stiffness

We determine the macroscopic superfluid stiffness at low frequencies via the relation Θ = (ℏ/2e) 2 ρ S d, where d is the film thickness, and ρ S is the superfluid density entering the London's equation:

j = - 1 c ρ S A ext , divA ext = 0, (17)
where A ext and j are, respectively, the external vector potential and the current density averaged over a macroscopically large region. The current density is, in turn, determined from the response equation at vanishing frequency ω → 0:

j α = -ˆd3 r ′ Q αβ 0 (r, r ′ ) A β (r ′ ) , (18)
with Q ω being the Fourier transform of the microscopic current response to external vector potential

Q αβ (t -t ′ ; r, r ′ ) = -δ ⟨j α (t, r)⟩ /δA β (t ′ , r ′ ).
Due to the microscopic disorder, the Q kernel depends on both coordinates r, r ′ rather than on their difference rr ′ and has a nontrivial tensor structure. As a result, the current density induced by the external field A ext does not automatically satisfy the charge conservation,

divj = 0, (19)
so the current j induces additional electromagnetic field A int such that the current response [START_REF] Coumou | Electrodynamic response and local tunneling spectroscopy of strongly disordered superconducting tin films[END_REF] to the total field A ext + A int satisfies Eq. ( 19). Consider first the case of weak disorder with diffusive transport in the normal state characterized by k F l ≫ 1, where l is the mean free path, and k F is the Fermi wave number. The average current response j to a smooth external field A ext obeying divA ext = 0 already satisfies the charge conservation [START_REF] Mattis | Theory of the anomalous skin effect in normal and superconducting metals[END_REF], while the disorder-induced deviation δj = jj turns out to be small [START_REF] Spivak | Mesoscopic fluctuations of the superfluid current density in disordered superconductors[END_REF]:

δj 2 /j 2 ∼ √ ξ0/l k F l ≪ 1
, where ξ 0 is the zero-temperature coherence length. This allows one to neglect the contribution of A int to current response and calculate ρ S simply as ρ S δ αβ = c ´d3 r ′ Q αβ 0 (r, r ′ ), where • denotes average over disorder, which in this case is carried out by means of the impurity technique [START_REF] Abrikosov | Methods of quantum field theory in statistical physics[END_REF].

In our model, on the other hand, the statistical and spatial fluctuations jj of the response to a smooth external field are much larger than j itself, necessitating a consistent account of the induced field A int . However, we can still assume that the spatial scale at which the total response j convergence to its average value is much smaller than the London's penetration length λ L = c 2 /4πρ S . In this case, one can neglect the induced magnetic field B = rotA int at the relevant length scales, so A int is dominated by its potential component, viz.

A int ≈ -c iω ∇ϕ. The current distribution is then described by Eq. ( 18), where the induced component A int of total electromagnetic field A ext + A int is found self-consistently from the following system of equations:

divj = 0, rotA int = 0, (20)
In particular, due to the assumed division of scales, one can set A ext = A ext = const. A more detailed derivation of ( 20) is presented in Appendix A. Due to Eqs. ( 11) and ( 16), the problem (18-20) is equivalent to the following discrete system of equations on the values of the onsite electric potential ϕ i :

j∈∂i

I i→j = 0, E i→j = -(ϕ j -ϕ i ) , I i→j = c iω Q ij (ω = 0) E i→j , (21)
with the sum in the first equation going over all neighbors of site i, and

Q ij (ω) = ´dt e iωt Q ij (t)
being the Fourier transform of the current response of the graph edge ⟨ij⟩:

Q ij (t) = - 4e 2 c [⟨[V i→j (t) , V i→j (0)]⟩ -⟨N ij ⟩ δ (t)] . (22)
While one would expect non-local current response ∝ ⟨⟨I e I e ′ ⟩⟩ with e ̸ = e ′ to be also present in Eq. ( 21), all such contributions vanish due to locally tree-like structure of the underlying graph, as explained in detail in Appendix B. Eq. ( 21) should be solved for the values of ϕ i for all sites i of the system, which then also yields the values of the edge currents I i→j .

The distribution of the electric potential in the real space ϕ (r) and the associated electric field E (r) = -∇ϕ (r) are then restored by inverting the following relation:

ϕ i = ˆd3 r ϕ (r) |ψ i (r)| 2 , (23)
where ψ i (r) are the single-electron wave functions, see Subsection III A. Upon also computing the current density j (r) with the help of Eq. ( 11), one is able to calculate the true superfluid density by means of the last relation in Eq. [START_REF] Dubouchet | Collective energy gap of preformed cooper pairs in disordered superconductors[END_REF]. One can, however, avoid the procedure of calculating the fields ϕ (r) and j (r) in real space altogether by noting that Eq. ( 21) is structurally identical to the Kirchhoff's law and the Ohm's law for the interaction graph, with Q ij c/iω playing the role of graph edge conductance and ρ S /iω corresponding to macroscopic conductance. Macroscopic Ohm's law then suggests that for a system with the geometry of a brick with macroscopical sizes L×w×d the total "conductance" along the L direction is equal to ∆ϕ/I total = (ρ S /iω) (L/wd), where I total is the total current through all boundary sites and ∆ϕ is the potential difference between the two boundaries. This constitutes a way to calculate the superfluid stiffness numerically for a given realization of disorder, as done in Appendix C.

The resemblance of Eq. ( 21) to the local Ohm's law also bears certain physical meaning: according to the second Josephson's relation, the quantity φ i = -2e ℏ ϕ i /iω is precisely the superconducting phase of a given site, whereas the last two relations in Eq. ( 21) express the first Josephson's relation linearized w.r.t the phase difference. The local structure of Eq. ( 21) thus becomes a direct consequence of the fact that the system conducts via coherent tunneling of the Cooper pairs. We are interested in estimating the macroscopic superfluid stiffness analytically, and there are currently no means to do this for a general setting, especially given that the values of the current response Q ij are randomly distributed across many orders of magnitude. However, our numerical experiments on the solution of ( 21) for 2D systems have shown (see Appendix C for details) that the following approximate expression is applicable to our problem for κ ≥ 1:

Θ (T ) ≈ C exp {⟨ln Q ij (ω = 0)⟩} = C • Q typ (T ) , ( 24
)
where C is a prefactor that depends only on the details of the graph structure (such as K, ξ loc and concentration n), but not on temperature T . As it is explained in Appendix E, one can estimate C ∼ nL 2 0 dK, where L 0 ∼ ξ loc is the interaction length from Subsection III A, n is site concentration, and d is the thickness of the film. The exact value of C is additionally modified by the presence of short-range correlations in Q ij , but this does not change the order of magnitude of the answer, as explained in Appendix C. We therefore approximate the temperature change of the superfluid stiffness by that of the typical current response Q typ .

The qualitative origin of the result ( 24) could be understood by noting that both the original problem and its discrete reformulation [START_REF] Feigel | Microwave properties of superconductors close to the superconductor-insulator transition[END_REF][START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF] are similar to the problem of macroscopic conductivity of disordered media studied in the seminal paper [START_REF] Dykhne | Conductivity of a two-dimensional two-phase system[END_REF] by A.M. Dykhne. In Ref. [START_REF] Dykhne | Conductivity of a two-dimensional two-phase system[END_REF] it was shown that the macroscopic conductance σ eff of a 2D random medium is equal to the typical value of the microscopic conductance σ typ = exp {⟨ln σ (r)⟩} provided that ln σ (r) is distributed symmetrically around its mean value. At first glance, our problem is different: it is three-dimensional, and there are no a priori reasons to assume that either the actual current response Q αβ ω (r, r ′ ) or its discrete counterpart Q ij satisfy the requirement on the symmetry of the distribution (although for large enough κ the distribution turns out to be sufficiently close to the symmetric log-normal one, as shown in Appendix C). However, there are no qualitative reasons for Θ to differ substantially from the typical current response (at least, in temperature dependence). Indeed, while the calculation of Ref. [START_REF] Dykhne | Conductivity of a two-dimensional two-phase system[END_REF] formally relies on certain duality properties of the 2D problem, it still illustrates the general intuition for the conductance of random media: regions with large conductance can be short-circuited, while regions with small conductance do not conduct at all, rendering the typical conductance to be the only relevant scale of the problem. In the future, we plan to address the validity of the approximate relation [START_REF] Day | A broadband superconducting detector suitable for use in large arrays[END_REF] in more detail.

D. Current response in a given disorder realization

The next step is to calculate the local current response Q ij . This is done by means of the extension of the Method of Belief Propagation (BP) to quantum problems. The standard BP scheme is applicable to classical Hamiltonians with discrete local degrees of freedom (e.g. Ising spins). These include several prominent examples in the theory of spin glasses, where the BP method is also known as the cavity method [START_REF] Mézard | The bethe lattice spin glass revisited[END_REF][START_REF] Mézard | The cavity method at zero temperature[END_REF]. In essence, BP method tries to accurately describe the structure of local two-point correlations in the Gibbs ensemble, resulting in a system of algebraic self-consistency equations that describe the conditional distributions of one site with respect to the other [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF]. Among other things, the BP approach is known to be exact for the systems with tree-like topology, i.e. the ones with no loops in the interaction graph [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF].

Remarkably, this latter property is not sensitive to the particular structure of the local degrees of freedom, and a generalization is available to thermal averages of quantum Hamiltonians, such as Eq. ( 9). The quantum nature of the problem implies replacing the algebraic equations with functional ones, with the latter being analytically intractable in general case. In some cases [START_REF] Biroli | Quantum thouless-andersonpalmer equations for glassy systems[END_REF] it is sufficient to characterize the whole problem, in addition to the set of local static fields, by the common intensity of effective quantum noise acting on the spins. In our problem, corrections to the static approximation can be shown to be proportional to small interaction constant λ ≪ 1 and definitely negligible at κ ≤ 1. Although it is not evident that these quantum corrections are irrelevant at large κ as well, the numerical analysis provided in Sec. IIIG of Ref. [START_REF] Feigel'man | Superconductorinsulator transition and energy localization[END_REF] indicates they are weak even very close to SIT. In the present paper we will use static approximation, leaving the account of quantum noise for future work. One can interpret the proposed Approximate Quantum Belief Propagation (AQBP) scheme as the mean field theory of the same type as that of Ref. [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF] but with the Onsager reaction terms taken into account, essentially representing the analogue of the classical Thouless-Anderson-Palmer equations [START_REF] Thouless | Solution of 'solvable model of a spin glass[END_REF] for our system.

The central object of the AQBP scheme is the set of local fields h α k→i that describe the contribution of spin k to the action of spin i, so the latter is described by the effective Hamiltonian:

H eff = -2ξ i S z i + 1 Z i -1 j∈∂i 2h α j→i S α i , (25)
where Z i is the number of neighbors of the given site i (equal to K + 1 = const within our model). The AQBP provides the self-consistency equations on the values of those fields:

h α k→i = j∈∂i\{k} J ij 2 Ŝα j j(i) , α = x, y, (26)
where the summation now goes over all neighbors of i except k (i.e. K summation terms), and 2 Ŝα

j j(i)
is the thermal expectation of yet another single-spin Hamiltonian:

H j(i) = -2ξ i Ŝz i -2h x j→i Ŝx j -2h y j→i Ŝy j , (27)
2 Ŝα j j(i)

:= h α i→j B i→j tanh βB i→j , (28)
with

B 2 i→j = ξ 2 j + h x i→j 2 + h y i→j 2 .
In order to compute any observable in a given disorder realization, one has to solve this system of algebraic equations for the h fields. In particular, the values of Q ij at frequency ω are computed as quantum averages over two-spin Hamiltonians that depend on the h fields:

Q ij (ω) (2e) 2 /c = 3 m,n=0 |⟨m |V i→j | n⟩| 2 1 Z (λ n -λ m ) e -βλn -e -βλm (λ n -λ m ) 2 -(ω + i0) 2 + 3 n=0 e -βλn Z ⟨n |N ij | n⟩ , Z = 3 n=0 e -βλn . (29)
Here V i→j and N ij are given by ( 13) and ( 14), and {λ n , |n⟩} , n = 0, 1, 2, 3 are the eigensystem of the following two-spin Hamiltonian:

H ⟨ij⟩ = - n=i,j 2ξ n S z n -4J ij S x i S x j + S y i S y j - α=x,y 2h α i→j S α j + 2h α j→i S α i . (30)
At this point we can identify the main reason behind the temperature dependence of Q typ in Eq. [START_REF] Day | A broadband superconducting detector suitable for use in large arrays[END_REF]. Almost by definition, this quantity is contributed by the most probable disorder configurations. Now, the majority of the disorder configurations have |ξ i | , |ξ j | ≫ J, h, T , hence the spectral problem for Hamiltonian [START_REF] Weißl | Kerr coefficients of plasma resonances in josephson junction chains[END_REF] can be addressed perturbatively, rendering λ ≈ ±ξ 1 ± ξ 2 , with corrections being small as J/ξ, h/ξ. As a result, the value of Q ij for such disorder configurations at low temperatures can be estimated as

Q ij (2e) 2 /c ≈ -2 3 n=1 |⟨n |V i→j | 0⟩| 2 E n0 -(ω + i0) 2 /E n0 + ⟨0 |N ij | 0⟩ , (31)
with n = 0 denoting the ground state, n > 0 enumerating the excited states, and

E n0 = λ n -λ 1 ∼ E F ≫ ∆ 0 , J.
Generally, temperature influences the values of h, thus modifying both the matrix elements and the spectral gaps E n0 . The latter, however, are only shifted by a quantity of the order ∆ 0 , J ≪ E F , hence this effect can be discarded along with the frequency dependence, as we are interested in ω ∼ ∆ 0 . The main temperature dependence is thus given by the matrix elements, and by means of the perturbation theory it can be estimated as

Q ij (2e) 2 /c ≈ 4Jh 1 h 2 E 10 E 20 , (32)
from which it immediately follows that

Q typ (T ) Q typ (T = 0) ≈ h 2 typ (T ) h 2 typ (T = 0) . ( 33
)
To obtain this latter Eq. ( 33), we used the fact that h i→j , h j→i are all uncorrelated (see below), so the r.h.s of Eq. ( 33) is expressed via the typical value of the h field. Essentially, this implies that the relative change of the superfluid stiffness mirrors that of the typical order parameter at low temperatures.

E. Statistical properties of the current response

According to the Eq. ( 24), the information about the temperature dependence of the superfluid stiffness is encoded in the statistical distribution of the local current response Q ij in the form of the typical value Q typ . It is thus our aim to compute the average of this quantity over various realizations of disorder. As Eq. ( 29) suggests, the value of Q ij is expressed via the values of ξ i , ξ j on the neighboring sites and on the pair of values h i→j , h j→i . Crucially, all four quantities are statistically independent. Indeed, for region of parameters in question, the solution to Eq. ( 26) is short-correlated due to the large number of summation terms in the right-hand side, as discussed in [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF]. As a result, the statistics of h on a given edge does not depend on whether this edge is a part of a tree or a locally tree-like interaction graph that is realized in our case. Iteratively expanding the r.h.s of Eq. ( 26) then reveals, that these equations possess a directed structure, in contrast to similar equations of Ref. [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF]. In other words, the value of h on edge i → j gathers statistical information only from the the finite branch rooted at i and not containing edge i → j or sites i or j. From this, it follows straightforwardly that all four quantities h i→j , h j→i , ξ i , ξ j are uncorrelated.

As a result, we can substantially simplify the procedure of calculating the value of Q typ : instead of solving the system (26) on a finite size instance of the interaction graph, we can only track the distribution of the h fields with subsequent averaging of the Q ij value (29) over the the distribution of h and ξ. The distribution of h can, in turn, be found using the Method of Population Dynamics (MPD), which boils down to claiming that the two sides of Eq. ( 26) are equal in distribution, which follows from the aforementioned directed structure of Eq. ( 26) and the short range of correlations in the solution. MPD then allows one to efficiently prepare large number of samples from this distribution: given a large initial pool of h values, one updates each value by replacing it with the value of the r.h.s of Eq. ( 26), where the values of ξ are sampled randomly for each term, and the values of h are randomly selected from the current pool. After a number of such updates, the pool of h values converges to a large sample from the target distribution. The required number of iterations is not large because the corresponding distribution does not have fat tails in the given range of parameters [22, Sec. III C]. In our implementation, we terminated the process once the mean value ⟨h⟩ has converged. The pool sizes should simply be sufficient to capture the rare events responsible for the formation of the relevant parts of the distribution. In particular, for all simulations presented below, pools of size 10 7 or larger were used to guarantee the convergence of low-value tails of the distribution.

It is worth noting that the assumption of short correlation distance is crucial for this procedure and is not guaranteed in general. In particular, the work [START_REF] Feigel'man | Superconductorinsulator transition and energy localization[END_REF] analyzes Eq. ( 26) in the context of the disorder-driven Superconductor-Insulator Transition (SIT) and finds that for K ≤ K RSB = λ exp {1/2λ} the distribution of the h fields becomes fat-tailed, and the associated spatial configuration is not at all short-correlated. We are thus interested in the interval K > K RSB , which still contains rich physics, as shown in [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF].

IV. THEORETICAL RESULTS

A. Numerical results

We start by discussing the qualitative shape of the temperature dependence of the local current response presented on Figure 3. There are three temperature regimes: i) at very low temperatures, one observes exponentially small change in the superfluid stiffness, with the numerical method being unable to properly resolve these values. ii) At moderately low temperatures T ∼ 0.05 ∆ 0 -0.3 ∆ 0 the dependence resembles a power law, although the exponent decreases with temperature. Crucially, the apparent exponent of the power law also decreases with disorder, as shown on Figure 4. iii) At higher temperatures 0.3 ∆ 0 ≤ T < ∆ 0 the local current response continues to decrease until it vanishes at the transition point. We do not analyze the region of higher temperatures T ∼ ∆ 0 , which might also be influenced by quaisparticles due to finite value of the single-particle pseudogap ∆ P .

Another important observation from the numerical results is that the relative change in the typical current response δQ typ /Q typ goes in line with that of the mean order parameter δ ⟨∆⟩ / ⟨∆⟩. In particular, it is true that for all points highlighted on Figure 3, regardless of the model parameters. This result is in agreement with the qualitative prediction (33), as our numerical experiments show that ∆ typ and ⟨∆⟩ differ by a temperature-independent factor. In addition, this numerical relation is roughly consistent with the result of Ref. [START_REF] Feigel | Superfluid density of a pseudogapped superconductor near the superconductor-insulator transition[END_REF], where the relation ρ S ∝ ∆ 2 ∼ h 2 is claimed. While our numerical experiments clearly imply that such a relation does not hold literally, it is apparently correct as far as the relative temperature variations are concerned. Eq. ( 34) then allows us to calculate the temperature dependence of Q typ from that of ⟨∆⟩. The latter, as it turns out, is amenable for analytical description.

δ ⟨∆⟩ / ⟨∆⟩ δQ typ /Q typ = 0.5 ± 0.05, (34)

B. Analytical analysis of low-T behavior

While calculating Q ij requires diagonalization of the twospin Hamiltonian [START_REF] Weißl | Kerr coefficients of plasma resonances in josephson junction chains[END_REF], the order parameter ∆ is almost directly accessible from the statistics of the h fields. Indeed, from [START_REF] Kulik | Surface-charge oscillations in superconductors[END_REF] and (26) one obtains:

∆ i = j∈∂i J ij 2 Ŝα j j(i) , ( 35
)
which coincides with the Eq. ( 26) for h itself for up to the difference in the summation number of terms in the r.h.s (K +1 instead of K).

In what follows, we will compute the expectation ⟨h⟩ of the h field, which coincides with the average order parameter ⟨∆⟩ up to an insignificant factor (K + 1) /K, as can be seen from Eq. [START_REF] Bollobás | Random graphs. Number 73 in Camridge studies in advanced mathematics[END_REF]. At the cost of additional technical effort, our theory also allows us to calculate the typical value h typ = exp {⟨ln h⟩}, which is more relevant for the value of the superfluid stiffness according to Eq. ( 33). However, this appears to be unnecessary: the statistical data for the h field obtained by the MPD (see Subsection III E) demonstrates that the relation ⟨h⟩ /h typ is nearly temperature-independent up until the transition region T ∼ ∆ 0 . This also explains why we chose to use the mean value of ⟨∆⟩ rather than the typical value ∆ typ in the numerical relation [START_REF] Feigel'man | Superconductorinsulator transition and energy localization[END_REF]: the former is easier to calculate analytically while being proportional to the latter with a nearly temperature-independent coefficient in the relevant range of temperatures.

The distribution of h fields can be computed straightforwardly by means presented earlier in [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF], although the calculations are considerably simpler due to the absence of the self-action. Indeed, the work [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF] dealt with the necessity to disentangle the mutual correlation between the values of ∆ on neighboring sites due to their equivalence in the selfconsistency equation. On the other hand, the value of h j→i only collects information about the values of ξ k and h l→m in one direction (safe for the present of loops, that, however, do not influence the value of h due to rapid decay of correlations with distance on the graph), This feature has made this problem amenable to the MPD in the first place, and the latter is equivalent to claiming that the two sides of the Eq. ( 26) are equal in distributional sense, with all random variables in the r.h.s being statistically independent:

P h (h) = K k=1 ∞ 0 P h (h k ) dh k ˆPξ (ξ k ) dξ k • δ h - λ K K k=1 f (h k , ξ k ) , (36)
where f (h k , ξ k ) is the spin average in the r.h.s of Eq. [START_REF] Mooij | Propagating plasma mode in thin superconducting filaments[END_REF]. As soon as the equation on distribution of h is obtained in such a way, the approach described below is similar to that of [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF].

In particular, at zero temperature one finds for the cumulant generating function m (s):

m (s) := ln e ish = is ⟨h⟩ + λ ⟨h⟩ κ F (κs) , Ims ≪ e 1/λ , (37)
F (σ) = 1 0 dw e iσw -1 -iσw w 2 √ 1 -w 2 , (38)
with ⟨h⟩ determined self-consistently from

⟨h ln h⟩ = 0. (39)
Here, the l.h.s is calculated by using the expression (37) for m (s), thus representing a function of ⟨h⟩, and Eq. ( 39) has to be solved for ⟨h⟩. In the same way, one can show that the statistics of the order parameter is described by the following generating function:

ln ⟨exp {is∆}⟩ = K + 1 K m (s) , (40)
and we once again remind that both h and ∆ are now measured in units of ∆ 0 . We are interested in K ≫ 1, hence the prefactor in front of m, arising from the different number of summation terms in Eqs. ( 26) and [START_REF] Bollobás | Random graphs. Number 73 in Camridge studies in advanced mathematics[END_REF], can be discarded. It is expected that any effect that is larger than the BCS-like exponential dependence of the form e -2⟨∆⟩/T will be produced by anomalously low values of the order parameter (the unusual factor 2 in the exponent is due to the absence of single-electron quasiparticles in our model). In other words, we are interested in the probability of the order parameter ∆ and the h field to attain values of the order of temperature, with latter being much smaller than mean order parameter ⟨∆⟩. Quantitatively, the extreme value statistics of h is encoded in the asymptotic of the m function at s = it, t ≫ 1, as shown in [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF]. This asymptotic is given by m (it) ≈ -At + Bt ln t, 1 ≪ t ≪ t * ∼ e A/B-1 , (41)

A = ⟨h⟩ {1 -λ [γ -1 + ln 2κ]} , B = λ ⟨h⟩ , (42)
where γ = 0.577... is the Euler-Mascheroni constant, and the estimation for t * follows from the fact that expression [START_REF] Mézard | The cavity method at zero temperature[END_REF] is only valid until it no longer represents a decreasing function. The distribution function P (h) for h < ⟨h⟩ behaves as follows:

P (h) ≈ ζ (h) 2πB 2 exp {-ζ (h)} , (43)
ζ (h) = B exp A -h B -1 . ( 44
)
For instance, at T = 0, one uses Eq. ( 42) to obtain

ζ (h) = λ⟨h⟩e -γ 2κ
exp 1 λ (1 -h/ ⟨h⟩) . Note that while the exponent in ζ (h) features a large factor ∼ 1/λ ≫ 1, the prefactor ∼ λ ⟨h⟩ /2κ in front of this exponent is actually small in the region of interest, rendering the probability density (43) a complicated function of h.

Eq. ( 43) ceases to work at sufficiently small h, where the relevant value of t exceeds t * , rendering Eq. ( 41) inapplicable. This corresponds to ζ (h) reaching the value of Bt * , which happens at

h ∼ h min = B = λ ⟨h⟩ . (45)
Physically, h min plays the role of the minimum value of the h field in the sense that the average e -2h/T for T ≪ 2h min is essentially given by e -2hmin/T , which can also be expressed by the following description of the m function for t ≥ t * :

m (it) ∼ -h min t, t ≫ t * . ( 46
)
The value of h min is then consistent with the continuity of m, i.e. the values of m (it * ) given by the two asymptotic expressions ( 41) and ( 46) coincide.

We also note that the expression ( 41) is only valid for sufficiently weak fluctuations of the interaction matrix elements J ij of the original Hamiltonian (9), viz. δJ/J ≪ λ/κ [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF], whereas our model assumes J ij = const for all connected sites (see approximation iii) in Subsection III A).

The same analysis that lead to Eq. [START_REF] Spivak | Mesoscopic fluctuations of the superfluid current density in disordered superconductors[END_REF] shows that for finite temperatures the value of m (it) acquires a correction

δm (it) ≈ +λt ∂ ∂β 1 0 dw e m(2iβ/w) √ 1 -w 2 , ( 47
)
which is produced by the leading term in the expansion of tanh β ξ 2 + ∆ 2 in Eq. ( 28) in powers of

exp -2β ξ 2 + ∆ 2 .
As already anticipated, the integral in the r.h.s is only sensitive to the asymptotic of the m function at s = it, t ≫ 1. The correction itself amounts to renormalization of the A coefficient in Eq. [START_REF] Mézard | The cavity method at zero temperature[END_REF]:

A (T ) -A (0) = -λ ∂ ∂β 1 0 dw e m(2iβ/w) √ 1 -w 2 . ( 48
)
Eq. ( 39) for the value of ⟨h⟩ also changes to

⟨h ln h⟩ ≈ + ∂ ∂β 1 0 dw e m(2iβ/w) √ 1 -w 2 , ( 49
)
and contains exactly the same w integral as the one in [START_REF] Gao | Experimental evidence for a surface distribution of two-level systems in superconducting lithographed microwave resonators[END_REF].

One has to solve the new equation for ⟨h⟩, similarly to the zero-temperature case.

It is then natural to replace m (it) in Eqs. [START_REF] Gao | Experimental evidence for a surface distribution of two-level systems in superconducting lithographed microwave resonators[END_REF][START_REF] De Graaf | Two-level systems in superconducting quantum devices due to trapped quasiparticles[END_REF] with its intermediate asymptotic expression [START_REF] Mézard | The cavity method at zero temperature[END_REF] to obtain a closed system of equations on A and B = λ ⟨h⟩. The emerging divergence of the integrals at small values of w should be cut at w ∼ w * = 2β/t * , in accordance with the limit of applicability of the intermediate asymptotic [START_REF] Mézard | The cavity method at zero temperature[END_REF]. The remaining region w ≤ w * produces an exponentially small contribution because for t ≥ t * the value of e m(it) is exponentially small: e m(it) ∼ e -hmint ≪ 1. If the temperatures are also exponentially small, i.e. T ≤ T * = 2/t * = e γ e -1/λ /κ ≪ 1 one obtains w * ≥ 1, so a more accurate calculation of the w integral has to be performed in this case. However, the same Eq. ( 46) still implies that the value of the w integrals is of the order of e -2hmin/T ≤ e -hmint * ≪ 1, implying an exponentially small change in all physical quantities, which is consistent with the numerical data on Figure 3.

We also note that Eqs. [START_REF] Sacépé | The fate of the superfluid density near the SIT in amorphous superconductors[END_REF][START_REF] Gao | Experimental evidence for a surface distribution of two-level systems in superconducting lithographed microwave resonators[END_REF][START_REF] De Graaf | Two-level systems in superconducting quantum devices due to trapped quasiparticles[END_REF] are applicable while the approximation 1-tanh x ≈ 2e -2x is applicable for the typical values of x = β ξ 2 + h 2 , which is true for T ≤ T max ∼ 2 ⟨h⟩.

The qualitative behavior of the solution to Eqs. (48-49) can be understood by employing the direct perturbation theory in T . Namely, one treats the change in A and ⟨h⟩ as a perturbation, and the leading order of the latter is given by:

δh ∝ I (T ) = ∂ ∂β 1 ŵ * dw e -A 2β w +B 2β w ln 2β w √ 1 -w 2 , w * = 2 T e 1-A/B , (50)
where the values A, B are taken at T = 0, i.e. from [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF], and the omitted positive coefficient of proportionality is temperature-independent (but does depend on κ, λ). The r.h.s of this equation is negative, and the plot of its absolute value for a reasonable choice of the parameters is shown on Figure 5. It clearly indicates the same type of power-law-like behavior as the one shown on Figure 3, although the power varies with T , as can be see e.g. by the values of the log-derivative d ln I/d ln T .

The comparison of the approximate theory with the numerical calculation by means of the MPD is shown on Figure 6 and is rather satisfying. Discrepancies are only visible at high temperature, where the expressions [START_REF] Sacépé | The fate of the superfluid density near the SIT in amorphous superconductors[END_REF][START_REF] Gao | Experimental evidence for a surface distribution of two-level systems in superconducting lithographed microwave resonators[END_REF][START_REF] De Graaf | Two-level systems in superconducting quantum devices due to trapped quasiparticles[END_REF] are no longer applicable due to higher powers of e -βx in the low-temperature expansion of tanh βx. This theory allows us to give more quantitative description to the claims above:

1. At low temperatures T ≤ T max ∼ ⟨∆⟩ the temperature dependence of the superfluid stiffness δΘ (T ) = Θ (T = 0) -Θ (T ) follows that of the mean order pa- rameter:

δΘ (T ) Θ (T = 0) ≈ 2 δ ⟨∆⟩ (T ) ⟨∆⟩ (T = 0) (51)
2. For λ ⟨∆⟩ ∼ T min ≤ T ≤ T max ∼ ⟨∆⟩ a complicated profile of the temperature dependence of the superfluid stiffness is observed. It can be roughly described by a power law

δΘ Θ ≈ δQ typ Q typ ≈ 2 δ ⟨∆⟩ ⟨∆⟩ ∼ (T /T 0 ) b , (52)
although the power b gradually decreases with T , which also allows other descriptions of the data (see Appendix D). The common trend is that b decreases with disorder strength κ. The qualitative shape of the dependence is described by the integral of the form [START_REF] Smith | Debye mechanism of giant microwave absorption in superconductors[END_REF]. We also note that our analysis suggests that the shape of this dependence is sensitive to certain details of the microscopic model as the latter determine the statistics of low values of the order parameter.

3. At T ≤ T min the dependence of all physical quantities on temperature roughly follows an activation profile:

δΘ Θ ≈ δQ typ Q typ ≈ 2 δ ⟨∆⟩ ⟨∆⟩ ∼ exp - 2T min T . (53)
We also note that the proposed description naturally incorporates the small disorder limit κ ≪ 1. In this case, the distribution of the order parameter is accurately described by a narrow Gaussian distribution, leading to a BCS-like behavior of the form e -2⟨∆⟩/T for both δQ typ and δ ⟨∆⟩, as demonstrated on Figure 7. In terms of our analysis, this corresponds to replacing the m function by a Gaussian one:

m (it) ≈ -⟨h⟩ t - πκλ 4 t 2 , (54)
which corresponds to small-t expansion of the exact expression [START_REF] Spivak | Mesoscopic fluctuations of the superfluid current density in disordered superconductors[END_REF]. This further reinforces the general conclusion: the nontrivial power-law-like behavior of physical quantities is the direct consequence of the nontrivial distribution of the order parameter ∆. We also emphasize that the BCS-like result on Figure 7 is obtained in the model with no quasi-particles (due to large pseudo-gap); as a result, the activation energy is equal to twice superconducting gap, instead of the gap itself. One may wonder if some modified activation-type behavior could describe the data and our simulation results even at large κ. The corresponding analysis is presented in Appendix D. In particular, Figure 11 shows that the activation-law fit δΘ(T ) ∼ exp {-T 1 /T } works in a limited range of temperature, and the latter shrinks as the dimensionless disorder strength κ increases. However, for a larger value of the coordination number, K = 20, we cannot numerically resolve the distinction between power-law and exponential fits, see Figure 12.

V. DISCUSSION AND CONCLUSIONS

We presented both experimental data and analytical theory of the low-temperature suppression of superfluid stiffness, δΘ(T ) = Θ(0) -Θ(T ), in strongly disordered pseudogapped superconductors.

The direct measurement of the superfluid stiffness at low temperatures T ≪ T c in amorphous InO x films revealed (see Fig. 1 c) a strong deviation from the activation-like behavior expected within the semiclassical theory [START_REF] Mattis | Theory of the anomalous skin effect in normal and superconducting metals[END_REF]. We observed instead a power-law

δΘ(T )/Θ(0) ∼ (T /T 0 ) b , (55)
with an exponent b ∼ 1.6 that is weakly sensitive to the disorder strength and a characteristic scale T 0 that is roughly proportional to the superfluid stiffness at low temperatures, viz., T 0 ∼ 5.5 Θ(0), see Fig. 2. A systematic departure from (55) towards higher values of δΘ is observed at higher temperatures (T 0.2T c ≈ 0.5 K) for the less disordered films, and might be attributed to the conventional Mattis-Bardeen contribution due to quasiparticles.

To describe the experimental data, we developed an analytical approach based on a microscopic model of strongly disordered superconductor previously proposed in the literature [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF][START_REF] Feigel'man | Fractal superconductivity near localization threshold[END_REF][START_REF] Ma | Localized superconductors[END_REF]. The key ingredient of the model is the presence of a broad distribution of the local pairing amplitude ∆(r) ̸ = const. The latter is known to exist [1] in a relatively wide range of normal-state resistances and was described theoretically in Ref. [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF]. The qualitative behavior of the whole model is characterized primarily by the dimensionless disorder strength κ. In particular, an estimation for κ can be extracted from the observed shape of the pairing amplitude distribution P (∆) at large values [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF]. The main challenge of the theoretical description is to connect theoretically computable distribution P (∆) with physically observable quantities (e.g., the superfluid stiffness Θ), as the latter can no longer be expressed via the pairing amplitude in a simple fashion. By combining analytical and numerical methods, we have shown that this model also exhibits near-power-law suppression (55) of the superfluid stiffness with temperature (see Figure 3), with disorder-dependent exponent b ∼ 1.6 -3 in the region of large disorder, κ ≫ 1. In particular, the exponent b decreases with disorder and may become less than 2, as shown on Figure 4. The characteristic temperature T 0 in Eq. ( 55) is found to be of the order of the mean pairing amplitude ⟨∆⟩. The dependence (55) takes place in the temperature range λ ⟨∆⟩ ≤ T ≤ ⟨∆⟩ /2, and the temperature range itself reaches low temperatures T ≪ ⟨∆⟩ due to the smallness of the dimensionless Cooper coupling constant λ. However, even within this temperature interval the power-law profile is only approximate, rendering the effective values of b and T 0 temperaturedependent. In particular, Figure 5 demonstrates that the "local" value of b decrease with the growth of temperature. At very low temperatures T λ ⟨∆⟩, the value of δΘ is better described by an activation law δΘ/Θ ∼ exp {-2λ ⟨∆⟩ /T } with a diminished value of the gap in comparison to its semiclassical value. This behavior results in lower values of δΘ than those suggested by Eq. (55). At high temperatures T ∼ ∆, our theory is not applicable to the complete lack of quasiparticles within the proposed model.

However, even in the absence of the latter, we find within the same formalism for the case of small disorder κ ≪ 1 the activation behavior δΘ/Θ ∼ exp {-2 ⟨∆⟩ /T }, but with twice the standard value of the gap, see Figure 7. The effective gap is doubled since the energy 2∆ corresponds in our case to a single excitation instead of two quasi-particles in the usual BCS theory.

As a result, one can claim qualitative agreement between theory and experiment, while a quantitative match is currently beyond reach due to numerous simplifications employed in the theoretical model. The quantitative result of our theory (summarized on Figure 6) is only valid in the model ( 9) with weak statistical fluctuations of the magnitude of interaction matrix elements J ij . As discussed in Ref. [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF], strong fluctuations of the matrix elements render the low-value tail of the distribution of the order parameter even more pronounced and can thus increase the reported suppression of Θ as well as change the overall shape of the temperature dependence. As a result, the latter will still resemble a power-law, probably with a smaller exponent b, due to a substantial fraction of system sites with ∆ i < ⟨∆⟩. An accurate description of the matrix element J ij and the resulting effects on the superfluid density and other physical observables are subjects of future work. We also emphasize that these unusual low-T properties are predicted to exist even at relatively large K ≥ λ exp(1/2λ), where the replica-symmetry-breaking theory of the SIT [START_REF] Feigel'man | Superconductorinsulator transition and energy localization[END_REF] is not yet relevant. An additional element missed in our simplified theory is the implication of the energy dependence of elements J ij , which is present due to Mott hybridization mechanism, as discussed in Ref. [START_REF] Feigel'man | Fractal superconductivity near localization threshold[END_REF]Sec. 2.2.5]; this issue we also leave for future studies.

Among other things to be considered in the future is the issue of the relation between the macroscopic superfluid stiffness Θ and the local current response Q αβ (ω; r, r ′ ) discussed in Subsection III C. The approximate Eq. ( 24) is currently supported by certain preliminary simulations as well as similarity to the problem of the macroscopic response of disordered 2D media discussed in Ref. [START_REF] Dykhne | Conductivity of a two-dimensional two-phase system[END_REF], while a qualitative solution to this problem is to be addressed in the near future. Importantly, the employed approximation does not give access to the superfluid stiffness itself, providing only its relative suppression with temperature δΘ/Θ. This difficulty can be traced down to the fact that kinetic quantities, such as Θ, include information about the embedding of the model in real space. For instance, the semiclassical expression Θ (T = 0) = ℏπ∆/(2e) 2 R □ valid for weak disorder contains this information via the normal-state resistivity per square R □ , with the latter typically determined experimentally. However, sufficiently disordered a:InO films as well as the theoretical model used in the present paper demonstrate insulating behavior in the normal state [START_REF] Charpentier | Quantum circuits and the superconductorinsulator transition in a strongly disordered superconductor[END_REF], so experimentally measured R □ at T > T c might be delivered by a different mechanism than the one responsible for the formation of finite superfluid stiffness Θ for T < T c , rendering R □ irrelevant for the value of Θ. A consistent treatment of this issue would include a more elaborate account of the microscopic properties of the underlying single-particle Anderson localization problem, as the latter unavoidably enter the kinetic quantities.

Another important point worth examining more thoroughly is the validity of description of the local current response proposed in Subsection III D. This comprises neglecting i) the quantum fluctuation and ii) the role of short loops that might be present in the interaction graph of the system. There does exist certain evidence in favor of the first approximation [START_REF] Feigel'man | Superconductorinsulator transition and energy localization[END_REF], but a more detailed analysis is required. The second of these approximations is at least partially controllable [START_REF] Bollobás | Random graphs. Number 73 in Camridge studies in advanced mathematics[END_REF], but because our description involves dealing with extreme value statistics, small loops might be important for certain physical quantities. However, estimations of Ref. [START_REF] Khvalyuk | Distribution of the order parameter in strongly disordered superconductors: An analytic theory[END_REF] indicate that finite concentration of short loops does not change the distribution of the order parameter at low temperature, leaving the effect on physical quantities for future studies.

On purely phenomenological grounds, the presence of nearly power-law behavior (55) indicates the existence of low-energy excitations with energies much below ⟨∆⟩. The nature of these excitations can be revealed by the analysis of dissipation properties of the same microwave resonances used to determine the spectrum of plasmon waves [START_REF] Charpentier | Quantum circuits and the superconductorinsulator transition in a strongly disordered superconductor[END_REF][START_REF] Sacépé | The fate of the superfluid density near the SIT in amorphous superconductors[END_REF]. We expect that the majority of these modes are localized at the scale of the typical interaction length between preformed Cooper pairs (up to few localization lengths). Therefore these modes may represent a specific kind of two level quantum systems (TLS). While those are usually related with atomic degrees of freedom in amorphous solids [START_REF] Gao | Experimental evidence for a surface distribution of two-level systems in superconducting lithographed microwave resonators[END_REF] or localized quasiparticles [START_REF] De Graaf | Two-level systems in superconducting quantum devices due to trapped quasiparticles[END_REF], in our case those TLS are formed by the coherent hopping of the Cooper pairs across several neighboring sites. An indication of the presence of TLS is the non-monotonic temperature dependence of the microwave resonance quality factor Q(T ), with a maximum at some nonzero temperature [START_REF] Charpentier | Quantum circuits and the superconductorinsulator transition in a strongly disordered superconductor[END_REF]. However, the observed data does not rule out the existence of additional mechanisms of dissipation at non-zero frequency. One such mechanism could be an analogue of the Debye relaxation recently discussed in relation with superconductors in the mixed state [START_REF] Smith | Debye mechanism of giant microwave absorption in superconductors[END_REF][START_REF] Pashinsky | Microwave response of type-II superconductors at weak pinning[END_REF]. We are interested in the low-frequency limit of the system (A5). Because the system is superconducting, the zerofrequency response function Q αβ (ω = 0, r, r ′ ) does not vanish. Eqs. A5 then suggest the following substitution to reproduce the correct low-frequency behavior of the response:

j = u iω , A = - c iω E (A6)
which then leads to the following equivalent system:

               divu = (iω) 2 ∆φ 4πε(ω) , -∆ + iω c 2 (-cE) + ∇ (iω) 2 c φ = 4π c u, divE = 0, E = -∇φ r + E, u α (r) = c ´dr ′ Q αβ (ω; r, r ′ ) E β full (r ′ ) . (A7)
One can then take the zero-frequency limit of these equations by simply putting ω = 0, which renders

               divu = 0, ∆E = 4π c 2 u, divE = 0, E = -∇φ r + E, u α (r) = c ´dr ′ Q αβ (0; r, r ′ ) E β full (r ′ ) .
(A8)

As we can see, even in the limit of vanishing frequency, the medium generally responses both by electric and magnetic fields, as rotE = rotE ̸ = 0. However, there's a natural division of scales in this equations, which allows one to simplify the system significantly as far as the value of the superfluid stiffness is concerned. Let's decompose all fields into a sum of fast and slow components, with the distinction made along the scale l m at which one can consider the superfluid stiffness as a self-averaging quantity. The order of l m can roughly be estimated by the size of the tree-like structure of the underlying interaction graph (see Subsection III A of the main text), so l m ∼ Cξ loc , where C is a number of order 2. In this case, all equations except the last in system (A8) independently describe slow and fast components, while the last equation in system (A8) mixes both slow and fast components of the field E due to small spatial scale of the change in the Q kernel. However, from the second equation in system (A8) the fast component of E can be estimated as

E (fast) ∼ 4π c 2 l 2 m u ∼ 4π c l 2 m ρ S E ∼ l 2 m λ 2 L E, (A9)
where λ L = c 2 /4πρ S is the London penetration depth. Because we expect λ L ≫ l m , this contribution to the electric field is small, so we can neglect it and obtain

     divu = 0, E = -∇φ, u α (r) = c ´dr ′ Q αβ (0; r, r ′ ) E β full (r ′ ) , (A10)
which is equivalent to Eq. ( 20) presented in the main text. Note that this only works for the response at small scales, while at large scales one still has to solve the exact system (A8). However, because we expect that the superfluid responses averages over scales l m , one can replace the last equation with a simple London-type relation at large scales:

               divu = 0, ∆E = 4π c 2 u, divE = 0, E = -∇φ + E, u α (r) = ρ S E β full (r) , (A11)
which can be cast in the form of the standard London equations because for ρ S = const the solution for φ is trivial, viz. φ (r) = 0:

         divj = 0, -∆A = 4π c j, divA = 0, j = -1 c ρ S A full . (A12)
Another case where the approximate Eqs. (A10) are applicable is that of a response of a thin film, for which the magnetic part of the response is small regardless of the scale, so one neglect it and put E = 0.

As a result, the system (A7) accurately describes the full electromagnetic response of a disordered superconductor, with the system (A8) being the corresponding low-frequency limit. This includes both the microscopic effects arising from strongly inhomogeneous superconducting state and the macroscopic effects such as the Meissner effect. At the same time, the division of scales or smallness of the absolute value of the currents in the material (such as in case of a thin film) allow one to determine the superfluid density of such a superconductor by using the potential approximation, Eqs. (A10), i.e., neglecting the magnetic part of the response.

Appendix B: Spatial structure of the current response to a potential electric field

In this Appendix, we analyze the spatial structure of the low-frequency current response to a potential electric field. In particular, we derive Eq. ( 21) used to determine the superfluid density in the system via the characteristics of the problem on a graph. We start from verifying the discrete current conservation identity. According to Eq. ( 11), the charge conservation condition divj = 0 corresponds to

0 = 1 2 e divD e (r) I e , (B1)
where I e is the current along directed edge e. Due to identity [START_REF] Mondal | Enhancement of the finite-frequency superfluid response in the pseudogap regime of strongly disordered superconducting films[END_REF], this equation reduces to

0 = i |ψ i (r)| 2 j∈∂i I i→j , (B2)
where we have used current anti-symmetry I i→j = -I j→i . This should be true for any point r, which is only possible if the coefficient in front of each |ψ i (r)| 2 vanishes, rendering the discrete charge conservation, Eq. [START_REF] Feigel | Microwave properties of superconductors close to the superconductor-insulator transition[END_REF].

We should now establish the connection between the current response function in real space and that on the graph. The vector potential is coupled to the system as H EM = -1 c ´dr ′ (j (r) , A (r)), which, together with Eq. ( 11), implies that the response of I e (t) to A (t, r) reads

I e (t) = - 1 2 e ′ ˆdt ′ Q ee ′ (t, t ′ ) A e ′ (t ′ ) , (B3)
where the discrete edge potential A e ′ is given by ( 15), and Q ee ′ (t, t ′ ) is the nonlocal current response on the graph, including the diamagnetic term due to the explicit dependence of the current operator on the vector potential, Eq [START_REF] Zhang | Microresonators fabricated from high-kinetic-inductance aluminum films[END_REF]. According to Eq. ( 20), we are interested in the current response to a potential vector field A = c/iω ∇ϕ, for which Eqs. ( 15) and ( 16) allow to express the corresponding discrete edge potential as

A e = c iω ϕ end(e) -ϕ beg(e) 2 , (B4)
where for a directed edge e = i → j we denoted end (e) = j and beg (e) = i, and the discrete field ϕ i is defined according to Eq. ( 23) of the main text. This expression can further be reformulated in terms of the edge potential induced by a scalar potential applied to a given site j:

A e = j A (j) e , A (j) e =      c iω ϕj 2 , end (e) = j, -c iω ϕj 2 , beg (e) = j, 0,
otherwise.

(B5)

Our next step will be to compute the response of the current along a given edge e to the edge potential of the form (B5) within the AQBP scheme described in Subsection III D of the main text. Together with Eqs. (B3) and [START_REF] Nierce | High kinetic inductance nbn nanowire superconductors[END_REF], this allows us to calculate the density current j (r) in the response to the potential field at low frequencies and eventually derive Eq. ( 21) the main text. The first step is to figure out the way to calculate the nonlocal response, by which we refer to the case when the edge e in Eq. (B3) does not belong the to the nearest neighborhood of site j from Eq. (B5). Within the AQBP scheme, any change in the operator on edge e by a perturbation on a different edge e ′ ̸ = ±e is produced by the response of the order parameter δh e , so one has

Q ee ′ (t, t ′ ) = -ˆdt ′′ δ ⟨I e (t)⟩ δh e (t ′′ ) • δ ⟨h e (t ′′ )⟩ δA e ′ (t ′ ) + δ ⟨I e (t)⟩ δh -e (t ′′ ) • δ ⟨h -e (t ′′ )⟩ δA e ′ (t ′ ) , (B6)
where the two terms in the sum correspond to contribution from the two directions of the target edge. The change of the order parameter is, in turn, described by the time-dependent generalization of the self-consistency equation ( 26) of the main text:

δh k→i (t) = j∈∂i\{k} J ij ˆdt Pi→j (t, t ′ ) δh i→j (t ′ ) + 2e c A i→j (t) ⟨2S x ⟩ i→j e y , (B7)
where e y is the unit vector along y direction, ⟨2S x ⟩ i→j is the static spin average given by Eq. ( 28) of the main text, and

P αβ i→j (t, t ′ ) = 2S α (t) 2S β (t ′
) is the polarization operator of a single spin with the Hamiltonian given by Eq. ( 27), whose Fourier transform reads

Pi→j (ω) = tanh βB i→j B i→j B 2 i→j -ω 2 /4 ξ 2 j iωξj 2 - iωξj 2 B 2 i→j , (B8) with B i→j = ξ 2 j + h 2 i→j ,
and ω shifted by +i0 to restore the retarded structure of the response. In Eqs. (B7-B8) we have set the direction of the unperturbed order parameter h along the x axis, analogous to the real order parameter configuration in conventional superconductors.

The first important thing to note about Eq. (B7) for the edge potential A

e given by Eq. (B5) is that due to the locally treelike structure of the graph and directed nature of this equation, δh e = 0 for all edges that are directed away from the site j where the potential is applied. For instance, one has δh j→i = 0 for any neighboring site i ∈ ∂j, whereas δh i→j ̸ = 0, as we will see later. This is apparent from the fact that for all edges e directed away from site j, the homogeneous part of Eq. (B7) contains only edges e ′ that originate from end (e) and point away from j as well, as explained on Figure 8. The same figure also illustrates that the actual value of δh j→i in the system that contains loops of length at least l can be estimated as δh j→i ∼ (JP ) l/2 δh m→j ∼ (1/Z) l/2 δh i→j because each additional edge traversed from the source of the perturbation to the target edge brings in additional power of J • P ∼ 1 Z , i.e. it is indeed negligible. to the presence of a large loop (red dashed line), the other branch of the graph also has some influence over the value of h, but according to Eq. (B7), each edge in this loop contributes with an additional factor of JP ∼ 1/Z to the magnitude of this influence. As a result, the change of h practically vanishes for all edges pointing away from the source site j. c) The change of h on the edge ending at the source site j is thus only sensitive to the source term itself d) The change of h on the edge at distance greater than 1 away from the source site j is sensitive only to the change of h on the edge that precedes it on the path to the source.

Upon application of the edge potential of the form (B5), δh i→j for any i ∈ ∂j is described only by the source term of Eq. (B7) and reads

δh i→j = - 2e iω s∈∂j\{i} J js ⟨2S x ⟩ j→s ϕ j 2 e y . (B9)
We then immediately recognize the r.h.s of the static selfconsistency equation ( 26) for h i→j , which leads to

δh i→j = - 2e iω h i→j ϕ j 2 e y , (B10)
i.e. the time-dependent change h i→j of the order parameter for edges pointing at j is transverse to the direction of the order parameter itself and directly proportional to the external field φ j . This result is otherwise expected from the gauge invariance: one can move the potential from site j to all sites of the system that are contained in branch rooted at i by a time-depend gauge transform. Because there's no direct response of h i→j to this new field due to the tree-like structure, as described above, the only thing that remains is the gauge change of h i→j , which is given exactly by Eq. (B10).

Consider now an edge e = k → i that is within distance 2 from site j, i.e. end (e) ∈ ∂j but beg (e) ̸ = j (see also Figure 8). This time around, both terms in Eq. (B7) contribute to the result, but the sum contains only one term corresponding to edge i → j that connects the target edge e to the source site j:

δh k→i = J ij Pi→j (ω) δh i→j + 2e c A i→j ⟨2S x ⟩ i→j e y = 2e iω J ij -h i→j Pi→j (ω) e y + ⟨2S x ⟩ i→j e y ϕ j 2 , (B11)
where we have used the Fourier representation of Eq. (B7).

Remarkably, in the ω → 0 limit, we observe from Eq. (B8) that h i→j Pi→j (0) e y ≡ ⟨2S x ⟩ i→j e y , (B12) so at zero frequency we arrive at δh k→i (ω = 0) = 0.

(B13)

Moreover, due to the fact that Eq. (B7) contains only the first term for any edge that is not in the nearest neighborhood of site j, we arrive to

δh e (ω = 0) = 0, (B14)
for any edge e that is not in the nearest neighborhood of site j. This and Eq. (B6) finally implies that in response to a potential on site j the current is induced only along edges that are adjacent to j. Finally, let's compute the current response directly along the adjacent edge, i.e e = i → j for some i ∈ ∂j. This time round, the response of I e to the edge potential (B5) contains both the nonlocal term from the edge potential on adjacent edges as well as the explicit current response to the vector potential on the same edge: We now have to compute the current response to a change of the order parameter along the same edge. The current itself is described by the time-dependent version of the two-spin Hamiltonian [START_REF] Weißl | Kerr coefficients of plasma resonances in josephson junction chains[END_REF], with one of the two h fields now being the function of time. We first note that the current response to the change of the longitudinal component of the δh ∝ e x field vanishes due to the symmetry of the two-spin Hamiltonian. Indeed, consider the transformation S x → -S x , S y → -S y followed by h x (t) → -h x (t), with both transformations applied to both sites. In the absence of the transverse field h y , this transformation leaves both the two-spin Hamiltonian and the current intact, so we conclude that δ ⟨I e (t)⟩ δh x e (t ′′ )

I e = ˆdt ′   -Q ee (t, t ′ ) A (j
≡ - δ ⟨I e (t)⟩ δh x e (t ′′ ) = 0. (B16)
The response to the transverse component δh ∝ e y can be calculated by performing a gauge transform on the corresponding spin, which delivers the following identity:

δ ⟨I e (t)⟩ δh e (t ′′ ) = 1 (2e) /c Q ee (t, t ′ ) h e e T y , (B17)
This and Eq. (B15) the renders the following identity for the value of the response to the potential field at the adjacent edge:

I e = - 1 iω cQ ee (ω) ϕ j . (B18)
According to Eqs. (B3)-(B5), to compute the current response from arbitrary configuration of the scalar potential, one should sum the responses over all sites j, while taking into account the anti-symmetry property I j→i = -I i→j . This finally renders for low frequencies

I i→j = - 1 iω cQ ee (0) (ϕ j -ϕ i ) , (B19)
where we neglected the frequency dependence of the Q kernel. As a result, the current response to a potential field at zero frequency is present only along the edges adjacent to the site were the potential is applied, and the magnitude of this response coincides with the one calculated from the direct effect of the applied potential on this adjacent edge. Together with Eq. ( 11) this yields Eq. ( 21) of the main text.

Appendix C: Verification of Dykhne's law for a discrete system on a locally tree-like graph

The set of numerical experiments presented below demonstrates the approximate validity of the relation [START_REF] Day | A broadband superconducting detector suitable for use in large arrays[END_REF] in the following sense: the macroscopic response ρ S and the typical local response Q typ obey this relation as disorder strength κ is varied within a fixed topology of the graph.

a. Generating the interaction graph.

For a given length L and width w, N ≫ 1 points were spread universally over a rectangle of size L × w, with L ≫ w. For each point, Z neighbors were selected in the neighborhood of radius ξ ≪ L, w (periodic boundary conditions along the w directions where used). To ensure the locally tree-like structure, the condition Z ≪ nξ D was also ensured (here D = 2 is the dimensionality of the problem, and n = N/Lw is the concentration). In case the condition of exactly Z neighbors could not be satisfied for all N sites, the number of such sites was controlled to be negligible (to maintain the average number of neighbors close enough to Z). In this way, the instance of a locally-tree like graph described in Subsection III A of the main text was obtained.

b. Determining the macroscopic response in a given disorder realization.

The next step was to determine the configuration of the order parameter according to Eq. ( 26) and calculate the local edge responses according to Eq. [START_REF] Maki | The behavior of superconducting thin films in the presence of magnetic fields and currents[END_REF]. Note that this procedure keeps track of the configuration of the h fields and Q values in a given disorder realization, in contrast to the MPD described in Subsection III E of the main text that addresses only the distribution of those quantities. Figure 9 shows the distribution of the local responses Q ij in a typical disorder realization.

Eqs. ( 21) (discarding the irrelevant -c/iω factor in front of the current) were then solved numerically with the following boundary conditions: φ i = 0 for a strip of size d × w with n -1/D ≪ d ≪ L near the one of the system's boundaries, and φ i = φ 0 for a strip of the same size situated at the other boundary, thus arranging the geometry of the two-contact measurement. The resulting distribution of the potential was used to determine the distribution of current I i→j . The total current I total was then calculated as a sum of the currents along all edges that connected a site inside the strip of fixed potential with a site outside of this strip. The resulting value of the superfluid density ρ S = φ 0 w/I total L was then recorded for a given disorder realization. is estimated as w0 ∼ 0.8. For these values of the concentration n = N/Lw, Z and ξ, the approximate depth of the tree like structure is m ∼ 6, i.e., the neighborhood of each site of depth ≤ 6 is nearly tree-like, while for larger depths numerous loops start to emerge, indicating that the graph is embedded in the 2D space. c. Verification of macroscopic limit.

To ensure that system size L was sufficient, the procedure described in Step b. was repeated for a set of systems with lengths L ′ < L obtained by truncating the original system, while still controlling the relation L ′ ≫ w. The resulting function ρ S (L ′ ) was then fitted by the dependence A + B/ (L ′ -2d), with the values of the parameters A, B used to ensure good convergence of the estimation ρ S (L) to the macroscopic quantity. In particular, the length scale w 0 = B/A was interpreted as the typical scale of the inhomogeneity of the solution, so the condition w > 2w 0 has been verified as well. Steps a. and b. described above were repeated for 10 disorder realizations. The average value of ρ S and its statistical error represent a single point on Figure 10. Other points where obtained by changing the value of the coupling constant λ (while preserving the structure of the graph) and performing

Step b.

e. The effect of correlations.

The distribution of Q ij possesses short-range correlations that might influence the resulting macroscopic response. To estimate this influence, Step b. was repeated for a system where the values of the local responses Q ij had been randomly reshuffled, so the short-range correlations in the values of Q ij had been destroyed while preserving the distribution of Q ij itself. The resulting macroscopic response ρ shuffled S turned out to be smaller than the correct one ρ S by a disorder-dependent factor of the order of unity: for κ = 7 we obtained ρ S /ρ shuffled S ∼ 4.5, while κ = 0.5 rendered ρ S /ρ shuffled S ∼ 1.1.

f. Conclusion.

The numerical experiment shown above is a satisfactory demonstration of the following statement: the Dykhne's law [START_REF] Day | A broadband superconducting detector suitable for use in large arrays[END_REF] is approximately correct as a function of dimensionless disorder strength κ in 2D systems at zero temperature for a certain locally tree-like graph. Because neither the graph's topology nor the qualitative shape of the Q distribution change at low but finite temperatures, we thus also expect this dependence to correctly describe the connection between ρ S and Q typ as a function of temperature. It is not clear whether the result holds for the 3D problem, which is relevant for the actual experimental setup, as the thickness of the films used in the experiments [? ] is much larger than, e.g., the localization length corresponding to the ξ value at Step a.

Appendix D: Activation-law representation of the numerical and analytical data

The reported dependence of the superfluid density on temperature can be interpreted in multiple ways depending on the presentation. Figure 3 of the main text plots the dependence δQ typ /Q typ in log-log scale, echoing the "power law" point of view on the experimental data of Ref. [? ]. The theoretical description challenges this picture by showing that the apparent power law actually depends on temperature and makes sense only in a certain temperature interval (see Figure 5). In this Appendix, we also present the data of the main text in the form that illustrates the degree of applicability of the activation-law fit δQ typ (T ) ∝ exp {-T 1 /T }.

We start by plotting on Figure 11 the numerical data of Figure 3 in logarithmic scale along the vertical axis and reciprocal scale 1/T along the horizontal axis. Similarly to the data processing of the main text, we then select a subset of points to apply the activation-law fit δQ typ /Q typ = A exp {-T 1 /T } and plot the resulting curve. On Figure 12 we also plot the data of Figure 6 (both experimental and theoretical) with reciprocal scale of temperatures. The plot also features the activation-law fit of the numerical data. The latter appears to describe the data fairly well in a certain range of temperatures, similarly to what is seen on Figure 3 in the main text.

The corresponding value of T 1 for various parameters is shown on Figure 13. As κ decreases, the value of T 1 increases and eventually approaches 2∆ 0 at small κ, as evident from Figure 7 in the main text. At large κ, the value of T 1 can roughly be described by the typical order parameter ∆ typ = exp {⟨ln ∆⟩} at zero temperature, as suggested by Figure 14. Note, however, that such a relation is not directly supported by the analytical approach of Subsection IV B. Indeed, both δQ typ and ∆ typ are complicated functions of the parameters, and the temperature dependence of δQ typ is not, in general, described by an activation law, so the relation T 1 ∼ ∆ typ is only empirical. The order of magnitude of the C constant in Eq. ( 24) can be estimated by examining the limit of weak disorder κ ≪ 1 corresponding to large K. We first note that Q scales as 1/K, as can be seen e.g. from Eq. ( 32), so in order for Θ to remain finite as K → ∞, the C value should contain a factor of K. Moreover, in the limit of weak disorder the value of Q typ in Eq. ( 24) should be replaced by a simple average ⟨Q⟩. Indeed, in the limit of weak disorder one expects Θ to depend only on a local self-averaging quantity, and a natural choice for the latter with the correct dependence on Q and K is the sum of Qs over all neighbors of a given site.

The value of ⟨Q⟩ can then be evaluated directly, as it gains its value on configurations with ξ 1 , ξ 2 ≫ ∆, so Eq. ( 32) is directly applicable with E i = |ξ i |. By cutting the logarithmic integrals over ξ at ∆ 0 from below and E F from above one obtains

K ⟨Q⟩ = ∆ 0 (2e) 2 c 2ν 0 ∆ 0 λn , ( E1
)
where n is the concentration of sites. This results in the following expression for the superfluid density

ρ S = C K K ⟨Q⟩ = C K ∆ 0 (2e) 2 c 2ν 0 ∆ 0 λn . (E2)
On the other hand, Ref. [START_REF] Feigel | Superfluid density of a pseudogapped superconductor near the superconductor-insulator transition[END_REF] provides an estimation of ρ S in a similar model:

ρ S = 2ν 0 e 2 R 2 0 ℏ ⟨∆⟩ 2 , ( E3
)
where ⟨∆⟩ ≈ ∆ 0 is the mean order parameter, R 0 ∼ L 0 ln δloc ∆0 ∼ L 0 /λ, with L 0 being the interaction radius from Subsection III A, and δ loc = ν 0 L 3 0 -1 being the level spacing in the localization volume. Comparing the two expressions, one arrives at

C ∼ nR 2 0 dK. ( E4
)
Appendix F: The cross-Kerr effect in moderately disordered superconductors for a strip geometry

In this Appendix, we address the cross-Kerr coefficients for the plasmonic excitations in a moderately disordered superconductor (obeying semiclassical description). The latter is to be contrasted with the case of strongly disordered superconductors with a pseudogap discussed in Sections III and IV of the main text. We consider a strip of horizontal sizes L × w and thickness d with the following hierarchy of scales L ≫ λ ≫ w ≫ d ≫ ξ 0 ≫ l, with λ being the wavelength of the 1D plasmon in question, ξ 0 being the zero-temperature superconducting coherence length in the dirty limit, and l being the mean free path. In particular, we assume both the charge and current density to be uniform across the thickness of the film. We also neglect the light retardation, i.e. the plasmon frequencies satisfies ω ≪ c/λ. For simplicity, we restrict our derivation to the case T = 0.

The main source of non-harmonicity of the plasmonic modes is the nonlinearity of the current density response j to the vector potential A. In the Coulomb gauge divA = 0, the latter can be written as [START_REF] Maki | The behavior of superconducting thin films in the presence of magnetic fields and currents[END_REF]

j (A) = - πσ∆ c 1 - 4 3π ξ 2 2e ℏc A 2 A, ( F1
)
where σ is the normal-state conductance, and ξ = ℏD/2∆ is the zero-temperature coherence length in terms of the diffusion constant D. The order parameter ∆ also experiences a correction:

∆ = ∆ 0 exp - π 4 ξ 2 2e ℏc A 2 , (F2)
with ∆ 0 being the order parameter at A = 0 and T = 0. The density of the kinetic energy is then obtained by integrating the current density (F1) w.r.t the vector potential:

f kin (A) = - 1 c A 0 (j, dA) = ρ S A 2 2c 2 1 - αξ 2 2 2e ℏc A 2 + O A 4 , (F3)
where in the last line we have used the standard relation ρ S = πσ∆ 0 for the dirty superconductor and expanded the result up to the leading powers of A 2 , yielding α = 4 3π + π 4 , as a result of combining Eqs. (F1) and (F2). For the case of a thin film, the corresponding free energy reads

F kin = volume d 3 r f kin = d film d 2 r f kin , (F4)
where in the second line we assumed that A is uniform across the film's thickness, and the integrations in the two expressions go over the film's volume and surface, respectively. The potential part of the plasmonic Hamiltonian is described by the electrostatic energy for a given charge density:

F pot = 1 2 film d 2 rd 2 r ′ ρ (r) V (r -r ′ ) ρ (r ′ ) , (F5)
where the integration goes over the surface of the film, ρ is the 2D charge density, and V (r) is the electrostatic potential from a point charge (including the effects of screening from the ground plate and the substrate if present). The corresponding electrostatic potential reads [START_REF] Yu | Kerr nonlinearity in a superconducting josephson metamaterial[END_REF]:

V (x, y) = 4π 1 + ε 1 x 2 + y 2 - 4π 1 + ε 2ε 1 + ε ∞ j=1 1-ε 1+ε j-1 x 2 + y 2 + (2hj) 2 (F6)
where h ≫ d is the distance from the film to the ground plane, and ε is the dielectric permittivity of the substrate. In what follows, we will neglect the screening, for simplicity, which amounts to assuming h ∼ L ≫ w.

The Hamiltonian of the plasmonic modes in the absence of retardation then reads

H = 1 2 ˆd2 rd 2 r ′ ρ (r) V (r -r ′ ) ρ (r ′ ) + d ˆd2 r f kin ℏc 2e ∇φ , ( F7
)
where d is the film thickness, and the superconducting phase φ is connected to the vector potential as A = ℏc 2e ∇φ (and also assumed to be uniform across the thickness of the film). The canonical commutation relations for the fields φ, ρ take place:

[φ (r) , ρ (r ′ )] = (2e) i δ (r -r ′ ) .

(F8)

At the boundary of the film, the field φ should satisfy the condition (n, j) = 0, corresponding to the absence of the current through the edges of the film:

(n, ∇φ) r∈boundary = 0, (F9)

where n is the normal vector to the boundary of the film. The Heisenberg's equations of motion associated with Eq. (F7) correspond to the charge conservation law and the Josephson's relation:

ρ = -divj, φ (r) = 2e ℏ ˆd2 r ′ V (r -r ′ ) ρ (r ′ ) . (F10)
By using the expansion (F3) of f kin in powers of A, we can rewrite the Hamiltonian (F7) as

H = 1 2 ˆd2 rd 2 r ′ ρ (r) V (r -r ′ ) ρ (r ′ ) + 1 2 Θ ˆd2 r (∇φ) 2 1 - α 2 (ξ∇φ) 2 + O (ξ∇φ) 4 , (F11)
where Θ = ρ S d (ℏ/2e) 2 is the 2D superfluid stiffness. Strictly speaking, the energy of phase fluctuations in Eqs. (F7) and (F11) also contains terms with higher derivatives of φ, such as, e.g., ξ∇ 2 φ 2 , but those are responsible for the relative corrections of the order O (kξ) 2 to the plasmonic dispersion law, whereas we are interested in the leading nonlinearity w.r.t φ that produces interaction of plasmons. As as a result, such terms can be neglected as far as the small frequencies and wave numbers are concerned.

Low-lying normal modes

To calculate the cross-Kerr coefficients, one first has to extract the normal modes of the Hamiltonian (F11) in the absence of nonlinearity. The latter can be done by finding the oscillating solutions to the classical equations of motions (F10): ρ = e iωt ρ ω (x, y) , φ = e iωt φ ω (x, y) , (F12) rendering

iωρ ω = ρ S d ℏ 2e ∆φ, iωφ ω = 2e ℏ V ρ ω , (F13)
where we introduced the shorthand V ρ ω for the r.h.s. of the second relation in Eq. (F10) for brevity. For a given set of solutions {ω, ρ ω , φ ω } with positive angular frequencies ω > 0, the harmonic part of the plasmonic Hamiltonian (F11) is rewritten as

H = ω ℏω α † ω α ω + 1 2 , ( F14
)
where α † ω , α ω are the creation and annihilation operators of the plasmonic modes satisfying the standard bosonic commutation relations, and the field operators are then expressed as ρ (x, y) = where denoted φ ω ρ ω = ´d2 r φ ω (r) ρ ω (r) and ρ ω V ρ ω = ´d2 r d 2 r ′ ρ ω (r) V (r -r ′ ) ρ ω (r) for brevity.

For the low-lying modes in the strip geometry, one can perform the plain wave ansatz:

ρ ω = 2 L cos kx ρ k (y) , φ ω = 2 L cos kx φ k (y) , (F17)
with ∂ y φ k (±w/2) = 0 as a result of the boundary conditions (F9). This leads to the following system of equations for the profile of the plasmon in the perpendicular direction:

iωρ k (y) = ρ S d ℏ 2e -k 2 + ∂ 2 y φ k , (F18) iωφ k (y) = 2e ℏ w/2 -w/2 dy ′ V k (y -y ′ ) ρ k (y ′ ) , (F19)
V k (y) = ∞ -∞ dx e ikx V (x, y) .

(F20)

The system (F18-F20) can be solved analytically for the case of the potential (F6) in the limit kw ≪ 1 and w ≪ h. We start by discussing the static limit k → 0, as it formally corresponds to the electrostatic problem for the same geometry. The charge distribution is then found from the electroneutrality condition E = 0:

∂ y w/2 -w/2
dy ′ V 0 (y -y ′ ) ρ 0 (y ′ ) = 0, (F21)

where we have used only the y-component of the electric field, as the x component vanishes automatically in the static limit: ∂ x V ρ ∼ k V ρ → 0. In the limit k → 0, the Fourier transform (F20) of the electrostatic potential (F6) reads:

V 0 (y) = 4π 1 + ε   - 2ε 1 + ε ∞ j=1 1 -ε 1 + ε j-1 2 ln w/2 2hj   + 4π 1 + ε 2 ln w/2 |y| + O y h 2 , ( F22
)
where in the second line we have taken into account that w ≪ h.

In the absence of screening, i.e., h → ∞, Eq. (F21) can be solved exactly:

ρ 0 (y) = C (w/2) 2 -y 2 , ( F23
)
where C is a normalization constant. The corrections from finite h can be retrieved by formal expansion in even Chebyshev polynomials: At finite k, the charge distribution ρ k (y) perturbatively deviates from its static profile ρ 0 (y), so we can find the leading approximation for the phase profile φ k by using the Josephson's relation F19 with the unperturbed value of the charge distribution:

ρ 0 (y) = 1 (w/2) 2 -y 2 ∞ n=0 ρ (n) T 2n (2y/w) , (
φ k = 1 iω k 2e ℏ w/2 -w/2 dy ′ V k (y -y ′ ) ρ k (y ′ ) ≈ 1 iω k 2e ℏ w/2
-w/2 dy ′ V k (y -y ′ ) ρ 0 (y ′ ) , (F26)

By construction, the last expression is constant at k → 0, so the leading approximation the phase is distributed uniformly. The exact value of φ k can then be found by integrating the charge conservation (F18):

φ k ≈ - iω k k 2 w w/2 -w/2 dy ρ 0 (y) = - iω k πC ρ S d ℏ 2e k 2 w . ( F27
)
It is instructive to follow how the charge conservation (F18) is respected given the singular profile (F23) of the charge density and nearly constant profile of the phase. The singularity in the charge profile is compensated by a small but singular corrections to φ k at finite k:

φ k = C k + ω k 1 - y w/2 2 3/2 f k (y) , (F28)
where the constant term C k is given by Eq. (F27), and f (y) is a regular function of y that weakly depends on k at k → 0.

Plugging this expression in the r.h.s Eq. (F18) renders a singular term identical to the one present in the l.h.s. Demanding the two Eqs. (F26) and (F27) to be consistent then delivers the plasmonic spectrum:

ω 2 k = ρ S d k 2
´w/2 -w/2 dy ′ V k (y -y ′ ) ρ 0 (y ′ ) ´w/2 -w/2 dy ρ 0 (y)

, ( F29
)
where we have retained the k-dependence in V k as it contains the weak logarithmic correction to the leading approximation.

Using the exact profile (F6) of the electrostatic potential, one obtains asymptotic behavior of the result in the two limiting cases:

ω 2 k = 8π 1+ε ρ S dk 2 w ln 8h we δ(ε) , k ≪ h -1 , 8π 1+ε ρ S dk 2 w ln 8 kwe γ , w -1 ≫ k ≫ h -1 , (F30)
where γ = 0.577... is the Euler-Mascheroni constant, and δ (ε) is defined as

δ (ε) = 2ε 1 + ε ∞ j=1 1 -ε 1 + ε j-1 ln 1 j . (F31)
Note, in particular, that the logarithmic correction to the 1D plasmon dispersion law at w -1 ≫ k ≫ h -1 has been predicted theoretically [START_REF] Mooij | Propagating plasma mode in thin superconducting filaments[END_REF] and observed experimentally [START_REF] Camarota | Experimental evidence of one-dimensional plasma modes in superconducting thin wires[END_REF][START_REF] Charpentier | Quantum circuits and the superconductorinsulator transition in a strongly disordered superconductor[END_REF], with the latter being in perfect agreement, including the numeric coefficient inside the logarithm. Corrections to ω k , ρ k , φ k for h ∼ w and/or kw ∼ 1 can be studied numerically, as has been done to verify the presented results.

Cross-Kerr coefficients

We can now expand the leading anharmonic part of the Hamiltonian in harmonic modes, according to Eq. (F27): (F33)

The calculation of the cross-Kerr effect is then similar to that of Ref. [START_REF] Weißl | Kerr coefficients of plasma resonances in josephson junction chains[END_REF]Appendix B]. We substitute the particular form, Eqs. (F17), (F23), and (F27), of the φ, ρ vectors, rendering a simple expression:

η {ωi} = α 4 Θ 2ξ 2 Lw 4 i=1 ℏω i Θ 1 L ˆdx i sin k i x. (F34)

Figure 1 .

 1 Figure 1. Temperature frequency shift of amorphous indium oxide microwave resonators. a. Transmission resonance as a function of frequency of sample TC007-3 for different temperatures. Upon increasing temperature, the resonance shifts towards lower frequencies. b. Relative frequency shift [f (20 mK) -f (T )] /f (20 mK) for the five samples of Table I. The low-temperature part of the data departs from the standard Mattis-Bardeen behavior but can be fitted by adding a power law contribution as in Eq. (1), as shown by the black dashed lines. The power law exponent is 1.57 ≤ b ≤ 1.69. c. Frequency shift in log-log scale. At the lowest temperatures all curves show a power law with exponent b ∼ 1.6. Above 0.5 K ∼ 0.2 Tc the Mattis-Bardeen mechanism kicks in and the frequency shift retrieves an exponential behavior for the moderately disordered samples (TC002, TC003 and TC001).

Figure 2 .

 2 Figure 2. Dependence of the energy scale T0 obtained from the fit of the frequency shift with Eq. (1) as a function of the zero-temperature superfluid stiffness Θ(0). Dashed line is a linear fit indicating that T0 = 5.5 Θ(0).

Figure 3 .

 3 Figure 3. The temperature dependence of the relative change of typical local current response 1-Qtyp (T ) /Qtyp (T = 0) in logarithmic scale along both axes for various dimensionless disorder strength κ and K = 15. The temperature is measured in units of ∆0. The points with solid filling are selected according to the empirically chosen criteria T ∈ [λ ⟨∆⟩ , λ ⟨∆⟩ + 0.2∆0] with ⟨∆⟩ being the mean order parameter at zero temperature. These highlighted points are used for power-law fitting on Figure 4.

Figure 4 .

 4 Figure 4. Dependence of the parameters b (main panel) and T0 (inset) of the power-law fit (T /T0) b of the portion of numerical data for δQtyp (T ) /Qtyp (0). For each set of parameters, the points for fitting are selected according to the empirically chosen criteria T ∈ [λ ⟨∆⟩ , λ ⟨∆⟩ + 0.2∆0], with ⟨∆⟩ being the mean order parameter at zero temperature, so the dataset for K = 15 corresponds to highlighted points on Figure 3. The value of κ was varied by changing the value of λ while keeping K constant, and various colors correspond to various values of K.

Figure 5 .

 5 Figure 5. Temperature dependence of the r.h.s of (50) for ⟨h⟩ = 0.46, κ = 10, λ = 0.123. The dashed straight lines visualize the log-derivative b = d ln I (T ) /d ln T , green: T = Tmin = 2λ ⟨h⟩, b = 2.48, black: T = 0.4, b = 1.46.

Figure 6 .

 6 Figure 6. Low-temperature dependence of the relative changes in the typical local current response δQ = 1 -Qtyp (T ) /Qtyp (T = 0) and the mean order parameter δ∆ = 1 -⟨∆⟩ (T ) / ⟨∆⟩ (T = 0) according to the numerical MPD (solid dots) and theoretical description given by Eqs. (48-49) (solid line) in logarithmic scale along both axes. The values of δQ are divided by two, which makes relation (34) apparent. The dashed line corresponds to the power-law fit (T /T0) b of numerical data for δQ in the range T ∈ [λ ⟨∆⟩ , λ ⟨∆⟩ + 0.2], with ⟨∆⟩ = 0.49 ∆0 (the values of the fitting curve are also divided by two). The parameters of the model are K = 20, κ = 10. The green region delineates the limits of applicability of the theoretical approach corresponding to λ ⟨∆⟩ ≤ T ≤ ⟨∆⟩ /2.

Figure 7 .

 7 Figure 7. The temperature dependence of the relative change of typical local current response 1 -Qtyp (T ) /Qtyp (T = 0) for κ = 0.25, K = 10 (small disorder). The temperature is measured in units of ∆0. The scale along the horizontal axis is reciprocal (1/T ), and the scale of the vertical axis is logarithmic. The solid line corresponds to fit of the points with T ≤ 0.45 ∆0 with A exp {-2Eg/T }, resulting in A = 2.41, Eg = 0.947 ∆0, which is consistent with the standard BCS theory. Attempting to fit the same points with a power law (T /T0) b results in b ∼ 5.9 (shown with the dashed line).

ACKNOWLEDGMENTS

  Authors would like to thank Denis Basko for numerous fruitful discussions. A.V.K. is grateful for the support by Laboratoire d'excellence LANEF in Grenoble (ANR-10-LABX-51-01). B.S. has received funding from the European Union's Horizon 2020 research and innovation program under the ERC grant SUPERGRAPH No. 866365. N.R. and B.S. acknowledge funding from the ANR agency under the 'France 2030 plan', with reference ANR-22-PETQ-0003. Th.C. and B.S. acknowledge funding from the ANR project ANR-19-CE30-0014 -CP-Insulators. [1] B. Sacépé, T. Dubouchet, C. Chapelier, M. Sanquer, M. Ovadia, D. Shahar, M.V. Feigel'man, and L.B. Ioffe. Localization of preformed cooper pairs in disordered superconductors. Nature Physics, 7(3):239-244, 2011.

Figure 8 .

 8 Figure 8. Visualization of the structure of dependencies in the self-consistency equation (B7) on a locally tree-like graph with coordination number Z = 3. a) the value of the h associated with a given directed edge is influenced by the branch originating from the end vertex of this directed edge (colored in green), and only the values of h along the edges oriented in the same direction contribute to the equation. b) Due to the presence of a large loop (red dashed line), the other branch of the graph also has some influence over the value of h, but according to Eq. (B7), each edge in this loop contributes with an additional factor of JP ∼ 1/Z to the magnitude of this influence. As a result, the change of h practically vanishes for all edges pointing away from the source site j. c) The change of h on the edge ending at the source site j is thus only sensitive to the source term itself d) The change of h on the edge at distance greater than 1 away from the source site j is sensitive only to the change of h on the edge that precedes it on the path to the source.

Figure 9 .

 9 Figure 9. The distribution of the decimal logarithm of the current response functions, as found from the Eq. (29) of the main text. The parameters of the model are T = 0, κ = 7, λ = 0.165, Z = 6. The orange curve corresponds to the Gaussian distribution with the parameters inferred from actual distribution: ⟨log 10 Q⟩ = -4.66 and (log 10 Q) 2 -⟨log 10 Q⟩ 2 = 1.31.

Figure 10 .

 10 Figure 10. Set of points corresponding to pairs (Qtyp, ρS) for various values of the coupling constant λ, while maintaining Z = 6 and T = 0. The corresponding values of κ are (from left to right): 17, 15, 13, 11, 9, 7, 5, 3, 2, 1, and 0.5. The parameters of the numerical experiment of Steps a.-b. are L = 20, w = 2, ξ = 0.05, and N = 960000. The inhomogeneity scale mentioned in Step c.is estimated as w0 ∼ 0.8. For these values of the concentration n = N/Lw, Z and ξ, the approximate depth of the tree like structure is m ∼ 6, i.e., the neighborhood of each site of depth ≤ 6 is nearly tree-like, while for larger depths numerous loops start to emerge, indicating that the graph is embedded in the 2D space.

  d. Averaging over disorder realizations.

Figure 11 .

 11 Figure 11. The temperature dependence of the relative change of typical local current response 1 -Qtyp (T ) /Qtyp (T = 0) in logarithmic scale along the vertical axis and reciprocal (1/T ) scale along the horizontal axis. Various colors correspond to various dimensionless disorder strength κ and K = 15. The temperature is measured in units of ∆0. The points with solid filling are selected according to the criteria T ∈ [λ ⟨∆⟩ , λ ⟨∆⟩ + 0.2∆0] with ⟨∆⟩ being the mean order parameter at zero temperature. These highlighted points are used for exponential fitting A exp {-T1/T }, with the latter plotted with dashed lines.

Figure 12 .

 12 Figure 12. Low-temperature dependence of the relative changes in the typical local current response 1 -Qtyp (T ) /Qtyp (T = 0) and the mean order parameter 1 -⟨∆⟩ (T ) / ⟨∆⟩ (T = 0) according to the numerical MPD (solid dots) and theoretical description given by Eqs. (48-49) (solid line). The values of Qtyp are divided by two, according to relation (34). The dashed lines correspond to the power-law (T /T0) b (blue) and activation-law A exp {-T1/T } (green) fits and of numerical data for Qtyp in the range T ∈ [λ ⟨∆⟩ , λ ⟨∆⟩ + 0.2∆0], with ⟨∆⟩ = 0.49 ∆0. The parameters of the model are K = 20, κ = 10.

Figure 13 .

 13 Figure 13. Dependence of the T1 parameter of the activationlaw fit A exp {-T1/T } of the portion of numerical data for δQtyp (T ) /Qtyp (0). For each set of parameters, the points for fitting are selected according to the criteria T ∈ [λ ⟨∆⟩ , λ ⟨∆⟩ + 0.2∆0],with ⟨∆⟩ being the mean order parameter at zero temperature, so the dataset for K = 15 corresponds to highlighted points on Figure11. The value of κ was varied by changing the value of λ while keeping K constant, and various colors correspond to various values of K.

Figure 14 .

 14 Figure 14. Comparison of the T1 scale of the activation-law fit A exp {-T1/T } from Figure 13 with the typical order parameter ∆typ = exp {⟨ln ∆⟩} at zero temperature for the same set of parameters κ, K and T = 0. The data points are the same as those presented on Figure 13. The value of ∆typ decreases with the increase of κ, thus giving the measure of disorder. The low-disorder case of Figure 7 corresponds to (∆typ/∆0, T1/∆0) ≈ (1, 2) and thus does not fall into the apparent linear tendency.

  F24) with ρ (n) = O w 2h 2n . For instance, ρ (0) = C, ρ (1) = -C w 2h where Li n (z) = ∞ k=1 z k /k nis the the polylogarithm function, and C is the normalization constant. In what follows, we will ignore those corrections.

  φ i ρ i ˆd2 r (∇φ 1 , ∇φ 2 ) (∇φ 3 , ∇φ 4 ) .

  ) e (t ′ ) + ˆdt ′′ δ ⟨I e (t)⟩ δh e (t ′′ ) •

	e ′ ̸ =-e:beg(e ′ )=j	δ ⟨h e (t ′′ )⟩ δA e ′ (t ′ )	A	 e ′ (t ′ ) (j)  .	(B15)

Appendix A: Electrodynamics of a disordered superconductor

In this Appendix, we formulate the low-frequency description of the electromagnetic response of a disordered superconductor, with the particular aim of deriving Eq. [START_REF] Dubouchet | Collective energy gap of preformed cooper pairs in disordered superconductors[END_REF].

Consider the Fourier transform of the Maxwell equations for the electromagnetic potentials (φ, A) in the Coulomb gauge divA = 0:

where ε (ω) is the dielectric permittivity due to the media surrounding the superconductor and the electrons deep within the Fermi surface (for simplicity, we neglect its spatial dependence, it can straightforwardly be restored). All fields here depend on the frequency ω and on the coordinate r. The material response fields ρ, j are to be determined from

and are bound to satisfy the charge conservation law

Here, the kernels Q and R are the material responses of the medium to the vector potential computed as direct variational derivatives:

As per usual, the system (A1-A3) is over-complete because the second equation in Eq. (A2) describing the charge density response is actually satisfied automatically by the true configuration of the electromagnetic fields due to the charge conservation law, Eq. (A3), so we can discard this equation. We can then exclude the charge variable with the help of Eq. A1 and obtain

We now want to calculate the superconducting response to a given external electromagnetic field E ext (r, t), which is created by currents other than those in the material in question. Certainly, this external field automatically satisfies the first four relations in system (A5), but with the current density j produced by some external sources (whose currents are obviously unaltered by the response of the medium). Therefore, to analyze the response of the target medium we have to the substitute the full electric field E full = E + E ext in the last expression describing the current response of the medium. In this way, the Eqs. (A5) are rendered inhomogeneous, so the response problem is well-formulated.

The main difference of the resulting system of equations (A5) with the standard semiclassical description is that due to the inhomogeneity of the superconducting state the current conservation law is not satisfied "automatically" by an appropriate choice of gauge. Instead, the system adjusts the values of microscopic currents and induced electromagnetic field in order to satisfy the charge conservation, making the resulting value of the macroscopic response much more complicated than simply the mean value of the miscroscopic response. In this way, the problem is similar to the one of the macroscopic conductance of a disordered media, where the potential distribution is determined from the charge conservation (or, equivalently, the Kirchhoff's rule for a discrete system).

In particular, η i has complete symmetry w.r.t permutations. We can then expand the sum in Eq. (F32) according to the number of coinciding frequencies:

To describe the cross-Kerr effect by this expression, we pick out only terms that map any two-excitation state to itself. The result reads

We can finally rewrite this as

where

The relevant values of η then read

rendering the final result for the cross-Kerr coefficients K ων :

where the value of α = 3 4π + π 4 ≈ 1.02 originates from the particular form of the nonlinearity of the current-phase relation (F3) in moderately disordered superconductors.

The result (F41) can be qualitatively compared to that for a chain of Josephson junctions with short-range capacitive coupling from Ref. [START_REF] Yu | Kerr nonlinearity in a superconducting josephson metamaterial[END_REF]. In the latter case, the nonlinearity in the current-phase relation is described by sin x ≈ x -x 3 /6 ⇒ α = 1/6, and the factor N in [31, Eq. ( 14)] corresponds to Lw/2ξ 2 .

We finally remind the reader that this result is only applicable to the case of moderately disordered superconductors described by semiclassical approximation. In contrast, the current-phase nonlinearity in strongly disordered superconductors with a pseudogap is yet to be described theoretically.