Topological properties of finite-size heterostructures of magnetic topological insulators and superconductors
Julian Legendre Eduárd Zsurka Daniele Di Miceli Llorenç Serra Kristof Moors Thomas L. Schmidt






Department of Physics and Materials Science, University of Luxembourg, 1511 Luxembourg, Luxembourg


Department of Physics and Materials Science, University of Luxembourg, 1511 Luxembourg, Luxembourg
Peter Grünberg Institute (PGI-9), Forschungszentrum Jülich, 52425 Jülich, Germany
JARA-Fundamentals of Future Information Technology, Jülich-Aachen Research Alliance, Forschungszentrum Jülich and RWTH Aachen University, Germany


Institute for Cross-Disciplinary Physics and Complex Systems IFISC (CSIC-UIB), E-07122 Palma, Spain
Department of Physics and Materials Science, University of Luxembourg, 1511 Luxembourg, Luxembourg


Institute for Cross-Disciplinary Physics and Complex Systems IFISC (CSIC-UIB), E-07122 Palma, Spain
Department of Physics, University of the Balearic Islands, E-07122 Palma, Spain


Peter Grünberg Institute (PGI-9), Forschungszentrum Jülich, 52425 Jülich, Germany
JARA-Fundamentals of Future Information Technology, Jülich-Aachen Research Alliance, Forschungszentrum Jülich and RWTH Aachen University, Germany


Department of Physics and Materials Science, University of Luxembourg, 1511 Luxembourg, Luxembourg



Heterostructures of magnetic topological insulators (MTIs) and superconductors (SCs) in two-dimensional (2D) slab and one-dimensional (1D) nanoribbon geometries have been predicted to host, respectively, chiral Majorana edge states (CMESs) and Majorana bound states (MBSs). We study the topological properties of such MTI/SC heterostructures upon variation of the geometry from wide slabs to quasi-1D nanoribbon systems and as a function of the chemical potential, the magnetic doping, and the induced superconducting pairing potential. To do so, we construct effective symmetry-constrained low-energy Hamiltonians accounting for the real-space confinement. For a nanoribbon geometry with finite width and length, we observe different phases characterized by CMESs, MBSs, as well as coexisting CMESs and MBSs, as the chemical potential, the magnetic doping and/or the width are varied.


Topological properties of finite-size heterostructures of magnetic topological insulators and superconductors
    Thomas L. Schmidt
    October 2023
=============================================================================================================







§ INTRODUCTION


Topological superconductors are fascinating phases of matter which have stirred significant interest in the scientific community <cit.>. These phases display gapped bulk states with superconducting pairing as well as topologically protected gapless surface states, which have been predicted to be Majorana states. The search for these quasi-particles has stimulated an intense research activity, primarily owing to their potential for quantum computing <cit.>. Nevertheless, proposed realizations of topological superconductors with large enough bulk gaps remain rare and Majorana states remain elusive and controversial.

Bringing together superconducting pairing and spin-orbit (SO) or SO-like interactions is a promising avenue for creating topological superconductivity. It has for instance been studied in superconductors with strong SO interactions <cit.>, in topological materials where doping with Nb, Sr or Cu yields a superconducting gap in the bulk <cit.>, or in heterostructures combining strong SO interactions or SO-like interactions induced by a magnetic texture with a conventional superconductor <cit.>. Breaking time-reversal symmetry (TRS) is also often necessary for the emergence of low-dimensional surface states such as chiral edge states or bound states. 






TRS can also be broken by an external magnetic field, but this may not be compatible with superconductivity. It is thus desirable to explore intrinsic magnetism (e.g. via magnetic doping) of heterostructures as an alternative <cit.>. In our work, we study MTI/SC heterostructures where the interplay between superconducting pairing, spin-orbit coupling, and TRS breaking leads to the appearance of topological superconducting phases. Our study will apply to heterostructures consisting of s-wave SCs and magnetically doped compounds of the Bi_2Se_3 family. The effective realization of such heterostructures has recently shown promising progress <cit.>. We will consider thin MTI films, which have become experimentally realizable over the past decade <cit.>.
A comprehensive description of MTI thin films is achieved by the construction of a symmetry-constrained k·p Hamiltonian <cit.>. Here, we introduce this model as a basis for the subsequent finite-size calculations and to relate our study to concrete material systems. We discuss its applicability for the system we investigate and we comment on the limits of such a model <cit.>.

In recent works, planar translation-invariant geometries and confined quasi-1D geometries of MTI/SC heterostructures have been suggested as a potential platform for Majorana physics. Such systems have been studied for specific values of the chemical potential μ, of the strength λ of the magnetic exchange interaction with the magnetic dopants, and for specific sizes <cit.>. Here, we provide a more comprehensive treatment, studying a wider region of the phase diagram in (μ,λ) space and the nature of the topological states when the in-plane size of the heterostructure is varied from infinitely large to finite. 
Moreover, we discuss the transition from chiral Majorana edge states (CMES) to Majorana bound states (MBS) and their respective localization. In the following, we will refer to planar translation-invariant geometries as “slab geometries" and to in-plane confined geometries as “nanoribbons". For nanoribbons with intermediate width, and depending on the magnitude of the magnetic exchange term in the MTI, we observe regions in the phase diagram where the low-energy sector hosts coexisting CMESs and MBSs. Our understanding of the topological properties associated with these finite-size systems, from two-dimensional to quasi-one-dimensional, is based on symmetry-constrained analytical low-energy models that we construct throughout this paper. These characterize the appearance of gapless surface (edge or end) states through the bulk-edge correspondence. 

This paper is organized as follows: In Sec. <ref>, we introduce the effective models we use for the description of the MTI. First, we review the construction of a symmetry-constrained k·p Hamiltonian which characterizes a 3-dimensional (3D) topological insulator (TI). Then we construct effective models for a thin MTI slab system and for a nanoribbon geometry by considering the low-energy states arising from quantum confinement in the 3D model. We study the occurence of chiral edge modes as a function of the magnetic exchange term and as a function of the width of the nanoribbon.
In Sec. <ref>, we study the topological properties of an MTI/SC slab and then we investigate the occurrence of CMES and MBS in MTI/SC with nanoribbon geometries.




§ EFFECTIVE MODEL FOR THE MTI
 



 §.§ Symmetry-constrained k dot p Hamiltonian


A convenient way of describing the topological properties of a (3D) compound of the Bi_2 Se_3 family is via the following Hamiltonian <cit.>,

    H_TI^3D( k)
       = H_0(k_z) + H_1(k_x,k_y), 
    
    H_0(k_z)    = D_1 k_z^2  𝕀 + (M-B_1 k_z^2 ) τ_z + A_1 k_z τ_x σ_z, 
     
    H_1(k_x,k_y)
       = ( C +D_2k^2)  𝕀 -B_2k^2 τ_z 
       + A_2 τ_x ( k_x σ_x + k_y σ_y ),

