Dia- and paramagnetic Meissner effect from odd-frequency pairing in multi-orbital superconductors

The Meissner effect is one of the defining properties of superconductivity, with a conventional superconductor completely repelling an external magnetic field. In contrast to this diamagnetic behavior, odd-frequency superconducting pairing has often been seen to produce a paramagnetic Meissner effect, which instead makes the superconductor unstable due to the attraction of magnetic field. In this work we study how both even- and odd-frequency superconducting pairing contributes to the Meissner effect in a generic two-orbital superconductor with a tunable odd-frequency pairing component. By dividing the contributions to the Meissner effect into intra- and inter-band processes, we find that the odd-frequency pairing actually generates both dia- and paramagnetic Meissner responses, determined by the normal-state band structure. More specifically, for materials with two electron-like (hole-like) low-energy bands we find that the odd-frequency inter-band contribution is paramagnetic but nearly canceled by a diamagnetic odd-frequency intra-band contribution. Combined with a diamagnetic even-frequency contribution, such superconductors thus always display a large diamagnetic Meissner response to an external magnetic field, even in the presence of large odd-frequency pairing. For materials with an inverted, or topological, band structure, we find the odd-frequency inter-band contribution to instead be diamagnetic and even the dominating contribution to the Meissner effect in the near-metallic regime. Taken together, our results show that odd-frequency pairing in multi-orbital superconductors does not generate a destabilizing paramagnetic Meissner effect and can even generate a diamagnetic response in topological materials.


I. INTRODUCTION
Odd-frequency superconductivity is an unusual superconducting state where the two electrons forming a Cooper pair join each other at different times, such that the resulting pair amplitude is odd under the interchange of the internal time coordinate, or equivalently, has an odd frequency dependence [1][2][3][4]. The oddness under time interchange has direct consequences for the other symmetry properties of the superconducting state, since the fermionic nature of the pair amplitude requires it to be odd under a full interchange of all the electron quantum numbers. For example, the conventional superconducting spin-singlet s-wave state can be converted into a spintriplet s-wave state if the frequency dependence becomes odd. This particular situation has been shown to arise when a conventional superconductor is in proximity with a ferromagnet, where the spin-triplet pairing allows for a long-range proximity effect into the ferromagnet [4][5][6].
More recently, odd-frequency superconductivity was also shown to occur in materials with an additional lowenergy degree of freedom [7], such as multi-orbital or band [8,9], double dot [10,11], or bilayer structures [13]. In these systems the orbital (or similar) degree of freedom adds an additional index to the pair amplitude, such that it is, e.g., possible to have odd-frequency pairing in the form of spin-singlet s-wave odd inter-orbital pairing. Notably, this odd-frequency pairing easily occurs in the bulk and does not require any spatial inhomogeneity in form of junctions or similar. The condition to generate odd-frequency pairing in a multi-orbital superconductor [7,8] is in fact easily fulfilled as it only requires finite inter-orbital (single particle) hybridization and that the strength of the superconducting pairing is different in different orbitals, also known as a pairing asymmetry or orbital selectivity [19]. As such, many known multiband superconductors have recently been shown to host odd-frequency pairing, including doped topological insulators [9,14], Sr 2 RuO 4 [15], UPt 3 [16], and superconducting Weyl semimetals [17,24]. Based on these results, even the iron-based superconductors with orbital selective pairing [19] likely host odd-frequency pairing.
As odd-frequency pairing is intrinsically time dependent it has so far been hard to detect directly, becoming essentially a hidden dynamic order. However, for specific systems signatures of odd-frequency pairing have been found in experimentally accessible quantities, such as the existence of finite density of states at zero energy [21], Kerr effect [15,16], or Josephson current in junctions where even-frequency superconductivity becomes forbidden [17,18,22,23]. Most prominent of the odd-frequency signatures is probably however the prediction of a paramagnetic Meissner effect [25][26][27][28][29][30][31][32].
The Meissner effect is the response of a superconductor to a weak magnetic field, with a conventional superconductor always completely repelling the magnetic field, known as a diamagnetic Meissner effect [33]. In contrast to this conventional diamagnetic response, odd-frequency pairing has instead in many circumstances been predicted to cause a paramagnetic Meissner effect, where the superconductor instead attracts a magnetic field, which then also easily destabilizes the whole superconducting state [25][26][27][28][29][30][31][32]. This paramagnetic response has recently been experimentally confirmed in regions where odd-frequency superconductivity is proximity-induced in various heterostructures [29,32]. In these heterostructures the paramagnetic response is not a problem, since superconductivity is originating from a conventional bulk supercon-ductor that still has a diamagnetic Meissner effect and is thus stable in the presence of a weak magnetic field. However, finding a paramagnetic Meissner effect in the bulk of a superconductor directly raises issues of stability of the superconducting phase. In fact, the existence of a paramagnetic Meissner effect has been used as a key argument against the existence of an intrinsic bulk oddfrequency superconducting state [30,34,35].
The issue with a possible paramagnetic Meissner effect from odd-frequency pairing causes a clear conundrum for odd-frequency superconductivity in multi-orbital superconductors, as the odd-frequency pairing is here clearly a bulk effect and therefore directly raises the question of how stable these superconductors actually are. In this work we tackle this problem by investigating the relation between odd-and even-frequency pairing and Meissner effect in a generic two-orbital bulk superconductor, where the amount of odd-frequency pairing is directly tunable by the asymmetry in the superconducting paring between the two orbitals. By dividing up the contributions to the Meissner effect from the even-and odd-frequency pairing and also from intra-and inter-band processes, we are able to derive simple analytical expressions that show that odd-frequency pairing generates both diamagnetic (positive) and paramagnetic (negative) Meissner contributions. In particular, the intra-and inter-band oddfrequency contributions always have different signs, with the signs determined by the normal-state band structure.
More specifically our results show that for a material with two electron-like bands (or hole-like), see Fig. 1(a), we find that the odd-frequency inter-band Meissner contribution is paramagnetic, but it is always compensated by an almost equally large but diamagnetic intra-band contribution. Adding to that an overall diamagnetic contribution from the even-frequency pairing, the total Meissner effect of the superconductor is always diamagnetic, even in the limit of large odd-frequency pairing. On the other hand, for a material with an inverted band structure, i.e. one electron-like and one hole-like band as found in topological materials, see Fig. 1(b), oddfrequency pairing instead generates a diamagnetic interband contribution to the Meissner effect that is always larger than the intra-band paramagnetic contribution. Very interestingly, in the near-metallic regime, the oddfrequency inter-band contribution becomes dominating and thus the total diamagnetic Meissner effect of the superconductor is even driven by the odd-frequency pairing. These results show that odd-frequency pairing in multi-orbital superconductors does not produce a destabilizing paramagnetic Meissner effect but can even generate a diamagnetic response in topological materials.
The organization of the rest of the article is as follows. In section II we first present the Hamiltonian for a generic two-orbital superconductor and then derive its superconducting pair amplitudes in subsection II A and the general theory and calculational details for the Meissner effect in subsection II B. We then present our results in section III, both analytical and numerical, divided up into two different cases: when the resulting energy bands have the same signs of their curvatures, both being electron-like in subsection III A, and when they are inverted, or topological, with different signs of their curvatures in subsection III B). Finally, we summarize our results in section IV. We also provide some more details of the calculations in Appendix A and B.

