Odd-frequency superconductivity and Meissner effect in the doped topological insulator Bi$_2$Se$_3$

Doped Bi$_2$Se$_3$ is proposed to be a nematic superconductor originating from unusual inter-orbital pairing. We classify all induced superconducting pair correlations in Bi$_2$Se$_3$ and discover that intra-orbital odd-frequency pairing dominates over a significant range of frequencies. Moreover, we explore the contributions of even- and odd-frequency pairing to the Meissner effect, including separating intra- and inter-band processes in the response function. Contrary to expectations, and due to inter-band contributions, we find a diamagnetic odd-frequency Meissner effect.

Superconductivity in Bi 2 Se 3 has been a focal point of experimental and theoretical research in the past few years. Bi 2 Se 3 had already gained prominence as a strong topological insulator [1,2], when it was discovered that electron doping by intercalation of Cu (and later Nb and Sr) leads to the appearance of superconductivity [3]. Due to the strong spin-orbit coupling, the pairing was already early on proposed to be of odd parity, making doped Bi 2 Se 3 an intrinsic topological superconductor [4]. More recently, a series of experiments have additionally discovered a surprising breaking of the rotation symmetry of Bi 2 Se 3 in the superconducting state [5][6][7][8][9]. This so-called nematic superconducting state has been found to only be consistent with an exotic inter-orbital order parameter, meaning the Cooper pairs consist of two electrons from different orbitals in the Bi 2 Se 3 low-energy structure [10], in contrast to more standard intra-orbital pairing.
In the presence of an additional electronic degree of freedom, like the orbital index in Bi 2 Se 3 , the symmetry classification of superconducting pairing is extended [11,12]. In particular, in the presence of spin, spatial parity, orbital (or band, layer, valley etc.), and time/frequency quantum numbers a total eight classes are possible, of which four are odd in frequency [11,12]. Odd-frequency superconductivity has the electrons paired at different times, with an odd relative time difference. It was originally envisaged for 3 He [13] and has played a major role in the understanding of proximity-induced superconductivity in heterostructures [14][15][16][17][18][19][20][21]. More recently, oddfrequency pair correlations have also been shown to be common in bulk systems in the presence of inter-orbital terms in the Hamiltonian [22][23][24][25][26].
As odd-frequency pairing joins electrons at different times, direct probes of this intrinsically dynamical state have proven to be challenging. One of the few easily accessible physical probes has been the Meissner effect. In contrast to the usual diamagnetic Meissner response of even-frequency pairing, which expels an external magnetic field from the superconductor, simple quasiclassical approaches have predicted a paramagnetic Meissner effect for the proximity-induced odd-frequency pairing in heterostructures [21,[27][28][29], which was re-cently also experimentally observed in superconductorferromagnet junctions [30]. A paramagnetic response for odd-frequency pairing has also been calculated for generic multiband Hamiltonians [31]. Because a paramagnetic Meissner effect means the superconductor attracts instead of repels magnetic fields, the superconducting state should become unstable [32]. It would therefore be very intriguing to discover odd-frequency pairing with a diamagnetic Meissner response.
In this letter, we classify all pair correlations induced in the inter-orbital nematic superconductor Bi 2 Se 3 and discover prominent intra-orbital pairing terms, which are odd in frequency, even dominating the even-frequency pairing for a wide range of frequencies. Moreover, we calculate the Meissner response of both the even-and odd-frequency pair correlations and analyze inter-and intra-band processes separately. Surprisingly, we obtain a diamagnetic Meissner effect for the odd-frequency components, due to a diamagnetic inter-band response. This is both important for the stability of the superconducting phase in Bi 2 Se 3 and opens a road for designing stable odd-frequency superconductors. Finally, we find a clear two-fold rotational symmetry in the Meissner response signaling the nematic state.
Model-The low-energy physics in the normal state of Bi 2 Se 3 is well-captured by a linear momentum model with two orbitals [4,33] where σ i and s i indicate the Pauli matrices in orbital and spin spaces, respectively, v and v z are the Fermi velocities of the electrons in-and out-of-plane, m hybridizes the different orbitals, and µ is the chemical potential. We also set = 1. The eigenvalues ofĤ 0 are given by forming a gapped three-dimensional (3D) Dirac dispersion, with the twofold degenerate valence and conduction bands separated by a 2m energy gap.
We introduce superconductivity through the pairing matrix∆ and write the full Hamiltonian in Nambu space , and ∆ = 0∆ ∆ † 0 , where we use... (...) to signal 4×4 (8×8) matrices. Out of the four possible k-independent pairing symmetries identified for doped Bi 2 Se 3 [4], the experimental discovery of nematic superconductivity [5,6] singles out an unconventional inter-orbital spin-triplet state. Specifically, this state transforms according to the 2D E u irreducible representation of the D 3 point group of Bi 2 Se 3 and has the form (∆ x ,∆ y ) = ∆(s 0 ⊗ iσ y , s z ⊗ iσ y ) [10]. Below the transition temperature T c , the order parameter can form any linear combination, parameterized bŷ ∆ = A x∆x + A y∆y with coefficients A x,y . The nematic state is formed by (A x , A y ) = (cos(θ), sin(θ)), where θ corresponds to the in-plane angle of the nematic director. With complex coefficients (A x , A y ) = 1 √ 2 (1, i) the order parameter instead describes a chiral superconductor.
For all numerical evaluations we set ∆ = 0.3 meV, which is similar to values measured in scanning tunneling experiments [34]. The other parameters are obtained from comparison to DFT calculations as m = −0.28 eV, v = 0.434 eV, and v z = −0.248 eV [33]. Superconductivity in Bi 2 Se 3 is observed for electron doping [35] and we therefore study 0.27 eV < µ < 0.4 eV, ranging from the Fermi level close to the bottom of the bulk conduction band to the doped metallic regime.
Odd-frequency pairing-The presence of odd-frequency pairing is unveiled by a complete classification of the superconducting pair correlations, or amplitudes, present in the anomalous Green's functionF , which is directly obtained from the Matsubara Green's functioň Here,Ĝ andĜ are the normal Green's functions in particle and hole spaces, respectively. To obtain compact analytical expressions for all types of pair correlations, we treat, in a first step, the superconducting order parameter ∆ as a small quantity relative to the other energy scales and expand the Green's function to first order in ∆. The anomalous part then reduces toF (1) =Ĝ 0∆Ĝ0 , whereĜ 0 is the Green's functions of the bare Hamiltonian H 0 . An overview of all pairing terms inF (1) for any choice of (A x , A y ) and their classification according to their spin, parity, orbital, and frequency symmetries is shown in Table I Table I). One of the eight different terms constitutes spin-triplet, s-wave, inter-orbital, even-frequency pairing, i.e. the same symmetry as the order parameter (blue in Table I). Apart from this, there is only one other s-wave term, i.e. k-independent: a spin-triplet, intraorbital, odd-frequency pairing (red). Thus, even though the superconducting order parameter has unconventional ↑↑, ↓↓ pz even-inter even 2(A · k)vµ∆ ↑↓ + ↓↑ px,y even-inter even  (1) for generic (Ax, Ay) according to their spin, parity, orbital, and frequency symmetries. Blue term has the same symmetries as the order parameter, while the red highlights the s-wave intra-orbital odd-frequency pairing. Here inter -orbital symmetry, there also exists odd-frequency spin-triplet, but otherwise conventional intra-orbital swave pairing, explicitly induced by the hybridization between orbitals m. We also obtain another odd-frequency term, with spin-singlet, p-wave, even inter-orbital symmetry. Both of these odd-frequency (and all evenfrequency) pair correlations are proportional to A x and A y in such a way, that they are present in both the nematic and chiral superconducting phases. The induced odd-frequency pair amplitudes exceed the even-frequency amplitudes with the same spatial parity for a wide range of frequencies. Comparing the even-and odd-frequency s-wave pair amplitudes inF (1) , the oddfrequency pairing dominates for m + µ < |ω| < −m + µ, when the system is doped in the bulk regime, i.e. µ > |m|. The odd-frequency p-wave pairing becomes larger than the two even-frequency in-plane p-wave pair amplitudes for |ω| > m, and |ω| > µ, respectively. These analytical findings are corroborated by a numerical analysis to infinite order in ∆. We compare pair amplitudes as function of frequency by integrating the absolute value |F (iω)| over k making the usual replacement iω → ω ± i0 + . To study p-wave symmetries we multiply by the corresponding form factors before the integration. The window of frequencies in which the odd-frequency pairing dominates is clearly visible in Fig. 1 and agrees well with the analytical prediction from first order pertubation theory (shaded areas), which also further strengthens the confidence in our pertubative analysis.
Meissner effect-Having established the presence of significant odd-frequency pairing in doped Bi 2 Se 3 , we turn to its influence on the Meissner effect. Tradition- ally, odd-frequency pairing has been assumed to manifest itself in a paramagnetic Meissner response, in contrast to the traditional diamagnetic Meissner effect of evenfrequency superconductivity [21,[27][28][29][30][31]. The Meissner effect is the response of a superconductor to an external magnetic field. Within linear response theory the current response j µ (q, ω e ) to an external vector potential A µ (q, ω e ) is governed by the current-current response function through j µ (q, ω e ) = −K µν (q, ω e )A ν (q, ω e ). Here, q and ω e are the wave vector and angular frequency of the external vector potential, and µ and ν spatial indices x, y, z. The Meissner response is obtained in the limit of a static, uniform magnetic field ω e → 0, q → 0 (in that order) [37,38].
We introduce the vector potential in our Hamiltonian Eq. (1) by replacing k → k− A (setting e = 1) and calculate the current density operators by taking a variational derivative with respect to A µ [39]. Usually, the current density consists of a paramagnetic (j P ∝ A 0 µ ) and a diamagnetic part (j D ∝ A µ ). However, due to the linear spectrum of the HamiltonianĤ 0 , the diamagnetic current vanishes and the current-current response function reduces to K µν (q, ω e ) = j P µ (q, ω e )j P ν (−q, ω e ) [40], see SM for full derivation [36]. We express this expectation value to infinite order in ∆ with the help of the Green's function as where Tr e is a trace over the particle part of the matrix, only, T is the temperature andĵ P µ = −∂Ĥ 0 (k − A)/∂A µ andĵ P µ = −ĵ P * µ are the paramagnetic current operators in particle and hole space, respectively. We have also suppressed the momentum and frequency dependence for legibility.
The even-and odd-frequency pair correlations are naturally contained only in the anomalous Green's function F =F e +F o , such that we can work out their influence on the Meissner response by focusing only on the last term in Eq. (3) [31]. This term, K (S) µν , captures the superconducting contribution to the Meissner effect, which we also checked by numerically subtracting the normal state contribution of K µν from Eq. (3) [40], see SM [36]. It can be expressed as a sum of even-and odd-frequency contributions K (S)  [31]. To further simplify the results, we use the bands obtained from diagonalizingȞ to split the response function K (S) µν into intra-and inter-band contributions. The intra-band contributions are dominated by quasiparticles excited just above the superconducting gap, while the inter-band contribution, on the other hand, becomes enhanced when two bands approach each other in the vicinity of, but not necessarily at, the Fermi surface. After carrying out the Matsubara summation analytically, we integrate numerically over k at T = 2 K, see SM [36] for further details.
The resulting even-and odd-frequency contributions to the Meissner effect in the nematic state are presented in Figs. 2 (a,b) as a function of chemical potential. In total, the Meissner response is diamagnetic and dominated by the even-frequency pairing, especially for µ ≫ m. Contrary to expectation, however, the odd-frequency contribution is not para-, but diamagnetic for most field directions and parameters. This surprising response can be understood by separately studying the intra-and interband processes, which we display for the K (S) xx component at doping levels around the onset of the conduction band in Figs. 2 (c,d). For intra-band processes, the even-frequency pairing gives rise to a diamagnetic and the odd-frequency pairing to a paramagnetic Meissner response, as expected. The inter-band processes, however, contribute with the opposite sign. The total response can then be para-or diamagnetic depending on the balance between intra-and inter-band processes. For the odd-frequency pairing in doped Bi 2 Se 3 we find intraand inter-band processes of the same order of magnitude and also notably larger than the corresponding evenfrequency processes below the onset of the conduction band. The odd-frequency contributions thus nearly cancel in an intricate way, and for the realistic parameters used in Fig. 2 (c), in total yield a reduced but still diamagnetic odd-frequency Meissner response. However, depending on the ratio of intra-to inter-band contributions, the Meissner response from the odd-frequency pairing can also be paramagnetic, such as for the K (S) yy component in Fig. 2 (b), or even strongly diamagnetic, if the inter-band processes become more dominant. The evenfrequency pairing, on the other hand, yields a dominating diamagnetic Meissner response, because the intra-band contribution is much larger than the inter-band one for µ > |m|, as clearly seen in Fig. 2 (c).
Nematicity-Finally, we discuss the nematicity of the Meissner response, which is already apparent in Fig. 2 zz . A different response in the z-direction is expected from the different Fermi velocities, |v| > |v z |. The in-plane Fermi surface is, however, isotropic, such that the observed nematicity in the inplane Meissner response, with K The nematicity also appears in the response for different nematic angles θ as we display in Fig. 3, which is equivalent to the Meissner response for a fixed nematic angle but with the magnetic field rotated in the x-y plane. The K Before concluding, we note that the response presented in Fig. 3 only contains the superconducting contribution to the Meissner effect. The total Meissner response K µν requires also including the first term in Eq. (3), as well as the contribution from the diamagnetic current operator for corrections to the low-energy spectrum beyond linear order in k. However, adding any of these neglected terms will only overlay additional two-or six-fold symmetric quantities on top of the very pronounced two-fold rotation symmetry in Fig. 3. Thus, we still expect a clear signature of nematic superconductivity in the Meissner effect, experimentally detectable e.g. in the London penetration depth.
Conclusions -In conclusion, we identify large and dominating odd-frequency intra-orbital s-wave pair correlations in nematic superconducting state in doped Bi 2 Se 3 , even though the order parameter itself constitutes only even-frequency inter-orbital s-wave pairing. We further discover strong odd-frequency intra-and inter-band contributions to the Meissner response of opposite signs that cancel in an intricate way to give rise to a small, but unexpectedly diamagnetic Meissner response. This both stabilizes the superconducting state in Bi 2 Se 3 and opens the possibility to find or engineer other odd-frequency superconductors where enhanced interband contributions give a stable diamagnetic Meissner response. Finally, we calculate a strong two-fold symmetric Meissner response as a function of the direction of an in-plane magnetic field, giving a measurable signature of nematic superconductivity in Bi 2 Se 3 .
We are grateful to D. Kuzmanovski  In this supplementary material we first present the complete form of the anomalous Green's function. Then we explain the derivation of the Meissner kernel in more detail and finally describe how we perform the Matsubara summation with focus on the splitting into the intra-and inter-band contributions.

I. ANOMALOUS GREEN'S FUNCTION TO INFINITE ORDER IN ∆
In the main text, we perform a perturbation expansion of the anomalous Green's function up to first order in the order parameter ∆, with all induced pair correlations are presented in Table I of the main text. Here we report the anomalous Green's function obtained completely non-perturbatively from Eq. 2 in the main text for the nematic state. We express the anomalous Green's function in the basis (1 ↑, 1 ↓, 2 ↑, 2 ↓), where 1, 2 indicate the different p 1z , p 2z -orbitals that form the basis of the Hamiltonian of the normal-state doped Bi 2 Se 3 (Eq. 1 in the main text) and ↑, ↓ show different spins [1]. We also separate the odd-and even-frequency parts F = F o + F e , which take the form and with the shorthand notation k x,y ≡ vk x,y , k z ≡ v z k z , γ ± = ((m ± ik z ) 2 − µ 2 − ω 2 − ∆ 2 ), A ± = A x ± iA y , A × k = A x k y − A y k x , and A·k = A x k x + A y k y . The denominator terms are given by D ± = ((iω) 2 −ξ 2 ± ), where ξ ± = ǫ 2 + ∆ 2 + µ 2 ± ∆ 2 (m 2 + (A × k) 2 ) + ǫ 2 µ 2 are the eigenvalues of the full HamiltonianȞ. Here ǫ = |ǫ 0 ± + µ| where, ǫ 0 ± are the eigenvalues of the normal Hamiltonian, ǫ 0 ± = ± m 2 + k 2 x + k 2 y + k 2 z − µ. It should be noted that the denominator is even in frequency and does not influence the symmetry classification.
Analyzing the non-perturbative results in Eqs. (1)-(2), we find that they are slightly more complicated but fully consistent with the results in the main text calculated perturbatively to first order in ∆. The first-order results are obtained by replacing γ ± → ((m ± ik z ) 2 − µ 2 − ω 2 ) and ξ ± → ǫ 0 ± , which are notably also valid for the chiral state, not just the real combinations of A x,y which generate the nematic state.

II. DETAILED DERIVATION OF THE MEISSNER EFFECT
In this section we present a more detailed derivation of the Meissner effect for doped Bi 2 Se 3 . The Meissner effect describes the response of a superconductor to an external magnetic field. It can phenomenologically be derived from the London equation, which relates the superconducting current density j to an external vector potential A by j = − n s e 2 m A., where n s is the superconducting density, and m and e are the electron mass and charge, respectively. Together with the definition of the vector potential ∇ × A = B and Ampere's law ∇ × B = µ 0 j, it can be used to derive the differential equation for the magnetic field inside a superconductor: ∇ 2 B = 1 λ 2 B. This differential equation yields an exponentially decaying solution and thus the magnetic field is expelled in the interior of a superconductor. Here, λ = m n s e 2 µ 0 is the London penetration depth, the characteristic length scale at which the magnetic field is suppressed.
In the London picture, the suppression of the magnetic field occurs because the coefficient n s e 2 m is assumed to always be positive. Within a microscopic theory the phenomenological London equation is replaced by the more general relation presented in the main text, where q and ω e are the wave vector and angular frequency of the external vector potential, and µ and ν represent the spatial indices x, y, z. The current-current correlation function K µν can in principle take any sign, such that the magnetic field can either be suppressed (for K > 0) or increased (K < 0) in the superconductor. The traditional Meissner effect, suppressing the magnetic field, is then called the diamagnetic Meissner effect in contrast to a paramagnetic Meissner effect. Within linear response theory, we can treat the vector potential as a perturbation such that K µν (q, ω e ) takes the general form [2][3][4] where all expectation values are taken with respect to the unperturbed Hamiltonian and j P and j D are the paramagnetic and diamagnetic current operators, respectively. To obtain the current operators, we introduce the vector potential A in the Hamiltonian through the substitution k → k − A, expand to second order in A, and differentiate with respect to A [4,5]. The resulting current can be split into the sum j P µ + j D µ,ν A ν . The Hamiltonian for Bi 2 Se 3 presented in Eq. (1) in the main text is linear in the crystal momentum k. Thus, when we introduce the vector potential through the substitution k → k − A, there is no quadratic term in A, and the hence j D µ automatically vanishes. As a result, the current-current response function in Eq. (5) reduces to the expectation value K µν (q, ω e ) = j P µ (q, ω e ) j P ν (−q, ω e ) . The Meissner effect is the response to a static, uniform external vector potential, corresponding to the limit ω e → 0, q → 0. Then, expressing the expectation value of the second quantized current operators in the expression for K µν with the help of the Green's functions, we arrive at Eq. (3) in the main text: HereǦ is the full 8 × 8 Green's function,Ĝ andF are the 4 × 4 normal and anomalous Green's functions, respectively, T is the temperature, and Tr e is a trace over the particle part only of the full 8 × 8 matrix. Moreover,ĵ P are the first quantized paramagnetic current operators, which we easily obtained from such thatĵ x = −ĵ * x = vs y ⊗ σ z ,ĵ y = −ĵ * y = −vs x ⊗ σ z , andĵ z = −ĵ * z = v z σ y . The vanishing diamagnetic current due to the linear dispersion also means that the normal Green's function contribution T r[Gĵ µ Gĵ ν ] technically contributes with a finite Meissner response even in the limit ∆ → 0. Usually, this contribution is cancelled by the diamagnetic current from higher order terms in k in the full band structure [6], but such terms are missing here. Because we are primarily interested in the contributions due to the odd-and even-frequency superconducting pair correlations in the anomalous Green's function, we can safely consider only the superconducting contribution to the Meissner response [7], given by the terms involving the anomalous Green's functions, Moreover, the normal Green's function contribution T r[Gĵ µ Gĵ ν ] not only attains a non-zero value in the absence of superconductivity, it also diverges at high energies. It is possible to regularize it by subtracting the normal-state response T r[Ĝ 0ĵµĜ0ĵν ] from Eq. (5), which provides another way of guaranteeing a vanishing Meissner response when ∆ = 0 [6,8]. We have numerically performed this regularization and find that, T r[Ĝĵ µĜĵν ] − T r[Ĝ 0ĵµĜ0ĵν ] vanishes indeed for ∆ = 0, but even for finite and realistic values of ∆. Hence, K (S ) µν in Eq. (7) represents the full superconducting contribution to the Meissner response, and is the quantity used in the calculations reported in the main text.
Finally, in order to analyze the Meissner kernel we split it into a summation of the even-and odd-frequency pairing contributions: