Non-topological parafermions in a one-dimensional fermionic model with even multiplet pairing
Leonardo Mazza Fernando Iemini Marcello Dalmonte Christophe Mora
cond-mat.str-el





LPTMS, CNRS, Université Paris-Sud, Université Paris-Saclay, 15 rue Georges Clémenceau, 91405 Orsay, France
Département de Physique, École Normale Supérieure / PSL Research University, CNRS, 24 rue Lhomond, 75005 Paris, France


ICTP, Strada Costiera 11, I-34151 Trieste, Italy
Instituto de Física, Universidade Federal Fluminense, 24210-346 Niterói, Brazil


ICTP, Strada Costiera 11, I-34151 Trieste, Italy


Laboratoire  Pierre  Aigrain,
École  Normale  Supérieure / PSL  Research  University,
CNRS,  Université  Pierre  et  Marie  Curie-Sorbonne  Universités,
Université  Paris  Diderot-Sorbonne  Paris  Cité,  24  rue  Lhomond,  75231  Paris  Cedex  05,  France



We discuss a one-dimensional fermionic model with a generalized ℤ_N even multiplet pairing extending Kitaev ℤ_2 chain. The system shares many features with  models believed to host localized edge parafermions, the most prominent being a similar bosonized Hamiltonian and a ℤ_N symmetry enforcing an N-fold degenerate ground state robust to certain disorder. Interestingly, we show that the system supports a pair of parafermions but they are non-local instead of being boundary operators. As a result, the degeneracy of the ground state is only partly topological and coexists with spontaneous symmetry breaking by a (two-particle) pairing field. Each symmetry-breaking sector is shown to possess a pair of Majorana edge modes encoding the topological twofold degeneracy. Surrounded by two band insulators, the model exhibits for N=4 the dual of an 8 π fractional Josephson effect highlighting the presence of parafermions.


    Christophe Mora
    Received: date / Revised version: date
==========================================



Introduction —
The prospect of localizing Majorana fermions at the two edges of a one-dimensional topological 
superconductor <cit.> has spurred an intense theoretical and experimental activity focused on semiconductor
nanowires <cit.> and chains of magnetic adatoms <cit.> coupled to conventional superconductors.
Not only do the Majorana fermions reveal the topological nature of the bulk groundstate, but their
non-Abelian statistics also holds promises for realizing topologically protected qubit gates
and quantum computation <cit.>.
Parafermions are ℤ_N fractional generalizations of ℤ_2 Majorana fermions <cit.>. Local
parafermionic operators have been proposed in hydrid systems combining fractional quantum Hall states
with superconductivity <cit.>. Their existence, topologically protecting the ground
state (GS) degeneracy, originates from the fractional  edge excitations in quantum Hall states. Improving over Majorana fermions, braiding parafermions
provides access to a richer set of gate operations in quantum computation processing <cit.>.

The existence of local parafermions in truly one-dimensional fermionic systems is still elusive. Symmetry-protected
topological phases of fermions have been classified in one dimension <cit.> according to
the most common symmetries, and none of them  is found to support edge parafermions.
Despite this classification, parafermions were identified using bosonization in some one-dimensional 
fermionic models <cit.>, thus raising two  important questions. The first one is about locality: bosonization
being a non-local transformation, the edge character of these parafermions is not clearly settled. The second
and somewhat related question is whether these parafermions generate topological protection.

In this letter, we delve more deeply into these issues by extending Kitaev's chain <cit.> to a lattice fermionic model with ℤ_N pairing, where N is an even integer. 
Interestingly, the model shares the same bosonization description  with many previous studies where parafermions were identified, and compared to them, it allows for a more microscopic and complete solution through exact results and  numerics.





We find that the GS N-fold degeneracy is a mixture of topology and spontaneous symmetry breaking (SSB). We demonstrate
that it is not generated by localized edge parafermions but rather by N/2 pairs of Majorana fermions acting in each
symmetry-breaking sector. The system also hosts a pair of non-topological parafermions, called poor man's parafermions in Ref. kane2015, that we microscopically demonstrate to be non-local. In fact, it is straightforward to prove that operators, local in terms of fermions, cannot exhibit parafermionic commutation relations (the proof is detailed in [See Supplemental Material for more details on bosonization and the exact solution to the GS, the construction of non-local parafermions and edge Majorana modes and the absence of local parafermions in one-dimension fermionic models]), challenging the existence of strong edge parafermions in truly one-dimensional fermionic systems.
Nevertheless, the N-fold degeneracy of the GS persists with exponential accuracy for ℤ_N-preserving disorder and the model exhibits the dual of the 8 π fractional Josephson effect <cit.>, thus shading a positive light on the possibility of using non-local and non-topological parafermions in future applications.