II. THEORY
A. Generic two-orbital superconductor To investigate odd-frequency superconductivity and Meissner effect in multi-orbital superconductors, we assume a bulk material with a two generic low-energy orbitals in the normal state. This is the simplest model that leads to the presence of bulk odd-frequency superconducting pairing. The normal state is described in the orbital basis as h(k) = ξ + τ 0 + ξ − τ z + ξ 12 τ x , where ξ 1(2) = ξ + ± ξ − and ξ 12 represent the intra-and interorbital kinetic energy dispersions, respectively. Here, τ represents the Pauli matrices in the orbital basis and we assume the Hamiltonian to be spin-degenerate. Diagonalizing this normal part results in the two energy bands: .
To make our work general, yet simple enough to allow a partial analytical treatment, we assume the dispersion relations in two orbitals to be of the form ξ 1(2) = −µ + t 1(2) k 2 − (−1) 1(2) δµ/2, with chemical potential µ, orbital energy difference δµ, electron wavevector k, and effective curvature (the inverse of the effective mass) t 1,2 . For the numerical results we limit ourselves to a twodimensional wave vector to keep the momentum integrations feasible, but we expect no significant changes for three dimensions. Moreover, we mainly assume a kindependent hybridization ξ 12 = t 12 between the two orbitals, but we report complementary results for kdependent hybridization ξ 12 (k) in Appendix A. Throughout this work, we set the energy scale through t 1 = 1 and in order to investigate the effects of the band structure, we separately treat two cases, t 2 > 0 and t 2 < 0. In Fig. 1 we plot the band structure given by Eq. (1) for the two different cases, where in (a) t 2 = 0.5 and (b) t 2 = −0.5. Thus, in (a) we have two parabolic bands with the similar electron-like dispersions but different curvatures with their energy minima at ± t 2 12 + (δµ/2) 2 . In contrast, in (b) we get a so-called inverted band structure, consisting originally of overlapping electron-and hole-like bands, which hybridize and form a gapped band structure, such that the Fermi level never crosses both bands. This second model has a band structure closely resembling that of topological insulators and other non-trivial topological materials. We include superconductivity in the system by introducing a conventional spin-singlet intra-orbital pairing in each orbital, written as∆ = (δ + τ 0 + δ − τ z ) ⊗ iσ y , with σ being the Pauli matrices in spin space. Here the spin-singlet order parameter in each orbital is δ 1(2) = δ + + (−)δ − , which is the simplest possible implementation of superconductivity in a two-orbital model. As we will explicitly see later, odd-frequency pairing is always present in this model for a finite δ − , as long as also ξ 12 is finite. The other alternative route to generate oddfrequency pairing is to assume a finite inter-orbital pairing δ 12 [7,15], but such term can always be described within a change of basis of our current model, see Appendix B, and thus it cannot qualitatively change our results. Finally, writing the full Hamiltonian in the orbital basis ψ † = (c † 1↑ c † 2↑ c 1↓ c 2↓ ) we arrive at where γs are Pauli matrices in the particle-hole basis. The energy eigenvalues of the full Hamiltonian, H are equal to

