Effect of interorbital scattering on superconductivity in doped Dirac semimetals

Unconventional superconductivity has been discovered in a variety of doped quantum materials, including topological insulators and semimetals. A unifying property of these systems is strong orbital hybridization, which leads to pairing of states with nontrivial Bloch wave functions. In contrast to naive expectation, however, many of these superconductors are relatively resilient to disorder. Here we study the interplay of superconductivity and disorder in doped three-dimensional Dirac systems, which serve as a paradigmatic dispersion in quantum materials, using Abrikosov-Gor’kov theory. In this way, the role of disorder is captured by a single parameter (cid:2) , the pair scattering rate. In contrast to previous studies, we argue that interorbital scattering cannot be neglected due to the strong orbital hybridization in Dirac systems. We ﬁnd that the robustness of different pairing states highly depends on the relative strength of the different interorbital scattering channels. In particular, we ﬁnd that the “nematic” superconducting state, which is argued to be the ground state in many Bi 2 Se 3 -related compounds, is not protected from disorder in any way. The pair scattering rate in this case is at best smaller by a factor of 3 compared to systems without spin-orbit coupling. We also ﬁnd that the odd-parity pairing state with total angular momentum zero (the B phase of superﬂuid 3 He) is protected against certain types of disorder, which include a family of time-reversal odd (magnetic) impurities. Namely, this odd-parity state is a singlet of partners under CT symmetry (rather than T symmetry in the standard Anderson’s theory), where C and T are chiral and time-reversal symmetries, respectively. As a result, it is protected against any disorder potential that respects CT symmetry. Our procedure is very general and can be readily applied to different band structures and disorder conﬁgurations.


I. INTRODUCTION
Anderson's theory explains why conventional s-wave superconductors are weakly affected by nonmagnetic disorder [1][2][3].It is based on two essential conditions.The first one is that the Cooper pairs in these superconductors form singlet states of time-reversed partners.The second one is that the phase of the pair wave function is featureless over the entire Fermi surface.These two ensure that the pairing interaction, written in the basis that diagonalizes the disorder potential, remains the same as in the clean limit.Consequently, one can always pair time-reversed partners with the same interaction and the same transition temperature [4].In contrast, the Cooper pairs in unconventional superconductors violate one of these conditions and, as a result, are not protected [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22].
Mechanisms based on the huge spin-orbit coupling characterizing the topological materials have been suggested to explain this robustness.The authors of Ref. [48] studied the effect of disorder on an odd-parity paring state with zero total angular momentum (equivalent to the B phase in superfluid 3 He).They found that an additional chiral symmetry can protect this state from certain types of disorder when it is present.It was later suggested that the pair wave function in doped Bi 2 Se 3 is a multicomponent nematic state which breaks the rotational symmetry of the crystal [39,49].References [46,50] studied the effect of disorder on the nodal nematic state, also arguing for some robustness; however, their results were based on less generic grounds.Studies of the effects of scalar disorder on unconventional paring in topological materials were recently discussed for materials other than Bi 2 Se 3 [47,51] and also in the context of the surface of topological materials [52][53][54][55].
The studies mentioned above, however, focus only on the effects of intraorbital scattering, which has equal weight on all orbitals (i.e., density disorder).Because topological materials FIG. 1. Different types of nontrivial time-reversal-symmetric intra-and interorbital disorder potentials in Bi 2 Se 3 (as a prototypical Dirac system).The solid and dashed lines represent the Se and Bi layers in the quintuple unit cell, respectively.The purple p z orbitals represent the itinerant states on the top and bottom Se layers, which disperse as Dirac fermions close to the point.(a) Disorder that breaks the symmetry between top and bottom layers induces γ 5 and γ 3 potentials [our convention for the γ matrices is given below Eq. (1)].This can be caused by a polar impurity or a charged impurity closer to one layer than the other (e.g., due to the intercalation of Cu).(b) Disorder that modifies the tunneling between top and bottom (either by causing a barrier or by stretching and/or squeezing the z axis) induces disorder in the mass term of the form γ 0 .(c) An in-plane polar impurity induces γ 1 and γ 2 potentials.We note that in all cases we also anticipate an intraorbital density potential (proportional to the identity matrix).
are multiorbital systems with huge orbital hybridization, interorbital scattering is not expected to be particularly weaker than density disorder.In Fig. 1, we schematically depict three types of time-reversal-symmetric disorder potentials, which are expected to be present in Bi 2 Se 3 and lead to interorbital scattering.Thus, it is important to understand the effects of interorbital scattering in the superconducting topological materials.
The influence of interorbital scattering on pairing was first emphasized by Golubov and Mazin [56] and was later studied in the context of systems with multiple Fermi surfaces (see, for example, Refs.[57][58][59][60][61]).We emphasize that in topological materials the multiorbital nature is embedded in the Bloch wave functions rather than the presence of multiple Fermi surfaces, making them somewhat different.
In this paper, we study the effect of short-ranged intraand interorbital scattering on superconductivity in threedimensional materials with Dirac dispersion.The Dirac dispersion is a paradigmatic example of dispersion relations in topological materials, which also naturally have large spinorbit coupling [62].
We provide an extensive picture of how the transition temperature T c in different pairing channels is affected by all possible types of short-ranged scattering potentials.Our results are expressed in terms of the pair-breaking rate , which also enters the Abrikosov-Gor'kov theory of superconductors with magnetic impurities [2].We first discuss the case of a massless Dirac dispersion where the pair-breaking rates are proportional to the single-particle scattering rates through universal rational numbers (see Table III).We then compute how these scattering rates vary with mass, which allows us to interpolate to the well-known results for the systems with no spin-orbit coupling [2,5].
We find that the only state robust to time-reversal symmetric (TRS) disorder is the s wave, in agreement with the Anderson's theorem [1].Another fully gapped isotropic state is the odd-parity state with zero total angular momentum (analogous to the B phase in superfluid 3 He), which does not exhibit such robustness, in contrast to previous expectations [47,48].Despite being protected against certain types of disorder, it turns out to be very sensitive to some other types of defects, such as mass and polar impurities, which we expect to be generally present in topological materials.
On the other hand, we find that this fully isotropic oddparity state is protected from any disorder that respects CT (assuming this symmetry is also present in the clean system), where C and T are chiral and time-reversal symmetries, respectively.The reason for such robustness is that this state corresponds to a singlet pairing state of CT partners, in perfect analogy to pairing of T partners in Anderson's original argument [1].Interestingly, this result implies that the fully gapped odd-parity state is protected against certain disorder potentials that are odd under time reversal T .
We also study the multicomponent states, which are the O(3)-symmetry group related to the nematic pairing states in doped Bi 2 Se 3 .The relation between these states is obtained when the symmetry group is reduced from the fully isotropic O(3) to trigonal D 3d group of Bi 2 Se 3 .We find that the pair breaking rate for these states can be smaller than the corresponding rate in systems without spin-orbit coupling.However, they are still significantly influenced by disorder in contrast to previous studies [46,50].We also consider the effect of magnetic impurities, where the O(3)-symmetry pairing state can also be slightly more protected than in systems without spin-orbit coupling [5] (depending on the microscopic nature of the magnetic impurities).It should be added that when time reversal is broken, the chiral state might be preferred over the nematic one [63][64][65][66][67].
The rest of this paper is organized as follows.In Sec.II, we present the basic ingredients of our model, namely, an action for a Dirac fermion subjected to a generic disorder potential and an attractive pairing interaction.We project the disorder onto the Bloch basis of the conduction electrons near Fermi surface, which we will use throughout this paper.In Sec.III, we calculate the pair-breaking rate in topological materials, which is the main parameter that affects T c and the whole thermodynamics.Our results are presented in Sec.IV, where we discuss the effect of different types of disorder on the various pairing channels.Finally, Sec.V provides a summary of our main results and a discussion of the application of our method for other systems.Multiple technical details of our calculation are delegated to the Appendices.