where k^2 ≡ k_x^2+k_y^2, k_±≡ k_x± i k_y, and C, M, D_1, D_2, B_1, B_2, A_1, A_2, are real coefficients. This Hamiltonian acts in the (atomic) low-energy basis of states {|+,↑⟩,|-,↑⟩,|+,↓⟩,|-,↓⟩} through the Pauli matrices σ_x,y,z and τ_x,y,z which act, respectively, in the spin {↑,↓} and parity {+,-} spaces. The parity operation, equivalent to inversion symmetry in 3D, acts on {|+,↑⟩,|-,↑⟩,|+,↓⟩,|-,↓⟩} according to the matrix representation τ^z and flips the momentum vector, leaving the Hamiltonian invariant: τ^z H_TI^3D(-  k) τ^z = H_TI^3D( k).
This Hamiltonian is valid in the vicinity of the Γ point, where it describes the low-energy properties of the system. In Fig. <ref>, we show the eigenenergies of this Hamiltonian.

This Hamiltonian relies on a long-wavelength approximation and is a priori valid only for finite systems above a certain size. The coefficients appearing in the Hamiltonian can be determined by fitting the energy spectrum to experimental data or extracting the parameters from ab-initio calculations (see Ref. <cit.>). The distance from the Γ point up to which this model Hamiltonian is a good approximation allows an estimate of the possible size of the system.
The band structure of the Hamiltonian (<ref>) is in good agreement with ab-initio calculations for |k_x| ≲ 0.02 ^-1, |k_y| ≲ 0.02 ^-1, and |k_z| ≲ 0.05 ^-1 <cit.>. This implies that the system thickness d in the z direction should be larger than 20  while the width W and the length L, respectively, in the y direction and x-directions should be larger than 50 .

It is also important to note that H_TI^3D( k) is constrained by the symmetries of the system, specifically TRS, inversion symmetry, and the three-fold rotation symmetry around the z axis. The atomic states used to construct the Hamiltonian are based on quintuple layer unit cells which respect the aforementioned symmetries <cit.>. For the materials we consider, the thickness (along the z direction) of a quintuple layer is approximately 1 nm. In the following, we will consider systems with widths W and lengths L larger than 40 nm, which is much larger than the length of a unit cell in the x and y-directions (below 1 nm). This allows us to study a continuum of lengths above 40 nm along these directions. We notice that W ∼ 40 nm or L ∼ 40 nm is a reasonable lower limit of what is currently achieved experimentally. We will consider slabs and nanoribbons of Bi_2 Se_3 family compounds consisting of two or more quintuple layers (thickness d larger than 20 ).









 §.§ Slab geometry
 
Let us consider a TI slab (a planar translation-invariant geometry) with a finite thickness d along the z direction, with bottom and top surfaces of the slab located, respectively, at z=0 and z=d. We replace k_z → -i ∂_z, impose vanishing wave functions as boundary conditions at z= 0 and z= d, and denote by E_+ and E_- the two lowest eigenenergies of H_0(-i∂_z), each of which is twofold degenerate because of TRS. Magnetic doping can be accounted for by adding a TRS-breaking Zeeman term λσ_z; the system is now an MTI. The respective eigenstates of the Hamiltonian will be denoted by |φ^σ⟩ and |χ^σ⟩, where σ∈{↑,↓}.
Projecting H_TI^3D( k) on the basis {|φ^↑⟩, |χ^↑⟩, |φ^↓⟩, |χ^↓⟩} gives rise to a 4 × 4 Hamiltonian which describes an MTI slab <cit.> (we use ħ = 1 in the following),

    H_MTI^2D( k)    =
    [   h_+( k)         0;         0 h_-^*( k) ], 
    
    h_±( k)    ≡ - μ - D k^2 + v_F( k_y σ_x -  k_x σ_y ) 
       +  (m_0 ±λ +m_1 k^2 )σ_z,

where the coefficients are given by

    -μ   ≡(E_- + E_+)/2 + C, 
    
        D    ≡(B_2/2) ( ⟨φ^↑|τ_z |φ^↑⟩ + ⟨χ^↑|τ_z |χ^↑⟩) - D_2, 
    
        v_F   ≡ -i A_2 ⟨φ^↑|τ_x |χ^↓⟩∈ℝ, 
    
        2 m_0    ≡ E_- - E_+, 
    
        m_1    ≡- (B_2/2) ( ⟨φ^↑|τ_z |φ^↑⟩ - ⟨χ^↑|τ_z |χ^↑⟩).

In Table <ref> we list the values of the coefficients m_0, m_1, D and v_F for several values of the thickness, starting from two quintuple layers, and calculated from Eq. (<ref>), using the parameters of the 3D k· p Hamiltonian (<ref>) for bulk Bi_2Se_3 <cit.>. A diagonalization of H_MTI^2D( k) yields four energy bands with a finite gap at the Γ point. This is shown in Fig. <ref>(a) where, for simplicity, we considered the chemical potential of the system to be tuned to μ =0. 






The eigenstates |φ^σ⟩ and |χ^σ⟩ allow the construction of states |t,σ⟩ and |b,σ⟩ localized near the top and bottom surfaces, respectively. These are indeed given by |t,σ⟩ = (|φ^σ⟩ + |χ^σ⟩)/√(2) and |b,σ⟩ = (|φ^σ⟩ - |χ^σ⟩)/√(2) and have localization lengths ∼ 10  (see Fig. <ref>(b)).
In this basis of states, the Hamiltonian reads

    H̃_MTI^2D( k) =     -μ - D k^2 +  v_F( k_y σ_x -  k_x σ_y ) τ̃_z 
       + λσ_z + m( k) τ̃_x,

where m( k) = m_0 +m_1 k^2 and the Pauli matrices τ̃_x,y,z and σ_x,y,z act, respectively, on the top/bottom (t/s) degree of freedom and the spin degree of freedom. Eq. (<ref>) is the starting point of Sec. <ref>, where we study an MTI/SC heterostructure.


In case of magnetic doping (λ≠ 0), the symmetry class of the Hamiltonian is denoted A (unitary symmetry class) and the topological phase is characterized by an integer Chern number <cit.>. In contrast, at λ = 0, the Hamiltonian has TRS and fits in the symmetry class AII (symplectic symmetry class) where the topological sector is characterized by a ℤ_2 topological invariant.

Determining the value of the topological invariant requires not only the low-energy eigenstates but also information at large momenta. However, determining the difference between the topological invariants characterizing two phases separated by the closing of the gap is possible: in this case, we only need the eigenstates around the Γ point, where the gap closes, which we determine from the Hamiltonian in Eq. (<ref>). These eigenstates and thus the topological invariants do not depend on μ or D since those enter the Hamiltonian with an identity matrix. Hence, let us consider the case μ =0.

