Self-consistent study of Abelian and non-Abelian order in a two-dimensional topological superconductor

# Self-consistent study of Abelian and non-Abelian order in a two-dimensional topological superconductor

## Abstract

We perform self-consistent studies of two-dimensional (2D) -wave topological superconductivity (TSC) with Rashba spin-orbit coupling and Zeeman field by solving the Bogoliubov-de Gennes equations. In particular, we examine the effects of a nonmagnetic impurity in detail and show that the nature of the spin-polarised midgap bound state varies significantly depending on the material parameters. Most notably, a nonmagnetic impurity in a 2D -wave topological superconductor can act like a magnetic impurity in a conventional -wave superconductor, leading to phase transitions of the ground state as the impurity potential is varied. Furthermore, by solving for the spin-dependent Hartree potential self-consistently along with the superconducting order parameter, we demonstrate that topological charge density waves can coexist with TSC at half filling just as in a conventional -wave superconductor.

###### pacs:
74.20.-z, 74.20.Rp, 74.81.-g, 75.70.Tj

## I Introduction

The Rashba spin-orbit (SO) interaction(1) that occurs in materials without inversion symmetry is a key ingredient for topological insulating or superconducting states and for creating Majorana fermions as elementary excitation in a wide variety of systems.(2) Two-dimensional (2D) -wave topological superconductivity (TSC) with Rashba SO coupling and Zeeman field has been proposed to be realised in an ultracold Fermi gas, where -wave superfluidity and SO coupling can be generated, respectively, via the -wave Feshbach resonance and spatially varying laser fields.(3); (4) 2D -wave TSC has also been proposed to be achievable in a solid device made of more conventional materials,(4); (5) such as a semiconductor with SO coupling sandwiched between a conventional -wave superconductor and a ferromagnetic insulator,(6) or a semiconductor with Rashba and Dresselhaus coupling in proximity to an -wave superconductor only.(7) Furthermore, recently developed, one-atom-layer Tl-Pb compounds on Si(111) with giant Rashba effects(8) and ionic-liquid based double-layer transistors(9) are promising candidate systems for realising 2D -wave TSC.(10) In the latter system, for example, the layered structure can consist of an -wave superconductor and a ferromagnetic insulator, while the electric field controls the charge-carrier density as well as the Rashba SO interaction.(10)

In a 2D -wave topological superconductor, vortices host Majorana fermions as zero-energy bound states(6); (11) and hence obey non-Abelian exchange statistics. Therefore, creating 2D -wave TSC in a condensed matter system, where one has control over key parameters such as the filling factor and the strength of Rashba SO coupling and/or Zeeman field, potentially leads to realisation of fault-tolerant topological quantum computation.(5); (12) For setups where superconductivity is induced by the proximity effect, in general there have been several arguments(13); (14); (5) that coupling to the superconductor can renormalise the original parameters, e.g., on the surface of a three-dimensional topological insulator and can even be detrimental to the other ingredients necessary for TSC; although a self-consistent study of the proximity effect(15) has contradicted some results in this regard in Ref. (13). It is ideal to have an intrinsic pairing interaction in the system that drives -wave superconductivity or superfluidity, as in -wave superfluids of fermionic atoms(3), one-atom-layer superconducting compounds,(8) or electric-field double-layer transistors.(9); (10)

Despite the rapid advancement in device fabrication techniques and the 2D -wave TSC model(3); (4); (6); (5); (7) being one of the most promising models for platforms for topological quantum computation, self-consistent studies of this model are severely lacking. Few studies made so far, in which the superconducting order parameter is solved self-consistently, are a momentum-space study of average impurity effects(16) and studies of a single vortex(17); (18) and the vortex lattice(19) in terms of the tight-binding model;(3); (4) and a study of the effects of a nonmagnetic impurity in an -wave superfluid(20) using the continuum model.(6) To the best of our knowledge, there has been no work where the Hartree potential is solved self-consistently as well.

The purpose of the present work is to perform self-consistent studies of 2D -wave TSC by solving the Bogoliubov-de Gennes (BdG) equations(21) on the tight-binding model of Sato, Takahashi, and Fujimoto(3); (4) directly and numerically. As a tight-binding model, it is versatile in that band structure and the filling factor can easily be adjusted to model real systems in terms of hopping amplitudes and the chemical potential. Moreover, unlike the continuum model,(6) in addition to non-Abelian phase this model can host Abelian phase where topological order can be realised not only in superconducting states, but also in density wave states.(4) In this work, we examine the effects of a nonmagnetic impurity and illustrate how the properties of the midgap bound state vary with the material parameters – in stark contrast to the claim made by Hu et al.(20) on “universal” midgap bound states. In particular, we find that the midgap excitation bound to a nonmagnetic impurity is spin-polarised(22) and its spin can flip as the impurity potential is varied. Most notably, we find that a nonmagnetic impurity in a 2D -wave topological superconductor can act exactly like a magnetic impurity (classical spin) in a conventional -wave superconductor, resulting in phase transitions of the ground state. We also demonstrate that the spin-dependent Hartree potential effectively reduces the Zeeman field. Furthermore, by solving for the Hartree potential as well as the superconducting order parameter self-consistently, we show the existence of topological charge density waves (TCDW) as the ground state that is degenerate with TSC at half filling, just as in a conventional -wave superconductor.

