Exploring Interacting Topological Insulators with Ultracold Atoms: the Synthetic Creutz-Hubbard Model

Understanding the robustness of topological phases of matter in the presence of strong interactions, and synthesising novel strongly-correlated topological materials, lie among the most important and difficult challenges of modern theoretical and experimental physics. In this work, we present a complete theoretical analysis of the synthetic Creutz-Hubbard ladder, which is a paradigmatic model that provides a neat playground to address these challenges. We put special attention to the competition of correlated topological phases and orbital quantum magnetism in the regime of strong interactions. These results are furthermore confirmed and extended by extensive numerical simulations. Moreover we propose how to experimentally realize this model in a synthetic ladder, made of two internal states of ultracold fermionic atoms in a one-dimensional optical lattice. Our work paves the way towards quantum simulators of interacting topological insulators with cold atoms.

Topological features of quantum many-body systems provide a new paradigm in our understanding of the phases of matter [1], and give rise to a promising avenue towards faulttolerant quantum computation [2]. From a condensed-matter perspective, such features lead to exotic groundstates beyond the conventional phases of matter, which are typically understood by the principle of symmetry breaking and the notion of a local order parameter. On the contrary, these exotic states can only be characterised by certain topological properties.
The integer quantum Hall effect, which is a paradigmatic example of such peculiar phases [3], requires the introduction of a topological invariant to describe the different plateaus and their associated transverse conductivities [4]. Another interesting property of this state of matter is the bulk-boundary correspondence, which relates such a topological conductivity, a bulk property, to the existence of current-carrying edge states localised within the boundaries of the system [5]. Although the bulk of an integer quantum Hall sample appears as a trivial band insulator, its boundary corresponds to a holographic chiral liquid where interactions merely renormalise the edge-state Fermi velocity [6].
As realised in a series of seminal works [7,9,10], these remarkable properties are not unique to quantum Hall samples subjected to strong magnetic fields. Instead, they arise in various models with different symmetries and in different dimensions [11], the so-called topological insulators and superconductors [12], which also lead to the notion of symmetryprotected topological phases in the context of topological order [1]. Remarkably enough, some of these models have turned out to be accurate descriptions of real insulating materials [12][13][14], and promising candidates to account for observations in proximitized superconducting materials [15]. This has positioned the subject of topological insulators and superconductors not only at the forefront of academic research, but also at the focus of technological applications, such as topo-logical quantum computation with Majorana fermions [2].
Despite this success, (i) there are still several paradigmatic models of topological models whose connection to real materials still remains unknown, or even seems quite unlikely, as it occurs for the Hofstadter model with magnetic fluxes on the order of the flux quantum [16], or the Haldane model [7]. Moreover, (ii) most of the topological materials explored in the laboratory so far do not display important electronic correlation effects [17]. This is rather unfortunate in view of the richness of the fractional quantum Hall effect [18], where such correlations are responsible for a plethora of exotic topological phases of matter.
In the present work, rather than considering real materials, we shall be concerned with the so-called synthetic quantum matter in atomic, molecular and optical (AMO) platforms, more particularly with ultracold gases of neutral atoms trapped in periodic potentials made of light, i.e. optical lattices [19]. The ever-improving experimental control over these quantum many-body systems has already allowed to design their microscopic Hamiltonian to a great extent. In this way, it is possible to target interesting condensed-matter models, some of which are still lacking unambiguous experimental feedback from experiments with real materials [20,21], as occurs for the bosonic [22][23][24] and fermionic [25][26][27] Hubbard models. The initial interest in optical-lattice implementations of integer quantum Hall phases [28], and other time-reversal invariant topological insulators [29,30], has risen considerably in the recent years due to the experimental progress [31,32]. In particular, elusive topological systems such as the Hofstadter [16] or the Haldane [7] models, have been realised in experiments with bosonic [33] and fermionic [34] ultracold atoms, respectively. The possibility of including two or more different atomic states/species, controlling their interactions via Feshbach resonances [35], may eventually lead to the experimental test of correlation effects in these topological states.
Here, we introduce a variant of the quasi-one-dimensional (quasi-1D) Creutz topological insulator [36], and study the effect of repulsive Hubbard-type interactions on the topological phase. This model, which shall be referred to as the imbalanced Creutz-Hubbard ladder, is readily implementable with ultracold fermions in intensity-modulated optical lattices. We argue that such a model has all the required ingredients to become a workhorse in the study of strongly-correlated topological phases in AMO setups. Indeed, two accessible AMO ingredients as (i) a simple Zeeman shift between the atomic internal states, which yields a leg imbalance in the ladder, and (ii) on-site repulsive contact interactions tuned by some Feshbach resonances, which lead to Hubbard-type rung interactions in the ladder, can be employed to access the rich phase diagram of this model (see Fig. 1), which was widely uncharted prior to our study. Starting from a flat-band regime, we show that the imbalance and the interactions lead to a competition between a topological phase and two different phases of orbital quantum magnetism. At large interaction strength, a long-range in-plane ferromagnetic order arises, related to the symmetry-broken phase of an orbital quantum Ising model; while the Zeeman imbalance then drives a standard quantum Phase diagram displaying a topological insulator (TI) phase, and a pair of non-topological phases: an orbital phase with long-range ferromagnetic Ising order (oFM), and an orbital paramagnetic phase (oPM). The horizontal axis represents the ratio of the inter-particle interactions to the tunneling strength, whereas the vertical axis corresponds to the ratio of the energy imbalance to the tun- phase transition in the Ising universality class towards an orbital paramagnetic phase. In order to understand the model away from this limit, we introduce two new methods based on mappings onto models of quantum magnetism and quantum impurity physics. These methods allow us to locate exactly certain critical points/lines, and to predict topological quantum phase transitions for weak and intermediate interactions with different underlying conformal field theories (CFTs), i.e. Dirac versus Majorana CFTs, which then fit very well with numerical results based on Matrix-Product-States (MPS) [37]. We also provide suggestions for experimental observables to pinpoint these three phases.
This manuscript is organized as follows: in Sec. II the standard Creutz ladder is introduced, and some previous studies are briefly accounted for, which show that this model leads to a topological insulator in the BDI symmetry class [11]. Therefore, the standard Creutz ladder lies in the same symmetry class as the Su-Schrieffer-Heeger model [38], which has already been implemented in optical lattices [39]. We then introduce the imbalanced Creutz ladder, which provides an instance of a topological insulator in the AIII class (chiral unitary), which still lacks an AMO implementation. Moreover, as argued in Sec. III, the Hubbard interactions lead to a very neat interplay of strongly-correlated and topological effects. We identify the different phases of the model, together with standard and topological quantum phase transitions that connect them. Finally, in Sec. IV, we lay out a detailed proposal to implement the imbalanced Creutz-Hubbard ladder in an ultracold Fermi gas with two different internal states trapped in a standard one-dimensional (1D) optical lattice, enriched by intensity-shaking and Raman laser-assisted tunnelling. We present our conclusions and outlook in Sec. V.