Superconducting pair amplitude
To understand superconducting pairing in our generic two-orbital superconductor, Eq. (2), we extract the Green's function of the system: Here, G andḠ represent the electron and hole propagators, respectively, while F andF † = F are the anomalous Green's functions, or equivalently, the superconducting pair amplitude. The anomalous Green's function can further be decomposed into its even-and odd-frequency components F = F e + F o . For the Hamiltonian given in Eq.
The even-frequency pair amplitude in Eq. (5) include both intra-orbital (diagonal components) and interorbital (off-diagonal components) components, where the latter is directly proportional to the inter-orbital hybridization, ξ 12 . The odd-frequency pair amplitude Eq. (6) has only inter-orbital components and is a consequence of the existence of superconducting orbital selectivity, i.e. a finite δ − , and finite inter-orbital hybridization, ξ 12 , i.e. there must exist an asymmetry in the superconducting pairing between the two orbitals, or an orbital selectivity, together with a finite hybridization between the two orbitals for odd-frequency pairs to exist. In addition to generating odd-frequency pairs, the orbital selectivity term δ − causes the superconducting part of the Hamiltonian to not commute with the normal part. This incompatibility between the normal and superconducting parts has recently also been used to define the concept of superconducting fitness [20], which has been shown to be directly linked to the existence of odd-frequency pairs [7]. We can here explicitly verify that all pair amplitudes satisfy the fermionic nature of superconductivity: the inter-orbital pairs are even (odd) under the interchange of orbital index when they have an even-(odd-)frequency dependence, as required for spin-singlet pairs with no k-dependence (s-wave symmetry). Above we explicitly worked in the orbital basis, for the form of the Hamiltonian in the band basis, where the kinetic energy is diagonal, see Appendix B.

B. Meissner effect
The Meissner effect is the response of a superconductor to an externally applied magnetic field. Within linear response theory, the current response function j to the vector potential of the magnetic field A can be obtained via j µ (q, ω) = −K µν (q, ω)A ν (q, ω). Here, K is the currentcurrent correlation function and ν, µ are the directions of the applied vector potential and current response, respectively, while q, ω are the wavevector and frequency of the response function. The Meissner response is obtained for a static and uniform magnetic field, thus taking the limits q → 0 and ω → 0. To calculate K we need the current operator matrix J = J p µ + J d µν A ν , divided into its paramagnetic and diamagnetic part and which we find in the standard way by first introducing the vector potential A into the Hamiltonian through the substitution k → k−A and then taking the derivative with respect to A ν . The Meissner response can then be written as [9,36] where G is the Green's function and Tr e represents the trace over the electron part of the matrix. We here work in natural units, such that we set = c = e = m = 1. Moreover, since our system is isotropic in space, we can without loss of generality focus on a specific direction, say the x-direction, of the Meissner response function and thus drop the µ, ν indices. Finally, using the fact that the current operator J p(d) in the particle-hole basis takes the , the Meissner response in Eq. (7), can be simplified as Since the focus of this study is to obtain the contributions of the odd-and even-frequency pairing to the Meissner effect, we only need to focus on the second term of the above equation: because only this term depends on the pair amplitudes F,F . Moreover, the first and the last terms usually cancel each other in conventional and also multi-band superconductors [9,37]. Our interest is thus focused on the Meissner Kernel of Eq. (9): K = Tr{FJ pF J p }. This term can always be decomposed into K = K e + K o , where K e = Tr{F eJ pF e J p } is entirely due to the evenfrequency pairing, while K o = Tr{F oJ pF o J p } comes exclusively from the odd-frequency pairing. Technically there also exists terms of the form F oJ pF e J p in K, but, since they are always odd in frequency, they automatically cancel during the final frequency summation in Eq. (8).