Sau and Demler(23) have studied the effects of a nonmagnetic impurity in a semiconductor nanowire with Rashba SO coupling and Zeeman field in proximity to an -wave superconductor,(24); (25) where one-dimensional (1D) TSC as in Kitaev’s model can be realised.(26); (5) By searching for poles in the Green function, they found no midgap excitation bound to a nonmagnetic impurity when either one of the Rashba SO interaction and the Zeeman field – the two ingredients required for TSC – was missing. In TSC states, when the Zeeman field is larger than the assumed uniform order parameter, they have found midgap bound states whose energy can cross zero as the strength of the impurity potential is increased. Our self-consistent results for zero-energy crossing of the midgap excitation in a 2D -wave topological superconductor are similar to their findings; except that in the 1D system, the impurity cuts the wire in half and creates Majorana edge modes in the limit of infinitely strong potential.

We employ the Chebyshev polynomial expansion method(27); (28) for solving the BdG equations self-consistently for the mean fields as well as calculating the local density of states (LDOS) after self-consistency has been achieved. This method allows one to obtain self-consistent mean fields without diagonalization of the BdG Hamiltonian matrix and it can also gain significant speed-up from parallel computation. It is thus much more efficient than the traditional way of solving the BdG equations by direct diagonalization, especially when the BdG matrix is spin-dependent and required to be solved for relatively large system size. We also circumvent the high numerical demand of diagonalizing the BdG matrix by utilising the efficient algorithm of Sakurai and Sugiura (SS) (29); (30) to obtain the quasiparticle spectrum within a selected energy window. This numerical technique also benefits greatly from parallelism and if desired, the entire spectrum can be obtained readily by dividing the energy range into smaller subranges and applying the SS method to each subrange.

The paper is organised as follows. The model is described in Sec. II, results are presented and discussed in Sec. III, and the work is summarised in Sec. IV.

## Ii Model

We solve the BdG equations on the tight-binding model for 2D -wave TSC(3); (4) with the Hamiltonian,

 H = ∑⟨ij⟩σtijc†iσcjσ+∑iσ(ϵi−μ+hσ+V(H)i¯σ)c†iσciσ (1) + α2[∑i(c†i−^x↓ci↑−c†i+^x↓ci↑) + i(c†i−^y↓ci↑−c†i+^y↓ci↑)+H.c.] + ∑i(Δic†i↑c†i↓+H.c.),

with the Zeeman energy and for and , respectively, and . We assume a uniform pairing interaction that results in the -wave order parameter and the Hartree potential :

 Δi = U⟨ci↓ci↑⟩, (2) V(H)iσ = U⟨c†iσciσ⟩, (3)

where the electron creation and annihilation operators at site with spin are denoted as and , respectively, and is the Hartree potential created by the electrons with spin (and felt by those with the opposite spin ) at site . In the Hamiltonian (1), we consider hopping among nearest-neighbour lattice sites only with the probability amplitude , is the chemical potential, is the single-particle potential due to a nonmagnetic impurity at site , is the Rashba spin-orbit coupling strength, and H.c. stands for the Hermitian conjugate. We set the lattice constant to be unity, and and are the unit vectors in the - and -directions.

Defining the average Hartree potential for each spin component ,

 ¯V(H)σ=1N∑iV(H)iσ, (4)

where is the total number of lattice sites, the diagonal terms of the Hamiltonian in Eq. (1) can be written as

 ∑iσ(ϵi−~μ+~hσ+V(H)i¯σ−¯V(H)¯σ)c†iσciσ (5)

with () for (). We have defined

 ~μ = μ−¯V(H)↑+¯V(H)↓2, (6) ~h = h+¯V(H)↑−¯V(H)↓2. (7)

For , typically there are more spin-up electrons than spin-down electrons and the Hartree potential effectively reduces the Zeeman field. Intuitively, the average energy gain by the pairing interaction with spin-up electrons makes the electron to have its spin down less costly in terms of the Zeeman energy. When the system has translational symmetry, , and the Hamiltonian (1) can be expressed in momentum space as

 H=12∑kΨ†kH(k)Ψk, (8)

where and

 H(k)=(ϵ(k)−~hσz+αL(k)⋅σiΔ(k)σy−iΔ(k)∗σy−ϵ(k)+~hσz+αL(k)⋅σ∗), (9)

