Proximity-induced gapless superconductivity in two-dimensional Rashba semiconductor in magnetic field

Two-dimensional semiconductor-superconductor heterostructures form the foundation of numerous nanoscale physical systems. However, measuring the properties of such heterostructures, and characterizing the semiconductor in-situ is challenging. A recent experimental study [arXiv:2107.03695] was able to probe the semiconductor within the heterostructure using microwave measurements of the superfluid density. This work revealed a rapid depletion of superfluid density in semiconductor, caused by the in-plane magnetic field which in presence of spin-orbit coupling creates so-called Bogoliubov Fermi surfaces. The experimental work used a simplified theoretical model that neglected the presence of non-magnetic disorder in the semiconductor, hence describing the data only qualitatively. Motivated by experiments, we introduce a theoretical model describing a disordered semiconductor with strong spin-orbit coupling that is proximitized by a superconductor. Our model provides specific predictions for the density of states and superfluid density. Presence of disorder leads to the emergence of a gapless superconducting phase, that may be viewed as a manifestation of Bogoliubov Fermi surface. When applied to real experimental data, our model showcases excellent quantitative agreement, enabling the extraction of material parameters such as mean free path and mobility, and estimating $g$-tensor after taking into account the orbital contribution of magnetic field. Our model can be used to probe in-situ parameters of other superconductor-semiconductor heterostructures and can be further extended to give access to transport properties.

Two-dimensional superconductor-semiconductor heterostructures, where the semiconductor possesses Rashba spin-orbit coupling (SOC) [2], have attracted a lot of attention recently due to several promising applications.These include superconducting qubits [3,4], engineered pwave superconductivity [5,6], and, when subjected to a magnetic field, they present a promising platform for hosting Majorana zero-modes [7][8][9][10].However, experimental investigations of such heterostructures pose challenges.In transport measurements, for instance, the superconductor acts as a shunt.As a result, in order to study properties of the buried interface, some additional experimental probes are required.To address this requirement, recent experiments have probed superfluid density [1,11], vortex inductance [12], and terahertz cyclotron resonance [13].
The measurement of the superfluid density in superconductor-semiconductor heterostructure made out of aluminum deposited on top of spin-orbit coupled two dimensional electron gas (2DEG) using a resonant microwave circuit was implemented by some of the present authors and collaborators [1].The key insight used in the experiment [1] was the qualitatively different response of superfluid density in the conventional superconductor (SC) and proximitized 2DEG to in-plane magnetic field.Specifically, at a certain value of magnetic field Bogoliubov bands in 2DEG touch the Fermi energy, causing an abrupt depletion of superfluid density that was confirmed by the experiment.
Although the superfluid density depletion observed in experiment confirmed the proximityinduced order parameter in spin orbit coupled 2DEG, numerous questions remained open.In particular, to extract specific parameters of 2DEG from experimental data, g-factor, mobility, or carrier density, Ref. [1] relied on the simplified theoretical model that neglected disorder in 2DEG.This model resulted only in a qualitative agreement with the data, thus calling for more realistic theoretical model.Another outstanding question was the fate of the system in the regime when Bogoliubov bands touch the Fermi energy.On the one hand, at this point so-called Bogoliubov Fermi surfaces were predicted to emerge [6].On the other hand, in disorder-free theories the superfluid density contribution from 2DEG becomes negative in that regime, signaling potential instability, that was previously considered in a number of different contexts [14][15][16].Thus, reliable theoretical description of the system in the regime with Bogoliubov Fermi surfaces remains missing.
In this paper we develop a more realistic model for the proximitized 2DEG with strong spinorbit coupling that incorporates the presence of non-magnetic disorder and in-plane magnetic field.To describe the proximitized 2DEG with disorder, we use a Green's function formalism and perform a self-consistent calculation of the self-energy in the 2DEG, that takes into account impurity scattering.Using Green's function, we calculate the density of states (DOS) in the 2DEG and its superfluid density.Our main finding is that disorder leads to a regime of stable gapless superconductivity for a range of magnetic fields that expands with disorder strength.Although in the presence of impurities momentum is not a good quantum number, the gapless regime in presence of disorder may be viewed as a stable generalization of the state with Bogoliubov Fermi surfaces to the case of a disordered material.
We note that the 2DEG with spin-orbit coupling and pairing was previously considered in the literature.In particular, in early works [17,18], the 2D superconductor with spinorbit coupling was studied.In these papers, in contrast to ours, superconductivity is an internal property of 2D layers and does not come from an external superconductor through the proximity effect.Besides, in Ref. [17] the non-magnetic disorder was not incorporated.Later works [5][6][7][8] studied the superconductor-2DEG heterostructure without incorporating non-magnetic disorder and suppression of order parameter by magnetic field.Several studies have incorporated the effects of disorder and band-bending to more realistically model the superconductor-semiconductor interface [19][20][21][22][23]. Also, proximity induced triplet pairing was considered in Ref. [24] but for low carrier concentrations when spin-orbit splitting and Fermi energies are comparable.
In addition to 2DEG proximitized by conventional s-wave superconductors, the physics of Bogoliubov Fermi surfaces and interplay between disorder, spin-orbit coupling and superconductivity is actively studied in other material systems.In particular, Bogoliubov Fermi surfaces in presence of in-plane magnetic field were probed by the scanning tunneling microscopy experiments performed on the proximitized surface states of topological insulator [25].Also, Bogoliubov Fermi surfaces were studied in the multiband superconductors with broken timereversal symmetry [26][27][28][29][30].However, typically in such systems the pairing has higher angular momentum and is expected to be rapidly suppressed by disorder.In a different direction, Refs.[31][32][33][34][35] considered the phase diagram and effect of disorder and magnetic field on the so-called Ising superconductivity, where strong spin-orbit coupling locks spin to a particular direction.
Application of our theoretical framework to the real experimental data from Ref. [1] results in a much better quantitative agreement with the data, compared to the oversimplified theoretical model in [1].Using simultaneous multi-fitting of the model parameters, we extract the values of scattering time, mobility, and carrier density of the 2DEG, and estimate values of the g-tensor anisotropy and g-factor.The mobility is qualitatively consistent with Hall measurements [1], and extracted scattering time suggests that semiconductor has disor-der strength that puts it in between clean and dirty regimes.Moreover, we identify the fairly broad range of magnetic fields where the 2DEG realizes gapless proximity-induced superconductivity, and predict the shape of the density of states that could be potentially probed in tunneling experiments.
Although our model results in a quantitative agreement with experimental data, it still has a number of limitations that are related to our treatment of proximity effect.Specifically, we ignore an inverse proximity effect, and also treat the orbital effect of the magnetic field related to the motion of carriers between superconductor and 2DEG only qualitatively.Incorporation of these effects is important for understanding the ground state of the heterostructure from Ref. [1] for magnetic fields larger than 0.75 T, where our model still predicts that the superfluid turns negative signaling a potential instability.Also, self-consistent treatment of orbital effect of magnetic field is required for a more reliable extraction of the g-factor, that is currently done using phenomenological considerations.We hope that these shortcomings of our model can be addressed in the future work.The present theoretical model may be used as a guide for the values of magnetic field required for realizing the potential instability of the gapless superconducting state and studying it in future experiments.
More broadly, our theoretical model with some modifications will be applicable to a much broader family of materials, such as hybrid semiconductor [36] and topological insulator [37,38] nanowires or proximitized surface states [25], Germanium based 2DEGs [39], ferromagnetic hybrids [11], and two-dimensional materials with strong spin orbit coupling such as transition metal dichalcogenides.Also, it can be extended to predict the dissipative electromagnetic response, spin susceptibility and other characteristics accessible in the future experiments.Such extension of our work could allow a comprehensive in-situ characterization of heterostructures before fabrication of more complicated devices that i.e. attempts to realize Majorana physics [9,10].
The rest of the paper is structured as follows: In Section 2 we introduce a theoretical model based on the Green's function formalism.This section covers the calculation of the density of states (DOS) and superfluid density using Green functions.We also discuss how the magnetic field affects the order parameter and extend our model to consider cases with magnetic anisotropy.Next, in Section 3 we discuss theoretical predictions for the DOS and superfluid density, both with and without suppression of the order parameter by the magnetic field.In addition, we identify the parameter range of magnetic field where the semiconductor has gapless proximity-induced superconductivity as a function of disorder strengths.In Section 4 we introduce and implement a fitting procedure for experimental data for Al-InAs heterostructure [1], discuss the resulting material parameters and show predictions for the density of states.Finally, we conclude in Section 5 with the summary of our results and outstanding open questions.

