{'fraction_non_alphanumeric': 0.03921330581522399, 'fraction_numerical': 0.01015741977257092, 'mean_word_length': 4.454475223485266, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Majorana Zero Modes in Fermionic Wires coupled by Aharonov-Bohm Cages

We devise a number-conserving scheme for the realization of Majorana Zero Modes in an interacting fermionic ladder coupled by Aharonov-Bohm cages. The latter provide an efficient mechanism to cancel single-particle hopping by destructive interference. The crucial parity symmetry in each wire is thus encoded in the geometry of the setup, in particular, its translation invariance. A generic nearest-neighbor interaction generates the desired correlated hopping of pairs. We exhibit the presence of an extended topological region in parameter space, first in a simplified effective model via bosonization techniques, and subsequently in a larger parameter regime with matrix-product-states numerical simulations. We demonstrate the adiabatic connection to previous models, including exactly-solvable ones, and we briefly comment on possible experimental realizations in synthetic quantum platforms, like cold atomic samples.

We devise a number-conserving scheme for the realization of Majorana Zero Modes in an interacting fermionic ladder coupled by Aharonov-Bohm cages. The latter provide an efficient mechanism to cancel single-particle hopping by destructive interference. The crucial parity symmetry in each wire is thus encoded in the geometry of the setup, in particular, its translation invariance. A generic nearest-neighbor interaction generates the desired correlated hopping of pairs. We exhibit the presence of an extended topological region in parameter space, first in a simplified effective model via bosonization techniques, and subsequently in a larger parameter regime with matrix-product-states numerical simulations. We demonstrate the adiabatic connection to previous models, including exactly-solvable ones, and we briefly comment on possible experimental realizations in synthetic quantum platforms, like cold atomic samples.

I. MOTIVATION
In the last decade, the quest for topological states of matter has arguably been one of the driving directions in condensed matter physics [1][2][3], partially motivated also by their envisioned usage as platforms for quantum computation [4]. Among all possible topological states, Majorana Zero Modes (MZM) are one of the simplest examples realizing an anyonic excitation, i.e., not obeying either fermionic or bosonic statistics [5][6][7]. Despite them being conceptually quite simple, an ongoing quest is being pursued towards an unambiguous measurement of their existence. While they do not enable a universal set of gates, they have been proposed to serve as a topological protected quantum memory [8,9].
Stemming from the cornerstone paper by Kitaev [10], most efforts have focused on a one-dimensional realisation of Majorana Zero Modes via coupling some semiconducting nanowire to a bulk superconductor [11][12][13][14][15][16][17]. The superconductor serves as a reservoir inducing p-wave superconductivity into the nanowire via the proximity effect resulting in an effective breaking of the U(1) symmetry of number conserving down to a residual fermionic parity symmetry Z 2 .
In recent years, alternative proposals for realizing MZMs without breaking the number conservation were put forward [18,19]. These schemes are based on the field-theoretical observation that a minimal U(1) × Z 2 model of two (fermionic) Luttinger liquids coupled exclusively by a pair-hopping term indeed leads to the same topological signatures [20]. Noticeably, even some exactly solvable instances were found [21,22], giving deep insights into the nature of this phase. Moreover, a number-conserving scheme is particularly appealing for synthetic quantum matter platforms like cold atoms [23]. However, a perfect cancellation of single-particle tunnelings between the chains is needed to ensure the Z 2 protection of MZMs. In previous works, this was only achievable in an approximate fashion via some perturbative suppression. Here, we present an alternative scheme which makes use of exact interference terms of complex hopping amplitudes, also known as Aharonov-Bohm Cages [24,25], and perfectly cancels all single-particle poisoning. These cages are arranged in a translation invariant sequence across the two target fermionic chains, and a generic nearest-neighbor interaction term enables the sought-after correlated hopping of particles.
The paper is structured as follows: First in section II we introduce our model which involves four spinless fermionic species. After discussing the basic properties of that model, we integrate out two of these spinless fermions using a Schrieffer-Wolff (SW) transformation [26], and show that the obtained effective Hamiltonian falls in the same class as those of previous proposals. This effective Hamiltonian is first investigated in section III by using bosonization [27,28], in order to find the most favorable parameter regime for realizing the MZM phase. In this section we also review the basic indicators used for detecting the MZM phase: The non-local behavior of the end-to-end correlation function together with a relative sign between the ground-states of the two parity sectors and the exact double degeneracy of the entanglement spectrum [29]. Next, in section IV, we show numerical results using tensor network techniques [30], exhibiting all defining features of the MZM phase, not only for the effective model, but also for the full four-flavor setup in regimes very far from the perturbative expansion conducted before. Finally, in section V we summarize our findings and give a short outlook of open questions.

II. MODEL INTRODUCTION
Let us consider two (lattice) wires a and b, populated by spinless fermions, and connected to each other via additional sites c and d in a rhomboidal configuration pierced by a π-phase, as described by the following Hamiltonian: where α ( †) j annihilates (creates) a fermion in the site of kind α = a, b, c, d at the j-th lattice position, and p = superposition of the intermediate modes. For the sake of simplicity, we picked up a gauge where the whole πphase has been collected on a single link: however, any redistribution along the rhombi would, of course, lead to the same Aharonov-Bohm caging effect in the end. As made explicit by Eq. (2), indeed, destructive interference prohibits single-particle motion between the (decorated) wires, i.e., the two charges with N α = j n α j = j α † j α j , are separately conserved and denoted as the U (1) ± symmetries in the following.
In order to let Majorana physics emerge, we need to partially break the U (1) − associated to the N − charge into a residual Z 2 symmetry for each dressed chain. We could naturally achieve it by considering nearestneighbour density-density interaction terms of the kind: which does not preserve N p and N m (but still preserves their parity) away from the fine-tuned point The residual unbroken symmetry is actually (Z 4 ) − /(Z 2 ) + , since the parity of the conserved overall population fixes the parity of the relative population, too: as a convention, we decide to look at the parity in the dressed a chain, i.e., Actually, we show in App. C that our model is adiabatically connected to a regime where the U (1) + symmetry is further broken and the residual group is an even simpler Z 2 × Z 2 , i.e., the same symmetry class as two individual Majorana chains. Later we will abuse this relation to derive the signatures of the topological regime in the single particle correlation functions. Any additional generic intra-wire Hamiltonian H α=a,b and any Hamiltonian of the kind acting on the intermediate sites would still fall in the same U (1) + ×(Z 2 ) − symmetry class. We initially set t = t ⊥ = 0 for the sake keeping most calculations analytically feasible, but in App. A we provide some estimate on their utility for making the desired topological signatures even more evident. A pictorial sketch of the generic Hamiltonian, is given in Fig. 1b. Once we integrate out the intermediate sites (c, d) via a Schrieffer-Wolff Transformation along the lines of Ref. [26], we are left with a low-energy description of the dressed wires (a, b), illustrated in Fig. 1c: The form of Eq. (8) allows for a direct comparison with the model of Ref. [18] and the exactly solvable one of Ref. [21], as discussed in App. C. The two pair-hopping terms have rather different effects: The inter-chain one, W 1 , embodies the original Kitaev-chain model per each wire separately, and it is indeed the one responsible for the desired topological effect [18,20,21]. The intrachain one, W 2 and W 3 , instead promotes a (pseudo-)spindensity wave ordering in the wire-label degree of freedom, as we will discuss after considering the bosonized version of the Hamlitonian (8) below Eq. (12). The attainable  7), with intra-chain hopping elements (black), interchain ABC hoppings (blue) as in panel (a), and bubbles indicating the density-density interactions; (c) Effective low energy model of Eq. (8), with the blue arrows standing for the correlated inter-chain pair hopping, and the red/orange ones for the cross-correlated hopping terms. From the microscopic derivation one finds W1 = W2 and W3 = 0, however we considered all three couplings for comparison to the exact solvable model of [21], see App. C. couplings with the presented microscopic derivation are: and W 3 = 0. Due to the similar effect of W 2 and W 3 , we set W 3 = 0 for the rest of the main text. We will show in the following, via a combination of field-theory calculations and numerical simulations, that W 2 = W 1 does not impair the formation of the wished topological order, at least in an extended region of the parameter space. The dressed HamiltoniansH α (α = a, b) read where we used the abbreviation for the commutator , it amounts to a simple renormalization of the bare parameter in H α , along with some three-body interactions, density-supported and next-nearest-neighbour hoppings. From renormalization point of view these terms can be safely considered to be less relevant. Therefore we decide to drop them henceforth to keep the model simpler, and leave details for the interested reader in App. B.