with and . and are the creation and annihilation operators of the electron with momentum and spin , and and are the Pauli matrices. above reduces to the momentum-space Hamiltonian given in Ref. (4) when the Hartree potential is neglected, with and . Thus, various topological phases as classified in Ref. (4) according to the first Chern number or the Thouless-Kohmoto-Nightingale-Nijs (TKNN) number(31) can be achieved by replacing the chemical potential and Zeeman field by and , respectively, when the Hartree potential is taken into account. Topological phase transitions between topologically distinct phases occur when the energy gap of the bulk quasiparticle spectrum closes.(4) Assuming an isotropic -wave order parameter, , diagonalization of in Eq. (9) yields

 E±(k)=√ϵ(k)2+α2|L(k)|2+~h2+|Δ0|2±2ξ(k), (10)

where , and the minimum value of is the bulk spectral gap . For example, as is varied for a given set of , , and , the system transitions from one topological phase (trivial, Abelian, or non-Abelian) to another every time vanishes. The topological invariant that classifies each phase is the TKNN number(31) and can be calculated by(32); (33); (34)

 ν= 18π2∫dkdω[Tr(G∂kxG−1G∂kyG−1G∂ωG−1) (11) −Tr(G∂kyG−1G∂kxG−1G∂ωG−1)],

where and .(35) The system is in trivial phase when , and in Abelian and non-Abelian phase, respectively, when is even and odd ( and in this model).

In this 2D -wave TSC model, there are two underlying chiralities, ,(36); (37) as can be seen by expressing in the “chirality basis” that diagonalizes the normal-state Hamiltonian :(4); (19)

 Unknown environment '% (12)

where and

 ^Δ=1Δϵ(k)(−α|L(k)|η+Δ(k)~hΔ(k)−~hΔ(k)−α|L(k)|η−Δ(k)). (13)

The intraband pairing in the band () has the chirality of (), while the interband pairing is purely -wave. This is analogous to the two chiralities present in the nontrivial (non-Abelian) phase of the continuum model.(38); (7); (39) The two chiralities associated with the two Fermi surfaces are always mixed in Abelian phase. In contrast, the chirality associated with the single Fermi surface can be dominant over the other in non-Abelian phase for relatively weak SO coupling, and which chirality is more manifest is determined by the sign of as well as the sign of .(4); (40); (19) When SO coupling is strong, the two chiralities can be mixed strongly also in non-Abelian phase.(41); (19)

To obtain the mean fields, we perform self-consistent iterations up to the -th iteration step, where, e.g., the order parameter as a complex vector of length satisfies

 ∥→Δ(l)−→Δ(l−1)∥∥→Δ(l−1)∥<10−6 (14)

and similarly for each spin component of the Hartree potential that can be regarded as a real vector of length . All the calculation presented below has been performed for zero temperature.

## Iii Results

### iii.1 Self-consistency

We first illustrate how the order parameter in a uniform system varies with the strength of the pairing interaction . In our self-consistent calculation, the order parameter is defined in terms of Eq. (2) such that is positive. Figure 1 presents the order parameter as a function of both in units of in a uniform 5050-site lattice with the periodic boundary condition (PBC), with (green squares) and without (red circles) solving for the Hartree potential. Without (with) the Hartree potential, (), , and (), , for the systems shown in Fig. 1(a) and (b), respectively. When the Hartree potential is neglected, for () the non-Abelian (Abelian) TSC state with the TKNN number () is realised according to the conditions,(4)

 (4t−|μ|)2+Δ20

and the system is in the trivial phase with when (). In both systems presented in Fig. 1 the transition to the trivial phase will occur for . This limits the range of the pairing interaction to and , respectively, for the systems shown in Fig. 1(a) and (b) to stay in the nontrivial phase. As can be seen above, the phase boundaries do not depend on the signs of and . The self-consistently obtained is not affected by the sign of or either.

With the inclusion of the Hartree potential, and in Eqs. (15) and (16) are replaced by and . For , , and , the system now transitions from the non-Abelian to trivial phase at around , where we have found and . The system with , and has a very narrow region of Abelian TSC for , where at we have obtained and . Thus, in general, due to the effective reduction of the Zeeman field by the Hartree potential, one needs not assume too strong a coupling constant for the pairing and/or Rashba SO interaction when solving for the Hartree potential self-consistently. What this implies in terms of real materials is that the Pauli depairing effect of the Zeeman field is overcome partially by the Hartree potential.

### iii.2 Effects of a single nonmagnetic impurity

In a conventional -wave superconductor, a nonmagnetic impurity can locally suppress the order parameter and cause Friedel-like oscillations around it;(42); (43) however, a single nonmagnetic impurity does not bind a quasiparticle. On the other hand, in a spin-triplet chiral -wave superconductor, a nonmagnetic impurity can create midgap quasiparticle excitation and spontaneous supercurrent around it.(44) Furthermore, in the domain where the order is suppressed around a nonmagnetic impurity, the order is induced, and vice versa.(44); (45) It is an intriguing question as to how nonmagnetic impurities affect 2D -wave TSC due to the presence of the two underlying chiralities ( in the continuum model).(20); (16); (22); (39)

