# From fractionally charged solitons to Majorana bound states in a one-dimensional interacting model

###### Abstract

We consider one-dimensional topological insulators hosting fractionally charged midgap states in the presence and absence of induced superconductivity pairing. Under the protection of a discrete symmetry, relating positive and negative energy states, the solitonic midgap states remain pinned at zero energy when superconducting correlations are induced by proximity effect. When the superconducting pairing dominates the initial insulating gap, Majorana fermion phases develop for a class of insulators. As a concrete example, we study the Creutz model with induced -wave superconductivity and repulsive Hubbard-type interactions. For a finite wire, without interactions, the solitonic modes originating from the nonsuperconducting model survive at zero energy, revealing a fourfold-degenerate ground state. However, interactions break the aforementioned discrete symmetry and completely remove this degeneracy, thereby producing a unique ground state which is characterized by a topological bulk invariant with respect to the product of fermion parity and bond inversion. In contrast, the Majorana edge modes are globally robust to interactions. Moreover, the parameter range for which a topological Majorana phase is stabilized expands when increasing the repulsive Hubbard interaction. The topological phase diagram of the interacting model is obtained using a combination of mean-field theory and density matrix renormalization group techniques.

## I Introduction

In the context of quantum field theory, Jackiw and Rebbi Jackiw and Rebbi (1976) introduced a general mechanism to generate zero modes with fractional charges. They considered a one-dimensional massless fermion coupled to a bosonic scalar field. Owing to symmetry breaking, the bosonic field can acquire a finite uniform expectation value , which provides a mass to the fermion. When the scalar field has a kink interpolating between two opposite values of the mass, the one-particle Dirac equation exhibits a nondegenerate solution which is isolated in the middle of the gap (at zero energy ), and spatially localized around the mass inversion point where vanishes. This nondegenerate solution is protected by a unitary symmetry that connects each single-electron state with energy to its partner located at the opposite energy. The zero-energy state is self-conjugated under this discrete symmetry and thereby protected. In the many-particle description, there are two degenerate many-body ground states corresponding to the state being empty or filled, while negative energy states of the Fermi sea are all filled. Besides, the existence of a U(1) symmetry ensures that the electric charge is a good quantum number. Then, the ground state with empty (occupied) carries the charge (), in units of the original fermion charge. This has been generalized to arbitrary fractional charges Su and Schrieffer (1981); Goldstone and Wilczek (1981) and more complicated bosonic kinks.Jackiw and Rossi (1981); Volovik (2003) The first condensed matter realization of the Jackiw-Rebbi mechanism came with the study of conducting polymers, such as polyacetylene. Su et al. (1979); Heeger et al. (1988)

More recently, another type of zero-energy excitation, the Majorana fermion, has attracted a lot of attention from the condensed matter community.Alicea (2012); *Beenakker2013 A Majorana fermion is its own antiparticle,Majorana (1937) and can appear at defects of topological superconductors, such as -wave superconductorsRead and Green (2000) or superconducting hybrid systems mimicking them.Fu and Kane (2008); Lutchyn et al. (2010); Oreg et al. (2010) In presence of superconductivity, the U(1) symmetry is broken, such that the electrical charge is no longer conserved. Moreover, a discrete antiunitary particle-hole symmetry emerges which differs from the one discussed above in the case of fractionally charged solitons. Nevertheless, a Majorana excitation can also be described within the Jackiw-Rebbi scheme and its existence is signaled by a nondegenerate and self-conjugated solution of the Bogoliubov–de Gennes (BdG) equations. Read and Green (2000) The crucial difference appears in the many-particle description of the system: two solutions of the BdG equations are needed to construct a fermionic doublet.

This paper aims at understanding the physics of systems that host fractionally charged solitons in their normal state and Majorana modes in some of their superconducting phases. In particular, we study how the zero-energy solitonic modes survive under the addition of proximity-induced superconductivity and/or Hubbard repulsive interactions. Several phases are a priori possible, including a fully gapped superconducting phase without midgap states, a topological superconductor with Majorana end states, and finally a superconducting phase with midgap states which are not Majorana states. The general phase diagram of such a system is expected to be ruled by at least two symmetries: the chiral symmetry of the normal system, and the particle-hole symmetry associated with superconductivity. Although the U(1) symmetry is broken, zero-energy states can still be protected by the discrete symmetry between the positive and negative energies. The Majorana bound states (MBS) are end states described by particle-hole self-conjugated wave functions. The solitonic modes which do not have this Majorana self-conjugation property (and eventually survive in the superconducting phase) are henceforth called, interchangeably, chiral bound states (CBS) or, simply, solitons. Recently, the coexistence of CBS and MBS was studied in Rashba nanowires where the localized end states are realized by applying a non uniform magnetic field.Klinovaja et al. (2012) In this nanowire system, the end states occur usually at finite energy and a fine tuning of the parameters of the model is required to bring these states to zero energy. The BDI symmetry classAltland and Zirnbauer (1997) can provide a family of models that exhibit both CBS and MBS, over a wide range of their phase diagram, without fine tuning. Indeed, the BDI class harbors both topological insulators and superconductors associated with integer topological numbers ().Schnyder et al. (2008); Kitaev (2001) On one hand, this implies that BDI superconductors host multiple Majorana particles which are spatially close but do not hybridize, thereby remaining at zero energy.Niu et al. (2012); Tewari and Sau (2012); Sticlet et al. (2013a) On the other hand, the BDI class also encompasses insulating (nonsuperconducting) systems with zero-energy fractionally charged solitonic end states. The most famous representative example of such BDI nontrivial insulators is the polyacetylene chain described by the Su-Schrieffer-Heeger (SSH) model.Su et al. (1979)

Among the possible BDI systems, this paper investigates in detail the Creutz model,Creutz and Horváth (1994); Creutz (1999) and also some aspects of the SSH modelSu et al. (1979) in the Appendix. Initially introduced in the context of lattice quantum chromodynamics, the Creutz model has gained a foothold in condensed matter physics as a versatile model to test different ideas, such as the topological production of defects when crossing a quantum transition point,Bermudez et al. (2009) the dynamics of Dirac points under interaction,Dóra et al. (2013) the decay of the edge states in the presence of a thermal bath,Viyuela et al. (2012) or the persistent currents in Dirac fermion rings.Sticlet et al. (2013b) Proposals for the realization of the Creutz model in cold-atom systems have been recently advanced.Mazza et al. (2012) Moreover, phases of the bosonic Creutz model were also recently investigated.Takayoshi et al. (2013); Tovmasyan et al. (2013)

In the absence of superconductivity, the Creutz model exhibits gapped phases with and without solitonic modes.Creutz and Horváth (1994); Creutz (1999) The transition between these inequivalent gapped phases occurs by closing the bulk gap, thereby leading to a semi-metallic phase with a single Dirac cone at the transition. Our purpose is to enrich the Creutz system by adding spin-singlet superconducting pairing, induced by a proximity effect. When the superconducting order parameter becomes larger than the initial insulating gap, a topological superconductor (TSC) phase develops, hosting Majorana edge modes. In view of this, the superconducting Creutz model is henceforth called the Creutz-Majorana (CM) model. This phenomenology is also derived in the SSH model in the presence of proximity-induced superconducting pairing. Moreover, the superconducting phase in which the solitons survive proves to be in the bulk a topological phase protected by spatial inversion at a bond and fermion parity. It is topologically distinct from the Majorana phase and from the trivial phase without in-gap states.

The second thread followed by this paper concerns the effect of the repulsive Hubbard interaction on the different phases of the CM model and on its edge states. The effect of Coulomb interactions on MBS has been investigated recently in nanowires with strong spin-orbit in the presence of superconducting proximity effect which are relevant for experiments.Mourik et al. (2012) Such nanowire models Oreg et al. (2010); Lutchyn et al. (2010) break the time-reversal symmetry, but they nevertheless map onto a BDI system at low energy, the Kitaev Majorana chain,Kitaev (2001) and it was shown theoretically that Majorana bound states are robust with respect to repulsive Hubbard interactions.Gangadharaiah et al. (2011a)

Soon after, Stoudenmire et al. Stoudenmire et al. (2011) established that MBS are not only robust to interactions, but that, counterintuitively, the parameter regime where the MBS can be observed increases, at the price of a reduction in the bulk superconducting gap. Following the methodology of Ref. Stoudenmire et al., 2011, we study the effect of a repulsive Hubbard interaction on the Creutz-Majorana model using a combination of mean-field theory and density matrix renormalization group (DMRG). Globally, we obtain comparable results for the Majorana phase of the Creutz model: its parameter range grows upon increasing the Hubbard interaction. Nevertheless, there is a region in parameter space, where interactions prove detrimental to the Majorana states (cf. Fig. 8). In contrast to nanowire models, the models studied in this paper (Creutz and SSH models) host CBS in their normal state in a wide range of parameters. The Hubbard interaction breaks a chiral symmetry that protected the CBS modes, pushing them to finite energy, while the MBS modes remain pinned at zero energy. Even at finite interactions, the phases surrounding the TSC have a different topological character from each other. We describe these phases in terms of their entanglement properties and excitations.