The model —
We consider a linear lattice of spinless fermions subject to a multiplet proximity effect involving N fermions 
(the case N=2 is simply Kitaev's chain),




    Ĥ_P = - t ∑_j  ĉ_j+1^†ĉ_j  + Δ∑_j  ( ∏_k=0^N-1ĉ_j+k^†)  + h.c.

t is the hopping amplitude, Δ promotes condensation of clusters of N fermions.

 The model breaks the usual U(1) symmetry but possesses a ℤ_N symmetry, ĉ_j →  e^2 i π/Nĉ_j, associated to the number of fermions N̂_e being conserved modulo N. The corresponding charge is the generalized parity operator Q̂ = ω^N̂_e, with  ω = e^2 i π/N, which commutes with the Hamiltonian. 
 
 The low-energy physics is  conveniently described using bosonization <cit.>.  Following standard conventions, we introduce the phase and displacement fields θ̂ and φ̂, where ∂_x θ̂/π is the momentum conjugated to φ̂. The Hamiltonian is written in terms of these bosonic fields Ĥ_P ∼Ĥ_0 + Ĥ_1, where

    Ĥ_0 =  v_F/2 π∫ dx [ K (∂_x θ̂)^2 + 1/K (∂_x φ̂)^2 ]

describes the gapless excitations with velocity v_F close to the left and right Fermi points, and K=1 in the absence of interactions. The pairing term (a_0 is the lattice spacing)

    Ĥ_1 = - Δ/a_0∫ d x  cos (N  θ̂)

tends to pin the field θ, increasing its stifness, and thus the Luttinger parameter K in Eq. (<ref>), which makes the system more attractive with increasing Δ. Specifically, the perturbative renormalization group (RG) equations 

    d Δ/d ℓ = (2 - N^2/4 K) Δ,     d K/d ℓ = ( a_0 Δ/v_F)^2 N^2 ,

with the short length cutoff a_0 e^ℓ, predict a Kosterlitz-Thouless fixed point <cit.> at K = N^2/8, Δ=0, separating a gapless Luttinger phase from a gapped phase. Δ flows to strong coupling in the gapped phase and θ̂ gets locked semiclassically to 2 m π/N, m being an integer. As θ̂ and θ̂+ 2 π are identified, there are N distincts ground states following the ℤ_N symmetry of the model (see Fig. <ref>).





The properties of the gapped phase can be examined more thoroughly with density-matrix renormalization-group (DMRG) calculations for open chains; in Fig. <ref> we present data for a representative choice N=4, Δ = t.
We find N degenerate GSs protected by a gap (top panels); the entanglement spectrum organizes in degenerate multiplets of N levels  (bottom-left panel).
The exponential energy separation is lifted by disorder terms breaking the  ℤ_N symmetry and is immune to others. We tested explicitly the cases of an on-site potential ∑_j δ_j n̂_j and of a local interaction term ∑_j δ'_j n̂_j n̂_j+1  (bottom-right panel) preserving ℤ_N. Despite the partial resilience to disorder and the gap induced by Eq. (<ref>) as  in Refs. <cit.>, the GS degeneracy is not entirely topological and the model is not hosting edge parafermions as discussed below.



Exact solution — A better understanding of the gapped phase can be achieved by seeking an exact solution to the model <cit.>. This is done by modifying the Hamiltonian H_P, i.e. adding the term

    Ĥ_A =    -t  ∑_j [2 n̂_j n̂_j+1 - (n̂_j + n̂_j+1) ]  
       + Δ∑_j  [ ∏_k=0^N-1n̂_j+k + ∏_k=0^N-1(1 - n̂_j+k) ],

such that its ground states can be explicitly determined. One first verifies that the Hamiltonian H_P + H_A is a non-negative operator. We then observe that for an open chain of size L, the following states are zero-energy eigenstates
  
    |ψ_p⟩ = √(N/2^L)∑_j_1 < … <j_N_e N_e = p[N]|{ j_N_e}⟩.

They  correspond to the coherent sum <cit.> of all Fock states |{ j_N_e}⟩ = ĉ_j_1^†ĉ_j_2^†…ĉ_j_N_e^†|0⟩  such that the number of fermions N_e is equal to  p  modulo N. 
Since Ĥ_P + Ĥ_A is non-negative, the N states |ψ_p⟩, with 0 ≤ p ≤ N-1, belong to the GS manifold (we have verified numerically that they exhaust it). As they are eigenvalues of the ℤ_N parity operator,  Q̂|ψ_p⟩ = ω^p |ψ_p⟩, each parity sector has a unique zero-energy GS.



 

Before studying the properties of these states, we need
to ensure that they describe the same gapped phase of Ĥ_P (or, equivalently, of Ĥ_0 + Ĥ_1). 

We use a concrete numerical strategy by building an adiabatic path relating the two models Ĥ_P and Ĥ_P + λĤ_A. 
In Fig. <ref>  we monitor the gap of the model along the path λ∈ [0,1] and test the GS degeneracy and resilience to disorder for λ = 0.5. This demonstrates by continuity the uniqueness of the gapped phase and thus the relevance of the analytical states |ψ_p⟩ to describe it.

Topological properties of the exact model —
Topological phases typically host boundary zero modes in open geometries which induce circular permutations among the degenerate GS. The presence of these modes on the edges reflects the insensitivity of the bulk to perturbations. Closing the chain lifts entirely the GS degeneracy by coupling the edge states <cit.>. This situation contrasts with a non-topological phase with SSB, in which case closing the chain does not change the degeneracy. 

In order to test the topological properties of our model, we consider the Hamiltonian Ĥ_P + Ĥ_A on a  chain with periodic boundary conditions. The states |ψ_p⟩ with odd p remain ground states with zero energy; those with even p frustrate at least one link between lattice sites and thus have energy larger than zero. 
The N-fold GS degeneracy is thus reduced to N/2-fold  when closing the chain (see Fig. <ref>).
This partial reduction unambiguously indicates that the GS degeneracy is only partially topological <cit.>.  

As we show below, the partial lift of degeneracy can be accounted for by the presence of edge Majorana fermions. The residual N/2 degeneracy, instead, comes from SSB  induced by the operator


    Ô_2^(j) = 2 e^-i φ ĉ_j^†ĉ_j+1^† + h.c.,

where φ is an arbitrary phase. 
The peculiar mix of topology and SSB clearly appears once the matrix elements ⟨ψ_p+1|ĉ_j^†|ψ_p⟩, exponentially small in the bulk, and ⟨ψ_p+2|ĉ_j^†ĉ_j+1^†|ψ_p⟩=1/4 are analytically computed.

In parafermionic chains such mixing of topology and SSB is resulting from pair condensation of parafermions <cit.>.

Restricting operators (<ref>) to the GS manifold spanned by the states (<ref>), they are diagonalized by the states

    |n,e⟩ =  ∑_q=0^N/2-1ω^2qn/√(N/2)|ψ_2q⟩,
        |n,o⟩ = ∑_q=0^N/2-1ω^ (2q+1)n/√(N/2)|ψ_2q+1⟩;





σ=e,o denoting the fermionic parity, which spontaneously break the ℤ_N symmetry down to a topological ℤ_2 symmetry. This is shown by the eigenvalues O_2^(j)|n,σ⟩ =  cos( φ - 4 π n/N )  |n,σ⟩ independent on σ.



A further insight into the fact that the pairs |n,e⟩ and |n,o⟩ are protected by topology is given by the observation that
they are also the exact ground states of the Kitaev models (n=1,…,N/2)

    Ĥ_K,n = -t_0  ∑_j  ( ĉ_j^†ĉ_j+1 +  ω^2 n ĉ_j^†ĉ_j+1^† + h.c. )

related to each other by the circular gauge transform ĉ_j →ω ĉ_j. A direct consequence is that each of the following pairs of edge Majorana operators

    α̂_1,n = ω^-nĉ_1 + ω^nĉ_1^†    α̂_2,n = i ( ω^-nĉ_L - ω^nĉ_L^† ),

generates the algebra of Pauli matrices when restricted to the states |n,e⟩ and |n,o⟩.  

To summarize the above discussion, the GS manifold splits into N/2 symmetry-breaking sectors. In each sector, a pair of Majorana fermions enforces a topological degeneracy which is lifted when closing the chain while a residual N/2 degeneracy is preserved.

Topological properties through bosonization —
To confirm the GS structure beyond the fine-tuned Hamiltonian Ĥ_P + Ĥ_A, we return to the bosonized form Ĥ_0+ Ĥ_1 with a chain of finite length L and periodic boundary conditions. We determine the size of the GS manifold through

 
  the mode expansion

    θ̂(x)    =  θ̂_0 + π x/LĴ  + θ̂_k  0 (x) 
    φ̂(x)    =  φ̂_0 + π x/LN̂_e + φ̂_k  0  (x)

that decomposes the bosonic fields into zero modes and finite momentum excitations. The phase variables φ̂_0 and θ̂_0 are conjugate variables to the (integer-valued) current Ĵ and particle number N̂_e [More precisely, the number of fermions is N_0 + N̂_e, where N_0 is a reference value.],

    [φ̂_0, Ĵ] = - i,          [θ̂_0, N̂_e] = - i,

with the gluing (parity) condition <cit.> (-1)^N̂_e =  (-1)^Ĵ. 
Substituting the expansion (<ref>)  into the Hamiltonian Ĥ_0+ Ĥ_1, we get that in the limit of very large pairing Δ the current must vanish Ĵ = 0 for θ̂ to be spatially uniform. 
The gluing condition then enforces that N̂_e is an even integer restricting the GS manifold to N/2 states in agreement with the exact-solution analysis.

Next, we consider the gapped phase in the open geometry with a vanishing particle current on both ends, ∂_x θ̂(x=0,L)=0. In this case, the mode expansion is simplified and Ĵ and φ̂_0 are removed from Eq. (<ref>). The minima of the cosine in Ĥ_1  are the eigenstates |θ_m⟩ of θ̂_0 with eigenvalues 2 m π/N (m=0,1 …,N-1). These states are not eigenstates of the generalized parity operator Q̂,  Q̂|θ_m⟩ = |θ_m+1⟩ from Eq. (<ref>). They are coupled at finite size L by tunneling of instantons between the different minima of the cosine potential in Eq. (<ref>). Such tunnelings transfer a charge Δθ/2 π = 1/N between the two chiral left-moving and right-moving channels <cit.>.
The GS eigenstates of Q̂ are given by the linear combinations <cit.> |Ψ_p⟩ = ∑_m=0^N-1ω^-m p|θ_m⟩, in agreement with (<ref>) and (<ref>) [a tight analogy requires the identification  | n,e  ⟩  + | n,o  ⟩→ | θ_n ⟩ and  | n,e  ⟩  - | n,o  ⟩→ | θ_n+N/2⟩].
Their energy splitting is exponentially small with L as a result of instanton tunneling <cit.>.







Having set the bosonization framework, we reconsider the resilience of the gapped topological phase to disorder made of local ℤ_N-preserving operators, as displayed in Fig. <ref>. It amounts to test whether the different GS can be distinguished by such operators. ℤ_N-preserving operators have a bosonized form involving either the field ϕ̂, derivatives of θ̂ or  e^i N θ. They cannot measure the zero mode θ̂_0 and thus discriminate the different GS <cit.>. By contrast, the bosonized form of the operators (<ref>)

    Ô_2^(j)∼cos (φ- 2 θ̂),   Ô_2^(j)|θ_m⟩ =  cos( φ - 4 π m/N ) |θ_m⟩

clearly discriminate the GS manifold into N/2 sectors in full agreement with the exact solution. The N-fold degeneracy of the model (<ref>) is therefore symmetry-protected <cit.> by number conservation modulo N. Even when pairing terms lift the GS degeneracy, a twofold degeneracy protected by edge Majorana fermions is nonetheless guaranteed.



Parafermions —
At this stage, it might appear surprising to discover that the gapped phase hosts parafermionic zero modes

    γ̂_1 = e^i θ̂_0,     γ̂_2 = e^i ( θ̂_0 + 2 πN̂_e/N).

θ̂_0 is a phase variable with the identification θ̂_0 ∼θ̂_0 + 2 π so that
e^i θ̂_0 is a legitimate operator. θ̂_0 is pinned at low energy in the minima of Eq. (<ref>) and can be written as  θ̂_0 = 2 πn̂_θ/N. n̂_θ is an integer-valued operator which does not commute with N̂_e, as Eq. (<ref>) indicates. As a result, γ̂_1 and γ̂_2 satisfy the properties expected for parafermions, namely γ̂_1^N = γ̂_2^N = 1 and γ̂_1  γ̂_2 = ωγ̂_2 γ̂_1. They commute with the Hamiltonian Ĥ_0 + Ĥ_1 and therefore define zero-energy modes permuting the whole GS manifold with N states. Corresponding parafermionic operators are also constructed from the exact solution <cit.>.   Both constructions highlight the fact that these operators are non-local, i.e. not edge properties of the chain. For this reason, their existence does not contradict the partially topological nature of the gapped phase.
As such, the degeneracy that they describe is not topological: they correspond to the poor man's parafermions discussed in Ref. <cit.>.


These non-local parafermions are beyond a nice mathematical construction: the system exhibits an anomaly very similar to the 8 π Josephson effect predicted for quantum spin Hall edge states <cit.>, as the two models are dual to each other in bosonization. Let us consider the open chain for N=4 connected on each side to a 1D band (Peierls) insulator. The insulators are described in bosonization by the gapping term cos(2 ϕ̂- φ_1/2)  in addition to the quadratic part (<ref>), preserving the ℤ_4 symmetry. φ_1 (φ_2=0) is the phase acquired in the left (right) insulator upon backscattering a right-moving to a left-moving single fermion. 
  As θ_m →θ_m+1 tunneling events in the multiplet pairing region involve fractional fermion backscattering of charge 1/4, their amplitudes are dressed by e^i φ_1/4. Therefore, although the Hamiltonian is invariant upon shifting φ_1 by 2 π, the spectral flow of the splitted GS manifold is 8 π periodic in  φ_1 and level crossings are protected by ℤ_4 symmetry, i.e. fermion number conservation modulo 4. The case N=2, in the context of the magneto-Josephson effet, was solved exactly <cit.>.



Following Ref. <cit.>, it is also possible to introduce local boundary zero-mode operators, i.e. commuting with the Hamiltonian, sitting on the edges of the open chain. They take the form

    β̂_1 ∼ e^i θ̂_0,     β̂_2 ∼ e^i ( θ̂_0 +πN̂_e ).

Although β̂_1^N = β̂_2^N =1,  the fermionic statistics β̂_1 β̂_2 = - β̂_2 β̂_1 implies that these operators permute the pair of states |n,e⟩ and |n,o⟩, see Eq. (<ref>), but do not couple the SSB sectors with different n. In each SSB sector, β̂_1 and β̂_2 are identified as the Majorana fermions α̂_1,n and α̂_2,n (Eq. (<ref>)) introduced with the exact solution. The comparison between the non-local parafermions (<ref>) and the local edge (Majorana) modes (<ref>) reveals that boundary operators, such as β̂_2, can distinguish <cit.> at most the fermionic parity e^i πN̂_e. This contrasts with operators sensitive to the  Q̂-parity e^i  2 πN̂_e /N, such as γ̂_2, but necessarily delocalized on the whole chain.

Conclusions —
 We  discussed a microscopic one-dimensional fermionic model characterized by even multiplet pairing, in many aspects similar to models put forward as candidates for realizing parafermions in one-dimensional fermionic systems.
 We explicitly constructed the parafermions of our model and shown that they are neither local nor topological. The peculiar discussed phenomenology is instead a mixture of topological Majorana physics and SSB. We believe that these results extend beyond our model. However, the eccentricity of the specific SSB that we found might still leave non-topological parafermions as a viable route for technological applications.

  Note added.
 After the upload of this article on the arXiv website, two other independent preprints appeared that, discussing different models and using related techniques, analyzed the phenomenology that is the object of this paper, namely that of topological Majorana physics endowed with spontaneous symmetry breaking <cit.>.
 

Acknowledgements —
We thank A. Bernevig, J. Budich, M. Burrello, R. Fazio, D. Loss, N. Regnault, E. Sela, A. Stern and P. Zoller for enlightening discussions.
L.M was supported by LabEX ENS-ICFP: 
ANR-10-LABX-0010/ANR-10-IDEX-0001-02 PSL*. C. M. acknowledges support from Idex PSL Research University (ANR-10-IDEX-0001-02 PSL).






















  Supplemental Material for: 
 



 Leonardo Mazza, Fernando Iemini, Marcello Dalmonte and Christophe Mora














The different analytical techniques mentionned in the main text are presented with more details in this Supplementary information.
Three topics are covered: the formulation of an exact solution for the ground state at a fine-tuned point, the low-energy bosonisation
approach and the generic demonstration that parafermionic statistics cannot be obtained from spatially separated operators
in one-dimensional fermionic systems.



§ EXACT SOLUTION FOR THE GROUND STATES




 §.§ Number-conserving Hamiltonian


The special Hamiltonian advertised in the main text, H_P + H_A can be decomposed into two contributions,  H = H_P + H_A = H_Hei + H_pairing. With ĉ_j the operator annihilating a spinless fermion at site j and n̂_j = ĉ_j^†ĉ_j the corresponding on-site density, the first term in H,

    H_Hei  = -t  ∑_j=1^L-1[ ( ĉ_j^†ĉ_j+1 +  H.c.) - (n̂_j + n̂_j+1
    ) +  2 n̂_j n̂_j+1],

maps onto the ferromagnetic Heisenberg model upon a Jordan-Wigner transform. The number of fermions N̂_e is conserved with this Hamiltonian. Also by construction H_Hei≥ 0 and a zero-energy ground state is found in each sector with fixed N̂_e = N_e. The ground state wavefunctions take the form
 
    |ψ_N_e⟩ = 1/√(𝒩_N_e)∑_ 1 ≤ j_1 < j_2 < … <j_N_e≤ L|{ j_N_e}⟩,
         𝒩_N_e = LN_e

corresponding to a coherent sum of Fock states 

    |{ j_N_e}⟩ = ĉ_j_1^†ĉ_j_2^†…ĉ_j_N_e^†|0⟩,
 |0⟩ denoting the vacuum state. Already at this level, it is possible to distinguish different boundary conditions. For open boundary conditions, the states (<ref>) are genuine zero-energy ground states of the positive Hamiltonian H_Hei. However, for periodic or antiperiodic boundary conditions ĉ_j+L = ±ĉ_j, only the states (<ref>) with N_e odd/even are ground states of H_Hei.



 §.§ Adding multiplet pairing


We now add the multiplet pairing contribution to the Hamiltonian corresponding to the proximity hopping of N fermions,

    H_pairing = Δ∑_j=1^L-N+1[ ( ∏_k=0^N-1ĉ_j+k^†  + h.c. ) + ∏_k=0^N-1n̂_j+k + ∏_k=0^N-1(1 - n̂_j+k) ]

H_pairing is a  positive operator. Consider A_j = B_j + B_j^2 with B_j = ∏_k=0^N-1ĉ_j+k^†  + h.c.. Noticing that B_j^3 = B_j, one easily obtains  A_j^2 = 2 A_j  which implies that the eigenvalues of A_j are either 0 and 2. H_pairing = Δ∑_j A_j is therefore a sum of positive operators and itself a positive operator. The total Hamiltonian H = H_Hei + H_pairing does not conserve the number of fermions,  [Ĥ,N̂_e]  0, but only the Z_N parity, [Ĥ,Q̂] = 0, where Q̂ = ω^N̂_e and ω = e^2 i π/N.

The Hamiltonian has a complicated form but, most importantly, it is local in terms of the fermionic operators. It has analytical eigenstates at zero energy. Among the states (<ref>), there are N linear combinations which minimize both H_Hei and H_pairing, namely

    |ψ_p⟩ = √(N/2^L)∑_j_1 < … <j_N_e N_e = p[N]|{ j_N_e}⟩,          p=0,…, N-1

with H_Hei|ψ_p⟩ = H_pairing|ψ_p⟩ = 0. These states form the ground state manifold of H. Following the above analysis, for periodic or antiperiodic boundary conditions, only the states with p odd/even are ground states of H.



 §.§ Symmetry breaking operator and degeneracy lifting


Part of the degeneracy is the result of symmetry breaking just as in the quantum Ising model. The corresponding local operator involves the hopping of a pair of fermions (measuring p-wave pairing),

    O_2^(j) = e^i φ ĉ_j^†ĉ_j+1^† + h.c.,

decoupling the odd and even particle number sectors. φ here is an arbitrary phase. Restricting O_2^(j) to the ground state manifold with p odd, it takes the form, for instance when N=6,

    O_2^(j) = 1/4[      0 e^-i φ  e^i φ;  e^i φ      0 e^-i φ; e^-i φ  e^i φ      0 ],

and the same reduced form for p even. Hence, it is easily diagonalized as shown in the main text. Adding a coupling to the local operator O_2^(j) to the Hamiltonian would lift the degeneracy between the different symmetry-breaking sectors. Each of these sectors still preserves a two-fold degeneracy, i.e. a ground state manifold of dimension two, associated to the even/odd parity and to the presence of Majorana edge modes (see the corresponding discussion in the main text). This residual degeneracy is topologically protected.




 §.§ Non-local Parafermions


With the knowledge of the exact wavefunctions, it is possible to define operators which behave as parafermions when restricted to the ground state manifold. To construct these operators, we first introduce the notation

    c̃_j  = ĉ_j e^i π∑_k<jn̂_k,

corresponding to an annihilation operator dressed by the standard Jordan-Wigner string [In the spin language, c̃_j  corresponds to the operator S_j^- flipping the local spin from up to down.]. The string implies that c̃_j commutes with all fermionic operators located on its left,  i.e. [c̃_j,ĉ_k] = [c̃_j,ĉ_k^†] = 0 for k<j.

Using this commutation property, we introduce the set of all operators removing one fermion (modulor N)  from the chain with the ordered products

    K_ℓ⃗,k⃗ = c̃_ℓ_1^†…c̃_ℓ_N_+^†c̃_k_1…c̃_k_N_-    ℓ_1 < … < ℓ_N_+     k_1 < … < k_N_-

and N_+ - N_- = -1 [N]. For any pair of Fock states |{ j_N_e}⟩ and |{ j_N_e'}⟩ with N_e' = N_e-1 [N], there exists a unique operator K_ℓ⃗,k⃗ such that 
|{ j_N_e'}⟩ =  K_ℓ⃗,k⃗|{ j_N_e}⟩. Hence, summing over the whole set of K_ℓ⃗,k⃗ operators, we define a pair of non-local operators 

    γ̂_1 =  N/2^L∑_ℓ⃗,k⃗ K_ℓ⃗,k⃗,         γ̂_2 = γ̂_1 e^2 i π  (N̂_e-1/2)/N,         γ̂_1 γ̂_2 =  ω γ̂_2 γ̂_1.

In terms of the fermionic degrees of freedom, the domain of existence of γ̂_1 and γ̂_2 coincides with the full chain. Acting on the ground states, one finds

    γ̂_1  |ψ_p⟩ = |ψ_p-1⟩        γ̂_2  |ψ_p⟩ = ω^p-1/2|ψ_p-1⟩

This implies γ̂_1^N =  γ̂_2^N=1 within the ground state manifold which, together with the commutation relation (<ref>), shows that γ̂_1 and γ̂_2 behave as parafermions encoding the ground state degeneracy.




§ BOSONIZATION


We discuss below the details of the bosonization method applied to the generalized pairing model.



 §.§ Notations


We assume that interactions and pairing are weak enough and consider only the kinetic energy. We  take the continuum limit of the lattice model and project the fermion field operator around the two Fermi points

    ĉ_j ≃√(a_0)( ψ̂_R (x) e^i k_F x + ψ̂_L (x) e^-i k_F x),

at x = j a_0. We bosonize the two fields 

    ψ̂_R/L (x) = 1/√(2 π a_0) e^i [θ̂(x) ±φ̂(x)],

 a_0 is the short-distance cutoff or lattice spacing (Klein factors are discarded here for simplicity as they play no role in the forthcoming discussion). The two bosonic fields θ̂ and φ̂ satisfy the canonical commutation relation

    [ ∂_x φ̂(x) , θ̂(x') ] = i πδ (x-x').

We now add interactions and pairing, and analyze the model by the perturbative renormalization group method. If we suppose that the fermion interactions are overall attractive, the Luttinger parameter takes a value K >1. Cosine terms involving the field φ̂ flow to zero under RG and are thus irrelevant. Rescaling the coupling Δ to absorb numerical coefficients of order 1, we obtain

    H_bos = v_F/2 π∫ dx [ K (∂_x θ̂)^2 + 1/K (∂_x φ̂)^2 ] - Δ/a_0∫ d x  cos (N  θ̂),

written as H_0+H_1 in the main text. The RG equations are given as Eq. (<ref>) in the main text.



 §.§ Boundary conditions


We first consider periodic boundary conditions and expand the bosonic fields over the mode expansions, written as Eqs (<ref>) in the main text,

    φ̂(x)      = φ̂_0 + π x/LN̂_e + i π/L∑_k  0( L |k|/2 π)^1/21/k e^- a_0 |k|/2 - i k x ( â_k^† + â_-k ), 
    θ̂(x)      = θ̂_0 + π x/LĴ + i π/L∑_k  0( L |k|/2 π)^1/21/|k| e^- a_0 |k|/2 - i k x ( â_k^† - â_-k ),

where k are multiples of 2 π /L, and the phase variables φ̂_0 and θ̂_0 are conjugate variables to the current Ĵ and particle number N̂_e,

    [φ̂_0, Ĵ] = - i,          [θ̂_0, N̂_e] = - i.

In addition, the periodicities of the fields ψ̂_R/L imply the selection rule (-1)^N̂_e = (-1)^Ĵ. When relevant, the pinning of the field θ̂ in the minima of the cosine in Eq. (<ref>) imposes that θ̂ is spatially homogeneous, thus discarding a non-zero value for Ĵ in Eq. (<ref>). From the gluing condition (-1)^N̂_e = (-1)^Ĵ, we obtain that the ground state must have an even N̂_e.

If instead we consider an open chain with a vanishing current at both ends,  ∂_x θ̂(x=0,L)=0, the mode expansion takes a simplier form

    φ̂(x)      = φ̂_0 + π x/LN̂_e + i ∑_k > 0( π/L k)^1/2sin (k x)   ( â_k^† + â_k ), 
    θ̂(x)      = θ̂_0  + i ∑_k > 0( π/L k)^1/2cos( k x)   ( â_k^† - â_k ),

k are multiples of π /L and the canonical commutation relations [â_k,â_k'^†] = δ_k,k', [θ̂_0, N̂_e] = - i apply.
Here all values of the particle number N̂_e are admissible and the ground state manifold has dimension N in contrast to the closed chain.



 §.§ Parafermionic zero modes


Substituting the expansion (<ref>) into the Hamiltonian (<ref>) (or H_0+H_1 in the main text), one finds that the operator θ̂_0 is pinned in the gapped phase to the minima of the cosine. This is taking into account by writting  θ̂_0 = 2 πn̂_θ/N where n̂_θ is an integer-valued operator. Assuming at low energy that the modes â_k have a perturbative effect, we expand the Hamiltonian to isolate the part which involves the zero modes

    H_bos = π v_F/2 K LN̂_e^2 - Δ/a_0 L cos ( N θ̂_0 ) + …

where the cosine can also be written cos ( N θ̂_0 ) = cos(2 πn̂_θ). The Z_N parity operator Q̂ = ω^N̂_e shifts the value of n̂_θ by one units, Q̂n̂_θQ̂^-1 = n̂_θ -1, and therefore it commutes with the Hamiltonian [H_bos,Q̂]=0.

For a system of sufficient large length L ≫ 1, the second term in Eq. (<ref>) dominates and e^i θ̂_0 = e^2 i πn̂_θ/N commutes with the Hamiltonian  H_bos (it also commutes at finite length when K →∞ which is the case for the exact solution of Sec. <ref>). The pair of operators 
γ̂_1 and γ̂_2 given by (also Eq. (<ref>) in the main text),

    γ̂_1 = e^i θ̂_0,         γ̂_2 = γ̂_1 e^2 i π  (N̂_e-1/2)/N,         Q̂ = ω^1/2γ̂_1^†γ̂_2,

commutes with the Hamiltonian, [γ̂_1/2,H_bos]=0, and satisfies commutation relations (see main text) characterizing them as parafermions.

The expressions for the parafermions can also be obtained from the analysis of the exact solution in Sec. <ref>. We bosonize the operator

    c̃_j  =  1/√(2 π) e^i θ̂(j a_0) ,

which only involves the field θ̂(x) due to the attached Jordan-Wigner string. We then make use of the expansion (<ref>) in Eq. (<ref>) and neglect the modes â_k by assuming proximity to the ground state. The result coincides with Eq. (<ref>) which strongly points towards to the fact that these parafermionic modes are non-local.




 §.§ 8 π-periodic spectral flow dual to the fractional Josephson effect


Electron-electron interactions in a quantum spin Hall edge state have been shown <cit.> to mediate a fractional Josephson effect with 8 π periodicity protected by time-reversal symmetry and fermion parity. As we now discuss, our model of one-dimensional fermions with multiplet pairing is dual to model of Ref. zhang2014 and therefore predicts a similar ground state anomaly.

We focus on the case  N=4 for simplicity.
At low energy (short-distance variations ∼ 1/k_F are discarded), the model exhibits an emergent antiunitary symmetry Ξ

    i →  -i     ψ̂_R (x) →ψ̂_L^† (x)     ψ̂_L (x) →  - ψ̂_R^† (x)     Ξ^2= -1

commuting with the Hamiltonian [Ξ,H]=0. In terms of bosonic fields, the symmetry acts as θ̂→θ̂+ π/2, ϕ̂→  - ϕ̂- π/2. Exchanging the transformations on θ̂ and φ̂ precisely recovers time-reversal symmetry, namely ψ_R →ψ_L, ψ_L → -ψ_R.   Whereas a local Zeeman term ∼cos 2 ϕ̂ breaks time-reversal symmetry in Ref. zhang2014, here the symmetry Ξ is broken by a pairing term ĉ_j  ĉ_j+1∼ψ_R (x) ψ_L (x) ∼cos 2 θ̂ distinguishing the states  | θ_m ⟩.

To reveal the ground state anomaly dual to the 8 π Josephson effect, we sandwich the multiplet pairing model (Eq. (<ref>) in the main text) between two band insulators with commensurate potentials cos (2 k_F x + φ_1/2). The phase φ_1 (φ_2) parametrizes the positions of the potential minima in the insulator located on the left (right). We set φ_2=0 for simplicity. The corresponding bosonized form of this backscattering term is ψ_R^†ψ_L e^i φ_1/2 + h.c. ∼cos(2 ϕ̂- φ_1/2). The commensurate potentials induce backscattering of single electrons transferred between the right-moving and left-moving channels with the phases φ_1/2.

The setup thus alternates two regions dominated by cos(2 ϕ̂) with a central one gapped by cos(4 θ̂) and is thus dual to the model of Refs. <cit.>. It is known <cit.> that a magnetic domain wall in a quantum spin Hall edge state binds half an electronic charge e/2 and that the same charge e/2 is also pumped when the magnetization is reversed. This is due to the pinning of the field ϕ̂ by the cos(2 ϕ̂) term which shifts its minima under these processes. In our setup, a tunneling event between the states θ_m and  θ_m+1 involves a change in the conjugated field Δθ =  π /2 corresponding to the transfer of a fractional charge  Δθ/2 π = 1/4 between the right-moving and left-moving channels. Such instanton tunneling is accompanied by the phase term e^i φ_1/4 and the Hamiltonian governing the low-energy subspace is given by (m=4 and m=0 are identified) <cit.>

    H = Γ∑_m=0^3 e^i φ_1/4 | θ_m ⟩⟨θ_m+1 | + h.c.

where the energy scale Γ for instanton tunneling is exponentially small with the size L. The resulting spectral flow is 8 π periodic in φ_1. The energy crossings are protected by ℤ_4 symmetry, i.e. each eigenstate has a well-defined number of fermion modulo 4. Similarly to its dual model <cit.>, the crossings at φ_1=0   [2 π] are also protected by the symmetry Ξ while the crossings at φ_1= π  [2 π] are protected by standard fermion number parity. The former are splitted in the presence of a local pairing term ∼cos 2 θ̂ breaking  Ξ whereas the latter are topologically protected by the edge Majorana zero modes.




 §.§ Majorana edge modes


In fractional quantum Hall states with proximity s-wave superconductivity, parafermionic operators have been explicitely constructed by Clarke, Alicea and Shtengel <cit.> and rigorously shown to be localized at edges. We adapt below their analysis to our model and show that edge Majorana fermions rather than edge parafermions  thus emerge. We first add two regions of small size ℓ≪ L where Δ =0. They are localed at the left - ℓ <x<0 and right  L <x<L+ℓ sides of the system. In the left region, the bosonic fields φ̂_R/L (x) = φ̂(x) ±θ (x) are given by

    φ̂_R/L (x)   =  ±2 πn̂_θ/N±√(2)∑_q=0^∞( e^± i λ_q (x)b̂_q + h.c. )/√(2 q + 1)        φ̂_L (- ℓ) =  - φ̂_R ( - ℓ) 
    φ̂_L (0) =  φ̂_R (0) - 4 πn̂_θ/N

with λ_q (x) = (2 q+1) π (x + ℓ)/2 ℓ and the bosonic modes b̂_q satisfying [b̂_q,b̂_q'^†] = δ_q,q'. Using the Hamiltonian (<ref>) (with Δ=0), the following identity can be derived (we consider K=1 for simplicity)

    [ H_bos, e^i  φ̂_R/L (x)  ] = ± i v_F ∂_x e^i  φ̂_R/L (x) ,

such that the operator

    β̂_1 = e^2 i πn̂_θ/N∫_- ℓ^0 d x [  e^-2 i πn̂_θ/N e^i φ̂_R (x) +  e^2 i πn̂_θ/N e^i φ̂_L (x) + h.c. ]

is a zero-mode commuting with the Hamiltonian H_bos. Noting that  φ̂_L (0) =  φ̂_R (0) - 4 πn̂_θ/N, we observe that the expression of  β̂_1 involves only the local fermionic fields e^± i  φ̂_R/L (x,0) localized inside the small region left to the system. At low energy and for ℓ→ 0, the modes b̂_q are essentially gapped <cit.> and the boundary mode simplifies as β̂_1 ∼ e^2 i πn̂_θ/N as given by Eq. (<ref>) in the main text. A similar construction gives the right boundary operator β̂_2 also given Eq. (<ref>).




§ ABSENCE OF EDGE PARAFERMIONS IN ONE DIMENSIONAL FERMIONIC MODELS


In this section we present details for the demonstration 
that parafermionic statistics cannot be obtained from spatially separated operators
in one-dimensional fermionic systems.

 Let us first introduce 2L Majorana fermions α̂_j for the L sites of model,
 satisfying the usual anticommutation relations:
 
    {α̂_j,α̂_ℓ} = 2 δ_j,ℓ,    α̂_j^†  = α̂_j,

with j=1,..,2L. A complete basis for the full Hilbert space ℋ, of dimension 2^L, can be constructed 
 in terms of products of such Majorana fermions, as follows <cit.>,
 
    Γ̂_j⃗≡α̂_j_1α̂_j_2 ... α̂_j_m,

 where j⃗ = (j_1, j_2, ..., j_m) is an ordered vector of size m
 parametrizing a list of occupied modes, with j_i = 1,...,L and m=0,...,L. 
 For the special case m=0 ones defines Γ̂_j⃗≡𝕀.
The operators Γ̂_j⃗ are linearly independent
and any operator Ô in the Hilbert space has the decomposition,
 
    Ô = ∑_{j⃗| j_i ∈ S(O) } f_j⃗ Γ̂_j⃗

 where S(O) = (s_1,s_2,...,s_L_O)  denotes the support of the operator along the sites of model,
  with s_i=1,..,L representing the sites the operator occupies, and L_O = 1,...,L its effective length.
 
 Consider now two operators O_A and O_B, with local supports S(A) and S(B) in exclusive regions
  of space: S(A) ∩ S(B) = ∅. Using the decomposition of Eq.(<ref>), 
  the commutation relation Ô_A Ô_B + εÔ_B Ô_A, with ε an a priori arbitrary phase, 
  can be written as,
 
    Ô_A Ô_B + εÔ_B Ô_A    =   ∑_{j⃗_A | j_A_i∈ S(A)}{j⃗_B | j_B_i∈ S(B)}
     f_j⃗_A f_j⃗_B  (Γ̂_j⃗_AΓ̂_j⃗_B  + εΓ̂_j⃗_BΓ̂_j⃗_A ) 
       =   ∑_{j⃗_A | j_A_i∈ S(A)}{j⃗_B | j_B_i∈ S(B)}
     f_j⃗_Af_j⃗_B Γ̂_j⃗_AΓ̂_j⃗_B  (1 +εϕ_j⃗_A,j⃗_B )
       =   ∑_{j⃗_A | j_A_i∈ S(A)}{j⃗_B | j_B_i∈ S(B)}
     f_j⃗_Af_j⃗_B Γ̂_j⃗_AB≡j⃗_A∪j⃗_B   (1 +εϕ_j⃗_A,j⃗_B )

where in the second line we used the commutation relation Γ̂_j⃗_BΓ̂_j⃗_A  = 
ϕ_j⃗_A,j⃗_BΓ̂_j⃗_AΓ̂_j⃗_B, with ϕ_j⃗_A,j⃗_B a phase 
which can assume values ± 1, depending only on the length of the vectors j⃗_A and j⃗_B.
 In the third line we used the fact that 
 Γ̂_j⃗_AB≡j⃗_A∪j⃗_B =  Γ̂_j⃗_AΓ̂_j⃗_B,
 since the support of the operators Γ̂_j⃗_A and Γ̂_j⃗_B are in exclusive regions 
 of space with their respective vectors j⃗_A and j⃗_B not overlapping between each other.
 
 Notice that the operators Γ̂_j⃗_AB in the sum of Eq.(<ref>) are orthogonal to each other.
 Thus, in order to Ô_A and Ô_B satisfy the commutation relation Ô_A Ô_B + εÔ_B Ô_A =0 
 all of the terms in the sum 
 must vanish,
 
    (1 +εϕ_j⃗_A,j⃗_B ) = 0,    ∀j⃗_A,j⃗_B,

 a condition that can only be satisfied if ε = ± 1, corresponding in this way 
 to fermionic or bosonic commutation relations.





