Multifarious topological quantum phase transitions in two-dimensional topological superconductors

We study the two-dimensional topological superconductors of spinless fermions in a checkerboard-lattice Chern-insulator model. With the short-range p-wave superconducting pairing, multifarious topological quantum phase transitions have been found and several phases with high Chern numbers have been observed. We have established a rich phase diagram for these topological superconducting states. A finite-size checkerboard-lattice cylinder with a harmonic trap potential has been further investigated. Based upon the self-consistent numerical calculations of the Bogoliubov-de Gennes equations, various phase transitions have also been identified at different regions of the system. Multiple pairs of Majorana fermions are found to be well-separated and localized at the phase boundaries between the phases characterized by different Chern numbers.

Topological superconductors (TSCs) have been a hot topic in recent a few years due to their unique properties and potential applications 1 . One of the most significant characteristics of TSCs is Majorana fermions (MFs) which are found to exist in vortices [2][3][4][5] or on the boundaries 6-15 of the TSC system. MFs are their own anti-particles 16,17 . They obey non-Abelian statistics 4 and could be used in topological quantum computation 18,19 , which is also a significant reason for the recent wide concern on TSCs.
Generally, two proposals are suggested to realize TSCs and MFs, either with the spin-triplet p-wave superconductivity 6,7,[20][21][22] or with the conventional s-wave superconductivity by proximity effect in some materials [8][9][10][11][12][13][14][15][23][24][25][26][27][28][29][30][31][32] . The simplest and most straightforward model systems are the one-dimensional spinless p-wave TSCs or two-dimensional p x ± ip y TSCs. For example, Read and Green 20 considered a chiral TSC state with a Chern number 33 = 1  for p x ± ip y paired spinless fermions. In a one-dimensional quantum wire model by Kitaev 6 , MFs were found to locate at the ends of the chain. Possible topological quantum phase transitions (TQPTs) among TSCs and conventional superconducting or insulating states are also very intriguing, and it is possible to identify MFs at the phase boundaries of TQPTs. However, previous works 6,7,[20][21][22] are mostly done based upon models with topological pairing of fermions in a single topologically trivial band and varieties of TSC phases seem more or less limited, and additional condition such as long-range superconducting pairing 21,22 is required in some cases to achieve other TSC phases. Recently, Qi et al. proposed that when a quantum Hall (QH) state or a quantum anomalous Hall (QAH) state near the plateau transition of N to N − 1 is coupled to a conventional s-wave superconductor through the proximity effect, a new TSC phase with a Chern number = − N 2 1  will appear between the phase  = N 2 and 2N − 2. Their work illustrates that a nontrivial TSC phase with a Chern number = 1  or 2 is achievable through TQPTs from trivial phases 30 .
In this paper, we address a two-dimensional TSC model of paired spinless fermions (or spin-polarized electrons) on a Chern insulator (CI) with QAH states and focus on the TQPTs that will possibly occur. The checkerboard-lattice CI/QAH model which has two nontrivial topological energy bands is our preferred candidate system. We only consider the short-range nearest-neighbor superconducting pairing. The possible phase diagram has been explored, and rich and interesting TQPTs have been identified. Various TSC phases with the Here, c i † c ( ) i annihilates (creates) a spinless fermion on the site i. Staggered fluxes are superimposed in the plaquettes and induce additional phase factors ± φ on the nearest-neighbor hopping t. This means the C 4 rotational symmetry is broken, leading to two sublattices A (red) and B (blue). The second-nearest-neighbor hopping takes different value ′ t ij ( ′ t 1 , and ′ t 2 for solid, and dotted connection, respectively). μ is the chemical potential. ∆ = − 〈 〉 V c c ij j i is the pairing potential with V the strength of the attractive interaction. We set t = 1 as unit and ′ = − ′ = ′ = . t t t 0 5 1 2 in the calculations. In fact, our main results are not sensitive to the selected parameters. In the normal state, the system is a CI/QAH state with two topological bands of Chern numbers N = ± 1 when φ ≠ nπ 35 . Particularly, one of the two bands can be tuned to be very flat with the third-nearest-neighbor hopping, which hosts new fractional quantum Hall states 36,37 . TQPTs and phase diagram of TSCs. In order to explore the possible TQPTs, we solve Eq. (2) (See Methods) with different parameters at fixed p x + ip y pairing order parameter Δ = 0.1. The Hamiltonian can be diagonalized under a rational transformation and further a Bogoliubov transformation. In principle, the distinct TSC phases can be characterized by the Chern numbers of two lower BdG bands (C 1 , C 2 ) (also two higher BdG bands with C 3 = − C 2 and C 4 = − C 1 due to particle-hole symmetry). However, the individual Chern number is not well defined in some cases (especially for small |μ|) due to overlap of the two lower BdG bands. We use the sum of the two Chern numbers  = + C C 1 2 to characterize the TSC phases. The Chern number of each BdG band can be calculated with those BdG wave functions numerically 33 . A sophisticated way is to find out the condition under which the two middle BdG bands touch and re-open. The phase boundaries are μ = ± 4 sin φ with a single Dirac point (π, π), μ = ± 4 cos φ with a single Dirac point (0, 0), and µ = ± ′ − ∆ t 16( ) 8 2 2 with two Dirac points (π, 0) and (0, π). Numerical calculations show that each Dirac point carries a Berry phase ± π and the total Berry phase is ± 2π for the case of two symmetric Dirac points with the same sign Berry phases. Therefore, when the system crosses a phase boundary, a TQPT occurs with the changed Chern number  ∆ = ±1, or ± 2 for single, and double Dirac points, respectively. The overall phase diagram of TSCs with the characteristic Chern numbers  ranging from − 3 to 3 is depicted in Fig. 2, showing multifarious TQPTs. In fact, for small enough μ (uncolored region), the system is the same as the single band p + ip superconductor with strong pairing suggested by Read and Green 20 irrespective of the topology of the two normal-state bands. Here, the high Chern number TSC is naturally produced without Figure 1. The checkerboard-lattice CI/QAH model. Sites in red and blue constitute two sublattices respectively. Staggered fluxes are superimposed on the plaquettes, resulting in an additional phase factor ± φ on the nearest-neighbor hopping (± denoted by the arrow direction). Solid and dotted lines, connecting the second-nearest-neighbor sites (same sublattice), denote the hopping parameters ′ t 1 and ′ t 2 .
Scientific RepoRts | 6:28471 | DOI: 10.1038/srep28471 applying additional condition, for example the long-range superconducting pairing 21,22 . The topological phase transition in present CI with + p i p x y superconductivity is much richer than that in the QH/QAH states coupled to a conventional s-wave superconductor through the proximity effect 30,31 . Besides the previous suggested  = 1 TSC phase between the  = 2 TSC phase and = 0  superconducting phase, TSC phases with  = 3 (or − 3) also emerge via a series of TQPTs.
Phase diagrams with ∆ < 2 /2 are almost the same except that the phase boundaries of µ = ± ′ − ∆ t 16( ) 8 2 2 (green lines in Fig. 2) have slight shifts. These two phase boundaries disappear when ∆ > 2 /2, resulting in the limited TQPTs with Chern numbers  ranging from − 1 to 1. We will not illustrate these cases in detail here. We also remark that the corresponding phase diagram for p x − ip y pairing with Δ = 0.1 is almost the same as Both the x and y lattice directions are treated with periodic boundary conditions. Due to the puny size of the lattice in y direction compared with the large size in x direction, in the numerical calculations, we will neglect the y dependence and only consider the x dependence of physical quantities. The trap center is located at L 0 = 200. By solving the BdG equations self-consistently, some interesting results are revealed below.
We select the parameters φ = 0.42π, V = 1.75, μ = − 0.9, U trap = 0.00018 to study the TQPTs and the corresponding MFs. The fermion density n i [red line in Fig. 3(a)] decreases monotonically to zero from the center to the left/right edge, which is well confined by the trap potential. In comparison, the superconducting order parameter Δ (green line) exhibits interesting site dependence. Its amplitude increases from a relative smaller value at the center site (L i = 200), and reaches a maximum value at a critical length (L i = 90, 310). The superconducting gap then weakens down to zero upon further moving outwards. The non-monotonic behavior can be understood as follows. The energy scales of the two normal-state bands are ascertained to , resulting in the near zero value of Δ . Here, we emphasize that the superconducting gap Δ near the edge in the p + ip TSC is always zero regardless whether the trap potential is applied due to lack of counterpart edge state. This differs from the on-site s-wave pairing in HgTe quantum well 40 and QAH state with a conventional s-wave superconductivity proximity effect 30 , where the strongest superconductivity locates near the edge. , respectively. The TSC Chern numbers  are specified. Green arrows mark the slightly shifted positions of the phase boundaries in green for a fixed Δ = 0.35. Black dots with serial numbers will be discussed in the next subsection.
In fact, the possible TSC phase transitions are well illustrated within the present trap potential. From the center to the left/right edge, the system undergoes the TSC phases with  = −2, − 3, − 1 respectively, and then the insulating phase  = 0 with decreasing µ i eff , consisting with the notations in Fig. 2. This process can also be signified by the first-order derivatives of the superconducting order parameter and the fermion density [ Fig. 3(b)]. A few peaks can be identified, corresponding to the above mentioned TSC phase transitions.
In general, the MFs emerge near the phase boundary between the different TSC or TSC and insulating phases, as well as near the edge of the TSC phases. In Fig. 3(c), we output all 6400 eigenvalues of this cylinder system. The middle energy spectrum and four pairs of zero-energy modes protected by an energy gap about 0.06 can be clearly observed in the insert. As we know, for the E = 0, one fermion can be separated into two MFs which are expressed as γ = + † c c ( )/ 2 with c † the zero-energy fermion 16,17 . According to zero modes, we plot spatial distributions of the corresponding four pairs of MFs in Fig. 4. Two MFs of each pair are well-separated and locally distributed on the lattice. Distributions of all the MFs are summarized in Fig. 4(e) where a subtraction |ρ i − ρ i+1 | (i = 1, 3, 5, 7) is adopted to eliminate the overlapping between the MFs in the same pair. The locations of the MF peaks in Fig. 4(e) are the same as those of the peaks in Fig. 3(b), further manifesting the TQPTs in TSC and MFs emerging near the phase boundaries. Moreover, the effective chemical potential where the MF peaks locate also corresponds to the phase boundaries μ = − 4 cos φ, µ = − ′ − ∆ t 16( ) 8 2 2 and μ = − 4 sin φ as we discussed in the above section. The TQPTs and MFs can also be verified by the zero-energy LDOS 11 along the lattice sites, where six peaks exist with the positions same as those of the MFs, and peaks of the

