Spin-Polarized Tunneling in Critically Disordered Be-Al Bilayers

We report spin-polarized tunneling density of states measurements of the proximity modulated superconductor-insulator transition in ultra thin Be-Al bilayers. The bilayer samples consisted of a Be film of varying thickness, $d_\mathrm{Be}=\,$0.8-4.5 nm, on which a 1 nm thick capping layer of Al was deposited. Detailed measurements of the Zeeman splitting of the BCS coherence peaks in samples with sheet resistances $R\sim h/4e^2$ revealed a super-linear Zeeman shift near the critical field. Our data suggests that critically disordered samples have a broad distribution of gap energies and that only the higher portion of the distribution survives as the Zeeman critical field is approached. This produces a counter-intuitive field dependence in which the gap apparently increases with increasing parallel field.


I. INTRODUCTION
The disorder driven superconductor-insulator transition (SIT) has been the subject of intense investigation for more than 30 years now [1][2][3][4][5][6]. Early studies suggested a relatively simple picture of the SIT in homogeneously disordered two-dimensional (2D) systems. As the disorder of the superconductor is increased, the underlying repulsive Coulomb correlations are enhanced until they eventually overwhelm the resident superconducting correlations. This was believed to occur at a relatively welldefined critical disorder characterized by the quantum resistance R Q = h/4e 2 [1,7]. However, more recent studies of the SIT have shown that the disorder-driven transition can be more complex than originally thought. Under the proper conditions, local variations in film disorder are amplified by Coulomb interactions thereby producing regions of varying gap strength. This leads to an intrinsic puddling of the superconducting condensate as the disorder is increased thought the SIT [8][9][10][11][12]. In this scenario the insulating side of the SIT is fundamentally bosonic in that it consists of phase-decoupled superconducting puddles each having a finite order parameter [9]. Recent studies utilizing multiply connected geometries [13,14] show that superconducting pair correlations exist well into the insulating phase of highly disordered, nominally homogeneous, Bi films [6,9,10]. This phase is commonly referred to as a Bosonic insulator. But the details of the topological structure of the phase and its corresponding gap distribution remains unclear. Here we present a spin-polarized tunneling study of the evolution of the superconducting gap in critically disordered Be-Al bilayers that are tuned through the SIT into a Bose insulator phase via a Zeeman field. Parallel magnetic field is used to induce a Zeeman splitting of the BCS density of states (DoS) spectrum while minimizing the orbital field broadening of the spectral features. As the critical field is approached from below, we find that the supragap spin band shifts super-linearly with increasing field. This suggests that at low Zeeman fields the tunneling conductance measures the average of a rather broad distribution of gap energies [15]. But as the critical field is approached from below, only the highest energy portion of the distribution remains intact thereby producing an apparent increase in average gap energy with increasing field. The resulting spectra represent a direct probe of the pairing amplitude distribution across the SIT.
Our primary goal in this study was to use the Zeeman splitting of the BCS DoS spectrum to probe the superconducting phase of quasi-homogeneously disordered films which are close to the disorder-driven zero field SIT. In order to resolve the Zeeman splitting, the superconductor must have a low spin-orbit (SO) scattering rate. We also require that the transition temperature of the film be well above the base temperature of our fridge. Thin Al films meet these conditions; they have a T c ∼ 2.7 K and an extremely low intrinsic SO rate. Unfortunately, however, Al forms granular films, even when they are quench-condensed at 84 K. The granularity makes high resistance Al films somewhat unstable in air. Furthermore, the transport properties of high resistance Al films can be completely dominated by the intragrain coupling and not necessarily by many-body effects. In order to circumvent this limitation, we deposit the Al on a thin layer of Be. Beryllium forms dense, adherent, amorphous films on glass substrates and, like Al, it has a low SO rate [3]. However, the transition temperature of Be films, T c ∼ 0.5 K, is too low for our purposes [16]. The bilayer arrangement allows us to partially mitigate the granularity of the Al by providing an underlying metallic coupling between the Al grains while still maintaining a reasonably high T c by virtue of the proximity effect [17].

