Topological states of multiband superconductors with interband pairing

We study the effects of interband pairing in two-band s-wave and d-wave superconductors with D4h symmetry in both time-reversal invariant as well as time-reversal symmetry breaking states. The presence of interband pairing qualitatively changes the nodal structure of the superconductor: nodes can (dis)appear, merge, and leave high-symmetry locations when interband pairing is tuned. Furthermore, in the d-wave case, we find that also the boundary modes change qualitatively when interband pairing increases: flat zero-energy Andreev bound states gap out and transition to helical edge states.

Theoretically, a two-band generalization of the Bardeed-Cooper-Schrieffer (BCS) model was introduced in Refs.[14,15].Under the assumption that the Cooper pairs are formed by the quasiparticles in the same band, i.e. intraband Cooper pairs, the order parameter in a one-dimensional (1D) pairing channel, such as s-wave or d-wave, has two components, η 1 and η 2 , which describe the pairing state in each of the two bands.The interband scattering of the Cooper pairs between the bands couples the two order parameters as η * 1 η 2 +c.c. in the lowest order within a Ginzburg-Landau (GL) expansion, analogous to the Josephson tunneling.Depending on the sign of the scattering matrix element (the coefficient of the GL coupling term), the relative phase between η 1 and η 2 in the uniform ground state is either 0 or π, corresponding to a time-reversal (TR) invariant combination.Subsequent studies have shown that the most significant qualitative differences from the single-band case are connected with the spatial and temporal variations of the relative phase, which produce such novel features as the Leggett modes [16], phase solitons [17], and fractional vortices [18] (see for a review Ref. [19]).
Recent experimental and theoretical developments have motivated a further extension of the standard theory of multiband superconductivity, by taking into account the pairing among quasiparticles from different bands, i.e. interband pairing.Within the BCS approach of pairing in the momentum space, interband pairing is * kirill.samokhin@brocku.cafeasible if the pairing interaction cutoff energy exceeds the band splitting.Alternatively, starting from a realspace pairing interaction involving different atomic orbitals in a crystalline lattice, we find interband pairing components after transformation into the band representation [20][21][22][23][24], or interband pairs arise by the proximity effect [25].Assuming that interband pairs are stabilized through a suitable microscopic mechanism, one can characterize their condensate by an additional order parameter component.Thus, a complete phenomenological description of a two-band SC involves a GL free energy which depends on three complex order parameters, two intraband ones, η 1 and η 2 , and one interband one, η.This increases the number of possible stable superconducting states, some of them breaking TR symmetry [26][27][28].
In this paper, we show how the interband pairing affects the topological properties of a two-band SC, which is manifested in a qualitative reconstruction of the energy gap of the Bogoliubov excitations.We focus on two 1D pairing channels, s-wave and d-wave, on a twodimensional (2D) square lattice and consider both TRinvariant and TR symmetry-breaking superconducting states.The gap functions corresponding to the intraband and interband pairing are introduced using a symmetrybased phenomenological approach.This approach allows one to determine the gap structure, in particular, the location of the gap nodes, even if the microscopic pairing mechanism is not known, and has proved to be very useful in the studies of unconventional fermionic superfuilds and superconductors [29,30].
According to the bulk-boundary correspondence principle, changes in the topology of the bulk state are reflected in the spectrum of the fermionic modes at the boundary [31,32].In particular, the boundary modes are expected to be different for nodeless (fully gapped) and nodal (gapless) superconducting states.These boundary modes, also known as the Andreev bound states (ABSs), have been extensively used in experimental probes to identify unconventional pairing states [33,34].In our study, we calculate the boundary mode spectrum by solving numerically the Bogoliubov-de Gennes (BdG) equations for a 2D lattice model of a two-band SC, and show how varying the strength of the interband pairing causes the system to undergo a series of topological phase tran-arXiv:2310.12002v1[cond-mat.supr-con]18 Oct 2023 sitions.
The paper is organized as follows: In Sec.II, we derive the possible interband pairing gap functions compatible with s-wave, d xy -wave, and d x 2 −y 2 -wave intraband pairing, respectively.In Sec.III, we discuss the bulk spectrum and, in particular, the movement of the gap nodes in the Brillouin zone when tuning the interband pairing.In Sec.IV, we numerically compute the edge spectrum of a d xy -wave superconductor with a strip geometry and find a topological phase transition driven by the interband pairing strength.Finally, in Sec.V, we analyse the topological phase found in Sec.IV and calculate the corresponding topological invariant(s).
Throughout the paper we use the units in which ℏ = 1, neglecting, in particular, the difference between the quasiparticle wave vector and momentum.Additionally, the lattice constant is set to unity.

II. GAP SYMMETRY: GENERAL CONSIDERATIONS
We focus on a quasi-2D centrosymmetric time-reversal (TR) invariant crystal described by the point group G = D 4h (however, our results can be straightforwardly generalized to other crystal symmetries); g ∈ G is either a proper rotation R ∈ SO(3) or an improper rotation IR, where I denotes spatial inversion.The electron Bloch states are twofold degenerate at each wave vector k = (k x , k y ) due to the combined symmetry KI, called conjugation [35].We use the index n to label the bands and an additional Kramers index s to distinguish two orthonormal conjugate states within the same band.
We further assume that only two bands n = 1, 2 cross the chemical potential and participate in superconductivity, and also that, despite the presence of the electron-lattice spin-orbit coupling, the Bloch states in both bands transform under the point-group operations and TR in the same way as the pure spin-1/2 states.Then, the conjugacy index s =↑, ↓ can be regarded as a pseudospin projection transforming under time-reversal as K |k, n ↑⟩ = |−k, n ↓⟩ and K |k, n ↓⟩ = − |−k, n ↑⟩, and we have Here, D(1/2) (R) is the spin-1/2 representation of R. In other words, we assume that both bands correspond to the double-valued irreducible representation (irrep) Γ + 6 of D 4h [36].The assumption (1), which is widely used in the theory of unconventional superconductivity [37], can be relaxed and the band symmetries corresponding to other, non-pseudospin, double-valued irreps of the point group can be considered, with important consequences for the superconducting gap structure [38,39].The superconducting system is described by the Hamiltonian where H 0 is the single-particle Hamiltonian and H sc the attractive two-particle interaction Hamiltonian within a mean-field approximation.The single-particle Hamiltonian is given by where ξ n (k) = ξ n (−k) are the band dispersions counted from the chemical potential, so that ξ 1 (k) < ξ 2 (k) for all k between the two Fermi surfaces.The superconducting mean-field pairing Hamiltonian can be represented in the following form: where are the creation operators in the time-reversed states.The intraband pairing in the nth band is described by ∆nn , whereas ∆12 and ∆21 describe the pairing of quasiparticles from different bands (the interband pairing).In order to have a non-vanishing interband pairing in a BCS-like model, one has to assume that the pairing interaction shells near the Fermi surfaces, which are defined by |ξ 1 |, |ξ 2 | ≤ ϵ c , overlap, i.e., the pairing energy cutoff ϵ c exceeds the typical band splitting E b .We do not attempt to derive the pairing Hamiltonian (3) from any microscopic model and regard the gap functions as phenomenological parameters.
Note that the gap functions ∆nn ′ (k) are defined in Eq. (3) as the measures of the pairing between the quasiparticles in the states |k, ns⟩ and K |k, n ′ s ′ ⟩, not between |k, ns⟩ and |−k, n ′ s ′ ⟩.The gap function matrices can be represented as where σ0 and σ are respectively the unit matrix and the Pauli matrices in the pseudospin space, then Eq. ( 3) takes the form Therefore, ψ nn ′ can be interpreted as the pseudospinsinglet component of the gap function and d nn ′ as the pseudospin-triplet component.It follows from the anticommutation of the fermionic operators that Additional symmetry constraints on the gap functions are obtained by looking at the transformation of the mean-field Hamiltonian (3) under the point-group operations and TR.Different pairing channels correspond to different single-valued irreps γ of G.For G = D 4h , we consider only three even-parity pairing channels: the 'swave' pairing which corresponds to the trivial irrep A 1g and the 'd-wave' pairing which corresponds to either B 1g (d x 2 −y 2 pairing) or B 2g (d xy pairing).It follows from Eq. (1) that if the pairing is described by a 1D irrep, then the gap functions satisfy the following constraints: (5) where χ γ (g) are the group characters of the irrep γ.In particular, setting g = I, we have ∆nn ′ (−k) = ∆nn ′ (k) because χ γ (I) = 1 in the even irreps.The response of the gap functions to TR is given by ∆nn ′ (k) → ∆ † n ′ n (k).Next, we introduce the order parameter components η nn ′ and represent the gap functions in the form ∆nn ′ (k) = η nn ′ φnn ′ (k).The basis functions φnn ′ which determine the momentum dependence of the gap-in particular, the location of the gap nodes-are 2 × 2 matrices in the pseudospin space, which satisfy the point-group constraint Eq. ( 5) and can have singlet and triplet components similar to Eq. ( 4).Note that φ21 (k) = σ2 φ⊤ 12 (−k)σ 2 due to the anticommutation of the fermionic operators.Regarding the constraint imposed by TR, one can prove that the basis functions can be chosen to satisfy φnn ′ (k) = φ † n ′ n (k).Denoting the intraband order parameters as η n ≡ η nn and observing that the interband gap functions ∆12 and ∆21 are not independent and characterized by the same order parameter η ≡ η 12 = η 21 , we finally obtain Here α 1 , α 2 , α and β are real even functions of k.The intraband pairing in the even channels is purely singlet; the interband pairing has both singlet and triplet components.The Pauli principle is not violated because the exchange of electrons in an interband pair involves not only the reversal of their pseudospins and momenta but also the exchange of the band indices.
Our system is characterized by three order parameter components η 1 , η 2 , and η which can be found by minimizing the Ginzburg-Landau free energy.It is easy to show that the action of TR on the order parameter components η 1 , η 2 and η is equivalent to complex conjugation, see Appendix A. One can always choose one of the components, say η, to be real and positive; then, η 1 and η 2 are either both real (positive or negative), which corresponds to a TR-invariant superconducting state, or have complex phases other than 0 or π, which corresponds to a TR symmetry-breaking superconducting state, see Appendix B.
The point-group constraints on the basis functions take the following form where R(g) denotes the rotational part of g.Below we are looking for real and even-in-k solutions of these equations, for g = C 4z and C 2y (the two rotational generators of the group D 4h ).To facilitate the numerical analysis later in the paper, the solutions are expressed in terms of the lattice-adapted basis functions of the even 1D irreps of D 4h , namely For analytical calculations, it is more convenient to use the expressions that depend only on the direction of the wave vector in the xy plane: where k = k(cos θ, sin θ).