The low-energy model can predict the occurrence of topological phase transitions via a sign change of the effective gap m_0 ±λ. 
At λ = 0, the sign change of the gap m_0, when increasing the thickness from 2 to 3 quintuple layers, shows a transition between a trivial insulating and a quantum spin Hall insulating phase (see Table <ref>). At 0 < |λ| < |m_0|, the Chern number vanishes, indicating a trivial phase. In contrast, a sufficiently large magnetic exchange term |λ|>|m_0| yields a quantum anomalous Hall (QAH) phase, independently of the sign and magnitude of m_0, as was also observed in Ref. <cit.>.






Whether or not the materials we consider indeed display a quantum spin Hall phase (without doping) below 6 quintuple layers is still under debate in the literature. For instance, for Bi_2Se_3 it has been argued that the Coulomb interaction leads to a significant hybridization of the surface states, which could open a trivial gap below 6 quintuple layers <cit.>. GW computations performed for Bi_2Se_3 also concluded that the gap below 6 quintuple layers is trivial, but a gap inversion remains for Bi_2Te_3 <cit.>. Our following study will apply to a thin slab in (or near) the QAH regime, which is attainable independently of the sign of m_0 for a large enough magnetic exchange magnitude λ. Experimental studies on such systems have already been performed, and the QAH phase has been characterized, e.g., in Cr- and V-doped (Bi,Sb)_2Te_3 <cit.>.




 §.§ One-dimensional model for an MTI nanoribbon
 

Next, we consider a 2D nanoribbon geometry with a finite width W along the y direction with left and right edges of the slab respectively located at y=0 and y=W. Moreover, we assume λ>|m_0| which would correspond to the QAH regime for W →∞. The goal of this section is to estimate the critical width above which the nanoribbon effectively enters a QAH phase, manifested in one chiral state at each edge of the slab. Below this width, the hybridization of edge states on opposite edges becomes significant. 


The Hamiltonian (<ref>) describes the 2D MTI slab and consists of the two independent blocks h_+( k) and h_-^*( k). For the description of the nanoribbon, we substitute k_y → -i ∂_y and we impose vanishing wave functions at the edges of the slab. In the following, we assume (m_0 - λ) m_1 <0 without loss of generality, since the case (m_0 - λ) m_1 >0 would only exchange the roles of h_-^*(k_x,-i ∂_y) and h_+(k_x,-i ∂_y). The block h_-^*(k_x,-i ∂_y) has an inverted mass gap (m_0 - λ) m_1<0, which results in a pair of edge states arising from the quantum confinement along the y direction. The associated energies at k_x=0 converge to -μ+(m_0-λ)D/m_1 when W or λ is increased as shown respectively in Fig. <ref>(a) or Fig. <ref>(b) (red line). 
In contrast, h_+(k_x,-i ∂_y) has a normal mass gap (m_0 + λ) m_1>0, so for the parameter regime we consider its low-energy states arising from the quantum confinement are topological trivial. They have energies greater than |m_0 + λ| or smaller than -|m_0 + λ|, respectively, as shown in Fig. <ref>. 

For non-vanishing D, the energies of both topological edge states are shifted compared to the energies of the bulk states. As a consequence, for large enough D, i.e., |(m_0-λ)D/m_1|>|m_0 + λ|, the low-energy states at the Γ point are outside the energy gap. As we would like to build a simple ribbon model, we will consider small values of D in the following, such that the low-energy states at the Γ point are in the gap. Hence, in Fig. <ref> we have assumed D=-1 eV. For a description of the effects of larger values of D, we refer to Ref. <cit.>.


At μ= (m_0-λ)D/m_1, where the Fermi energy corresponds to the red line in Fig. <ref>, we can construct a Hamiltonian which describes the edge states around the Fermi level by only considering the block h_-^*(k_x,-i ∂_y), as it describes the relevant low-energy physics. Substituting k_y → -i ∂_y and imposing vanishing wave functions at the edges of the slab, h_-^* becomes,

    h_-^* = h_0-(-i ∂_y) + h_1-(k_x), 
    
       h_0-(-i ∂_y) = D ∂_y^2 -i  v_Fσ_x ∂_y  +  (m_0 - λ -m_1 ∂_y^2 )σ_z,
       h_1-(k_x) = -μ - D k_x^2 +  v_F  k_x σ_y   +m_1 k_x^2 σ_z.

Let us denote the lowest eigenergies at k_x=0 by E_+^y and E_-^y, and the associated eigenstates by χ_y and φ_y. Similarly to the previous section, we project h_1-(k_x) on the low-energy eigenstates of h_0-(-i ∂_y) and thus obtain an effective 1D Hamiltonian describing the slab with finite width W at low energies (see App. <ref>)

    h_-^1D(k_x) =   -μ + E_1 -D̃ k_x^2  -  ṽ_F k_x σ̃_x + (m̃_0 + m̃_1 k_x^2 ) σ̃_z,

where σ̃_x,z are Pauli matrices acting in the {φ_y,χ_y} basis and the parameters are given by

    2 E_1    = E_+^y + E_-^y, 
    
        2m̃_0    ≡-E_+^y + E_-^y 
    
        2D̃   ≡ 2D-m_1(⟨φ_y|σ_z |φ_y⟩ +⟨χ_y|σ_z |χ_y⟩), 
    
        2m̃_1    ≡ m_1 (⟨φ_y|σ_z |φ_y⟩ -⟨χ_y|σ_z |χ_y⟩), 
    ṽ_F   ≡ v_F⟨φ_y|σ_y |χ_y⟩.





The coefficients D̃ and m̃_1 are plotted as functions of the width W in Fig. <ref>. The coefficient ṽ_F varies only very weakly: for the parameters considered in Fig. <ref>, ṽ_F = 1 eV and it has a maximum variation ∼ 10^-3eV over the range of W investigated.

Next, we introduce the states |L⟩ = (|φ_y⟩ + |χ_y⟩)/√(2) and |R⟩ = (|φ_y⟩ - |χ_y⟩)/√(2) which are localized on the left and right edges of the ribbon, respectively. The parameter m̃_0 represents the overlap between both edge states, so it determines the topological phase of the slab. For small overlap, in the QAH phase, one chiral state appears at each edge. For large m̃_0, the overlap of these edge states brings the system into a topologically trivial phase. In this 1D limit, a symmetry-class A Hamiltonian has indeed trivial topological properties <cit.>.



From Fig. <ref>, we see that for W ≲ 500 , m̃_0 cannot be neglected compared to the other terms in the Hamiltonian, so the 1D system becomes topologically trivial. For W ≳ 500 , on the other hand, m̃_0 is negligible, so the slab enters a QAH phase. 







§ MTI/SC HETEROSTRUCTURES
 



 §.§ Surface states of the MTI and proximity effect