II. EXPERIMENTAL METHODS
Superconducting Be-Al bilayer films were formed by first depositing a thin Be layer of varying thickness onto arXiv:2010.08510v2 [cond-mat.supr-con] 12 May 2021 a fire-polished glass substrate followed by the deposition of a 1 nm-thick Al film. The depositions were performed by e-beam evaporation from 99.9% Be and 99.999% Al targets at a rate of ∼ 0.2 nm/s. The glass substrates were maintained at 84 K during the deposition of both the Be and Al layers. The bilayers were deposited without breaking a vacuum of P < 3 × 10 −7 Torr. Planar tunnel junctions were formed between the upper Al layer of the samples and a counter-electrode composed of a nonsuperconducting Al alloy using a 1 nm layer of SiO as the tunnel barrier. Bilayers with an Al thickness of 1 nm and Be thicknesses ranging from 0.8 to 4.5 nm had normal state sheet resistances that ranged from R ≈ 100 Ω to 10 4 Ω at low temperature. Magnetotransport measurements were made on a Quantum Design Physical Properties Measurement System He3 probe. The maximum applied field was 9 T and the base temperature of the system was 400 mK. The tunneling measurements were carried out using a standard 27 Hz 4-wire lock-in amplifier technique.

III. RESULTS AND DISCUSSION
In general the critical field of a thin film superconductor has both an orbital and a Zeeman component [18]. If one makes a low atomic mass film [19], such as Al, sufficiently thin and orients the field parallel to the film surface then the orbital response will be suppressed and one can realize a purely Zeeman-mediated critical field transition [16,20], which is first order at low temperatures. The T = 0 parallel critical field is given by the Clogston-Chandrasekhar equation [21], H c = √ 2∆0 gµB , where ∆ 0 is the zero temperature -zero field gap energy, µ B is the Bohr magneton, and g is the Landé g-factor. In this series of experiments we have explored the Zeeman response of Be/Al bilayers with sheet resistances ranging from R R Q to R ∼ R Q . Shown in Fig. 1 are the bilayer transition temperatures T c as a function of the Be thickness. For Be thicknesses greater than d Be ∼ 1.5 nm the bilayer transition temperatures lie between that of pure Be film, T Be c ∼ 0.5 K, and that of a pure 2.5 nm Al film, ∼ 2.7 K. The fact the T c decreases with increasing d Be in this range suggests that the bilayer transition temperature was mediated by the proximity effect [17]. In the thickness range d Be < 1.5 nm, there is a local maximum in T c . In the limit d Be → 0, T c should be asymptotic to the transition temperature of the underlying Al layer. Unfortunately, the T c of the 1 nm Al layer is unknown because the bilayers become electrically discontinuous for d Be 0.5 nm. Nevertheless, it is possible that the Al layer has a transition temperature well above that of a 2.5 nm Al film. In fact, the data in Fig. 1 suggests T Al c ∼ 3.5 K. Bilayers with d Be 1 nm have sheet resistances approaching the quantum resistance R Q = h/4e 2 [7]. Since R Q represents the threshold for the SIT we believe that the local maximum in Fig. 1 arises from the preemptive effects of increasing disorder. In Fig. 2 we plot the temperature dependence of the sheet resistance of a d Be = 0.8 nm bilayer in a range of applied parallel magnetic fields. The "fan-like" structure of the data is associated with the field-tuned SIT and has been reported in a wide variety of critically disordered thin film superconductors [1,2,5,11]. Note the separatrix between the superconducting and insulating phases at a transport parallel critical field H tr c ∼ 5 T. These data demonstrate that our thinnest Be layer samples are near the zero field SIT and that their underlying normal state is insulating [22,23]. Furthermore, there is an obvious low temperature crossing of the 7 T and 9 T traces. This suggests that there is a possible superconducting contribution to the insulating behavior over a narrow range of fields. In other words, the data in Fig. 2 is consistent with the emergence of a Bose insulator phase at intermediate parallel fields. Below we present a tunneling density of states study of the Zeeman-tuned SIT via the application of parallel magnetic [24]. Shown in Fig. 3 are tunneling conductance spectra taken on the bilayer of Fig. 2. At low temperatures the tunneling conductance G is simply proportional to the single-particle density of states (DoS) [25]. The bias voltage is relative to the Fermi energy and the conductances have been normalized by the conductance at 2 mV. The solid trace in the upper panel represents the tunneling conductance in the superconducting phase of the bilayer and the dashed trace is the corresponding conductance in the high-field normal phase. We note that based on the DoS measuerements, the parallel critical field, H c ∼ 7 T, is higher than that estimated from transport measurements; we will return on this point later. The are three features of the tunneling spectra in Fig. 3 that are of particular importance to this study. The first is the Zeeman splitting of the BCS coherence peaks [18,26,27]. The second is a broad logarithmic suppression of the DoS in the normal phase of the bilayer. This suppression arises from e − e interaction effects [23] and is often referred to as the zero bias anomaly (ZBA) or the Coulomb anomaly [3]. It is a direct microscopic measure of the disorder-induced repulsive correlations and has been welldocumented in a wide variety of 2D systems. The third is the pairing resonance (PR) represented by the dips riding on top of the normal state spectra [28,29], as indicated by the arrows in Fig. 3. Details of the PR have been published elsewhere [30], but for our purposes the resonance, which arises from evanescent Cooper pairs, provides us with a direct probe of the spin properties of the normal state (see Appendix A). In particular, we can use the PR to extract the effective normal state g-factor, which may differ from the naive value g = 2 due to Fermi liquid (FL) effects [31].
As the beryllium thickness is decreased below 1 nm both the resistance of the bilayers and the magnitude of the ZBA increase precipitously. This is due, in part, to the fact that a 1 nm Al film deposited directly on the glass substrate would not be electrically continuous. Therefore, for our chosen geometry the Be thickness controls the level of disorder. In order to quantify the strength of the ZBA we define the dimensionless parameter α = [G(2mV ) − G(0)]/G(2mV ) which measures the relative depletion of electron states at the Fermi energy in the high-field normal phase. A plot of α as function of beryllium thickness is shown in the inset of Fig. 2. Note that α grows rapidly as d Be is lowered below 2 nm, indicative of the approach to the zero-field SIT. As is evident in Fig. 3 the depletion of single particle states due to the ZBA and the depletion due to opening of the superconducting gap are comparable in magnitude in our most disordered samples.
Shown in the lower panel of Fig. 3 is the superconducting spectrum of the upper panel after the ZBA spectrum has been divided out of the data. In order to sup- The value of∆ has been chosen to match the extrapolated zero-field peak position of the experimental data; this also fixes the field scale as∆/µB. press the PR the normal state ZBA spectrum was measured by applying a 9 T perpendicular field. The four peaks, the two outer supra-gap peaks and the two inner sub-gap peaks, represent the Zeeman splitting of the BCS spin-up and spin-down subbands [18,26,27]. The occupied and unoccupied subband peaks are located at eV = ±(∆ 0 ± 1 2 eV Z ) where eV Z = gµ B H is the Zeeman energy. Although these data were taken relatively close to the parallel critical field H c ∼ 7 T, the peaks are sharp and their positions can, in principle, be measured quite accurately through the parallel critical field transition. In Fig. 4 we plot the position of the supra-gap and sub-gap coherence peaks as a function of parallel field for a d Be = 0.8 nm bilayer having an α ∼ 0.4. However, near H c the spectrum begins to display characteristics of both the superconducting and normal phases. In particular, the PR dip emerges and is superimposed on the superconducting spectrum. The position of the PR, is quite close to, and partially overlaps with, the superconducting sub-gap peaks near H c . This makes it difficult to determine the sup-gap positions in the critical region. The supra-gap peaks, however, are positioned well away from the PR, so we have chosen to focus our analysis on their field dependence. The supra-gap peak positions from the bilayer of Fig. 4 are shown in Fig. 5 along with the corresponding peaks in a moderately disordered Al film, R ∼ 1000 Ω. The solid line in this plot represents V Z with g = 2. At fields well below H c the supra-gap peaks in both films exhibit the expected Zeeman shift. However, as the critical field is approached the Al data falls below the Zeeman line while the Be-Al data rises super-linearly. The sub-linear field dependence of the Al peak is due to both a downward FL renormalization of the quasiparticle spin [31] and to small, but finite, orbital pair-breaking effects. The former can be estimated by measuring the PR position as a function of field, see inset of Fig. 5. The dashed line in the inset represents a least-squares fit of Eq. (1) to the data in which g and ∆ 0 where varied. The fit clearly shows that the normal state g-factor in the bilayers is also well below 2. Therefore, the super-linear field dependence of V pk cannot be due to FL exchange effects.
The super-linear field dependence of V pk was only observed in highly disordered bilayers having α 0.2 (see Appendix B), which suggests that the effect is a consequence of disorder-enhanced e − e interactions as manifest in the ZBA. In the strong disorder limit the system begins to break up into weakly connected superconducting islands which have a broad range of local gap energies [15,32,33]. Thus, the super-linear field dependence of V pk reflects the portion of the gap distribution which survives as H → H c . At low fields the tunneling conductance samples the entire distribution of gaps but near the critical field the sampling is skewed towards high gap values. We have modeled this behavior by averaging over a Gaussian distribution of zero-field gaps (see Appendix C). This approach is oversimplified at low field, where different parts of the sample are all superconducting and can influence each other directly, but should be qualitatively correct at high fields, where the surviving superconducting regions are disconnected [32] (a fully self-consistent calculation as those presented in Refs. [15,32] is beyond our scope). The result of such a calculation displays a supra-linear behavior resembling the experimental data, see the inset in Fig. 4. The existence of isolated superconducting regions between 5 and 7 T could also explain the discrepancy between transport and DoS estimates for the parallel critical field. Interestingly, fitting the DoS with a phenomenological Dynes broadening, while inadequate, also supports our interpretation, see Appendix D.