Kernel decomposition into intra-and inter-band processes
To better understand the physics of the Meissner response, we decompose its Kernel K into parts coming from inter-and intra-band processes, respectively. Using the fact that the Meissner Kernel can generically be written as K = (aω 4 + 2bω 2 + c)/D 2 , we can form the following decomposition: Here K ± and K 12 are frequency-independent coefficients for processes involving only one individual band (i.e. intra-band processes) and both bands (i.e. inter-band processes), respectively. These coefficients are given by We here note that the intra-band term is the only term appearing in a single-orbital superconductor, or for a two-orbital superconductor in the trivial ξ 12 = 0 limit, while the inter-band term only appears when ξ 12 is finite.
We also emphasize here that these are processes between different bands, i.e. the eigenstates of the Hamiltonian Eq. (2), and not processes between the different orbitals.

Frequency summation
In order to analytically perform the summation over frequency in the Meissner Kernel in Eq. (9), we use analytical continuation from real frequency to Matsubara frequency, ω → iω n . Then, by applying Matsubara frequency summation, the Meissner response reads where n(ξ) is the Fermi-Dirac distribution function, n (ξ) its derivative with respect to energy, and C(ξ) = (n(−ξ) − n(ξ))/2ξ. Here we note that both the C(ε) and functions are always positive, while n (ε) is always negative. Moreover, assuming we always have a gapped superconductor, as is usually the case for the bulk of s-wave states, the n (ε) term is negligible, and hence it is the signs of K ± and K 12 that become the decisive factors for determining the sign of the Meissner effect.
We can further simplify the expression for the Meissner effect by noting that, if we have a k-independent orbital hybridization ξ 12 as we assume here in the main text, the current operator is diagonal in the orbital basis: J p = j 1 0 0 j 2 and similarly forJ p = −J p (very results for a k-dependent hybridization are reported in Appendix A). Moreover, considering that the anomalous Green's function takes the form can be written where the last equality is just iterating the statement from before used to parametrize the Meissner Kernel.
Here, the first two terms in the first expression are always positive, while the last term changes sign based on the sign of j 1 j 2 and f 12 f 21 . Since j 1(2) = −∂ξ 1(2) (k − A)/∂A = 2t 1(2) k, the sign of the product j 1 j 2 is equivalent to the sign of t 2 , since we set t 1 = 1. It is thus natural when studying the Meissner effect to distinguish between the two cases: t 2 > 0 giving j 1 j 2 > 0, which has two bands with the same curvature, and t 2 < 0 with j 1 j 2 < 0, which has an inverted band structure. This is particularly important for the case of the odd-frequency pairing contribution to the Meissner effect as this pairing contains only the interorbital terms f 12 , f 21 .