Theoretical model
In this section we introduce the main theoretical model for describing the 2DEG proximitized with a conventional superconductor.We begin by describing the physical system and the approximations used in our theoretical treatment, covered in Sec.2.1 and Sec.2.2.Following this, the Green's function formalism is introduced in Sec.2.3, where we also discuss selfconsistency conditions.The derivation of physical properties, such as the density of states and superfluid density, is detailed in Sec.2.4.Then in Sec.2.5 we discuss the generalization of our calculations to the case of anisotropic g-factor.Finally, we conclude this section with the review of the depairing theory of the conventional superconductor in Sec.2.6 required for the realistic modeling of the heterostructure.

System and physical assumptions
We study 2D heterostructure which consists of superconducting and semiconducting layers schematically depicted in Fig. 1(a).The superconductor is assumed to be isotropic material with a conventional s-wave order parameter.Moreover, motivated by the common use of aluminum, we assume that superconductor is strongly disordered.The 2DEG layer, in contrast is generally anisotropic and has a strong spin-orbit coupling.For simplicity, we choose the spin-orbit coupling to be of the Rashba form with potential generalizations discussed in Appendix D.Moreover, we also neglect the anisotropy in most of 2DEG parameters, and take into account only anisotropic g-factor that controls the coupling to the external in-plane magnetic field of varying direction and magnitude, see Fig. 1(a).Finally, the new crucial ingredient incorporated in the present work is the presence of disorder of arbitrary strength in the 2DEG layer.
Simultaneous presence of spin orbit coupling and magnetic field breaks both spin-rotation symmetry and time-reversal symmetries of the system.The effect of the in-plane magnetic field on the conventional superconductor is the suppression of the order parameter, that can be described by the conventional depairing theory [40,41] and will be reviewed in Sec.2.6.Likewise, we incorporate the coupling of the in-plane magnetic field to the charge carriers in the 2DEG.However, assuming the small thickness of the 2DEG, we take into account only Zeeman term [42].
The crucial assumption that underlies our treatment is that the 2DEG gets the superconducting order parameter by proximity effect but is not affecting the order parameter in the superconductor.Such absence of inverse proximity effect from the 2DEG onto superconductor (that would result in additional suppression of the order parameter in the superconductor) is justified if the density of carriers is much higher in the superconducting material.The absence of inverse proximity effect along with the assumption of transparent interface between semiconductor and superconductor (we note that our model can be easily extended to include suppression of the induced order parameter due to imperfect interface transparency, in practice it is also possible to characterize induced order parameter in the 2DEG via tunneling measurements) allows us to treat the 2DEG with the induced order parameter that is entirely determined by depairing intrinsic to superconductor.This approximation distinguishes our model from earlier Ref. [18], where the pairing mechanism was assumed to be intrinsic to the 2DEG.

Hamiltonian and band structure
In order to construct the Green's function, we first consider the system in the absence of disorder in the momentum space.To this end, we use the enlarged space that includes spin and Nambu (particle-hole) degrees of freedom, with operators c † k (c k ) that are Fourier transformation of real space creation (annihilation) operators, We write the Hamiltonian as: where the function ξ k = k 2 /(2m)−µ encodes the parabolic band structure with effective mass m and chemical potential µ.Parameter λ so determines the strength of the Rashba spin-orbit coupling, and we used approximation λ so k ≈ λ so k F (k F is a Fermi momentum), that assumes that spin orbit splitting is much smaller than the Fermi energy.Set of Pauli matrices σ x, y,z acts in the spin space and Pauli matrices τ x, y,z operate in the particle-hole space.We introduced φ k as the angle of momentum direction with respect to x-axis, cos φ k = k x /k.Finally, terms in the second line of Eq. ( 2) describe the Zeeman energy and the proximity-induced isotropic order parameter.
In the definition of Hamiltonian (2) we fix the direction of magnetic field to point along y-axis as in Fig. 1(a).This results in the Zeeman energy −σ y V with where g is the isotropic g-factor (generalization to anisotropic g-tensor is presented in Sec.2.5 below) and µ B is the Bohr magneton.We emphasize the presence of a factor 1/2 in Eq. ( 3) that was omitted in Refs.[1,6] but is essential for comparing the values of g-factor with the previous literature.Another consequence of the in-plane magnetic field, not included in Eq. (2), is its orbital effect, which is known to introduce an additional phase between carriers in the superconductor and 2DEG [43].While incorporating this effect into our self-consistent treatment is beyond the scope of the present work, we take it into account phenomenologically in Sec.4.3 in order to extract the realistic values of g-factor.
The left panel of Figure 1(b) shows a sketch of the spin-orbit locked Fermi surfaces corresponding to the Hamiltonian (2) without pairing and magnetic field.Non-zero induced order parameter, ∆, without magnetic field V = 0 results in the isotropic gap opening for each of the two spin-momentum locked Fermi surfaces, see the sketch in Fig 1(b) on the right.Non-zero magnetic field, V > 0, results in the tilting of the Bogoliubov's bands.Although formally the distance between top and bottom of quasiparticle band remains 2∆ (for now we neglect the depairing effect of magnetic field on superconductor), the separation from the hole branch and Fermi level decreases for k x > 0 and increases for the negative k x < 0, see Fig 1(b).At the special value of the field, when V = gµ B H/2 = ∆, the Bogoliubov's bands touch the Fermi level, signaling emergence of Bogoliubov's Fermi surfaces [6].
The qualitatively different form of the band structure and associated density of states (DOS) is illustrated in Fig. 1(c)-(d) for V < ∆ and in panels (e)-(f) when V > ∆.In the latter case, the DOS has no gap anymore -a direct manifestation of the emergence of Bogoliubov's Fermi surfaces.We note, that although the analytic expression for the band structure resulting from Eq. ( 2) was obtained in the previous literature [6], the associated DOS was studied only numerically.In the Appendix A we present the analytic derivation of the DOS for this band structure.We show that the standard square-root singularity in the DOS at ∆ is split due to presence of magnetic-field and spin-orbit coupling into two Van Hove singularities: at the band bottom the DOS has a discontinuous jump at energies ±|∆ − V | rather than square root divergence as in BCS case.In addition, at energies ±|∆ + V | the band structure develops a logarithmic Van Hove singularity, see Fig. 1(d) and (f).
We are interested in the case when the energy scale associated with the spin-orbit coupling much larger compared to Zeeman energy, λ so k F ≫ V .In this case, the band-splitting due to spin-orbit coupling is much larger compared to the Zeeman-induced anisotropic shift of energy bands.This allows us to neglect the inter-band terms induced by the magnetic field, as they are suppressed by the small parameter [V /(λ so k F )] 2 .More details on this approximation are presented in the Appendix A.

Green function formalism in presence of disorder
We incorporate the non-magnetic disorder using the self-energy in the Green's function formalism.The matrix form of disorder potential in coordinate representation reads h dis (r )=τ z U (r ), where U(r ) = i u 0 δ(r − r i ) and u 0 and r i are the impurity strength and coordinate location respectively.We calculate the average self-energy Σ due to scattering of electrons by the impurity potential, here and in what follows we use the following notation: k ...
.., and denote by n imp the concentration of non-magnetic impurities.The average self-energy enters into the expression for the Greens function in presence of disorder, where the first two terms correspond to the inverse of the Green's function without disorder defined with Hamiltonian (2), and ε n = 2πT (n + 1/2) correspond to fermionic Matsubara frequencies.