We now consider the superconducting proximity effect <cit.> originating from the contact with s-wave superconductors on the top and bottom surfaces of the MTI slab. We assume that they give rise to the pairing potentials Δ_1 and Δ_2, respectively, with two real parameters Δ_1 and Δ_2. This is described by extending the Hamiltonian (<ref>) to the Nambu basis,

    H_BdG^2D( k) 
       =  
         v_F k_y σ_x τ̃_z  + [ - v_F k_x σ_y τ̃_z + λσ_z 
       + 
         m( k)  τ̃_x  - μ - D k^2 ] γ_z 
       - 
        ( Δ_1 + Δ_22  + Δ_1 - Δ_22τ̃_z ) σ_y γ_y.

where the Pauli matrices γ_y,z act in the particle-hole space. Expressed in this form, it is straightforward to check that the Hamiltonian has particle-hole symmetry 𝒫 H_BdG^2D( k)  𝒫^-1 = -H_BdG^2D(- k) and at λ=0 it also has TRS Θ H_BdG^2D( k)  Θ^-1 = H_BdG^2D(- k), where Θ = i σ_y 𝒦, 𝒫 = γ_x 𝒦, and 𝒦 is the complex conjugation operator.



 §.§ Topological properties for the slab geometry
 

For a slab geometry, the topological phase transitions are signaled by the change of the 2D bulk invariant of the model (<ref>) and happen at the gap closing points. For the special case Δ_1=-Δ_2, μ=0 and D =0 the Hamiltonian can be simplified to (see App. <ref>)

    H_BdG^2D   = ∑_ k∑_η=±∑_κ=±ψ_κ, k^η† h_κ^η(  k) ψ_κ, k^η , 
    
        h_κ^η( k)    =   v_F(  k_y σ_x  - κ 
     k_x σ_y ) +  m_κ^η( k) σ_z,

with m_κ^η( k) = m_0,κ^η + m_1,κ^η k^2 and m_0,κ^η = κλ  + ηκ m_0   + |Δ|, and m_1,κ^η = ηκ m_1.
It is then straightforward to calculate the topological invariant <cit.>. The Hamiltonian is a sum of independent massive Dirac Hamiltonians, so the topological invariant is the sum N=∑_η,κ N_κ^η of the winding numbers N_κ^η associated with each of the independent Hamiltonians h_κ^η( k). At fixed m_1, the cases m_0>0 and m_0<0 both result in the same phase diagram for N, even though the decompositions in terms of N_κ^η are different. In Fig. <ref>, we show this phase diagram along with the details of the decomposition in terms of N_κ^η for m_0<0 and m_1>0. 





Building on this, we can extend the phase diagram for Δ_1 ≠ -Δ_2, μ≠ 0 and D ≠ 0 by considering the phase boundaries. As topological phase transitions happen only at the gap closing points, we find that they are described by the equation

    α^2 Δ^4  + Δ^2 [2 α m_0^2- (1+ α^2) (λ^2-μ^2)] 
       + (λ^2+μ^2 - m_0^2)^2-4λ^2 μ^2 =0,

with α≡Δ_2/Δ_1. This result generalizes the study performed in Ref. <cit.> where an equation for the phase boundaries was given for μ=0. 
The D k^2 term, being proportional to γ_z and involving only terms of second order in the momentum, does not influence the phase boundaries.



 §.§ Ribbon geometry for mu=0, Delta1=-Delta2
 

Next we investigate the effect of in-plane confinement on the topological properties of the MTI/SC heterostructure for the case where μ=0 and Δ_1=-Δ_2. We show that for a confined geometry, the topological properties depend on the decomposition of N in terms of N_κ^η, and on the mass values m_0,κ^η. This is in contrast to the translation-invariant geometry, where N alone determines the topological properties at λ≠ 0.

First, we consider a ribbon geometry with width W along the y direction. The calculation performed in Sec. <ref> can be easily adapted since the Hamiltonian (<ref>) is a sum of independent massive Dirac Hamiltonians.
For each term where the mass m_0,κ^η and the parameter m_1,κ^η correspond to the inverted regime, i.e., m_0,κ^η m_1,κ^η<0, a pair of low-energy states appears, which we denote by |φ_y,κ^η⟩ and |χ_y,κ^η⟩ and which have the respective energies m_0,κ^η and - m_0,κ^η. These energies are plotted in Fig. <ref> as a function of the width of the ribbon. Figure <ref>(a) shows the localization of |φ_y,κ^η⟩ and |χ_y,κ^η⟩ along the y direction in an MTI/SC nanoribbon with width W=2500 . One finds that m_0,κ^η decreases when m_0,κ^η increases and that the overlap between |φ_y,κ^η⟩ and |χ_y,κ^η⟩ is proportional to m_0,κ^η (see Fig. <ref>(a) and Fig. <ref>(a)). 

For large regions of the phase diagram the masses m_0,κ^η can differ significantly. At intermediate widths W, this causes the coexistence of chiral low-energy states strongly localized at the edges with states which overlap along the y direction <cit.>. For instance, in Fig. <ref>, for 1000 < W < 5000 , m̃_0,+^- (blue line) and m̃_0,-^- (red line) are small, so we expect that |φ_y,+^-⟩, |χ_y,+^-⟩, |φ_y,-^-⟩, and |χ_y,-^-⟩ are edge states. In contrast, m̃_0,-^+ (black line) is larger so the states |φ_y,-^+⟩, and |χ_y,-^+⟩ overlap along y. 
The states overlaping along the y direction are then described by a 1D bulk Hamiltonian, with parameters m̃_0,κ^η and m̃_1,κ^η = m_1,κ^η (⟨φ_y,κ^η|σ_z |φ_y,κ^η⟩ -⟨χ_y,κ^η|σ_z |χ_y,κ^η⟩ )/2 which depends on W.

Next, we also consider confinement along the x direction such that the length in the x direction satisfies L ≫ W. In this case, we find that the overlapping states |φ_y,κ^η⟩ and |χ_y,κ^η⟩ give rise to new low-energy states |φ_x,κ^η⟩ and |χ_x,κ^η⟩ which are localized at both ends of the ribbon if m̃_0,κ^ηm̃_1,κ^η<0 (see Fig. <ref>(b)). These states are Majorana bound states (MBS), which arise from the confinement of a 1D BdG D Hamiltonian with topologically non-trivial ℤ_2 number.