III. BOSONIZATION
We now proceed with a field-theoretical analysis of the toy Hamiltonian (8) via bosonization along the notation of Ref. [27]. After having introduced density and phase fluctuating fields, ϕ α and ϑ α , for both fermionic species, α = a, b, and their (anti-)symmetric combinations ϕ ± = 1/ √ 2(ϕ a ± ϕ b ) (same for ϑ ± ), we find the following Hamiltonian: where K τ and v τ denote the Luttinger parameter and Fermi velocity in the τ = ± sectors, and we already dropped less relevant terms (see App. B), including those becoming resonant only at half-filling. We notice that a similar bosonized Hamiltonian appears when considering spinful fermions with anisotropic spin interactions [31], and moreover with β 2 = 0 in other discussions of numberpreserving models for Majorana zero modes [18,20].
Besides an ordinary Tomonaga-Luttinger liquid in the symmetric (charge) sector, τ = +, which is therefore always gapless, the anti-symmetric (spin) sector, τ = −, exhibits a a double Sine-Gordon interaction with bare couplings and scaling dimensions: It is therefore apparent that a gap will arise whenever K − = 1 [32]: while for K − < 1 the ϕ − field is pinned and the phase is trivially a charge-density-wave or singletpairing, depending on the sign of the β 2 coupling [33], the β 1 term causes instead the appearance of unpaired Majorana edge modes for K − > 1 [20]. Furthermore, the refermionization argument for K − = 2 given in Ref. [20] also shows that this mode has to be identified with single particle transitions from one of the two chains to the other.
Therefore, one expect certain overlaps like Ω oo | a † j b j |Ω ee to be nonzero at the edges with a exponential decay to a possible non-zero bulk value. Otherwise said, the degenerate ground-states in the topological phase are related to the ground-states of two independent Kitaev-Majorana chains, once projected on a fixed total charge N + [21]. Thus, using the standard classification of topological insulators, one expect two pairs of MZM to be present in the system, for the unconstrained model, corresponding to the two independent Kitaev-Majorana chains. However, after fixing the total particle number we do not observe 4 = 2 × 2 Majorana modes, but -depending on the total parity P + -only two out of the four possible combinations.
The picture of two unconnected chains is particularly useful to get grasp of some fundamental behaviour of single-particle correlation functions, which we will employ as fingerprint of the desired edge physics. Let us consider the (four) ground-states to be connected to with |Ω being the vacuum of the theory, and c α,E the fermionic operator formed by two Majorana edge modes, with . . . denoting the (gapped) excitations of the system, the correlation length, and A being a normalization factor. For each single chain then holds where P = ± for the even/odd sector, andG(j) is the exponentially decaying correlation function coming from the residual (gapped) excitations of the spectrum. In section IV we will use this exponential revival of the end-toend correlation function together with this characteristic relative sign between the two parity sectors as one of the indicators for having a MZM phase. Closely related to that behavior of the correlation function is the vanishing of the energy gap between the two parity sectors: ∆E = |E − − E + | ∼ e −2L/l . A second indicator is provided by studying the entanglement spectrum [34], which should be exactly double degenerate in the case of being in a Majorana-like phase [29]. While working at fixed particle number might circumvent the formation of a charge gap by forbidding hybridisation of different fillings (as it is indeed the case in our setup), the spin sector remains instead gapped (see App. B).Therefore, we expect an exponentially decaying behavior to the middle of the system, followed by an exponential revival with a π phase difference between the two ground-states. The same holds true for matrix elements of inter-chain operators like the so-called Majorana wave-function [21], Ω ee | a j b † j |Ω oo . Since the characteristics is similar for both observables, we decide to only present results for the single particle correlation functions. This is also motivated from the fact that, in a generic interacting model, the overlap Ω ee | a j b † j |Ω oo may have a non-zero bulk value, making it harder to uniquely identify the edge contribution. This problem is absence for the single particle correlation function, since a j is fundamentally zero.
In order to determine the most favourable regime of the microscopic parameters for achieving the topological phase, we consider the perturbative RG equations (strictly valid only around K − ≈ 1) [33]: We now have to integrate these differential equations starting from the bare values of K − and β j on the original lattice couplings, Eq. B17 (see App. B for details). Thereby we get a rough estimation of the phase diagram, presented in Fig. 2: The exact position of the phase boundaries is (highly) depending on the non-universal constant A. Interestingly, anyway, both the very asymmetric role played by W 2 and the strikingly almost straight critical lines in the (W 1 , W 2 , U α = 0)-plane can be predicted by the equation with two non-universal constants D and C. This equation represents the linearized version of the criticality condition found in [31], see also appendix B 1 for more details. By inspecting Fig. 2, we notice that, in the absence of intra-wire interactions (U α = 0), the line W 1 = W 2 = W dictated by Eq. (9) is well inside the topological phase for W < 0, while no definite conclusion can be reached on the boundary for W > 0. Noticeably, for U α < 0 the bare parameters for the RG-flow are sensibly pushed away from the boundary, thus making the topological phase observable also for W > 0, see Fig. 2b.
We stress here that the extra terms in H c,d of Eq. (6) are foreseen to contribute to stabilize the topological phase, too. Including non-vanishing couplings t and t ⊥ , indeed, a regime with W 1 > W 2 can be achieved, which pushes the model deeper into the topological region even for U α = 0, as can be seen from Fig. 2a (see App. A for more details).

IV. NUMERICAL RESULTS
Next, we validate the cut-off dependent predictions of bosonization against unbiased numerical simulations on the lattice, performed via Matrix Product States (MPS) not only for the effective Hamiltonian (8), but also for the full model (7), i.e., without introducing any perturbative description. We focus on two pristine indicators of the emergence of unpaired Majorana modes at the edges, i.e., i) finite end-to-end single-particle correlations with an exponential decay in the bulk, with relative π-phase between the two parity sectors, and ii) double-degeneracy of the entanglement spectrum, dictated by the Z 2 protecting symmetry.
We conduct our numerical investigations at a fixed density of ν = N + /(2L) = 1/3, so that additional resonances arising at half-filling are avoided. As an exemplary parameter set for the effective model we choose with a chain of length L = 256 and N + = 170 fermions in the system. The specific choice of U nn = 0 was made to simplify the number of parameters to a minimum, without affecting the qualitative picture, as we verified for a a wide range of U nn . Indeed, from a RG point of view, the operator coupled to U nn only has a minimal influence by slightly detuning the bare Luttinger parameter K − and decreasing the bare coupling strength β 2 of the bosonized Hamiltonian, see App. B. Moreover, this choice is always reachable, at least in this fourth-order effective Hamiltonian description, by suitably tuning the bare interaction parameters V 1 and V 2 relative to µ and J.  : Single particle correlation function a1a † k between the leftmost site of the chain and the k-th one. The magenta and blue lines are representing the expectation value to the ground-state in the parity sector P = ±. The black line displays a simulation with neglecting the parity conservation allowwing a superposition between the two different sectors. The inset zooms around the right edge of the chain, revealing the relative π phase between the recovery of the correlation function in the two parity sectors.
First, Fig. 3 illustrates the decay of the single-particle terms a 1 a † x ± with correlation length ≈ 7.35 and their strong revival at the opposite edge r ≈ ∓0.5, with the sign depending on the parity sector, as discussed in Sec. III and predicted in Eq. (15). The quantity r is thereby defined as the amplitude of an exponential fit performed on both ends of the correlation function. Ignoring the underlying Z 2 symmetry in the DMRG simulation results into a pure exponential decay, indicating an equal weighted superposition of the two (nearly) degenerated ground-states of the two parity sectors, see black line in Fig. 3. Such scenario is confirmed by examining the parity expectation value in this setup, P ≈ 0. This is an expected behavior, since DMRG favors the least entangled ground-state [35][36][37]. Second, Fig. 4a shows the entanglement spectrum [34], i.e., − ln λ 2 j with λ j the Schmidt values of a L/2 bipartition of the system, as a function of the quantum numbers δN + = n + − N + /2 and P − = (−1) na . On one hand, the perfect double-degeneracy between the two parity sectors is a clear fingerprint of the symmetry protected topological nature of the anti-symmetric channel (τ = −) [29].
On the other hand, the parabolic shape (with particlehole symmetry) indicates the gapless nature of the symmetric channel (τ = +) [38,39], with the curvature giving back a Luttinger parameter K + 0.97 [40], pretty close to the bare value of Eq. (B17), K  The degeneracy split closes exponentially in system size, with a decay length roughly equal to twice the singleparticle correlation length l since this splitting originates in the exponential small overlap between the two Majorana wave-functions localized at the two ends of the chain [10].
On the other hand, the energy gap to the first excited state vanishes as 1/L originating from the discretization of the momentum in a finite size system as expected for a wellbehaved Luttinger Liquid [33] having a linear dispersion relation.
Noticeably, from Fig. 5 it can be seen that, even in the presence of this gapless channel, the energy difference between the even and the odd sector vanishes exponentially, as one would expect for a system with two topological ground-states. The finite-size gap inside each parity sector, instead, vanishes algebraically with L −1 as expected for a standard Luttinger liquid with a linear dispersion relation + (q) = v + |q|. Additionally, in App. C we show that an adiabatic path exists between our effective model and the exactly solvable one of Ref. [21]. In this article, a path is called adiabatic if the U(1) × Z 2 symmetry group is preserved all along that path and the single particle gap, as defined by the antisymmetric sector, stays finite. This is analogous to requiring that we observe a finite correlation length of the single particle correlation function smaller than the system size along that path. Constructing such a path requires the insertion of some extra operators, but its cartoon projection in the W 1 − W 2 − U nn parameter-space is illustrated as a dashed line in Fig. 8. The energy gap, extracted via fitting the exponential decay of single-particle correlation in the bulk, is plotted in Fig. 9a, and is evidently nonvanishing.
In the same Appendix C, we also demonstrate numeri-cally the adiabatic connection to the two uncoupled Majorana chains by further breaking the residual conservation of the total particle number down to a total parity and only requiring a residual global symmetry group of After reporting the topological fingerprints found in the effective model, we also studied the full model of Eq. (7) and showed that the topological phase persists once the full four-flavor Hilbert-space is taken into account. For this, we studied at the same indicators as reported above, namely the non-local correlation functions with the characteristic relative sign between the two parity sectors as predicted by Eq. (15) and the double degeneracy of the low lying entanglement spectrum. As a warm-up we treated the model with a set of parameters deep inside the perturbative regime of the effective Hamiltonian. The results are reported in the appendix E, where we indeed find the full model to have non-local correlation functions as expected. Nonetheless, the parameters in this regime are not so appropriate to experimental realizations with interaction strengths of several orders of magnitude in difference.
However, we can do better by going away from the perturbative regime. Indeed, we showed that the Majoranalike phase is not bounded to the perturbative regime, but is considerably extended to a more realistic parameters, where we find a double degenerated entanglement spectrum and non-local correlation functions for a large range of parameters.
As an example consider Fig. 6a. In this figure we computed the average degeneracy of neighboring Schmidt values of the low lying entanglement spectrum for a cut at half of the system of length L = 60: computed by fixing the following set of parameters: Noticeably, there seems to be an overall separatrix region of nearly vanishingλ between two regions having a finite splitting between the neighboring entanglement values. Along that set of parameters we expect the Majoranalike phase to appear: This is confirmed by the end-to-end correlation function showing the characteristic π-phase revival, see Fig. 6b as one example. Furthermore, we stress here that this result is not sensitive to the concrete choice of the chemical potential and the filling. We have explicitly checkt it for all µ/J ∈ [2,5] and also for n e = 40, which corresponds to the 1/3 filling in the effective model.
As a final remark we want to discuss possible experimental platforms. We emphasize that the crucial ingredient in realizing our proposal is the cylindrical like structure defined by the rhombi-Hamiltonian H ♦ , Eq. (1). The requirement of having periodic boundary conditions along one direction together with imprinting a effective phase is usually a hard task in physical set-ups. Recently this task was achieved by using the internal degrees of freedom of cold atoms as synthetic dimensions and imprinting arbitrary gauge fluxes to the atoms [41,42]. Together with the good controllability of hopping transisiont by loading the cloud of atoms to an optical lattice [43][44][45][46] and the reliability of species with sizable nearest-neighbor interactions such as polar atoms [47][48][49] suggest cold atoms as the perfect platform, but also other synthetic platforms could be valid However, a concrete realization goes beyond the scope of this article.