Density of states and superfluid density
The DOS can be calculated using its expression via Green function and explicit parametrization of the latter in Eq. ( 7): where τ is the mean scattering time and ν 0 is the density of states at each of the spin-orbit split Fermi surfaces.Here we used the retarded Green function which is analytical continuation of Matsubara Green function G from imaginary frequencies to real ones: Analogously, in all calculations pertaining to the DOS below, we also redefine the parameter iε (1) → ε (1) .The second relevant observable -superfluid density n s -is calculated using the electromagnetic response kernel.The response kernel Q αβ (q , iω n ) where ω n = 2πT n (we use units with k B = 1) determines the current in response to the applied electromagnetic field A β , where α and β are two-dimensional spatial indices.Superfluid density is obtained as a zero frequency and zero momentum limit of the response kernel [47,48]: In the Appendix C we derive the expression for the response kernel using Green functions, where the velocity tensors v α (α = x, y) are defined as vx = k x /m and vy = k y /m.Note, that here we neglected the contribution from spin-orbit coupling to velocity operator, since it is suppressed by the parameter λ 2 so /v 2 F , where v F is the Fermi velocity.Substituting the Eq. ( 13) into ( 12) and using the explicit form of the Green's function (7) we get that the tensor of superfluid densities is diagonal, with two non-zero components n s,x x and n s, y y .The superfluid density x x-component is denoted as n s,⊥ (since magnetic field points in y-direction) is given by the following expression The expression for the y y-component denoted as n s,∥ can be obtained by replacing the cos 2 φ k with sin 2 φ k .We note, that the spin-orbit interaction strength λ so does not explicitly enter the expressions for superfluid density and DOS due to change of variables of integration.We use ξ that is the relative energy difference to the respective spin-orbit split Fermi surfaces and neglect the difference in DOS between two Fermi surfaces assuming that spin-orbit splitting is much smaller compared to the Fermi energy.At the same time, presence of spin-orbit energy splitting that is much larger than induced pairing and Zeeman term (we assume that renormalization due to disorder does not change the order of magnitude of parameters V (1) , ∆ (1,2) ) is important for our treatment, since this allows to neglect inter-band terms.

Arbitrary magnetic field direction and anisotropic g -tensor
Until now we considered the magnetic field pointing along y-direction and did not take into account possible g-tensor anisotropy.However, g-tensor anisotropy is generically expected in the presence of spin-orbit coupling.Indeed, in the relevant case of spin-orbit coupled asymmetric quantum wells, in-plane g-factor anisotropy is known to be a large effect [49][50][51].
Here we generalize our previous results to include both of these additional ingredients.
Assuming general in-plane magnetic field and arbitrary g-tensor, ĝ = g x x g x y g y x g y y , the Zeeman contribution to energy is given by the convolution of Pauli matrices and magnetic field, µ B σ α ĝαβ H β /2.Explicit convolution gives the expression −σ x V x − σ y V y for the Zeeman energy that was previously given by −σ y V in Eq. ( 2).However, now the parameters V x, y are defined as: where χ is the angle between magnetic field H and x-axis.Using these notations, we obtain the more general form of the Hamiltonian (2) in momentum representation, We can reduce this Hamiltonian to the situation considered previously by unitary transformation acting in the spin and Nambu spaces.To this end, we define the angle θ as Rotating the Hamiltonian using unitary matrix U θ we bring it to the same form as before, Eq. ( 2), however with the shifted angle φ k and redefined where we introduced Using this insight, we reexamine our calculation of physical quantities in Sec.2.4.The angle rotation does not affect the calculation of DOS, since the shift of angle φ k does not change the integral over φ k .Hence, the DOS in presence of g-factor anisotropy and general direction of magnetic field can be obtained from Eq. ( 10) using the generalized expression for V from Eq. ( 19).
The extension of the calculation of the superfluid density is more involved, since the shift of the angle φ k affects Green's functions, but does not affect velocity operators in Eq. ( 13).Careful incorporation of such shift in Eq. ( 14) gives the following expression for the x xcomponent of the superfluid density that relies on the expressions n ⊥ (H) and n || (H) introduced in Eq. ( 14), Note, that in the case of anisotropic g-tensor, the tensor n αβ is not diagonal anymore.Its remaining components can be calculated in the same way as n x x by incorporation of the shift of angle φ k in the Green function.In detail, for x x-component we performed the shift of the angle φ k in the Eq. ( 14) everywhere except the factor cos φ k which corresponds to the velocity operators.For the y y and x y-components one may use the same equation with cos 2 φ k replaced by sin 2 φ k and sin φ k cos φ k , respectively.This gives non-vanishing n x y = n y x component of the superfluid density, In particular, for two orthogonal orientations of magnetic field, H∥x and H ⊥ x we get where Figure 2: (a-c) Superfluid density normalized by the total electron density, n s /n as a function of magnetic field for three progressively increasing values of disorder.With increasing disorder the purple areas, representing the regions of magnetic field where there is a gapless regime, expand, and the decline of the superfluid density with magnetic field becomes smoother.In addition, the point where superfluid density becomes negative, so that our treatment is not applicable anymore (shown as dashed lines) is close to the clean case value V = ∆ in panel (a) and is shifted to progressively larger values of magnetic field.(d-e) DOS normalized by its value in the normal state, ρ/ν 0 for the two different values of magnetic field indicated by vertical lines in the upper panels.The larger value of the magnetic field puts the 2DEG into the gapless regime, as is clear from the non-zero value of the DOS at zero energy.For the calculation of superfluid density, the temperature was set to T /∆ = 0.0564, and suppression of the order parameter by magnetic field is ignored in all plots.
From here we see that off-diagonal terms of g-tensor induce off-diagonal contributions to the superfluid density.Such non-diagonal g-tensors naturally arise in the anisotropic quantum wells of interest here [50].In spin-3/2 hole systems within low-symmetry quantum wells (e.g., GaAs), the g-tensor can exhibit asymmetry due to the interaction between the p-like character of hole wave functions and asymmetric bandedge profiles, alongside specific crystallographic orientations [52].Equations ( 23)- (24) suggest that measurements of the offdiagonal components of the superfluid density for different orientations of the magnetic field can be used to probe asymmetry between off-diagonal elements of the g-tensor.In particular, the extremely non-symmetric form of the g-tensor, with i.e. g y x ≫ g x y would lead to the signature n x y (H, 0) ≫ n x y (H, π/2) potentially observable without quantitative fitting.

Review of depairing theory in conventional superconductor
One of the main approximations of the current theoretical model is the assumption that the order parameter ∆ in the 2DEG is the same as in the superconductor.Thus although the main focus of our work is the theory of 2DEG, we need to know the superconducting order parameter at non-zero field and temperature.The in-plane magnetic field suppresses the order parameter in the superconductor according to the conventional depairing theory [40,41] that is reviewed below for the sake of completeness.
To find the dependence of the superconductor order parameter on the magnetic field (still parametrized by V ) and temperature, ∆(V, T ) we assume that the superconducting layer is thick enough (unlike the semiconducting one) to neglect Zeeman contribution and leave only the orbital contribution of an in-plane magnetic field [42].At the same time we note that the thickness of the superconducting layer still has to be much smaller than the London penetration depth since we neglect the Meissner effect.In this case, the general approach of pair-breaking treatment in dirty superconductors is applicable [40][41][42].This approach establishes the following connection between the critical temperature in absence of pair-breaking, T c0 , and the critical strength of pair-breaking α at fixed temperature T : Here ψ -is the digamma function, and the parameter α can incorporate contributions from various pair-breaking mechanisms, such as magnetic impurities, magnetic field, electric current, etc.Using the Eq. ( 26) for a fixed temperature T , we can find the value α c (T ) -the value of pair-breaking parameter at which order parameter becomes equal to zero at temperature T .Recalling that in our case the physical nature of pair-breaking parameter α is magnetic field and that orbital contribution from in-plane magnetic field in thin film (the thickness of film is much smaller than coherence length) results in quadratic contribution [40]: α ∝ V 2 , we can express the dependence of α on magnetic field in the following way: where V c is the critical magnetic field for the given temperature T and is determined by intrinsic properties of specific superconductor.We notice that at zero temperature α c (T = 0) has the following form: 2) , where ∆ 00 -order parameter at zero temperature and zero magnetic field, T c0 -critical temperature at zero magnetic field and we used the well-known from BCS theory relation [53].In the next Section we will consider a finite but very small temperature, T = 0.0564∆ 00 so that α c (T ) = 0.49∆ 00 is very close to zero-temperature critical field, but at the same time small finite temperature provides a natural regularization for numerical calculations.Knowing the explicit dependence (27) of pair-breaking parameter α on magnetic field, we can find the order parameter ∆(T, V ) from the following equations [41], Here f n is the Eilenberger Green's function determined using numerical solution of Eq. ( 28) and substituted into Eq.( 29) to obtain the self-consistent value of the gap at non-zero temperature and magnetic field.Knowing the self-consistent gap value and Eilenberger Green's function f n also allows for calculation of the superfluid density as [1,41] where τ SC is a mean free path in superconductor that is assumed to satisfy τ SC ∆ ≪ 1, and n is the density of electrons in superconductor.The last equation for superfluid density n s of SC is not used in the next Section, because we analyze only superfluid density in the proximitized 2DEG.However, this expression will be used in the Section 4 where fit the experimental data, since experiment is probing the total superfluid density of 2DEG and SC.