IV. SUMMARY
In summary, we have exploited the Zeeman splitting of the BCS DoS spectrum to obtain direct evidence for a distribution of pairing energies in a critically disordered BCS superconductor. The effects of disorder serve to both suppress T c and broaden the critical field transition. Our data suggest that the Zeeman-tuned SIT is dominated by the tail of a broad distribution of local gap energies. Presumably this high energy tail also plays a prominent role in the zero-field, disorder-driven SIT transition.

ACKNOWLEDGMENTS
The magneto-transport measurements were performed by P.W.A. and F.N.W with the support of the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award No. DE-FG02-07ER46420. The theoretical analysis was carried out by G.C.

Appendix A: Pairing resonance
A low spin-orbit BCS superconductor in the presence of a pure Zeeman field will undergo a low-temperature first order transition to the normal state at a critical field that is simply proportional to the gap energy [18], For fields above H c the system enters a paramagnetic phase and the mean-field BCS order parameter is quenched. Although there is no global superconducting order parameter in the paramagnetic phase, a wellpronounced superconducting fluctuation mode still exists [34]. Electrons that tunnel with the proper energy can, in fact, produce a resonant excitation of this mode and thereby cause a strong tunneling DOS singularity which we term the pairing resonance (PR). The bias voltage V * corresponding to the resonant mode and the corresponding DOS singularity is a function of the supercritical field [34], see Eq.
(1). The mode persists to fields well beyond H c . As the field increases the PR broadens and moves out to higher energy as per Eq. (1). Note that the implicit g-factor g dependence in Eq. (1) can be used to determine g in the paramagnetic phase by measuring V * (H ) and the then fitting the data to Eq. (1). The PR is extremely sensitive to both spin-orbit scattering and perpendicular field. Both of these strongly suppress the resonance. Shown in the upper panel of Fig. 6 is the normal state tunneling spectrum for the d Be = 1.0 nm bilayer in a 9 T parallel field. The satellite dips are the PR which are superimposed on the e−e interaction mediated ZBA [23]. The lower panel is the resulting spectrum when a 4 T perpendicular field is applied. Note that PR is completely suppressed but the background ZBA remains largely unaffected.
Appendix B: Parallel-Field SIT in a 1 nm -1 nm Be-Al bilayer