A. s-wave pairing
For γ = A 1g , the simplest singlet solutions of the symmetry constraints (7) are given by α Collecting everything together, we arrive at the following expressions for the gap functions: where ρ is a real parameter.One can say that the interband gap functions in the A 1g channel correspond to an s + ig pairing, with the understanding that the s and g components have a different pseudospin structure.

B. d-wave pairing
The singlet components of the gap functions can be chosen in the standard form: or d xy -wave pairing, respectively.For the same reason as in the s-wave case, β1 = β2 = 0, and, using B 2g ×A 2g = B 1g and B 2g ×A 2g = B 1g , we obtain: for the d x 2 −y 2 -wave pairing and for the d xy -wave pairing.In both cases, ρ is a real parameter.
We would like to add two comments about the structure of the interband gap functions.First, the momentum dependence of their singlet and triplet components corresponds to different even irreps of the point group (for instance, in the case of d xy pairing, it is B 2g for α and B 1g for β3 ).However, the pseudospin also transforms under the point group operations, in such a way that both α and β σ3 corresponds to the same pairing channel.It is in this sense that the interband gap functions in both B 1g and B 2g channels correspond to a d+id pairing.Second, the singlet components in the interband gap functions in Eqs.(10) and (11) have the same symmetry as in the intraband ones.As seen from Eq. ( 5), this is a consequence of our assumption that both bands correspond to the same double-valued irrep of D 4h .In general, i.e., for the bands corresponding to different irreps, the symmetry of α may be different from that of α.

III. BOGOLIUBOV SPECTRUM IN THE BULK
The mean-field Hamiltonian Eq. ( 2) can be written in the form where we introduced the Nambu spinor operator and the Bogoliubov-de Gennes (BdG) Hamiltonian which is an 8 × 8 matrix in the tensor product of the band, Nambu, and pseudospin spaces.The gap functions ∆nn ′ (k) are given by Eqs. ( 9), (10), or (11).The Hamiltonian ( 14) is even in k and has the built-in particle-hole symmetry: where and τ are the the Pauli matrices in the Nambu space.
Since Û⊤ C = ÛC , the Hamiltonian is generically in the tenfold class D (Refs.[40][41][42]).The TR action on the Nambu operators Eq. ( 13) is given by KC Therefore, which is equivalent to replacing If the superconducting state is TR invariant, i.e., η 1 , η 2 , and η are all real, then the BdG Hamiltonian satisfies an additional constraint: Since Û⊤ K = − ÛK , the TR invariant BdG Hamiltonian is in the tenfold class DIII.
For the pairing symmetries we consider, the Hamiltonian Eq. ( 2) is invariant under an arbitrary U (1) pseudospin rotation c † kn↑(↓) → e ∓iθ/2 c † kn↑(↓) , and we have ĤBdG , Σ3 = 0, where Σ3 = 1 4×4 ⊗ σ3 .Therefore, Eq. ( 14) can be represented in the form where the Hamiltonians ) act in the two four-dimensional eigenspaces of Σ3 , corresponding to the two pseudospin projections.In Eq. ( 16) and everywhere below, we use the notation β = β3 and the interband order parameter η is chosen to be real and positive, but the phases of η 1 and η 2 can be arbitrary.
It follows from Eq. ( 15) that the pseudospin-resolved Hamiltonians Ĥ↑ and Ĥ↓ are transformed into each other by TR: Also, they satisfy the relation , where Û = 1 2×2 ⊗ τ2 .Introducing the magnitude and the phase of the interband gap functions with | ∆(k)| = ηg(k) and g = α2 + β2 , one can see that Ĥ↑ and Ĥ↓ are particle-hole symmetric at each k, in the following sense: where Therefore, Ĥ↑ and Ĥ↓ have the same bulk spectrum, which consists of symmetric pairs of eigenstates E and , the pseudospin-resolved Hamiltonians Ĥ↑ and Ĥ↓ are in the tenfold class C.
The matrices Ĥ↑ and Ĥ↓ can be diagonalized analytically, see Appendix C, and we find that the bulk Bogoliubov spectrum consists of four branches ±E ± , where The notations are as follows: are the intraband gap functions and One can show that P > Q in the presence of interband pairing.Therefore, E + is strictly greater than E − at all k.Each of the four branches ±E ± is twofold degenerate due to pseudospin.In the absence of interband pairing, we set η = 0 and recover the usual expressions for a twoband superconductor: where is the excitation energy in the nth band.While the upper Bogoliubov excitation branch E + is fully gapped in the superconducting state, the lower branch E − vanishes at the wave vector k if in which case E − (k) and −E − (k) touch, producing a gap node.In two spatial dimensions, the three real functions r 1,2,3 cannot simultaneously vanish at the same k, unless forced to do so by additional symmetries.
Writing the intraband order parameters in the form we see that r 3 identically vanishes in the states in which φ 1 +φ 2 = 0 or 2π.This happens, in particular, in the TR invariant states in which η 1 and η 2 are both either real positive or real negative.As shown in Appendix B, the TR-symmetry breaking states with φ 1 + φ 2 = 0 or 2π are stable only if the system's parameters are fine-tuned, the possibility that can be neglected.In a generic state with the interband pairing, i.e., when η ̸ = 0 and In the s-wave case, this can only happen accidentally and is neglected.In contrast, the d-wave intraband gap functions, and therefore r 3 , vanish along the high-symmetry directions for symmetry reasons.Thus, there exist four classes of the stable bulk nodal structures, which are studied below.