II. THE MODEL
We start by describing the normal-state action of the model.We consider massless Dirac fermions where ) and ψ ω,k,± correspond to two orbitals, each consisting of a Kramers pair.The orbitals are related to each other through inversion.v is an isotropic velocity and F is the Fermi energy.In addition, the γ matrices are taken to be Hermitian, γ = τ 2 s, γ 0 = τ 1 s 0 , and  2) and (3).The table also lists the discrete symmetry properties of each scattering process under I = γ 0 = τ 1 s 0 , T = Kγ 1 γ 3 = Kτ 0 (−is 2 ), and C = γ 5 = −τ 3 s 0 , corresponding to inversion, time-reversal, and chiral symmetries, respectively.ζ = +1 (−1) corresponds to Fermi level residing in conduction (valence) band, accordingly.
, where τ and s are Pauli matrices in the orbital and spin basis, correspondingly.δ is the mass of the Dirac point.In the following, we set δ = 0 for simplicity unless specified explicitly otherwise.Note that we have neglected higher order corrections in momentum.
The action in Eq. ( 1) can be conveniently diagonalized in the manifestly covariant Bloch basis (MCBB), in which the electron spinor transforms as an ordinary SU(2) spin-1/2 [40,49,68,69] , where ζ = ±1 corresponds to conduction-valence band, respectively, k j = k j /k, and k± = kx ± i ky .Without loss of generality, we assume electron doping (ζ = 1), and therefore omit index ζ henceforth.The field operators are then approximated by their weight on the band operators Finally, the action in Eq. ( 1) possesses inversion, chiral, and time-reversal symmetries.The representation of these symmetry operations in orbital basis is given by I = γ 0 , C = γ 5 , and T = Kγ 1 γ 3 , respectively, where K is complex conjugation.
Next, we consider the disorder potential.The crucial element in our theory is the inclusion of interorbital scattering.Within the Dirac notations, such (momentum-independent) scatterings can be represented using the Dirac matrices introduced above and their products.For elastic short-ranged scattering (compared to k F ), we have where N m is the number of impurities in channel m and r l is the position of these short-ranged impurities.There are 16 different Hermitian matrices M m representing different types of disorder.The representation of these matrices in the orbital-spin τ ⊗ s basis and their discrete symmetry properties are given in Table I.We further clarify that m = 0 corresponds to simple density disorder 1, which was considered in Refs.[46,48,50,51].m = 1 corresponds to mass disorder γ 0 , m = 2 is odd-parity scalar disorder γ 5 , m = 3, 4, 5 correspond to odd-parity dipolar disorder, and the magnetic disorder m = 6, 7, 8 correspond to a fully symmetric local moment with spin oriented along the axis S α = −iε αβγ γ β γ γ .Some examples of nontrivial scattering matrices of this type, which naturally appear in disordered Bi 2 Se 3 , are shown schematically in Fig. 1.When projecting the disorder potential onto the MCBB, we obtain a set of scattering matrices Q m ( p, k) with nontrivial momentum dependence: where ), and we defined the matrices Before proceeding, we make a few important remarks regarding the choice of disorder potential in Eq. ( 2).First, we assume that the disorder is Gaussian correlated with zero mean, V l,m = 0. Next, we assume that there are no spatial correlations and that different types of disorder do not correlate.The latter assumption implies that disorder potential does not break any symmetry on average, leading to V l,m V l ,m = V 2 m δ mm δ ll .There is one exception, however, which requires clarification.In the absence of chiral symmetry, the density disorder (m = 0) and the mass disorder (m = 1) can, in principle, mix, since they belong to the same trivial representation.We note that this is the reason why mass disorder is always present in Bi 2 Se 3 , even if it respects time-reversal and inversion symmetries (as opposed to the claim made in Ref. [46]).As we show in Sec.III, however, the correlations between m = 0 and m = 1 disorder channels do not affect the scattering rate or superconductivity.
Second, a central assumption of our theory is that disorder naturally appears in the orbital basis.Indeed, the set of matrices Q m ( p, k) introduced in Eq. ( 3) and listed in Table I is the result of starting from the orbital basis and projecting disorder potential onto the MCBB on the Fermi surface.The additional momentum-dependent form factors in the scattering matrices could have been easily overlooked if we started directly from the band basis, constructing a phenomenological picture of disorder [70].To emphasize this fact, we point out that even the density channel obtains nontrivial momentum dependence, which, as we show below, plays a crucial role in protecting some unconventional pairing states from density disorder.
The last ingredient required to estimate the superconducting transition temperature is the attractive interaction which leads to the instability.We study the superconducting instability in the band basis, in the spirit of the Bardeen-Cooper-Schrieffer (BCS) theory.
We then decompose the interaction into the irreducible representations in the Cooper channel: where F J ( k) are form factors in the MCBB corresponding to different representations J of the relevant symmetry group, which are specified in Table II.A superconducting instability can occur in any one of the channels depending on the attractive strength of coefficients g J .We are mainly interested in systems with large spin-orbit coupling, characteristic for topological materials, which do not have spin-rotational symmetry.That is why in this paper we focus on the fully isotropic O(3) group of joint rotations of spin and momentum.Different representations are labeled by the total angular momentum J. (Note that within our notations J labels both different representations and different components within the same representation.)

