Spin-orbit driven superconducting proximity effects in Pt/Nb thin films

Manipulating the spin state of thin layers of superconducting material is a promising route to generate dissipationless spin currents in spintronic devices. Approaches typically focus on using thin ferromagnetic elements to perturb the spin state of the superconducting condensate to create spin-triplet correlations. We have investigated simple structures that generate spin-triplet correlations without using ferromagnetic elements. Scanning tunneling spectroscopy and muon-spin rotation are used to probe the local electronic and magnetic properties of our hybrid structures, demonstrating a paramagnetic contribution to the magnetization that partially cancels the Meissner screening. This spin-orbit generated magnetization is shown to derive from the spin of the equal-spin pairs rather than from their orbital motion and is an important development in the field of superconducting spintronics.

This supplementary materials contains additional information about 1) the theoretical framework used to describe our data and 2) the experimental techniques used: LEM, STS, TEM and transport. While not essential to the understanding or conclusions of the manuscript, specialist readers may find the additional information useful and informative.

A. Quasiclassical Theory
To compare our experimental results with theory we make use of the quasiclassical framework in the dirty limit (coherence length much longer than the mean free path), where the Green functions obey the Usadel equation [1]. To simulate our scanning tunneling spectroscopy data we require the electronic density of states, which follows directly from the obtained Green functions. For our muon experiments, where we measure local flux, we apply the linear response theory to calculate the screening current in response to a small externally applied field.
We use realistic material parameters for our sample layouts, which where possible we determine directly from independent measurements. The parameters required to describe the interfaces between two materials are tuned such that the predicted superconducting transition temperature matches the experimentally observed value. Furthermore, we initially omit spin-orbit interaction as well as additional pair-breaking mechanisms. We find this describes our system well as long as the Pt layers are very thin (≤ 10 nm), but for our sample using ∼ 50 nm thick Pt layers a clear discrepancy between the flux expulsion observed in our muon experiments and the prediction from theory emerges. For this system with thick Pt, we use a detailed theoretical approach including spin-orbit interactions which is able to explain the muon data. Below follows first the theory modeling omitting spin-orbit interactions, and its application to the STS data (which requires the self-consistent solution for ∆(x)) and the LEM data (which requires the solution for B(x)). The theory including spin-orbit is tailored for sample SM3 and is detailed under the section: Triplet correlations induced by impurity spin-orbit interaction.
We take the x-axis normal to the metallic layers and assume translational invariance in the y,z plane. The Usadel equation for s-wave superconductivity then takes the form i D∂ x (ǧ∂ xǧ ) = [Ȟ ,ǧ ] withǧ the 4×4 matrix Green function in the Nambu ⊗ spin-space, the reduced Planck constant and D the diffusion constant (see e.g. [2]). The Hamiltonian can be described byȞ = i ω n (τ 3 ⊗σ 0 )+∆ with ω n the Matsubara frequencies defined by ω n = πk B T (2n+1) with k B the Boltzmann constant, n integer and the maximum allowed frequency given by the Debye frequency. Furthermore, σ i and τ i are the Pauli matrices of the spin space and Nambu space respectively. The matrix Green function anď ∆ only have non-zero elements on their main and anti-diagonals with diag(ǧ) = ( and anti-diag(∆) = (−∆, ∆, −∆ * , ∆ * ), where G and F are the quasiclassical normal and anomalous Green functions respectively, both being functions of (x, ω n ), ∆(x) is the order parameter and where up/down arrows indicate spin-up/spin-down. The matrix Green function satisfies the normalization conditionǧ 2 =1 and the order parameter must be solved self-consistently satisfying the gap equation: with T c0 the bulk critical temperature. We use the interface boundary conditions as formulated by Nazarov [3], which for the interface between two materials with labels l, r for the layer on the left and right side of the interface respectively can be written as: , with σ i the conductivity of layer i, 0 ≤ Γ ≤ 1 the interface transparency and R b the interface resistance times the interface area.
When the Green functions are known one can calculate the response of the superconductor to a (small) external field. Within the linear response theory, the shielding current density j y (x) in response to the vector potential A y (x) is written as [4]: with N 0 the normal state density of states near the Fermi energy. The vector potential is defined by B = ∇ × A and using Maxwell's equation ∇ × B = µ 0 j, the current density must thus satisfy µ 0 j = ∇ × (∇ × A). Using both expressions for the current density the vector potential (and thus the magnetic flux B) can be solved.
To obtain the electronic density of states, one needs to determine the real-time, energy (ϵ) dependent Green functions rather than the imaginary-time ones. Using the Keldysh technique, the matrix Green function and Hamiltonian (used in the Usadel equation) are now written as: Here,ǧ R andǧ A are the retarded and advanced components describing equilibrium properties andǧ K is the Keldysh component which describes the non-equilibrium properties. Each component is a 4 × 4 matrix Green function in the Nambu ⊗ spin-space and a function of (x, ϵ). As long as the system is in equilibrium it is fully determined by the retarded component and it is thus sufficient to calculate the retarded equation, which effectively replacesǧ(x, ω) by g R (x, ϵ) and i ω n by ϵ. Once the retarded Green functions are obtained, the electronic density of states (N (x, ϵ)) follows from taking the real part of the (retarded) normal Green function.

B. STS
The scanning tunneling microscopy/spectroscopy (STM/S) measurements were performed using a home-built, cryopumped low temperature STM machine operating at a base temperature below 2 K. The samples were placed in the entry-lock chamber of the machine and pumped at a vacuum pressure of 8×10 −6 mbar overnight, before they were introduced to the STM image. The STM tip was made from PtIr wire, which before use was conditioned by field emission on a Au single crystal target. Bias voltage (V ) was applied to the sample with the STM tip at virtual ground. Differential conductance (dI/dV ) versus V spectroscopy data were recorded using the standard lock-in technique, with the frequency and amplitude of voltage modulation set at 413 Hz and 0.2 mV, respectively. The differential conductance of the STS spectra are modeled as [5]: with G 0 the differential conductance at eV >> ∆, N the (normalized) density of states (DOS) at the sample surface, f the quasiparitcle distribution function and V the applied voltage bias (with respect to the grounded sample) to the STM tip. For the N we use the results obtained from the quasiclassical theory and also make a comparison with the phenomenological Dynes theory [6] which approximates the DOS by: where ∆ is a measure of the (induced) superconducting gap and Γ the effective pair-breaking strength. Figure S1 shows in more detail the analysis presented in Fig. 1 of the manuscript and includes a direct comparison to fitting the data with the Dynes model. The latter predicts stronger smearing of the peak features than what is actually measured and gives less compelling fits compared to using a more realistic DOS obtained from quasiclassical simulations. Figure S2 shows ∆(T ) obtained from both approaches, where for the quasiclassical approach (where ∆ is solved as function of position) ∆(T ) is shown both at the sample surface and at the Nb/Pt interface. Note that ∆ as used in the quasiclassical theory itself is zero outside the superconductor and what is actually shown is the imaginary part of ∑ F ↑↓ , scaled such that for bulk Nb it gives the correct bulk value for the gap (∆ 0 ). The obtained full spatial dependence of ∆ (normalized to ∆ 0 ) for each temperature is shown in Fig. S3. Figure S4 shows the STM topographic image for the Au(5 nm)/Pt(2 nm)/Nb(50 nm) sample.

C. LEM
All low-energy muon-spin rotation (LE-µSR) measurements were performed on the µE4 beamline at the Paul Scherrer Institut [7]. The positive muon is an unstable spin-1 2 lepton of charge +e with a lifetime τ µ = 2.197µs and on decay emits a positron preferentially along its momentary spin direction. Upon implantation into a material a muon will rapidly thermalize, while maintaining its spin direction, after which its spin precesses around the local field making it a local magnetic probe. By monitoring the decay positrons of implanted, 100% spin polarized muons, information about the precession frequency (and thus the local field) can be obtained. The implantation depth profile (or stopping profile) of the muon is energy dependent and can be calculated by a well-proven Monte-Carlo simulation [8,9]. This allows the average probing depth to be tuned from about 10 to 100 nm below the surface. For a typical measurement several million counting events are collected at a rate of about 1 k/s and errors as small as 0.1 G can be achieved [10]. All our muon measurements were undertaken in transverse field geometry (applied field orthogonal to the muon spin direction) with the applied field direction in the plane of the sample. From the measurement data taken at a particular muon energy E (and thus a particular probing depth profile) one can determine the average flux ⟨B⟩(E), which can also be presented as ⟨B⟩(⟨x⟩) with ⟨x⟩ the average probing depth of the muons at energy E. However, in cases where the shape of the flux profile is known, either analytically including free fit parameters or numerically, one can treat the measurement data imposing this B(x) to find best fits to the data. To compare this to the conventional approach it is necessary to convolute B(x) with the muon stopping profiles in order to estimate which ⟨B⟩(E) the profile corresponds to. An example measurement of a muon experiment is shown in Fig. S5 (including the corresponding muon stopping profiles) where N L (t) and N R (t) are the measured positron detector histograms of the left and right detector respectively and t is the time at which a positron is detected. The two detectors are placed on opposite sides of the sample and the implanted muons have their initial spin direction pointing towards the left detector, with the result that from all the muons that decay at t = 0 a bigger fraction hits the left detector compared to the right detector. As a function of time, the muon-spin precesses around the local field which results in the muon-spin direction alternating between the left and right detector. The main exponential decay of the measurement signal reflects the muon life-time and the superimposed damped oscillations contain the information about the local field and define the asymmetry A(t) of the signal. The histograms can be fitted using N L, with N 0 the amplitude of the signal, K the time-independent background contribution and A(t) the asymmetry of the signal which carries all the information concerning the field distributions. It can generally be modeled as , with A 0 the (setup dependent) maximum asymmetry that can be measured, p(x) the muon stopping profile, γ µ = 851 Mrad·s −1 the gyromagnetic ratio for the muon, B(x) the local flux density, ϕ the starting angle of the muon-spin direction and G(t) the depolarization function (G ≤ 1), where the integral runs over the full width of the stopping profile In Fig. S6 we present muon data obtained on a Pt(60nm)/Nb(50nm). The results are similar to those obtained on sample SM3 of the main manuscript, but due to the superconductor being much thinner (50 nm instead of 96 nm for SM3), the effect (e.g. difference between the measurement in the normal state compared to the superconducting state) is also much reduced and only just visible. In Fig. S7 we present muon data obtained on a Pt(94nm)/Nb(96nm), which compared to sample SM3 has a thicker Pt layer. Up to the highest muon energy used, where muons reach up to some 20-25 nm into the Nb layer (see Fig. S7), no signs of a net Meissner screening is observed. This apparent lack of Meissner screening is consistent with our experimental findings on SM3 and related theory modeling. Due to the Pt layer being thicker (94 nm compared to 56 nm for SM3), higher muon energies are required in order for them to reach the Pt/Nb interface, which also results in a more spread-out stopping profile.

D. Triplet correlations induced by impurity spin-orbit interaction
For sample SM3 we use a detailed quasi-classical simulation including spin-orbit interaction in order to find a realistic solution for the local magnetic flux density, and from that, using the muon stopping profiles, the average flux density that would be observed during a muon experiment.
For the theoretical explanation of the experimental results we make use of the quasiclassical framework with extrinsic spin-orbit coupling (SOC) in the dirty limit [11][12][13] (coherence length much longer than mean free path), where the Green functions obey the Usadel equation. It is known that both intrinsic and extrinsic SOC induce singlet-triplet conversion at hybrid interfaces of usual singlet superconductors with nonsuperconducting materials, even in the absence of the exchange field [11,12,[14][15][16][17][18][19]. In the framework of the considered model triplets are generated by the Meissner currents in the presence of the impurity SOC, quantified by the spin Hall angle θ. The triplets are small and their amplitude relative to the singlet component is controlled by the parameter p s ξθ, where p s is the condensate momentum and ξ is the superconducting coherence length. This parameter is very small and our estimates give p s ξ ∼ 0.05 for the system under consideration (at the applied field B = 300 G and the Nb thickness 96 nm). Therefore, we can consider the triplet correction as a small perturbation. Then the general Usadel equation [11][12][13] can be simplified and up to the first order with respect to the parameter p s ξθ takes the form: [ whereǧ is the quasiclassical Green's function, which is a 4 × 4 matrix in the direct product of particle-hole and spin spaces. We use Pauli matricesτ i for the particle-hole space andσ i for the spin space. Eq. (6) is written in the Matsubara representation and ω n is the Matsubara frequency.∆ = |∆|τ 3 e iτ3χ τ 1 , where |∆| is the absolute value of the superconducting order parameter and χ is its phase, which in our case can be set equal to zero.
where A k is the k-coordinate component of the vector potential of the magnetic field.
where D is the diffusion constant and θ is the spin Hall angle of the material. The main generator of the triplets is the second term in the round brackets in Eq. (8). The applied magnetic field generates Meissner currents in the superconductor, which lead to the appearance of the triplets via this term. Please note that the sign of θ is changed with respect to Ref. [12]. In this case the parameter θ exactly coincides with the spin Hall angle, defined as in the limit of the nonsuperconducting system. Below we consider the process in detail.
If the magnetic field is applied along the z axis, the only nonzero component of the vector potential is A y and nonzero∂ yǧ appears accounting for the Meissner currents along the y direction. Then the second term in the round brackets in Eq. (8) results in nonzero spin component ∝ σ z of the matrix currentJ x . Further it serves as a generator of the opposite-spin triplets with S z = 0, which are described by the anomalous Green's function (f s + f t σ z ).
For the problem under consideration it is convenient to calculate the Green's function via the Θ-parametrization [20]. The matrix structure in the particle-hole space is written explicitly: The singlet contribution to the anomalous Green's function is contained in the first term of Eq. (9) and is represented by sinh Θ. The triplet contribution is contained in the second term and is represented by δΘ cosh Θ. The both quantities depend on the coordinate x along the direction perpendicular to the S/N interface. δΘ ∝ p s ξθ, as it is demonstrated below, and, therefore, small. Further we substitute the ansatz Eq. (9) into the Usadel equation Eq. (6) and up to the first order with respect to the parameter p s ξ obtain the following equations for Θ(x) and δΘ(x): Eq. (11) contains the term depending on the spin Hall angle θ, which could be a generator of the triplet correlations. However this term is typically small for type-II superconductors being of the order of ξ/λ with respect to all the other terms. The main generator of the triplets appears via the boundary conditions at the S/N interface. The boundary conditions [12] represent the generalized Kupriyanov-Lukichev boundary conditions and take the form: Here p s = −2(e/c)A y is the condensate momentum due to the Meissner currents and θ N,S are the spin Hall angles in the normal and superconducting parts of the structure. Eqs. (16) and (17), which represent boundary conditions for δΘ, contain terms ∝ p s θ N,S , which are generators for the triplets. The superconducting order parameter should be calculated self-consistently according to: The spatial distribution and the amplitude of the triplet anomalous Green's function f t = δΘ cosh Θ is presented in Fig. S9 for the Pt/Nb system for several Matsubara frequencies. As it should be for dirty samples, f t is an example of odd-frequency s-wave superconducting correlations. The lengths of layers are chosen to fit the experimental sample SM3. The resistance of the Pt/Nb interface is R b = 2 · 10 −15 Ωm 2 , the diffusion constants are D N b = 3.6 · 10 −4 m 2 /c and D P t = 1.1 · 10 −3 m 2 /c, the conductivities are σ N b = [1.5 · 10 −7 Ωm] −1 and σ P t = [1.1 · 10 −7 Ωm] −1 . In agreement with the structure of Eqs. (16) and (17) the amplitude of the triplet component at the interface should be of the order of f t (x = 0) ∼ θp s ξf s (x = 0). In our case f s (ω 1 ) ∼ 1 and p s ξ ∼ 0.05 at B = 300 G. Taking θ P t = 0.096 (see below) we obtain that the numerical value of f t,P t (ω 1 , x = 0) is in reasonable agreement with these rough estimates.
In principle, the triplet component contributes to the measured magnetic field in two different ways: (i) it results in the spin splitting of the quasiparticle DOS caused by the spin-triplet pairs. Although the triplet pairs have zero spin projection onto the z-axis, they are inevitably accompanied by the appearance of the spin dependence of the normal Green's function, as it is suggested by Eq. (9). In its turn, the spin-dependent correction to the normal Green's function leads to the spin splitting of the DOS, giving rise to the magnetization; (ii) it leads to the paramagnetic contribution to the Meissner currents. The distribution of the magnetic field in the system can be found making use of the Maxwell's equation: In the considered geometry Eq. (19) can be rewritten in the form: The first term in the right-hand side of Eq. (20) accounts for the the magnetization created by the triplet pairs and the second one describes paramagnetic contribution to the Meissner currents. We consider the both effects and have to conclude that the second one is negligibly small. The paramagnetic contribution to j s can be estimated in terms of the paramagnetic correction to the kernel Q = Q s (1 + δQ), where Q s is the singlet contribution to the Meissner response. The magnitude of δQ ∝ (p s ξ) 2 θ 2 , in fact it is even smaller because of additional reducing factors coming from the exact solution. From the numerical calculation it is obtained that in the framework of the considered mechanism it cannot account for the observed paramagnetic response and even does not make a visible contribution to it. Therefore, we focus only on the magnetization, induced by the triplet pairs. It can be calculated as where F a is the Landau Fermi liquid parameter [21] and can be expressed via the experimentally accessible normal state susceptibility χ = 2µ 2 B N F /(1 + F a 0 ). Then the magnetization takes the form The direct measurements of the low temperature normal state susceptibility of a Pt(40nm) film grown on a Si substrate give the result 4πχ = 0.104/(T + 0.652) in Gaussian units (dimensionless). The spatial profile of the corresponding paramagnetic field contribution δB = 4πM z , calculated according to Eqs. (20) and (22), is presented in Fig.4(c) of the manuscript. To compare the calculated profile with the muon results it is necessary to convolute δB(x) with the energy-dependent muon stopping profile. In order to compare δB(x = 0) at the interface with the experimental results presented in Fig.4(b) of the manuscript, we convolute δB(x) with the muon stopping distribution corresponding to 19(54)keV (nm) (approximately at the interface). Fitting the result of the convolution, which is denoted by ⟨B + ⟩ to the experimental value we obtain θ P t = 0.096. We also assume nonzero value of the impurity SOC in Nb θ N b = −0.2θ P t . Under this assumption we obtain small nonzero value of the triplet component in the depth of Nb, which is maintained by the Meissner currents, that is decays over the length scale of the magnetic field penetration depth. It results in the nonzero negative magnetization deep in Nb in Fig.4(c) of the manuscript. However, the small spin Hall effect in Nb does not influence the result qualitatively.
In order to compare the temperature of δB(T ) at a given x with the experimental results presented in Fig.5(a)-(b) of the manuscript we convolute B(x) with the muon stopping distribution corresponding to 19(54) keV(nm) (approximately at the interface) and 27(82) keV(nm) (in the depth of Nb). The result of the convolution of δB, which is denoted by ⟨B + ⟩, is presented in Fig. S10 as a function of temperature. The green curve corresponds to the convolution of δB with the muon distribution 19(54) keV(nm) and the black curve is for the convolution with 27(82) keV(nm) distribution. It is seen that ⟨B + ⟩ is strongly suppressed with increasing temperature. We see three reasons for this suppression. The strongest effect comes from the temperature suppression of the triplet correlations, the second reason is the temperature dependence of the normal state susceptibility χ. The third reason is associated to the fact that the decay length of δB in Pt, which is, roughly speaking, nothing but the normal state coherence length ξ N = √ D pt /2πT declines with increasing temperature. After the convolution procedure it leads to smaller ⟨B + ⟩. Fig.5(c) of the manuscript represents the result of the convolution of the total magnetic field, which is a sum of the diamagnetic response of the singlet component and the paramagnetic contribution of the triplet component.
The paramagnetic contribution becomes smaller with shortening the Pt layer. It is connected to the fact that the triplets generated at the opposite interfaces of the Pt layer have opposite signs if the directions of p s are the same at those interfaces. The same effect has been discussed in Ref. [11]. The paramagnetic contribution should become stronger in case of further increase of the Nb length. Roughly speaking, the effect is proportional to the Nb length unless it reaches 2 magnetic field penetrating depths of Nb. The point is that the triplet f t is proportional to the value p s (x = 0) at the Nb/Pt interface. For the structures under consideration the dependence p s (x) is nearly linear with the zero crossing point corresponding approximately to the middle of the Nb layer. Then p s (x = 0) = (2e/c)BL/2, where L is the Nb length.

E. Transport
Transport measurements were performed in a commercial 7 T helium flow cryostat using a He-3 insert (Cryogenic Mini Cryogen Free System). A standard four-point geometry was used to measure resistance (R) as function of temperature (T ) and applied field (H). Figure S11 shows the results obtained on the Cu/Pt/Nb (SM2) sample. The R(H) curves are measured both for decreasing and increasing field but no significant hysteresis was observed. From the field decreasing branches, at fixed temperature, and always starting from fields above the critical field to ensure normal state, the upper (perpendicular) critical fields are determined (H c2,perp ). The obtained H c2,perp (T ) are fitted to µ 0 H c2,perp (T = 0)(1 − (T /T c (H = 0)) α ) β and the Ginzburg-Landau coherence length is determined from ξ GL = √ /(2eµ 0 H c2,perp (T = 0)).

F. SEM/TEM
Samples for transmission electron microscopy (TEM) were prepared by a focused ion beam technique using a FEI Scios Dualbeam FIB. Sites of interest were identified by SEM imaging and protected from the damaging effects of the gallium ion beam by the electron-deposition of a carbon layer followed by the application of a thicker gallium iondeposited platinum layer. Thin foils were then prepared by ion milling and transferred in-situ to a copper TEM grid for final thinning and ion-polishing. TEM and scanning transmission electron microscopy (STEM) were carried out using a FEI Titan Themis S/TEM operated at 200 kV equipped with a Super-X windowless x-ray energy dispersive spectrometer (EDS). Bruker Esprit software was used for analysis and quantification of EDS spectra. Figure S12 shows the EDS analysis of the Pt/Nb sample (SM3).