A. Generic s-wave pairing
In this case, the phases of η 1 and η 2 take any values, except φ 1 = φ 2 = 0 or π.The TR invariant states in which η 1 and η 2 are real but have opposite signs are also included here.Since r 3 is nonzero at all k, the s-wave superconducting state is fully gapped, regardless of the strength of the interband pairing.

B. Generic d-wave pairing
For concreteness, let us consider the evolution of the nodal structure in the d xy -wave case (for the d x 2 −y 2 -wave pairing, the nodes are just rotated by π/4).In the absence of the interband pairing, r 3 identically vanishes and the point gap nodes are located where ξ 1 = ψ 1 = 0 or ξ 2 = ψ 2 = 0, i.e., at the intersections of the axes of the 2D Brillouin zone with the Fermi surfaces.
In the presence of the interband pairing and for generic phases of η 1 and η 2 , r 3 only vanishes along the axes k x = 0 and k y = 0 for symmetry reasons.Moreover, r 2 is also zero there and the only remaining gap node condition, Eq. ( 20), takes the following form: along the k x = 0 or k y = 0 lines.As η increases, the nodes remain on the high-symmetry axes, but move into the "interband space", where ξ 1 < 0 and ξ 2 > 0 (recall that we assume ξ 1 < ξ 2 ).Eventually, at a sufficiently strong interband pairing, the nodes merge and annihilate each other, which marks the transition into a fully gapped bulk phase, as shown in Fig. 1.Annihilating nodes were also found in Ref. [43], in a model of a TRinvariant SC with a d x 2 −y 2 -wave intraband pairing and d xy -wave interband pairing.The evolution of the gap structure can be studied analytically using a simple model with two parabolic electron-like bands: where E b > 0 is the band splitting, in which the two Fermi surfaces are circles of radii k F,1(2) = k 0 1 ± mE b /k 2 0 .We use the following gap symmetry factors: see Eqs. (11) and (8).The gap node equation Eq. ( 22) becomes along the axes of the momentum space.Taking, for instance, the θ = 0 axis, at η = 0 the two nodes are located on the Fermi surfaces, at k = k F,1 x and k = k F,2 x.As η increases, the nodes move towards each other, to k = k 1 x and k = k 2 x, where Figure 2. Schematic illustration of the nodal behavior in the toy model of the exceptional s-wave case (Sec.III C).Empty red dots: stray nodes at interband pairing η = ηc,1, see Eq. ( 29).Filled red dots: stray nodes at interband pairing ηc,1 < η < ηc,2, see Eq. (30).Red crosses: annihilation of stray nodes at interband pairing η = ηc,2.
Finally, when the interband order parameter η reaches the critical value the nodes merge at k = k 0 x and "annihilate" each other.At stronger interband pairing, our d xy -wave superconductor is fully gapped.Note that the disappearance of the nodes happens only if ρ ̸ = 0, i.e., when the interband gap functions contain the triplet component with the d x 2 −y 2 -wave-like momentum dependence.
C. s-wave pairing, φ1 = φ2 = 0 or π In this TR invariant state, r 3 is identically zero everywhere, but we still need to solve the remaining equations r 1 = 0 and r 2 = 0. Assuming, without loss of generality, the same intraband symmetry factors in both bands: One way to satisfy this is to put α = 0, but that does not happen in the s-wave case, whereas in the d-wave case that can happen only along the high-symmetry lines, which was already considered in Sec.III B. Therefore, if we look for the gap nodes away from the high-symmetry axes, then we need to solve the following two equations: Note that the second equation can have solutions only between the two Fermi surfaces, where ξ 1 < 0 and ξ 2 > 0, i.e., ξ If the interband pairing is sufficiently strong, so that the right-hand side of the first equation is positive, then Eq. ( 27) defines two lines between the Fermi surfaces.The intersections of these lines, if they exist, correspond to accidental point nodes in the excitation spectrum.
To illustrate these points for the s-wave pairing, we use the band structure model Eq. ( 23), with the following angular dependence of the gap functions: see Eqs. ( 9) and (8).Solving Eq. ( 27), we obtain that the accidental nodes are located on the circle of the radius where at the angles found from the equation In the absence of the interband pairing, this last equation does not have any solutions.As η increases and reaches the nodes emerge in pairs, first at the angles given by θ = π/8, 3π/8, . . . .As the interband pairing further increases, the nodes split and move along the circle defined by Eq. ( 28) towards the angles θ = 0, π/4, . . ., where they finally merge and disappear at These changes in the gap structure are shown in Fig. 2. The superconducting state is fully gapped at η < ηc,1 and at η > ηc,2 .Note that the gap nodes appear in this model only if ρ ̸ = 0, i.e., when the interband gap functions contain the "triplet" component with an anisotropic, g-wave-like momentum dependence.Although these nodes are topologically unstable, since any deviation from the condition φ 1 = φ 2 = 0 or π will remove them, they are protected by TR symmetry.The difference from the generic d-wave case is that r 3 now vanishes everywhere, which makes it possible for additional gap nodes to appear away from the highsymmetry lines.Repeating the reasoning from Sec. III C, we find that there are two types of nodes: the "highsymmetry" ones, which are located where α(k) = 0, and (Top) Empty red dots: high-symmetry nodes without interband pairing η = 0. Filled red dots: high-symmetry nodes at small interband pairing 0 < η < ηc,1, see Eq. (32).(Bottom) Empty red dots: high-symmetry nodes at interband pairing η ≲ ηc,1.Filled red dots: high-symmetry nodes (along the main axes) and stray nodes (off the main axes) at interband pairing ηc,1 < η < ηc,2, see Eq. (33).Red crosses: annihilation of high-symmetry nodes (along the main axes) at interband pairing η = ηc and annihilation of stray nodes (off the main axes) at interband pairing η = ηc,2 (ηc < ηc,2).also the "stray" ones, which correspond to the solutions of Eq. (27).
To develop some analytical insight, we again assume a d xy -wave pairing and use the parabolic bands (23), with the symmetry factors given by Eq. (24).We obtain that the stray nodes are located on the circle defined by Eq. ( 28), at the angles determined by the equation At η = 0, this equation has no solutions.To illustrate the different scenarios of how the stray nodes are created and destroyed by increasing the interband pairing strength, we solve Eq. ( 31) in three cases, for ρ = 0, |ρ| ≫ 1, and ρ = 1.At ρ = 0, which corresponds to the absence of the triplet d x 2 −y 2 component in the interband gap functions, the stray nodes appear in pairs at θ = π/4, 3π/4, . . .when η reaches the critical value As the interband pairing strength further increases, the nodes split and move away from each other, staying on the circle (28) and asymptotically approaching the axes θ = 0, π/2, . . . .Since at ρ = 0 both the intraband and interband gap functions vanish along the axes, the highsymmetry nodes are not affected by η, i.e., always remain at the intersections of the two Fermi surfaces with the lines k x = 0 and k y = 0.
To see what happens in the general case, when both the singlet (d xy ) and the triplet (d x 2 −y 2 ) components are present in the interband gap functions, we set ρ = 1.Then, the solutions of Eq. ( 31) exist only if ηc,1 ≤ η ≤ ηc,2 , where and As η increases, the stray nodes first appear on the axes, i.e., at θ = 0, π/2, . . ., where they peel off in pairs from the high-symmetry nodes, then move towards θ = π/4, 3π/4, . . ., where they eventually merge and annihilate each other.Note that ηc,1 is less than the critical strength of the interband pairing at which the high-symmetry nodes disappear, see Eq. ( 26).Therefore, there is an interval of η, in which the stray nodes co-exist with the high-symmetry ones, so that there are sixteen nodes altogether (eight of each type), all located between the Fermi surfaces, as shown in Fig. 3. Similar behaviour of the nodes was also found in a different model in Ref. [44], in which the interband pairing in a TR invariant d-wave state is controlled by the inter-orbital SO coupling.In contrast to the limiting cases of ρ = 0 and ρ ≫ 1, in which ηc,2 = ∞ and the stray nodes survive the strong interband pairing, in the general case all nodes eventually disappear as η increases.
The stray nodes are accidental, in the sense that they are not protected by the crystal symmetry.In order to destroy them, one has to tune the intraband order parameter phases out of the TR invariance condition φ 1 = φ 2 = 0 or π.However, this condition always corresponds to a critical point of the free energy, see Appendix B, and if this critical point is a minimum, then the state is stable.Therefore, the stray nodes are protected by the TR symmetry.