We are interested only in states with negligible energy, i.e., with energy below a small energy threshold E_ths. If m_0,κ^η m_1,κ^η<0, the Hamiltonian h_κ^η(  k) has a non-zero BdG D Chern number <cit.>, and if m̃_0,κ^η<E_ths, the states |φ_y,κ^η⟩ and |χ_y,κ^η⟩ represent a chiral Majorana edge state (CMES) at each edge of the system. In this description, CMESs with different chiralities can coexist at the same edge if m_0,κ^η m_1,κ^η<0 and m̃_0,κ^η<E_ths for several {η,κ}. This description is valid at vanishing disorder in the system. At finite disorder, CMESs of opposite chirality will hybridize, resulting in 0, 1 or several copropagating CMESs at each edge. It is worth noting that a pair of copropagating CMESs at each edge is topologically equivalent to a QAH chiral edge state <cit.>. If m_0,κ^η m_1,κ^η<0, m̃_0,κ^η>E_ths, and m̃_0,κ^ηm̃_1,κ^η<0, then the states |φ_x,κ^η⟩ and |χ_x,κ^η⟩ form a pair of MBS with negligible energy, localized at the end of the nanoribbon. If disorder were included in our description, we would either find 0 or 1 pairs of MBS in the system. 

As an example, let us consider a Bi_2Se_3 nanoribbon with thickness d=30  (see Table <ref>) and with m_0<0, m_1>0 (the Bi_2Se_3 nanoribbon is a quantum spin hall insulator at λ=0 and Δ_1=Δ_2=0). We consider the regime where |m_0| ≫ |Δ_1|, which corresponds to the experimentally relevant range of superconducting pairings. We arbitrarily set E_ths to be one tenth of the superconducting pairing magnitude |Δ_1|.
In Fig. <ref>, we display the number of CMES and the number of MBS appearing below E_ths, as a function of λ and W.

In the limit W →∞, the number of CMES is in agreement with the phase diagram of Fig. <ref>. Moreover, for the region of Fig. <ref> where N=1, and when W is small enough such that there is no edge states in the y direction, the system simply hosts a pair of Majorana bound states <cit.>. For intermediate range of widths, we observe a richer phase diagram. Indeed, at fixed λ, varying the width up to ∼ 10^4 , we cross different phases. Namely, we observe a region characterized by either one CMES or two CMES with the same chirality and MBS.


In Fig. <ref>, we display the number of CMESs and the number of MBSs in the case m_0>0, m_1>0 (the Bi_2Se_3 nanoribbon is a trivial insulator at λ=0 and Δ_1=Δ_2=0). For a straightforward comparison with Fig. <ref>, we consider m_0 with opposite sign and otherwise identical parameters. For wide enough nanoribbons (W ≳ 2 × 10^4), we observe either 0, 1 or 2 CMESs in agreement with Fig. <ref>. Similarly to the results shown in Fig. <ref>, when W is small enough, we observe 0, 1, and 0 end states, respectively, for the regions of Fig. <ref> where N=0, N=1, and N=2. The region of intermediate width is qualitatively different from what we observe in Fig. <ref>. Only in the region of Fig. <ref> where N=2, we observe a phase characterized by one CMES and MBS. 






















 §.§ Ribbon geometry for general parameters
 

In the previous section, the choice of parameters D=0, Δ_1 = - Δ_2, and μ = 0 made it possible to calculate the low-energy states for the in-plane confined geometry without further approximations. Here, we discuss the more general parameter regime D≠0, μ≠ 0, μ≪ m_0, and Δ_1 ≠ - Δ_2. 

In the following, we consider Δ_1=2 meV and we write Δ_2=αΔ_1. Moreover, we restrict our study to the case -1≤α <1, and μ≪ E_3,4 where E_3,4 are the two highest energies associated to the Hamiltonian (<ref>) for the specific MTI thickness d=30 . For this regime of parameters, we assume that each pair of energy bands around the Γ point (see Eq. (<ref>)) can be described to a good approximation by a massive Dirac Hamiltonian. Therefore, close enough to the Γ point, we express H_BdG^2D( k) in Eq. (<ref>) by H( k)=∑_i=1^4H_i( k) where

    H_i( k) = v_i( k_y σ_x - k_x σ_y) + ( M_i  - B_i k^2 ) σ_z,

where the Pauli matrices act in a transformed basis of states which coincides with the low-energy states of H_BdG^2D( k) at k=0. The parameters v_i, M_i, and B_i are obtained by fitting the resulting energies around the Γ point with the corresponding energy band E_i of H_BdG^2D( k). We further impose that the effective parameters should reduce to m_0,κ^η and m_1,κ^η at μ=0, D=0, Δ_2 = - Δ_1 and that they evolve adiabatically in (λ, D, α, μ) parameter space at constant topological invariant N, which is determined by the phase boundaries evolution according to Eq. (<ref>).

For sufficiently small width W and M_i B_i> 0, the low-energy edge states associated to H_i( k_x, -i ∂_y) hybridize. Similarly to Sec. <ref>, a 1D Hamiltonian H̃_i( k_x) can be derived from H_i( k_x, -i ∂_y),

    H̃_i(k_x) =  -  ṽ_i k_x σ̃_x + (M̃_i - B̃_i k_x^2 ) σ̃_z,

where σ̃_x and σ̃_z are Pauli matrices acting in the {|φ_y⟩,|χ_y⟩} basis where |φ_y⟩ and |χ_y⟩ are the eigenstates of H(k_x=0, k_y → -i ∂_y ) with the two lowest eigenenergies
E_±. The coefficients appearing in the previous equation are given by ṽ_F= v  ⟨φ_y|σ_y |χ_y⟩, 2 m̃_0 = -E_+ + E_-, 2B̃ = B (⟨φ_y|σ_z |φ_y⟩ -⟨χ_y|σ_z |χ_y⟩).




  §.§.§ D ≠ 0
 
First, we consider -1 eV^2 ≤ D ≤ 0 (such that |(m_0-λ)D/m_1|<|m_0 + λ| as we already consider in Sec. <ref>), Δ_1 = - Δ_2 and μ = 0. This case is simplified by the facts that (i) the parameter D appears in the Hamiltonian with a γ_z matrix, in contrast to the m_1 term proportional to τ̃_x γ_z, and (ii) the parameter D multiplies a factor k^2. Therefore the D term does not change the value of the masses M_i and has an important effect only if D k^2 ≳Δ_1. Since here we consider at most D = -1 eV^2, for k ≲ 0.02 ^-1 D has no strong qualitative effect on the low-energy states at the system widths and lengths we consider. 

As an illustration in Fig. <ref>, we show the change of B_i with D for a specific value of λ, and we show how Fig. <ref>(a) is modified for D=-1 eV^2 due to the evolution of the low-energy states ±M_i which appear due to the confinement along y when M_iB_i>0.









  §.§.§ Δ_2 ≠ -Δ_1

At Δ_2 ≠ -Δ_1, the phase boundaries are shifted according to Eq. (<ref>). However, for each phase with a fixed topological invariant, we checked from H( k) that each pair of bulk bands retains the same topological character. Therefore, no strong qualitative changes happen for the phase diagram, as we show in Fig. <ref> for the case α=0.








  §.§.§ μ≠ 0