Results for proximity-induced superconductivity in disordered 2DEG
In this section we study the influence of disorder on the proximity effect in the 2DEG using methods introduces in the previous section.First, we describe the numerical procedure used for the calculation of the self-consistent Green's function.Then we present our results for DOS and superfluid density.Finally, we focus on the region of gapless superconductivity and discuss the range of magnetic fields and disorder strengths where it occurs.We note that in the first two subsections we do not consider the suppression of order parameter ∆ due to the magnetic field, and include it only in Sec.3.3.

Numerical solution for Green's function
In order to find the renormalized parameters and DOS we numerically solve the system of self-consistent equations ( 8) using iterations: we substitute an iteration n − 1 in the right side of an equation to obtain the value of the corresponding renormalized parameter at the next iteration, n, on the left side.For the calculation of DOS we need Green's function for real values of energies, which is why we make the analytic continuation in Eqs. ( 8) iε n → ε+i0 and iε (1) → ε (1) in the same way we did in Eq. (10).As an initialization of the iterative procedure, we use the following values of parameters: ε (0) = ∆, and ∆ (2) (0) is motivated by the fact that it reproduces DOS in normal state after substitution into Eq.( 10)].We use from 100 to 200 iterations in order to get convergent results, the stronger the disorder is, the more iterations were required.We note that iterative procedure can be performed independently at each energy separately.As a maximum value of the considered energy we use: ε max = 6∆, so, ε ∈ [−ε max , ε max ], since no additional feature were observed beyond this range.In addition, we discretize the values of ε and V using the grid with the step of δ = 0.1∆, and obtain functions at intermediate values of these parameters by interpolation.
For the calculation of superfluid density, we need to find Green function for imaginary frequencies iε n using the Eqs.( 8) and the same iterative procedure.In that case, we use the following values of parameters to seed the iteration procedure: iε (0) = V and ∆ (2) (0) = 0.For imaginary energies, numerical calculations converge faster, hence we use from 20 to 50 iterations.The temperature at which calculations are performed is T = 0.0564∆.Since temperature is finite, Matsubara energies are discrete, so, discretization is introduced only for the Zeeman energy V with the step δ = 0.1∆.The maximal value for the considered energy is ε max = 40 • 2πT .

DOS and superfluid density
Using the iteration procedure described above, we calculate the self-consistent Green's function across a range of disorder strengths.Figure 2 illustrates how the superfluid density is suppressed by magnetic field for three different values of disorder (zero-field values of superfluid density are in agreement with well-known analytical results [47]).At the same time, bottom panels of this figure illustrate the behavior of DOS for two typical values of magnetic field, the first being in the regular and second in the gapless regime (see discussion below).
At the smallest value of disorder, shown in Fig. 2(a) and (d), it has a weak effect on the superfluid density and behavior of DOS.Indeed, comparing Fig. 2(d) with Fig. 1(d) and (f) for the DOS in clean case we see qualitatively similar behavior.Moreover, the superfluid density n ⊥ in Fig. 2(a) shows a sharp decline to zero near the value of magnetic field that is only slightly larger then V = ∆, predicted to give rise to Bogoliubov's Fermi surfaces in the clean case [6].We note that nearly at the same time as the DOS becomes non-zero at ε = 0, the n ⊥ component of the superfluid density turns negative.This makes the theory used here inapplicable, as is shown by dashed lines in Fig. 2(a).Thus, the gapless regime -situation where 2DEG has non-zero proximity-induced order parameter, but no gap exists in the density of states -may be realized only in the extremely narrow range of magnetic fields at weak disorder.This will be discussed in more details in Sec.3.3 below.When disorder strength is increased so that scattering time becomes twice shorter than inverse gap, its effect becomes more pronounced.Figure 2(b) and (e) reveal that even at zero magnetic field, V = 0 the superfluid density is already suppressed by about 50% compared to the total carrier density.In addition, the decline of the superfluid density becomes more gradual and peaks in the DOS at energies where clean system has Van Hove singularities become much smoother.Most importantly, the gapless superconductivity is now realized in a narrow window of magnetic field near the value of V ≈ 1.5∆.Remarkably, the disorder not only "smears" all the singularities in the energy space, but also pushes the onset of rapid decline of superfluid density to significantly larger values of magnetic field and creates a gapless regime.
Finally, increasing disorder by additional factor of five in Fig. 2(c) and (f) results in a broad gapless regime that occurs for V ⪆ 2.5∆ and extends beyond the maximal value of V = 3∆ considered in this work.DOS in Fig. 2(f) has extremely broadened peaks at the energy of Van Hove singularities in the clean case and develops a broad dip at zero energy upon entering the gapless regime.Remarkably, the onset of negative superfluid density is pushed beyond the considered range of magnetic fields.This may be understood as an effect of renormalization of effective magnetic field by disorder as we discuss in Appendix F.
The effect of disorder on the system is encoded in the self-consistent solution for the Green's function.In Appendix F we consider the results for renormalized singlet and triplet components of the order parameter, ∆ (1) and ∆ (2) , that become now non-trivial functions of energy.In particular, we show that the disorder-induced triplet component has real part that is odd in frequency [54].In addition, we also calculate the data for superfluid density and DOS including the order parameter suppression in superconductor in Appendix G. Upon comparing Fig. 10 from this appendix with the data in Fig. 2, we conclude that while the order parameter's suppression does not introduce new qualitative features, it does result in the expected vanishing of the superfluid density at the superconductor's critical field.

Region of the gapless superconductivity and potential instability
It is clear from the Fig. 2 that the extent of the gapless regime increases with disorder strength.In order to quantify this phenomenon, we numerically calculated DOS ρ(ε) at ε = 0 for various values of disorder and magnetic field using the discrete grid with the step of δ = 0.02∆ (or δ = 0.02∆ 00 when we incorporate the order parameter suppression) for magnetic field.For a specific value of disorder strength, with the increase of the magnetic field from 0, initially, ρ(0) takes values much smaller than 1.At a certain (discrete) value of the magnetic field, ρ(0) abruptly increases and becomes of the order of 1.We identify this value as the magnetic field at which a gapless regime emerges.For the numerical calculation we use the same approach discussed in Sec.3.1.However, at the boundary between gapless and standard proximity induced superconducting regime more iteration -from 100 to 10000 -are required for convergence.
In Fig. 3(a) and (b) we present the area of parameters (in 1/τ and V axes) which corresponds to the gapless region.In Fig. 3(a) the suppression of order parameter by magnetic field is not taken into account while in Fig. 3(b) we apply the standard depairing theory reviewed in Sec.2.6 for fixed low temperature T = 0.0564∆ 00 and critical field of superconductor being V c = 2.5∆ 00 , where ∆ 00 is the order parameter in superconductor at zero field and zero temperature.The solid purple line on the bottom of the region represents the values of τ and V at which the superconducting gap disappears, the dashed line on the top of the regionvalues of τ and V beyond which superfluid density becomes negative.As we discuss below in greater detail, the negative superfluid density signals that the current theoretical model is inapplicable.Nevertheless, the region of gapless superconductivity is vastly expanding with increasing disorder.
The prediction of the negative superfluid density may be viewed as a results of oversimplification that neglects inverse proximity effect while assuming perfectly transparent interface between 2DEG and superconductor.A more complete model should account for the inverse proximity effect, possibly leading to a reduction in critical magnetic field and change of order parameter suppression, thereby addressing the issue of negative superfluid density.To this end, incorporation of the free energy stemming from both 2DEG and superconductor is necessary.The self-consistent equations for the order parameter resulting from such total free energy would incorporate the inverse proximity effect.In addition, it would also allow to check if the state with spatially modulated order parameter of the form ∆ (r) = ∆ s + ∆e iq•r with |∆ s | ≫ | ∆| is more favorable compared to the uniform state.We should specify, that the leading contribution to the order parameter cannot depend on coordinate because it comes from the superconducting layer which is strongly disordered, therefore, FFLO [55,56] phase is suppressed [18], besides, there is no spin-orbit coupling in the superconducting layer, thus, long-wave helical phase [18] also is not expected to appear.Meanwhile, the small correction comes from the 2DEG due to an inversed proximity effect and can depend on the spatial coordinate because the disorder is not assumed to be strong and, in addition, spin-orbit coupling is present.This model significantly differs from the model in Ref. [57] where the dependence of order parameter on coordinate was defined entirely by the 2DEG properties and as a result, there was no constant contribution from the superconducting layer.
The development of a more comprehensive theory discussed above is beyond the scope of this paper and could be an avenue for future research.At present, it is unclear if developing more realistic model will just quantitatively alter the predictions in vicinity of point where superfluid density turns negative, or if a non-uniform superconducting state may emerge above the dashed line in Fig. 3. Nevertheless, the prevalence of negative superfluid density diminishes with increasing disorder (as is seen on Fig. 3), making the current model applicable in a wider range of magnetic field.Moreover, despite the mentioned limitation, the current theory offers significant predictive power, as will be demonstrated in the next section by fitting experimental data.
Figure 5: The theoretical model that incorporates disorder (solid lines) results in a much better agreement with the experimental data (dots), compared to the model [1] that does not incorporate non-magnetic disorder (dashed line).The data in panel (a) corresponds to the resonance frequency suppression with field H x at fixed T = 0.1 K. Panel (b) shows the temperature dependence of the resonant frequency at H = 0 T.