III. COMPUTATION OF THE SCATTERING RATE
We now turn to the computation of the pair scattering rate , which enters the Abrikosov-Gor'kov theory and dictates the thermodynamics of superconductors.The procedure we employ consists of three main steps, which are described diagrammatically in Fig. 2.

A. Single-particle lifetime
We start with computing the single-particle lifetime.The bare electronic Green's function in the band basis is given by (ζ = 1) TABLE II.Different representations of the time-reversalinvariant order parameters F J .The labels L, S, and J correspond to orbital, spin, and total angular momentum, respectively.Note that for every total angular momentum J there are 2J + 1 states |J, J z , where J z gets integer values between −J and J. Thus, the two different states with J = 0 both have J z = 0 and are distinguished by their transformation properties under inversion, either even (g) or odd (u). ) Summation of the diagrams in Fig. 2(a) leads to a selfenergy correction to the Green's function, (3) and ( 5), we find that the self-energy is given by (iω) = m m (iω), with A diagrammatic representation of the summation of the Gor'kov ladder.(a) The summation over the scattering processes from all types of impurities within the first Born approximation, which leads to the self-energy correction to the full Green's function, Eq. ( 6).(b) The Cooperon vertex correction B J (iω), Eq. ( 12), which results from the summation over the bare pairing propagators A(iω), Eq. ( 9).(c) The summation of the Gor'kov ladder.The building block of the ladder, P J , is given by the sum of B J (iω) over Matsubara frequencies; see Eq. ( 14). where Here n m = N m /L 3 is the density of impurities in channel m, which arises after averaging over the positions of the impurities, and ν 0 = k 2 F /2π 2 v is the density of states at the Fermi level per pseudospin.We note a factor of 2 difference in the definition of the scattering time compared to a parabolic band (see Appendix A), which is a feature of topological touching points of two bands.Thus, we have obtained that the single-particle scattering rate decomposes into a sum over the different scattering channels: We recall that disorder is uncorrelated among different channels, which results from the assumption that disorder does not break any symmetry on average.Mass and density disorder are an exception, since they both belong to the trivial representation.However, we note that even if cross correlations between mass and density are present, for δ = 0 the cross term in Eq. ( 6) vanishes.To see this, we write the product where we have used the identities d k = d (cos θ k )dφ k , kα|γ 0 | kγ = 0, and d p |pβ pβ| = 2π 1, and the summation over repeated index β is implied.