Nagai, Ota, and Machida have examined the average effects of nonmagnetic impurities in the 2D -wave TSC model(3); (4) by solving for the impurity self-energy and the uniform order parameter self-consistently.(16) They have shown that while the system in trivial phase () is robust against nonmagnetic impurities, the Anderson theorem(46) can break down in a nontrivial phase () and that the superconducting transition temperature becoming sensitive to impurity concentration is accompanied by the appearance of midgap states. They have further studied quasiparticle bound states around a single nonmagnetic impurity in non-Abelian phase in terms of the BdG equations (without self-consistency)(22) and have found that the bound-state energy decreases with respect to the bulk spectral gap as the Zeeman field is increased. This is consistent with their finding of being more sensitive to nonmagnetic impurities for stronger Zeeman field in Ref. (16) and with the picture that the intraband chiral -wave pairing (with chirality or in non-Abelian phase) dominates over the interband -wave pairing for large .(7); (22)

Hu et al.(20) have obtained “universal” midgap bound states around a nonmagnetic point (delta-function) impurity in a 2D -wave superfluid with Rashba SO coupling and Zeeman field in a harmonic trap by solving the radial BdG equations self-consistently. They have found for one set of parameters in non-Abelian phase that a strong point impurity (nonmagnetic or magnetic) results in a midgap state with the “universal” bound-state energy , where and are the bulk order parameter and the Fermi energy, respectively. As pointed out by Shitade and Nagai,(39) this bound-state energy is not universal as it was derived by substituting the -wave order parameter in the formula for Majorana edge modes in a chiral -wave superfluid in confinement.(47) The argument made by the authors of Ref. (20) is that the impurity site where the order parameter is suppressed to zero acts as a singular point of vacuum and its “interface” with the nontrivial region hosts a Majorana edge state. We note, however, that the self-consistently solved order parameter increases continuously from zero away from the impurity site and there is no clear boundary with which trivial and nontrivial regions can be defined locally.

In this subsection, we demonstrate that the effects of a nonmagnetic impurity on 2D -wave TSC are far from being “universal” and depend significantly on the material parameters (, , and ). For relatively weak SO coupling in non-Abelian phase, where one of the two chiralities can be dominant over the other, we find in general that the smaller the , or the larger the , the more -wave-like the system reacts to a nonmagnetic impurity. Reversing the sign of while keeping all the other parameters the same does not change the order parameter nor the excitation spectrum: only the spin components of quasiparticle excitation are interchanged. In case where one chirality is dominant over the other in non-Abelian phase, or switches the major chirality.(40); (4); (19) Under , changing the sign of the impurity potential at the same time leaves the order parameter and excitation spectrum virtually unchanged. In the results presented below, and the Hartree potential was neglected.

We first illustrate the effects of the SO coupling strength . We place a single nonmagnetic impurity with potential at the centre of the lattice with PBC. The order parameter acquires a small imaginary part in the vicinity of a nonmagnetic impurity, while it remains real at the impurity site. In Fig. 2 we show (a) the quasiparticle bound-state energy in units of the bulk spectral gap and (b) the ratio of the real part of the order parameter at the impurity site to the bulk order parameter as a function of for and , for , , and . The system size is 5151 lattice sites for and , and 6464 for . We find that typically the disturbance of the order parameter by the impurity is contained well within 5151 lattice sites. The coupling constant for the pairing interaction has been chosen to be , and for , and , respectively, such that (). The spin of the quasiparticle excitation has been determined by identifying the peak at the eigenvalue in the spin-resolved LDOS at the impurity site. Consistently with the non-self-consistent results of Ref. (22), we find that the midgap excitation is always spin-polarised at the impurity site.

In these non-Abelian systems the chirality is dominant over , as can be apparent in the core of a vortex pinned by a nonmagnetic impurity.(19) One can see in Fig. 2(a) that the nonmagnetic impurity binds spin-down quasiparticle excitation (i.e., composed of spin-down particle and spin-up hole), which becomes more strongly bound to the impurity as decreases when . This can be understood by the fact that for a given , the smaller the , the more manifest the major chirality becomes and the more chiral -wave-like the system behaves.(39); (41); (19) For the bound-state energy is similar and is suppressed at the impurity site to a similar value for the three values of . In the range where the bound-state energy is substantially different for , , and , is peaked at the impurity site, more sharply for smaller , as can be seen in Fig. 2(b). This is illustrated further in Fig. 3 for , where (a) and (b) the average number of spin-up and spin-down electrons at the impurity site are plotted as a function of . The average number of electrons decreases monotonically as increases, while reaches its maximum value for in this system.