III. RESULTS
Having established the underlying theory we can now focus on calculating the Meissner effect in generic twoorbital superconductors described by Eq. (2). Thanks to the derivations in the preceding section we can proceed analytically to a large degree. To calculate the Meissner Kernel in terms of its intra-and inter-band contributions in Eq. (10), we need the coefficients a, b, c in Eq. (11). For a k-independent orbital hybridization ξ 12 these are directly accessible using Eq. (13). For the contribution to the Meissner Kernel from the even-frequency pairing, we use Eq. (5) for the anomalous Green's function components and arrive at the coefficients: These, plugged in Eq. (11), directly give the evenfrequency Meissner Kernel contributions K e ± and K e 12 . For the contributions to the Meissner Kernel from the odd-frequency pairing we instead use Eq. (6) and find that a o = c o = 0, which means we can straightforwardly simplify the expressions to arrive at These equations show that for odd-frequency pairing, the sign of the intra-and inter-band Meissner Kernels, and thus their contributions to the Meissner effect, are opposite to each other and their sign only depends on the sign of j 1 j 2 . Thus, for a band structure with two bands with the same curvature, odd-frequency pairing always has a positive or diamagnetic Meissner effect from intra-band processes, while inter-band processes always give a negative or paramagnetic Meissner effect. On the other hand, for an inverted band structure, where then j 1 j 2 < 0, we get the opposite behavior: intra-band processes always generate a paramagnetic Meissner effect, while inter-band processes always give a diamagnetic effect.
To further analyze the results, especially establishing the relative sizes of the even-and odd-frequency and intra-and inter-band contributions to the Meissner effect, we need to perform the summation over reciprocal space in Eq. (12). Having established the importance of the sign of the product j 1 j 2 in Eq. (13), we divide our numerical results into two distinct cases: two electronlike bands where j 1 j 2 > 0 (or equivalently two hole-like bands) and an inverted band structure with j 1 j 2 < 0. In particular, below we choose t 2 = ±0.5 and also set δ + = 0.03 and ξ 12 to an k-independent constant t 12 . We then vary all other parameters; µ, δµ, and δ − , to arrive at a comprehensive picture of the Meissner effect in generic two-orbital superconductors. In terms of determining the impact of odd-frequency pairing on the Meissner effect, it is particularly important to vary δ − as the odd-frequency pair amplitude is directly proportional to this pairing asymmetry, or orbital selectivity.
A. Two electron-like bands, j1j2 > 0 We begin with a material consisting of two low-energy electron-like bands, here generically modeled with t 1 = 1 and t 2 = 0.5, such that j 1 j 2 > 0. In this case the oddfrequency pairing generates an intra-band Meissner effect that is always diamagnetic, while the inter-band part is always paramagnetic, as given by Eq. (15). However, as we in the plots see below, even-frequency pairing generate terms which can change sign depending on parameter values.
We start by investigating in Fig. 2 the even-and oddfrequency contributions to the Meissner effect divided up into intra-and inter-band processes and as a function of the energy difference between orbitals δµ, for different values of orbital selective superconducting order parameter δ − . As we fix δ + = 0.03, we here investigate three different regimes. For δ + > δ − we are in a so-called s ++ phase where the phase of the superconducting order parameter on each orbital is the same: δ 1,2 > 0, while for δ + < δ − we are in a s +− phase where δ 1 > 0 but δ 2 < 0. At the boundary, δ − = δ + = 0.03, we find δ 2 = 0, i.e. only orbital 1 is intrinsically superconducting. Moreover, we set µ = 0.3 which makes the Fermi level cross both bands, see Fig. 1(a), and thus both bands give large contributions to the Meissner effect.
In Figs. 2(a,b) we plot the contributions from the evenfrequency pairing to the Meissner response. The contribution from intra-band processes (a) is completely dominating and remain positive for the full range of parameters values, i.e. a diamagnetic Meissner response. The behavior is more or less the constant for both s ++ and s +− orbital pairing. The only exception is the special case δ + = δ − , which developes a notable decrease when δµ is increasing. We can understand this behavior by noting that, in this case there is no intrinsic supercon- ductivity in orbital 2 (δ 2 = 0) and increasing δµ shifts the band bottom of orbital 1 to higher energies, thus making its Fermi surface smaller. As a consequence, the Meissner intra-band contributions goes down with increasing δµ. While the even-frequency inter-band contribution in (b) is much smaller, we note that it interestingly has no fixed sign. Turning to the odd-frequency contributions, we see in Figs. 2(c,d) that the intra-(inter-)band part is dia-(para)magnetic, as also established by Eq. (15). Furthermore, we find that the inter-band contribution is generally larger than the intra-band contribution, in contrast to the dominating intra-band process for the even-frequency pairing. This leaves the total Meissner effect from odd-frequency pairing paramagnetic, in agreement with the traditional expectation for the oddfrequency response [25][26][27][28][29][30][31][32], but we note that the overall effect is partly canceled due to a relatively large diamagnetic intra-band contribution. Moreover, we find that the Meissner response from the odd-frequency pairing is significantly increasing with increasing orbital asymmetry δ − , in agreement with the effectively linear dependence on δ − for the odd-frequency pair amplitude.
In Fig. 3 we continue to analyze the Meissner response, focusing especially on how the different components compare to each other, with even-frequency intra-band (blue), even-frequency inter-band (red), oddfrequency intra-band (yellow), and odd-frequency interband (purple) together with the total response (green). Based on the distinct behaviors observed in Fig. 2, we here focus on the overall dependence on the chemical po- tential (a,b), and the superconducting asymmetry δ − , as they determine the Fermi surface size and magnitude of the odd-frequency pairing, respectively. For the Meissner effect as a function of chemical potential we choose to report values for both the s ++ (a) and s +− (b) regimes. We here see how all terms increase as a function of µ, which is due to the Fermi surface increasing with µ and thus creating a stronger superconducting state, which is then reflected in an increase in the total Meissner effect. The only minor exception to this increasing trend is the evenfrequency intra-band contribution around µ ∼ 0.22 for s +− pairing. Analyzing the superconducting order parameters in the band basis, we find that one of the intraband order parameters δ + − δ − cos(θ) becomes zero at this particular point, thus causing the drop in the intraband Meissner effect (for details and definition of θ, see Appendix B). This drop does not exist in (a) because there δ + > δ − and hence δ + − δ − cos(θ) never becomes zero. We also find in Fig. 3(a,b) that both the intra-and inter-band odd-frequency contributions become larger when δ − increases. This increase is further verified in Fig. 3(c,d), where we plot the different contributions to the Meissner effect as a function of the superconducting asymmetry δ − . Here, we choose two different values of the inter-orbital hybridization t 12 = 0.05 (c) and 0.1 (d), with all other parameters similar to Fig. 3(a,b). All Meissner contributions increase as a function of δ − , but notably the odd-frequency terms increases faster. This is to be expected as the odd-frequency pair amplitude is linearly proportional to δ − , see Eq. (6). Note, however, that we do not find that the odd-frequency Meissner contributions increase in the same manner for the orbital hybridization t 12 , despite this term also being necessary to generate odd-frequency pair amplitude F o , see Eq. (6). This is due to the term t 12 also changing the whole band structure and thus the pair amplitude (through the denominator D) and the Meissner response are influenced rather substantially. We find numerically a maximum odd-frequency Meissner response around t 12 ∼ 0.05, see Appendix A. With δ − being a much smaller energy scale, the same band structure effects are not seen for realistic δ − . Interestingly, we find that for the odd-frequency part, the inter-band contributions are larger the intraband ones for all parameters in Fig. 3. However, the difference is always small, so, although the inter-band contribution is paramagnetic, the diamagnetic intra-band contribution makes the total paramagnetic Meissner effect from the odd-frequency pairing quite small. As a consequence, odd-frequency pairing has only a minor effect on the Meissner effect in two-orbital superconductors with two electron-like (or hole-like) bands. In fact, we find that the total Meissner effect is essentially constant when δ − increases, despite this causing a strong increase in both odd-frequency pair amplitude and its influence on the Meissner effect. Hence we conclude that multiorbital superconductors are highly stable in a magnetic field, even when substantial bulk odd-frequency pairing is present. We emphasize that this is not due to small contributions to the Meissner effect from the odd-frequency pairing, but the fact that the intra-and inter-band oddfrequency processes nearly cancel each other.