E. Summary
The effect of the interband pairing on the energy gap nodes in the bulk is fundamentally different in the four cases discussed in this section.The least interesting case is the generic s-wave state, which is fully gapped at η = 0 and remains so as η increases.This state is topologically trivial and does not support zero-energy boundary modes.
In the TR-invariant s-wave state with φ 1 = φ 2 = 0 or π, a sufficiently strong interband pairing can create and then destroy again point nodes between the Fermi surfaces.The critical values of η separating the gapped and gapless phases, as well as the locations of the nodes, are model-dependent.
In the generic d-wave state, the nodes are located only along the high-symmetry lines.As η increases, these nodes leave the Fermi surfaces, move towards each other into the interband space, and merge and disappear at η = ηc .The gapless and gapped phases separated by ηc are expected to be topologically different, which is confirmed in Secs.IV and V below.
The TR-invariant d-wave state with φ 1 = φ 2 = 0 or π exhibits the most complex behaviour.In this case, the nodes along the high-symmetry lines can co-exist with the additional (stray) nodes in the interband space, whose number and locations are model-dependent.As the interband pairing increases, the two families of nodes evolve as the system passes through a series of transitions characterized by the creation and destruction of the pairs of nodes.Across these transitions, the topological charges of the nodes are conserved, see Sec.V A. Eventually, at a sufficiently large η, all the nodes will have pairwise collided and annihilated each other, so that the superconducting state will be fully gapped.