II. THE IMBALANCED CREUTZ-HUBBARD LADDER
The standard Creutz model describes a system of spinless fermions on a two-leg ladder (see Fig. 2 (a)), which are created-annihilated by the fermionic operators c † i, , c i, , where i ∈ {1, . . . , N} labels the lattice sites within each leg of the ladder ∈ {u, d}. Fermions are allowed to hop vertically along the rungs of the ladder with tunnelling strength t v , and horizontally along the legs of the ladder with a complex tunnelling where t h is a tunnelling strength, and δ a,b is the Kronecker delta. The arrangement of complex phases in the horizontal links leads to a net 2θ -flux gained by a fermion hopping around a square unit cell, playing thus the role of the so-called Peierls phases of a magnetic field piercing the ladder. In addition, the kinetic part of the Hamiltonian also includes a diagonal tunnelling of strength t diag , yielding altogether This quadratic lattice model was put forth in Ref. [36] as a simple toy model to understand some of the key properties of higher-dimensional domain-wall fermions [40,41], which were introduced in the context of lattice gauge theories to bypass the fermion-doubling problem [42]. For periodic boundary conditions, this model leads to a couple of bands that display a pair of massive Dirac fermions with different Wilson masses m 0 , m π at momenta k D ∈ {0, π} [43]. For open boundary conditions, a pair of zero-energy modes exponentially localised to the left/right edges of the ladder appear as one of the Wilson masses gets inverted (m π < 0) when t v < 2t diag . Considering the bulk-boundary correspondence discussed in the introduction, these edge states resemble the holographic liquid of the higher-dimensional topological insulators. In fact, the change in polarisation of the system can be characterised by a topological invariant [44], the so-called Zak's phase [45], such that the appearance of these zero-energy modes coincides with a non-vanishing topological invariant, and the Creutz ladder yields a symmetry-protected topological phase in this regime. As discussed below, for θ = π/2, this topological phase corresponds to a BDI topological insulator.
Since the objective of this work is to study correlation effects, we now consider the simplest possible Hubbard interactions between the spinless fermions where V h (V v ) are the density-density interaction strengths between fermions residing in neighbouring sites along horizontal (vertical) bonds of the ladder, and we have introduced the fermion number operators Figure 2: Standard and imbalanced Creutz ladder: Two-leg ladder where fermions tunnel along the black links enclosing a net flux 2θ along a closed plaquete: (a) standard Creutz ladder (1), and (b) imbalanced Creutz ladder, which leads to Eq. (4) in the π-flux limit.
For reasons that will become clear in the cold-atom implementation discussed in Sec. IV, in the following we will deal with a variant of the Creutz Hamiltonian: (i) we shall substitute the vertical tunnelling by an energy imbalance between the legs of the ladder ε u = ∆ε/2 = −ε d , which changes the symmetry class of the topological insulator for θ = π/2 from BDI to AIII; (ii) we limit the interaction terms (2) to the anisotropic regime V h = 0; (iii) we set the amplitude of the diagonal hopping equal to the one along the legs (|t diag | = |t | = t) and finally (iv) we fix the phases in order to get a net π-flux through the plaquettes. The resulting Hamiltonian (see Fig. 2 (b)), which we will refer to is the imbalanced Creutz-Hubbard Hamiltonian, is where we have introduced the kinetic term with s u/d = ±1. This term leads to a pair of flat bands, and a couple of zero-energy topological edge states. The remaining terms contain the Hubbard interactions and the energy imbalance, which can have non-trivial effects on the flat-band physics, and induce a phase transition to other non-topological phases of matter. In the following section, we present a new formalism to understand such transitions.
We are encouraged to study this imbalanced Creutz-Hubbard ladder both for experimental and theoretical reasons that will be discussed in detail in the forthcoming sections. In particular, let us advance that this model shares topological properties with the original Hamiltonian (1). In particular, the Dirac fermions now occur at k D ∈ {−π/2, π/2} with different Wilson masses m ±π/2 . Provided that ∆ε < 4t h = 4t d , the mass m π/2 < 0 gets inverted, and we obtain analogous topological features and exponentially-localised edge states. The choice of the anisotropic regime V h = 0 is motivated by the use of ultracold fermionic atoms in optical lattices with contact interactions. Exploring also the regime of V h > 0 would require situations where the atoms have longer-range interactions.

Previous studies on related models
With the model under study been defined, let us comment on some relevant literature, and advance some comments in relation to our results. The standard Creutz ladder (1) filled with bosons has been studied in Refs. [46,47], which includes on-site Hubbard interactions instead of the nearest-neighbour terms of Eq. (2). The focus of these two articles was the appearance of pair superfluids due to the interplay of interactions and frustration flat-band effects [48]. A useful formalism was introduced in these works, which consisted on the projection of the Hamiltonian onto the lowest-energy flat band by introducing highly-localised Wannier functions. Moreover, using a unitarily-equivalent formulation of the Creutzladder Hamiltonian introduced in Ref. [47] would lead to a spin-orbit coupled Hubbard model with staggered energy imbalance in our case (3), which might broaden the relevance of our results beyond the ladder compound of Fig. 2. The standard Creutz ladder (1) populated by spinful fermions with on-site interactions between opposite spins has been studied recently in Ref. [49], which focused on an exact Bardeen-Cooper-Schrieffer description for attractive Hubbard interactions. To the best of our knowledge, the standard Creutz ladder with spinless fermions has only been studied previously in Ref. [50]. Here, the emphasis was placed on an additional superconducting s-wave pairing, and its interplay with the repulsive Hubbard interactions, which leads to an interesting Creutz-Majorana-Hubbard ladder. The limit of vanishing pairing, which corresponds to the standard Creutz-Hubbard model with the vertical tunneling (1) instead of the imbalance (5), and leads to the aforementioned BDI topological insulator, was only touched upon briefly by numerical meanfield and density-matrix renormalisation group studies. These results pointed towards the possibility of a phase with in-plane ferromagnetic order as the interactions are increased.
To the best of our knowledge, the effect of Hubbard interactions in the imbalanced Creutz ladder (3), and more generally on AIII topological insulators, remains largely unexplored. In relation to the above studies on related models, we show below that one cannot project onto a single flat band for the fermionic model, but must instead retain all flat bands and edge modes to account for correlation and topological effects accurately. In this way, we shall make interesting connections between the existence of topological phase transitions and the physics of quantum impurity models. Additionally, by going beyond a mean-field analysis in the strong-interaction limit of the imbalanced model (3), we show analytically that in-plane ferromagnetic order also arises in our model, and moreover corresponds to the symmetry-broken phase of an orbital quantum Ising model. In addition, we illustrate that the imbalance drives such Ising model into a orbital paramagnetic phase.

III. TOPOLOGICAL QUANTUM PHASE TRANSITIONS
The following subsections are devoted to the construction of the phase diagram shown in Fig. 1. We start by discussing the solution of the non-interacting imbalanced Creutz ladder, and the appearance of flat bands and fully-localised edge states in the Hamiltonian (4) (see Sec. III A). This corresponds to the vertical axis of the phase diagram. In Sec. III B, we examine the weakly-interacting regime and show that the model maps onto a pair of weakly-coupled Ising chains, which can be studied through a mean-field analysis (i.e. the region in the vicinity of the vertical axis of Fig. 1). In Sec. III C, we study the opposite limit of very strong interactions (i.e. rightmost region of Fig. 1), and discuss the possible non-topological orbital magnetic phases that can arise. In Sec. III D, we explore the intermediate regime, and show that the effect of the interactions and imbalance leads to edge-bulk couplings that can be mapped onto quantum impurity models. This new perspective yields a neat picture underlying the destabilisation of the topological phase in favour of the orbital magnets. These different methods allow us to build an analytical prediction of the phase diagram of the model. Finally, in Sec. III E, we test numerically the above predictions, and provide a detailed study of the phase diagram by means of Matrix-Product-State numerical simulations. Scattered through these sections, we shall also introduce comments on possible experimental tools to measure the relevant observables for the different phases of the model, which will become relevant for the specific experimental cold-atom proposal discussed in Sec. IV.