B. Two inverted bands, j1j2 < 0
We next turn to the case of an inverted band structure, such that t 1 t 2 < 0, and equivalently j 1 j 2 < 0, as schematically illustrated in Fig. 1(b) and typical for topological insulators [38]. This band structure never has both bands crossing the Fermi level, and there also needs to be a finite doping to reach beyond the insulating state. In this case, Eq. (15) shows that the odd-frequency pairing always gives a paramagnetic Meissner effect for intra-band processes and a diamagnetic Meissner effect for the interband processes, i.e. opposite to the behavior in the previous section with a non-inverted band structure. This analytical result is in agreement with what has been found earlier for nematic inter-orbital spin-triplet superconductivity in doped topological insulators [9], indicating that these signs of the Meissner effect contributions are likely stable for a range of different symmetries for∆ even though we here concentrate on intra-orbital spin-singlet s-wave symmetry. Proceeding with numerical results, we next report similar plots for the inverted band structure as in for two electron-like bands in Figs. 2-3. Thus, in Fig. 4 we show the even-and odd-frequency (rows) and intra-and interband (columns) contributions to the Meissner effect as a function of δµ. We here initially fix the chemical potential µ = 0 such that both bands contribute similarly at δµ = 0, but we check the µ dependence in the next figure. Again we choose three increasing values of δ − creating δ 1 δ 2 > 0 (s ++ ), δ 2 = 0, and δ 1 δ 2 > 0 (s +− ) orbital pairing, respectively. Here we find that all Meissner contributions increases rapidly when δµ becomes more negative until around δµ ∼ −0.3. This is the value when the Fermi level crosses into the valence band, and thus this increase is due to more low-energy states being available. For even more negative values of δµ we see different trends. The even-frequency intra-band contribution (a) continues to increase, as expected for a metallic state with an increasing Fermi surface. The only exception is the case of δ − = 0.04, where, in a very limited regime, we surprisingly find a negative, paramagnetic contribution. This is the regime where superconductivity in one of the bands, δ + − δ − cos(θ), approaches zero, see Appendix B. The even-frequency inter-band contribution (b) gives both diamagnetic and paramagnetic responses, dependent on the choice of parameters, and is now notably larger than in the previous case with j 1 j 2 > 0 in Fig. 2. However, for the even-frequency contribution the intra-band processes are still generally dominating. When it comes to the odd-frequency contributions to the Meissner effect, we see how they form a peak structure quite similarly to that of Fig. 2, indicating a general behavior of the odd-frequency contributions. Most importantly, here the diamagnetic inter-band contribution is much larger than the paramagnetic intra-band contribution, which leaves the odd-frequency pairing giving a diamagnetic contribution to the total Meissner response. Thus, odd-frequency pairing actually helps stabilizing the superconducting state through the generation of an increased diamagnetic Meissner effect, in complete contrast to many other systems where an odd-frequency paramagnetic Meissner effect has often been discussed [25][26][27][28][29][30][31][32]. Finally, in Fig. 5 we plot in the same subpanels all contributions to the Meissner effect as a function of chemical potential µ (a,b) and superconducting orbital asymmetry δ − (c,d), i.e. similarly to Fig. 3 but now for the inverted band structure. For small values of the chemical potential the system is an insulator in the normal state and we find that all contributions to the Meissner effect are often of similar size. Interestingly, as we increase δ − by comparing (a,b), we get a notably higher contribution from the odd-frequency pairing and the odd-frequency inter-band part even becomes the clearly dominant contribution in this low µ regime. However, by increasing the chemical potential to larger values, the conventional even-frequency intra-band part increases, and eventually becomes dominant, with the total Meissner response now instead approaching its value. The transition between these two regimes is marked by the system going into a metallic regime around µ ∼ 0.08. It is thus not surprising that inter-band processes are dominating in the insulating and near-metallic regime, µ < 0.08. Thus we find that the odd-frequency contributions to the Meissner effect is completely dominating in the insulating and nearmetallic regimes. Notably, the total Meissner effect is always diamagnetic, even with dominating odd-frequency contributions as they are also diamagnetic.
To further investigate the behavior of the Meissner effect when it is dominantly coming from the oddfrequency pairing, we plot in Fig. 5(c,d) the Meissner contributions as a function of δ − for µ = 0. These subpanels clearly show how, by increasing the δ − term both the intra-and inter-band contributions to the odd-frequency Meissner effect become large, especially for moderate t 12 ∼ 0.05 (c). While the intra-band contribution is paramagnetic, it is much smaller than the diamagnetic interband contribution, which makes the total Meissner response follow closely that of the odd-frequency inter-band contribution. Based on these results, we conclude that for inverted, or topological, band structures, odd-frequency pairing can easily dominate the Meissner response, especially at low doping levels in the near-metallic regime. It is inter-band odd-frequency processes that dominate in this regime and they always contribute diamagnetically to the Meissner effect. As a consequence, it is only the large odd-frequency pairing that stabilizes the superconducting state in a magnetic field.