The enhancement of the order parameter at the impurity site can be understood roughly by the fact that the chemical potential is shifted to locally and this effectively makes the site transition into the trivial phase. Although the phase boundaries that separate trivial and nontrivial phases are defined in terms of bulk , , and , and do not directly apply to a single lattice site, we have checked that with no impurity and the other parameters fixed, e.g., for , the BdG equations converge to a trivial state with uniform order parameter similar in value to for shown in Fig. 2(b), for and .

In Fig. 4 the order parameter is plotted as a function of spatial coordinates and for with at the centre of a 5151-site system for and , for (upper graphs) and (lower graphs). The left (right) panel shows the real (imaginary) part of the order parameter. This figure illustrates the chiral -wave nature of the system being enhanced for smaller . As can be seen in Fig. 4, the order parameter is suppressed to nearly zero at the impurity site for both and . For , however, the real part of the order parameter is enhanced at the nearest- and next-nearest-neighbour sites to the impurity (). This is reminiscent of the order parameter of the opposite chirality induced around a nonmagnetic impurity in a chiral -wave superconductor.(44); (45) Moreover, the small, oscillatory imaginary part in the order parameter indicates the presence of induced supercurrent around the impurity, and the oscillation amplitudes are increased slightly for compared to those for .

We have found curious effects of a nonmagnetic impurity by increasing the Zeeman field : the spin of the impurity bound state flips as the impurity potential is varied. Two different examples are presented below. Shown in Fig. 5 are (a) the quasiparticle bound-state energy and (b) the order parameter at the impurity site as a function of for and on a 5151 lattice, for , , and . , and for , and , respectively, that yield (). The energy of the spin-down bound state varies in a way similar to the results shown in Fig. 2(a) as is decreased from the positive side, reaching a maximum at some negative value of . The bound state, however, disappears into the gap edge rather abruptly as is decreased further and then reappears with spin up; although for it remains very close to the gap edge. In this region the order parameter is peaked at the impurity site, once again due to the local shift of the chemical potential to that effectively transforms the site into another phase – either trivial or Abelian depending on the value of .

Figure 6 demonstrates a striking example where a nonmagnetic impurity in a 2D -wave topological superconductor acts like a magnetic impurity (classical spin) in a conventional -wave superconductor. The system presented in this figure is a 5151 lattice with , , , and (). A magnetic impurity in a conventional -wave superconductor locally breaks time-reversal symmetry and creates spin-polarised quasiparticle excitation (the spin direction can be defined with respect to that of the impurity) in the energy gap(48); (49); (50) and as the potential strength of the impurity is increased, the excitation energy crosses zero, after which the ground state contains a spin-polarised quasiparticle (spin- up or down, depending on the definition) bound to the impurity and the midgap excitation has the opposite spin as it would now remove this extra spin from the ground state.(51) When this phase transition of the ground state happens, the order parameter becomes negative at the impurity site(52); (53) – caused by the sign change in the contribution from the impurity bound state(53); (54) – and vanishes in the limit of infinite potential strength.(52)

It can be seen in Fig. 6(a) that as is increased in magnitude (either from the positive or negative side), the energy of the spin-down quasiparticle excitation decreases and crosses zero, at which point the quasiparticle excitation becomes spin up and changes sign from positive to negative, as can be seen in Fig. 6(b). As is increased further, approaches zero. Since the two limits () should result in the same uniform state, the impurity bound-state energy must cross zero an even number of times:(23) this is indeed the case for the quasiparticle excitation shown in Fig. 6(a). The zero-energy crossing occurs at and in this system. We have found the excitation (spin up at the impurity site) with energy for , which is not a Majorana fermion. The spin-resolved LDOS at the bound-state energy is plotted as a function of and coordinates for the entire lattice in Fig. 7 for (upper panel) and (lower panel); where the spin of the quasiparticle excitation is down and up, respectively, at the impurity site. As can be seen in Fig. 7, the spin-up (spin-down) LDOS is zero at the impurity site for (). The diagonal extension of the LDOS reflects in the Fermi wave vector on large portions of the squarish Fermi surface.(19) The LDOS was calculated by using the eigenfunctions obtained by the SS method.(29); (30)

We have also found the same phase transition of the ground state in the presence of a nonmagnetic impurity in the system with , , , and (), where the variation of the bound-state energy as a function of is very similar to that shown in Fig. 6(a) for , but the zero crossing on the negative side occurs at , in contrast to for the system shown in Fig. 6 (). It is interesting to note that for , , and , the spin of the impurity bound state flips as the impurity potential is varied with and without crossing zero energy, respectively, for () and () [Fig. 5(a)].