Transport
In Fig. 7 we plot the temperature dependence of the sheet resistance of a bilayer having layer thicknesses d Be = 1.0 nm and d Al = 1.0 nm in a range of applied parallel magnetic fields. The "fan-like" structure of the data is very similar to that of the higher resistance d Be = 0.8 nm bilayer shown in Fig. 2. The transport data suggests a that the T = 0 SIT critical field for this sample is H c = 7 T. Note that the 7 T traces separates the superconducting-like and insulating-like temperature dependencies in Fig. 7. Indeed the 7 T trace suggests that an intermediate metallic phase exists at the boundary of the SIT. The insulating behavior of this bilayer is not as pronounced as that of the one in Fig. 2, but this is expected given that the normal state resistance of the former is less than half that of the latter. Interestingly, the crossing of the 8 T and 9 T traces in Fig. 7 is similar to what is seen in Fig. 2.

Tunneling
Here we present tunneling density of states measurements across the SIT shown in Fig. 7. Tunneling allows us to establish the relationship between the measured gap and the transport characteristics at a given parallel field. In panel (A) of Fig. 8 we plot the normalized tunnel conductance as a function of bias resistance at H = 4.5 T. As can be seen in Fig. 7 the system is well in the superconducting phase at this field which is reflected in the corresponding tunneling spectrum. The classic BCS coherence peaks are clearly evident but have been Zeemansplit into spin-up and spin-down sub-bands by the applied field. Interestingly, in panel (B) we see a similar superconducting spectrum at a substantially higher field of H = 7.0 T although the corresponding 7.0 T transport trace in Fig. 7 suggests that the system is in a metallic phase. At 7.5 T the transport clearly exhibits insulating behavior but the corresponding tunneling spectra in panel (C) of Fig. 8 indicates that a finite superconducting condensate persists. Although the coherence peaks have been suppressed in this spectrum, a sharp break in the background (indicated by the arrows) marks superconducting contribution to the tunneling conductance. This incommensurability between transport and tunneling is similar to what we observed in the d Be = 0.8 nm bilayer in the main text. Specifically, these data suggest that local superconductivity plays a role the insulating behavior observed in the transport. Finally in panel (D) we show what we believe to be a normal state spectrum at 9.0 T. All superconducting signatures are absent save for the pairing resonance, which appears as satellite dips su- perimposed on the zero bias anomaly (ZBA) background [23]. Note that the corresponding 9.0 T transport trace in Fig. 7 appears to be less insulating than that of the 8 T suggesting that there may be a bosonic contribution to the insulating behavior in the field range of 7 to 9 T.