The case μ≠ 0 is more subtle because the energies around the Γ point can change significantly. For simplicity, we only consider values for λ and μ corresponding to the region with topological phase N=1 in the infinite 2D geometry with max(|Δ_1|,|Δ_2|) ≲μ≪ E_3,4 where E_3,4 are the two highest energies of the Hamiltonian (<ref>). Moreover, we focus on the case W >W_min=10^3, unless otherwise stated. 
From our study in the previous section, we expect this region to be interesting since at μ = 0 it shows a coexistence of two CMES and one MBS (or 1 CMES at very large W). How do these topological states change for μ≠ 0?


Firstly, both high energies E_3,4 remain similar as in the case μ=0, since here μ≪ E_3,4. 
Therefore, above W = 1000 , the low-energy states arising from confinement of the bulk states associated to E_3,4  
are copropagating CMES with negligible energy, as it is the case for μ=0 (blue and red lines in Fig. <ref>(a)).
Secondly, we describe the evolution of the low energy bands E_1,2 at μ≠ 0 by H( k)=∑_i=1^2H_i( k) with H_i( k) given in Eq. (<ref>).  

For the N=1 phase we consider, we know that at W →∞ the topological low-energy states correspond to only one CMES at each edge. This means that the total number of inverted bulk bands cannot change when μ changes. Taking this constraint into account, we observe from our fit that the topological characters of the first and second energy bands are exchanged at small μ with concomitant sign changes of B_1 and B_2. 
Although this has no impact on the phase diagram at W →∞, this is important for the topological properties of the low-energy states for nanoribbons of intermediate width W.






In the example of Fig. <ref>, we observe that M̃_1 takes non-negligible values for W ≲ 15 × 10^3 , with an oscillating behavior as a function of W. This oscillating behavior is also observed as a function of μ and λ when the width W is fixed, as it is shown in Fig. <ref>. Moreover, the sign of B̃_1 and the resulting sign of M̃_1 B̃_1 also oscillate as function of W, μ and λ. For the parameters μ and λ considered in Fig. <ref>, B̃_1 is positive, so the sign changes of M̃_1 B̃_1 are determined by M̃_1. From our theory, this signals topological transitions as W, μ and λ are varied. 

Let us also consider confinement along the x direction, such that the length of the nanoribbon satisfies L ≫ W. Here again, for concreteness, we are interested only in states with energy below a certain threshold E_ths which we arbitrarily set to be one tenth of the superconducting pairing magnitude |Δ_1|. If M̃_1<E_ths, the state associated to the energy M̃_1 is a CMES with negligible energy and with chirality opposite to the copropagating states associated to E_3,4. This description is valid at vanishing disorder in the system. At finite disorder, CMES of opposite chirality hybridize, resulting in only one CMES in the system. If M̃_1>E_ths and M̃_1B̃_1>0 (see, e.g., the red phase in Fig. <ref>), then a pair of MBS with negligible energy, localized at the ends of the nanoribbon, appears. In this case, the low-energy states in the system are both copropagating CMES associated with E_3,4 and a pair of MBS. In Fig. <ref>, we display the number of CMES and the number of MBSs appearing at low energies as a function of μ and λ for W=5200. 

Finally, let us comment the situation in which the nanoribbon is very thin, here meaning W≲500, where the states arising from confinement of the bulk states associated to E_3,4 hybridize in the bulk and do not result in topological low-energy states. In this case, M̃_1>E_ths and M̃_1B̃_1>0 (see Fig. <ref>), which yields a topological phase with a pair of MBS, localized at the ends of the nanoribbon, and without any CMES, reflecting the scenario proposed in Refs. <cit.>. 
Note that this topological phase can arise for λ both greater and smaller than |m_0|, hence not requiring the MTI system (with Δ_1 = Δ_2 = 0) to be in the QAH phase for very large widths. While this phase only appears in a narrow window of λ near |m_0
|, this region is enlarged when μ is increased, similar to the topological phase with 2 CMES and 1 MBS in Fig. <ref>.












§ CONCLUSION

We have studied the topological properties of finite MTI/SC heterostructures using symmetry-constrained low-energy models. We started by developing analytical models for MTI slabs with a finite thickness as well as MTI nanoribbons with finite thickness and width. We investigated the appearance of low-energy states as a function of the magnetic doping, the chemical potential and the system size. Next, we considered such finite geometries subject to superconducting pairing induced by two superconductors at the top and bottom surfaces. For very wide nanoribbons the low-energy states are the chiral edge states as predicted by the 2D bulk topological invariant. For finite width nanoribbons, we constructed and studied low-dimensional models describing the low-energy properties of our system. In a nanoribbon geometry with finite width and length, we observed regions where the low-energy states can host coexisting chiral edge states and Majorana bound states depending on the strength of the magnetic exchange term. Finally, we investigated the effect of a finite chemical potential on the topology of our system. We have studied how the bulk invariant is modified and we have built low-energy models to study the modifications in the low-energy states which appear at the boundaries of the system. When varying the magnetic doping, the chemical potential and the size of the system, we observed topological transitions between two phases which differ by the presence or absence of MBS.


This project is financially supported by the QuantERA grant MAGMA, by the National Research Fund Luxembourg under the grant INTER/QUANTERA21/16447820/MAGMA, by the German Research Foundation under grant 491798118, by MCIN/AEI/10.13039/501100011033 under project PCI2022-132927, and by the European Union NextGenerationEU/PRTR.
L.S. acknowledges support from Grants No. PID2020-117347GB-I00 funded by MCIN/AEI/10.13039/501100011033 and No. PDR2020-12 funded by GOIB.
K.M. acknowledges the financial support by the Bavarian Ministry of Economic Affairs, Regional Development and Energy under Grant No. 07 02/686 58/1/21 1/22 2/23 and by the German Federal Ministry of Education and Research (BMBF) via the Quantum Future project ‘MajoranaChips’ (Grant No. 13N15264) within the funding program Photonic Research Germany.







§ EFFECTIVE RIBBON HAMILTONIAN VIA A PROJECTION ONTO THE LOW-ENERGY STATES
 

Here we determine the lowest eigenergies and the associated eigenstates of the Hamiltonian block,

    h_0-(-i ∂_y) = D ∂_y^2 -i  v_Fσ_x ∂_y  +  (m_0 - λ -m_1 ∂_y^2 )σ_z

The energies E are given by the transcendental equation

    λ_2[ ξ (m_0 - λ) -E - (-D + ξ m_1 ) λ_1^2 ]λ_1 [ ξ (m_0 - λ) -E - (-D + ξ m_1 ) λ_2^2 ] 
          =  tanhλ_2 wtanhλ_1 w ,

with w ≡ W/2,

    λ_α =  √(-F + (-1)^α-1√(R)2(m_1^2+D^2)^2),  α={1,2},