We summarise the parameter values used for Figs. 2, 5 and 6 in the phase diagram as a function of and shown in Fig. 8, where different topological phases with , , and are separated by the phase boundaries, and , with .(4) The system is in the trivial phase with when and . Although there is no direct correspondence of the phase diagram to the single impurity site, for the systems presented in Fig. 2 and for in Fig. 5, the range of where the order parameter is peaked at the impurity site loosely corresponds to ranges of the local chemical potential where or . In contrast, for in Fig. 5 () and for the system shown in Fig. 6 (), the TKNN number associated with the local chemical potential varies among , and and reaches only for . Interestingly, the spin reversal of the midgap excitation occurs for () for in Fig. 5 and for and ( and ) for the system shown in Fig. 6, which correspond to and , respectively. Namely, the “local” TKNN number at the impurity site has the opposite sign to that of the bulk TKNN number, which is analogous to locally changing the sign of (19) and thus flipping the spin of the quasiparticle excitation. This does not explain, however, the rather striking difference in the behaviour of the impurity bound state between Figs. 5 and 6, which warrants further investigation.

Sau and Demler have found similar zero-energy crossing of midgap excitation bound to a nonmagnetic impurity in a 1D semiconductor with Rashba SO coupling and Zeeman field in proximity to an -wave superconductor.(23) In this system, the bound-state energy also crosses zero in the limit of infinitely strong impurity potential, as it corresponds to the system effectively cut in half by the impurity and the zero-energy states are the Majorana fermions bound at the two edges. Thus, in the case of a 1D topological superconductor, there is an odd number of zero-energy crossings for finite impurity potential.

### iii.3 Coexistence of TCDW and TSC

The 2D s-wave TSC model can also host topological charge density waves (TCDW).(4) With the order parameter for charge density waves (CDW) with wavevector , the mean-field Hamiltonian can be written as(4)

 HCDW=12∑kΨ†CkHCDW(k)ΨCk, (17)

where and

 (ϵ(k)−~hσz+αL(k)⋅σiΔCσy−iΔCσyϵ(k+Q)+~hσz−αL(k+Q)⋅σ∗). (18)

We note that the term in given in Ref. (4) () has a wrong sign: it should be minus as above. If satisfies and , is equivalent to in Eq. (9) for TSC. In our tight-binding model, these two conditions are satisfied for and . For , the only nontrivial phase possible is Abelian with . For the same , , and the uniform order parameter , the TCDW and TSC states then have the same quasiparticle spectrum, with two zero-energy bound states per surface. In the TCDW state, however, the two zero-energy surface states are equivalent to each other due to folding of the Brillouin zone by the CDW order and hence there is only one zero mode per surface.(4) Furthermore, while the zero modes are each a Majorana fermion in the TSC state, quasiparticle excitation in the TCDW state is not and possesses (1) charge.(4)

We have solved the BdG equations with the Hamiltonian in Eq. (1) self-consistently for each spin component of the Hartree potential with PBC, no impurity, and the superconducting (SC) order parameter set to zero, and have found uniform TCDW states with and . For example, , , and on a 6464 lattice yield the TCDW state with , the effective Zeeman field , and . In real space, the CDW order parameter is the deviation of the Hartree potential from its average value, which alternates between and from one site to its nearest-neighbour site, in each spin component. The average electron density (number of electrons per site) is exactly 1.0. This confirms the fact that in the presence of Rashba SO coupling and Zeeman field corresponds to half filling, as for the tight-binding system for conventional -wave superconductivity (with nearest-neighbour hopping only).

We show in Fig. 9 the (a) spin-up and (b) spin-down components of the LDOS as a function of excitation energy at lattice sites = (1,1), (2,1), (1,2), and (2,2) in the above TCDW state for , , and on a 6464 lattice. The LDOS was calculated using the Lorentz kernel(55); (27) with the corresponding Lorentzian smoothing width of . The CDW order is such that the electron density (0.465904) and (0.196576) at sites (1,1) and (2,2) ((2,1) and (1,2)). These density modulations are reflected in the LDOS, showing more hole and electron excitations available for (1,1) and (2,1), respectively, for both spin components.

By solving for the SC order parameter self-consistently along with the Hartree potential for the same set of input parameters, we have also obtained the TSC state with the uniform electron density of 1.0 with () and a mixed TSC+TCDW state with and , both with . The three (pure TCDW, pure TSC, and mixed) states have exactly the same ground-state energy and excitation spectrum with the spectral gap . Thus, the TSC and TCDW states are degenerate ground states at half filling, as in the attractive Hubbard model for conventional -wave superconductivity.(42) Hence any linear superposition of the TSC and TCDW states is also a ground state with and the two topological orders can coexist.

By introducing surface boundaries, we have also confirmed the existence of zero-energy bound states in our self-consistent solutions. We have obtained the pure TCDW state at for , , and on a 8080 lattice with surface edges at and and PBC in the direction. A cross section of the (a) spin-up and (b) spin-down electron density is plotted as a function of at and in Fig. 10. The presence of an edge affects only few lattice sites close to the edge and we observe perfect CDW in the bulk of the system, where . The effective Zeeman field is . The excitation spectrum is presented in Fig. 11, where the abscissa is an index numbering the eigenvalues, clearly showing the existence of zero-energy states (). We have also found the pure TSC state in which the SC order parameter is enhanced significantly at and very near and , but with in the bulk and two Majorana bound states per surface with . Moreover, a mixed state has been found where the SC order parameter has the same overall structure as in the pure TSC state, but with along with in the bulk. For both pure TSC and mixed states, as for the pure TCDW state. The mixed state also hosts two zero modes per surface with .