B. Vertex correction
In addition to the single-particle processes, it is also important to take into account the effect of pair scattering.Namely, now we calculate the correction to the BCS vertex due to intermediate scattering on disorder.In the limit of weak disorder, F τ 1, the most important correction to the Gor'kov ladder comes from the diagrams with nonintersecting impurity lines (so-called Cooperon), as shown in Fig. 2(b) [3].
To compute the disorder contribution to the vertex, we need two ingredients.First, we calculate the bare propagator of a Cooper pair: This propagator links between the scattering events on the Gor'kov ladder and corresponds diagrammatically to the first term on the right-hand side of Fig. 2(b) [i.e., it forms the legs of the ladder].When decomposing a generic pairing interaction into the irreducible representations [as in Eq. ( 4)], the Gor'kov ladder decomposes into scattering channels of the orthogonal basis functions F J ( k), which are labeled by FIG. 3. A diagrammatic representation of the decomposition from particle-hole to particle-particle channels; see Eq. (10).Matrices M, Q, F, and b are defined in Eqs. ( 2)-( 4) and (11) [see also Tables I and II and Eq J. Thus, when writing Eq. ( 9), we assume that the Cooper pair propagator is contracted on both sides with the interaction lines in the corresponding channel J.
The next important ingredient for calculating the Gor'kov ladder is the scattering amplitude from a single impurity in the particle-particle basis, which is the building element of a Cooperon and shown diagrammatically as the second term on the right-hand side of Fig. 2(b).Thus, we need the scattering amplitude of a Cooper pair with any momenta k and −k into a pair with any other momenta p and −p due to an impurity of type m.This amplitude is given by the product of two singleparticle events: where is a matrix of weights corresponding to the conversion from the particle-hole to particle-particle basis (similar to the Fierz identity [42]) and d k = d (cos θ k )dφ k is the solid angle element.The matrix in Eq. ( 11) is given explicitly in Appendix E.Here we only note that |b Jm | 1/2.The procedure of decomposition from particle-hole to particle-particle basis is shown schematically in Fig. 3. Contracting the two ingredients, Eqs. ( 9) and (10), and using the orthogonality of the superconducting form factors F J ( k), we obtain for the disorder corrected block of the Gor'kov ladder, shown schematically in Fig. 2(b): where is the pair scattering rate, which is a sum of independent scattering rates Jm originating from the different intra-and interorbital disorder channels m.The values for the partial pair scattering rates from Eq. ( 13) are the main result of this paper and are listed in Table III.
As mentioned below Eq. ( 1), up to this point we have focused on the case of zero mass, δ = 0.In this case, the product of τ m Jm takes the universal rational values presented in Table III.To discuss how these numbers vary with a finite TABLE III.Top: Indicates different types of intra-and interorbital scattering matrices and their properties under the discrete symmetries.The scattering matrices are labeled by m = 0, . . ., 15 and appear as γ matrices and their products.T = Kγ 1 γ 3 , I = γ 0 , and C = γ 5 correspond to time-reversal, inversion, and chiral symmetries, respectively.Bottom: The values of the dimensionless pair scattering rate τ m Jm = 1 − 2b Jm [where the matrix b Jm is defined in Eq. ( 11)] for different superconducting pairing states.The labels L, S, and J correspond to orbital, spin, and total angular momentum, respectively (we use same convention as in Ref. [68]).It is clear that the F 0g state is protected from disorder that respects T symmetry (Anderson's theorem), while the F 0u state is protected from disorder that respects CT symmetry.In the last line, F 2 implies all the pairing states with J = 2, i.e., F 21 − F 25 .The result for all of these states is the same.
mass δ, we introduce a parameter α = δ/ F , which ranges between 0 (no mass) and 1 (infinite mass).As a consequence, the single-particle rates in Eq. ( 7) and the pair-breaking rates in Eq. ( 13) are modified to Eqs. (A5) and (A11), respectively (for details, see Appendix A).We can distinguish three cases.The first one is the case of s-wave pairing (F 0g ), in which the mass does not affect the results in Table III.The second case is J = 0g (non-s-wave) and the disorder matrices are inversion symmetric (i.e., I −1 M m I = M m ).In this case, the values of Jm τ m continuously interpolate between those in Table III and the asymptotic value without spin-orbit coupling Jm τ m = 1 which was computed by Larkin [5].The functional form of this interpolation is given by Eq. (A11) and plotted in Fig. 4. Finally, the third case is J = 0g and disorder matrices are odd under inversion (i.e., I −1 M m I = −M m ).In this case, the matrices M m act as purely interband operators within the Bloch wave functions when taking the limit α → 1.Consequently, Jm and 1/τ m go to zero as the conduction and valence bands become infinitely separated, i.e., when α increases.It is interesting, however, that their ratio remains the same for all α and is given in Table III.

C. Computation of T c
The final step of the calculation is to use the disordermodified interaction in the superconducting channel J to compute the renormalized pairing vertex.To do that, we insert a Cooperon in each block of the Gor'kov ladder [this step gives us factor B J (iω)], perform the summation over intermediate Matsubara frequencies ω n = π T (2n + 1), and sum up all the blocks [as shown in Fig. 2(c)].The result reads as gJ = g J /(1 + g J P J ), (14) where P J = −T ω n B J (iω n ) and T is temperature.The transition temperature T c,J in the channel J is determined as a singularity (vanishing denominator) in Eq. ( 14), leading to the result that has the Abrikosov-Gor'kov form where (x) is the digamma function and T c,J,0 is the transition temperature in the absence of any disorder.For the detailed calculation of the transition temperature, see Appendix B. An alternative derivation of this result via the method of Gor'kov Green's functions (which also gives a solution below T c ) is presented in Appendix C. The solution of Eq. ( 15) predicts that superconductivity is completely suppressed at Jcr. = π e −γ e T c,J,0 ≈ 1.76 T c,J,0 (16) (γ e ≈ 0.577... is the Euler's constant), as shown in Fig. 5.
FIG. 4. The dependence of the pair scattering rate τ m Jm on Dirac mass α = δ/ F (where F = (vk F ) 2 + δ 2 ), for odd-parity pairings F 0u (blue) and F 11 (red).The upper half is the effect of mass disorder γ 0 , whereas the lower half shows magnetic disorder iγ 1 γ 5 .While the former has the most severe effect on F 0u , this pairing channel is protected against the latter due to the CT symmetry.Nevertheless, they all approach the Larkin's result τ m Jm = 1 (black) as the mass goes to infinity (α = 1).FIG. 5. T c /T c,0 as a function of /T c,0 obtained from the solution of Eq. ( 15).

IV. RESULTS
In Sec.III, we described how the pair scattering rates in Eq. ( 13) affect the transition temperature.As we show in Appendix C, the effect of these rates is actually much more general as they dictate the entire low-temperature thermodynamics of these superconductors (up to phase fluctuation effects) [2,3].For example, we recall that the gap may close in the superconducting state when disorder is sufficiently strong.
Having established the importance of the pair scattering rates in Eq. ( 13) for superconductivity in doped Dirac systems, we now turn to discuss their value for different pairing states and different interorbital disorder potentials.The results are summarized in Table III, which includes both nonmagnetic (m = 0, . . ., 5) and magnetic (m = 6, . . ., 15) impurities.
As mentioned above, the elements of the matrix b Jm in Eq. ( 11) range between 1/2 and −1/2 (see Appendix E).Consequently, the rates Jm appearing in Table III, which shows the values of 1 − 2b Jm , range between 0 and twice the single-particle scattering rate 1/τ m .The former implies that disorder in channel m does not affect superconductivity in channel J, while the latter corresponds to the most severe effect possible.
Indeed, for the s-wave channel (J = 0g), we find that the pair scattering rate vanishes, 0g = 0, for all the T -even disorder matrices m = 0, . . ., 5, which is manifestly the Anderson's theorem for nonmagnetic impurities.
Additionally, we recover the well-known Abrikosov-Gor'kov result for magnetic impurities, that the pair scattering rate is twice that of single particles [2], such that overall the pair-breaking rate is given by Another limit of interest is the odd-parity state with total angular momentum zero (J = 0u), which is equivalent to the B phase of superfluid 3 He [71].We notice that, similar to the swave state F 0g , this pairing state is completely protected from certain types of disorder, namely m = 0, 2 and 12, . . ., 15. Inspecting Table III, we identify that the common symmetry of these disorder potentials is that they are all even under the product of chiral and time reversal CT .What makes this result even more interesting is that some of the CT -even matrices are T odd.Thus, the F 0u is protected from certain types of magnetic impurities.We identify this protection with similar results for superconductors with multiple Fermi surfaces [56].
To understand this protection, we now show that the F 0u is essentially a singlet pairing state between partners related to each other by CT symmetry.Thus, in complete equivalence to the Anderson's original argument [1], it follows that as long as the disorder potential does not violate CT symmetry, we can always pair CT partners in the basis that diagonalizes the disorder potential (see Appendix D).
Let us show that the F 0u pairing state is indeed a singlet state of CT partners.This is most easily seen in the orbital basis.We find it convenient to rotate the orbital basis by π/2 about the τ 2 axis first.This transforms from the basis of chirality to the basis of parity (i.e., the orbitals are labeled by their parity τ = ±).Note that this does not affect the operation of time-reversal T .Then the action of chiral symmetry is implemented by C = γ5 τ s 0 and the corresponding pairing state where ψ sτ (k) is a field operator in the rotated basis.Inspecting this pairing state, it is evident that it is fully antisymmetric and that each term consists of a pair of operators related to each other by CT symmetry.However, the F 0u state becomes vulnerable to disorder when CT symmetry is not present, such as in doped Bi 2 Se 3 .In that case, mass belongs to the same representation as density and is always present.Moreover, we argue that Dirac materials are often polar ionic crystals (e.g., Bi 2 Se 3 , SnTe, PdTe, etc.), and therefore it is likely that the disorder potential also induces dipolar moments of type m = 3, 4, 5.This argument should be contrasted with the claim made in Ref. [72], where it was stated that only density disorder should be present.Overall, we find that the pair-breaking rate in the F 0u channel equals It should also be noted that the authors of Ref. [48] were the first ones to identify that the F 0u state can be protected from disorder in the massless limit.However, they concluded that it is protected by C symmetry.As we show here, it is actually protected by CT symmetry.To emphasize this distinction between the two, we point out that the F 0u state is immune to some disorder potentials that are odd under C, such as m = 12, . . ., 15.
Next, we consider the odd-parity pairing states with total angular momentum J = 1.These states (more accurately, the nematic E u states of the D 3d symmetry group of Bi 2 Se 3 , which derive from the J = 1 representation by breaking the rotational symmetry from spherical to trigonal) are of special interest experimentally, since they are considered to be the pairing state in doped Bi 2 Se 3 [28][29][30][31][32][33][34][35]39,49].We find that all disorder channels affect superconductivity with this pairing symmetry, and the dimensionless rate τ m m takes two possible values, 5/3 and 1/3.Thus, depending on the relative weight in these channels, the scattering rate can vary significantly.In particular, for density disorder, the rate 1 j,0 = 1/3τ 0 ( j = 1 − 3 here; see Tables II and III) is much smaller than in systems without spin-orbit coupling, where it is expected to be 1/τ 0 [5,11].However, our results are not consistent with experiments that find this state to be protected [45,46] and also not consistent with previous theoretical work where density disorder was considered [46,50], and even more so when considering the m > 0 channels, which are more harmful.For example, if polar disorder is present, then an average over all possible directions gives 1 j,3−5 = 11/9τ 3 , where we assume τ 3 = τ 4 = τ 5 by symmetry and, again, j = 1 − 3.

V. DISCUSSION
In this paper, we construct rigorous method to evaluate the effect of both intra-and interorbital disorder on superconductivity in doped Dirac materials.We argue that generic disorder potential always induces interorbital scattering processes given by Eq. ( 2) and listed in Table I.We compute the contribution of each type of intra-and interorbital scattering channels to the pair-breaking rate for a given pairing potential, which dictates the entire thermodynamics of a superconductor.This result is summarized in Table III.Our main conclusions from this analysis are as follows: (1) We have found a version of the Anderson's theorem which is based on the pairing of CT partners rather than T partners (C is chiral and T is time-reversal symmetry).The odd-parity state with total angular momentum zero (J = 0u) is such a pairing state.Consequently, it is protected from disorder that respects CT symmetry, which includes certain types of T -odd impurities.This also generalizes the results of Ref. [48].It is interesting to understand in the future if such a symmetry can exist (or nearly exist) in a solid-state material.
(2) As expected from the Anderson's theorem, the s-wave (J = 0g) pairing state is protected from all nonmagnetic scattering processes, including interorbital ones.
We emphasize that the results presented in Table III are based on a model where the mass term in the single-particle Hamiltonian, Eq. ( 1), was assumed to be zero.It should be noted, however, that in the majority of doped topological materials which become superconducting such a mass term exists.The dependence of these numbers on the mass is discussed in Appendix A and plotted in Fig. 4. For inversionsymmetric disorder potentials, the inclusion of this term modifies these numbers toward their known values without spin-orbit coupling [2,5].In particular, it removes the protection of the F 0u state.For inversion-odd disorder, on the other hand, the values in Table III remain unchanged.For more details, we refer the reader to Appendix A.
The analysis performed in this work assumes a shortranged disorder potential.However, in doped materials, one may also anticipate a correlated potential emerging from charged impurities [73].Therefore, it is important to also understand the influence of a soft potential on superconductivity.
Another question which was not addressed in this paper is the microscopic origin of pairing and how it is affected by disorder [74].In particular, doped topological materials are characterized by small electronic density and small density of states.As a result, the pairing interaction must be more singular [69].Such an interaction is expected to be sensitive to the presence of disorder [75].
Finally, our results hold only for the case of a finite Fermi energy and weak disorder, implying F and F τ 1.It would be interesting to consider the limit of low density, where both conduction and valence bands are important.In this limit, however, a number of issues arise.First, the assumption k F l 1 breaks down and therefore the self-energy and Born approximations are invalid and other approaches (e.g., the replica approach) must be used [76][77][78].Second, the omission of one of the bands is invalid and all four bands must be taken into account.Finally, we note that this scenario is, however, very exotic and, as discussed above, requires long-ranged interactions [69].
Looking forward, we argue that our theory is useful to many other systems with strong orbital hybridization.In particular, Eq. ( 11) is easily generalizable to different Hamiltonians and reduced dimensions.Of special interest are semimetallic systems, including a quadratic band touching point relevant to the half-Heusler compounds [79], line-node semimetals, Weyl semimetals that emerge when inversion is broken in a Dirac material [69], and higher-order band touching points [80].
We also note that the results in Table III suggest that the effects of disorder on unconventional pairing states depend strongly on the microscopic nature of disorder (we note that a similar disorder dependent pair-breaking rate has been argued to exist in superfluid 3 He in a nematic aerogel [81]).This opens an interesting avenue to manipulate the superconducting ground state by selectively inducing specific types of disorder potentials.
Before concluding this paper, we note that arguments for robustness of unconventional superconductivity to disorder were recently cast in terms of the so-called superconducting fitness [46,47,51], which was first discussed in Ref. [82] in the context of clean systems.An intuitive understanding of the fitness can actually be obtained based on Ref. [83], where the affects of disorder on superconductivity were assessed by looking at the minimal excitation of a system and comparing it with the clean limit.As shown in Ref. [47], this translates to the condition that the Hamiltonian including disorder commutes with the gap function, [ Ĥ + V , ˆ ] ± = 0, where the ± stands for commutation-anticommutation for TR even and TR odd disorder, respectively.Our results are consistent with this picture.
for m = 2-5, 9-12) does not affect superconductivity in the case of an infinite mass, α → 1.Indeed, all matrices Q m are independent of momentum in the limit α = 1, implying that inversion-odd disorder only scatters between the conduction and valence bands.Such transitions are obviously suppressed in the limit of large band gap, leading to the divergent single-particle scattering time τ * m → ∞ and, consequently, vanishing pair-breaking rate * Jm → 0. Finally, we comment on how finite mass affects the robustness of J = 0u channel.While finite mass breaks CT symmetry and the channel becomes susceptible to most of the types of disorder, it is still unaffected by γ 5 (m = 2) and iγ 0 γ 5 (m = 12), which not only respect CT symmetry, but also odd under inversion.

APPENDIX B: CALCULATION OF T c
Using the condition for the superconducting instability, Eq. ( 14) of the main text, we obtain that P ≡ −T c ω B(iω) = − 1 g .Note that the derivation in this Ap-pendix is independent of the pairing channel, so we omit the subscript J for brevity.We now plug in the result for B(iω) as obtained in Eq. ( 12) and find Using the definition of Matsubara frequencies ω n = 2π T (n + 1/2), this equation can be rewritten as Following Abrikosov and Gor'kov [2], we make use of the fact that in the clean limit one has where γ e is the Euler's constant.Thus, we can rewrite Eq. (B2) as We can identify the two terms inside the square brackets as digamma functions and we know from Eq. (B3) that 1 gν 0 = log ( 4e γe π ω D 2T c,0 ), where T c,0 is the critical temperature for a clean system.Consequently, we obtain which coincides with Eq. ( 15) of the main text.

APPENDIX C: ABRIKOSOV-GOR'KOV EQUATIONS AT ARBITRARY TEMPERATURE: GAPLESS SUPERCONDUCTIVITY
Now we present the complementary approach to derive the effect of disorder on superconductivity, which exploits the formalism of Gor'kov Green's functions [84].The advantage of this method is that it allows us to treat the problem at arbitrary temperature and study thermodynamic and electromagnetic properties of a disordered superconductor at temperatures down to T = 0.
To start with, we introduce the Nambu space (N) for the MCBB electron operators according to In this basis, the bare (without disorder) Gor'kov Green's function takes form [70,85] Ĝ0 where (C5) The self-energy due to disorder is then given by (C6) We notice that the self-consistency requires us to use full Green's function Ĝ, instead of the bare one Ĝ0 .Following Ref. [2], we look for a solution of the form Ĝ Performing integration over ξ k first, we obtain the very general expression which applies to any superconducting state with unitary pairing: The above expression can be easily used to reproduce the result for the transition temperature T c,J , Eq. (15).Neglecting ˜ 2 n,k in the denominator and performing integration over k and summation over m, we find for the channel J where J is given by Eq. ( 13).In deriving the last equation, we also used Eqs.( 10) and (11).Utilizing further Dyson equation Ĝ−1 = Ĝ−1 0 − ˆ , we easily obtain ωn = ω n + sgn( ωn ) 2τ , which can be readily resolved yielding Finally, using the gap equation in the channel J ) we obtain after summation over p which is identical to Eq. (B1) and leads eventually to Eq. ( 15).We emphasize that Eq. ( C10) is very general and can be used to study thermodynamic properties of any (unitary) superconducting state.As an example, we focus on the fully isotropic pairing functions F 0g and F 0u .Performing integration over k in Eq. (C10), we find a set of coupled equations for ωn and ˜ n : accompanied with the gap equation These equations, in principle, can be solved numerically to find the value of the pairing gap at arbitrary temperature and study the thermodynamic properties of a superconductor.For instance, in the case of nonmagnetic (T -even) disorder for the s-wave F 0g pairing or CT -even disorder for the pwave F 0u pairing, we have J = 0, and Eq.(C16) admits simple solution ωn / ˜ n = ω n / .In these cases, the latter result implies that the gap equation (C17) is not modified by disorder at all, consequently, the transition temperature and all the thermodynamic properties below T c remain unchanged compared to the clean case.
Abrikosov and Gor'kov have analyzed Eqs.(C16) and (C17) in detail (with J = 0) for the most interesting limiting cases in Ref. [2].In particular, they found that there is a range of the impurity concentration where superconductivity is not entirely suppressed, while becoming gapless.In our language, this corresponds to the threshold value of the pair-breaking rate J , above which the gap in the spectrum of elementary excitations vanishes: Jcr. ≈ 0.91 Jcr., (C18) where the critical value Jcr. is given by Eq. ( 16).As a result, the low-temperature behavior of the specific heat changes from exponential to T linear in the range J < J < Jcr. .The analysis of this Appendix can be straightforwardly generalized to study the effect of different types of disorder on the anisotropic nematic pairing states in doped Bi 2 Se 3 compounds [39,49] or nonunitary chiral pairing in Majorana superconductors [49,86].We leave these and related interesting questions to a future publication.

APPENDIX D: GENERALIZATION OF ANDERSON'S ARGUMENT TO A GENERIC ANTIUNITARY SYMMETRY
In the main text, we have seen that the F 0u state is protected from any disorder respecting the product of time-reversal and chiral symmetries.In this Appendix, we will show how to (trivially) generalize Anderson's original argument [1] to any antiunitary discrete symmetry that squares to −1 and acts within a two-orbital basis.To this end, we consider a discrete unitary symmetry C, which acts within the space of the two orbitals, and construct the antiunitary symmetry by multiplying it by TRS, T .We assume a basis of operators ψ τ σ , where τ = ± is the orbital and σ = ± is spin ("+" and "−" correspond to spin up and spin down, respectively), such that the symmetries act as follows: where the convention is that −α has both spin and orbit flipped compared to α and that the sign of α is given only by the spin component, such that α = 1 for (+, τ ) and α = −1 for (−, τ ).
We will now show that the transition temperature of the pairing state To show that the pairing state is not affected by disorder, we follow Anderson's original argument [1].We consider a generic dispersion Hamiltonian and a disorder Hamiltonian where, as before, the sign of the orbital indices a and α is given by the spin component alone.The first equality in Eq. (D11) stems from the definition of CT symmetry, while the second one is true because the Hamiltonian preserves this symmetry.
We are now ready for the final step.We transform the operators in Eq. (D6) to the basis that diagonalizes the sum of H 0 + H d according to Eq. (D10): where in the second line we used Eq.(D11).The final result is that this Hamiltonian has exactly same form as Eq.(D6) and the interaction weight g remains unchanged in the new basis.Consequently, T c is not modified as long as the density of the energy states na is the same as in the clean Hamiltonian.

APPENDIX E: THE CONVERSION MATRIX b Jm
The matrix elements b Jm from Eq. (11), where m = 0-15 numerates different types of disorder, equal to FIG.3.A diagrammatic representation of the decomposition from particle-hole to particle-particle channels; see Eq.(10).Matrices M, Q, F, and b are defined in Eqs.(2)-(4) and(11) [see also Tables I and II and Eq.(E1)].

) and 2 k ≡ 2
Tr F † ( k)F ( k). (C4) Matrices τ0 and τ3 here are the corresponding Pauli matrices in the Nambu space (not to be confused with the Pauli matrices in the orbital basis).The factor √ 2 in Eq. (C3) is introduced for convenience only, and simply reflects the normalization condition for functions F ( k).The form of the Gor'kov Green's function (C2) holds for the states with unitary pairing, i.e., satisfying the relation ˆ † k ˆ k ∝ 1.The nonunitary states[85], which do not satisfy this relation and can be realized in multicomponent superconductors, will be considered in future works.The matrices Q αβ m (p, k) ≡ pα|M m | kβ , describing the scattering of electrons on the impurities of type m in the MCBB basis, in the Nambu space become D1)Cψ τ σ (k) = ψ −τ σ (k), (D2)and thusCT ψ σ τ (k) = σ ψ −τ −σ (−k).(D3)For the sake of brevity, we will denote orbit and spin under the same index α = (σ, τ ), such that CT ψ kα = αψ −k−α , (D4)