Summary and Discussion
In conclusion, we illustrate the multifarious TQPTs of two-dimensional spinless fermions with p-wave superconducting pairing on a checkerboard-lattice CI/QAH model. Various TSC phases, especially with high Chern numbers, are established with just short-range nearest-neighbor pairing. A rich TSC phase diagram is revealed by tuning the chemical potential and staggered-flux phase factor. This is in sharp contrast to the previous theoretical works based on single-band p + ip TSC or QAH states coupling to a conventional s-wave superconductor through the proximity effect. Furthermore, with the recent development on generating artificial gauge fields 41 and especially realizing artificial staggered fluxes in optical lattices [42][43][44][45] , such a checkerboard-lattice CI/QAH model could be potentially realized in the near future. A finite-size checkerboard-lattice cylinder with a confining harmonic potential trap has been further explored using the self-consistent BdG method. Well-separated multiple pairs of MFs emerge at topological phase boundaries, identifying the TQPTs. These MFs are well controllable and robust against perturbations, and might be used for implementing the non-Abelian statistics and topological quantum computation 18,19 . In addition, such a checkerboard-lattice model might carry a topological flat band at selected model parameters which is beneficial to the superconducting pairing due to the large density of states. Future work might be very interesting to explore the competition of these TSC phases with other possible phases involved with the topological flat band.

Methods
For the system with the periodic boundary conditions, the Hamiltonian can be rewritten in momentum space as z x y . σ and I are the Pauli matrix and unit matrix, respectively. The p x + ip y pairing in momentum space is with Δ the superconducting order parameter. When the lattice system has a finite size L x × L y , we solve the BdG Hamiltonian self-consistently in real space where , ij is the single particle Hamiltonian with τ and τ′ vectors linked to the nearest-neighbor and second-nearest neighbor sites. u v ( , ) i n i n T is the quasiparticle wave function corresponding to the eigenvalue E n . Due to the particle-hole symmetry of the BdG equations, the wave vector ⁎ ⁎ v u ( , ) i n i n T is also an eigenvector corresponding to eigenvalue − E n . The superconducting order parameter, on-site particle number n i , and the LDOS ρ i (ω), are determined self-consistently by ij n i n j n n B ∑ = n u f , i n i n n 2 respectively with f n as the Fermi distribution function.