A. Non-interacting limit: Flat bands and edge states
We start by solving the kinetic part (4) of the π-flux Creutz-Hubbard Hamiltonian (3). For periodic boundary conditions, and after introducing the spinor Ψ(q) = (c u (q), c d (q)) t for the fermion operators in momentum space c (q) = ∑ i e −iqai c i, / √ N, one finds where σ σ σ = (σ x , σ y , σ z ) is the vector of Pauli matrices, and B B B(q) = 2t(− cos(qa), 0, sin(qa)). By direct diagonalization, one finds that the system develops two flat bands ε ± := ε ± (q) = ±2t, where q ∈ BZ = [−π/a, π/a) is the quasimomentum, and the lattice constant shall be set to a = 1 henceforth. The vanishing group velocity associated with these bands, v g = ∂ q ε ± (q) = 0, indicates that the groundstate must be insulating, regardless of the particular filling. This can be considered as a new type of insulator, namely a flatband insulator, which corresponds neither to the usual band insulator, nor to the Mott insulators. It shares some properties with the former (i.e. no correlations), and with the latter (i.e. localized fermions), but it differs from both insulators in the large degeneracy of the ground state, except for half-filling conditions. On top of this, the flat bands are also topological: the diagonalization of the discrete chiral symmetry σ y (s.t. σ y H(q)σ y = −H(q), with H(q) = B B B(q) · σ σ σ ) puts the Hamiltonian (6) in a purely off-diagonal form, with elements B x ± iB z ; its complex phase gets a non-trivial winding number W = sgn(t) = 0 [52]. Equivalently, we could consider the eigenvectors |ε ± (q) ∝ (B x + iB z ) 1/2 , ±(B x − iB z ) 1/2 t and realize that they exhibit a uniform Berry connection A ± (q) = i ε ± (q)|∂ q |ε ± (q) = 1 2 . The uniform Berry connection leads to a finite Zak's phase, which is defined [45] as and equals ϕ Zak,± = π for our topological flat bands. This Zak's phase pinpoints the topological properties of the bands, and can be connected to a macroscopic observable: the polarization of the system [44]. Interestingly enough, the Creutz ladder displays an infinite flatness parameter without requiring long-range tunnelings, as necessary in higher-dimensional models of topological flat bands [53]. From this perspective, switching on the leg imbalance ∆ε > 0 in Eq. (5) leads to some curvature in the energy bands where we have introduced the flatness parameter f = 4t/∆ε, which becomes infinite for vanishing imbalance. The presence of the imbalance also drives the system out of the BDI class and into the AIII class [11], since B z (q) = ∆ε/2+2t sin q acquires a mixed parity under k ↔ −k and therefore no effective time-reversal (U T H(−q) * U † T = +H(q)) or particle-hole (U C H(−q) * U † C = −H(q)) operators can be found. Anyway, the discrete chiral symmetry is still described by σ y and the whole procedure described above can be employed. The Berry connection of the bands becomes non-uniform due to their curvature which leads to the following Zak's phase ϕ Zak,± = πθ (f − 1), where θ (x) is the Heaviside step function. Hence, the Zak's phase yields topological effects provided that the bands are sufficiently flat, i.e. f > 1. Conversely, when ∆ε > 4t, the curvature of the bands is large, i.e. f < 1, and no topological phenomenon occurs. We note that this point ∆ε = 4t is exactly the regime where one of the fermionic Wilson masses vanishes, as mentioned in the previous section in connection to the domain-wall fermions of lattice gauge theories. This marks a quantum phase transition between the AIII topological insulator and a trivial band insulator, as shown in the vertical axis of Fig. 1.
Regarding cold-atom experiments, we note that the Zak's phase has been measured for another paradigmatic topological insulator in 1D: the Su-Schrieffer-Heeger model of polyacetilene [38]. In this non-interacting limit, the Zak's phase is a single-particle property of the topological bands, and can be thus measured by a Ramsey interferometric protocol that exploits Bloch oscillations of single-particle initial states [39]. This measurement scheme has been generalised to other topological insulators [54], and can also be applied to the cold-atom implementation of the Creutz ladder considered in Sec. IV.
To have an alternative view on these topological features, let us introduce the so-called Aharonov-Bohm cages [55], which will also become very useful once interactions are switched on. In the π-flux Creutz ladder (4), the fermions cannot tunnel two sites apart due to the Aharonov-Bohm effect [36] (see Fig. 3(a)). One can thus find single-particle eigenstates strictly localised in cages formed by simple square plaquettes (see Fig. 3(c)). In second quantisation, such Aharonov-Bohm cages (AB-c) with energies ε ± = ±2t are For these particular weights, the fermions cannot tunnel one site apart due to the Aharonov-Bohm effect (see Fig. 3(a)), and are thus localised within a boundary AB-c (see Fig. 3(b)), which corresponds to an edge state within the bulk-edge correspondence of the topological insulator. We can finally express the π-flux Creutz Hamiltonian (3) for open boundary conditions as Although generic fillings can lead to a variety of interesting phases in the presence of interactions, potentially connected to fractional topological insulators [53], we shall only be concerned in this work with half-filling, i.e. N fermions, where the groundstate of Eq. (12) is two-fold degenerate and the groundstate energy is We thus see that the groundstate degeneracy corresponds to the two possible choices in populating either of the zeroenergy edge modes, and it is related to the topology of the ladder (i.e. ring versus line with open edges).  (4), one can identify two type of rhombic plaquettes that enclose a synthetic π-flux (left). Therefore, a particle trying to tunnel two sites apart (middle) will be subjected to a destructive Aharonov-Bohm interference that forbids this process. Accordingly, particles are confined to the so-called Aharonov-Bohm cages and cannot spread through the entire lattice. These cages correspond to square plaquettes, except at the edges (b), where destructive interference can also be found for a particle trying to tunnel one site apart (right). (c) Aharonov Bohm cages with the relative amplitudes for a single-particle state in two possible flat bands, and in the two possible zero-energy edge modes.
The effects of the leg imbalance ∆ε > 0 in (5) can be understood from this edge perspective by writing where t imb = −i∆ε/4 is an effective tunnelling induced by the imbalance, and has two relevant effects. The first line describes the hopping of fermions in neighboring AB-c, which leads to the aforementioned curvature of the bulk energy bands (8). The second line represents the hopping between the topological edge modes and the bulk AB-c. As discussed in more detail in Sec. III D, both terms conspire to induce a broadening of the edge modes in the regime 4t ≤ ∆ε, which is the regime where the topological Zak's phase (7) vanishes, signaling an imbalance-induced topological phase transition.
We also note that, precisely at the point 4t = ∆ε, the bulk bands become ε ± (q) = ±2t(1 + sin q), and the gap vanishes exactly at q = −π/2, which coincides with the momentum of the massless Wilson fermion. Regarding cold atoms, previous experiments on synthetic two-leg ladders subjected to artificial gauge fluxes, where the diagonal inter-leg tunneling in Eq. (1) is substituted by a vertical one, have measured non-vanishing chiral currents circulating in opposite directions along each leg of the ladder [56]. These states are connected to the edge states of the Hofstadter model as the number of legs is increased [56,58], but the bulk and edges coincide in the limiting case of the two-leg ladder. The situation is different for the two-leg Creutz ladder, as the edge states are not extended along the legs of the ladder, but exponentially localised to the left-and right-most rungs (11). We note that observing the particular localisation in Eq. (11) would require implementing a box confining potential, which have already been achieved in several cold-atom experiments for bosonic gases [57]. However, we note that for milder confining potentials, one expects that the edge states will remain to be confined in the boundaries of the system, provided that the confining potential increases at least quadratically with the distance to the center of the trap [59].
According to the above discussion, no additional legs would be required to identify the difference between edge and bulk excitations in the Creutz ladder, such that detecting the presence of zero-energy modes localised at the edges would be a proof of the topological nature of the phase. Using a pair of laser beams with a frequency difference that can be scanned within the edge-bulk gap |ω L,1 − ω L,2 | < 2t, and a spatial profile that can be localised to a few sites from the leftand right-most rungs of the ladder, the associated Bragg signal due to Raman excitations can detect the presence of these localised zero-modes, provided that the atomic gas is confined in a box potential. This method is analogous to the situation in higher-dimensional topological insulators [60], where additional angular momentum of the laser beams can be exploited to probe the chirality of the edge modes and improve the measurement resolution [59]. For milder confinement potentials, non-topological states might also happen to be localised within the edges of the system [60]. In Reference [59], the authors show that the Bragg scheme can be adapted to distinguish the chirality of the localized edge states of a higherdimensional topological insulator, and thus differentiate them from the non-topological states, probing in this way the topological nature of the phase. It would be interesting to explore analog techniques to differentiate the edge states (11) of the Creutz ladder from spurious states that would arise if no box potential is implemented.

B. Weak interactions: Quantum Ising ladder
Let us now address the fate of this topological phase as the Hubbard repulsion is switched on. We start by exploring the weakly-interacting regimet V v , and by noting that the two representations of the non-interacting imbalanced model , both in the original (Eqs. (1) and (2)) and in the plaquette (Eqs. (12) and (15)) bases, are dual to each other. In particular, the mapping between spinors induces a duality transformation H πC (t, ∆ε) → ∆ε 4t H πC t, 16t 2 ∆ε with a self-dual point ∆ε = (4t) 2 /∆ε corresponding to the previous critical point ∆ε = 4t. Such a duality is reminiscent of the one occurring in the quantum Ising model (QIM) [65], and suggests a possible equivalence between both models. Indeed, a formal equivalence is found by introducing the following rung operators Under this canonical transformation, the Hamiltonian is transformed onto In this particle-hole rung basis (17), we identify two independent subsystems which no longer display particle-number conservation, but instead have parity conservation. A Jordan-Wigner transformation [66], namely reveals the Ising nature of the two subsystems, and leads to a Hamiltonian that can be understood as a two-leg quantum Ising ladder For open boundary conditions, this description allows an alternative interpretation of the edge-state behaviour by writing the Ising model as a Kitaev-Majorana chain [9]. The ferromagnetic regime is associated with two uncoupled Majorana zero-energy modes on the opposing edges of the system for each of the legs of the Ising ladder. Combining the two free Majoranas on either edge of the Ising-ladder then yields the local fermionic edge modes defined in Eq. (11). This contrasts the case of a single Majorana chain, where the fermionic zero mode is highly non-local since it can only be built from the Majoranas at the opposite edges of the chain. Let us turn on the interactions, and express the weak Hubbard repulsion in terms of these rung operators These terms introduce a coupling between the two legs of the quantum Ising ladder. If we restrict to half-filling in the original model, the tunneling term vanishes, and we are left with the quartic interactions. These in turn can be expressed in terms of spin-spin interactions that couple the two legs of the ladder For weak interactions, V v ∆ε,t, we can treat the influence of one Ising chain on the remaining one with a mean-field approximation where we have introduced the transverse magnetization mn(∆ε,V v ,t) = σ z j,n for each leg of the ladder, andn = 2, 1 for n = 1, 2. We thus observe a renormalization of the imbalance parameter that controls the transverse field of the Ising model, and thus leads to a shift of the critical point as the interaction strength V v increases.
Accordingly, the topological phase of Sec. III A survives in a finite region of parameter space as the interactions are switched on. We find this region by solving the selfconsistency mean-field equations by iterating to convergence. In this equation, we have used the exact expression of the ground-state transverse magnetization in the QIM [67], namely where we have introduced the complete elliptic integrals of the first and third kind, namely The self-consistent mean-field approximation reproduces very precisely the numerical results for the density imbalance of the Creutz ladder In Fig. 4, we compare the groundstate density imbalance ∆n with the mean-field transverse magnetization m 1 (∆ε,V v ,t) = m 2 (∆ε,V v ,t), and find a very good agreement even for considerable interactions V v ∼t.
If we solve the self-consistency equation to first order in V v /t, we find that the critical point ∆ε/t = 4 flows towards smaller values of the imbalance as the interactions increase This weak-coupling expansion defines a critical line in parameter space that separates the topological and non-topological phases, and agrees well with our numerical findings for the phase diagram of the model, as discussed below (see the red dashed line in Fig. 1).
Regarding the cold-atom experiment, we note that the definition of the Zak phase (7), which could be used to probe the topological nature of the phase via interferometric methods [54], is only valid in the non-interacting limit. For weak interactions, the single-particle excitations will be described instead by quasi-particles, and the interferometric protocol to measure a many-body Zak's phase can be obtained by coupling additional impurities to the interacting Fermi gas [68]. We also note that the method to detect the zero-energy edge modes based on Bragg scattering [60], should also hold in the interacting regime. Finally, we advance that a measurement of the leg imbalance density (28) could also be used to experimentally test the validity of our prediction of the critical line (29) (see Sec. III E). Concerning the mapping in Eq. (78), this leg imbalance can be inferred from spin-resolved density measurements obtained by optical imaging, either in-situ (dispersive imaging) or after a time-of-flight expansion (absorption imaging). Exploiting light polarization and selection rules for the particular internal states of the atoms, one could measure separately the density of each of the two components, and even spin structure factors [61] that give access to the leg-imbalance susceptibility. Note that this susceptibility can be used to extract the position of the critical line (see top panel of Fig. 6). Regarding in-situ absorption imaging, spin-selective density measurements and correlations can also be accomplished [62]. As alternative possibilities, the spin structure factor could be measured in a non-demolition manner by exploiting a quantum Faraday effect [63], while spin-resolved density measurements can be achieved by using recent quantum gas microscopes (see e.g. [64]).