= 1 2 ψ k CT ψ k = α α 2 ψ
kα ψ −k−α (D5)is not affected by disorder that respects the product CT .Notice that when C anticommutes with inversion, this state is an odd-parity pairing state.The pairing state (D5) is driven by the interaction HamiltonianH I = −g kp ψ † p (CT ) −1 ψ † p ψ k CT ψ k = −g αβ kp αβ ψ † −p−α ψ † pα ψ kβ ψ −k−β .(D6) dpk = lm e i(k−p)•r l M m .It is assumed that both Hamiltonians (D7) and (D8) respect CT .Notice that the disorder potential or dispersion Hamiltonian need not respect T or C individually, but only the product of the two.Following the Anderson argument, we now diagonalize the sum of dispersion and disorder Hamiltonians pkαβ U aα np * ĥαβ p δ kp + dαβ pk U bβ n k = n,a δ nn δ ab (D9) and the corresponding field operators c na = U aα np * ψ pα , (D10) where n denotes the spatial states that diagonalize the above Hamiltonian.Because the sum H 0 + H d possesses CT symmetry, we may assume that the new operators c na come in pairs related to each other under CT .Thus, we can denote the CT partner of c na by c −n−a (given the symmetry, we can always label states in this manner).It then follows that (CT ) −1 U aα nk CT = aα U −a−α −n,−k * = U aα nk , (D11)

kp n 1 n 2 n 3 n 4 αβabel αβ c † n 1 e U e−α n 1 −p * U lα n 2 p * c † n 2 l c n 3 b U bβ n 3 k U a−β n 4 −k c n 4 a = −g kp n 1 n 2 n 3 n 4 αβabel ae c † n 1 e U −eα −n 1 p U lα n 2 p * c † n 2 l c n 3 b U bβ n 3 k
U −aβ −n 4 k * c n 4 a = −g nn ae ae c † −n−e c † ne c n a c −n −a , (D12)

TABLE I .
Table of the impurity scattering matrices in the orbital and band (MCBB) bases appearing in Eqs. (