Fitting experimental data
In this section, we apply the theoretical model described above to the experimental data on Al-InAs heterostructure reported in Ref. [1].Since the experiment probes the total superfluid density, first in Sec.4.1 we fit the dependence of pair breaking parameter in aluminum.Next, in Sec.4.2 we discuss the fitting procedure that we use to fit the experimental data and extract material parameters.Finally, in Sec.4.3 we discuss the extracted values of parameters and show predictions for the DOS that can be checked in future experiments.

Pair-breaking effects in Al
In order to understand the dependence of the order parameter ∆ in Al on temperature T and magnetic field H we apply the depairing theory reviewed in Sec.2.6.Specifically, we use the dependence of the critical temperature on magnetic field measured in Ref. [1] to determine the pair-breaking parameter α as a function of magnetic field.In contrast to Ref. [1] that extracted the depairing parameter as a phenomenological function of magnetic field, here we aim to capture the dependence of depairing parameter by a simple function thereby providing insights into depairing mechanisms in aluminum.Moreover, we emphasize that although experimental data from Ref. [1] measures the critical temperature of the 2DEG-superconductor heterostructure, we ignore the effect of the 2DEG onto superconductor.This is consistent with the approximation of neglecting the inverse proximity effect adopted in this work.In Sec.2.6 we discussed that the orbital contribution from in-plane magnetic field results in the quadratic dependence α ∝ H 2 .However, assuming purely quadratic dependence of depairing on the field does not result in a good quantitative agreement with the experimental data (see dashed line in the Fig. 4).This discrepancy may be consistent with the tiny outof-plane component of magnetic field that, despite careful alignment may be present in the experiment, see Appendix E. Alternatively it may potentially emerge from the influence of magnetic field on the inverse proximity effect that is not considered in the present work.The out of plane component of magnetic field contributes linearly to the depairing parameter α [42,58,59].Thus we use the following ansatz for the depairing parameter where two dimensionless constants ζ 1,2 parametrize the entire dependence.Here, T c0 is the critical temperature in the absence of the field, and H cr is the zero-temperature critical magnetic field.
We can establish additional relation between constants ζ 1,2 introduced above.For this we consider the vicinity of critical field H → H cr .For such values of magnetic field, the critical temperature tends to zero T c → 0 and we can simplify the relation (26) from Sec. 2.6 by expanding it assuming that T c → 0. This expansion gives us the following relation Comparing this with Eq. ( 31 Using specific dependence of α on magnetic field in Eq. ( 31), we calculate the dependence of T c on applied magnetic field according to Eq. ( 26).The resulting curve is fitted to experimental data on the dependence of critical temperature on magnetic field using the least square method to determine the best value of ζ 1 .The fit has excellent agreement with experimental data, as shown in Fig. 4, yielding the following value of ζ 1 , Knowledge of the explicit dependence of pair-breaking parameter α on magnetic field H allows us to find the order parameter ∆ (T, H) and superfluid density in Al using Eq. ( 28)- (30).Note, that we take τ SC ∆ 00 = 0.001 for Al, which corresponds to the dirty limit, but the functional form of all dependences is not affected by specific value of τ SC as soon as it is sufficiently small.

Fitting procedure
After determining the depairing in Al, we turn to fitting the data on the superfluid density measured by the experiment [1].The experiment probes the Al-InAs heterostructure by placing it in resonator and sensitively measuring the resonance frequency f r .This frequency has a constant "geometric" contribution denoted as f geo , and also receives a contribution from superconducting condensate in aluminum and 2DEG, denoted as f kin .The total resonance frequency measured in experiment reads [1]: Assuming the distributed circuit model, Ref. [1] suggested that the kinetic contribution to inductance may be viewed as a sum of 2DEG and superconductor contributions, Here we wrote a combination of two individual contributions from 2DEG and superconductor in a slightly different form, and emphasized the magnetic field and temperature dependence.Both contributions depend on the constants c p,s that have dimension of squared frequency and will be our fitting parameters.For the 2DEG the constant c p multiplies the dimensionless ratio between the superfluid density in 2DEG and normal carrier density.In contrast, for superconductor, the constant c s multiplies the superfluid density at a given field and temperature, normalized by the superfluid density at zero temperature and field, n (Al) s (0, 0) = π(τ SC ∆)n (Al) , that is suppressed by a small factor τ SC ∆ compared to the full carrier density n (Al) in the dirty limit.In both cases, the normalized superfluid densities will be provided by the theoretical model developed above.
After we obtain the value of parameters c p and c s from fitting, we may use them to extract material parameters.To this end we use on the following relation between these constants and material parameters, developed in Ref. [1]: Here l denotes the resonator length, C is the resonator capacitance, γ is a geometric factor characterizing the resonator and given in Eq. (S4) in the supplement of Ref. [1], m (InAs) is the effective mass of a quasiparticle in InAs, e is the electron charge, ∆ 00 is zero-field zerotemperature order parameter in Al, and □ is the aluminum sheet resistance in normal state.From here it is clear that knowledge of constants c p,s will yield an estimate of the carrier  36).
The fitting procedure used to match experimental data relies on Eqs. ( 34)- (35).These equations explicitly contain three fit parameters, c p , c s , and f geo .The theoretical model for the ratio n (Al)  s (H, T )/n (Al) s (0, 0) is entirely fixed by Eqs. ( 31) and ( 33) in previous section and contains no additional free parameters.In contrast, the normalized superfluid density of 2DEG additionally depends on the scattering time τ intrinsic to 2DEG and components of g-tensor, g eff x x and g eff y y (remaining components are assumed to be zero).Thus in total the fit between theoretical model and experimental data is performed over six parameters c p , c s , f geo , τ, g eff x x , and g eff y y .Note, that we denote the g-factors in the fitting as effective ones, since in what follows we phenomenologically incorporate the orbital effect to estimate their physical values.
In practice, to achieve such six-parameter fit we simultaneously fit three experimental dependencies as illustrated in Fig. 5 and Fig. 6. the dependence of f r on T at H = 0 T, the dependence of f r on H at T = 0.1 K and the dependence of f r on the angle between magnetic field and x-axis for six different values of magnetic field listed in Fig. 6.The simultaneous fit was conducted in the following manner: we constructed the cost function for each observable F α as S α = i F (th)  α , where F (th) α (x) is theoretical prediction for the observable F α and {F  x} ≡ { f r (H, χ), (H, χ)} for magnetic field and its direction sweep in Fig. 6.Then, to derive the total cost function, we sum these functions with specific weights w α .The weights are selected in such a way that at the minimum of the total cost function, contributions from different cost functions were approximately of the same order: S tot (c p , c s , f geo , g x x , g y y , τ) = w 1 S 1 + w 2 S 2 + w 3 S 3 .By minimizing this cost function, we extracted the best values of fitting parameters.
To estimate errors for these fitting parameters, we use a procedure similar to bootstrapping [60].Specifically, we construct a set of cost functions {S ( j) tot } in the following way: for each observable, we select only a portion of the experimental points in a random manner {F (exp) α (x), x} ( j) and then proceeded with the fitting in the same way as described above.This yields a set of different values for the fitting parameters {c p , c s , f geo , g eff x x , g eff y y , τ} ( j) .We check the resulting distribution of fit parameters and estimate error bars from its variance.Figure 7: The material parameters extracted from the best fits to experimental data are used to generate predictions for superfluid density of 2DEG and DOS for two different orientations of magnetic field.Purple areas correspond to the regions of the magnetic field where the gapless regime is realized.Vertical lines on the plots for superfluid densities in (a)-(b) show the values of the magnetic field at which DOS is plotted on the respective panels (c)-(d), with green line illustrating the regime with a gap, while the blue line corresponding to the field when system is in the gapless regime.

Material parameters and prediction for DOS and gapless regime
The fitting of the theoretical model to experiment using the procedure described above gives the value of parameters that are summarized in Table 1.Figures 5-6 demonstrate the good quantitative agreement of best fits from the theory model with the experimental data.In particular, the agreement is much better compared to the model where disorder is not incorporated from Ref. [1], compare dashed and solid lines in Fig. 5. Also, due to presence of 1/2 in the definition of V , Eq. ( 3), the values of g-factor are now approximately two times larger.
Out of six fit parameters in Table 1, g-factors have the most immediate physical interpretation.However, the g-factors values of g eff x x = 30.8± 0.8 and g eff y y = 17.8 ± 0.8 are considerably larger compared to the ones used in the literature, where a range of values 3-11 was reported [61][62][63][64].We note, that comparable values of g-factor can be obtained from the qualitative model in Ref. [1] that ignores disorder, provided one uses definition Eq. (3) for Zeeman energy.Such overestimation of g-factor points out that the orbital effect of the in-plane field cannot be neglected.For instance, recent work [43] reported similar overestimation of g-factor due to ignorance of orbital effects in a Josephson junction made out of similar superconductor-2DEG heterostructures.
While rigorous incorporation of the orbital effect of the magnetic field is reserved for the future work, we take this effect into account phenomenologically.Specifically, we assume that 2DEG is separated by a distance d 2DEG from the superconducting layer.Then taking into account the phase accumulated to in-plane magnetic field results in the following induced order parameter in InAs, ∆ InAs (r) = ∆ Al • e i2φ (r) , where φ(r) = q • r, with |q | = πH d InAs /Φ 0 , where Φ 0 is the magnetic flux quantum [57].Proceeding in the same way as in Ref. [57], we derive that effect of such non-zero momentum q can be incorporated by changing the value of the parameter V into V (eff) with where ± sign corresponds to different Rashba-split bands and we neglect the contribution λ so q since λ so ≪ v F is assumed throughout this work.We can refine our estimate of the g-factor by using the right hand side of Eq. (37).Equating expression in parenthesis to the fitted value of g-factor, we can extract its physical value according to expression: The extraction of g-factor in such phenomenological way strongly depends on the distance between 2DEG and superconductor.For instance, assuming distance d 2DEG ∼ 0.2 nm we obtain 40% smaller value of g x x ∼ 18 and even more suppressed value of g y y ∼ 5. Larger distance of d 2DEG ∼ 0.5 nm would result in a nearly vanishing value of g x x .While these values of thickness are approximately an order of magnitude smaller compared to the estimated distance between aluminum and 2DEG layer, the band bending can decrease the effective distance.Also, our phenomenological considerations do not treat the orbital effect self-consistently, thus overestimating its magnitude.Although orbital effect for such thickness is (within our phenomenology) nearly equivalent to the Zeeman energy contribution, it does not involve the electron spin and hence is expected to be isotropic.Thus, we expect that quantitative incorporation of the orbital effect may allow for extraction of effective thickness as well as g-factors using the fitting of the data rather than phenomenological considerations.
We note that treatment discussed above is phenomenological, and proper incorporation of such orbital effect and disorder requires revisiting our self-consistent treatment presented above.In particular, we expect the renormalization due to disorder to suppress the orbital effect.Additional subtlety is related to the sign of the g-factor that is know to be negative in InAs, and to the relative ± sign in Eq. (37).Since our treatment is not sensitive to the relative sign of the g-factor, in order to obtain Eq. ( 38), we assume g > 0 and also choose the plus sign in Eq. ( 37) corresponding to the Fermi surface where gap closes earlier.
Next, we discuss the scattering time in 2DEG that is also extracted from the fitting.The value τ∆ 00 = 0.23 suggests that InAs heterostructure is in intermediate between clean and dirty regimes, with its mean free path being of the order of 1.3 µm, assuming λ F = 10 nm.Estimating mobility from the scattering time we obtain where for the effective mass of electron in InAs we used m (InAs) = 0.04m e [64], m e is electron mass, and e is the charge of electron.The InAs mobility, measured in the absence of Al, is 1.3 • 10 4 cm 2 /(V • s) [1], suggesting that the presence of Al does not drastically alter the InAs carrier mobility.Finally, the parameters c s , c p , and f geo provide consistency check.Parameter f geo is consistent with the expected value based on electromagnetic simulations in the Ref. [1].Knowing the parameters c s and c p we can calculate the aluminum sheet resistance R (Al) □ and carrier density in InAs, n (InAs) , using Eq. ( 36): The resulting value of the aluminum sheet resistance is approximately the same as in Ref. [1].However, the InAs carrier density is twice as large.The different estimate arises due to suppression of superfluid density by disorder that was not incorporated in the theoretical model used in Ref. [1].This omission resulted in a lower value of c p compared to ours, which in turn led to lower values of the carrier density.Comparing the fit carrier density with the Hall value n (InAs) = 1.06 • 10 12 cm −2 , measured in the absence of Al, we find that the InAs carrier density is dramatically increased by the presence of Al.
Using the material parameters from the Table 1 we generate predictions for the DOS and superfluid density in proximitized InAs and check whether gapless regime is accessible.Fig. 7 shows our predictions for two orientations of magnetic field.When magnetic field is parallel to the x-axis we see that the system is predicted to have a wide gapless regime.In contrast, for the field aligned with y-axis we see the sharper decline of the superfluid density, gapless regime occurs at larger values of magnetic field, and has smaller extent.This is explained by the anisotropic response of superfluid density in 2DEG due to presence of spin orbit coupling.In principle, the dependence of DOS on energy can be studied experimentally using a scanning tunneling microscope.Thus, the predictions presented in Fig. 7 can be experimentally verified.
Finally, let us discuss the uncertainties in determining the parameters.We note, that although error bars estimated by bootstrapping in Table 1 are relatively small, the systematic uncertainties that come from potential further simplifications present in the model, such as phenomenological treatment of orbital effect, assumption of perfectly isotropic mass, consideration of only Rashba spin orbit coupling, diagonal form of g-tensor, and assumption of fully transparent interface between 2DEG and superconductor are potentially much larger.While relaxing these assumptions is straightforward, additional data is needed to prevent overfitting and enable reliable extraction of additional material parameters.