V. CONCLUSION & OUTLOOK
Motivated by the ongoing search for an unambiguous detection of topological Majorana zero modes, we have put forward a new number conserving realization of a Majorana-like phase. Our proposal uses the geometry of the underlying lattice together with the Aharonov-Bohm effect to achieve exact cancellation of all possible single particle processes and only allow for pair transitions.
Starting from a perturbative analysis, we found clear fingerprints of a Majorana-like phase with Majorana zero modes being present in an extended parameter regime. We also showed that this is still true in the full model far away from the perturbative regime, thus rising the hope for realizations using synthetic dimensions in a cold atom platform. Due to the large amount of possible parameters to tune, we postpone the development of a concrete scheme together with an experimentally reachable parameter space to future work.
Among the important open questions for all possible quasi-one dimensional number conserving set-ups, the influence of finite temperature on the Majorana Zero Modes plays an important role. To be concrete, it is unclear how possible higher order terms in the bosonization may couple the symmetric gapless sector to the antisymemtric gapped sector hosting the Majorana Zero Modes. Such effects could lead to a much smaller lifetime of the MZM than expected from the limit of two unconnected Kitaev-Majorana chains. This will be the subject for future investigations.

VI. ACKNOWLEDGMENT
We thank M. Burello, R. Egger, A. Haller and M. Köhl for inspiring discussions about theoretical and experimental aspects of this work. We also want to warmly acknowledge a fruitful exchange with F. Lisandrini and C. Kollath, working on a different scheme for the same scope [50]. This work has been funded by the Deutsche