IV. BOUNDARY MODES
Having discussed the bulk properties of the different superconducting phases in Sec.III, we turn our attention towards the boundaries of the material.We consider a strip geometry: the system is infinitely extending along  the x direction, but has a finite width along the y direction.To formulate the lattice model, we assume N x (N y ) lattice sites with periodic (open) boundary conditions along the x (y) direction.As a result, the system effectively possesses two infinitely extended edges parallel to the x axis.
The bulk band structure is described by where µ is the chemical potential, t n is the nearest, and t ′ n -the next-nearest neighbor hopping amplitude in the nth band.Superconductivity is either of sor d xy -wave type, given by Eq. ( 9) and Eq. ( 11), respectively; for d x 2 −y 2 -wave superconductivity, we would consider edges rotated by π/4.The total bulk Hamiltonian of the lattice model is described by Eq. ( 2).
The formal description of the superconducting strip system follows Ref. [45].The momentum component k y is not a good quantum number, because translation sym- metry is broken along the y direction.To account for this, we only consider the momentum representation k ≡ k x along the x direction, while we keep the real space representation i ≡ i y along the y direction.Assuming a sufficiently wide strip, the superconducting gap is approximately constant along the y direction and we neglect any potential surface effects causing the order parameter to be spatially deformed close to the edges.The order parameters η 1 , η 2 , and η are not computed self-consistently, but set to their respective bulk values.
The BdG Hamiltonian ĤBdG is diagonal in the pseudospin space, with the blocks Ĥ↑ and Ĥ↓ being related by TR, see Eq. ( 17).Furthermore, it is diagonal in k space.As a consequence, the problem is reduced to the diagonalization of a 4N y × 4N y matrix.Here, we employ an exact diagonalization procedure to solve for the eigenvalues as well as the corresponding eigenstates of Ĥ↑ (k).
We begin the discussion of our results with the s-wave case.For a weak interband pairing, the spectrum is fully gapped in both the generic (cf.Sec.III A) and the exceptional cases (cf.Sec.III C).In the generic case, this situation remains true regardless of the strength of the interband pairing.In contrast, in the exceptional case, there exist two critical values, ηc,1 and ηc,2 , between which the spectrum is gapless for four distinct k values (three k values at η = ηc,1 and two k values at η = ηc,2 ).As the interband pairing strength increases from ηc,1 to ηc,2 , these nodes move towards each other until they annihilate eventually, so that the spectrum is fully gapped again.There are no edge states present-regardless of the interband pairing strength.
Next, we discuss the generic d xy -wave case (cf.Sec.III B).The results are summarized in Fig. 4. In the absence of interband pairing (the top panels), the energy spectrum shows five zeros at These correspond to the bulk nodes along the main axes and are located exactly on the two Fermi surfaces.The node at k 0 is four-fold degenerate corresponding to the four nodes of the system along the y direction.Between k LO and k RO the spectrum shows flat ABSs.As soon as the interband pairing is turned on, the ABSs between k LI and k RI gap out, while they remain intact between k LO and k LI as well as between k RI and k RO .
As the interband pairing strength increases, the nodes move away from the Fermi surfaces until they meet each other and annihilate.After this point, the bulk is completely gapped but eight zero-energy crossing ABS branches remain.They are singly degenerate (doubly degenerate for Ĥ↑ ⊕ Ĥ↓ if both edges are taken into account, see Fig. 5) and the corresponding eigenstates are localized near the edges of the strip, see Appendix D. They mark a different, topologically non-trivial, superconducting phase.
The edge states are schematically illustrated in Fig. 5.For Ĥ↑ , four states are located close to the left edge of the strip, while the other four are located close to the right edge of the strip.Depending on their slope they move either along the positive or negative x direction.Furthermore, they mix electrons from one band with holes from the other band, see Appendix D.
Finally, we turn our attention towards the exceptional d xy -wave case (cf.Sec.III D).Similarly to the generic case, we observe a topological phase transition when the interband pairing strength increases.However, apart from the phase transition, the behavior of the spectrum significantly differs from the generic case.Weak interband pairing does not immediately and fully gap out the ABSs.Indeed, there are ABSs present until the stray nodes annihilate along the diagonals of the Brillouin zone (cf.Sec.III D).
We make two final remarks.First, the topological phase transition occurs regardless of the presence of TR symmetry.The characteristics of the edge states in the strip spectrum (cf.Fig. 4) are the same regardless of the intraband order parameter phases φ n .Only the behavior of the bulk nodes and the ABSs in the gapless regime differs between the TR symmetry-breaking generic states and the TR invariant exceptional states.Second, the interband pairing strength required to reach the topologically non-trivial superconducting phase strongly depends on the distance between the two Fermi surfaces.The closer the Fermi surfaces along the main axes of the Brillouin zone, the weaker the required interband pairing in order to annihilate the gap nodes.

V. TOPOLOGICAL ARGUMENTS
The results of the previous two sections show that the effects of the interband pairing are most profound in the d-wave states.The evolution of the bulk gap structure, which is reflected in the changes of the ABS spectrum, can be interpreted in terms of a series of transitions between topologically distinct superconducting phases.In this section, we discuss the relevant bulk topological invariants, focusing as before on the d xy -wave states.