IV. CONCLUSIONS
Odd-frequency pairing has often been discussed to give a paramagnetic Meissner effect, with only conventional even-frequency pairing assume to give the diamagnetic Meissner effect necessary to stabilize a superconductor in an external magnetic field [25][26][27][28][29][30][31][32]. As a result, the existence of large odd-frequency pair amplitudes has been thought to destabilize superconductivity, possibly even to the degree where it cannot exist as a bulk effect [30]. Here we disprove this simplistic picture for multi-orbital superconductors where odd-frequency bulk superconducting pairing is ubiquitous. We show that odd-frequency pairing actually generates both dia-and paramagnetic contributions, and that they usually either nearly cancel each other or even generate an overall diamagnetic Meissner effect. As a result, even large odd-frequency pairing in the bulk does not destabilize superconductivity in multi-orbital superconductors.
In more detail, we are able to derive, using only a few non-restrictive assumptions, simple analytical results for the Meissner response in a generic two-orbital superconductor, divided up into its contributions from evenand odd-frequency pairing, as well as coming from intraand inter-band processes. In a two-orbital superconduc-tor, odd-frequency pairing is always present as soon as there exists a finite orbital hybridization and an asymmetry in the superconducting pairing between the two orbitals, with the odd-frequency pair amplitude growing linearly with the latter. We show analytically that for odd-frequency pairing, intra-and inter-band processes contribute with different signs to the Meissner effect, with the signs directly tied to the character of the normal state band structure.
For two electron-like (or hole-like) bands (i.e. same sign of the curvature in both bands) we always find that the odd-frequency pairing gives diamagnetic intra-band but paramagnetic inter-band contributions to the Meissner effect. Numerically we further show that, while this inter-band part is usually slightly larger, the two contributions nearly cancel in all relevant parameter regimes. Together with a relatively large even-frequency diamagnetic Meissner response, this results in an overall stabilizing diamagnetic Meissner effect, even in the limit of large odd-frequency pairing in the superconductor. In contrast, for an inverted band structure (i.e. with different signs of the curvatures in the two bands, as found in topological materials), we instead always find that the odd-frequency pairing gives diamagnetic inter-band but paramagnetic intra-band contributions. Very interestingly, the odd-frequency inter-band contribution becomes even the dominating contribution to the total Meissner effect in the near-metallic regime where the Fermi level is close to the valence or conduction band bottom. Thus, the Meissner effect can even be completely driven by the odd-frequency pairing, but still be diamagnetic. Moreover, for both types of band structures we find that the even-frequency contributions to the Meissner effect is often diamagnetic but can actually become paramagnetic for some parameter ranges. Still, we always find a stable diamagnetic Meissner effect, even in superconductors with large odd-frequency bulk pairing.
To summarize, our results show that bulk oddfrequency bulk pairing does not cause a destabilizing paramagnetic Meissner effect in multi-orbital superconductors. Instead, odd-frequency pairing either does not contribute significantly to the Meissner effect or, very interestingly, gives a diamagnetic contribution in topological band structures. Based on this diamagnetic oddfrequency Meissner response in topological materials, we speculate that odd-frequency pairing might be quite common in a wide range of topological materials. In Fig. 6 we plot the behavior of the even-and oddfrequency and intra-and inter-band contributions to the Meissner effect as a function of the k-independent orbital hybridization t 12 . We further set δ + = 0.03, δ − = 0.04, δµ = −0.2 and use in (a) an electron-like band structure with t 2 = 0.5 and µ = 0.3 and in (b) an inverted band structure with t 2 = −0.5 and µ = 0.0 in order to comply with the parameter choices in the main text. These superconducting order parameters renders a s +− state with substantial odd-frequency pairing due to the large size of δ − . The results show that for the region t 12 0.1 all different contributions to the Meissner effects are reasonably large, while for larger t 12 the even-frequency intraband contribution becomes dominant. This motivates our choice of t 12 ≤ 0.1 in the main text in order to capture the most interesting regime. Still, all conclusions drawn in the main text, in terms of both magnitudes and signs of the Meissner effect, hold true for the full range of t 12 . Next we investigate how the Meissner response changes if we instead use a momentum dependent orbital hybridization, ξ 12 (k). In Fig. 7 we show the different contributions to the Meissner effect as a function of the orbital energy difference δµ for a kdependent hybridization between orbitals of the form ξ 12 (k) = t 12 k 2 . We further set δ + = 0.03, δ − = 0.04, creating a s +− state and report results for both an electron-like band structure (a,c) with t 2 = 0.5, µ = 0.3, and an inverted band structure (b,d) with t 2 = −0.5, µ = 0.0. We also vary the hybridization strength, with t 12 = 0.05 (a,b) and t 12 = 0.2 (c,d). The results in these figures are in very good agreement with figures in the main text, Figs. 2 and 4, regarding the sign and magnitude of different contributions to the Meissner effect, and thus do not qualitatively change any of our conclusions. Exactly the same as in the main text, the odd-frequency intra-and inter-band contributions to the Meissner effect have opposite signs with respect to each other and also inverting the band structure flips the signs of these two contributions. Due to the opposite signs of the intra-and inter-band contributions, the total Meissner contribution from odd-frequency pairing is always small and hence even large odd-frequency pairing does never produce a destabilizing Meissner response, even for a momentum dependent orbital hybridization. Additionally, we find for the case t 12 = 0.05, that the odd-frequency contribution is relatively small, but it increases when we use t 12 = 0.2 in panels (c,d). This can easily be understood from the fact that the t 12 -term is now multiplied by a k 2 dispersion, which reduces its effect at low momenta. Hence, in order to get a similar behavior to that of the k-independent case with ξ 12 = t 12 , we need to use larger values of pre-factor t 12 .