Appendix A: Derivation of the Effective Hamiltonian
In this Appendix, we briefly recall the basics of the Schrieffer-Wolff (SW) transformation in the version of Ref. [26], which we used to derive the effective model of Eqs. (8)- (9). Furthermore, we discuss the effect of the extra terms (t , t ⊥ ) of Eq. (6), showing that they may lead to an even more convenient regime for the topological phase to arise.
In the SW formalism, the Hamiltonian is considered to be divided between a block-diagonal non-interacting part H 0 with a clear energy-scale separation, and some small interaction V , i.e.: where H = P ⊕ Q denotes the division of the Hilbert space such that in the sense that all eigenvalues from the P subspace are much smaller than the eigenvalues from the Q subspace, and that the matrix elements of the V operator are much smaller than the energy separation between the two subspaces. We recall that a possible term P P V P P can be set to 0 w.l.o.g.. The target SW transformation is a rotation X of the Hilbert space, such that the Hamilton operator is brought back to a block-diagonal form under its action: The desired effective low-energy Hamilton operator is then given by While X, and thus H eff , are only known exactly for a few special cases, there exists a perturbative solution in terms of powers of the interaction V , with terms in the typical form: and variations thereof, especially in case the original lowenergy subspace is not exactly degenerate, i.e., if not all states in P share the same eigenvalue E 0 under H 0 . Anyway, we can easily identify the Green operatorĜ(ω) restricted to the high energy space Q and evaluated at E 0 : a fact which will come handy in the following. In our specific setup, we chose i.e., the Hamiltonian acting on the auxiliary sites. Thereby, it is easy to identify the low-energy space P as the one containing all states with empty c and d sites, while the high-energy configurations Q are all the remaining ones with at least one fermion placed on these auxiliary sites. As a consequence, E 0 = 0 and the energy separation is of the order of µ. As long as t = t ⊥ = 0, H 0 is already diagonal in the Fock basis, and this allows for an exact evaluation ofĜ Q (E 0 ), leading to the compact expressions of Eqs. (8)-(9).
If we now include such terms, i.e. t , t ⊥ = 0, this is not true anymore since a fermionic excitation on the c and d states is now allowed to delocalize across the auxiliary sites. Fortunately, we can still evaluateĜ Q (E 0 ) if we restrict for a moment to the case in which at most one fermionic state in the c and d sites can be occupied. The matrix elements of this operator decay exponentially in real space. This in turn leads to an additional exponentially decaying hopping in the a, b chains of the form: where the −(+) holds for α = a(b). However, the fourthorder term in the Schrieffer-Wolff transformation, which generates the desired pair hopping term, requires to deal with two fermions on the c, d sites, and therefore to solve the full interacting problem. However an analytic solution is not that easy any more. Anyway, in the limit of t , t ⊥ µ (consistently with all other energy-scales of the setup), one can treat them as small perturbations and compute the Green operator perturbatively. The sizeable diagonal elements are responsible for the generation of the pair hopping terms, while the off-diagonal contributions are again exponentially suppressed. We finally arrive at the expression: where Γ j are 2 × 2 matrices, depending on all other parameters, with Γ 1 = Γ 2 , which allows for detuning In Fig 7 some examples for W j (t , t ⊥ ) are shown. Moreover, a third pair coupling operator is generated: This was also one of the reasons to consider the slightly more general toy model of equation (8) where one finds the relation W 3 = (t , t ⊥ )Γ 3 (t , t ⊥ ) T . The possibility of detuning W 1 relative to W 2 and also the generation of W 3 shifts the effective model near to the vicinity of the exactly solvable model [21], for which the relation figure 7, indeed, we see that introducing t increases W 1 relative to W 2 and also introduces a positive W 3 .