Summary and Outlook
To conclude, we examined a heterostructure consisting of 2DEG with strong spin-orbit coupling proximitized by superconductor and subjected to an in-plane magnetic field.Our theoretical model incorporates non-magnetic disorder of arbitrary strength in the semiconducting layer, a feature that was not addressed in previous studies.Our model gives predictions for the density of states and the superfluid density as a function of the magnetic field and varying disorder strength.We observe that increasing disorder stabilizes a gapless superconducting phase within a progressively increasing range of magnetic fields.This regime may be viewed as an extension of the phase with Bogoliubov Fermi surfaces that incorporates disorder in the system.
We applied our theoretical model to experimental data for an Al-InAs heterostructure obtained in Ref. [1].Our model that incorporates disorder was able to quantitatively describe the data, also enabling us to extract parameters of the InAs layer, such as the anisotropic g-tensor, scattering time due to disorder, and mobility.Using extracted parameters of the heterostructure, we identified the range of magnetic fields, where the gapless regime of superconductivity is realized and generated predictions for the density of states that may be tested in future experiments.
Although our model was able to describe the experimental data and generate predictions, a number of questions remains open.First, while our theoretical model incorporates disorder compared to earlier theoretical studies, it still relies on a number of approximations.In particular, it may be desirable to relax the assumption of the perfectly transparent interface between 2DEG and superconductor, incorporate orbital effect of the magnetic field and inverse proximity effect -phenomena related to more realistic description of the motion of electrons between 2DEG and superconductor. 1Incorporation of orbital effects may be particularly important for Germanium that has small value of g-factor and also for surface states of topological insulator proximitized by the bulk superconductor [25] where orbital contribution dominates the physics.In addition, to make the model of 2DEG more realistic, one may incorporate the Dresselhaus spin-orbit coupling along with Rashba spin-orbit considered here, and take into account more realistic band structure of 2DEG.These ingredients may be particularly important for treatment of heterostructures beyond Al-InAs, such as Al-Ge [39], Al-InSbAs [65], and other materials.
Bringing additional ingredients into our theoretical model, that already relied on a sixparameter fitting for describing experimental data, would most likely require additional experimental probes.In particular, it is easy to relax the assumption that the induced order parameter in the 2DEG is equal to its value in the superconductor, and relying on the tunneling measurements use more accurate relation.Using the framework for Green's function calculation set up in this work, one can easily calculate other experimental observables such as spin susceptibility and finite-frequency electromagnetic response kernel.Including more observables measured in situ, such as spin susceptibility [66], optical conductivity [48], density of states, or noise spectra [67] would enable even more detailed material characterization and allow for independent verification of our model.
Another set of questions that remain open is the eventual fate of the gapless proximityinduced superconducting state in 2DEG at sufficiently weak disorder and in the clean limit [1].While the presence of superconducting film suggests that this instability will not destroy pairing in the system, it may cause reconstruction of the low energy band structure leading to reduced density of states.Understanding the instability requires incorporation of the inverse proximity effect and self-consistent treatment of the superconductor that takes into account the 2DEG, that are beyond the present work.Construction of such model may assist the future experiments in the identification and characterization of the phase diagram.
Finally, the experimental control available in the microwave cavity setups akin to Refs.
[1] call for the theoretical development of nonequilibrium probes of 2DEG-superconductor heterostructures.In particular, the cavity enables to excite the system in the regime beyond liner response, and potentially probe it in the time-resolved fashion.In addition, electric contacts may enable to use current as an alternative pair-breaking mechanism to the in-plane magnetic field.These capabilities invite the theory of such 2DEG-superconducting system subject to nonequilibrium effects and currents, that would allow to gain further insights into the material properties and potentially uncover new phases.
The DOS is even function, ρ(ω) = ρ(−ω) so we defined it only for positive values of energy.For V > ∆ we have, In these equations K(n) is the complete elliptic integral of the first kind and Π(n, m) is the complete elliptic integral of the third kind.These results for DOS are presented in the Fig. 1.We note, that these expressions have logarithmic singularities at energies ω = ±(∆ + V ) and discontinuous jumps at ω = ±(∆ − V ), that replace the conventional square-root BCS divergence in the DOS at ω = ±∆.

C Derivation of a response kernel Q αβ
In this appendix we present the derivation of the electromagnetic response kernel Q αβ .To construct the kernel, we use the explicit expression for current operator that is given by the sum of three terms when written in the real space, ĵ = ĵ (grad) + ĵ(so) + ĵ(dia) , (C.1) Figure 8: Dependence of the normalized superfluid density n s,x x /n on the angle between the magnetic field and x-axis for different forms of g-tensor.Dashed blue line in both panels serves as a reference and corresponds to isotropic and diagonal g-tensor with g x x = g y y = 10 and g x y = g y x = 0. Panel (a) shows the effect of the anisotropy in the diagonal terms of g-tensor, when off-diagonal terms are vanishing.Panel (b) demonstrates that the off-diagonal terms in the g-tensor tilt principal symmetry axes.Radial divisions start at n s,x x /n = 0 and are in increments of 0.2; µ B H = 0.13∆ 00 and τ∆ = 0.5.
that represent the standard gradient contribution, spin-orbit contribution, and, finally the diamagnetic term, Separating the diamagnetic contribution, that results in the term proportional to the total carrier density, n, the kernel Q αβ (q, iω n ) can be expressed as follows where we indicated the coordinate and time dependence of operators.Substituting explicit expressions for electric current, rewriting the correlator using Green's function, and expressing the diamagnetic contribution as the product of Green's functions with ∆ set to zero we obtain Eq. ( 13) in the main text.We note, that there is no vertex correction from disorder due to the following facts: 1) scattering potential is local; 2) the Green function G(p) is an even function of p.As a result, the vertex correction vanishes because it is proportional to the following integral: p p α G(p)G(p), which is equal to 0 since the integrand is an odd function of p.More details on how vertex correction is calculated can be found in Refs.[44,45], where disorder potential is not assumed to be local and as a result, the vertex correction is non-trivial.

D Anisotropic effects D.1 Superfluid density in the presence of magnetic anisotropy
In this Appendix we examine the different effects of magnetic anisotropy, encoded into gtensor on the superfluid density, n s,x x (in the absence of order parameter suppression).In particular, we show the effect of the off-diagonal g-tensor and discuss potential for observing its effect experimentally.We use Eq. ( 21) to plot the dependence of the superfluid density, n s,x x , on the angle between a magnetic field and x-axis for different forms of g-tensor.First, we consider the case when off-diagonal elements of g-tensor are absent, g x y = g y x = 0 and diagonal terms are different, see Fig. 8(a).Note, that even in this case the superfluid density depends on the direction of the magnetic field -a property that stems from the spin-momentum locking.Upon including the anisotropy g x x ̸ = g y y , the dependence of n s,x x on the angle of the field becomes either more symmetric when g x x > g y y .In the opposite case, g x x < g y y , the dependence of n s,x x on the field angle develops more pronounced asymmetry and assumes an hourglass shape, with the smallest superfluid density corresponding to the case when magnetic field is parallel to the y-axis.Next, we explore the case when off-diagonal elements of g-tensor are present and diagonal are equal, see Fig. 8(b).This type of anisotropy gives the dependence also an hourglass shape, and tilts its principal symmetry axes away from x and y directions.
The qualitative difference in the shape of n s,x x as a function of the angle may be used to detect and quantify the g-tensor anisotropy in the experiment.At the same time, addition of the off-diagonal elements in the g-tensor as fitting parameters may potentially result in the overfitting, especially if other six fitting parameters used in the main text are still present.Hence, in order to unambiguously determine the form of the g-tensor using this framework, it is desirable to fix at least some of the other material properties using different or additional experimental data, or directly measure off-diagonal components of the superfluid density as discussed in the main text.

D.2 Generalization of spin-orbit coupling
In the main text we assumed that spin-orbit coupling is of a Rashba type.In this Appendix we show that our results are applicable not only to Rashba spin-orbit coupled system and can be generalized by a simple redefinition of the g-tensor.In particular, they can be also applied to the system that has only Dresselhaus spin-orbit coupling.
Rashba spin-orbit coupling has the following contribution to the Hamiltonian in momentum representation, Instead of 2 × 2 antisymmetric tensor above, we consider more general case -an arbitrary orthogonal matrix U. Any orthogonal matrix U can be expressed as: where U rot = cos θ − sin θ sin θ cos θ is an arbitrary rotation matrix and d = ±1 corresponds to reflections.Then, we make the changes of variables, that corresponds to rotation in the spin space and potentially reflection in the momentum space.After this transformation, the Hamiltonian takes the same form as for the case of Rashba spin-orbit coupling in Eq. (D.1).
Let us now analyze how this change of variables affects other terms in the Hamiltonian.To begin with, we consider the Zeeman contribution, here, ĝ = U rot ĝ.It is clear that the reflection in momentum space does not affect the kinetic term H k because it is quadratic in momentum, the same applies to the calculation of diagonal components of superfluid density tensor; off-diagonal components are multiplied by factor d since they are proportional to k x k y .Thus, we have reduced the problem to the original one and all results for Rashba spin-orbit coupling are applicable to this generalization of spin-orbit coupling, the only change is the redefinition of g-tensor: ĝ → U rot ĝ.In particular, if we put θ = π/2 in U rot and d = −1 in reflection matrix, we get the Dresselhaus spin-orbit coupling.