Finally, we notice that recent papers discuss the possibility of having time-reversal-invariant topological superconductors supporting two MBS at each end of a one-dimensional (1D) system.Zhang et al. (2013); Keselman et al. (2013); Klinovaja and Loss (2013) These Kramers pairs of MBS are protected by coexisting time-reversal symmetry and chiral symmetry. We emphasize that the Creutz model belongs to a completely different symmetry class (BDI rather than DIII) and has at best one MBS at each end. Nevertheless, there are also two symmetries at play, the pseudo-time-reversal symmetry and the chiral symmetry, that allow us to understand the behavior of the solitonic edge modes and of the MBS.

In brief, this paper investigates the mechanism of competition between fractionally charged solitons and Majorana modes in the large class of BDI models, thereby complementing the recent proposal in Ref. Klinovaja and Loss, 2013. Furthermore, we demonstrate that the stabilization of Majoranas by repulsive interactions, first obtained in Ref. Stoudenmire et al., 2011, is also globally observed in a distinct topological insulator, the Creutz model, although its phase diagram also presents regions where the Majoranas are destabilized.

The paper is organized as follows: Section II reviews the Creutz model and its topological properties. Section III investigates the Creutz-Majorana model, including the survival of the CBS in the superconducting region, and their transition to MBS upon increasing the pairing. Section IV takes into account the effects of repulsive interactions in the Creutz-Majorana-Hubbard model (CMH) using self-consistent Hartree-Fock approximation and density matrix renormalization group (DMRG). Section V examines a particular case of the interacting, but nonsuperconducting, Creutz model where a bulk metallic phase develops. Section VI holds the conclusions. The Appendix comments on the generality of the results by treating the noninteracting Su-Schrieffer-Heeger (SSH) polyacetylene modelSu et al. (1979) enriched with a -wave superconducting pairing.

## Ii Creutz model

Here, we present the noninteracting part of the lattice model. Historically, this type of model was introduced by WilsonWilson (1977) to solve the fermion doubling problem and simulate Dirac fermions on lattices. Later, Creutz investigated the existence of edge states in finite chains or ladders.Creutz and Horváth (1994) We shall refer to this model simply as the Creutz model in the remainder of this paper. This section provides a review of the spectral and topological properties of the Creutz model.

### ii.1 Model

In the absence of superconductivity and interactions, our starting point is the lattice Hamiltonian

(1) |

where the sum runs over all sites indexed by (see illustration in Fig. 1). In this paper, the electronic spin is represented by the standard Pauli matrices (). The spin indices for the electron annihilation operators and spin matrices are implicit. The electrons can jump from one site to the nearest-neighboring site while conserving their spin: this process has a complex amplitude , i.e., the electrons gain or lose a phase when hopping between the same spin states. The electrons can also hop between sites with amplitude while flipping their spin, which mimics a spin-orbit coupling. Finally, there is an onsite mass term which favors the polarization of the electronic spin along the direction. This Hamiltonian can be seen as an adaptation from Ref. Creutz, 1999, where the chain index of the Creutz ladder has been replaced by the electronic spin.

The Hamiltonian is diagonalized by Fourier transforming it into momentum space :

(2) |

Consequently, there are two bands with the energy dispersion

(3) |

shown in Fig. 1. From Eq. (3) it follows that the two energy bands can touch either at , or at . At these momenta, the energy dispersion is linear and there is a Dirac cone band touching. For , the band-touching takes place at , and for , at . When both and vanish, the system exhibits two Dirac cones.

When and are finite and close in absolute value, the low-energy physics is given at or , corresponding to one gapped fermion. Without loss of generality, the parameters in the Creutz model () can be considered positive. In this case, the Dirac cone forms at , and it is gapped with a mass , which can be either positive or negative. If becomes larger than , new minima (maxima) for the upper (lower) band in Eq. (3) develop near .

### ii.2 Topological characterization

The Creutz model is classified in the BDI class of topological insulators.Schnyder et al. (2008); Kitaev (2001) Although it is a normal (nonsuperconducting) system, it is possible to find an antiunitary operator , which anticommutes with the Hamiltonian, and thus realizes a particle-hole symmetry . The operator is the complex-conjugation operator. There is also a pseudo-time-reversal symmetry represented by antiunitary operator, , such that . Both pseudo-time-reversal and particle-hole symmetry operators square to one. Finally, there is a chiral symmetry represented by , which is proportional to the product of the other two symmetries. The chiral symmetry anticommutes with the Hamiltonian

(4) |

thereby allowing us to associate any eigenstate with energy and momentum to another eigenstate with opposite energy and the same momentum.

When (resp. ), the system is gapless with a Dirac cone at momentum (resp. ). A particular point is the case of vanishing and , when the system recovers the time-reversal symmetry, and exhibits two Dirac cones at and . In all these gapless phases, there is no well-defined global topological winding number.

For all other values of the parameters , the chain is a fully gapped insulator, whose topological properties are encoded in its associated winding number

(5) |

From the expression of the winding number it follows that when , the system is a topologically trivial insulator () and when , it is a topologically nontrivial insulator characterized by .Sticlet et al. (2013b)

In the topologically nontrivial state there are chiral zero-energy bound states (CBS) at the edges of an open system. These localized states are protected by the chiral symmetry and the bulk gap. For a finite chain, the two edge CBS overlap through the insulating bulk. This overlap is exponentially small for chains longer than the spatial extension of the CBS.Creutz (1999)

## Iii Creutz-Majorana model

A natural question is whether Majorana fermions can be induced in the Creutz model by considering a superconducting pairing potential. This would not be very surprising because the system has the necessary ingredients to realize Majorana fermions. Crucially the system has the possibility to realize a single Dirac cone in the Brillouin zone (for ). At the Dirac point, the electron wave functions are eigenstates of with zero energy. Then a spin-singlet superconducting pairing could open a gap in the spectrum and Cooper pairs of spin-polarized fermions can be formed. In this way, -wave pairing plus a topological insulating state could generate a -wave superconducting regime where MBS become possible at the ends of a finite system. A similar mechanism was thoroughly investigated in different proposals to realize MBS using the proximity effect to an -wave superconductor.Fu and Kane (2008); Lutchyn et al. (2010); Oreg et al. (2010)

In addition, the Creutz model displays CBS in the inverted gap regime, (where all the parameters in the model were considered real positive numbers). Therefore a key question is whether induced superconductivity immediately removes these modes, and leads to the formation of Majorana fermions. A central finding of this section is that the solitonic (CBS) modes are robust with respect to the superconducting pairing. These modes are protected by the chiral symmetry and they will subsist at finite superconducting pairing. It will be necessary that the bulk gap closes at finite , in order to subsequently have Majorana fermions forming (see Fig. 2).

It is noteworthy to remark that the CBS will prove fragile in comparison to the Majorana fermions. The reason is that the CBS crucially depend on the chiral symmetry for staying at zero energy, and they are removed even by scalar potentials that are created by nonmagnetic disorder. In the strong pairing regime (large , with respect to the other system parameters), it is expected that the topological Majorana phase is destroyed Read and Green (2000).

### iii.1 Bulk structure

#### iii.1.1 Model

Let us modify the Creutz Hamiltonian (1) by adding an -wave singlet superconducting pairing with amplitude . This leads to the Creutz-Majorana model described by the Hamiltonian

(6) |

where the sum is taken over all the lattice sites, and the order parameter can be considered real without any loss of generality.

The momentum-space Hamiltonian is quadratic in the standard BdG form. Let us choose the basis . In this basis, the Hamiltonian is written as with

(7) |

The are the spin Pauli matrices and are Pauli matrices in particle-hole space. The products of two Pauli matrices from different spaces is understood as a tensor product.

By diagonalizing the BdG Hamiltonian , it follows that there are four energy bands satisfying

(8) |

There are four possible gap closings in the system at , for , and at , for . As it will be shown, all these lines mark topological transitions between different gapped topological phases.

#### iii.1.2 Symmetries