Appendix B: Bosonization Details
In this Appendix, we provide the details of the derivation of the bosonized low-energy theory and its bare couplings in terms of the microscopic parameters. The starting point is the effective Hamiltonian of Eq. (8). In order to keep the equations simple, we start here by using the simplified single chain Hamiltonian from the main text: However, towards the end of this appendix we will shortly discuss additional operators which are generated by the Schrieffer-Wolff transformation. The first step is to rewrite the lattice annihilation/creation operators α ( †) in terms of two (slowly-varying) envelope functions ψ R/L,α defining the right/left moving fields: The "continuum" position x j is defined as x j := aj, with a the lattice spacing, in terms of which we will express all quantities (like lengths, energies, etc.) in the following. The Fermi momentum is defined by k F = N π 2La , with N being the number of fermions and L the number of lattice sites in the system. In order to formulate a theory in the continuum, we send a → 0 while keeping constant the product ak F := δ = πν, as well as all energies ta, U α a, . . . This also amounts to replacing sums by integrals according to the rule a j → dx. Henceforth, we will also remove any residual explicit dependence on the lattice spacing by appropriately rescaling the fields and the coupling constants for better readability.
The next step is the assumption of linearity for the fermionic dispersion relation near the (two) Fermi points. The right and left moving fields become thus indepen-dent from each other, and the free fermion Hamiltonian is written as H 0 = α=a,b dx h 0,α (x): where v F = 2t sin(δ) defines the Fermi velocity. The independent left/right moving fields are then rewritten in terms of vertex operators of continuous bosonic fields ϕ α (x) and ϑ α (x), describing respectively the density and phase fluctuations: with s σ = +1 if σ = R (and −1 for L), and the latter relation defined in the Heisenberg picture with explicit time-dependent operators. It is particularly useful to recall also the expression for the current densities: The Klein factors η α,σ , forming a Clifford algebra (i.e., {η α,σ , η β,ρ } = 2δ α,β δ σ,ρ ), are essential to obtain the correct anti-commuting behavior of the fermionic operators [28]. However, thanks to the particle-number preserving character of the Hamiltonian, we can simply treat them as simple hermitian matrices and reorder their strings to be the same in all terms: henceforth, we consider to have already performed such a reordering and drop all Klein factors from our formulas. Subsequently, we rewrite the different lattice operators in terms of the bosonic fields, according to the dictionary presented in Eqs.(B4)-(B5), and taking care of normal ordered products along standard procedures [33]. One major consequence is that, in most cases, the algebra amounts to directly summing the exponents appearing in equation (B4) when dealing with products of ψ ( †) fields. Since all terms turns out to be diagonal in the bosonic fields, and we are considering identical (a, b) species, it is convenient to resort to symmetric and anti-symmetric combinations of the fields: The free Hamiltonian and the intra-chain interactions are thereby well-known to be mapped to a quadratic form: with v τ and K τ the Fermi velocity and the so-called Luttinger parameter in each sector. These are equal to v F and 1 in the free case, and get renormalized by the interactions. Indeed, the representation of the lattice number operator reads xj accounting for charge density waves. By integrating their product on neighbouring sites over the whole lattice, all oscillating terms will average out unless we are at half-filling, and we are left with: with the following coefficients: By applying the same procedure to the total unit-cell interaction, the inter-chain terms give rise to an extra Sine-Gordon interaction involving the field ϕ − , originating from scattering terms of the form: with coefficients: g nn,ϑ,+ = −g nn,ϕ,− = g nn,ϑ,− = U nn cos(2δ) , g nn,ϕ,+ = U nn (4 − cos(2δ)) (B13) Finally, we can also translate the pair-hopping terms by similar algebra, and obtain with coefficients: Noticeably, the two kinds of pair-hoppings give rise to Sine-Gordon terms for the two conjugate bosonic fields: it will be the one in ϑ − which will be responsible for the topological phase, while the one in φ − is already present with other types of density interactions between the two chains. The Sine-Gordon term involving φ − favors the formation of either a spin density wave or charge density wave, depending on the sign of the coupling [33]. By putting all these contributions together, we get to Eq. (11) with the bare couplings of the low-energy theory expressed in terms of those of the microscopic lattice Hamiltonian: (B17) In order to determine the actual phase the system will end up into, such bare couplings should be analysed from the renormalization group (RG) perspective, i.e., by integrating out short-distance degrees of freedom and retaining the long-distance ones only, thus moving from a full quantum action to a low-energy effective one.

