Proximity-induced odd-frequency superconductivity in a topological insulator

At an interface between a topological insulator (TI) and a conventional superconductor (SC), superconductivity has been predicted to change dramatically and exhibit novel correlations. In particular, the induced superconductivity by an $s$-wave SC in a TI can develop an order parameter with a $p$-wave component. Here we present experimental evidence for an unexpected proximity-induced novel superconducting state in a thin layer of the prototypical TI, Bi$_2$Se$_3$, proximity-coupled to Nb. From depth-resolved magnetic field measurements below the superconducting transition temperature of Nb, we observe a local enhancement of the magnetic field in Bi$_2$Se$_3$ that exceeds the externally applied field, thus supporting the existence of an intrinsic paramagnetic Meissner effect arising from an odd-frequency superconducting state. Our experimental results are complemented by theoretical calculations supporting the appearance of an odd-frequency component at the interface which extends into the TI. This state is topologically distinct from the conventional Bardeen-Cooper-Schrieffer (BCS) state it originates from. To the best of our knowledge, these findings present a first observation of bulk odd-frequency superconductivity in a TI. We thus reaffirm the potential of the TI/SC interface as a versatile platform to produce novel superconducting states.

At an interface between a topological insulator (TI) and a conventional superconductor (SC), superconductivity has been predicted to change dramatically and exhibit novel correlations. In particular, the induced superconductivity by an s-wave SC in a TI can develop an order parameter with a p-wave component. Here we present experimental evidence for an unexpected proximity-induced novel superconducting state in a thin layer of the prototypical TI, Bi2Se3 proximity-coupled to Nb. From depth-resolved magnetic field measurements below the superconducting transition temperature of Nb, we observe a local enhancement of the magnetic field in Bi2Se3 that exceeds the externally applied field, thus supporting the existence of an intrinsic paramagnetic Meissner effect arising from an odd-frequency superconducting state. Our experimental results are complemented by theoretical calculations supporting the appearance of an odd-frequency component at the interface which extends into the TI. This state is topologically distinct from the conventional Bardeen-Cooper-Schrieffer (BCS) state it originates from. To the best of our knowledge, these findings present a first observation of bulk odd-frequency superconductivity in a TI. We thus reaffirm the potential of the TI/SC interface as a versatile platform to produce novel superconducting states.
In a conventional superconductor (SC), the electronic excitations are usually described as condensation of Cooper pairs [1]. Fermi statistics implies symmetry constraints on permutation properties of the pair wave function, thus limiting the possible SC classes [2]. According to conventional classification, states with even parity (s-, d-wave) must be in a spin-singlet configuration while states with odd parity (p-, f -wave) must be in a spin-triplet configuration. However, additional classes are possible when permutation with respect to time and, if present, orbital degrees of freedom are included. This general classification allows for the oddfrequency, or Berezinskii, state characterized by superconducting pairing which is nonlocal and odd in time [2][3][4]. Odd-frequency pairing gives rise to SC states with symmetries different from conventional states, for example triplet s-wave and singlet p-wave states. The Berezinskii state is currently recognized as an inherently dynamical order that can be realized in a variety of systems, including bulk SCs, heterostructures, and dynamically driven systems. It is especially relevant at interfaces, where locally broken symmetries can influence the type of pairing [4,5].
One of the peculiarities of such a state is that the sign of the Meissner screening can be reversed in some cases, causing an attraction instead of a repulsion of external magnetic fields. Paramagnetic Meissner response has been predicted in proximity structures [4,[7][8][9][10][11] and in multiband SCs [12]. While such a paramagnetic Meissner screening cannot be stable in the bulk, it has been observed at interfaces, where the superconductivity is induced in a non-superconducting layer by the proximity effect [13]. One particular superconducting interface that has attracted significant attention in recent years is between a conventional s-wave SC and a 3D topological insulator (TI). Fu and Kane have predicted that induced superconductivity in the topological surface state (TSS) may have a p x +ip y even-frequency order parameter, that might allow stabilization of Majorana bound states in vortex cores [14]. The latter are the key ingredient in a proposal for fault-tolerant quantum computation. Furthermore, the occurrence of Majorana zero modes can be related to the presence of odd-frequency superconducting components [15][16][17][18]. Aside from the conventional evenfrequency superconductivity one should be on a lookout for odd-frequency correlations induced at such interfaces.
In this work we provide experimental evidence of proximity induced odd-frequency superconductivity in a heterostructure of Bi 2 Se 3 (∼ 110 nm) on Nb (∼ 80 nm). In particular, we measure the depth dependence of the magnetic field parallel to the interface in the Meissner state using low energy muon spin rotation [19][20][21][22] (LE-µSR) at the µE4 beam line [23]  µSR geometry where the applied field is parallel to the surface and perpendicular to the muon spin and momentum directions. The lines depict the expected magnetic field depth profile due to screening of the applied field inside the heterostructure with (solid gray) and without (dotted black) proximity induced superconductivity in the Bi2Se3 layer. (bottom) Calculated muon implantation profiles at different implantation energies [6]. (b) Measured local mean field as a function of implantation energy above (black squares) and below (red circles) the superconducting transition temperature. The black horizontal line represents the applied field value, Bext, and the blue line is a fit to the theoretical model (see text).
sion characterization of the magnetic field profile by measuring the average Larmor precession frequency of the muons' spins as a function of their implantation energy (and corresponding implantation depth), as illustrated in Fig. 1(a). Below the superconducting transition temperature of Nb we observe conventional diamagnetic Meissner screening in the Nb layer. In contrast, a paramagnetic Meissner screening is observed in the Bi 2 Se 3 layer, indicating a proximity-induced odd-frequency superconducting component in the TI layer. Our experimental results are complemented by theoretical calculations supporting the appearance of an odd-frequency component at the interface, extending deep into the TI.
The local magnetic field as a function of implantation energy, E, is shown in Fig. 1(b). At a temperature of 20 K, i.e., in the normal state, the muons probe a depthindependent magnetic field (black squares). At temperatures below the superconducting transition of Nb, we observe a strong variation in the measured field as a function of depth, featuring a typical Meissner screening inside the Nb layer (red circles). In a conventional metal/SC proximity structure, the induced superfluid density in the metal results is a decreased mean field which increases monotonically as a function of distance from the metal/SC interface, reaching the applied field value far inside the metal [24][25][26]. In contrast, the behaviour in Bi 2 Se 3 /Nb is reversed; the field in the TI is enhanced compared to the applied field value, B ext (i.e., its value at 20 K). This can be clearly seen in the fast Fourier transform spectra of the muon spin ploarization spectra which represents the local field distribution sensed by the muons [ Fig. 2(a)]. A careful inspection of the temperature dependence of this effect shows that the paramagnetic field shift in Bi 2 Se 3 occurs below the superconducting transition temperature, T c ≈ 9 K, of Nb [ Fig. 2(b)]. This is in agreement with T c obtained from transport measurements which show a sharp superconducting transition below 9 K (see Supplementary Material (SM) [27]). We have also measured the magnetic field as a function of depth for various B ext values in the range 7 mT to 20 mT. We find that the local field at low temperatures is proportional to the applied field, always exhibiting a paramagnetic shift inside the TI. This indicates that the induced superfluid density is almost field-independent within this field range (see SM [27]).
The observed paramagnetic shift of about 0.2 % is too large to be attributed to the demagnetization fields of the Nb layer, which in a thin film are only relevant close to the edges [27]. Furthermore, a misalignment of the applied field with respect to the interface could only reduce the mean field measured in the sample. Another possible source for a positive shift are microscopic demagnetization fields caused by the roughness of the SC/TI interface. However, such stray fields should decrease exponentially with distance away from the interface [28,29], which does not agree with our observation. We can also exclude a temperature dependence of the field in bulk Bi 2 Se 3 or any systematic deviations caused by the ex- (a) Local field distribution calculated from fast Fourier transform of the muon polarization measured in the heterostructure, above and below the superconducting transition temperature. Different energies are offset for clarity. The solid lines represent results from the fitting procedure (see SI [27]). (b) Temperature dependence of the field shift in the Bi2Se3 and Nb layers.
Odd-frequency components in superconductors are expected at any ballistic interface [5] and in TIs, they may also occur in the presence of in-plane gap gradients [31], exchange fields [32] and multigap odd-orbital coupling [33]. However, they should be a subdominant component compared to the even-frequency pairs. The total shielding current given by j = −e 2 (n e − n o )A/mc, where n e/o is the superfluid density of the even/odd pairs and A the vector-potential, should therefore still be diamagnetic [10]. The observed extension of the positive shift deep into the Bi 2 Se 3 layer, shown in Fig. 1(b), is a clear indication that not only the TSS but also the bulk conduction band of Bi 2 Se 3 is relevant for the proximity effect, as has also been pointed out by previous studies [34]. Therefore, we conclude that the observed paramagnetic Meissner screening is due to supercarriers induced into the bulk conduction band of Bi 2 Se 3 by the proximity to superconducting Nb. Furthermore, this observation suggests that the proximity-induced superconducting state is the unconventional odd-frequency (Berezinskii) state. Such a nontrivial pairing state in TIs has not been observed before.
To make a better comparison and illustrate the veracity of our conclusion we provide a summary of theoretical calculations that support our experimental results. We have developed a theoretical description of proximity-induced superconductivity in the TI/SC heterostructure based on a well established two-orbital tightbinding model for Bi 2 Se 3 [33,35] (see SM [27]). We show that the induced odd-frequency pairing persists in the bulk of the TI, and that it dominates the Meissner effect near the TI/SC interface. Based on the insights from this microscopic model, we propose a theoretical depth profile of the magnetic field, derived from Maxwell's equations and linear response theory, that quantitatively fits our experimental data.
The band inversion in Bi 2 Se 3 implies that the bulk conduction and valence bands are formed by two orbitals with different parity, originating from hybridized Se and Bi p z states. This orbital degree of freedom in the TI allows the generation of odd-frequency pairing components, in addition to the even-frequency ones [36]. For a proximity-coupled TI to a singlet s-wave SC, the symmetry allowed odd(even)-frequency components are odd(even) in the orbital index [33]. In addition to the dominant s-wave singlet pairing, even-and oddfrequency p-wave triplet components are present. Our calculations show strong odd-frequency SC correlations that propagate away from the interface due to the coupling between Nb and the bulk electronic states of Bi 2 Se 3 . The two orbitals in the tight-binding model correspond to top and bottom Se p z states in a quintuple-layer, each hybridized with the neighbouring Bi atom [37]. We assume stronger tunneling from the SC into the orbital closest to the interface. This gives rise to odd-frequency components that can be comparable or larger in magnitude than the even-frequency ones over a wide range of frequencies.
The penetration depth of the induced SC pairing in the TI layer depends on the position of the chemical potential. When the chemical potential is at the Dirac point of the TI, the induced SC pairing amplitudes decay within approximately two quintuple layers from the interface, i.e., a typical penetration depth of the TI surface states. However, when the chemical potential is in the bulk conduction band, as is the case for an (intrinsically) n-doped Bi 2 Se 3 samples used here, we find finite SC pairing amplitudes far away from the interface. We attribute this to coupling between the SC and bulk TI states. This is consistent with the experimental observation of the Meissner screening inside the TI layer shown in Fig. 1(b). The role of bulk states in the TI/SC heterostructure has been previously pointed out for trivial even-frequency pairing [34]. We extend this approach to allow for odd-frequency pairing induced in the TI.
We will now focus on the interpretation of the positive field shift in the observed Meissner screening. In ordinary metals or semiconductors, odd-frequency pairing leads to the paramagnetic Meissner effect [12]. Recently, it has been pointed out that for materials with Dirac dispersion both even-and odd-frequency components can contribute to dia-and paramagnetic Meissner effect [38,39]. This complication, however, is less relevant for the system considered here since the coupling between SC and TI bulk states is dominant. Therefore, we attribute the paramagnetic shift to the odd-frequency components. Based on these considerations, we develop a phenomenological model of the Meissner effect.
We consider a TI/SC heterostructure with the interface at z = 0, SC (Nb) extending from −d SC to 0 and the TI (Bi 2 Se 3 ) from 0 to d TI . Similarly to Ref. 13, we solve a differential equation for the vector potential where K xx (z) is the current-current correlation function, or the Meissner kernel, which determines the magnetic response of the system. The local magnetic field is calculated as B(z) = dA(z)/dz. We use matching boundary conditions for A(z) and B(z) at the TI/SC interface and set B(z) = B ext outside the heterostructure. Furthermore, we propose the following form of the Meissner kernel in the TI and SC layers where λ SC is the London penetration depth of the SC while λ TI and z 0 are the characteristic lengthscales in the TI layer. Thus, we assume a conventional diamagnetic Meissner screening in the SC which leads to exponential suppression of the local magnetic field extending into the SC layer. On the TI side, the kernel describes the proximity-induced odd-frequency screening which exponentially decays over a persistence length z 0 from the TI/SC interface. We solve Eq. (S25) with the kernel K xx given in Eq.
(2) and with boundary conditions at all three interfaces to obtain the theoretical depth profile of the local magnetic field, B theo loc (z). To fit the experimental data, we use B theo loc (z) to calculate the muon precession signal averaged over its stopping distribution at a given E (see details in SM [27]). The minimization of the fit was performed for all E values (at T = 3.7 K) simultaneously, using the same theoretical field profile. The parameters extracted from this global fit [gray line in Fig. 1(a)] give λ TI = 1.62(4) µm and λ SC = 117.8(6) nm. In this fit z 0 is chosen to reflect the fact that superconductivity is induced in the Bi 2 Se 3 bulk conductance states above the TI gap and, therefore, has a long decay range inside the TI. We checked that the specific choice of z 0 does not affect the quality of the fit. However, the exact value of of λ TI systematically depends on the choice of z 0 and varies between 0.90(2) µm and 1.74(5) µm for z 0 within [d TI , ∞). The penetration depth λ SC is set by the penetration depth in Nb, which we assume is not modified significantly by the interface. Note that the penetration depth obtained for the Nb is much larger than that in the clean limit (27(3) nm [20]). We also fixed the thickness of the layers obtained from RBS and XRR measurements, d TI = 111.1 nm and d SC = 81 nm, together with the applied field obtained from the measurements at T = 20 K, B ext = 19.818(2) mT. Finally, we calculated the mean local field as a function of E by averaging B theo loc (z) over the corresponding muons' stopping distribution. These are plotted in Fig. 1(b) (blue line) and exhibit a very good agreement with the mean field obtained independently, confirming again the validity of the obtained B theo loc (z). Hence, the proposed theoretical field profile, based on the assumption of a large odd-frequency pairing amplitude in the bulk of the TI supported by the microscopic model, explains qualitatively and quantitatively the observed paramagnetic Meissner shift. The measured local field profile shows an almost constant paramagnetic shift in the TI extending at least 30 nm from the TI/SC interface and a conventional screening on the superconducting Nb side. The paramagnetic shift decreases gradually towards the opposite TI surface. Deviations of the experimental data from the predicted profile, for example a small positive shift of the mean magnetic field above B ext on the SC side just below the interface may be due to the limited accuracy in calculating the stopping profiles using Trim.SP [ Fig. 1(a)]. Furthermore, some of the experimental details that may affect the local field profile, such as interface roughness and inhomogeneous thickness of the layers, are not included in the theoretical modelling. However, we except this to have a negligible effect, since the roughness is on an atomic length scale. Therefore, the magnetic field response is set by the macroscopic length scales, such as the penetration depth. We also point out that both, the induced dominant odd-frequency superconductivity and the underlying Nb superconductivity, are of s-wave type. Therefore the effects of atomic disorder at the interface are not likely to degrade the induced odd-frequency state.
In conclusion, we observe an intrinsic paramagnetic Meissner shift in proximity-induced superconductivity in a Bi 2 Se 3 /Nb heterostructure. We attribute this effect to odd-frequency superconductivity which persists up to tens of nanometers away from the Bi 2 Se 3 /Nb interface. This finding, which is supported by our theoretical calculations, is the first observation of a bulk induced odd-frequency superconducting state in a TI. Our results demonstrate that the experimental phenomenology of superconductivity at TI interfaces is richer than previously thought, and it highlights the potential of TI/SC heterostructures for realizing novel electronic states. Odd-frequency superconducting compo-nents may be of particular importance for the theoretical description of TI and semiconductor-based Majorana heterostructures which operate in a non-zero magnetic field.
Supplemental Material: Proximity-induced odd-frequency superconductivity in a topological insulator

EXPERIMENTAL METHODS
The sample was prepared by sputtering a Nb layer on top of sapphire substrate and transferred ex situ to a molecular beam epitaxy (MBE) chamber (Createc GmbH). The Nb surface was then cleaned from any oxide impurities, followed by MBE growth of a Bi2Se3 layer on top of Nb at 200 • C, and finally a cool-down under Se flow to 25 • C [S1]. A standard effusion cell was used for the deposition of Bi (99.9999% purity) while Se (99.9999% purity) was sublimated out of a cracker source (Createc GmbH). A Bi:Se flux ratio of 1:10 was used to reduces the formation Bi2Se3 anti-site defects and Se vacancies.
Low energy muon spin spectroscopy measurements were performed on the µE4 beamline of the Swiss Muon Source at Paul Scherrer Institute in Villigen, Switzerland [S2]. In these measurements, fully polarized muons are implanted at a given implantation energy, E, and the evolution of their polarization is monitored via their anisotropic beta decay. The probing depth inside the sample was tuned by varying E in the range 3 keV to 25 keV. The decay positrons were detected by a set of detectors placed around the sample region. The sample was cooled with a 4 He flow cold finger cryostat. To suppress the signal from muons missing the sample, the backing plate of the sample was sputtered with a thin layer of Ni [S3]. The muon stopping profiles in the sample were modeled with the Trim.SP code [S4]. As input to these calculations we used the area density measured with Rutherford backscattering (RBS) and the thicknesses determined with X-ray reflectometry (XRR), as shown below. The measured µSR asymmetry (polarization) spectra were analyzed using the Musrfit software [S5]. In a magnetic field applied perpendicular to the initial muon spin direction, those time spectra can be understood as the Fourier transform of the local field distribution averaged over all muon stopping sites [S6].

µSR DATA ANALYSIS
To confirm that the observed behavior is not an artifact of the analysis procedure, we have analyzed the raw data using two independent approaches. In the first approach, we parametrize each individual spectrum separately to calculate the mean field sensed by the muons at a given E and temperature, T . In the second approach, we use the theoretical field profile to model the muon spin polarization for all values of E and T simultaneously, using the corresponding stopping distributions. Both methods are described in detail below and a schematic overview is given in Fig. S1.

Parametrization of individual µSR spectra
The asymmetry of the individual spectra can be fitted to an exponentially damped oscillation, where A0 is the initial asymmetry, σ is the depolarization rate, γµ = 2π × 135.5 MHz/T the muon's gyromagnetic ratio and ϕ is a detector-dependent geometrical phase factor. The Nb nuclear magnetic moments are significantly larger and have higher natural abundance than those of Bi and Se. This leads to a different muon spin depolarization rates in the Nb and Bi2Se3 layers. This effect is of little importance above the superconducting transition temperature where the measured asymmetry can still be described by Eq. (S1) at all implantation energies. However, in the superconducting state, Eq. (S1) fails to fit the data due to the strong depth dependence of the field across the sample, in particular at E > 15 keV, as can be clearly seen in Fig. 2(a) in the main text. Therefore, we used a sum of two exponentially damped signals to fit spectra measured with E ≥ 15 keV and T ≤ 8 K, Here the subscripts TI and Nb indicate which layer the signal originates from. The mean field is then given by Note that the exact details of the fit do not affect the obtained value of Bmean significantly as long as the fit mimics the spectra accurately. The values ob Bmean as a function of E are plotted in Fig. 1(b) in the main text.

Fit of the phenomenological field distribution
Instead of parametrizing individual µSR spectra, it is possible to directly model them using a single depth profile of the magnetic field B(z), see e.g. Ref. [S7, S8]. This requires knowledge of the muon stopping distribution n(z, E), which is obtained from the Trim.SP simulations. For this analysis we have taken the experimentally determined layer thickness from RBS and XRR (see below), dTI = 111.1 nm and d Nb = 81.0 nm. At a given depth, z, we can write the time evolution of the muon spin polarization as Therefore, the µSR asymmetry spectrum at a given E is calculated by performing the ensemble average over z, where A0 is the initial asymmetry. To fit the experimental data, we evaluated this integral discreetly on a grid of 0.5 nm. At 20 K we assumed a constant field inside the sample, i.e. B(z) = Bext. In the superconducting state, we assumed B(z) to be the solution to the phenomenological Meissner kernel from Eq. (2) in the main text, using the externally applied field Bext as a boundary condition. The details of the functional form of B(z) are given below in Eqs. (S32) and (S35). For consistency, we fixed Bext, σTI and σ Nb to their values obtained from the fits at 20 K (Table I). Furthermore, we assumed a persistence length z0 = 1 µm in all fits, where z0 describes the propagation of the Nb electrons into the TI. This large z0 value corresponds to proximity effect with the chemical potential of Nb lying within the bulk conduction band of Bi2Se3. On the other hand, if the chemical potential were within the TI band gap, there should only be proximity to the surface states and z0 would be on the scale of the lattice parameter. Our fitting procedure gives reasonable parameters only for large z0 values, although the quality of the fit is not affected by the specific choice of z0. The results from these global fits of all spectra simultaneously at 20 K and at 3.7 K for different Bext are listed in Table I. The errors quoted for these parameters reflects only the statistical errors. For example, the specific choice of z0 does not affect the quality of the fit but changes the value of λTI. z0 between dTI and ∞ results in λTI between 0.90(2) µm and 1.74(5) µm.  In order to compare the measured Bmean (Eq. S3) with that calculated based B(z) from the global fit, we use,

K 3.7 K
This is shown as solid lines in Fig. S2 and Fig. 1(b) of the main text.

Rutherford backscattering
The structure and stoichiometry of the sample was verified with RBS at the Tandem-accelerator of ETH Zurich. Measurements were performed with 5 MeV 4 He 2+ primary ions and a Si PIN diode detector placed at 168 • . Figure S4 shows the RBS yield as a function of the final He energy. The spectrum was analyzed with the RUMP software [S9]. We found that the area densities of Bi2Se3 and Nb are 3.75(15) × 10 17 at /cm 2 and 4.28(17) × 10 17 at /cm 2 , respectively.

X-ray reflectometry
To study the thickness of the layers we performed XRR on a Seifert four-circle diffractometer, using the Cu-Kα line. The angle dependence of the reflected intensity is shown in Fig. S5. The corresponding scattering length densities were calculated with GenX from the interface, on a length-scale of ∼1 nm [S11]. This is an order of magnitude shorter than the observed superconducting proximity.

Tight-binding model
We consider a tight-binding model on a three-dimensional (3D) cubic lattice which describes a 3D topological insulator (TI) coupled to a superconducting layer (SC). The Hamiltonian of the system is written as where H TI(SC) is the Hamiltonian of the TI(SC) layer and Ht is the coupling Hamiltonian. To model the electronic structure of Bi2Se3, we include two orbitals per atoms and spin-orbit interaction [S12-S14] . The two orbitals represent two Se pz-orbitals, originating from top and bottom Se layers in a quintuple-layer unit cell of Bi2Se3, each hybridized with its neighbouring Bi pz-orbital. The hybridized Se and Bi orbitals give rise to the inverted valence and conduction bands of opposite parity. The conduction(valence) band has a predominant contribution from Se(Bi) derived states. In the bulk, the two orbitals can be transformed into each other by inversion with respect to the central Se atom. However, in the TI/SC heterostructure the inversion symmetry is broken and the two orbitals are shifted along the z-axis with respect to the inversion center. In order to reflect this asymmetry in the tight-binding model for a heterostructure, we assume a larger tunnelling amplitude between the atoms of the SC layer and one of the orbitals, corresponding the hybridized orbital closest to the interface. The Hamiltonian of the TI layer can be written in the spin and orbital basis as where Φ † = (a ↑ † 1k , a ↓ † 1k , a ↑ † 2k , a ↓ † 2k ), a σ † τ k (a σ τ k ) is the creation (annihilation) operator for an electron with spin σ, orbital τ = 1, 2 and momentum k inside the TI layer. The Hamiltonian matrix HTI(k) is given by, where the diagonal elements hαα(k) = γ0 − j γi(e ik j a + e −ik j a ) (α = 1, ..., 4, j = x, y, z) represent same-orbital (intra-orbital) same-spin hopping, off-diagonal elements in spin index h12(k) = −h34(k) = −λx(e ikxa − e −ikxa ) + iλy(e iky a − e −iky a ) represent intra-orbital spin flips, and off-diagonal elements in orbital index h13(k) = h24(k) = − j tj(e ik j a +e −ik j a )+λz(e ikz a −e −ikz a ) represent inter-orbital same-spin hopping. Here, a is the lattice constant and k ≡ (k || , kz), with k || = (kx, ky) is the twodimensional (2D) (in-plane) momentum. The Hamiltonian can be rewritten in concise form as, whereÎ is the identity matrix in spin and orbital space, d4 = γ0 −2 j γjcos(kja), d0 = ε−2 j tjcos(kja), dj = −2λj sin(kja), Γ0 = τx ⊗ σ0, Γx = −τz ⊗ σy, Γy = τz ⊗ σx, and Γz = τy ⊗ σ0, where τ and σ are Pauli matrices. The parameters of the model, γ0, γj, tj, λj and (with j = x, y,z) are fitted to point [S13]. We use the following parameters: γ0 = 0.3391, γx,y = 0.0506, γz = 0.0717, ε = 1.6912, tx,y = 0.3892, tz = 0.2072, λx,y = 0.2170, and λz = 0.1240 eV [S13]. The eigenstates of HTI can be calculated analytically and are given by E± = d4 ± µ d 2 µ . For the Nb layer, we consider a single layer of a singlet s-wave SC, with the Hamiltonian given by where k ≡ k || = (kx, ky) is the in-plane momentum, d σ † k (d σ k ) is the creation (annihilation) operator for an electron with spin σ and momentum k, ε k = −2[cos(kxa) + cos(kya)] is the dispersion, and ∆0 is the superconducting order parameter in the SC layer.
We assume nearest neighbors hopping between the atoms of the SC and the TI interfacial layer, with different tunneling amplitudes for the two orbitals of the TI atoms. The tunneling Hamiltonian is given by, where tτ , with τ = 1, 2, are the orbital-dependent tunneling matrix elements. Here, the creation (annihilation) operators for TI are labeled by the in-plane momentum k ≡ k || and the atomic index iz, where iz = 0 for the interfacial layer. As described above, in order to reflect the spatial asymmetry between the two TI orbitals, we allow the tunnelling amplitudes t1 and t2 to be different (for the sake of definiteness, we assume t1 < t2).

Proximity-induced Berezinskii paring at TI/SC interface
Due to proximity to the SC layer, superconductivity is induced inside the TI. Based on the symmetry analysis for an s-wave SC and a multiband model of a TI with tetragonal symmetry [S13], the possible induced pairing in the TI are even-frequency even-orbital and odd-frequency (Berezinskii) odd-orbital spin-singlet and spin-triplet pairings. In addition to the dominant s-wave spin-singlet pairing, even-and odd-frequency p-wave triplet components are present. In order to find the pairing amplitudes, we calculate the proximity-induced anomalous Green's function at the interface. We define the electron, hole and anomalous Green's functions as functions of momentum k and imaginary time t, From the spin-and orbital-resolved anomalous Green's function, we present the singlet and triplet pairing amplitudes that are either odd or even in the orbital index. We introduce the following singlet/triplet odd-frequency odd-orbital (F τ τ (o) s/t ) and even-frequency even-orbital (F τ τ (e) s/t ) pairing amplitudes, where τ, τ = 1, 2 andŜ s/t is the spin-dependent symmetry factor. Here,Ŝs(k) = iσy for the singlet component andŜt(k) = 2i(d · σ)σy for the triplet component, with d a vector in spin space which specifies the symmetry of the triplet state. For instance, d = (kx, ky, 0) for the triplet state with A1u symmetry and d = (ky, −kx, 0) for a triplet state with A2u symmetry. N is the normalization factor, equal to the number of k-points in the Brillouin zone. The only non-zero odd-frequency component is F

12(o)
s/t , i.e. the inter -orbital component which is odd in the orbital index. Both even-frequency intra-and inter-orbital components, F τ τ (e) s/t (τ = 1, 2) and F 12(e) s/t respectively, are present. We find numerically, that the even-frequency intra-orbital component is dominant, therefore we compare F in order to estimate the magnitude of the odd-frequency pairing. To the leading order in the tunneling matrix elements, the anomalous Green's function in the TI is given by, whereĜ TI 0 is the Green's function of TI in the normal state, i.e. decoupled from the SC layer, andF SC 0 is the anomalous Green's function of an isolated SC layer. Here all Green's functions are matrices in spin (σ), orbital (τ ) and atomic indexes (i ≡ ix, iy, iz), e.g.F TI ≡ [F T I ] σσ τ τ ,ii . For the inter-orbital pairing amplitude, we get whereTτ | = tτÎ andÎ is the identity matrix. After performing the Fourier transformation, the inter-orbital pairing amplitude at the interface becomeŝ where we used shorthand notation for the in-plane atomic indices, i = (ix, iy). Note that hereF TI 12 is a matrix in spin basis. Other pairing amplitudes in the orbital basis are calculated in a similar way.
We also performed calculations using the tight-binding Hamiltonian of a TI as described above, for a finite slab of TI coupled to a single layer of SC. The induced pairing amplitudes calculated with this approach are in qualitative agreement with the semi-analytical approach based on the bulk model and second-order perturbation theory [see Eq. (S21)]. However, the semi-analytical model is more convenient for calculating the current-current correlation function that determine the magnetic response of the system. The odd-and even-frequency pairing amplitudes at the interface (interfacial TI layer), calculated using this approach, are shown in Fig. S6. The odd-frequency component is comparable to or larger in magnitude than the evenfrequency ones over a wide range of frequencies. Figure S7 shows the odd-frequency pairing amplitude for different tunneling ratios t1/t2. The pairing amplitude decreases with increasing t1/t2, i.e. with decreasing the asymmetry in the coupling between the atoms of the SC and the two TI orbitals, which is the source of the odd-frequency pairing.
Enhancement of Berezinskii pairing due to coupling to bulk TI The depth dependence of the even-and odd-pairing amplitudes for different positions of the chemical potential in the TI is presented in Fig. S8. We find that SC pairing extends deeper into the bulk TI when µTI is shifted away from the Dirac point (µTI = 0), i.e., into the conduction or valence band. This is due to the fact that for a chemical potential larger than the bulk TI bandgap, both the topological surface states and the bulk TI states are present at the chemical potential and can therefore couple to the SC. This supports the experimental result that SC exists not only at the interface where the TI surface states reside, but also in the bulk with a large persistence length, z0, such as that used in our fits. This behavior has previously been seen for the usual even-frequency pairing in the TI/SC heterostructures using momentum-resolved photoemission spectroscopy [S15]. Here, we extend this result to odd-frequency pairing. The details of the interplay between the Dirac node and the Berezinskii state will be explored more in a separate paper [S16]. Here we quote our finding that the Berezinskii order-parameter defined as ∆Ber ≡ ∂F (o) (ω) ∂ω |ω=0 grows into the bulk TI states as the doping pushes the chemical potential away from the node.

Current-current correlation function at the TI/SC interface
To justify the phenomenological description of the paramagnetic Meissner effect in the TI/SC heterostructure, we calculate the current-current correlation function from linear-response theory. In the case of static and uniform magnetic field, the Meissner kernel, Kxx(ω, q), is given by whereĵ(k) = ∂Ĥ T I (k) ∂k . Here,ĵx,F TI andĜ TI are matrices in spin and orbital basis. The Green's function for the TI in proximity to the SC layer is calculated using second order perturbation theory as explained in the previous section.
For the calculation based on the bulk model of the TI, Kxx(0, 0) gives the current-current correlation function at the interface and is a constant dependent on the model parameters. Taking the corresponding parameters for the TI Bi2Se3, we find Kxx(0, 0) = −κ 2 < 0, which is consistent with the possible paramagnetic Meissner effect.

THEORETICAL PARAMETRIZATION OF THE LOCAL MAGNETIC FIELD PROFILE
We consider the TI/SC heterostructure with interface at z = 0, Nb extending from −dSC to 0 and TI from 0 to dTI. The external magnetic field Bext is applied along the y-axis, Bext = (0, Bext, 0); z-axis is the TI quintuple layer growth direction. For a phenomenological description of the depth profile of the local magnetic field across the TI/SC heterostructure, we use a combination of Maxwell's equations and a linear response theory [S17].
From Maxwell's equations in transverse gauge, namely ∇ 2 A = −J, B = ∇ × A for ∇ · A = 0, where A is the vector potential, B is the local magnetic field and J is the current density, we obtain the equation for the vector potential From linear response theory, the vector potential and the current density are related by where Kxx(z) is the current-current correlation function, or the Meissner kernel. Substituting Eq. (S24) into Eq. (S23), we get The local magnetic field inside the TI is given by B(z) = Bext + M (z) and from Maxwell's equations B(z) = dAx(z)/dz. We consider the following boundary conditions at the interfaces (S28)

Sign of the Meissner effect
For now, we will neglect the depth dependence of the current-current correlation function and assume that Kxx(z) =const. For even-frequency superconducting pairing, Kxx(z) = κ 2 > 0, which gives a diamagnetic Meissner effect. The solution of Eq. (S25) gives, Here we have normalized the z-coordinate to the TI layer thickness, z → z/dTI.
It can be shown analytically for systems with quadratic dispersion that for any odd-frequency superconducting pairing, regardless of its origin, Kxx(z) = −κ 2 < 0 [S18]. The magnetization is given by, Hence, the magnetization and, as a result, the induced magnetic field B(z) inside the TI are an oscillatory decaying function of the distance. Therefore, B(z) can be larger or smaller than Bext and the exact shape of the field profile is determined by material parameters. However, in order to have a paramagnetic shift, the sign of Kxx must be opposite to that of the conventional SC. As shown in the previous section, the calculation based on microscopic tight-binding model and linear response theory yields Kxx(z = 0) < 0 at the interface. This further supports the dominant odd-frequency pairing assumed in the phenomenological model. It should be mentioned that recent theoretical studies on odd-frequency pairing in Dirac materials showed that due to inter-band terms in the Dirac Hamiltonian, both even-and odd-frequency components can contribute to diamagnetic and paramagnetic Meissner effect [S19, S20]. This means that the paramagnetic Meissner effect in Dirac materials is not necessarily due solely to the odd-frequency components. However, in the system studied here, the induced SC pairing has bulk nature due to the electron doping in the sample which shifts the chemical potential away from the Dirac node and into the bulk TI conduction band. Therefore, we attribute the paramagnetic Meissner shift to the odd-frequency pairing.
For the TI/SC heterostructure, the induced superconducting pairing has both the even and odd-frequency components, as shown in the theoretical calculations. In order to have the overall paramagnetic response, the odd-frequency component should be dominant. In addition, the oscillating odd-frequency component should produce a positive shift of the local magnetic field with respect to Bext. Below we construct a theoretical parametrization of magnetic field depth profile that describes well the experimental data.

Phenomenological Meissner kernel Kxx in TI/SC heterostructure
We assume exponential decay of the magnetic field inside the Nb layer, from both interfaces (TI/SC and SC/vacuum or substrate). The Meissner kernel is given by where λSC is the London penetration depth of the SC layer. This leads to the following solutions for the field and the vector potential, The kernel inside the TI layer can be written as Kxx(z) = −a 2 e −2bz , z > 0 (S34) where a ≡ 1/λTI and b ≡ 1/z0 are fit parameters (see Eq. (2) of the main text). This leads to the following solutions for the field and the vector potential BTI(z) = ae −bz d1J1 a b e −bz + d2Y1 a b e −bz , z > 0, (S35) where J0,1 and Y0,1 are 0 th or 1 st order Bessel functions of the first and second kind, respectively. The values of c1, c2, d1 and d2 are determined by the boundary conditions. The results of the fit of the experimental field profile using the Bessel-function based solution for the vector potential are presented in Fig. 1 of the main text.

FINITE ELEMENT MODELLING OF TI/SC HETEROSTRUCTURE
Here we model the stray fields from a slab of thin superconducting layer to rule out such effects as a source of the observed magnetic field depth dependence in the TI/SC heterostructure. Although trivial, we present these calculations for completeness.