C. Strong interactions: Orbital Ising ferromagnet
In this section, we explore the opposite limit of the Creutz-Hubbard ladder (3), namely the strongly-interacting regimẽ t V v . In the limitt = 0, the groundstate of the Creutz-Hubbard Hamiltonian corresponds to a Mott insulator where the N fermions are distributed in the ladder avoiding simultaneous occupancies of two sites within the same rung.
By switching on the tunnellingt V v , an orbital analog of the well-known super-exchange [69][70][71][72] takes place, which reduces the energy of the groundstate by allowing for processes where a rung of the ladder is virtually occupied by two fermions. Instead of the usual Heisenberg model associated to is reasonable even for considerably large interaction strengths. The green curve shows as a reference the magnetization of the non-interacting model M(∆ε/4t).
the strongly-interacting Hubbard model [69], a different spin model arises in the present case due to the combination of intra-and inter-leg tunnellings. To second order of perturbation theory [73], and considering the half-filling regime, the relevant Hamiltonian describing the low-energy physics corresponds to an orbital quantum Ising model, namely where the super-exchange coupling is J = −8t 2 /V v , and the leg imbalance ∆ε plays the role of an effective transverse field. The above spin operators are defined as the orbital analogue of the usual spin operators for electrons Finally, P r = Π i (1 − n i,u n i,d ) is a Gutzwiller projector onto the subspace of singly-occupied vertical rungs.
The 1D quantum Ising model can be exactly solved [74] by introducing a Jordan-Wigner transformation [66], followed by a fermionic Bogoliubov transformation [75]. In comparison with the non-interacting groundstate (13), which displays a topological two-fold degeneracy, the strongly-interacting groundstate for ∆ε < |J|/2 has a non-topological degeneracy related to the Z 2 symmetry of the Ising model. In this regime, the groundstate develops long-range order as a consequence of spontaneous symmetry breaking lim r→∞ T y i T y i+r SI = 1 where h = 2∆ε/|J| < 1. This defines a critical line that separates the phase of long-range order, i.e. an orbital ferromagnet (oFM), from the disordered phase, i.e. an orbital paramagnet (oPM), and is depicted by a yellow dashed line in Fig. 1. As the leg imbalance is increased above a critical value ∆ε t | c = 4t V v , a standard quantum phase transition occurs between the long-range ordered Ising ferromagnet and a disordered orbital paramagnet, where all fermions tend to occupy the lower leg of the ladder. We note that this transition is not of a topological origin, as it can be understood by a local order parameter: the orbital magnetisation T y i T y i+r SI → m 2 y . We also note that the long-range ferromagnetic order is totally absent in the non-interacting topological groundstate (13), where one finds According to this discussion, it is clear that one cannot connect the non-interacting topological and strongly-interacting ferromagnetic phases adiabatically. Therefore, there should be an interaction-induced topological quantum phase transition between the symmetry-protected topological phase, and a state with magneto-orbital long-range order, for intermediate interactions V v /t. Our analytical treatment of the stronglyinteracting regime also points to the possible origin of inplane ferromagnetic order in the Creutz-Hubbard ladder with vertical tunnelings (1) instead of the imbalance, as recently found through a mean-field and numerical analysis [50]. Similar topological quantum phase transitions to magneticallyordered phases have also been found numerically in higherdimensional models [17].
Concerning the cold-atom realization, we note that the orbital spin operators (31) developing the long-range order (32) in the orbital ferromagnet, correspond to the standard spin operators for the two internal states of the atoms (78). Therefore, the magnetic correlations can be inferred from the spin structure factor, which can be measured for instance via Bragg scattering by playing with the polarization of the laser beams [61]. This can allow to test experimentally the validity of the critical line (33).

D. Intermediate interactions: Extended Hubbard models and quantum impurity physics
In this section, we elaborate on a theory that allows us to predict this interaction-induced topological quantum phase transition by making an interesting connection to the physics of quantum impurity models. We start by expressing the Hubbard part of the Hamiltonian (5) in the basis of the bulk (10) and edge (11) Aharonov-Bohm cages (AB-c). A similar step has been performed by rewriting the bosonic Creutz-Hubbard model in terms of the bulk negative-energy AB-c [46,47]. This corresponds to a projection onto the low-energy flat band, which is justified fort V v , and is also customary in the theory of fractional topological insulators [53]. In our case, however, it will be crucial to retain the positive-and zero-energy AB-c to capture the topological phases, and the possible quantum phase transitions that take place in the regimet ∼ V v .
In this basis, we find an extended Hubbard model with a nearest-neighbour interaction between fermions confined in adjacent AB-c where the AB-c number operators are n i,α = w † i,α w i,α for α ∈ {+, −}, and n η = η † η for η ∈ { , r}. This effective repulsion between the fermions can be depicted in a new ladder scheme, where the legs of the ladder correspond to the positive/negative energy flat bands, the sites correspond to the labels of the different bulk cages, except at the boundaries of the ladder, where the sites correspond to the edge cages (see Fig. 5). All the possible interactions, including the ones at the edges, are represented by curly lines. We also include the nearest-neighbour intra-and inter-leg tunnelling, which arises from interpreting the imbalance (15) within this virtual ladder.
In addition to this effective tunnelling and repulsion, we also find correlated tunnelling processes of two types. We obtain a pair-tunnelling Hamiltonian whereJ = −V v /4. The first term in the parenthesis describes how fermions confined in adjacent AB-c of opposite legs tunnel simultaneously to empty neighbouring cages, and can be understood as the anti-correlated pair-tunnelling of Fig. 5. The second term describes how fermions confined in adjacent Aharonov-Bohm cages on the same leg tunnel simultaneously to empty neighbouring cages, and can be understood as the correlated pair-tunnelling of Fig. 5. In addition, we also obtain an inter-leg density-dependent tunnelling where T d = V v /4 is the tunneling strength. The first line describes an inter-leg tunnelling within the bulk of the ladder that depends on the density difference of the neighbouring Aharonov-Bohm cages, and will be negligible for a groundstate with translationally-invariant bulk properties. On the other hand, the second line describes inter-leg tunnelings that occur at the boundaries of the Creutz-Hubbard ladder, and depend on the density of the edge AB-c (see Fig. 5). These terms are not negligible for translational-invariant bulks, and will play a key role in the topological phase transitions of the model. As shown in this section, the flat-band physics of the Creutz-Hubbard ladder provides an alternative route to the e + e e l = e r t imb Figure 5: Non-standard Hubbard terms in the Aharonov-Bohm-cage basis: Virtual ladder structure to represent the relevant processes in the Creutz-Hubbard ladder. The two legs of the ladder correspond to the two possible flat bands at energies ε ± = ±2t, with sites representing the bulk Aharonov-Bohm cages (10). The two boundary sites represent the two zero-energy modes ε l = ε r = 0, and their corresponding edge cages (11). The imbalance of the original Creutz-Hubbard ladder (15) leads to standard intra-and inter-leg tunnelling, but also to edge-bulk tunnelling processes represented by curved arrows and with strengths proportional to t imb . The Hubbard interaction (35) can be decomposed into nearest-neighbour interactions of strength proportional to V v that are represented by curly lines, and involve two adjacent bulk cages in any possible leg, or adjacent edge-bulk ones. Non-standard Hubbard terms involving a pair tunnelling also arise, where pairs of fermions tunnel simultaneously either from the same leg (correlated tunnelling) or from opposite legs (anti-correlated tunnelling) of the ladder with amplitudẽ J. Finally, we also describe a density-dependent inter-leg tunnelling of amplitude amplitude T d , where the fermion tunnelling depends on the density of fermions populating the Aharonov-Bohm cage at the neighbouring edge.
physics of non-standard Hubbard models, which typically arise in optical lattices as a consequence of dipolar interactions or higher orbitals [76]. In the present case, these nonstandard terms arise due to the interplay of interactions and the flatness of the bands, which yield with the terms introduced in Eqs. (12), (15), and (36)- (38). At first sight, this formulation leads to a Hamiltonian (39) that seems more complicated than the original one (3). However, as shown below, it is an ideal starting point to calculate an effective boundary theory for the topological edge states, based on which one can understand possible topological phase transitions as the imbalance and/or interactions are increased.