Appendix B: Band basis Hamiltonian
In this appendix, we transform the two-orbital Hamiltonian in Eq. (2) into the band basis where the kinetic energy is fully diagonal. While the original orbital basis of Eq. (2) is convenient for a lot of the derivations in the main text, the band basis gives a better intuition for some of our results.
The normal-state Hamiltonian in the main text, h(k) = ξ + τ 0 + ξ − τ z + ξ 12 τ x , is diagonalized and thus transformed into the band basis with the rotation ma-trixR θ/2 = cos(θ/2) − sin(θ/2) sin(θ/2) cos(θ/2) , where θ = sin −1 ξ 12 / ξ 2 12 + ξ − . Applying this rotation matrix also to the superconducting part of Eq. (2) we retrieve the Hamiltonian in the band basis as: The Hamiltonian H b in the band basis nicely illustrates two features. One, the only coupling between the two bands is through the inter-band pairing terms, proportional to δ − . Thus, in the δ − → 0 limit, the superconducting part of the Hamiltonian H b would commute with the normal part. As a result, odd-frequency pairing is a direct consequence of this inter-band term [7]. This shows that odd-frequency pairing from the interband pairing is fundamentally the same, and only differ with a basis change, from the odd-frequency pairing generated by orbital hybridization ξ 12 and orbital pairing asymmetry δ − used in the main text. Two, the intraband pairing terms becomes explicit as δ + ± δ − cos(θ) and it is easy to see when one of the bands become nonsuperconducting, a result used in the main text to explain some of the even-frequency intra-band results.