as well as F= v_F^2 - 2 D E + 2 (m_0-λ) m_1 and R=F^2-4(D^2-m_1^2) [ E^2 - (m_0-λ)^2]. The values ξ=± 1 denote, respectively, the solutions for the energy E_+ with eigenstate |χ_y⟩ and the energy E_- with eigenstate |φ_y⟩. We restrict the study to both low-lying energies. Their dependency with respect to the width W is shown in Fig. <ref>.
The eigenvectors are given by

    |χ_y⟩ =  N_1 [ -(-D -m_1 ) η_2 f_+(y);        i  v_F   f_-(y) ],

and 

    |φ_y⟩=  N_2 [ -(-D -m_1 ) η_1 f_-(y);        i  v_F   f_+(y) ],

with 

    f_+(y)   =coshλ_1 ycoshλ_1 w  - coshλ_2 ycoshλ_2 w, 
     f_-(y)   =sinhλ_1 ysinhλ_1 w - sinhλ_2 ysinhλ_2 w, 
    η_2    = λ_1^2 -λ_2^2 λ_1 tanhλ_1 w -λ_2 tanhλ_2 w, 
    η_1    = [ λ_1^2 -λ_2^2  ]λ_1 λ_1 w  -λ_2 λ_2 w ,

and normalization constants N_1 and N_2.

Next we project h_1-(k_x),

    h_1-(k_x) = -μ - D k_x^2 +  v_F  k_x σ_y   +m_1 k_x^2 σ_z,

on the low-energy eigenstates of h_0-(-i ∂_y) and we obtain an effective Hamiltonian for the system initially described by h_-^*, valid around k_x=0. Here we consider only the two low-energy eigenstates of h_0-(-i ∂_y) and we assume that these eigenstates are well separated in energy from the other eigenstates. This assumption is valid for the mass inverted regime (m_0- λ)m_1 <0 we are considering here.
Then we find the Hamiltonian given by the Eq. (<ref>).



§ TRANSFORMING THE HAMILTONIAN INTO A SUM OF DIRAC TERMS IN SPIN SPACE
 

The matrix σ_z τ̃_x commutes with the Hamiltonian at μ=0 and D=0 and Δ_1 =- Δ_2,

    H_BdG^2D( k) 
    =  
            v_F k_y σ_x τ̃_z  + [ - v_F k_x σ_y τ̃_z + λσ_z 
       + 
         m( k)  τ̃_x  ] γ_z - 
        Δ_1 τ̃_z σ_y γ_y.

Therefore it is possible to diagonalize H_BdG^2D( k) and σ_z τ̃_x using a common basis transformation defined via a unitary matrix U which diagonalizes σ_z τ̃_x. For U = ( τ̃_z + τ̃_x )/√(2) one finds Uσ_z τ̃_x U^†= σ_z τ̃_z. This diagonal matrix has only two different eigenvalues ± 1. It is convenient to define P as the projector on the eigenspace with eigenvalue 1 whereas P=1-P projects on the eigenspace with eigenvalue -1. Then we have, 

    H_BdG   = ∑_ kψ_ k^† U^† P U H_BdG( k) U^† P U ψ_ k
       + ∑_ kψ_ k^† U^†P U H_BdG( k) U^†P U ψ_ k.

Next, we have

    U H_BdG( k) U^†   = v_F k_y σ_x τ̃_x - |Δ| τ̃_x σ_y γ_y 
       + [ - v_F k_x σ_y τ̃_x + λσ_z  +  m( k)  τ̃_z ] γ_z.

Choosing P =( σ_z τ̃_z + 1 )/ 2 and using τ̃_z P=σ_z P and τ̃_z P=-σ_z P, we write 

    H_BdG   = ∑_ kψ_ k^† U^† P H_BdG^P( k)  P U ψ_ k
       + ∑_ kψ_ k^† U^†P H_BdG^P( k) P U ψ_ k,

with 

    H_BdG^P,P( k)    =  A_2 k_y σ_x τ̃_x  -A_2 k_x σ_y τ̃_x γ_z 
        + [ λ± m( k)  ] σ_z  γ_z - |Δ| τ̃_x σ_y γ_y,


Now we see that applying a basis transformation which diagonalizes τ̃_x makes the Hamiltonian diagonal in {top,bottom} space. Therefore we perform another unitary transformation using U,

    H_BdG=    ∑_ kψ_ k^P † U H_BdG^P( k) U^†ψ_ k^P 
       + ∑_ kψ_ k^P† U H_BdG^P( k) U^†ψ_ k^P,

with ψ_ k^P =U P U ψ_ k and ψ_ k^P =U P U ψ_ k and we have

    U H_BdG^P,P( k) U^†   =  A_2 k_y σ_x τ̃_z  -A_2 k_x σ_y τ̃_z γ_z 
       + [ λ± m( k)  ] σ_z  γ_z - |Δ|  σ_y τ̃_z γ_y,

Now we almost have a Hamiltonian which is diagonal in {top,bottom} space and in Nambu space. Only the last term in the two previous equations remains off-diagonal. Another unitary transformation, which we denote by U_1 = (σ_x+iγ_y )/√(2), makes the total Hamiltonian diagonal in {top,bottom} space and in Nambu space, 

    H_BdG= H_BdG,1 + H_BdG,2,

with H_BdG,1 = ∑_ kψ_ k^P † U_1^† H_1^P( k) U_1 ψ_ k^P, H_BdG,2 = ∑_ kψ_ k^P† U_1^† H_1^P U_1 ψ_ k^P
and H_1^P( k) = U_1 U H_BdG^P( k) U^† U_1^†, H_1^P=U_1 U H_BdG^P( k) U^† U_1^†. Then, we indeed have

    H_1^P,P( k)    =  A_2 k_y σ_x τ̃_z  +A_2 k_x σ_y τ̃_z γ_z 
       - [ λ±  m( k)  ] σ_z  γ_z + |Δ|  σ_z τ̃_z ,

The Hamiltonian is now diagonal in {top,bottom} space and in Nambu space. We now project it over the eigensubspace of γ_z and of τ̃_z.

First we define P_γ_z = (γ_z + 1 )/ 2 as the projector on the eigenspace with eigenvalue +1 of γ_z and P_γ_z = 1 - P_γ_z projects on the eigenspace with eigenvalue -1 of γ_z. Then we obtain

    H_BdG,1   =∑_ kψ_ k^0,P† H_1^0,P( k) ψ_ k^0,P
       +∑_ kψ_ k^1,P† H_1^1,P( k) ψ_ k^1,P,

with ψ_ k^0,P = P_γ_z U_1 ψ_ k^P, and ψ_ k^1,P = P_γ_z U_1 ψ_ k^P, and

    H_1^(0,1),P( k)    =  A_2 k_y σ_x  τ̃_z +A_2 k_x σ_y  τ̃_z 
       ∓[ λ  +  m( k)  ] σ_z + |Δ|  σ_z  τ̃_z.

Moreover we notice that the unitary transformation σ_y τ̃_x performed in the eigenspace with eigenvalue +1 of γ_z gives

    σ_y τ̃_x H_1^0,P( k) σ_y τ̃_x = H_1^1,P( k)  ,