Bulk-mediated broadening and edge-to-edge tunnelling
As a warm-up to the discussion of the model for intermediate interactions, let us discuss how to derive an effective boundary theory in the non-interacting limit V v = 0. In Sec. III A, we have already discussed how the imbalance (15) may induce a quantum phase transition at ∆ε t | c = 4 that can be understood through a topological invariant (7) of the bulk bands. We now discuss how an effective boundary theory, derived via the above non-standard Hubbard Hamiltonian (39), yields an alternative description of such a phase transition.
The Hamiltonian (3) can be rearranged as follows stands for the Hamiltonian of the zero-energy edge modes. The bulk part of the Hamiltonian can be expressed as follows which is readily diagonalised by a Fourier transform considering periodic boundary conditions w N,α = w 1,α . This leads to the energy bands in Eq. (8), where the quasi-momentum lies in the Brillouin zone q ∈ (−π, π], and there is a total of 2N single-particle modes. To approximately solve the two-leg Creutz ladder with open boundary conditions, we build on the solution of an open single chain. We note that the bands in our ladder fulfill ε α ( π 2 + q) = ε α ( π 2 − q), which allows us to combine modes with momenta π 2 ± q to construct AB-c momentum operators via standing waves that respect the open boundary conditions of a finite ladder w N,α = w 0,α = 0. Note that the −q solutions of the problem with periodic boundary conditions are implicitly considered in the +q solutions. Together with the fact that the q = 0, π modes of the periodic solution yield trivial standing waves, this forces us to enlarge the number of sites in the problem with periodic boundary conditions N − 1 → 2N, such that the allowed momenta become q n = 2π N−1 n → q n = π N n. The transformation that approximately diagonalises the bulk Hamiltonian (41) is

Using these operators, the bulk Hamiltonian for the ladder with open ends approximately becomes
where B B B(q n ) = (0, ∆ε 2 cos q n , 2t + ∆ε 2 sin q n ), and we have introduced new spinors Ψ(q n ) = (w + (q n ), w − (q n )) t . Diagonalising this Hamiltonian leads to the curved energy bands (8), such that H bulk = ∑ n,α αε(q n )γ † α (q n )γ α (q n ) [77], where we have introduced γ + (q n ) = +u * q n w + (q n ) + v q n w − (q n ), γ − (q n ) = −v * q n w + (q n ) + u q n w − (q n ), together with the following constants and ζ (q n ) = 2t + ∆ε 2 sin q n , e iφ n = −i sgn(cos(q n )). We note that by using the periodic bulk energies (8), our approximation is essentially introducing an intensive contribution coming from the added bonds at the ends of the ladder.
Finally, the imbalance (15) also induces a bulk-edge coupling that can be understood as a hybridization between the edge and bulk orbitals To obtain the correct bulk-edge couplings g η n,α , note that the connection between the solutions of the periodic and open chains required enlarging the number of sites of the periodic chain. To preserve the distance between the edges, one must change the site indexes of operators that form the edge-bulk coupling (15), and leads to Interestingly, the Hamiltonian of the Creutz ladder in this formulation H edge + H bulk + H b−e is similar to a non-interacting two-impurity Fano-Anderson model [78], where the impurities correspond to the modes localised at the edges of the Creutz ladder, and the gapless metallic bands are substituted by gapped bulk bands also described by free fermions. This analogy yields an insightful interpretation of the topological quantum phase transition at ∆ε t | c = 4. Interpreting the dispersive bulk bands as reservoirs, the hybridisation can have two effects: (i) the zero-energy edge modes η = l, r may get shifted and broadened ε η = 0 → ∆ε η − iΓ η /2, (ii) the bulk fermions mediate an edge-to-edge tunnelling. Due to the particle-hole symmetry of the bands, the level shifts cancel ∆ε η = 0, such that the poles representing the edge excitations always correspond to zero energy. On the other hand, the broadening can be expressed as Γ η = ∑ α J η,α (ε η ), where we have introduced the spectral density for the coupling of the edge states to the bulk bands J ηα (ω) = 2π ∑ n |g η n,α | 2 δ (ω − αε(q n )).
Accordingly, the broadening of the levels depends on the value of the spectral function at zero energy Γ η = ∑ α J ηα (0). For ∆ε < 4t, the bulk bands (8) remain gapped, and thus J ηα (0) = 0, such that Γ l = Γ r = 0, and the edge modes thus remain in-gap bound states. As announced before, there will be a bulk-mediated tunnelling between these in-gap modes H e−e = t ee r † l + H.c., (49) where we have introduced the tunnelling strength t ee = ∑ n,α g r n,α g l n,α * αε(q n ) .
Provided that the bulk bands (8) remain gapped, this edgeedge tunnelling will decrease exponentially with the edge-toedge distance, such that the topological degeneracy of the groundstates (13) is preserved in the thermodynamic limit. Conversely, when ∆ε = 4t, the gap vanishes and the exponential localisation of the edge states disappears. Moreover, the bulk density of states at zero energy does not vanish anymore, such that J ηα (0) = 0, and the edge modes become broadened resonances Γ η > 0 instead of the previous bound states. It is very instructive to understand how the onset of a topological phase transition can be predicted by looking at the effective theory for the edges. This will become very useful in the presence of interactions, where a simple interpretation in terms of observables associated to a topological invariant (7) for non-interacting systems cannot be applied. Additionally, this result also underlies the importance of preserving all the orbitals in the effective description in terms of AB-c, and not simply projecting to the low-energy flat band when trying to describe flat-band effects in fermionic topological insulators.