In the presence of translational invariance, the following symmetries catalog the Creutz-Majorana Hamiltonian in the BDI topological class.Altland and Zirnbauer (1997); Schnyder et al. (2008); Kitaev (2001) There are two symmetries which are represented by anti-unitary operators squaring to one. First, the Bogoliubov–de Gennes (BdG) Hamiltonian has an intrinsic built-in particle-hole symmetry

(9) |

represented by . This anti-unitary operator anticommutes with the BdG Hamiltonian and squares to one.

Second, the system also has a pseudo-time-reversal symmetry represented by the anti-unitary operator satisfying

(10) |

and commuting with the Hamiltonian.

Finally, there is a chiral symmetry proportional to the product of the two symmetries and :

(11) |

This symmetry associates to each state with energy , a partner at , with the same momentum. These three symmetries are sufficient to place the CM model in the BDI class.

Nevertheless, it is crucial to note a second set of symmetries that also place the model in the same BDI class. There is a second time-reversal symmetry, , and a second chiral symmetry, :

(12) |

The new chiral symmetry can be considered as inherited from the original Creutz model, since the nonsuperconducting system already exhibits a chiral symmetry represented by . Therefore, under the BdG doubling of degrees of freedom, the symmetry can be extended to the hole sector, such that an electron of a definite chirality is reflected in a hole of the same chirality. Hence, the addition of superconductivity does not destroy this chiral symmetry, but it simply promotes it in the CM model to .

Nevertheless, there is an important difference between the two chiral symmetries of the CM model, and . The latter, , remains (as in the Creutz model) fragile in the sense that it can be destroyed by a scalar (nonmagnetic) on-site potential (which multiplies a matrix) because . In contrast, the chiral symmetry is preserved as it anticommutes with the scalar potential, . Therefore, unless the bulk gap closes or a perturbation breaks the chiral symmetry, it is guaranteed that the zero-energy modes of the topological insulating phase of the Creutz model remain protected in the Creutz-Majorana model.

The fragility of the chiral symmetry to (local or global) changes in the chemical potential will prove central in explaining the phenomenology of the topological phases in the Creutz-Majorana model with and without interactions.

#### iii.1.3 Phase diagram

A one-dimensional BdG Hamiltonian from the BDI class always has an associated topological invariant which is a winding number.Schnyder et al. (2008) Using either one of the two chiral symmetries, it is possible to block off-diagonalize the CM model and to obtain

(13) |

The topological invariant can be extracted from the matrix by defining a winding number as

(14) |

This procedure can be further simplified by using both chiral symmetries to completely off-diagonalize the Hamiltonian. Subsequently, the winding number can be determined algebraically.

In the present case of the Creutz model, the winding number does not exceed 1 in absolute value, so one can expect at most a single MBS at a given edge. This allows us to use a perhaps more transparent determination of the invariant by using the Majorana number .Kitaev (2001) This invariant is sensitive only to the winding number parity. But since the Creutz-Majorana model does not realize higher winding number phases, , can decide unequivocally whether the system is in a topological trivial or non-trivial state. In the trivial phase, without Majorana states, the invariant is , while in the nontrivial phase (topological superconductor), .

Following Ref. Kitaev, 2001, the Majorana number is defined at the time-reversal-invariant momenta, and

(15) |

where Pf denotes the Pfaffian of a matrix. It follows immediately that the topological invariant reads as

(16) |

where we have assumed a nonzero on-site .