A. Gapless bulk
According to Sec.III, the gap structure of a generic nodal d xy -wave state is insensitive to the phases of the intraband order parameters: the nodes move on the highsymmetry axes as the interband pairing strength varies.In order to study the topological properties, one can focus on the TR invariant states, in which the order parameter components η 1 , η 2 , and η are all real.If η 1 η 2 < 0, then there exist only the high-symmetry nodes (Sec.III B), whereas at η 1 η 2 > 0 the stray nodes are also possible (Sec.III D).
For the real order parameters, the Hamiltonians Ĥ↑ and Ĥ↓ , see Eq. ( 16), have a 'chiral' symmetry: In the basis in which ÛS is diagonal, the Hamiltonians can be brought to a block off-diagonal form, e.g.,  where The υ-matrix for Ĥ↓ is obtained from Eq. ( 36) by replacing φ → − φ.
The positions of the gap nodes are determined by the zeros of |det υ|, whereas the topological charges of the nodes are given by the winding number of the phase of det υ: see Refs.[46,47].The integration here is performed around an infinitesimally small circular contour wrapping counterclockwise around the node.From Eq. ( 36), we have where ψ n (k) = η n α(k), assuming the same intraband symmetry factors in both bands.In agreement with the results of Sec.III, we see that the zeros of | det υ| are located in the interband space, either where α = 0 and ξ 1 ξ 2 + ∆2 = 0 (the high-symmetry nodes), or away from the symmetry axes (the stray nodes), the latter being possible only if η 1 and η 2 have the same sign.The topological charges of the nodes can be easily calculated by expanding det υ(k) in the vicinity of the nodes.In the case η 1 > 0, η 2 < 0, we find (see Appendix E) that the gap nodes located on the same axis are oppositely "charged", as shown in Figs. 6, 7, and 8, which makes it possible for the nodes to "annihilate" each other, as discussed in Sec.III B.
According to Refs.[46,47], the number of the zeroenergy edge modes at given momentum k x along the boundary, per one pseudospin projection, is equal to |N (k x )|, where ∇ ky ln det υ(k).
The integral here is taken along a straight line which runs between the opposite edges of the Brillouin zone perpendicular to the boundary.In a continuum model, the limits are extended to infinity.Assuming that all gap functions vanish far from the Fermi surfaces, one can integrate along a closed contour C shown in Figs. 6,  7, and 8. Using Stokes' theorem to contract the contour without crossing any gap nodes, we find that N (k x ) is equal to the total charge of the nodes enclosed by C.
In this way, we obtain: where k 1,2 are the positions of the bulk nodes, see Eq. (25).Taking into account the pseudospin degeneracy, the total number of the zero-energy ABS localized near one edge of the sample is equal to 2|N (k x )|.We see that the momentum range in which the topologically protected zero-energy boundary modes exist shrinks with increasing the interband pairing and eventually disappears, in agreement with the numerical results of Sec.IV.

B. Gapped bulk
We have seen in the previous sections that the bulk Bogoliubov spectrum becomes fully gapped when the interband pairing exceeds certain value.Moreover, if ∆ ̸ = 0, then the Bogoliubov branches E + and E − , see Eq. (19), are always separated.In the absence of any level crossings, the intraband pairing can be adiabatically turned off without affecting the bulk topology.Therefore, in order to study the topology of the mappings k → Ĥ↑ (k) and k → Ĥ↓ (k) in the nodeless regime, we can set η 1 = η 2 = 0, which considerably simplifies the calculations.The pseudospin-resolved Hamiltonians (16) are then reduced to direct sums of 2 × 2 matrices: where and ν ± = (η α, ∓η β, ξ) and ξ = (ξ 1 + ξ 2 )/2.Since the Hamiltonians Ĥ↑ and Ĥ↓ are in the class C, the gapped bulk states are characterized by an even (2Z) topological invariant.[42] Therefore, we expect an even number of zero-energy boundary modes.
In the case of the d xy -wave pairing, we have see Eqs. (11) and (8).Therefore, the matrices ( 40) have the same form as the BdG Hamiltonians for the chiral d ± id states, shifted up or down in energy.It is well known [48] that the d + id and d − id superconductors can support chiral boundary modes, which are protected by the bulk topology.Diagonalizing Eq. ( 40) and using the band model ( 23), we find that the bulk spectra of Ĥ↑ and Ĥ↓ are the same and, in agreement with Eq. ( 19), are given by four particle-hole symmetric branches ±E ± , where and ∆ has the form (18).Note that the branch indices in this last expression have nothing to do with the "chirality" index ± in Eq. (40).At a sufficiently strong interband pairing, | ∆| > E b /2, the bulk spectrum is fully gapped.According to Ref. [31], the topological invariant characterizing a gapped chiral d-wave state has the following form: where ν = ν/|ν| and ν = ν + or ν − .Here we integrate over the 2D momentum space, which can be compactified into an S 2 sphere, because the gap functions vanish outside the overlapping BCS pairing shells, see Sec.II, so that ν = ẑ sign ξ and the integrand in Eq. ( 42) is equal to zero far from the Fermi surfaces.The expression ( 42) is nothing but the degree of the mapping k → ν(k), which takes integer values and can be used to enumerate different equivalence classes of the Hamiltonians Eq. (40).
Writing the interband gap functions in the form Eq. ( 18), with | ∆| nonvanishing only inside the pairing shells of thickness ϵ c , Eq. ( 42) takes the form Finally, neglecting the ξ-dependence of | ∆| and φ inside the pairing shell, sending ϵ c → ∞, and integrating with respect to ξ, we obtain: where the integration is performed along the ξ = 0 line.Since we neglect the ξ-dependence of φ, one could integrate along either of the two Fermi surfaces, with the same result (note that the lines ξ = 0, ξ 1 = 0, and ξ 2 = 0 all lie within the BCS pairing shell, which encompasses both Fermi surfaces).Thus, the invariant ( 42) is equal to the phase winding number of the interband gap function.
For the d xy -wave pairing, see Eq. ( 41), we obtain the following winding numbers for ĥ± ↑ and ĥ± ↓ : Therefore, each of the four Hamiltonians ĥ± ↑ and ĥ± ↓ has two chiral zero modes near each edge of the sample.These modes have opposite slopes for opposite chiralities, and are also shifted up and down in energy by ±E b /2, as shown in Fig. 9 for ĥ+ ↑ and ρ > 0. The 4 × 4 Hamiltonian Ĥs corresponding to one pseudospin channel has four helical modes composed of two pairs of the counterpropagating chiral modes from ĥ+ s and ĥ− s , as shown in Fig. 5.This result is in agreement with the numerical solution of the BdG equations, see Fig. 4.