Bulk-mediated dephasing and edge-edge interactions
Let us now derive an effective boundary theory in the balanced interacting model (∆ε = 0) for intermediate interactions V v /t. Therefore, we must analyse the following part H πCH = H FB + V nn + V pt + V dt of the Hamiltonian (39). Although the pair and density-dependent tunnellings in Eqs. (37)-(38) modify the distribution of particle-hole pairs in the rungs of the virtual ladder (see Fig. 5), the nearest-neighbour interactions (36) do not change, since the number of neighbouring AB-c that are occupied is preserved under such processes. Therefore, in this limit, the nearest-neighbour interactions (36) can be substituted by a c-number, and only the flat-band (12) and the correlated-tunnelling terms (37)- (38) have an important effect on the non-interacting groundstates (13). Moreover, as argued below Eq. (38), only the density-dependent tunnelling at the edges of the ladder will play a role to determine the order of the translationally-invariant groundstates.
According to this discussion, we can rearrange these terms as H πCH = H edge + H bulk + H b−e , where the edge Hamiltonian coincides with Eq. (40) for the non-interacting imbalanced case. In contrast, the bulk Hamiltonian is no longer described by a quadratic fermionic model, but instead by quartic terms that can be understood as some effective spin exchange. This becomes apparent after introducing the su(2) algebrã which should not be confused with the strong-coupling orbital spin operators of Eq. (31). With this notation, the bulk Hamiltonian becomes , which corresponds to a ferromagnetic Ising model in a transverse field. This model can be solved exactly for periodic boundary conditions [74] by means of a Jordan-Wigner transformation [66], namelỹ where f † i = ( f i ) † are spinless fermionic operators. Considering periodic boundary conditions f N = f 1 , we would obtain the energy bands for single-particle excitations where we have introduced the flatness parameterf = 8t/V v , such that we recover perfect flat bands in the non-interacting regime V v = 0. Here, the quasi-momentum is defined within a halved Brillouin zone q ∈ (0, π], such that there is a total of N single-particle modes. Let us remark that, just by looking at this effective theory of the bulk, one would predict a standard quantum phase transition at |J| = 2t onto an orbital ferromagnet that can be described by a local order parameter. However, this theory would predict a featureless disordered phase for |J| < 2t, and thus completely miss the topological features of the model. It is thus of paramount importance to include the boundary AB-c in the full description. To prepare for that, we use a similar approximation to treat the bulk as in the previous section. We note that the above energies have the symmetryε α (+q) =ε α (−q), such that we can move back to the entire Brillouin zone, and construct ABc momentum operators that respect the open boundary conditions f N = f 0 = 0 of a finite ladder by combining the ±q solutions of the periodic chain. In this case, in order to recover the N single-particle modes, and taking into account that the q = π mode yields a trivial sanding wave, we need to enlarge the number of sites in the problem with periodic boundary conditions N − 1 → N, such that the allowed momenta become q n = 2π N−1 n → q n = 2π N n with n ∈ {1, · · · , 1 2 N}. Accordingly, the mapping transforms the bulk Hamiltonian (52) approximately into whereB B B(q n ) = (0, −2J sin q n , 4t + 2J cos q n ), andΨ(q n ) = ( f (q n ), f † (−q n )) t are the Nambu spinors. Diagonalising this Hamiltonian leads to the bulk curved energy bands (54), such that H bulk = ∑ n,α αε(q n )γ † α (q n )γ α (q n ), where we have introducedγ together with the following constants andζ (q n ) = 2t +J cos q n , e iφ n = isgn(sin(q n )). In contrast to the single-particle bulk excitations in the imbalanced case (44), the fermionic operators (57) now describe Bogoliubov fermionic excitations [79]. We note again that by using the periodic bulk energies, our approximation is essentially introducing an intensive contribution coming from the added bonds at the ends of the ladder. Finally, the bulk-edge coupling corresponds to a spindensity interaction We can express this interaction in terms of the Bogoliubov excitations, after using the above Jordan-Wigner transformation, together with Eqs. (55) and (57). In a first approximation, we neglect the Jordan-Wigner string, which allows to derive a bulk-edge coupling that is analogous to the imbalanceinduced hybridization (46), namely (g l n,α l † l −g r n,α r † r )γ α (q n ) + H.c..
Once again, in order to obtain the correct bulk-edge couplings g η n,α , we have to consider that the connection between the modes of the periodic and open chain require enlarging the number of lattice sites, and thus changing the indexing of the operators in the edge-bulk coupling. This leads tõ Accordingly, the Creutz-Hubbard ladder in this formulation H πCH = H edge + H bulk + H b−e becomes again similar to a noninteracting two-impurity Fano-Anderson model [78]. In this case, the gapless metallic bands in the Fano-Anderson model are substituted by gapped bands of Bogoliubov fermions (57) due to the lack of conservation of the number of Jordan-Wigner fermions. The bulk-edge hybridisation (46) is substituted by a bulk-edge coupling (60) that does not conserve the number of Jordan-Wigner fermions either, and depends on the population of the edge modes. This analogy shall allow us to predict the onset of a topological quantum phase transition in the Creutz ladder caused solely by the Hubbard interactions.
Interpreting once more the dispersive bulk bands as reservoirs, the bulk-edge coupling can have two effects: (i) the zero-energy edge modes η = l, r may get shifted ε η = 0 → ε η = ∆ε η , and dephased with a rateΓ η /2, (ii) the bulk Bogoliubov excitations can mediate an edge-edge interaction (i.e. density-density interactions between the edge fermions mediated by spin-wave-type excitations). In analogy to the imbalance-induced broadening (48), the dephasing rate can be expressed asΓ η = ∑ αJηα (ε η ), where we have introduced the spectral density for the new bulk-edge coupling Accordingly, the dephasing of the superposition states depends on the value of this spectral function at zero energỹ Γ η = ∑ αJηα (0). For V v /4 < 2t, the bulk bands (54) remain gapped, and thus J ηα (0) = 0, such that there is no dephasing. In this regime, there are only level shifts and edge-edge interactions described by H e−e = ∆ε l n l + ∆ε r n r +U ee n l n r , where we have introduced Let us note that, in contrast to the imbalanced case, the particle-hole symmetry does not impose the vanishing of the level shifts. Due to the different bulk-edge coupling, these shifts depend on the filling of the Bogoliubov levels for the bulk groundstate |g bulk =γ † − (q 1 )γ † − (q 2 ) · · ·γ † − (q N−1 )|0 B , where |0 B is the Bogoliubov vacuum, and thus do not vanish.
If we now use a similar argument as for the imbalanceinduced quantum phase transition, we find that (i) the bulkmediated edge-edge interaction has no effect on the degenerate groundstate (13), since both edge modes are not populated simultaneously for a half-filled system. (ii) Although the zero-modes are shifted in energy, the topological degeneracy is preserved ∆ε l = ∆ε r , and no dephasing occurs provided that the Bogoliubov excitations remain gapped. Conversely, when |J| = 2t, or equivalently V v = 8t, dephasing takes place signalling that the edge modes are not well-defined singleparticle excitations. This argument thus locates the critical point of the interaction-induced quantum phase transition at V ṽ t | c = 8, which has been represented by a red circle in the phase diagram of Fig. 1.

E. Phase diagram of the Creutz-Hubbard ladder
Our considerations in the above sections already allowed us to determine the possible phases of the model and their phase boundaries in certain parameter regimes. In the following, we lay out the full phase-diagram of the model in the ( ∆ε t , V ṽ t ) plane using DMRG calculations. Our DMRG-code is based on matrix-product states and uses a two-site optimization scheme. It is built on the Abelian Symmetric Tensor Networks Library (developed in collaboration with the group of S. Montangero in Ulm) and capable of implementing multiple Abelian symmetries, such as particlenumber conservation. We consider lattices up to N = 256 sites and virtual bond-dimensions of up to m = 200.
The analytic and numerical results for the phase diagram, collected in Fig. 1, are described in the following subsections: 1. Topological insulator to orbital paramagnet phase transition As shown in Sec. III B, the mapping of the Creutz ladder onto a Quantum Ising ladder allows us to predict a critical line (29) separating the topological insulator (TI) and the orbital paramagnet (oPM) for sufficiently weak interactions. This critical line is represented by a red dashed line that starts form the point ∆ε = 4t,V v = 0 in Fig. 1.
The Ising transverse magnetization, or equivalently the density imbalance between the legs of the Creutz ladder (28), can serve as a good indicator for this quantum phase transition that can be easily calculated numerically using the DMRG code. We can determine the critical line by studying the divergence of the imbalance susceptibility χ ∆n = ∂ ∆n /∂ (∆ε/t) (Fig. 6(top)).
As an alternative means of identifying the TI phase, we can study the behaviour of the ground-state degeneracy in the Creutz-Hubbard model with variable filling. We therefore introduce the single-and two-particle energy gaps where E(x) is the ground-state energy of a system with x particles. It can be shown that the two quantities coincide for gapless systems (∆ = δ = 0) and conventional insulators (∆ = δ = 0). In a topological insulator, however, δ = 0 due to the presence of zero-energy edge modes while ∆ = 0 measures (half) the band-gap. In Fig. 6 (bottom), we show that the predicted behaviour is indeed observed. The critical points obtained through these different observables are represented by yellow stars in the left part of Fig. 1. As can be seen from these results, the analytical prediction of the phase boundary (29) is reasonably accurate even for quite large interactions, where the exact critical line given by DMRG departs from a straight line, and bends up.

Orbital ferromagnet to orbital paramagnet phase transition
In Sec. III C, we introduced an effective orbital Ising model in the limit of very strong interactions, which allowed us to predict a critical line (33) separating the orbital ferromagnet (oFM) and orbital paramagnet (oPM). This critical line is represented by a yellow dashed line in Fig. 1.
Indeed by measuring the paramagnetic and ferromagnetic magnetization ( T z i resp. T y i in Eq. (31)), we confirm that these quantities scale equally, and identify the phase-transition point also for finite interactions (Fig. 7). Technically, we determine the paramagnetic magnetization by measuring the fermionic observable that defines T z i , which is proportional to the leg density imbalance discussed above (see Eq. (31)). In order to avoid problems due to incomplete symmetry breaking when studying the ferromagnetic order-parameter T y i (i.e. between the possible alignments in the ferromagnetic phase), we determine instead the zero-momentum component of the orbital magnetic structure factor S T y T y (k) = 1 which yields the desired ferromagnetic magnetization in the thermodynamic limit T y i = (S T y T y (0)) 1/2 . We observe for both quantities an Ising-like scaling, which differs from the strong-coupling prediction only by a renormalization of the critical point and of the maximum ferromagnetic magnetization (comp. Fig. 7).
The critical points obtained through these magnetizations are represented by yellow stars in the right part of Fig. 1. As can be seen from these results, the analytical prediction of the phase boundary (29) is reasonably accurate even for moderate interactions.

Topological insulator to orbital ferromagnet phase transition
In Sec. III D, we derived an interesting connection between the balanced Creutz-Hubbard ladder for intermediate interactions and a Fano-Anderson-type model, which allowed us to predict the extension of the topological phase along the (0, V ṽ t ) axis of the phase diagram until a critical point V ṽ t = 8. Beyond this point, the long-range ordered orbital Ising magnet sets in, and the topological edge modes disappear into the bulk. This critical point is represented by a red circle in Fig. 1.
The numerical analysis (Fig. 8) confirms the validity of the effective Ising-model derived in Eq. (52), and the exact location of this critical point. Moreover, in the case of finite imbalance (∆ε = 0), the divergence of the paramagnetic susceptibility serves as a criterion for the determination of the phase-boundary (see inset Fig. 8). The critical points obtained by these means are represented by yellow stars in the middle part of Fig. 1.

Conformal field theories for the critical lines and entanglement spectrum for the phases
So far, we have used a conventional condensed-matter approach to explore numerically the phase diagram of the model, which is based on exploiting energy gaps, susceptibilities, and correlation functions to identify phases with long-range order or symmetry-protected topological phases, and critical lines that separate them. An alternative approach, based on the groundstate entanglement, has recently become a complementary method to understand the phase diagram of quantum many-body models [80]. For instance, the pair-wise concurrence 'susceptibility' can serve as a probe to localize quantum critical points [81], displaying a scaling behavior similar to that of certain observables in the more conventional condensed-matter approach. Other entanglement measures can also serve as probes of quantum criticality, as occurs for the block entanglement entropy S(l) = −Tr{ρ l log ρ l }, where ρ l = Tr L−l {|ε gs ε gs |} is the reduced density matrix of the left block with l sites for a bipartition of a chain of L sites. Remarkably enough, not only does the block entanglement entropy serve as a probe of criticality due to its divergence at a phase transition, but its scaling with the system size also reveals the central charge c of the conformal field theory (CFT) underlying the critical behavior [82,83]. For a critical system with open boundary conditions, the block entanglement entropy scales as where we have introduced a non-universal constant a. Since such entanglement entropy can be easily recovered from our MPS numerical results, calculating the central charge of the different critical lines of our phase diagram can serve as an additional confirmation of our previous derivations. In Sec. III B, we argued that the synthetic Creutz ladder for sufficiently weak interactions can be understood as a couple of Ising models of length L = N with a renormalized transverse field. Accordingly, the corresponding CFT should have central charge of c = 1/2 + 1/2 = 1, such that we would expect the scaling S(l) = 1 6 ln 2N π sin πl N + a. This is a natural connection to the non-interacting regime of Sec. III, where we argued that the phase transition can be understood in terms of a Dirac fermion with a Wilson mass in lattice gauge theories. This Wilson fermion becomes massless at the critical point, and corresponds to the CFT of a single massless Dirac fermion with central charge c = 1.
For the strongly-interacting regime of Sec. III C, we showed that the oFM-oPM quantum phase transition can be predicted in terms of a single Ising model of length L = N in a transverse field. Accordingly, the corresponding CFT should have central charge of c = 1/2, and S(l) = 1 12 ln 2N π sin πl N +ã. In contrast to the previous case, the critical phenomena is governed by the CFT of a single Majorana fermion with central charge c = 1/2.
Finally, in the intermediate interacting regime of Sec. III D, we argued that the relevant physics to understand the TI-oFM phase transition is by approximating a complicated nonstandard Hubbard model with a simpler quantum impurity model for the edges coupled to a bulk of spins described by yet another Ising model of L = N − 1 sites in a transverse field. Technically, however, we determine the entanglement entropy between the N physical sites in the original basis and therefore measure S(l) = 1 12 ln 2N π sin πl N + a for a system of length L = N, corresponding to the CFT of a single Majorana fermion with a central charge of c = 1/2.
We confirm the above predictions through the numerical determination of the central charge along the critical lines in three representative cases (see Fig. 9). We find central charge values agreeing with c = 1/2 for the TI-oFM and the oFM-oPM-transition. The charge c = 1 along the TI-oPM-transition originates from the hybrid nature of the Ising-model describing it (see Eq. (20)). Building on these results, we depict the In the previous sections, we provided several indicators of the non-trivial topological nature of the TI phase, such as the resilience of the zero-energy edge modes as interactions are switched on (see Sec. III D), or the different single-and two-particle gaps (65) that distinguish topological and nontopological phases (see Sec. III E). Let us now provide further evidence by using entanglement properties of the groundstate. In particular, a strong signature of the presence of topological order can be extracted from the study of the entanglement spectrum (ES) [8]. Similarly to the case of the entanglement entropy, to define the ES we take a bipartition of the ladder and the resulting Schmidt decomposition of the ground state |ψ = ∑ i λ i |ψ i L ⊗ |ψ i L−l , where |ψ i L and |ψ i L−l are basis vectors of the two parts, satisfying the orthogonality condition ψ j l |ψ i L−l = δ i j , whereas the λ i are the Schmidt values. By convention, the ES is defined as a logarithmic rescaling of the Schmidt values, −2 log(λ i ) and it can be again extracted straightforwardly from MPS calculations. As originally pointed out in [51] in the context of the characterization Figure 10: Degeneracies of the entanglement spectrum for different phases. For a ladder of length L = 128 and for a bipartition in the half chain, the twenty lower eigenvalues of the ES are depicted for the three different phases. The dots represent the degeneracy of the corresponding eigenvalue. In the TI phase, for V v /t = 4 and ∆ ε /t = 1.5 the eigenvalues are all doubly degenerate. In the oPM phase, for V v /t = 4 and ∆ ε /t = 2, and in the oFM phase for V v /t = 9 and ∆ ε /t = 0.2 almost all the eigenvalues are not degenerate. We found the same behaviour elsewhere in phase space.
of the Haldane phase of Heisenberg-type magnets, the degeneracy of the ES robustly identifies the symmetries protecting the topological phase. As shown in Fig. 10 in the imbalanced Creutz-Hubbard model, the ES in the TI phase clearly shows doubly-degenerate eigenvalues, whereas in the oPM and oFM phases the ES is trivial and almost completely non-degenerate. This supports the topological nature of the wide region of the phase diagram labelled as TI (see Fig. 1), and demonstrates that the topological insulating phase survives to considerably strong interactions.

IV. SYNTHETIC CREUTZ-HUBBARD LADDER
In this section, we discuss an experimental setup capable of implementing the imbalanced Creutz-Hubbard Hamiltonian (3). We shall focus on a particular set of microscopic couplings that can be realised with ultracold atoms in optical lattices by exploiting the tools of laser-assisted tunnelling.

A. Cold-atom Creutz-Hubbard ladder
We consider a cubic state-independent optical lattice that traps a two-component atomic gas with hyperfine states |↑ = |F = I − 1 2 , M , |↓ = |F = I + 1 2 , M (see Fig. 11(a)). In the Wannier basis, the Hamiltonian corresponds to the standard Hubbard model [26], namely where the fermionic operators f † i,σ , f i,σ create and annihilate an atom in the electronic state σ at site i of the lattice, and n i,σ = f † i,σ f i,σ are the number operators. In this equation, t is the tunnelling strength of atoms between neighbouring potential wells along the x-axis, ε σ stand for the energies of the electronic levels, and U ↑↓ corresponds to the on-site interaction strength due to s-wave scattering. As customary [23], tunnelings and interactions of a longer range, as well as tunnelings along the y and z axes, have been neglected in this expression. This is justified for sufficiently deep optical lattices V 0,y ,V 0,z V 0,x E R , where V 0α are the corresponding amplitudes of the optical potential V (r r r) = ∑ α=x,y,z V 0α sin 2 (kr α ), E R = k 2 /2m is the recoil energy, and k is the wavevector of the retro-reflected laser beams forming the optical lattice. In this regime, the parameters of Eq. (69) can be expressed as where a ↑↓ is the s-wave scattering length for the collision of two atoms in the two different hyperfine states. The first step towards a possible implementation of the Creutz-Hubbard model based on Eq. (69) is to represent the legs of the ladder by the two hyperfine states [84,85] (see Fig. 11(b)). This might be interpreted as the smallest-possible synthetic dimension along which real [86] or complex [87] tunnelings can be implemented via Raman transitions, leading to recent experiments realising synthetic gauge fields [56,58]. Unfortunately, the Creutz ladder involves more complicated tunnelings (see Fig. 11(c)) that cannot be directly obtained using this scheme. One possibility to implement the required tunnelings could be to combine bi-chromatic optical lattices with additional Raman transitions and a staggered optical potential, which allows for a very flexible toolbox to realise spindependent tunnelings [85]. However, in view of the success of experimental schemes based on periodic modulation of the lattice [33,34] that can be understood in terms of Floquet engineering [88], we hereby present an alternative scheme to achieve the desired Hamiltonian combining periodic modulations with a generalisation of the Raman-assisted scheme [28] applied to a spin-independent optical lattice.

State-independent energy gradient
We consider a linear energy gradient that tilts the optical lattice independently of the hyperfine state. Such a linear gradient can be obtained by accelerating the lattice, or by exploiting the ac-Stark effect, and yields Note that we have considered the regime ∆ V 0x , such that the gradient does not modify the inter-site terms of the original Hubbard model (69), but only leads to the local on-site term (71) in the Wannier basis. Moreover, we impose t ∆, such that the tilt inhibits the original intra-leg tunnelling in Fig. 11(b). The goal now is to re-activate the tunneling against this gradient, by generalising the ideas of schemes based on Raman transitions [28], or on periodic modulations [89] with additional shallower optical lattices [33].