RG Flow
The RG-flow is controlled by the flow parameter l, representing the effective momentum cutoff in terms of the overall ultra-violet one via Λ UV /Λ 1 + dl. At firstorder, the equations for the Sine-Gordon couplings β k , with k = 1, 2, are determined by their scaling dimension ∆ k : If a coupling flows to ∞ for l → ∞, then the theory acquires a gap, and the coupling is dubbed relevant: this happens if ∆ k < 2. One should actually stop the flow when the value overcomes the cut-off, and could then estimate thereby the value of the gap. If instead the scaling dimension is large, ∆ k > 2, the coupling is irrelevant, since it flows to 0 and disappears from the effective theory. The limiting case, ∆ k = 2, the coupling is labeled as marginal, and higher orders are required to find out its actual behaviour. In our Eq. (11) we find the common result ∆ 1 = 2/K − and ∆ 2 = 2K − [33]: while K + does not flow at all, and the symmetric sector remains gapless in all cases, we should resort to higher orders of perturbation theory to inspect the flow of K − , at least around where A is some cutoff depending constant. No additional contribution to the flow of the β k couplings is generated at second order, and thus the set of equations (B18)-(B19) is consistent. Similar RG flow equations have been studied in the past, and it has been shown that all points on the plane defined by flow to a critical model [31]. Upon inserting the bare values of K − and β j in terms of the original lattice couplings, Eq. B17, and linearising the dependence of v − and K − on small W 2 values, the criticality condition can be recast, for U nn = 0, as with two non-universal constants D = π 2 /(2 sin 2 δ) and C = cos(2δ)/(π 2 √ A) (in the case of U α = 0). Therefore, we expect a pretty different behaviour depending on the sign of W 2 , while a symmetry in W 1 should appear. For example, inserting all numbers we would expect a slope α + = 1 (α − = −1/3) in the case of W 2 > 0 (W 2 < 0) which perfectly matches the observed slope in Fig. 2a.
The same equation can also be used to predict the behavior of the critical lines in the W 1 = W 2 = W and U α phase diagram of figure 2b. However, the equations are not compact and easy to write. But for the choice of the non-universal constant A used for the numerical integration one finds a leading-order linear behavior for W < 0 and a quadratic leading-order for W > 0 matching the numerical observations. We tested these predictions by numerically integrating the differential equations (B18)-(B19) starting from the bare values of the couplings, up to a point where one of the two β k coupling constants reach a certain cut-off value β k (l * ) = c. This indicates the formation of a spectral gap of one or the other kind, which can be estimated according to ∆ ∼ e −l * . The precise predictions depend on the non-universal constant A appearing in the flow equations, too: Nevertheless, we can use them for a rough estimation of the phase diagram, presented in Fig. 2. The (asymmetric) linearity of the boundaries is evidently kept up to fairly large values of the W couplings.

Two Kitaev Chains
In the context of the adiabatic connection between our model and the situation of two independent Kitaev chains (which will be deepened in App. C), it is worth briefly mentioning the bosonization of the pair-hopping operator acting equally on the two chains: By using the recipe exposed above, one finds: with g ∆ = 2∆ sin(δ)/π. In a situation where the ϑ − field is already locked, as in the topological phase through the H W1 operator, H ∆ is basically the operator cos( √ 2πϑ + (x)). The scaling dimension of this operator is given by ∆ = (2K + ) −1 , i.e., as long as K + > 1/4, this operator becomes relevant and gaps out the charge sector.

Irrelevance of Additional Effective Terms
As promised at the beginning of this long Appendix, we want to justify our assumption of dropping all the additional terms which appear in the effective single chain Hamiltonians H α . As a starting point, the single chain Hamiltonians of the full model are chosen as the usual spinless Fermi-Hubbard Hamiltonian consisting of a hopping term and a nearest-neighbor interaction term: The commutators K α,j = [H α , α j ]/t can be readily computed to be: and inserted into Eq. (10): while α † j K α,j + h.c. will only redefine the effective values of t and U , the product K α,j K α,j will also generate three-body terms (n j−1 n j n j+1 ).
In general, the lowest order of an operator consisting of N fermionic densities is given in bosonization by a power-N operator ∼ (∂ x ϕ α ) N . The scaling dimension of these operators can be shown to be ∆ N = N , i.e., their flow equations are of the form dβ N /dl = (2 − N )β N . Thus, these operators become surely irrelevant for all N ≥ 3. In addition to these N -power operators, also some higher harmonic cosine terms might appear, i.e. cos(n √ 4πϕ α ).
However, their scaling dimension is a monotonic increasing function of n, meaning that the most relevant operator is given by the first harmonic n = 1. Nevertheless, with increasing interaction strength these higher harmonics might become relevant, if K ± 1: however, this is by far not the scenario we are considering in this paper.