## Iv Conclusions

In summary, by solving for the superconducting order parameter self-consistently, we have found that the effects of a nonmagnetic impurity in a 2D -wave topological superconductor can vary significantly depending on the material parameters. In particular, the weaker the Rashba SO coupling, or the stronger the Zeeman field, the more -wave-like the system reacts to a nonmagnetic impurity. The midgap excitation bound to the impurity always carries down spin at the impurity site, i.e., spin antiparallel to the direction of the Zeeman field. For relatively strong Zeeman field, however, the spin of the midgap excitation can flip as the impurity potential is varied. We have found two such cases: in one case, the midgap excitation disappears into the gap edge and reappears with opposite spin, and in the other the system undergoes a phase transition of the ground state, whereby the spin-down quasiparticle is bound to the impurity in the ground state and thus the midgap excitation becomes spin up. In this case, a nonmagnetic impurity in a 2D -wave topological superconductor acts exactly like a magnetic impurity (classical spin) in a conventional -wave superconductor. The spin flip of the quasiparticle excitation can be understood partially in terms of the TKNN number that corresponds to the shifted chemical potential at the impurity site having the opposite sign to that of the bulk TKNN number.

We have also shown by solving for each spin component of the Hartree potential as well as the order parameter self-consistently that the Hartree potential effectively reduces the Zeeman field. Furthermore, we have demonstrated in terms of our self-consistent solutions the coexistence of TCDW and TSC in the Abelian phase at half filling. This is analogous to the CDW and uniform superconducting states being the degenerate ground states at half filling in the tight-binding model for conventional -wave superconductivity: the presence of a nonmagnetic impurity, however, lifts this degeneracy in favour of CDW.(42)

In this work, we did not solve for the Hartree potential self-consistently when studying the impurity effects. When the Hartree potential is included, the TSC states presented in Sec. III.2 will require weaker pairing interaction and/or SO coupling. Further studies of the effects of nonmagnetic impurities and their influence on possible interplay of TCDW and TSC in a 2D -wave topological superconductor will be presented in a future publication.(56)

## V Acknowledgements

We thank Evan Smith for useful discussions. K.T. is grateful to CCSE at Japan Atomic Energy Agency for hospitality, where part of the research was performed. This work was enabled in part by support provided by WestGrid (www.westgrid.ca) and Compute Canada Calcul Canada (www.computecanada.ca). Some preliminary calculation was also performed on the supercomputing system SGI ICE X at the Japan Atomic Energy Agency. The research was supported by the Natural Sciences and Engineering Research Council of Canada, and partially by JSPS KAKENHI Grant No. 26800197 and the “Topological Materials Science” (No. 16H00995) KAKENHI on Innovative Areas from JSPS of Japan.

### References