Raman-assisted tunnelling
We consider three additional laser beams. These lead to a couple of Raman two-photon excitations that assist the inter-leg tunnelings (see Fig. 12(a)). A first pair of laser beams with frequencies ω 1 , ω 2 assists the inter-leg tunnelings in Fig. 12(b) by virtually populating an excited state. The effect of the lasers is two-fold, they provide the energy to overcome the potential offset in the tunnelling process ω 1 − ω 2 = (ε ↓ − ε ↑ ) + ∆, and they also induce a recoil kick δ k = (k 1 − k 2 ) · e x along the x-axis, thus allowing the tunnelling to occur (i.e. overlap of neighbouring Wannier functions). The other pair of laser beams assists the inter-leg tunnelling in Fig. 12(c) provided that ω 1 − ω 3 = (ε ↑ − ε ↓ ) + ∆, also yielding a recoil kick δk = (k 1 − k 3 ) · e x . Altogether, these Raman-assisted terms lead to where we have introduced the Raman-assisted tunnelling strengths which are expressed in terms of the respective two-photon Rabi frequencies Ω 12 , Ω 13 corresponding to the particular Raman transition, and an exponential term due to the laserassisted overlap of the neighbouring Wannier orbitals, which has been calculated using a Gaussian approximation. We shall consider the limit where the Raman-assisted tunneling strengths are set to be equal Ω =Ω. To arrive at these expressions, we have considered again the regime of deep optical lattices, and that δ k · λ = 2πn, δk · λ = 2πñ, where n,ñ ∈ Z and λ is the wavelength of the retro-reflected laser beam that leads to the original optical lattice. This can be achieved by controlling the angle of the Raman lasers with respect to the xaxis. More importantly, we have assumed that the two-photon Rabi frequencies fulfil |Ω 12 |, |Ω 13 | ∆, such that on-site Raman transitions are highly off-resonant, and can be neglected. Let us note that, by avoiding a state-dependent optical lattice that traps each of the two components differently, our scheme is not subjected to the heating problems that become specially troublesome for fermionic atoms [90].

Intensity-modulated optical lattice
To obtain the complex horizontal hopping in Fig. 11(c), we must re-activate the intra-leg tunnelling inhibited by the lattice tilting (71), and dress it with the correct Peierls phases. This can be achieved by introducing a periodic driving provided by e " e # e " e # w 1 w 2 w 3 w 1 w 1 w 3 = (e " e # ) + D w 1 w 2 = (e " e # ) + D Figure 12: Raman-assisted tunnelling scheme for a binary hyperfine mixture of ultracold atoms in an optical lattice: (a) Atoms in a given hyperfine state from the groundstate manifold nS 1/2 can tunnel to the neighbouring potential well changing its hyperfine state by means of a Raman transition that virtually populates an excited state from the manifolds nP 1/2 , nP 3/2 . (b) The inter-leg tunnelling of |↑ (blue circles) to |↓ (red circles) is mediated by a pair of lasers providing the energy necessary to overcome the offset ω 1 − ω 2 = (ε ↓ − ε ↑ ) + ∆, and a recoil kick to assist the tunnelling. (c) The inter-leg tunnelling of |↓ (red circles) to |↑ (blue circles) is mediated by a pair of lasers providing the energy necessary to overcome the offset ω 1 − ω 3 = (ε ↑ − ε ↓ ) + ∆, and a recoil kick to assist the tunnelling. an additional optical lattice with a modulated intensity (see Fig. 13(a)), namely where the intensity of the ac-Stark shifts for each of the hyperfine states is modulated periodically in time, V d,σ (t) = V d,0 sin(ω d,σ t − φ σ ) with frequency ω d,σ and a phase φ σ . In analogy with the lattice tilting (71), we have assumed that the amplitude fulfills V d,0 V 0x , such that we can neglect any periodic modulation on the bare tunnelling of the original Hubbard model (70), and consider instead a periodic driving of the on-site energies (74). More importantly, since this lattice must be very shallow, one can minimise the scattering of photons from the excited state in the calculation of the two-photon ac-Stark shift, which can lead to the aforementioned problematic heating effects in the opposite regime of deep state-dependent optical lattices [90]. As an interesting alternative, we note that a time-dependent magnetic-field gradient [91] can also be exploited to obtain state-dependent time-periodic modulations.

Effective Creutz-Hubbard Hamiltonian
Once the ingredients have been introduced, let us show how the synthetic Creutz-Hubbard model is obtained. If the wave- t h e i2f # Figure 13: Photon-assisted tunnelling scheme for a binary hyperfine mixture of ultracold atoms in an optical lattice: (a) Atoms in a given hyperfine state can tunnel to the neighbouring potential well preserving its hyperfine state by absorbing energy from an additional optical lattice whose intensity is periodically modulated. (b) The intra-leg tunnelling of |↑ (blue circles) is activated by the weak modulated lattice with ω d,↑ = ∆/2, such that the tunnelling picks the modulation phase e −i2φ ↑ . (c) The intra-leg tunnelling of |↓ (red circles) is activated by the weak modulated lattice with ω d,↓ = ∆/2, such that the tunnelling picks the modulation phase e −i2φ ↓ .
length of this intensity-modulated lattice is twice that of the original lattice λ d = 2λ , only the even sites will be subjected to the periodic modulation (see Fig. 13(a)). Moreover, if the resonance condition is met, namely nω d,σ = ∆ for n ∈ Z, the nearest-neighbour tunnelling can be restored by absorbing energy quanta from the periodic drive. By setting ω d,σ = ∆/2, the intra-leg tunnelings acquire the phase of the corresponding drivings (see Fig. 13(b)-(c)), and lead to one can explore the properties of the imbalanced Creutz-Hubbard Hamiltonian in Eq. (3) by applying a laser-assistedtunnelling scheme to a two-component gas of fermionic atoms loaded 1D optical lattice, such that the components play the role of the two legs of the ladder (78). The Hamiltonian parameters can be controlled experimentally using the mappings in Eqs. (79), (80), and (82), together with the expressions (70) and (73), as a guiding principle. The particular Hamiltonian (3) witht = tJ 0 (V d,0 /∆), ∆ε = δ /2, and V v = U ↑↓ , is obtained by tuning the ratio of the bare and Raman tunnelings t/Ω = −2J 0 (V d,0 /∆)/J 2 (V d,0 /∆), such that the strength of the inter-and intra-leg tunnelings is equal, and by fixing the phases such that φ ↑ = π/4 = −φ ↓ leads to a net π-flux.

V. CONCLUSIONS AND OUTLOOK
In this work, we have advanced on the understanding of the competition between topological features and interaction effects in quantum many-body systems. By focussing on the paradigmatic Creutz-Hubbard ladder, we have developed a variety of analytical tools, and performed a thorough numerical analysis, unveiling characteristic features and general methods that could be applied to other strongly-correlated topological insulators/superconductors. Moreover, our predictions can be readily tested in table-top experiments with ultracold fermionic atoms in optical lattices, which may serve as quantum simulators of interacting topological phases.
By using a pair of electronic states of the atoms as synthetic legs of the ladder, and by exploiting laser-assisted tunnelling to control their dynamics, we have shown that the physics of a variant of the Creutz ladder can be implemented in a onedimensional optical lattice. Such a variant is obtained by applying a Zeeman splitting between the two legs of the ladder, which puts the system into the AIII class of topological insulators. To the best of our knowledge, topological insulators in such a symmetry class have not been realized in coldatom experiments yet. Moreover, by adding on-site repulsive interactions only, we can induce phase transitions of different universality classes into an orbital ferro-or paramagnetic phase. Surprisingly, we found that the expected Dirac CFT transition line at weak interactions splits into two Majorana CFT critical lines, once the Hubbard term dominates over the Zeeman term: these findings are also numerically supported by the scaling of entanglement entropy. In addition, we provided an understanding of these topological phase transitions in the language of quantum impurity physics, shedding new light on the hybridisation mechanism of the edge states. Finally, we have identified experimentally accessible signatures to test our theoretical predictions.
From the above results, we are convinced that the synthetic Creutz-Hubbard ladder can become a workhorse in the theoretical and experimental study of correlated topological insulators, as it allows for a very neat understanding of the underlying phase diagram, and a clear path to implement the model with ultracold atoms. This could be helpful in order to gain even deeper insight into the tri-critical point, serving as a guide to construct an effective quantum field theory that describes the mechanism underlying the splitting of the central charge into the two critical lines. It would be very interesting to study the imbalanced Creutz-Hubbard model at different fillings, and to explore the possibility of finding topological phases of matter that disappear for vanishing interactions. In this respect, the analytic and numerical methods hereby presented may be generalised to other fillings, allowing to go beyond mean-field arguments that support the existence of such interesting ground states.