E An effect of out-of-plane component of magnetic field on pairbreaking parameter
In Sec.4.1 we studied pair-breaking in Al and found a linear contribution from the magnetic field to the pair-breaking parameter α.Assuming that this linear contribution appears because of a small out-of-plane component of a magnetic field in the experimental measurements, let us estimate this component.The pair-breaking parameter for the case when a magnetic field is perpendicular to the layer is expressed as follows: α ⊥ (H) = DeH [58,59], while for the magnetic field parallel to the layer, it reads: α ∥ (H) = D (eH d) 2 /6, where D is the diffusion constant, d is the thickness of the superconducting layer, e is the charge of the electron, and H is a magnetic field.Assuming that pair-breaking is weak, α/∆ 00 ≪ 1, we can make the following approximation for resulting pair-braking parameter: α(H) ≈ α ∥ (H ∥ ) + α ⊥ (H ⊥ ) where H ⊥ = H sin φ and H || = H cos φ are out-of-plane and in-plane components of magnetic field and φ is an angle between a magnetic field and the plane.Then, taking into account Eq. (31)  Taking the ratio of these equations we obtain For the calculation of superfluid density and order parameter suppression, the following temperature and critical magnetic field are used: T /∆ 00 = 0.0564, V cr = 2.5∆ 00 .Dashed lines correspond to the areas of magnetic field where current theory is not applicable (superfluid density becomes negative).Purple areas represent the regions of magnetic field where there is a gapless regime.Vertical lines in the plots from the upper panel show the magnetic field at which the DOS is plotted in the lower panel; the green line represents the magnetic field at which the system is in the regime with a gap, while the blue one represents the gapless regime.
We should notice that this estimate is applicable only when pair-breaking is weak which corresponds to the small magnetic fields: H/H cr ≪ 1. Experimentally Ref. [1] performed very careful field alignment with the out of plane component being constrained to less than 0.1% at considerable values of the magnetic field, that is an order of magnitude smaller compared to our estimate above.This discrepancy may be attributed to the fact that our estimate is applicable only at very small fields, also the systematic contributions that are beyond our theoretical model, such as field-dependence of inverse proximity effect, may play a role here.

F Renormalization of the order parameter and magnetic field by disorder
In this Appendix, we present the dependencies of real and imaginary parts of parameters ∆ (1) , ∆ (2) , and V (1) on energy.The data is shown for three values of disorder strength: τ∆ = 0.1, τ∆ = 0.5, and τ∆ = 5, with each plot comparing the renormalized functions for two values of magnetic field that put the system into regime with the gap, and gapless regime.Figure 9 demonstrates that for energies far away from zero, real part of gaps, ∆ (1,2) and renormalized magnetic field V (1) tend to their true values, and imaginary parts approach zero.At the same time, for energies of order V and ∆, renormalized functions strongly depend on the energy.In particular, all renormalized functions have features at the energies of the order ±|V ± ∆|.For small disorder (large τ) these features are sharp, and they smoothen out with increasing disorder strength.
Another notable feature is the symmetry properties of the real and imaginary part of the renormalized functions.First, for small values of V when DOS has the gap, the Im ∆ (1) (ε) = 0 for ε that is sufficiently small.The same property also holds for imaginary parts of the Im ∆ (2) (ε) and Im V (1) (ε).Second, the renormalized s-wave gap component, ∆ (1) (ε) has even in energy real and odd in energy imaginary parts.
In contrast, the induced triplet component has odd in energy real and even in energy imaginary parts, highlighting that the presence of disorder induces odd frequency triplet component of the order parameter [69].The magnitude of the triplet component is increasing with disorder strength, and for large disorder it becomes in some energy intervals even larger that the order parameter in the superconductor.At the same time, in this regime both real and imaginary contributions to ∆ (2) (ε) are of comparable magnitude, intuitively suggesting that disorder-induced triplet pairs are short lived.

G Superfluid density and DOS in the presence of order parameter suppression
In this Appendix, we present the plots for superfluid density and DOS in the presence of order parameter suppression.Taking values of the temperature T = 0.0564∆ 00 and critical field V c = 2.5∆ 00 we show resulting superfluid density and DOS in Fig. 10.Data in this figure is qualitatively similar to the case where the order parameter suppression is not taken into account, see Fig. 2 in the main text.The only differences are that the superfluid density vanishes at the critical field of superconductor, as expected and the boundaries of the gapless regimes move to smaller values of magnetic field, which can also be observed in the plots in Fig. 3.

3 15 4 3 20 5and Outlook 22 A DOS in the clean case 24 B Approximations in Green function 25 C 27 E
Results for proximity-induced superconductivity in disordered 2DEG 13 3.1 Numerical solution for Green's function 13 3.2 DOS and superfluid density 13 3.3 Region of the gapless superconductivity and potential instability Material parameters and prediction for DOS and gapless regime Summary Derivation of a response kernel Q αβ 25 D Anisotropic effects 27 D.1 Superfluid density in the presence of magnetic anisotropy 27 D.2 Generalization of spin-orbit coupling An effect of out-of-plane component of magnetic field on pair-breaking parameter 28 1 Introduction

Figure 1 :
Figure 1: (a) Sketch showing 2D-heterostructure, where the superconductor (SC) with the s-wave pairing proximitizes the 2DEG.The 2DEG has strong Rashba spinorbit coupling (SOC), and the entire heterostructure is subject to the in-plane magnetic field pointing in y-direction.(b) Rashba-split bands of 2DEG spectrum without magnetic field lead to two spin-momentum locked Fermi surfaces (red and blue) shown on the left.The right part of panel (b) shows the red Fermi surface gapped due to proximity effect without (top right) and with (bottom right) the magnetic field.The presence of a magnetic field tilts the spectrum in an anisotropic manner with gapped bands touching the Fermi level first along the x-direction.(c) Cut of the band structure in presence of induced gap and weak magnetic field V < ∆ along the x-axis and (d) associated DOS that shows the decreased gap, the discontinuous jump in DOS at the gap edge, followed by the logarithmic Van Hove singularity.Panels (e)-(f) show similar data for larger value of magnetic field V > ∆, when Bogoliubov's Fermi surfaces form and DOS becomes gapless.Panels (c) and (e) use chemical potential µ = 20∆ and spin orbit coupling λ so k F = 5∆ in units of proximity-induced gap, ∆.

Figure 3 :
Figure3: The span of magnetic fields where gapless superconductivity is realized (filled region between solid and dashed violet lines) expands with increasing disorder strength.The solid violet line corresponds to the onset of the gapless regime, dashed violet line marks the point when superfluid density turns negative, signaling that our model becomes inapplicable.Data shown in panel (a) assumes that the order parameter is not suppressed by magnetic field, panel (b) takes the depairing into account assuming critical magnetic field: V c = 2.5∆ 00 , orange dotted line corresponds to this critical field V c .

Figure 4 :
Figure 4: Fitting the theoretical model (solid line) to the experimental dependence (dots) of critical temperature on a magnetic field gives an excellent agreement for the value of fitting parameter ζ 1 = 0.207.The dashed line corresponds to the fit in the absence of linear contribution from magnetic field to pair-braking parameter (ζ 1 = 0).

Figure 6 :
Figure 6: Polar plot for the dependence of resonance frequency on the angle between a magnetic field and x-axis, radial divisions start at f r = 4.8 GHz and are in 0.2 GHz increments.Different colors encode the values of the magnetic field, with dots corresponding to experimental data and solid lines showing the best fit of the theoretical model.

4 )Figure 10 :
Figure10: Results for superfluid density (top) and DOS (bottom) in the presence of order parameter suppression.For the calculation of superfluid density and order parameter suppression, the following temperature and critical magnetic field are used: T /∆ 00 = 0.0564, V cr = 2.5∆ 00 .Dashed lines correspond to the areas of magnetic field where current theory is not applicable (superfluid density becomes negative).Purple areas represent the regions of magnetic field where there is a gapless regime.Vertical lines in the plots from the upper panel show the magnetic field at which the DOS is plotted in the lower panel; the green line represents the magnetic field at which the system is in the regime with a gap, while the blue one represents the gapless regime.

Table 1 :
Result of the fitting procedure.Magnetic anisotropy is significant since g eff x x differs from g eff y y .The disorder is intermediate meaning that dirty-limit and cleanlimit approximations are not applicable to the specific Al-InAs heterostructure.
we obtain that Matching linear and quadratic in H terms with each other separately we get,