and because

    ψ_ k^0,P†σ_y τ̃_x H_1^1,P( k) ψ_ k^1,P = 0,

we obtain

    H_BdG,1= ∑_ kψ̃_ k^P† H_1^1,P( k) ψ̃_ k^P	 .

with ψ̃_ k^P† = ( ψ_ k^0,P†σ_y τ̃_x+ψ_ k^1,P†).
Similarly, defining ψ_ k^0,P = P_γ_z U_1 ψ_ k^P, ψ_ k^1,P = P_γ_z U_1 ψ_ k^P and 

    H_1^1,P( k) = A_2 k_y σ_x  τ̃_z -A_2 k_x σ_y  τ̃_z + [ λ  -  m( k)  ] σ_z + |Δ|  σ_z  τ̃_z ,

we obtain

    H_BdG,2=    ∑_ kψ̃_ k^P† H_1^1,P( k) ψ̃_ k^P,

with ψ̃_ k^P† = ( ψ_ k^0,P†σ_y τ̃_x+ψ_ k^1,P†). Hence, H_BdG now reads

    H_BdG   = ∑_ kψ̃_ k^P† H_1^1,P( k) ψ̃_ k^P
       + ∑_ kψ̃_ k^P† H_1^1,P( k) ψ̃_ k^P,

Now we project this Hamiltonian over the {top,bottom} space. We define P_τ̃_z = ( τ̃_z + 1 )/ 2 as the projector on the eigenspace with eigenvalue +1 of τ̃_z and P_τ̃_z = 1 - P_τ̃_z. Then we have
 
    ∑_ kψ̃_ k^P† H_1^1,P( k) ψ̃_ k^P=
       ∑_ kψ_ k^0,1,P † H_1^0,1,P( k) ψ_ k^0,1,P + ∑_ kψ_ k^1,1,P † H_1^1,1,P( k) ψ_ k^1,1,P,

with ψ_ k^0,1,P =  P_τ̃_zψ̃_ k^P
and ψ_ k^1,1,P =  P_τ̃_zψ̃_ k^P and

    H_1^0,1,P( k) = A_2 k_y σ_x  -A_2 k_x σ_y  + [ λ  +  m( k)  ] σ_z + |Δ|  σ_z , 
       H_1^1,1,P( k) = -A_2 k_y σ_x  +A_2 k_x σ_y  + [ λ  +  m( k)  ] σ_z - |Δ|  σ_z ,

Let us define

    h_±^P ( ks)=A_2 k_y σ_x  -A_2 k_x σ_y  + [ λ  +  m( k)   ± |Δ| ] σ_z

Then we have H_1^0,1,P( k) = h_+^P ( k) and σ_y  H_1^1,1,P( k) σ_y = - h_-^P *(- k) so we conclude that 
 
    ∑_ kψ_ k^P † U_1^† H_1^P( k) U_1 ψ_ k^P=   ∑_ kψ̃_ k^P †H̃_1^P( k) ψ̃_ k^P,

with ψ̃_ k^P = (ψ_ k^0,1,P, σ_y ψ_ k^1,1,P)^T and

    H̃_1^P( k) = [      h_+^P( k)              0;              0 - h_-^P *(- k) ]

The same transformations can be applied to the other part of the Hamiltonian and we find
 
    ψ̃_ k^P† H_1^1,P( k) ψ̃_ k^P=∑_ kψ̃_ k^P†H̃_1^P( k) ψ̃_ k^P,

with ψ̃_ k^P = (ψ_ k^0,1,P, σ_y ψ_ k^1,1,P)^T and

    H̃_1^P( k) = [      h_+^P( k)              0;              0 - h_-^P *(- k) ],

with ψ_ k^0,1,P =  P_τ̃_z (σ_y τ̃_x P_γ_z U_1 ψ_ k^P +  P_γ_z U_1 ψ_ k^P ), ψ_ k^1,1,P =  P_τ̃_z (σ_y τ̃_x P_γ_z U_1 ψ_ k^P +  P_γ_z U_1 ψ_ k^P ) and

    h_±^P( k)=A_2 k_y σ_x  -A_2 k_x σ_y  + [ λ  -  m( k)   ± |Δ| ] σ_z,

To sum up we have

    H_BdG= ∑_ kψ̃_ k^P †H̃_1^P( k) ψ̃_ k^P +  ∑_ kψ̃_ k^P†H̃_1^P( k) ψ̃_ k^P,

with 

    H̃_1^P( k) = [      h_+^P( k)              0;              0 - h_-^P *(- k) ], 
       H̃_1^P( k) = [      h_+^P( k)              0;              0 - h_-^P *(- k) ],

and

    h_±^P( k)=A_2 k_y σ_x  -A_2 k_x σ_y  + [ λ  +  m( k)   ± |Δ| ] σ_z, 
       h_±^P( k)=A_2 k_y σ_x  -A_2 k_x σ_y  + [ λ  -  m( k)   ± |Δ| ] σ_z.

We denote the energies of h_+^P( k), - h_-^P *(- k), h_+^P( k) and - h_-^P *(- k) by E_1, E_2, E_3, and E_4, respectively. We have

    E_1^±   = ±√(A_2^2 k_⊥^2 + ( λ  +  m( k) + |Δ|  )^2), 
    
        E_2^±   = ±√(A_2^2 k_⊥^2 + ( λ  +  m( k) - |Δ|  )^2), 
    
        E_3^±   = ±√(A_2^2 k_⊥^2 + ( λ  -  m( k) + |Δ|  )^2), 
    
        E_4^±   = ±√(A_2^2 k_⊥^2 + ( λ  -  m( k) - |Δ|  )^2).

Note that when the system has time-reversal symmetry at λ = 0, we have E_1^± = E_4^± and E_2^± = E_3^±. Generally speaking, in presence of time reversal symmetry, Kramers' theorem tells us that every (spin-1/2) Bloch state is degenerate with its time-reversal conjugate, i.e., a state |ψ_a ( k)⟩ with energy E_a( k) and a state |ψ_b(- k)⟩ = Θ|ψ_a( k)⟩, where Θ is the time-reversal operator, with energy E_b(- k) have the same energies, E_a( k) = E_b(- k). Additionally, we have E_i^±( k) = E_i^±(- k) (i∈{1,2,3,4}) due to the additional presence of inversion symmetry, so that E_a( k) = E_b( k), meaning that the energy bands which come in Kramers pairs E_a( k) = E_b(- k) are not only degenerate at the time-reversal invariant points k=- k but at each k point. From E_1^± = E_4^± and E_2^± = E_3^±, we find that E_1^- and E_4^- form a Kramers pair of energy bands and the same is true for E_1^+ and E_4^+, E_2^- and E_3^- and E_2^+ and E_3^+.