Supra-gap peak position
In Fig. 9 we plot the supra-gap peak position as a function of parallel field as obtained from the spectra represented in Fig. 8. This plot corresponds to Fig. 5. The solid line represents the linear Zeeman dependence assuming g = 2. Typically the supra-gap peak positions in a relatively low disorder film will fall sightly below the Zeeman line, see Fig. 5. This is primarily due to the fact that the parallel field causes a small suppression of the mean-field gap via orbital pair breaking. Although the bilayer thickness is much less than the coherence length, the orbital suppression of the gap cannot be completely attenuated. Furthermore, the g-factor in the normal state is somewhat less than 2 due to Fermi liquid exchange effects and as the critical field is approached the quasiparticle exchange turns on thereby effectively decreasing g : 2 → 1.6 [30]. Arrow 1 in Fig. 9 marks the initial break from the Zeeman line that is later followed by the emergence of super-linear field dependence in the region of arrow 2. We believe that the this super-linear field dependence arises from a distribution of gap energies as was the case for Fig. 5. However, since this bilayer is substantially less disordered the effect is somewhat less pronounced but nevertheless clearly evident.
Appendix C: Gaussian model of independent superconducting gaps Experimental evidence for a roughly Gaussian distribution of local superconducting gaps has been collected over the years in disparate superconducting material such as high-T c superconductors [35], iron-based superconductors [36], and low-T c materials [37]. The gap distribution standard deviation can be a significant fraction of the average gap: 20 %, 14 %, and 6 %, respectively, in the cited works. Since the critical fields increase with the gap, the local gap variation can lead to a local variation of the critical fields. In an extremely simplified model, we treat different parts of a sample as separate and independent superconducting puddles. We do not expect such model to be realistic at low field, where all the puddles are superconducting and therefore can influence each other via the proximity effect; however, if at high field a sufficient portion of the puddles have turned normal, such a treatment should be qualitatively correct.
In our numerical calculations, we proceed as follows: we start with a Gaussian distribution of gaps with the standard deviation σ being 10 % of the average zero field, zero temperature gap value∆ 0 and we truncate the distribution to ±2σ. We then divide the distribution into 39 bins and normalize the weights so that the sum of all weights is 1. From the zero-field, zero-temperature gap value ∆ 0 of each bin, we self-consistently calculate the order parameter as function of field, taking into account the orbital effect of the parallel field H as well as the temperature (0.4 K) at which the measurement is done. The approach used here is a simplified version of that in Ref. [31], since we are neglecting the possible effects of spin-orbit scattering and Fermi-liquid renormalization of the spin susceptibility. For the orbital parameter c = 0.02, we have used the value extrapolated from measurements in thicker Al films.
Once the dependence of the order parameter on parallel field is known, we can calculate the energy-and field-dependent superconducting density of state (DoS) ν ∆0 ( , H) for each gap value ∆ 0 ; finally, for a given field H we average over the Gaussian distribution by performing a weighted sum of the densities of states. If H exceed the maximum possible field [38] for a given ∆ 0 , the corresponding DoS is taken as normal-metal (i.e., energy independent) one. As an example, we show in Fig. 10 a plot of the average DoS calculated at the maximum field for the average gap (solid line) and for comparison the DoS calculated at the same field for a single value of the gap corresponding to the average gap (dashed line). The finite subgap DoS for the solid curve is due to the normalstate puddles (note that we have not included the effect of the zero-bias anomaly in our model). The average DoS clearly displays broader features than the non-averaged one, and the outer peak positions are at higher energy. Since its proposal by Dynes and collaborator in 1978 [39], fitting of (normalized) tunneling density of state data is often performed by using the following expression here generalized to allow for Zeeman splitting in the presence of parallel magnetic field, E Z = µ B H. In this formula, γ is the so-called Dynes broadening; in normal/superconducting bilayers with a tunnel contact (weak proximity effect), the Dynes formula correctly accounts for the finite time electrons spend in the superconducting layer before tunneling into the normal metal [40]. In contrast, for a good contact between two superconducting materials (strong proximity effect), a hard gap with a square root threshold replacing the square root singularity is theoretically expected. Nonetheless, we will use this formula to analyze the experimental data; assuming γ small compared to ∆, it predicts that the maxima of the outer peaks in density of states are located at eV pk ± ∆ + E Z + γ/ √ 3 (D2) In the above equations, ∆ is the self-consistently determined order parameter, which is suppressed in comparison with the order parameter ∆ 0 that one would obtain in the absence of broadening: To account phenomenologically for the orbital pairbreaking effect of the parallel magnetic field, we assume for γ a field dependence in the form Substitution of Eqs. (D3) and (D4) into Eq. (D2) leads to a sublinear evolution of the peak position with parallel field.
We have used Eq. (D1) to fit the outer peaks of tunneling density of states data for a d Be = 0.8 nm bilayer measured at 70 mK in increasing parallel field H up to 5.5 T. We treat γ and ∆ as fit parameters; as shown in Fig. 11, the outer peaks can be well captured by Eq. (D1), and the underestimation at energies further away from the Fermi energy is due to the logarithmic behavior of the zero-bias anomaly (see also Figs. 3 and 6). The density of states at the inner peaks and near the Fermi energy is overestimated, demonstrating that at finite field the simple Dynes phenomenology cannot accurately describe the data.
In Fig. 12 we present the extracted values of the broadening γ as function of field. The dashed line is obtained by fitting Eq. (D4) to the data up to 4.5 T, excluding the two data points at the highest fields. From this fit and the fitted value of ∆ at zero field, we find γ 0 ≈ 0.1 meV, c ≈ 0.19, and ∆ 0 ≈ 0.55 meV. Interestingly, these parameters also satisfactorily describe the behavior of ∆ as function of field, see Fig. 13. As the field increases above 4.5 T, we see that the broadening γ decreases and the gap parameter ∆ increases; this counterintuitive behavior is qualitatively in agreement with what one expects if some regions of the film with weaker superconductivity turn normal, since the surviving regions have on average a larger gap and the gap distribution of the surviving superconducting regions is narrower.