VI. CONCLUSION
Based on a symmetry analysis, we determined the possible interband pairing gap functions in the case of twoband superconductors with s-wave, d xy -wave, or d x 2 −y 2wave pairing, which can give rise to both TR invariant and TR symmetry-breaking superconducting states.As the interband pairing strength increases, the nodal structure changes fundamentally.Nodes leave the Fermi surfaces and eventually annihilate each other on the highsymmetry axes, whereas other nodes (stray nodes) appear, move, and merge in the interband space.
In the case of a d-wave superconductor with a strip geometry, the boundary modes exhibit qualitative changes when interband pairing increases.Starting from zeroenergy flat ABSs in the absence of interband pairing, these modes partially gap out as soon as interband pairing is turned on.In the limit of strong interband pairing, the system undergoes a topological phase transition to a fully gapped helical d ± id-wave superconducting state.The corresponding topological invariant is the phase winding number of the interband gap, which explains the existence of the eight gap-crossing zero-energy branches near one edge of the sample in the helical state.where ν n = (Re ψ n , − Im ψ n , ξ n ).It is manifestly particle-hole symmetric and one can find its spectrum either by a direct calculation of a 4 × 4 determinant [39] or by using the following trick [50].
Let us calculate the second and fourth powers of Eq. (C1): where One can see that the matrix M = Ĥ4 − (µ 1 + µ 2 ) Ĥ2 does not contain off-diagonal 2 × 2 blocks.Moreover, since we find that M is proportional to the 4 × 4 unit matrix.Therefore, the eigenvalues of Ĥ satisfy the following biquadratic equation: Solving it, we obtain the Bogoliubov energy branches given by Eq. (19).
If periodic boundary conditions are also considered along the y direction (i.e., the strip is closed to a torus), then additional off-diagonal terms appear in the corners of ξ↑(↓) and ∆↑(↓) .Fig. 10 illustrates that the gap-crossing energy branches in the strong interband-pairing regime belong to states which are localized near the edges of the strip.The respective eigenstates ψ of the Hamiltonian Ĥ↑ are of size 4N y and split into electron and hole contributions as well as band contributions n = 1 and n = 2, as shown in Fig. 11.Let us consider, for example, the high-symmetry nodes along the positive k x axis, i.e., at θ = 0.They are located between the two Fermi surfaces, at ξ = ±ξ 0 , where We expand Eq. (E1) near the nodes by setting ξ = ±ξ 0 + ξ 0 x, θ = y (|x|, |y| ≪ 1), and obtain det υ = ∓ξ 2 0 (x + iw ± y), where Therefore, the topological charges of the nodes, see Eq. ( 37), are given by q ± = sign (w ± ).
In particular, in the absence of the interband pairing, we have q + = − sign (η 1 ) and q − = − sign (η 2 ).It follows from Eq. (E2) that If η 1 and η 2 have opposite signs, then q + q − < 0, independently of the value of η.The high-symmetry nodes on the same axis have opposite charges and annihilate each other at the critical strength of the interband pairing, given by Eq. ( 26).
In contrast, if η 1 and η 2 have the same sign, then Eq. (E3) changes sign at At this value of the interband pairing, one of the highsymmetry nodes splits into two stray nodes of the same charge and one high-symmetry node of the opposite charge, see Sec.III D. All nodes eventually annihilate each other at a sufficiently strong interband pairing.

Figure 5 .
Figure 5. Schematic illustration of the edge states in a fully gapped dxy-wave SC, for strong interband pairing (top panel: the edge states for Ĥ↑ , bottom panel: the edge states for Ĥ↓ ).

Figure 9 .
Figure 9. Schematic illustration of the chiral edge modes for the Hamiltonian ĥ+ ↑ , see Eq. (40).Dashed lines: chiral modes without the energy shift.Solid lines: chiral modes shifted by −E b /2.

k 4 holes, k 1 holes, k 2 holes, k 3 Figure 11 . 2 2−
Figure 11.Eigenstate profiles of the states marked by red dots in Fig. 10 as a function of y position.Left column: electron components of the eigenstates.Right column: hole components of the eigenstates.