Appendix C: Adiabatic Connection to Exactly Solvable Models
In this Appendix, we provide details about the two paths in parameter space we chose for illustrating the adi- U nn = W 1 /2. However, since the latter lives in a higher dimensional space, this is to be understand as a cartoon, rather than an exact statement.
abatic connection between our effective model in Eq. (8) and i) the exactly solvable one of Ref. [21], or ii) the setup with two uncoupled non-interacting Majorana chains. In Fig. 8 we also show a cartoon picture of how all the different number-conserving models of this article are related. For the first one, whose non-vanishing gap is plotted in Fig. 9a, we chose where h.c. acts over every term which is not already explicitly hermitian. In the Hamiltonian above, the chemical potential µ half of the strength for the very first and very last site compared to the bulk sites and compared to eq. (8) an additional interaction term between the two chains U r is introduced. This is necessary in order to solve the model exactly at the special point considered in [21]. However, we explicitly checked that the existence of the Majorana-like phase does not depend on the lowering of the onsite potential for the first and last site, as one would expect for a topological phase.
For the second case, we define a simple linear interpolation, s ∈ [0, 1], between the two limiting cases: with H eff the one in Eq. (8) and H kitaev,α the Majorana chain for the α = a, b species defined at the sweet spot: The correlation length extracted from the single particle correlation functions stays finite along the interpolation path, as shown in Fig. 9b. The big drop in the beginning can be explained by the charge sector gapping out. This is supported by looking at the entanglement entropy (not shown here) which becomes asymptotically constant instead of following the logarithmic law of critical systems [51]. This behavior is expected since adding the pair potential terms allows for coupling of states with all possible particle numbers.
The existence of such an adiabatic connection, preserving the time-reversal symmetry and a (Z 2 ) + × (Z 2 ) − subgroup of the full symmetry group U(1) + × (Z 2 ) − of our model, is instrumental to understand and categorise the topological phase. Breaking the U(1) symmetry leads to a four fold degeneracy, differently from the case of preserving the U(1) symmetry. This can be understood by recognizing that the four ground-states split into two ground-states for each parity of the total particle number. By conserving the total paritcle number, and therefore fixing the parity, we restrict the model to one of the two subspaces having either a even partiy (N tot mod 2 = 0) or odd parity. This results in an effective two fold groundstate degeneracy as observed in the DMRG simulations, see Fig. 5.
Since the model we are dealing in this work is of interacting nature and the interaction is crucial to obtain the desired Majorana-like phase, it should be noted that the general Z classification of non-interacting fermionic systems in spatial dimension one breaks down to a Z 8 classification presence of interactions and time reversal symmetry, as was shown in Ref. [52]. As long as we only couple two chains by the effective Hamiltonian (8) such a distinction is not relevant.

Appendix D: Entanglement Spectrum Analysis
In this Section we briefly discuss the form of the entanglement spectrum within the topological phase, as depicted in Fig. 4a. From the bosonization analysis of our model we learned that the symmetric and antisymmetric sectors of the Hilbert space decouple. A similar decoupling is therefore also expected for the entanglement spectrum and the states corresponding to the Schmidt decomposition: δN,P,j = −2 log(λ δN,P,j ) Here, λ δN,P,j denotes the Schmidt values labeled by two quantum numbers associated to the symmetric (δN ) and antisymmetric sector (P ), namely the excess charge with respect to average filling and the parity of one of the dressed chains. Since the symmetric sector is gapless and the anti-symmetric one is in a gapped topological phase, we expect that the spectrum will show distinctive features of both. The gapless charge sector is indeed displaying the universal behavior with entanglement levels arranged in parabolas as a function of the quantum number δN , whose details are determined by the underlying conformal field theory [38,39].
This special form of the entanglement spectrum is also useful for extracting the Luttinger parameter K + , which was also used in the main text [40]: π sin(π/Ll) .
Moreover, the curvature of the parabolas is also mainly determinated by the Luttinger parameter K + by [39]: However, the exact numerical values of the entanglement spectrum and all correct degeneracies are hard to extract, since they are subject of strong finite bond dimension effects [38] and using the formula connecting the total number fluctuation to the Luttinger parameter is more stable.
The topological character of the gapped antisymmetric sector is dictating the presence of two copies of each level, transforming differently under the parity. This is in perfect agreement with the results of Ref. [29], and the degeneracy represents the fractionalization of the fermionic parity operator at the end of a finite subsystem, as more generally known for symmetry-protected topological phases [53]. Comparing the results for the system in the topological phase, Fig. 4a, with results form the system being in the trivial phase, Fig. 4b, we indeed find that the parabolas originating from the gapless charge sector are still present, however the non-trivial double degeneracy between the two different parity sectors is gone.

Appendix E: The Full Model in the Perturbative Regime
In this Section we report the results found for the full model deep inside the perturbative regime. The Schrieffer-Wolff transformation used for deriving the effective Hamiltonian (9) is valid in the limit of µ (the chemical potential on the c and d states) being the dominant energy scale. Together with requiring the resulting parameters of the effective model defined by the relation (9) being in the topological regime gives some additional constrains on the parameters of the full model.
For example, targeting the effective parameters U nn ≈ 0, W ≈ t/1.4 and U α = −0.5t while fixing µ/J = 5 , V 1 /J = −1 , V 2 /J = 1.25 , t = 0.01 leads to W/J ≈ 0.007 and U α /J ≈ −0.005. Thus, the intra-wire interaction U α is several orders of magnitude smaller than the interactions on the c and d states (V 1 and V 2 ), which are not constrained to be small since they act exclusively on the c and d subspace. Further, we used t = 0.05 together with t ⊥ = 0. This should lead to a slightly detuning W 1 > W 2 favoring the topological phase as expected from appendix A. The results for a simulation of a system with 60 sites and a filling of n e = 40 (ν = 1/3 in the effective model) are shown in Fig. 10. First looking at the density profiles in Fig 10b, one sees a slightly decrease of the population on the a and b sites compared to the effective model. However, this is expected due to the additional c and d states. Now looking at the end-to-end correlation function in Fig. 10a, one sees the same characteristic behavior as in the effective model. I.e., an exponential decay towards the middle of the system together with an exponential revival showing a relative sign between the two parity sectors. Remarkably, the energy gap (correlation length) of the full four-flavor model seems to be larger (smaller) than in the effective model. This can be explained by the additional terms in the Hamiltonian 8 which are discarded in the numerical simulations of the effective model, as discussed below Eq. 10.