The ensuing phase diagram is illustrated in two different representations in Figs. 2 and 3. Figure 2 uses Eq. (16) to present the phases as a function of system parameters squared, such that the relative sign difference is not an issue. Both the phase without in-gap states (gray) and the CBS hosting phase (yellow) are trivial, with respect to the Majorana number (). As we shall see in Sec. IV.2, these two phases are actually distinct with regards to a different topological invariant. This compact picture can be compared with Fig. 3, where it is shown more clearly that the topological transition lines in the phase diagram are in fact the lines where the bulk spectrum is gapless [Eq. (8]. Figure 3 refines Figs. 2 by labeling the respective topological phases, and also by illustrating the phase diagram for the original Creutz model ().

Note in both figures the effect of a finite superconducting pairing on the original Creutz model, which lies in the horizontal axis . When , an infinitesimal is sufficient to open a superconducting gap at the Dirac cone, and produces a topologically nontrivial superconductor. However, when the Creutz model is deep in a nontrivial insulating state with solitons at its end, the addition of superconducting pairing does not destroy these modes. Since the CBS are protected by the bulk gap and the chiral symmetry , they survive the induced superconducting pairing.

Because the solitons develop in a gapped phase without Majorana fermions, the topological invariants and remain impervious to the existence of the CBS phases. They can not differentiate between trivial phases and the phase with CBS states. As a consequence, the winding number of the Creutz model is not a limit case of the winding number in CM model, :

(17) |

Since the chiral bound states’ existence can only be inferred from the bulk properties, it is necessary to study their localization on defects in an infinite system, or at the edges of an open system.

### iii.2 From chiral bound states to Majorana fermions

Here, we further analyze the transition between the trivial superconductor with CBS and the topological superconductor with MBS. Let us first consider a single interface which can host either CBS or MBS. Using a continuum approximation of the lattice model, the single interface problem is solved analytically and the BdG wave functions are derived explicitly. In the second part of this section, we perform numerical diagonalization of a finite length system, which hosts CBS or MBS at its two end points. These studies yield a more precise account of the transition that takes place from two chiral bound states (in the redoubled BdG formalism) to a single Majorana edge mode that develops in the nontrivial superconducting phase.

#### iii.2.1 Edge state analysis at an interface in a continuum model

Let us consider an isolated interface between a trivial Creutz insulator and the superconducting Creutz-Majorana model (CM). Depending of the topological phase of CM, there will be zero-energy bound states at this interface. Then, by increasing the strength of the superconducting pairing , there is a transition induced from the region with CBS to the one with MBS (see Fig. 3). To further fix ideas, let us consider all the parameters in the original models () as positive real numbers. In this case, the low-energy physics, when and are comparable in magnitude, is given by that of a gapped Dirac fermion at , with mass . The sign change in the mass signals a transition from a topologically trivial () to a nontrivial phase (“inverted” gap; ) (Eq. 5).

Let us start with the Creutz model on an infinite line, and an additional spatial twist is introduced in the mass of the Dirac fermion

(18) |

where is a positive number, and is the spatial coordinate (Fig. 4). According to Eq. (5), the trivial phase is realized to the left (), and the nontrival phase, to the right (). At the interface there will be a single solitonic mode (CBS), eigenstate of the chiral symmetry operator .

Now, let us investigate the more interesting case of the Creutz-Majorana model, when a superconducting pairing is added to the nontrivial insulating phase. This will allow us to observe the effect of the superconducting gap on the bound state from the topological insulator. It is now possible to increase the pairing until the lower boundary to the Majorana phase is crossed (see Fig. 3). This provides a picture for the behavior of the edge states over the topological transition.

The continuum Hamiltonian that models the interface between the trivial Creutz insulator and the CM model is readily obtained by linearizing the lattice Hamiltonian near and neglecting second-order contributions in momentum

(19) |

with and . The function is the Heaviside step function. The Fermi velocity is determined from the lattice model , with the lattice constant.

It is possible to write analytically the wave functions at positive and negative , respectively. The solution is obtained by matching these wave functions at the interface . Moreover, the problem is simplified by looking specifically for zero-energy, bound solutions at the interface ( and ).

In the trivial insulating state to the left of the interface, the wave function reads as

(20) |

where and are space-independent coefficients, and transposes the row four-vectors. In the superconducting phase to the right, the solution will depend on whether it is the superconducting gap or the Dirac fermion mass which dominates. For , the right-side wave function reads as

(21) |

Matching the wave functions over the boundary, it follows that the coefficients of depend on and

(22) |

Crucially, the matching procedure has left two free coefficients, and . This indicates that there are two modes living at the interface, and they are eigenstates of the chiral symmetry operator . These modes are diagonal in electron-hole space, so one of them is electron-like while the other is a hole-like excitation (see Eq. 20). In the doubled (four-band) BdG formalism, they correspond to the solitonic mode that was already present in the Creutz model, in the absence of the superconducting pairing . Breaking the chiral symmetry that protects them allows removal of the CBS from zero energy. To sum up, when the “inverted” Dirac fermion mass dominates the superconducting gap, there are two modes at zero energy protected by the chiral symmetry inherited from the Creutz model.

When the two masses are equal (), the bulk gap closes. This signals the topological transition to the Majorana phase. Subsequently, for , the superconducting gap dominates, and the solution in Eq. (21) becomes non-normalizable at . The appropriate (normalizable) wave function in this new regime reads as

(23) |

Matching the right, , and left, , wave functions reveals that a single mode remains at the edge in the regime:

(24) |

At the wave-function level, one notices how the two CBS are combined into a Majorana fermion. Moreover, the resulting wave function has the property that it is self-particle-hole conjugate (for ), indicating that one true Majorana fermion lives now at the interface. As dictated by the bulk-boundary correspondence,Volovik (2003) the edge-state existence could be predicted from information about the bulk topology.

Adding a constant potential in the Hamiltonian, removes the CBS from zero energy, but it does not destroy the Majorana fermion. For a small with respect to both the superconducting gap and the insulating gap, it follows that the domain for Majorana fermions increases, such that the transition from the trivially gapped state to the nontrivial superconductor moves at lower pairing potential . Additionally, if the chemical potential is larger than the insulating gap , then any positive pairing potential opens a gap which hosts Majorana fermions. This is a case where the potential is in a band and the single-band electrons are subsequently coupled by a -wave-type pairing to realize a topological superconducting phase.

In conclusion, we have described the transition between the trivial superconducting phase with doubly-degenerate zero-energy states protected by the chiral symmetry to a phase with Majorana fermions. This section offers the proof that the local twofold degeneracy in the spectrum does not stand for multiple Majorana modes at an edge, but for chiral bound states protected by a fragile symmetry inherited from the Creutz insulator.

#### iii.2.2 Spectrum of the finite lattice model: Gapping the CBS

The CBS in the superconducting phase contiguous to the solitonic phase of the Creutz model, can be equally explored using exact diagonalization of the lattice Hamiltonian in Eq. (6). In an open wire, the CBS and the MBS form at the two edges of the system. Therefore there will be a doubling in the degeneracy of zero-energy modes. The numerical analysis of the lattice confirms the location of the bound states, reveals the fourfold energy degeneracy of the CBS, and their sensitivity to chiral symmetry breaking terms.

It is possible to take a crosscut through the phase diagram from Fig. 3. The system parameters are fixed, except the pairing which is varied in order to traverse the three distinct regimes of the model: trivial superconductor, superconductor with CBS, and nontrivial superconductor with MBS. In units where , the parameters are and . The pairing is varied as to encompass the topological transitions at and . The lowest eigenvalues resulting from exact diagonalization are plotted against in Fig. 5 (top). The chiral bound states are seen as fourfold-degenerate midgap states at small . The bulk spectrum closes at and a twofold degeneracy characteristic of Majorana fermions remains. In the strong-pairing regime , the system becomes a trivial superconductor devoid of midgap states.

The fragility of the chiral bound states can be tested by adding different types of scalar potentials. Because the chiral symmetry does not anticommute with the local on-site scalar potentials, these modes are not robust. Nonmagnetic disorder or even a constant chemical potential removes them from zero energy. In fact, it is sufficient to add a local potential on the first and last sites to lift them from zero energy. This chiral-symmetry-breaking term is implemented in Eq. (6) by adding

(25) |

The result is presented in the spectrum of the lowest eigenstates in Fig. 5(bottom), for and , deep in the CBS region of the phase diagram. Since the chemical potential is chosen identical on the first and last sites , the fourfold degeneracy is reduced to two degenerate pairs of states moving in opposite energy direction away from zero energy.

In the Creutz model, a chemical potential can move two degenerate states away from zero energy. In the CM superconducting model, a pair of electronic states is reflected in a pair of hole states that feel an opposite chemical potential. This observation accounts for the evolution of states seen in Fig. 5 (top). The remaining twofold degeneracy is naturally removed by rendering incommensurate the chemical potential magnitude on the first and last sites. That is the general case of a (nonmagnetically) disordered wire.

In contrast, in the nontrivial superconducting region there is a single Majorana fermion bound at each edge. Moderate disorder which does not close the bulk gap leaves them at zero energy.

Both kinds of edge modes are eigenstates of the chiral symmetry operator . Two CBS, the positive eigenstates are localized at one edge, while the negative eigenstates, at the other edge. If is a zero-energy eigenstate, one can represent the spatial distribution for the chiral operator expectation value . If is normalized to one, then the following equation holds:

(26) |

depending whether is a positive or negative eigenstate. The sum is carried over all lattice sites. This yields at once a representation for the localization of the modes, as well as a proof of them being chiral. The result is shown in Fig. 6. Note that the local weight of the expectation value between two eigenstates at two different edges is not correlated. There is a degree of freedom in distributing this weight for degenerate states, resting on the specifics of the diagonalization procedure. Nevertheless, the global Eq. (26) always holds.

In conclusion, we have found a trivially topological phase with fragile midgap edge states which are eigenstates of the chiral operator. These states are removed from zero energy by chiral symmetry breaking perturbation. Also they can disappear by driving the system in a Majorana phase.

### iii.3 Many-body ground-state degeneracy

Until now, we have only described the zero-energy edge modes from the perspective of single-particle physics, without touching on the question of the state filling. This subsection clarifies the way to count the ground states expected in the Creutz-Majorana model at the transition from fractionally charged solitons, to the CBS in the superconducting phase, and, subsequently, to Majorana bound states.

In the absence of superconductivity, the Creutz model provides nontrivial insulating phases with zero-energy states. These states can be either filled or empty with electrons. For long wires and at half filling, there will be no energy cost in having either of the states filled or empty. Therefore, while there are two states in the spectrum at zero energy, the ground state will be fourfold degenerate, with four possible configurations in the occupation of the zero modes (Fig. 7).

The states at positive energy can be thought as states carrying positive charges. Due to the chiral symmetry, each filled negative energy state has an empty positive energy partner. Therefore, when there is a single electron occupying one of the two midgap states, the overall charge of the system is zero. When both zero-energy states are filled, there will be an excess of one electron charge in the system. Because the charge of a midgap state can change only by one, it follows that each edge state will carry a fractional charge . Similarly, when both states are empty, the total energy does not change, and each edge state carries a charge .Heeger et al. (1988)

Moreover, the states are eigenstates of the chiral operator . Here, the left modes have spin projection up and the right states, spin projection down. Therefore, when both states are either filled or empty, the spin of the zero-energy pair is zero. When only the one state is filled, the total charge vanishes, but the ground state will have spin .

Under the addition of superconductivity, the charge symmetry is broken. However the system maintains its chiral symmetry and there will be CBS edge modes at zero energy. In the BdG formalism, the number of states is doubled. Therefore, in the BdG spectra of the CM model, there will be four zero-energy end states. Nevertheless, the degeneracy of the many-body ground state does not change. There are still four configurations, because the quasiparticle occupation numbers are not independent between particle-hole symmetric states.

While the states remain at zero energy, the ground-state is fourfold degenerate, either in the nontrivial insulating Creutz model, or in the CBS phases of the superconducting model. Breaking the chiral symmetry, allows states to be moved away from zero energy, while remaining localized in bulk energy gaps, as it was shown in different nanowire systems.Gangadharaiah et al. (2011b); Klinovaja et al. (2012) Nevertheless, the ground state becomes in this case unique as for the trivially gapped system.

When crossing the topological transition to the Majorana phase, the BdG spectrum will reveal the existence of two Majorana solutions bound at zero energy at the ends of the wire. The ground state is also twofold degenerate. Two pair of states can be either filled or empty with a single electron. In contrast to the fractionally charged solitons, the electron itself can be viewed as divided between the edge states of the system in the form of two Majoranas. This is pictorially represented in Fig. 7.

The upshot of this section is that the topological transition from CBS to Majorana fermions can be tracked either by looking at the spectrum of the BdG Hamiltonians, or from the many-body ground-state degeneracy. These two view points must be clearly distinguished in the following, when interactions are treated either in mean-field theory, or by using DMRG techniques. While in mean field one has access to BdG Hamiltonians, and to their spectral properties, the DMRG obtains the many-body ground state.

## Iv Creutz-Majorana-Hubbard Model

The Creutz-Majorana model on a finite-size chain hosts either zero-energy solitonic (CBS) states, or Majorana modes, depending on the strength of the superconducting pairing . The aim of this section is to investigate the effects of repulsive onsite interactions on these edge modes and to obtain the topological phase diagram of the model. To that end, we consider the following Hamiltonian of the Creutz-Majorana-Hubbard (CMH) model:

(27) |

The first term, , represents the Hamiltonian of the Creutz-Majorana model from Eq. (6), and contains the Hubbard interaction between onsite electronic densities of opposite spin

(28) |

with being the interaction strength.

In this section, the CMH model is studied using a combination of self-consistent Hartree-Fock theory and extensive DMRG simulations. Let us begin by summarizing the main results and defer the details to the following subsections.

Figure 8 shows the topological phase diagram for different values of interaction , combining results from the mean-field analysis and DMRG simulations. The agreement between both approaches is overall very good. The main standout feature of the phase diagram is that the topological phase and its associated Majorana bound states at zero energy are robust to interactions. In fact, the parameter range for which a topological Majorana phase is stabilized globally expands, upon increasing the Hubbard coupling. However, there is also a region in the phase diagram (), where the Hubbard interaction is detrimental to the Majoranas. Here their existence domain decreases with increasing interaction strength. The mean-field phase diagram captures the topological transitions of the model rather accurately for small and moderate interactions. At strong interactions, the mean field keeps a very good estimate of the topological transition at small , while at large it tends to overestimate the extension of the Majorana phase. The “spin-orbit” coupling tends to delocalize the electrons and leads to quantum fluctuations in the particle number. This explains, at a qualitative level, the divergence of the mean field results from the “exact” DMRG results at large “spin-orbit” coupling.

A priori, this overall very good agreement may seem surprising in a quantum 1D system, where fluctuations are expected to be pronounced. However, let us recall that U(1) charge and the SU(2) spin-rotation symmetries are broken in this model, and hence the system acquires a certain stiffness to fluctuations, as observed in nanowire models Stoudenmire et al. (2011) which also break these symmetries.

Crucially, the interactions have a different effect on the two types of zero-energy states localized at the edges of the wire. The chiral bound states from the large phase (cf. Fig. 2) are removed from zero energy. Recall that for , there is a fourfold zero-energy degeneracy in the BdG energy spectrum as well as a fourfold degenerate many-body ground state. At the mean-field level, upon increasing the interaction strength, the zero-energy CBS acquire a finite energy. Further increase of the interaction leads to a decrease of the bulk gap and pushes the end states into the bulk continuum. (Fig. 9). As soon as the CBS aquire a finite energy, the many-body ground-state degeneracy is lifted. DMRG results confirm that the ground state becomes unique, corresponding in the BdG picture to having all negative quasiparticle states filled (Fig. 10). Therefore, the degeneracy of the ground state in the CBS phase is not protected with respect to Hubbard interactions. However, the system remains in a “nontrivial” phase which is protected by a product of fermion parity and inversion symmetry. This phase is characterized by a bulk topological invariant.

In contrast, the Majorana end states generally remain pinned at zero energy as expected from the phase diagram. Consequently, the ground state of the system remains twofold degenerate. As already remarked, only in a limited region with a small superconducting gap () and small “spin-orbit” coupling () can Majorana fermions be destroyed by interactions (Fig. 8).

### iv.1 Mean-field study

At the mean-field level, the four-operator interaction is approximated by a sum of two-operator interactions. In the present case, we choose the following Hartree-Fock decoupling of density-density interaction term from Eq. (28):

where the site index was neglected, since the decoupling involves only local terms at a given site. The brackets denote the expectation value in the ground state. The site-independent parameters and are defined:

(30) |

A finite anomalous pairing indicates intrinsic superconducting correlation produced by repulsive interactions . Note that the Hamiltonian , already has an induced (extrinsic) superconducting pairing , which is independent on the interactions . The finite represents the tendency towards a polarization in the direction.

Let us also introduce the onsite electronic density and the onsite magnetization (in the direction),

(31) |

where again these quantities are uniform. Consequently, the ground state of the system is characterized by the expectation values . The Wick theorem is verified by evaluating both sides of Eq. (IV.1) in the ground state.

Because the specific decoupling in Eq. (IV.1) is site independent, the implicit assumption that we have made is that the ground state can develop only ferromagnetic instabilities. However, the DMRG simulations do not reveal evidence of antiferromagnetic instabilities, so in this approximation we neglect antiferromagnetic order parameters or more complicated nonlocal decoupling terms. A further simplification in the present mean-field theory is to consider only real order parameters.

Let us consider the mean-field Hamiltonian as describing a periodic chain, where the ground state is unique. Therefore, a Fourier transform allows us to determine the momentum-space mean-field Hamiltonian :

(32) |

The interaction energy reads as

(33) |

and the BdG Hamiltonian ,

(34) | |||||

where it becomes clear that interactions renormalize the onsite terms. The exchange term modifies the polarization and the anomalous pairing renormalizes the proximity-induced pairing . In addition, the interaction also introduces a chemical potential and a -polarization term . The term would act as a Zeeman term, breaking the time-reversal symmetry.

The mean-field Hamiltonian (34) is quadratic and can be diagonalized in a basis of quasiparticles as:

(35) |

with , and the Hamiltonian matrix, . The quasiparticle operators create quasiparticle excitations at positive energies . Also, note that the spectrum exhibits the particle-hole symmetry with a spectrum symmetric about the zero energy.

The system ground state is annihilated by the quasiparticle excitations. Therefore, it can be represented as

(36) |

where is the vacuum of conventional fermionic annihilation and creation operators. The usual creation and annihilation operators can be decomposed on the basis of the quasiparticle excitations with coefficients given by the eigenvectors of the Hamiltonian .

The order parameters for the mean-field Hamiltonian are found using a self-consistent procedure. At the start, the order parameters are assumed to be zero. After diagonalizing the mean-field Hamiltonian, the order parameters are evaluated in the ground state given by Eq. (36). This allows us to determine the interaction energy , and, consequently, the total Hartree-Fock energy, given by the occupied quasiparticle energy levels, plus the interaction energy (). The resulting order parameters are fed back into the mean-field Hamiltonian and a new ground state is obtained. Therefore, new order parameters can be calculated. The cycle of operations is repeated until self-consistency is reached (i.e., the ground-state energy settles to a minimum value among cycles). As a consequence, this procedure generates a final set of order parameters which completely characterize the mean-field Hamiltonian.

In the mean-field calculations, we consider wires of lengths at several interaction strengths. The self-consistent procedure is performed over a point grid in the parameter space to obtain the order parameters of the mean-field Hamiltonian.

Crucially, the self-consistent analysis (backed by DMRG results) reveals that there is no polarization developing in the direction. The vanishing Zeeman term implies that the mean-field Hamiltonian is invariant. Consequently, it can still be classified in the BDI class, as a noninteracting and translation-invariant Hamiltonian, which can be treated similarly to the Creutz-Majorana model in Sec. III. Therefore, information about its topological phases can be extracted from the topological invariant:

(37) |

where the renormalized on-site terms read as

(38) |

The topological phase diagram from Fig. 8 follows by plotting the transition lines where the Majorana number changes sign.

The properties of the edge states due to interactions can be explored at mean-field level. Retaining the interaction-renormalized order parameters obtained for the periodic system, one can study wire models with open edges. Hence, the existence of MBS and CBS can be tested in comparison with further DMRG results (Sec. IV.3).

At the mean-field level, interactions produce a finite uniform scalar potential , which breaks the chiral symmetry . Therefore, the zero-energy CBS inherited from the Creutz model are no longer protected. As shown in Sec. III.2, this leads to a removal of the zero-energy degeneracy for CBS states. Fig. 9 shows the lifting of the CBS modes from zero energy under increasing interaction strength .

### iv.2 DMRG Study

In this section, the DMRG study of the ground-state phases and phase transitions of the Creutz-Majorana-Hubbard model is discussed. We first employ the DMRG method White (1992) in an MPS formulation Schollwöck (2011) to obtain the three lowest-lying levels of the energy spectrum for open chains of length . This is done successively, by first obtaining the ground state in MPS form, and then performing a new DMRG calculation optimizing an MPS which is orthogonal to the previous ones.Schollwöck (2011) The resulting spectra, shown in Fig. 10, illustrate the fate of the MBS and CBS when interactions are considered. The results are consistent with mean-field theory of the open chain. When starting from a phase with CBS modes, interactions remove the fourfold degeneracy and yield a unique ground state (bottom-right panel in Fig. 10). It corresponds to having all negative quasiparticle states occupied. In contrast, the twofold degeneracy of the ground state in the Majorana phase is robust to interactions.

The above procedure on open chains is, however, rather cumbersome, especially when close to the critical points where the correlation length rises and finite-size effects from the coupling of edge modes can become relevant. For this reason, we employ the more efficient iDMRG algorithm McCulloch (2008); Schollwöck (2011); Kjäll et al. (2013) to obtain an MPS representation of the ground-state wave function for an infinite system (see Ref. Kjäll et al., 2013 for the details of the algorithm). We conserve the fermion number parity in our simulations, and use a maximum bond dimension up to when calculating ground states. Obtaining the energy spectrum of an infinite chain is problematic. However, the entanglement spectrum is readily available and also contains information about the degeneracies of the system.Li and Haldane (2008); Pollmann et al. (2010); Turner et al. (2011); Fidkowski and Kitaev (2011)

An MPS representation of a wave function can be cut anywhere in the chain (including across a physical site), yielding the Schmidt decomposition

(39) |

where represent the states to the left/right respectively of the cut. Taking , these are the eigenstates of the reduced density matrix for the right half of the system, . The bulk entanglement spectrum,

(40) |

reflects how the different Schmidt states contribute to the wave function . It can be interpreted as the spectrum of the so-called entanglement Hamiltonian which contains information about “artificial edges” of semi-infinite chains. The symmetries that protect the topological degeneracies in the energy spectrum associated with open edges also protect degeneracies in the entanglement spectrum.Pollmann et al. (2010); Turner et al. (2011); Fidkowski and Kitaev (2011) In the following, we work in units where .

The different phases found in Fig. 11 for can be distinguished by their entanglement spectrum, whose features survive the addition of interactions. We perform two different cuts of the chain, as illustrated in the insets of Fig. 11 : one of the cuts corresponds to separating two physical sites; and the other can be seen as the act of polarizing a single site.

The topological phase with Majorana end states (blue), is distinguished from the other phases by its twofold entanglement degeneracy for both kinds of cuts. This degeneracy is due to the fractionalization of the fermion parity symmetry and can not be lifted unless the system undergoes a phase transition, as discussed in detail in Ref. Turner et al., 2011.

The trivial phase at small is only degenerate for cuts inside a physical site. The phase at large (yellow) displays entanglement degeneracy only for a vertical cut, separating two physical sites. This degeneracy is protected by the product of the fermion parity and inversion symmetry, i.e., by with being the inversion of the chain at a bond. In order to distinguish this phase from the trivial phase at small , we employ a method similar to the one used to classify topological phases in spin chains.Pollmann and Turner (2012) The basic idea is to find the cohomology class to which the “fractionalized” representation of a symmetry belongs. The matrices are representations of the symmetry acting on the Schmidt states . In our case, assuming a two-site unit cell, we use the symmetry to calculate . We find that the topological invariant and thus the phase is a topological phase, while in the trivial phase for small we find . A direct consequence of is the observed even degeneracy in the entanglement spectrum (see Ref. Pollmann and Turner, 2012 for details). The spatial inversion symmetry is broken at the edges of an open chain and thus there will be no signatures of this phase in the energy spectrum of Fig. 10. The existence of the CBS end states at is stabilized by the presence of further symmetries, which are broken by interactions (as discussed above).

The boundaries between different phases correspond to gapless points in the thermodynamic limit. We extract the correlation length of the system from the spectrum of the transfer matrix of the MPS. A sharp increase near the boundaries between phases is observed in Fig. 12 (top), both with and without interactions. This reflects the diverging behavior of the correlation length close to the critical point, cutoff by the finite bond dimension used, as discussed in Ref. [Kjäll et al., 2013]. We therefore take the location of the maxima as an estimate for the position of the critical points.

The transition between the TSC phase and the topological phase stabilized by at , and in Fig. 12(bottom row) is a representative case of all the phase transitions reported here. In order to analyze the phase transitions, we employ a finite-entanglement scaling approach Pirvu et al. (2012); Kjäll et al. (2013) of the entanglement entropy and correlation length . The entanglement entropy, , also displays a divergent behavior when approaching the critical point (not shown). The spin-orbit term is first fine-tuned in order to give the largest possible for . This provides an estimate of and ensures that the system is in the finite-entanglement scaling region. At this point, the entanglement entropy scales with the calculated correlation length, setting an effective length scaleCalabrese and Cardy (2004),

(41) |

giving access to the central charge of the critical point. A standard scaling form for the correlation length is also used,

(42) |

The extracted central charge and correlation-length exponent are in very good agreement with those of the 2D Ising universality class, and . All the transitions reported between the different phases, with and without interactions, belong to the Ising universality class.

In order to connect with the mean-field results of Sec. IV.1, we analyze in Fig. 13 the average fermion number per site and associated fluctuations. For and (bottom), the fermion number deviates from half filling in order to reduce the Hubbard energy for on-site repulsions. In the mean-field description, this corresponds to an effective (uniform) chemical potential, that lifts the chiral symmetry which protected the zero-energy chiral bound states. Additionally, the increasing fluctuations of the fermion number with imply that the mean-field approach of Section IV.1 becomes less accurate for large .

The topological phase with MBS is extremely robust to strong interactions, as shown by the persistence of the twofold degeneracy of the entanglement spectrum in Fig 14(top) for , up to . The correlation length increases with increasing (bottom panel), but it is always found to converge to a finite value. Hence, the bulk energy gap remains finite, albeit reduced by increasing interactions. These conclusions are in good agreement with those reached by Stoudenmire et al. Stoudenmire et al. (2011) for the effects of Hubbard interactions on a topological phase hosted by a nanowire model. We have additionally checked that this phase is robust to a moderate chemical potential randomly distributed across a finite chain (not shown).

### iv.3 Spectral properties

The excitation properties of topological phases provide important clues to their nature, as well as insight into quantities which can be used experimentally to measure them. Alicea (2012); Beenakker (2013) For example, the spectral functions can be used to track the closing of the bulk energy gap, which may signal a topological transition. Moreover, the local spectral function describes the local single-particle density of states at the edges of an open chain. This information could be used to detect Majorana fermions via a quantized zero-bias signal in the differential conductance.Fu and Kane (2009); Akhmerov et al. (2009); Law et al. (2009); Flensberg (2010) In this section, the spectral functions are calculated within mean-field theory and using MPS techniques.

#### iv.3.1 Excitation spectrum

At the mean-field level, the excitation spectrum is obtained by solving for the eigenenergies of the BdG Hamiltonian from Eq. (34) on a periodic chain. The self-consistent Hartree-Fock calculation results are represented in comparison with the DMRG results in Fig. 15.

Spectral functions can be obtained with an MPS approach via the evolution of correlations in real time with the time-evolving block decimation (TEBD) method.Vidal (2004, 2007); White and Feiguin (2004); White and Affleck (2008); Schollwöck (2011) A fermionic operator or (the spin index is suppressed for clarity) is applied at site to the ground state MPS of a finite chain. The resulting excited MPS is evolved in real time, such that the hole and particle Green’s functions, defined as

(43) | |||

(44) |

are calculated separately. A Fourier transform to momentum and frequency space gives the retarded single-particle Green’s functionMahan (2000) =. The momentum-resolved spectral function is given by

(45) |

which probes the energy spectrum associated with single-particle excitations in the bulk. We let the matrix bond dimensions grow with time such that the truncation error is at most per step. A fourth-order Suzuki-Trotter decomposition with is employed, where the time scale is set by the inverse hopping . The maximum time is chosen such that the wave-front of the “light cone” of correlations does not reach the edges of the system, resulting in a cutoff at small frequencies. The sampled time is extended using linear predictionWhite and Affleck (2008) to improve the resolution in the limit: .

The bulk spectral functions for the different phases found are shown in Fig. 15. The color map represents the TEBD results for open chains of sizes up to and the black dashed lines are calculated within mean-field theory. The overall qualitative agreement is rather good. The single-particle spectra in Figs. 15 (a), (c), and (e) correspond exactly to the respective low-energy spectra of the noninteracting model at half filling, and the gapped nature of all phases is evident. The flat bands in Figs. 15(c) arise because there is no momentum dependence in a band in the parameter point (Eq. 8).

Turning on strong interactions (right column of Fig. 15) results in breaking of the particle-hole symmetry, seen in a corresponding shift of the spectral weight towards positive energies. Nevertheless, the mean-field Hamiltonian remains symmetric [Eq. (9)], because it captures only the single-particle physics and retains the symmetry between the positive and negative states. The band crossings of Figs. 15(c) and (e) are lost, with a subgap opening at small momentum. Increasing leads to a decreasing of the bulk gap, cf. Fig. 14, such that it quickly becomes too small to observe with our energy resolution. Inside the topological superconducting phase Fig. 15(d), the gap is not located at , but moves at some small momenta which depend on and . Also noteworthy is the appearance of a faint shoulder in for close to in Fig. 15(d) at , which may signal the presence of a scattering continuum. Unsurprisingly, a local probe in the bulk, such as , does not distinguish between phases with different topological character.

#### iv.3.2 Local spectral functions

We now turn our attention to the edge physics. In order to look at edge properties, we calculate the local spectral function for open chains of size for the same parameters as in Fig. 15. The edge spectral function is represented in Fig. 16 using both mean-field results and DMRG approaches. Once again, the qualitative agreement between the two methods is rather good.

At the mean-field level, the Hamiltonian for the open chain is still described by Eq. (34). However, the order parameters () on the open chain are determined from the Hartree-Fock calculation on the periodic system. The local density of states at site reads as

(46) |

The sum is carried over all eigenstates with respective eigenenergies . The Dirac delta function is numerically approximated using Gaussians of width inversely proportional to the system size.

In the context of DMRG, the local spectral functions are computed by adding a fermion (or removing) to an edge of the chain and then calculating only equal-space correlations, i.e., in Eqs. (43) and (44).

The signatures of the trivial phase are not affected much by the presence of interactions, [cf. Figs. 16(a) and 16(b)]. The absence of zero-energy quasiparticle peaks remains, but the energy gap decreases with interactions.

More interestingly, the topological phase shows a clear signature of the Majorana zero-energy mode for both and in Figs. 16(c) and 16(d). This can be interpreted as a direct observation of the Majorana bound state. Without interactions, this mode is well separated by a symmetric gap of width from the remainder of the excitations. The overall structure survives the addition of interactions, reflecting once again the robustness of the Majorana states. The survival of the zero-energy MBS peak was also recently observed in the Kitaev chain subject to nearest-neighbor repulsive interactions.Thomale et al. (2013)

The edge CBS at large are also clearly signalled by the presence of a zero-energy peak in the noninteracting regime [Fig. 16(e)]. In order to perform this simulation, the ground-state MPS is biased to one of the four degenerate states by applying a small pinning potential to an edge of the chain. Turning interactions to splits this peak into two peaks at nonzero energies, showing that the CBS have moved away from zero energy, as described in the preceding sections.

## V Creutz-Hubbard model

In this section, we focus on the original Creutz model without superconductivity, but with Hubbard interactions

(47) |

In the noninteracting Creutz model, there are two insulating phases meeting at a critical gapless point (Fig. 2). In contrast, at finite , there is an extended phase developing around along the nonsuperconducting line (Fig. 8). The present section seeks to elucidate the nature of this phase.

Close to the Dirac point (), adding a small interaction leads to moving the chemical potential from half filling into the bulk bands. Therefore, the system develops a metallic phase. However, when the Creutz model is in an insulating phase, small interactions with respect to this gap, keep the chemical potential inside the gap and leave the system insulating. This explanation accounts for the observed increase in the metallic phase with the interaction strength, near the original band touching.

At the mean-field level, we consider explicitly in this case possible antiferromagnetic ordering in and, respectively, directions

(48) |

However, both mean-field and DMRG indicate that there is only ferromagnetism in the direction developing for low to moderate (and vanishing ferromagnetic , and antiferromagnetic , ). Keeping the formalism where the degrees of freedom are doubled, the metallic region is identified with a closing of an energy gap. On Fig. 17 it is shown the system becoming metallic at precisely the predicted points on the phase diagram (Fig. 8).

This picture is fully supported by DMRG results (cf. Fig. 18 for a representative case and ). The system is in a gapless metallic phase, with a central charge as for free fermions. Superconducting -wave correlations are observed to decay exponentially fast with distance.

The interaction reinforces the term. For positive , is negative such that the renormalized on-site coupling is larger than the bare . Nevertheless, the “spin-orbit” term competes with the ferromagnetic ordering . It tends to delocalize the electrons and to spin-flip, thus leading to a reduction of the polarization at higher values of (Fig. 17).

At finite superconducting pairing , the metallic region is immediately gapped out. Consequently, Majorana fermions are hosted in the proximity-induced gap, leading to the topological superconducting region in Fig. 8.

## Vi Conclusions

This paper demonstrates that insulators hosting fractionally charged states (chiral bound states or CBS) continue to exhibit solitonic modes under the addition of proximity-induced superconducting pairing. It also explains in detail how the CBS can be transformed into Majorana bound states (MBS) when the pairing exceeds a certain threshold fixed by the gap of the insulator. Building on two representative examples, the Creutz model and the Su-Schrieffer-Heeger model (treated in the Appendix), we propose the following general mechanism for one-dimensional topological insulators in the BDI class.

For moderate pairing, the CBS remained pinned at zero energy owing to the protection of chiral symmetries. This is reflected in a fourfold-degenerate many-body ground state for a finite wire. Indeed, the two zero-energy end states of the wire can be either occupied or empty without changing the total energy of the electronic system. However, breaking the chiral symmetries, pushes the CBS at finite energy, thereby leading to a unique ground state (the chemical potential is fixed at zero energy) because the localized end states are then either completely filled or completely empty.

At stronger superconducting pairing, larger than the initial insulating gap, the Creutz, and the Su-Schrieffer-Heeger insulators are driven in a phase supporting Majorana fermions. In the absence of interaction, we have investigated this transition directly by finding the wave functions of the zero mode pinned at the interface between a topological and a trivial semi-infinite phase. This transition can be also seen in a halving of the many-body ground-state degeneracy for a finite wire: the two end states become then two spatially separated Majoranas forming a single ordinary fermion (topological quantum bit).

The behavior of the CBS and MBS is also considered under the effects of repulsive interactions. In the Creutz model, the interactions were shown to break the chiral symmetry, leading to a removal of the CBS from zero energy and a unique ground state. In contrast, the Majorana fermions are rather robust, such that the parameter regime for the existence of Majorana fermions increases, albeit the bulk gap is reduced by interactions.

Aside from the edge-state physics, we have also studied the bulk phase diagram of the Creutz-Majorana-Hubbard model using a combination of self-consistent Hartree-Fock theory and extensive DMRG simulations (cf. Fig. 8), the agreement between both approaches being overall very good. The main feature is that increasing the Hubbard repulsion tends to expand stability domain of the topological Majorana superconductor, in agreement with the phenomenology obtained in other nanowire models.Stoudenmire et al. (2011) In contrast, at weak pairing and small spin-orbit coupling , the interactions lead to a decrease of the Majorana phase (cf. Fig. 8). For large , the model exhibits a topological phase which is protected by a combination of fermion parity and inversion symmetry.

Besides the topological phases, the phase diagram (see Figure 8) presents various other interesting bulk phases. For example, in the absence of superconductivity (), a metallic phase develops for a finite range of values (which extends upon raising the interaction ) of the spin-orbit parameter owing to the effect of interactions. This gapless phase was restricted to in the absence of interactions. It would be interesting to study the possibility of intrinsic superconductivity in such a phase.

The authors thank E. Orignac and S. Huber for interesting discussions. D. S. thanks P. Simon for drawing his attention on the subject of fractionally charged solitons. J. C. acknowledges support from EU/FP7 under contract TEMSSOC. D. S. and J. C. were supported by the French ANR through projects ISOTOP and MASH.

*

## Appendix A Superconducting SSH model

The physics of survival of solitonic modes into the superconducting phase and their subsequent transformation into Majorana fermions is more general. This section comments on the topological phases in another model from the BDI class, the Su-Schrieffer-Heeger (SSH) model.Su et al. (1979)

The most simple superconducting SSH model is represented by a dimerized chain with a cell containing two atoms A and B:

(49) |

The alternating hopping strengths represent the “single bond”, , and the “double bond”, , of the polyacetylene chain. This model captures only the fermionic degrees of freedom from the Hamiltonian in Ref. Su et al., 1979, neglecting the bosonic degrees of freedom due to lattice distortion. Crucially, the model is enriched here by the addition of a homogeneous superconducting order parameter which couples near-neighbor sites (Fig. 19).

Let us consider first a finite one-dimensional wire in the absence of superconductivity () with an integer number of cells () and positive real values for and . If the inter-cell hopping strengths are larger than the intra-cell hopping terms, , then a zero-energy state is bound at the edge.

As in the CM model, the solitonic modes survive at zero energy under the addition of a finite superconducting pairing. In this case, they are protected by the insulating gap induced by the Peierls instability . The general mechanism for their survival depends on the relative magnitude between the insulating gap and the superconducting gap. The topology of the Hamiltonian is characterized using Eq. (16) and the resulting phase diagram is represented in Fig. 20.

For small superconducting pairing , the system has two phases without Majorana fermions. These are two distinct phases: one is an inversion symmetry-protected phase (only in the bulk), and with solitonic modes (CBS) that survive at finite until entering the Majorana region (); the other is a trivial phase devoid of midgap states ().

For large superconducting pairing , the system is in a topologically nontrivial phase which supports Majorana fermions at the ends. Note that Kitaev chain modelKitaev (2001) at zero chemical potential is recovered by making the system uniform in normal hopping, . In this case, any finite yields one Majorana bound state at the edge.

In the presence of repulsive interactions, we expect to see a similar removal of the ground-state degeneracy in the CBS phase. This was already investigated in the absence of superconductivity, revealing a split in energy between a charged soliton and the neutral solitons.Kivelson and Heim (1982) The addition of superconducting pairing should qualitatively change these results only by allowing a Majorana phase robust to interactions.

In this section, we have shown how another BDI system supporting fractionally charge modes can be driven in a Majorana phases under the effect of superconducting pairing. In a straight parallel to the Creutz model, the solitonic modes are protected by a chiral symmetry and remain at zero energy in a wide region of parameters, before entering a Majorana phase. The condition to see the transition to the Majorana phase comes to a competition between the insulating gap and the induced superconducting pairing.

## References

- Jackiw and Rebbi (1976) R. Jackiw and C. Rebbi, Phys. Rev. D 13, 3398 (1976).
- Su and Schrieffer (1981) W. P. Su and J. R. Schrieffer, Phys. Rev. Lett. 46, 738 (1981).
- Goldstone and Wilczek (1981) J. Goldstone and F. Wilczek, Phys. Rev. Lett. 47, 986 (1981).
- Jackiw and Rossi (1981) R. Jackiw and P. Rossi, Nuclear Physics B 190, 681 (1981).
- Volovik (2003) G. E. Volovik, The Universe in a Helium Droplet (The International Series of Monographs on Physics, 117) (Oxford University Press, USA, 2003).
- Su et al. (1979) W. P. Su, J. R. Schrieffer, and A. J. Heeger, Physical Review Letters 42, 1698 (1979).
- Heeger et al. (1988) A. J. Heeger, S. Kivelson, J. R. Schrieffer, and W. P. Su, Rev. Mod. Phys. 60, 781 (1988).
- Alicea (2012) J. Alicea, Reports on Progress in Physics 75, 076501 (2012).
- Beenakker (2013) C. Beenakker, Annual Review of Condensed Matter Physics 4, 113 (2013).
- Majorana (1937) E. Majorana, Il Nuovo Cimento (1924-1942) 14, 171 (1937).
- Read and Green (2000) N. Read and D. Green, Phys. Rev. B 61, 10267 (2000).
- Fu and Kane (2008) L. Fu and C. L. Kane, Phys. Rev. Lett. 100, 096407 (2008).
- Lutchyn et al. (2010) R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Phys. Rev. Lett. 105, 077001 (2010).
- Oreg et al. (2010) Y. Oreg, G. Refael, and F. von Oppen, Phys. Rev. Lett. 105, 177002 (2010).
- Klinovaja et al. (2012) J. Klinovaja, P. Stano, and D. Loss, Phys. Rev. Lett. 109, 236801 (2012).
- Altland and Zirnbauer (1997) A. Altland and M. R. Zirnbauer, Phys. Rev. B 55, 1142 (1997).
- Schnyder et al. (2008) A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Phys. Rev. B 78, 195125 (2008).
- Kitaev (2001) A. Y. Kitaev, Physics-Uspekhi 44, 131 (2001), arXiv:cond-mat/0010440 .
- Niu et al. (2012) Y. Niu, S. B. Chung, C.-H. Hsu, I. Mandal, S. Raghu, and S. Chakravarty, Phys. Rev. B 85, 035110 (2012).
- Tewari and Sau (2012) S. Tewari and J. D. Sau, Physical Review Letters 109, 150408 (2012).
- Sticlet et al. (2013a) D. Sticlet, C. Bena, and P. Simon, Phys. Rev. B 87, 104509 (2013a).
- Creutz and Horváth (1994) M. Creutz and I. Horváth, Phys. Rev. D 50, 2297 (1994).
- Creutz (1999) M. Creutz, Phys. Rev. Lett. 83, 2636 (1999).
- Bermudez et al. (2009) A. Bermudez, D. Patanè, L. Amico, and M. A. Martin-Delgado, Phys. Rev. Lett. 102, 135702 (2009).
- Dóra et al. (2013) B. Dóra, I. F. Herbut, and R. Moessner, Phys. Rev. B 88, 075126 (2013).
- Viyuela et al. (2012) O. Viyuela, A. Rivas, and M. A. Martin-Delgado, Phys. Rev. B 86, 155140 (2012).
- Sticlet et al. (2013b) D. Sticlet, B. Dóra, and J. Cayssol, Phys. Rev. B 88, 205401 (2013b).
- Mazza et al. (2012) L. Mazza, A. Bermudez, N. Goldman, M. Rizzi, M. A. Martin-Delgado, and M. Lewenstein, New Journal of Physics 14, 015007 (2012).
- Takayoshi et al. (2013) S. Takayoshi, H. Katsura, N. Watanabe, and H. Aoki, Phys. Rev. A 88, 063613 (2013).
- Tovmasyan et al. (2013) M. Tovmasyan, E. P. L. van Nieuwenburg, and S. D. Huber, Phys. Rev. B 88, 220510 (2013).
- Mourik et al. (2012) V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, and L. P. Kouwenhoven, Science 336, 1003 (2012).
- Gangadharaiah et al. (2011a) S. Gangadharaiah, B. Braunecker, P. Simon, and D. Loss, Phys. Rev. Lett. 107, 036801 (2011a).
- Stoudenmire et al. (2011) E. M. Stoudenmire, J. Alicea, O. A. Starykh, and M. P. Fisher, Phys. Rev. B 84, 014503 (2011).
- Zhang et al. (2013) F. Zhang, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 111, 056402 (2013).
- Keselman et al. (2013) A. Keselman, L. Fu, A. Stern, and E. Berg, Phys. Rev. Lett. 111, 116402 (2013).
- Klinovaja and Loss (2013) J. Klinovaja and D. Loss, (2013), arXiv:1312.1998 [cond-mat.mes-hall] .
- Wilson (1977) K. Wilson, in New Phenomena in Subnuclear Physics, edited by A. Zichichi (Plenum, 1977).
- Gangadharaiah et al. (2011b) S. Gangadharaiah, L. Trifunovic, and D. Loss, Phys. Rev. Lett. 108,, 136803 (2011b), 1111.5273 .
- White (1992) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
- Schollwöck (2011) U. Schollwöck, Annals of Physics 326, 96 (2011).
- McCulloch (2008) I. P. McCulloch, (2008), arXiv:0804.2509 [cond-mat.str-el] .
- Kjäll et al. (2013) J. A. Kjäll, M. P. Zaletel, R. S. K. Mong, J. H. Bardarson, and F. Pollmann, Phys. Rev. B 87, 235106 (2013).
- Li and Haldane (2008) H. Li and F. D. M. Haldane, Phys. Rev. Lett. 101, 010504 (2008).
- Pollmann et al. (2010) F. Pollmann, A. M. Turner, E. Berg, and M. Oshikawa, Phys. Rev. B 81, 064439 (2010).
- Turner et al. (2011) A. M. Turner, F. Pollmann, and E. Berg, Phys. Rev. B 83, 075102 (2011).
- Fidkowski and Kitaev (2011) L. Fidkowski and A. Kitaev, Phys. Rev. B 83, 075103 (2011).
- Pollmann and Turner (2012) F. Pollmann and A. M. Turner, Physical Review B 86, 125441 (2012).
- Pirvu et al. (2012) B. Pirvu, G. Vidal, F. Verstraete, and L. Tagliacozzo, Phys. Rev. B 86, 075117 (2012).
- Calabrese and Cardy (2004) P. Calabrese and J. Cardy, Journal of Statistical Mechanics: Theory and Experiment 2004, P06002 (2004).
- Fu and Kane (2009) L. Fu and C. L. Kane, Phys. Rev. Lett. 102, 216403 (2009).
- Akhmerov et al. (2009) A. R. Akhmerov, J. Nilsson, and C. W. J. Beenakker, Phys. Rev. Lett. 102, 216404 (2009).
- Law et al. (2009) K. T. Law, P. A. Lee, and T. K. Ng, Phys. Rev. Lett. 103, 237001 (2009).
- Flensberg (2010) K. Flensberg, Phys. Rev. B 82, 180516 (2010).
- Vidal (2004) G. Vidal, Phys. Rev. Lett. 93, 040502 (2004).
- Vidal (2007) G. Vidal, Phys. Rev. Lett. 98, 070201 (2007).
- White and Feiguin (2004) S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004).
- White and Affleck (2008) S. R. White and I. Affleck, Phys. Rev. B 77, 134437 (2008).
- Mahan (2000) G. D. Mahan, Many particle physics (Springer, 2000).
- Thomale et al. (2013) R. Thomale, S. Rachel, and P. Schmitteckert, Phys. Rev. B 88, 161103 (2013).
- Kivelson and Heim (1982) S. Kivelson and D. E. Heim, Phys. Rev. B 26, 4278 (1982).