1. E. I. Rashba, Sov. Phys. Solid State 2, 1109 (1960).
2. A. Manchon, H. C. Koo, J. Nitta, S. M. Frolove and R. A. Duine, Nat. Mat. 14, 871 (2015); see also references therein.
3. M. Sato, Y. Takahashi, and S. Fujimoto, Phys. Rev. Lett. 103, 020401 (2009).
4. M. Sato, Y. Takahashi, and S. Fujimoto, Phys. Rev. B 82, 134521 (2010).
5. J. Alicea, Rep. Prog. Phys. 75, 076501 (2012).
6. J. D. Sau, R. M. Lutchyn, S. Tewari, and S. Das Sarma, Phys. Rev. Lett. 104, 040502 (2010).
7. J. Alicea, Phys. Rev. B 81, 125318 (2010).
8. A. V. Matetskiy, S. Ichinokura, L. V. Bondarenko, A. Y. Tupchaya, D. V. Gruznev, A. V. Zotov, A. A. Saranin, R. Hobara, A. Takayama, and S. Hasegawa, Phys. Rev. Lett. 115, 147003 (2015).
9. L. J. Li, E. C. T. O’Farrell, K. P. Loh, G. Eda, B. Özyilmaz, and A. H. Castro Neto, Nature (London) 529, 185 (2016).
10. Y. Nagai, S. Hoshino, and Y. Ota, Phys. Rev. B 93, 220505(R) (2016).
11. S. Tewari, J. D. Sau, S. Das Sarma, Ann. Phys. 325, 219 (2010).
12. C. Nayak, S. H. Simon, A. Stern, M. Freedman, S. Das Sarma, Rev. Mod. Phys. 80, 1083 (2008).
13. T. D. Stanescu, J. D. Sau, R. M. Lutchyn, and S. Das Sarma, Phys. Rev. B 81, 241310(R) (2010).
14. A. C. Potter and P. A. Lee, Phys. Rev. B 83, 184520 (2011).
15. A. M. Black-Schaffer, Phys. Rev. B 83, 060504(R) (2011).
16. Y. Nagai, Y. Ota, and M. Machida, J. Phys. Soc. Jpn. 83, 094722 (2014).
17. K. Björnson and A. M. Black-Schaffer, Phys. Rev. B 88, 024501 (2013).
18. K. Björnson and A. M. Black-Schaffer, Phys. Rev. B 91, 214514 (2015).
19. E. D. B. Smith, K. Tanaka, Y. Nagai, Phys. Rev. B 94, 064515 (2016).
20. H. Hu, L. Jiang, H. Pu, Y. Chen, and X.-J. Liu, Phys. Rev. Lett. 110, 020401 (2013).
21. P. G. de Gennes, Superconductivity of Metals and Alloys (Westview Press, Boulder, 1999).
22. Y. Nagai, Y. Ota, and M. Machida, J. Phys. Soc. Jpn. 84, 034711 (2015).
23. J. D. Sau and E. Demler, Phys. Rev. B 88, 205402 (2013).
24. R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Phys. Rev. Lett. 105, 077001 (2010).
25. Y. Oreg, G. Refael, and F. von Oppen, Phys. Rev. Lett. 105, 177002 (2010).
26. A. Y. Kitaev, Phys.-Usp. 44, 131 (2001).
27. L. Covaci, F. M. Peeters, and M. Berciu, Phys. Rev. Lett. 105, 167006 (2010).
28. Y. Nagai, Y. Ota, and M. Machida, J. Phys. Soc. Jpn. 81, 024710 (2012).
29. T. Sakurai and H. Sugiura, J. Comput. Appl. Math. 159, 119 (2003).
30. Y. Nagai, Y. Shinohara, Y. Futamura, Y. Ota, and T. Sakurai, J. Phys. Soc. Jpn. 82, 094701 (2013).
31. D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Phys. Rev. Lett. 49, 405 (1982).
32. K. Ishikawa and T. Matsuyama, Nucl. Phys. B 280, 523 (1987).
33. G. E. Volovik, The Universe in a Helium Droplet (Oxford University Press, Oxford, 2009).
34. V. Gurarie, Phys. Rev. B 83, 085426 (2011).
35. A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Phys. Rev. B 78, 195125 (2008).
36. S. Fujimoto, Phys. Rev. B 77, 220501(R) (2008).
37. M. Sato and S. Fujimoto, Phys. Rev. B 79, 094504 (2009).
38. C. Zhang, S. Tewari, R. M. Lutchyn, and S. Das Sarma, Phys. Rev. Lett. 101, 160401 (2008).
39. A. Shitade and Y. Nagai, Phys. Rev. B 92, 024502 (2015).
40. L.-H. Wu, Q.-F. Liang, Z. Wang and X. Hu, J. Phys.: Conf. Ser. 393, 012018 (2012).
41. Y. Masaki and Y. Kato, J. Phys. Soc. Jpn. 84, 094701 (2015); J. Phys.: Conf. Ser. 568, 022028 (2014).
42. K. Tanaka and F. Marsiglio, Phys. Rev. B 62, 5345 (2000).
43. A. V. Balatsky, I. Vekhter, J.-X. Zhu, Rev. Mod. Phys. 78, 373 (2006).
44. Y. Okuno, J. Phys. Soc. Jpn. 69, 858 (2000); see also references therein.
45. M. Takigawa, M. Ichioka, K. Kuroki, Y. Tanaka, Phys. Rev. B 72, 224501 (2005).
46. P. W. Anderson, J. Phys. Chem. Solids 11, 26 (1959).
47. M. Stone and R. Roy, Phys. Rev. B 69, 184511 (2004).
48. L. Yu, Acta Phys. Sin. 21, 75 (1965).
49. H. Shiba, Prog. Theor. Phys. 40, 435 (1968).
50. A. I. Rusinov, Zh. Eksp. Teor. Fiz. 56, 2047 (1969) [Sov. Phys. JETP 29, 1101 (1969)].
51. A. Sakurai, Prog. Theor. Phys. 44, 1472 (1970).
52. M. I. Salkola, A. V. Balatsky, J. R. Schrieffer, Phys. Rev. B 55, 12648 (1997).
53. M. E. Flatté and J. M. Byers, Phys. Rev. Lett. 78, 3761 (1997).
54. M. E. Flatté, J. M. Byers, Phys. Rev. B 56, 11213 (1997).
55. A. Weiße, G. Wellein, A. Alvermann and H. Fehske, Rev. Mod. Phys. 78, 275 (2006).
56. S. L. Goertzen, Y. Nagai, E. D. B. Smith, K. Tanaka (unpublished).
101841