Methods for detecting charge fractionalization and winding numbers in an interacting fermionic ladder

The experimental engineering of topological phases of matter in ultracold atomic gases [1–6] lays the foundations for a deeper understanding of phase transitions that transcend the Landau paradigm of symmetry breaking. In these experiments, models displaying non-local forms of order are realized in highly controllable environments, where the parameters driving the system in and out the topological phases can be tuned with wide freedom, and where observables complementary to those of a typical solid-state experiment can be measured. These experiments have already revealed interesting properties associated with two-dimensional (2D) topological Bloch bands: the anomalous (Hall-like) velocity, which was detected in response to an external force [3, 4], the topologically invariant Chern number [4], and chiral edge currents [2, 5, 6].

Following these advances, an important objective would be to probe the edge modes of one-dimensional (1D) topological systems, which typically appear at zero energy and exhibit charge fractionalization. In particular, identifying an observable that unambiguously signals their presence in experiments constitutes a remarkable challenge. Detecting the properties of zero-energy edge modes would strongly complement the Zak-phase measurement, recently demonstrated with bosonic atoms in a 1D optical superlattice [1].

Several theoretical efforts have been devoted to the design of realistic platforms hosting topological superconducting phases with ultracold fermions [7–14], including number-conserving setups [13, 15–20]. Here, we envision an even simpler scenario, based on the fact that similar topological edge physics can be accessed without pairing mechanisms. Indeed, it is a generally overlooked fact that pairing interactions are not strictly necessary to mimic topological superconductors (TSCs). A fundamental example is offered by the Su–Schrieffer–Heeger (SSH) model [21], which presents particle–hole symmetry and belongs to a non-trivial topological class of chiral Hamiltonians, namely the class BDI of the Altland–Zirnbauer classification [22–26]. Similar 1D fermionic models without superconducting interactions display topologically protected edge modes localized at their boundary, which are Dirac-like [27] and feature remarkable properties, such as charge fractionalization [28].

The goal of this article is twofold. First, we propose a route to mimic the physics of TSCs using state-of-the-art ultracold fermionic experiments. We exploit a simple idea: employing a two-leg ladder to double the fermionic species, in such a way that the legs are respectively associated with effective holes and particles [29, 30], see figure 1. This geometry is indeed well-suited for experiments on ultracold gases, as it has already been realized for atoms trapped in optical lattices, either in physical ladder geometries [2], or exploiting internal degrees of freedom as an artificial dimension [5, 6, 31–33]. We characterize the topological properties of this model, both in the absence and presence of interactions, finding several affinities with previous studies of 1D topological interacting fermionic systems [33–40]. We argue that such systems constitute a useful playground, not only to examine the appearance of symmetry-protected edge modes, but also to study the role of contact interactions, which may be tuned to drive transitions between trivial and topological phases. This analysis of the Hubbard repulsion extends to spin 1/2 fermions previous studies about the effect of interaction in generalizations of the SSH model [41, 42].

Zoom In Zoom Out Reset image size

As a second objective, we analyze in detail how topological signatures might be directly observed in such interacting systems. Beside the detection of fractionalized edge-modes, we also focus on the winding number associated with the expectation value of the spin, which provides a good detection tool for topological phases also in the presence of interactions. In this way, we extend to a 1D model in the topological class BDI the techniques developed to reveal the topology of cold-atom realizations of the 2D Haldane model [43, 44] and other 2D topological systems [45–48]. For both these observables, we examine the effect of a trapping potential, which sets soft boundaries to the system (usually considered to alter the observation of edge physics). Our schemes are based on the direct observation of the atomic cloud or on time-of-flight measurements: in both cases they are extremely robust to such confinement.

This article is organized as follows. In section 2 we introduce the model and provide an intuitive description of its symmetries. In section 3 we focus on the non-interacting model at half-filling and thoroughly characterize the topological insulator (TI) that is reached for a certain range of parameters. Furthermore, methods to detect unambiguous signatures of the topological properties are proposed, based on the density profile and the momentum-distribution of the gas. In section 4 we study the role of interactions and characterize the related interacting TI. In section 5 we describe a possible physical realization of the model based on laser-assisted tunnelings and in section 6 we present our conclusions. Finally appendix presents a detailed analysis of the non-interacting model, its order parameter and spin winding number.

We consider a spinful fermionic ladder in the presence of external gauge potentials as depicted in figure 1. Here the two legs of the ladder are associated with a pseudo-spin ${\tau }_{z}$ and the lattice sites along the main axis (x direction) are labelled by $r=1,\ldots ,L.$ We introduce the four-component fermionic operator ${\hat{a}}_{r},$ defined on the lattice site r, which acts on both the pseudo-spin ${\tau }_{z}$ and the spin related to the two internal states of the fermion σ_z (two commuting sets of Pauli matrices ${\tau }_{\alpha }$ and σ_α are used to describe these degrees of freedom). The Hamiltonian describing the ladder system is taken to be of the form (${\hslash }=1$)

$$\begin{eqnarray}&&{\hat{H}}_{0}=\displaystyle \sum _{r}\left\{t\left({\hat{a}}_{r}^{\dagger }{\tau }_{z}{{\rm{e}}}^{{\rm{i}}\frac{B}{2}{\sigma }_{z}}{\hat{a}}_{r+1}+{\rm{h.c.}}\right)+{\hat{a}}_{r}^{\dagger }\left({\rm{\Omega }}{\sigma }_{x}+J{\tau }_{x}+\mu {\tau }_{z}+{\mu }_{0}\right){\hat{a}}_{r}\right\},\end{eqnarray} \tag{ 1 }$$

and is schematically represented in figure 1. The first term in equation (1) describes the intra-chain tunneling along the x direction, with hopping amplitude t and spin-dependent Peierls phase-factor $\mathrm{exp}({\rm{i}}B{\sigma }_{z}/2),$ which represents a 'spin–orbit coupling' analogous to those already realized in fermionic [5] and bosonic [6] chains. Note that, because of the ${\tau }_{z}$ factor, the two chains have opposite kinetic energy; consequently, the motion around each plaquette of the ladder acquires a π-phase independently of the σ_z-spin component. The second term describes an on-site spin-flip term with amplitude Ω, the inter-chain tunneling with amplitude J, the potential difference between the two chains μ, and the overall chemical potential μ₀.

In the following we will elaborate on the fact that this system effectively reproduces some of the physical features of TSCs without recurring to any physical pairing mechanism. This analogy is based on the idea that atoms in the first chain (${\tau }_{z}=\;+1$) can be identified with conduction electrons of a generic 1D superconducting model, ${\hat{a}}_{r;{\tau }_{z}=+1}^{\dagger }\equiv {\hat{c}}_{\mathrm{elect}}^{\dagger }(r),$ whereas those in the second chain (${\tau }_{z}=-1$) can be identified with its holes, ${\hat{a}}_{r;{\tau }_{z}=-1}^{\dagger }\equiv {\hat{c}}_{\mathrm{hole}}^{\dagger }(r).$ In this picture any tunneling from one chain to the other constitutes an effective pairing interaction, i.e. ${\hat{a}}_{r;{\tau }_{z}=+1}^{\dagger }{\hat{a}}_{r;{\tau }_{z}=-1}\equiv {\hat{c}}_{\mathrm{elect}}^{\dagger }{\hat{c}}_{\mathrm{hole}}\ \approx \ {\hat{c}}_{\mathrm{elect}}^{\dagger }{\hat{c}}_{\mathrm{elect}}^{\dagger },$ where the last equality is justified by the Bogoliubov–De Gennes treatment of the superconductor. Within this parallelism we interpret the four-band Hamiltonian (1) as a Bogoliubov–de Gennes Hamiltonian in the superconducting picture. Specifically, the particle–hole symmetry, which plays a key role in the physics of TSCs, is here represented by a swap of the two chains, $C={\tau }_{y}{\sigma }_{y},$ which have opposite kinetic energy in the same way as particles and holes do. Such mapping, though, must be seen only as an analogy, since the number of degrees of freedoms in the system (1) is doubled with respect to the superconducting wire and this has important physical consequences, as will be discussed in the following. Finally, note that Hamiltonian (1) is unitarily related to those considered in the four-wire setup of [29, 30] and may have a relevance also for the study of electronic gases.

The model in equation (1) can be realized using cold atoms trapped in an optical lattice. We present here an overview of the experimental proposal and refer the interested reader to section 5 for a detailed analysis of the implementation of the model. Let's start considering a 2D setup. The realization of a spin-dependent intra-chain tunneling, described by the first term in equation (1), is particularly challenging, as it requires a subtle control over the hopping amplitudes. This effect could be engineered by exploiting the laser-induced-tunneling methods implemented in recent experiments [2, 4, 49, 50]. Specifically, we propose to achieve this task by combining a spin-dependent staggered potential with large energy offset Δ between neighboring sites, inhibiting the bare hopping along the x direction, together with an onsite energy modulation set at the resonant frequency $\omega ={\rm{\Delta }}.$ The spin-dependent staggered potential is chosen to be opposite for the two internal states, i.e. ${V}_{\mathrm{stag}}(r)={(-1)}^{r}({\rm{\Delta }}/2){\sigma }_{z},$ which can be realized by considering an appropriate anti-magic wave-length [51]; this choice is motivated by the fact that the resonant modulation will then generate effective tunneling matrix elements of the desired form ${t}_{\mathrm{eff}}(r)=t\mathrm{exp}[{\rm{i}}\phi (r){\sigma }_{z}],$ see [49, 52, 53]. In order to make the Peierls phase-factors constant over the whole lattice, i.e. $\phi (r)=B/2,$ we propose to modulate the lattice with two pairs of lasers; such a configuration allows to address individual links independently [4], hence realizing the desired Peierls phase factors on all links, (see section 5). Using additional fields resonant with the energy difference between the two spin-states a tunable onsite spin-flip term Ωσ_x can be realized. Finally, the two-leg ladder can be isolated using an additional superlattice, or a light-intensity mask [54, 55].

Hamiltonian (1) is characterized by four energy bands; as shown in figure 2, by varying the filling of the ladder, and thus the chemical potential μ₀, the system is driven through a series of metallic and insulating quantum phases (we consider in this article only the case of zero temperature). For half filling, corresponding to the case where μ₀ = 0 and the particle density is $\rho \equiv N/L=2$ (N is the number of fermions), the single-particle Hamiltonian shows both the particle–hole symmetry we sought for, defined by the operator $C={\tau }_{y}{\sigma }_{y}$ and an additional time-reversal symmetry, $T={\sigma }_{x},$ which bring the system into the topological class BDI (see appendix for more details). This class includes, for example, the SSH model and, according to the periodic table of TIs and superconductors [23, 24], it may present topological phases with zero-energy modes. Specifically, our model displays a non-trivial topological insulating phase for ${{\rm{\Omega }}}_{c,1}\lt {\rm{\Omega }}\lt {{\rm{\Omega }}}_{c,2},$ where ${{\rm{\Omega }}}_{c,i}$ are defined, for $B\lt \pi ,$ as:

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{c,1}^{2}\equiv {J}^{2}+{\left(\mu -2\mathrm{cos}(B/2)\right)}^{2};\qquad {{\rm{\Omega }}}_{c,2}^{2}\equiv {J}^{2}+{\left(\mu +2\mathrm{cos}(B/2)\right)}^{2}.\end{eqnarray} \tag{ 2 }$$

We find that the topological phase is surrounded by two topologically trivial phases. For μ₀ = 0, trivial and topological phases are distinguished by a topological order parameter ${\mathcal{W}}$ which takes the respective values of $+1$ and $-1$ (see the appendix and in particular equation (A.7) for the definition which is based on the technique developed in [56]).

Zoom In Zoom Out Reset image size

3.1. Fractionalized edge modes

In the topological phase, two zero-energy fermionic modes appear in ladders with open boundary conditions, as showed in figure 2. These modes are exponentially localized at the ends of the system and have important consequences on the density distribution characterizing the topological insulating phase when N = 2L + 1 fermions are introduced in the system with hard-wall boundary conditions, as displayed in figure 3. We observe that such modes are described by Dirac operators and they are not Majorana modes as it would be expected in the superconducting analog wire. The figure shows that the density in the bulk of the system indeed corresponds to the expected value ρ = 2. Moreover, analogously to the SSH model, a charge 1/2 is exponentially localized at each boundary. This important signature of charge fractionalization can be suitably identified through the expectation value of the operator ${\hat{n}}_{j}^{*}={\displaystyle \sum }_{m=1}^{j}({\hat{n}}_{m}-2)$ where ${\hat{n}}_{m}={\hat{a}}_{m}^{\dagger }{\hat{a}}_{m},$ which detects the excess density with respect to the bulk value ρ = 2. As illustrated in the insets of figure 3, an overall excess density of 1/2 is localized within a few sites from the left and right edges of the sample.

Zoom In Zoom Out Reset image size

Such a signature can be observed even in the presence of a harmonic confinement, described by the following contribution to the Hamiltonian: ${\hat{H}}_{{\rm{tr}}}={\displaystyle \sum }_{r}{w}_{r}{\hat{a}}_{r}^{\dagger }{\hat{a}}_{r},$ with ${w}_{r}=\bar{w}{(r-L/2)}^{2}.$ The effect of an external potential can be understood in a Thomas–Fermi approach as a space-dependent chemical potential ${\mu }_{0}(r).$ Due to the four energy bands, the system has three insulating phases for intermediate fillings, and, in the presence of a harmonic trap, this yields a typical wedding-cake structure with integer density plateaus (see figure 4, first column). Remarkably, even in the presence of the trapping it is possible to identify the fractionalized modes, as we see next. These zero-energy modes extend in the intermediate metallic region between one trivial plateau $\left(\rho =1,3\right)$ and the topological one $\left(\rho =2\right),$ up to exponential corrections.

Zoom In Zoom Out Reset image size

In figure 4, second column, we show the expectation value of ${\hat{n}}_{j}^{**}={\displaystyle \sum }_{m=0}^{j}({\hat{n}}_{L/2+m}-1),$ which is particularly suited for the detection of fractionalized edge modes in cases where the density in the center of the trap is ρ = 2. Moving from the center of the trap to the next plateau (ρ = 1), this operator measures the excess density with respect to the ρ = 1 value. One can obtain either an integer (no fractional modes) or an half-integer value (presence of one fractional mode). This is an unambiguous signature of the non-trivial topological phase (third row): indeed, in this case, the quantity $\langle {\hat{n}}_{j}^{**}\rangle$ becomes half-integer for values of j corresponding to the distance of the ρ = 1 plateau from the center of the trap.

As a final remark, let us stress that the problem of detecting fractionalized edge modes through a density measurement was first addressed in [57, 58], where this detection method relies on the optical measurement of reflected light. The recent experimental advances, however, allow for the challenging method presented above, since the feasibility of a single-atom detection for ultracold fermions in optical lattices has indeed been demonstrated [59–61]. In particular, a combination of laser cooling and fluorescence detection enables an unambiguous measurement of the occupancy of single sites for both ⁴⁰K [59, 60] and ⁶Li [61] gases. This is of particular importance for the detection scheme that we are proposing, because it could suffer from the experimental inability to fix the total number of atoms which are used in the many experimental realizations necessary to reconstruct the signal $\langle {\hat{n}}_{j}\rangle .$ The novel single-atom microscopes will also allow a post-selection based on the global number of particles of the system, necessary to obtain an accurate measurement.

3.2. Spin winding number

Another interesting signature of the topological phase is offered by the expectation value of the spin operator:

$$\begin{eqnarray}&&\hat{{\boldsymbol{\Sigma }}}(k)=\displaystyle \frac{1}{2}{\hat{a}}_{k}^{\dagger }\left({\boldsymbol{\sigma }}\otimes {\tau }_{0}\right){\hat{a}}_{k},\end{eqnarray} \tag{ 3 }$$

where ${\hat{a}}_{k}$ is the four-component annihilation operator in momentum space. This approach is inspired by the techniques presented for 2D systems in [43–47] and for ladders in [62], where it was shown that the expectation value of the spin $\langle \hat{{\boldsymbol{\Sigma }}}(k)\rangle$ provides a good observable to identify the topological invariant (winding number) of certain TIs. In the following, we generalize this procedure to our quasi-1D ladder model and show that also here the topological invariant ${\mathcal{W}},$ which clearly identifies the non-trivial topological regime, can be extracted from $\langle \hat{{\boldsymbol{\Sigma }}}(k)\rangle .$ This is thus another example of the interesting concept that a time-of-flight measurement can detect topological order.

The Hamiltonian (1) can be written in a real form thanks to its symmetries; therefore $\langle {\hat{{\rm{\Sigma }}}}_{y}(k)\rangle =0$ for each eigenstate of the system and $\langle \hat{{\boldsymbol{\Sigma }}}(k)\rangle$ always lies in the $\hat{x}-\hat{z}$ plane. The 'spin winding number' ${\mathcal{S}}$ is defined as the number of times the vector $\langle \hat{{\boldsymbol{\Sigma }}}(k)\rangle$ encircles the origin for k going from 0 to $2\pi$ (the lattice spacing is set to 1). Let us denote with ${{\boldsymbol{S}}}_{2}(k)/2$ the expectation value of the spin operator (3) for a state of non-interacting fermions filling completely and solely the second energy band. Remarkably, the parity of ${{\mathcal{S}}}_{2},$ which is the winding number of ${{\boldsymbol{S}}}_{2}(k),$ coincides with the topological order parameter of the model:

$$\begin{eqnarray}&&{\mathcal{W}}={(-1)}^{{{\mathcal{S}}}_{2}}\end{eqnarray} \tag{ 4 }$$

(see appendix for a demonstration and for details on the analytical calculation of this topological index for this specific model). Equation (4) is analogous to those derived for several other 2D models [43–47]: it relates a topological invariant to a quantity, ${{\mathcal{S}}}_{2},$ to be extracted via time-of-flight imaging. A similar behavior was discussed in [62] for a two-band generalizations of the SSH model. In our case, we stress that the second band of the model is the lowest-energy band with non-trivial topological order; this is also related to the fact that the fractionalized edge modes appear in the second bulk gap.

We now describe how to measure ${{\mathcal{S}}}_{2},$ for realistic systems, even in the presence of a harmonic trap. The main problem is that the spin winding number has to be probed for the second band of the Hamiltonian only: in a physical realization of the topological phase, both the first and the second band are filled, and the acquired signal includes information of both. The right column of figure 3 shows the expectation value of the spin $\langle \hat{{\boldsymbol{\Sigma }}}(k)\rangle$ obtained when the system with hard-wall boundaries is filled with 2L fermions (top panel) and with only L fermions (middle panel). Both signals are not particularly interesting. If we consider the artificial situation where atoms populate the second band only (bottom panel), the spin expectation value is characterized by a winding number that reproduces the behavior of ${\mathcal{W}}$ and encircles the origin in the topological phase. In the ideal case of a hard-wall confining potential, the required value of the second band can be extracted by repeating the experiment twice, at densities ρ = 1 and ρ = 2: the difference of the measured distributions returns the sought information $\langle \hat{{\boldsymbol{\Sigma }}}(k){\rangle }_{\rho =2}-\langle \hat{{\boldsymbol{\Sigma }}}(k){\rangle }_{\rho =1}=\frac{1}{2}{{\boldsymbol{S}}}_{2}(k).$

In the presence of a harmonic trap, the wedding cake density profile suggests that the many-body wavefunction can be roughly thought as a state where each energy band α is uniformly populated by N_α atoms (${N}_{1}\geqslant {N}_{2}\geqslant {N}_{3}\geqslant {N}_{4}\geqslant 0$). In this case the measurement of the observable $\langle \hat{{\boldsymbol{\Sigma }}}(k)\rangle$ returns:

$$\begin{eqnarray}&&\langle \hat{{\boldsymbol{\Sigma }}}(k)\rangle =\displaystyle \frac{1}{2}\displaystyle \sum _{\alpha }\displaystyle \frac{{N}_{\alpha }}{{\mathcal{L}}}{{\boldsymbol{S}}}_{\alpha }(k),\end{eqnarray} \tag{ 5 }$$

where ${{\boldsymbol{S}}}_{\alpha }(k)$ is the expectation value of the spin calculated in the thermodynamic limit for the single particle eigenstate of the ${\alpha }^{{\rm{th}}}$ energy band (see appendix). ${\mathcal{L}}$ is the discretization adopted for the Brillouin zone in the time-of-flight imaging (see, for example, [43, 46]).

If we consider the case in which the density profile shows only two plateaus, the value of ${{\boldsymbol{S}}}_{2}(k)$ can be estimated by comparing the observed $\langle \hat{{\boldsymbol{\Sigma }}}(k){\rangle }_{\rho =2}$ with that of a realization with a single plateau only, $\langle \hat{{\boldsymbol{\Sigma }}}(k){\rangle }_{\rho =1}:$

$$\begin{eqnarray}&&{{\boldsymbol{S}}}_{2}(k)=\displaystyle \frac{2{\mathcal{L}}}{{N}_{2}}\left({\langle \hat{{\boldsymbol{\Sigma }}}(k)\rangle }_{\rho =2}-\displaystyle \frac{{N}_{1}}{{N}_{1}^{\prime }}{\langle \hat{{\boldsymbol{\Sigma }}}(k)\rangle }_{\rho =1}\right),\end{eqnarray} \tag{ 6 }$$

where N₁ and N₂ are the occupations of the two bands for the state with two plateaus, and ${N}_{1}^{\prime }$ is the total number of atoms in the reference state with a single plateau. All the quantities N₁, N₂ and ${N}_{1}^{\prime }$ can be experimentally accessed and we report in the right column of figure 4 the comparison of the data obtained for hard wall and harmonic potentials. Our numerical simulations confirm that even in the presence of the trap ${{\mathcal{S}}}_{2}$ is equal to ±1 in the topological phase, whereas in the trivial phases, it is either 0 or ±2 (see figure 4).

Let us conclude with some information on how to measure $\langle \hat{{\boldsymbol{\Sigma }}}(k)\rangle$ through spin-resolved time-of-flight imaging [43–47] in our setup. Special care is required in time-modulated systems with spin-dependent features [63], as considered in the specific proposal detailed in section 5 because spin-dependent observables can potentially undergo large and complicated micro-motion (rapid motion with a time-scale of the order of the driving period $2\pi /\omega$), which typically alters the accuracy of measurements. In such schemes, stroboscopic measurements performed at specific times, $2\pi /\omega \times n$, where n is integer, are generally required to extract relevant information relative to the spin-dependent quantities [63]. For the scheme detailed in section 5, the micro-motion can be estimated from the unitary operators $K(t)$ and R(t) defined in that section, through the method of [53]. We find that $\langle {\hat{{\rm{\Sigma }}}}_{z}(k)\rangle$ is unaffected by the micro-motion; in contrast, an accurate analysis of $\langle {\hat{{\rm{\Sigma }}}}_{x}\rangle$ does require a stroboscopic measurement. Moreover, we note that measuring the expectation value of ${\hat{{\rm{\Sigma }}}}_{x}$ also necessitates a $\pi /2$ pulse, which has to be short compared to the driving period in order to probe the system stroboscopically.

Let us now consider the role of interactions, with a special emphasis on the robustness of symmetry-protected topological order. It is experimentally relevant to consider an on-site Hubbard interaction in each leg:

$$\begin{eqnarray}&&{\hat{H}}_{{\rm{int}}}=U\displaystyle \sum _{r,{\tau }_{z}}{\hat{n}}_{r,{\tau }_{z},{\sigma }_{z}=\uparrow }{\hat{n}}_{r,{\tau }_{z},{\sigma }_{z}=\downarrow },\end{eqnarray} \tag{ 7 }$$

and to analyze the phase diagram of ${\hat{H}}_{0}+{\hat{H}}_{{\rm{int}}}$ at half filling, $\rho =2,$ which is characterized by the competition between the TI and Mott insulator (MI) occurring in the presence of a strong contact repulsion. We employ a density-matrix renormalization group algorithm based on a matrix-product state (MPS) ansatz [64, 65]. We will consider systems with open-boundary conditions with L = 72 and maximal bond dimension D = 200.

The transition between TI and MI can be located via the analysis of the charge gap at N = 2L. In particular, the single-particle gap is defined as:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{1}(N)=E(N+1)+E(N-1)-2E(N),\end{eqnarray} \tag{ 8 }$$

where E(N) is the ground-state energy of the system with N fermions. Clearly, ${{\rm{\Delta }}}_{1}(N)\gt 0$ for the MI because the system has a thermodynamic gap. On the other hand, the TI has zero-energy modes which ensure that $E(N-1)=E(N)=E(N+1)$ and thus ${{\rm{\Delta }}}_{1}(N)=0.$ Unfortunately, the mere calculation of Δ₁ does not permit to discriminate the TI from a generic gapless phase, for which ${{\rm{\Delta }}}_{1}(N)=0$ too. We thus consider also the two-particle gap:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{2}(N)=E(N+2)+E(N-2)-2E(N).\end{eqnarray} \tag{ 9 }$$

Whereas for a gapless phase ${{\rm{\Delta }}}_{2}(N)$ is also equal to zero, for a TI it is larger than zero, signaling the gap which is protecting the phase.

Based on this discussion, we now consider a systematic study of the Hamiltonian, focusing on the competition of the two terms which are responsible for a gap opening, namely, the interaction term proportional to U responsible for the MI, and the interchain tunneling proportional to J. Roughly speaking, we identify the pairing term as the one inducing the TI, since at U = 0 the system is in a topological phase for $0\lt {J}^{2}\lt {J}_{c,1}^{2}\equiv {{\rm{\Omega }}}^{2}-{(2\mathrm{cos}(B/2)-\mu )}^{2}$ (see appendix for more detail) and the two chains decouple at J = 0.

Figure 5 presents the numerical results for ${{\rm{\Delta }}}_{1}(2L)$ and ${{\rm{\Delta }}}_{2}(2L)$ in the parameter space spanned by U/t and $J/t.$ The other parameters are chosen such that at U = 0 there is a TI, and are listed here for completeness: ${\rm{\Omega }}/t=1.8,$ $B=\frac{\pi }{2}\frac{L}{L+1}$ and μ/t = 1. Calculations are limited to the size L = 72 and a systematic scaling to the thermodynamic limit, as well as the exact evaluation of the properties of the critical lines, is beyond our numerical possibilities; additionally, the 2D space is studied with a grid of 0.2 along both axis. Despite these limitations, the qualitative nature of the phase diagram emerges quite clearly. Indeed, through the study of Δ₁ and Δ₂ we are able to identify the TI, the MI and the critical regions which separate them, resulting in the schematic phase diagram presented in figure 6. The topological region is identified with the large region where Δ₁ = 0 but Δ₂ > 0 whereas for the MI both Δ₁ and Δ₂ are larger than zero. Critical regions with Δ₁ = Δ₂ = 0 separate the two insulators. The schematic phase diagram discriminates the MI, whose appearance is driven by the on-site repulsion, from the trivial band insulator (BI), which appears instead also at U = 0. Our investigation did not identify a phase transition between these two trivial insulating phases, which are adiabatically connected.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To better analyze the transition between TI and MI, in figure 7 we focus on the line at $J/t=1,$ which entails a phase transition for $U/t={u}_{{cr}}\sim 1.9.$ We show the behavior of the chemical potentials $E(2L+\alpha )-E(2L+\alpha -1)$ for $\alpha =\;+2,+1,0,-1$ as a function of U. Two qualitatively different behaviors are separated by u_cr. For $U/t\gt {u}_{{cr}}$ the energy cost for adding one particle to the states with N = 2L or N = 2L + 1 is approximately U (especially for large values of U/t). Conversely, subtracting one particle from the states with N = 2L or N = 2L − 1 does not yield any energy gain. Thus, ${{\rm{\Delta }}}_{1}(N=2L)$ and ${{\rm{\Delta }}}_{2}(N=2L)$ are larger than zero and are approximately equal to U and 2U, respectively: these are typical signatures of a MI.

Zoom In Zoom Out Reset image size

For $U/t\lt {u}_{{cr}}$ the energy cost and gain for adding and removing one particle to/from the state with N = 2L are both equal to U/2. We interpret this as a signature of the fractionalization of the zero-energy modes of the TI: the charge excess $n\sim 1/2$ on top of the density plateau ρ = 2 (one particle per site) does cost a repulsive energy U n. Since the zero-energy modes have fermionic nature, they cannot accommodate more than one particle: the energy cost for adding one additional particle to the state with $N=2L+1$ becomes significantly larger than U/2 (vice versa for removing one particle from the state with N = 2L − 1). Thus, ${{\rm{\Delta }}}_{1}(N=2L)=0$ but ${{\rm{\Delta }}}_{2}(N=2L)\gt 0,$ signaling a TI for $U/t\lt {u}_{{cr}}.$

The MI extends for $U/t\gg 1,$ where a perturbative expansion shows that the system can be described by a spin model in a paramagnetic phase: under the assumptions $U\gg J,t$ and $\rho =2,$ one atom is trapped in each site of the two legs. We thus introduce the Pauli operators ${\hat{\eta }}_{r,\tau }^{i}$ (i = x, y, z) acting on the local effective Hilbert space spanned by the two spin states $\sigma =\pm 1$ of the atom located at the site r of the chain $\tau =\pm 1.$ We obtain the following second-order perturbative spin Hamiltonian:

$$\begin{eqnarray}{\hat{H}}_{{\rm{pert}}} & = & \displaystyle \sum _{r,\tau }\left[{\rm{\Omega }}{\hat{\eta }}_{r,\tau }^{x}+{J}_{{\rm{eff}}}{\hat{\eta }}_{r,\tau }^{z}{\hat{\eta }}_{r+1,\tau }^{z}+{J}_{{\rm{eff}}}\mathrm{cos}B\left({\hat{\eta }}_{r,\tau }^{x}{\hat{\eta }}_{r+1,\tau }^{x}+{\hat{\eta }}_{r,\tau }^{y}{\hat{\eta }}_{r+1,\tau }^{y}\right)\right.\\ & & \left.+{J}_{{\rm{eff}}}\mathrm{sin}B\left({\hat{\eta }}_{r,\tau }^{x}{\hat{\eta }}_{r+1,\tau }^{y}-{\hat{\eta }}_{r,\tau }^{y}{\hat{\eta }}_{r+1,\tau }^{x}\right)\right]+\displaystyle \sum _{r,j}{K}_{{\rm{eff}}}{\hat{\eta }}_{r,\tau =1}^{j}{\hat{\eta }}_{r,\tau =-1}^{j},\end{eqnarray} \tag{ 10 }$$

where ${J}_{{\rm{eff}}}\propto {t}^{2}/U$ and ${K}_{{\rm{eff}}}\propto {J}^{2}/U.$ In this regime, the term proportional to Ω dominates and the ground state of ${H}_{{\rm{pert}}}$ is close to a trivial product state in which all the spins are oriented in the $\hat{x}$ direction. We expect that such state, characterizing the MI phase, might be adiabatically connected to the trivial BI at U = 0 and $J\gt {J}_{c,1}.$ Our numerics does not suggest the existence of a further phase transition between the Mott and the trivial band insulating phases.

The phase diagram in figure 5 shows that the topological region appears clearly as a thermodynamic region, within a well defined parameter regime. We emphasize that, for $1.8\lesssim J\lesssim 2.2,$ the system is in a topologically trivial phase for U = 0, and enters the symmetry-protected topological phase when the interaction parameter U is increased. Therefore the interaction is not necessarily obnoxious to the purpose of experimentally obtaining the topological phase but, on the contrary, it can also drive the system into it by shifting the position of the critical point. This has been verified also in the corresponding TSC systems [35–37], where the addition of repulsive interactions is proven to expand the topological phase for certain ranges of the physical parameters (see also [39, 40] for related models in terms of Majorana modes). This means that, for some particular value of $J\gt {J}_{c,1}$ the presence of a repulsive interaction allows the formation of edge modes otherwise absent. A similar behavior is also observed in 2D systems with time-reversal invariance [66]. Let us stress, however, that this has nothing to do with the physics of fractional Chern insulators, where interactions drive the system into distinct (strongly correlated) topological phases. As the phase diagram clearly shows, there is only one TI phase, which is strictly equivalent to that of the non-interacting system. Importantly, the phase diagram in figure 5 shows that interactions have a non-trivial role in tuning the system in and out the TI phase.

In order to further clarify this last point, we now investigate in more detail the properties of the topological phase in the interacting system. Numerical investigations reported in figure 8 show that the signatures of the non-interacting TI persist in the presence of interactions. First, the density profile of the gas allows for a clear identification of the presence of fractionalized edge modes located at the boundaries of the ladder via the computation of $\langle {\hat{n}}_{j}\rangle$ and $\langle {\hat{n}}_{j}^{*}\rangle .$ Indeed, for $U/t\leqslant {u}_{{cr}},$ figure 8 shows that $\langle {\hat{n}}_{j}^{*}\rangle$ saturates to 0.5 within few sites, which is strongly different from the behavior for $U/t\geqslant {u}_{{cr}}.$ It is interesting to observe that within the topological region the localization length of the edge modes has a weak dependence on $U/t.$ Second, the system displays also within the interacting region a non-zero winding number associated with the second band of the system. As in an interacting system bands are not well defined, the winding number is computed by subtraction of the data relative to $\rho =1$ to those relative to $\rho =2$ (see similar discussion in section 3). This robustness of the spin winding number against local interaction is consistent with similar results in 2D systems [48]. Finally, on a more abstract side, the analysis of the Schmidt spectrum presents the robust two-fold degeneracy of symmetry-protected topological phases [67].

Zoom In Zoom Out Reset image size

The physical realization of the ladder system in equation (1) can be obtained by extending the 2D setup elaborated and realized in [4]. The present proposal builds on a 2D optical superlattice subjected to a well-designed time-modulation as displayed in figure 9. Along the y direction a superlattice potential is used to partition the lattice into a 1D array of isolated ladders. Hopping between the two legs of the ladder corresponds to transitions ${\hat{\tau }}_{z}=-1\leftrightarrow +1,$ see figure 1. The main challenge in realizing the Hamiltonian in equation (1) consists in engineering the spin-dependent complex matrix elements for tunneling processes taking place along the legs of the ladder. In the following, we will show that this can be achieved by combining a spin-dependent superlattice potential x (figure 9), which introduces a spin-dependent energy offset Δ between neighboring sites, together with the space-dependent time-modulation of the lattice discussed in [4].

Zoom In Zoom Out Reset image size

We start by considering the time-independent part of the system, which can be described by the 2D tight-binding Hamiltonian

$$\begin{eqnarray}{\hat{H}}_{0} & = & -{J}_{x}\displaystyle \sum _{m,n}\left({\hat{a}}_{m+1,n}^{\dagger }{\hat{a}}_{m,n}+{\rm{h.c.}}\right)\\ & & -{J}_{y}\displaystyle \sum _{m,n}\left({\hat{a}}_{m,n+1}^{\dagger }{\hat{a}}_{m,n}+{\rm{h.c.}}\right)+\displaystyle \frac{{\rm{\Delta }}}{2}\displaystyle \sum _{m,n}{(-1)}^{m}{\hat{a}}_{m,n}^{\dagger }{\sigma }_{z}{\hat{a}}_{m,n},\end{eqnarray} \tag{ 11 }$$

where m and n label the horizontal and vertical integer coordinates. The spin-dependent staggered potential could be realized, for instance, by considering an appropriate anti-magic wavelength, for which the polarizability is opposite for the two spin species [51]. In order to keep the bare tunneling processes of strength J_{x, y} spin-independent the remaining lattice potentials need to be created using a magic wavelength, for which the polarizability is the same for the two spin-species.

The bare tunneling is suppressed along the legs due to the offset ${\rm{\Delta }}\gg {J}_{x,y},$ which allows for a complete control over induced-tunneling-matrix elements, such as those realized by modulating the lattice resonantly [53]. Following [4], the modulation is taken to be produced by two pairs of laser beams with frequency difference ${\omega }_{r/b}=\pm {\rm{\Delta }}$ (figure 9) in order to restore resonant tunneling. The corresponding time-dependent potential defined by these four lasers is then of the form

$$\begin{eqnarray}&&\hat{V}(t)=\kappa \displaystyle \sum _{m,n}{\hat{a}}_{m,n}^{\dagger }{\hat{a}}_{m,n}\left[v(m,n){{\rm{e}}}^{{\rm{i}}\omega t}+{v}^{*}(m,n){{\rm{e}}}^{-{\rm{i}}\omega t}\right],\end{eqnarray} \tag{ 12 }$$

with the resonance condition ω = Δ, and we choose the laser phases in such a way that

$$\begin{eqnarray}&&v(m,n)=\displaystyle \frac{1}{2}\left\{\mathrm{cos}\left(m\displaystyle \frac{\pi }{2}-\pi /4\right){{\rm{e}}}^{-{\rm{i}}\pi n-{\rm{i}}B/2}+\mathrm{cos}\left(m\displaystyle \frac{\pi }{2}+\pi /4\right){{\rm{e}}}^{{\rm{i}}\pi n+{\rm{i}}B/2}\right\}.\end{eqnarray} \tag{ 13 }$$

This requires a stabilization of the phase of the modulation relative to the static lattice potential, which is challenging and was not yet demonstrated in previous realizations [2, 4, 49, 50]. This specific choice of the potential $v(m,n)$ is made in order to independently address successive hopping terms along the x direction, which is generally required when engineering Peierls phase-factors in superlattice structures, see [4, 53] and below.

The time-evolution of the system is ruled by the Schrödinger equation ${\rm{i}}{\partial }_{t}\psi =\hat{H}(t)\psi ,$ where $\hat{H}(t)={\hat{H}}_{0}+\hat{V}(t)$ is defined by equations (11) and (12). The long-time dynamics of the system can be suitably described by an effective-Hamiltonian approach [53], which is valid in the high-frequency regime $\omega \to \infty .$ Since the static Hamiltonian H₀ contains a staggered-potential term that explicitly diverges linearly with Δ = ω, we first apply the unitary transformation [53]

$$\begin{eqnarray}&&\psi =\hat{R}(t)\tilde{\psi }=\mathrm{exp}\left(-{\rm{i}}\hat{W}t\right)\tilde{\psi },\qquad \hat{W}=\displaystyle \frac{{\rm{\Delta }}}{2}\displaystyle \sum _{m,n}{(-1)}^{m}{\hat{a}}_{m,n}^{\dagger }{\sigma }_{z}{\hat{a}}_{m,n},\end{eqnarray} \tag{ 14 }$$

which removes the diverging term. The effective Hamiltonian can then be derived in this moving frame, using the method of [53, 63] (see also [68]).

For the sake of simplicity, let us first consider the dynamics associated with the species σ_z = +1. For these atoms, the transformed Hamiltonian reads:

$$\begin{eqnarray}&&\tilde{H}(t)={\hat{R}}^{\dagger }(t)\left[{\hat{H}}_{0}+\hat{V}(t)\right]\hat{R}(t)-\hat{W}={\hat{V}}^{+}{{\rm{e}}}^{{\rm{i}}\omega t}+{\hat{V}}^{-}{{\rm{e}}}^{-{\rm{i}}\omega t},\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}&&{\hat{V}}^{+}=\kappa \displaystyle \sum _{m,n}{\hat{n}}_{m,n}v(m,n)-{J}_{x}\displaystyle \sum _{m\;{\rm{odd}},n}\left({\hat{a}}_{m+1,n}^{\dagger }{\hat{a}}_{m,n}+{\hat{a}}_{m-1,n}^{\dagger }{\hat{a}}_{m,n}\right),\\ &&{\hat{V}}^{-}=\kappa \displaystyle \sum _{m,n}{\hat{n}}_{m,n}{v}^{*}(m,n)-{J}_{x}\displaystyle \sum _{m\;{\rm{even}},n}\left({\hat{a}}_{m+1,n}^{\dagger }{\hat{a}}_{m,n}+{\hat{a}}_{m-1,n}^{\dagger }{\hat{a}}_{m,n}\right).\end{eqnarray} \tag{ 16 }$$

We describe the time-evolution of the system dictated by $\tilde{H}(t)$ through the evolution operator, which we partition as

$$\begin{eqnarray}&&\tilde{U}(t)={{\rm{e}}}^{-{\rm{i}}\hat{K}(t)}{{\rm{e}}}^{-{\rm{i}}t{\hat{H}}_{{\rm{eff}}}}{{\rm{e}}}^{{\rm{i}}\hat{K}(0)},\end{eqnarray} \tag{ 17 }$$

where the effective Hamiltonian ${\hat{H}}_{{\rm{eff}}}$ describes the long-time dynamics, and where the operator $\hat{K}(t)$ captures the micro-motion. Following [63], we find that the effective Hamiltonian associated with the general time-dependent Hamiltonian in equation (15) is given by:

$$\begin{eqnarray}{\hat{H}}_{{\rm{eff}}} & = & \displaystyle \frac{1}{\omega }[{\hat{V}}^{(+1)},{\hat{V}}^{(-1)}]+{\mathcal{O}}(1/{\omega }^{2})\\ & = & -\displaystyle \frac{{J}_{x}\kappa }{\omega }\left[\displaystyle \sum _{m\;{\rm{even}},n}{\hat{a}}_{m+1,n}^{\dagger }{\hat{a}}_{m,n}\left(v(m+1,n)-v(m,n)\right)\right.\\ & & \left.\phantom{\displaystyle \sum _{n}}+{\hat{a}}_{m-1,n}^{\dagger }{\hat{a}}_{m,n}\left(v(m-1,n)-v(m,n)\right)\right]+{\rm{h.c.}}\\ & = & \displaystyle \frac{{J}_{x}\kappa }{\sqrt{2}\omega }\left[\displaystyle \sum _{m\;{\rm{even}},n}{\hat{a}}_{m+1,n}^{\dagger }{\hat{a}}_{m,n}\mathrm{cos}\displaystyle \frac{m\pi }{2}{{\rm{e}}}^{{\rm{i}}n\pi +{\rm{i}}B/2}\right.\\ & & \left.+{\hat{a}}_{m-1,n}^{\dagger }{\hat{a}}_{m,n}\mathrm{cos}\displaystyle \frac{m\pi }{2}{{\rm{e}}}^{-{\rm{i}}n\pi -{\rm{i}}B/2},\phantom{\displaystyle \frac{{J}_{x}\kappa }{\sqrt{2}\omega }\displaystyle \sum _{m{\rm{even}},n}},\phantom{\rule{-4.8em}{0ex}}\right]+{\rm{h.c.}}\\ & = & \displaystyle \frac{{J}_{x}\kappa }{\sqrt{2}\omega }\displaystyle \sum _{n}{{\rm{e}}}^{{\rm{i}}n\pi +{\rm{i}}B/2}\left[\displaystyle \sum _{m\;{\rm{even}}}\left({\hat{a}}_{2m+1,n}^{\dagger }{\hat{a}}_{2\;m,n}+{\hat{a}}_{2\;m,n}^{\dagger }{\hat{a}}_{2m-1,n}\right)\right.\\ & & +\left.\displaystyle \sum _{m\;{\rm{odd}}}\left(-{\hat{a}}_{2m+1,n}^{\dagger }{\hat{a}}_{2\;m,n}-{\hat{a}}_{2\;m,n}^{\dagger }{\hat{a}}_{2m-1,n}\right)\right]+{\rm{h.c.}}.\end{eqnarray} \tag{ 18 }$$

The irrelevant sign change in the tunneling matrix elements (i.e. in the last line of equation (18)), can be removed by applying an additional gauge transformation

$$\begin{eqnarray}&&\hat{G}=\mathrm{exp}\left[{\rm{i}}\pi \displaystyle \sum _{m\;\mathrm{odd},n}{\hat{a}}_{2\;m,n}^{\dagger }{\hat{a}}_{2\;m,n}\right]={\hat{G}}^{\dagger },\quad {\hat{G}}^{2}=1,\end{eqnarray} \tag{ 19 }$$

which indeed reverses the sign of the tunneling terms ${\hat{a}}_{2m+1,n}^{\dagger }{\hat{a}}_{2\;m,n}$ and ${\hat{a}}_{2\;m,n}^{\dagger }{\hat{a}}_{2m-1,n}$ for m odd only. In this way, the final effective Hamiltonian describing the dynamics of the σ_z = +1 species yields

$$\begin{eqnarray}&&{\hat{H}}_{{\rm{eff}}}\to \hat{G}{\hat{H}}_{{\rm{eff}}}\hat{G}=\displaystyle \frac{{J}_{x}\kappa }{\sqrt{2}\omega }\displaystyle \sum _{m,n}{\hat{a}}_{m+1,n}^{\dagger }{\hat{a}}_{m,n}{{\rm{e}}}^{{\rm{i}}n\pi +{\rm{i}}B/2}+{\rm{h.c.}}\end{eqnarray} \tag{ 20 }$$

which is indeed the tunneling term in equation (A.1) for ${\sigma }_{z}\;=\;+1$ atoms. In the case of the ${\sigma }_{z}=-1$ species, the staggered potential is reversed, so that the even and odd sites must be inverted. This results in the final effective Hamiltonian:

$$\begin{eqnarray}&&{\hat{H}}_{{\rm{eff}}}=\displaystyle \frac{{J}_{x}\kappa }{\omega }\displaystyle \sum _{m,n}\left({(-1)}^{n}{\hat{a}}_{m+1,n}^{\dagger }{{\rm{e}}}^{{\rm{i}}B{\sigma }_{z}/2}{\hat{a}}_{m,n}+{\rm{h.c.}}\right),\end{eqnarray} \tag{ 21 }$$

where ${(-1)}^{n}$ is equivalent to the operator ${\tau }_{z}$ in equation (A.1).

The Zeeman term Ω σ_x present in the Hamiltonian (1) can be directly generated by two resonant coupling potentials ${\hat{V}}_{i}^{\mathrm{RC}}(t)={\displaystyle \sum }_{m,n}2{\rm{\Omega }}\mathrm{cos}({\nu }_{i}t){\hat{a}}_{m,n}{\sigma }_{x}{a}_{m,n}.$ Indeed, considering the bare atomic frequency ${\omega }_{z}\gg {\rm{\Delta }}$ between the two sublevels, they are effectively separated by the position dependent energy offset ${\omega }_{z}+{(-1)}^{m}{\rm{\Delta }}$ created by the spin-dependent potential in equation (11). Under the transformation $\hat{R}(t),$ the effect of the coupling becomes

$$\begin{eqnarray}&&{\hat{R}}^{\dagger }(t){\hat{V}}_{i}^{\mathrm{RC}}(t)\hat{R}(t)=\displaystyle \sum _{m,n}2{\rm{\Omega }}{{\rm{e}}}^{-{\rm{i}}{\rm{\Delta }}{(-1)}^{{mt}}}\mathrm{cos}\left({\nu }_{i}t\right){\hat{a}}_{m,n}^{\dagger }{\sigma }_{+}{\hat{a}}_{m,n}+{\rm{h.c.}},\end{eqnarray} \tag{ 22 }$$

where ${\sigma }_{+}=({\sigma }_{x}+{\rm{i}}{\sigma }_{y})/2.$ This terms commute with the gauge-transformation operator $\hat{G}$ and, by choosing the frequencies ${\nu }_{\mathrm{1,2}}={\omega }_{z}\pm {\rm{\Delta }},$ we recover the required Zeeman term ${\rm{\Omega }}{\sigma }_{x}$ in equation (1), both for even and odd sites, through the standard rotating-wave approximation.

We note that the effective Hamiltonian (1) was derived at first order in ${\omega }^{-1},$ in a basis provided by two commuting unitary operators: $\hat{R}(t)$ and $\hat{G}.$ It is important to notice that these latter operators commute with the τ operators, so that they neither affect the static hopping term $J{\tau }_{x},$ nor the static potential difference $\mu {\tau }_{z}:$ this indicates that these static terms can be directly included in the (effective) Hamiltonian (1). The latter remark is also valid for the spin-independent potential μ₀. Therefore, we conclude that the application of two pairs of Raman lasers and a radio-frequency field, combined with the spin-dependent staggered potential directed along the ladder, allows one to generate all the spin-dependent terms in the ladder Hamiltonian (1).

Importantly, we emphasize that the Hubbard interactions are also unaffected by the aforementioned transformations $\hat{G}$ and $\hat{R}(t).$ Thus, the effects of interactions can be directly incorporated into the effective Hamiltonian (1), at first order in ω⁻¹.

Finally, the time-evolution operator in equation (17) is then completely determined by computing the kick operator $\hat{K}(t),$ which is readily calculated using the expression [53, 63]

$$\begin{eqnarray}&&\hat{K}(t)=\displaystyle \frac{1}{{\rm{i}}\omega }\left[{\hat{V}}^{+}{{\rm{e}}}^{{\rm{i}}\omega t}-{\hat{V}}^{-}{{\rm{e}}}^{-{\rm{i}}\omega t}\right]\approx \displaystyle \frac{2\kappa }{\omega }\displaystyle \sum _{m,n}{\hat{n}}_{m,n}\;| v(m,n)| \mathrm{sin}(\omega t+{\theta }_{m,n}),\end{eqnarray} \tag{ 23 }$$

where we assumed that $\kappa \gg {J}_{x}$ and we defined ${\theta }_{m,n}=\mathrm{arg}[v(m,n)].$

In this work we presented a ladder setup for ultracold fermions subject to both the presence of an artificial π-flux magnetic potential and a spin–orbit coupling. Such a model may be seen as the synthesis of two accessible experimental techniques to realize synthetic gauge fields in optical lattices: on one side, the realization of complex tunneling matrix elements using time-modulated optical lattice [2–4, 49, 50], and, on the other, the implementation of spin–orbit terms through the control over internal atomic degrees of freedom [5, 6]. Analogously to models already discussed in the context of nanowires [29, 30], the combination of these two elements gives rise to a particle–hole symmetry which protects non-trivial topological phases within the Bogoliubov–de Gennes and chiral classes of TIs and superconductors. As the model conserves the number of particles and can be realized with state-of-the-art experimental techniques, our results might give a substantial advance towards the observation of topologically protected zero-energy modes in fermionic systems with and without interactions. In particular the physical realization that we present does not require the engineering of any pairing mechanism, neither the coupling to external molecular gases or superfluids, as, for example in [9, 11, 12, 14], nor an interchain pair-hopping, as exploited in the ladder model presented in [13].

No interaction is indeed necessary for the appearance of the symmetry-protected topological phase. Consequently, the realization of the particle–hole symmetry through the ladder geometry is not as robust as its counterpart in TIs due to the absence of a true superconducting gap, thus it should be considered an extrinsic feature. However, due to the absence of disorder and to the high degree of isolation in ultracold atom setups, we expect the topological features of the system to be experimentally detectable. For example, the introduction of a trapping potential does not spoil the observation of edge modes, despite breaking the particle–hole symmetry.

We analyzed two experimentally relevant signatures of the appearance of topological phases: the presence of fractionalized edge modes, detectable through site-resolved density measurements (as recently reported in [59–61]), and the winding behavior of the spin degree of freedom, which can be observed through spin-resolved time-of-flight imaging. In particular, we have shown how to detect these observables even in the presence of a trapping potential, which induces soft boundaries, often believed to be particularly disruptive for the detection of topological signatures.

Our study has also considered the effect of a contact repulsive interaction; this is possible thanks to the presence of the spin degree of freedom that differentiate the main features of our model from its spinless counterparts as the SSH and its interacting generalizations (see for example [41, 42]). Apart from mapping out the phase-diagram of the model, which entails two gapped phases, with and without topological properties, we have found that a Hubbard interaction can enhance the extension of the topological phase, instead of being detrimental. Furthermore, the spin winding number introduced in the article provides a good topological parameter also in the interacting case, where the usual order parameters based on single-particle wavefunctions fail.

Concluding, we mention that the model under scrutiny may be an interesting platform for the study of further fractionalization effects, based on a particular fine tuning of the parameters, which can be reminiscent of the physics of parafermionic zero-energy modes (see [69] for a recent review) in the spirit of [30].

We acknowledge enlightening discussions with Rosario Fazio, Fabian Hassler, Sylvain Nascimbene, Luca Taddia and Andrea Trombettoni and thank Simone Barbarino for careful reading of the manuscript. We thank Davide Rossini for providing the numerical algorithms based on Matrix-Product States employed in this work. LM acknowledges support from the Italian MIUR through FIRB project RBFR12NLNA and Regione Toscana POR FSE 2007–2013. MA, H-HT and MB acknowledge support from the EU grant SIQS. NG is financed by the FRS-FNRS Belgium and by the BSPO under the PAI project P7/18 DYGEST.

In this appendix we examine the non-interacting ladder model and we discuss, in particular, the relation between the order parameter ${\mathcal{W}},$ that distinguishes trivial and topological phases, and the observed spin winding number.

The Hamiltonian (1) is translationally invariant and can be also expressed in momentum space as $\hat{H}={\displaystyle \sum }_{k}{\hat{a}}_{k}^{\dagger }{\mathcal{H}}(k){\hat{a}}_{k}$ where:

$$\begin{eqnarray}&&{\mathcal{H}}(k)=2t{\tau }_{z}\mathrm{cos}(B/2)\mathrm{cos}k-2t{\tau }_{z}{\sigma }_{z}\mathrm{sin}(B/2)\mathrm{sin}k+{\rm{\Omega }}{\sigma }_{x}+J{\tau }_{x}+\mu {\tau }_{z}+{\mu }_{0}.\end{eqnarray} \tag{ A.1 }$$

We observe that its kinetic term corresponds to a Peierls substitution $k\to k+\frac{B{\sigma }_{z}}{2}.$

Its spectrum is symmetric for a tranformation mapping $B\to B+2\pi$ and $t\to -t.$ Therefore, even if the Hamiltonian is periodic in B with period $4\pi ,$ we can restrict our study to the case $0\leqslant B\lt 2\pi .$ Besides, we can consider only positive values of Ω and J because their sign trivially depends on the chosen basis for the spin and pseudospin. In particular ${\mathcal{H}}(k,-{\rm{\Omega }})={\sigma }_{z}{\mathcal{H}}(k,{\rm{\Omega }}){\sigma }_{z}$ and ${\mathcal{H}}(k,-J)={\tau }_{z}{\mathcal{H}}(k,J){\tau }_{z}.$ Analogously we consider only $t\gt 0,$ because of the symmetry between the two chains ${\mathcal{H}}(k,-t,-\mu )={\tau }_{x}{\mathcal{H}}(k,t,\mu ){\tau }_{x}.$ Hereafter we rescale all the energies in units of t in such a way that, below, we will always consider t = 1. Finally we observe that the expectation value of ${\sigma }_{y}$ is always null because the Hamiltonian (A.1) is real.

The system is characterized by an anti-unitary time-reversal-like symmetry $T={\sigma }_{x},$

$$\begin{eqnarray}&&T{\mathcal{H}}(k){T}^{\dagger }={{\mathcal{H}}}^{*}(-k),\end{eqnarray} \tag{ A.2 }$$

and, for ${\mu }_{0}=0$ (the system is exactly at half filling), we obtain the particle–hole symmetry $C={\tau }_{y}{\sigma }_{y}:$

$$\begin{eqnarray}&&C{\mathcal{H}}(k){C}^{\dagger }=-{{\mathcal{H}}}^{*}(-k).\end{eqnarray} \tag{ A.3 }$$

These non-unitary symmetries characterize the topological symmetry class BDI (see, for example, [23, 24]), which is also characterized by the unitary chiral symmetry $P={TC}={\tau }_{y}{\sigma }_{z}.$

Additional terms in the Hamiltonian may break the C symmetry which is indeed fragile, due to the lack of a physical pairing interaction; we emphasize however that such perturbations become significant only if their magnitude is comparable with the energy gap. To this concern, in a cold atom gas, the presence of noise and defects are negligible and the main effect we must consider is the trapping potential bringing to a space dependent chemical potential ${\mu }_{0}.$ As discussed in section 3, however, the local shift in energy provided by the trapping allows in general to isolate topological regions of the chain with the effect of binding fractionalized modes at the interface between these regions and the trivial ones.

An additional Zeeman term proportional to σ_z, which breaks only the time reversal symmetry, brings instead the system in the symmetry class D, which is still topologically non-trivial in one dimensions. Therefore, this sort of term does not alter in a fundamental way the properties of the system.

1D systems in the BDI class possess topological phases labelled by a topological invariant in ${\mathbb{Z}}$ [23]. This topological invariant, can be evaluated by exploiting the chiral basis defined by the symmetry P [56]. In this basis the Hamiltonian assumes the simple form:

$$\begin{eqnarray}{\mathcal{H}}(k) & = & \left(\begin{array}{cc}0 & A(k)\\ {A}^{\dagger }(k) & 0\end{array}\right);\\ A & = & \left(\begin{array}{cc}{\rm{i}}\mu +J+2\mathrm{icos}\left(\displaystyle \frac{B}{2}+k\right) & {\rm{i}}{\rm{\Omega }}\\ {\rm{i}}{\rm{\Omega }} & {\rm{i}}\mu -J+2\mathrm{icos}\left(\displaystyle \frac{B}{2}-k\right)\end{array}\right).\end{eqnarray} \tag{ A.4 }$$

If $\mathrm{det}(A(k))\ne 0$ for all k then the system is gapped. In this case, $\mathrm{det}A(k)=\left|\mathrm{det}A(k)\right|{{\rm{e}}}^{{\rm{i}}\xi (k)}$ and the winding number of $\xi (k)$ constitutes the topological invariant which distinguish topological and non-topological phases [56]. In this model, this winding number may assume only the values 0 or ±1 characterizing trivial and topological phases respectively. In more detail, one obtains:

$$\begin{eqnarray}\mathrm{det}A(k) & = & -{J}^{2}-{\mu }^{2}+{{\rm{\Omega }}}^{2}-2\mathrm{cos}(B)-2\mathrm{cos}(2k)\\ & & +-4\mu {cos}(B/2)\mathrm{cos}(k)+4{\rm{i}}J\mathrm{sin}(B/2)\mathrm{sin}(k).\end{eqnarray} \tag{ A.5 }$$

This determinant is purely real for $k=0,\pi$ or for B = 0. For B = 0, its phase $\xi (k)$ cannot change its value from 0 to π unless $\mathrm{det}A$ crosses zero and the gap closes; thus the case B = 0 is either trivial or gapless.

In all the other cases, for $0\lt B\lt 2\pi ,$ the topological invariant ${\mathcal{W}}$ is evaluated by considering the behavior of ξ between k = 0 and k = π. In particular, the tangent of the phase ξ is given by:

$$\begin{eqnarray}&&\mathrm{tan}\left[\xi (k)\right]=\displaystyle \frac{4J\mathrm{sin}\displaystyle \frac{B}{2}\mathrm{sin}k}{{{\rm{\Omega }}}^{2}-4\mu \mathrm{cos}\displaystyle \frac{B}{2}\mathrm{cos}k-2\mathrm{cos}B-{J}^{2}-2\mathrm{cos}(2k)-{\mu }^{2}}\end{eqnarray} \tag{ A.6 }$$

where the numerator is always positive for $B\in (0,2\pi )$ and $k\in (0,\pi ),$ and always negative for $B\in (0,2\pi )$ and $k\in (\pi ,2\pi ).$ Therefore, going from k = 0 to $k=\pi ,$ $\xi (k)$ must be always included in $(0,\pi )$ since $\mathrm{sin}(\xi )\gt 0$ in this regime, whereas for k which goes from π to $2\pi ,$ $\xi (k)$ must be either in $(-\pi ,0)$ (if $\xi (k=\pi )=0$), such that its winding number vanishes for $k\to 2\pi ,$ or $\xi (k)\in (\pi ,2\pi )$ which implies a final winding number equal to one (see figure A1 ). Thus the parity of the winding number results:

$$\begin{eqnarray}&&{\mathcal{W}}={\rm{sign}}\left[\displaystyle \frac{{{\rm{\Omega }}}^{2}-2-{J}^{2}-{\mu }^{2}-4\mu \mathrm{cos}(B/2)-2\mathrm{cos}(B)}{{{\rm{\Omega }}}^{2}-2-{J}^{2}-{\mu }^{2}+4\mu \mathrm{cos}(B/2)-2\mathrm{cos}(B)}\right].\end{eqnarray} \tag{ A.7 }$$

Zoom In Zoom Out Reset image size

${\mathcal{W}}=1$ when the winding number is 0, whereas ${\mathcal{W}}=-1$ for the winding number being ±1. Therefore, in terms of the parameter Ω, two phase transitions appear at:

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{c,1}^{2}\equiv {J}^{2}+\mathrm{min}\left[{\left(\mu -2\mathrm{cos}\displaystyle \frac{B}{2}\right)}^{2},{\left(\mu +2\mathrm{cos}\displaystyle \frac{B}{2}\right)}^{2}\right]\end{eqnarray} \tag{ A.8 }$$

and

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{c,2}^{2}\equiv {J}^{2}+\mathrm{max}\left[{\left(\mu -2\mathrm{cos}\displaystyle \frac{B}{2}\right)}^{2},{\left(\mu +2\mathrm{cos}\displaystyle \frac{B}{2}\right)}^{2}\right].\end{eqnarray} \tag{ A.9 }$$

The system is in a topological phase $({\mathcal{W}}=-1)$ for ${{\rm{\Omega }}}_{c,1}^{2}\lt {{\rm{\Omega }}}^{2}\lt {{\rm{\Omega }}}_{c,2}^{2}$ whereas for ${{\rm{\Omega }}}^{2}\lt {{\rm{\Omega }}}_{c,1}^{2}$ and ${{\rm{\Omega }}}^{2}\gt {{\rm{\Omega }}}_{c,1}^{2}$ we obtain trivial phases $({\mathcal{W}}=1).$ At the transition points ${{\rm{\Omega }}}^{2}={{\rm{\Omega }}}_{c,1}^{2},{{\rm{\Omega }}}_{c,2}^{2},$ the gap closes respectively for $k=\pi ,0,$ consistently with the results in [30, 70]. The corresponding critical values for J are given by:

$$\begin{eqnarray}&&{J}_{1,2}^{2}={{\rm{\Omega }}}^{2}-{\left(\mu \mp 2\mathrm{cos}(B/2)\right)}^{2}\end{eqnarray} \tag{ A.10 }$$

Depending on ${\rm{\Omega }},\mu$ and B there can be 0, 1 or 2 phase transitions as a function of J.

In order to evaluate the spin expectation value in the plane $\hat{x}-\hat{z},$ it is useful to adopt the Hamiltonian form in equation (A.4) which allows a simple diagonalization. In this basis the components of the physical spin, S^x and S^z, become:

$$\begin{eqnarray}&&{S}^{x}={\sigma }_{x}\otimes {\tau }_{y},\qquad {S}^{z}={\sigma }_{z}\otimes {\tau }_{0},\end{eqnarray} \tag{ A.11 }$$

where σ_i and ${\tau }_{i}$ are Pauli matrices in this basis. To obtain the eigenstates of the Hamiltonian we consider:

$$\begin{eqnarray}{{\mathcal{H}}}^{2}(k)=\left(\begin{array}{ll}A(k){A}^{\dagger }(k) & 0\\ 0 & {A}^{\dagger }(k)A(k)\end{array}\right),\end{eqnarray} \tag{ A.12 }$$

where, ${{AA}}^{\dagger }={({A}^{\dagger }A)}^{*}$ since $A={A}^{T}.$ Therefore the generic form of the eigenstates in this basis is:

$$\begin{eqnarray}&&{\rm{\Psi }}=\left(\begin{array}{l}\psi \\ \pm {{\rm{e}}}^{{\rm{i}}\alpha }{\psi }^{*}\end{array}\right),\end{eqnarray} \tag{ A.13 }$$

where ψ are the two eigenvectors of ${{AA}}^{\dagger }.$ In particular, we are interested in the eigenvector corresponding to the lowest eigenvalue of ${{AA}}^{\dagger },$ since it determines the second and the third bands of H which are topologically non-trivial. The value of α is fixed by the equation:

$$\begin{eqnarray}&&{\psi }^{\dagger }A{\psi }^{*}=\pm \varepsilon {{\rm{e}}}^{-{\rm{i}}\alpha },\end{eqnarray} \tag{ A.14 }$$

where $\pm \varepsilon$ are the energies of the two intermediate bands of the system.

The expectation value of S^z in the second band does not depend on the phase α and it is easily evaluated. It results:

$$\begin{eqnarray}&&{S}_{2}^{z}(k)\equiv \left\langle{S}^{z}\right\rangle=\displaystyle \frac{2\mathrm{sin}(k)\left(\mathrm{sin}(B)\mathrm{cos}(k)+\mu \mathrm{sin}\left(\displaystyle \frac{B}{2}\right)\right)}{N(k)},\end{eqnarray} \tag{ A.15 }$$

where N(k) is a positive quantity. In particular we must distinguish two cases: $| \mu | \gtrless \left|\mathrm{cos}B/2\right|.$

For $\mu \gt 2\left|\mathrm{cos}B/2\right|,$ ${S}_{2}^{z}\gt 0$ for $k\;\in \;(0,\pi ),$ thus S^z₂ behaves like the imaginary component of the phase $\mathrm{exp}[{\rm{i}}\xi (k)]$ and it can be showed that the sign of S^x₂ at k = 0 and $k=\pi$ is equal to the sign of the denominator of equation (A.6), therefore the spin winding number of the second band and the winding number of the phase $\xi (k)$ coincide.

For $0\lt \mu \lt 2\left|\mathrm{cos}B/2\right|,$ ${S}_{2}^{z}(k)$ has a further zero between k = 0 and $k=\pi ,$ therefore the spin vector is aligned along $\hat{x}$ in three different points in the interval $k\in [0,\pi ].$ This can be seen in the last column of figure 4 (obtained for $\mu =0.5$ and $B=\pi /2$), where the blue line, representing the average value of the spin in a system without the trapping, always crosses the horizontal axis three times. The spin winding number ${{\mathcal{S}}}_{2}$ is null if and only if ${S}_{2}^{x}(k)$ has the same sign in these three points (first and last rows). However if ${S}_{2}^{x}(k)$ assumes the same sign for k = 0 and $k=\pi ,$ then the total spin winding number for k going from 0 to $2\pi$ is even, due to the time-reversal symmetry T. Otherwise it is odd. At k = 0 and $k=\pi$ the sign of ${S}_{2}^{x}(k)$ assumes, respectively, the same value of the signs of the numerator and denominator of ${\mathcal{W}}$ in equation (A.7). Therefore the parity of the spin winding number, ${(-1)}^{{{\mathcal{S}}}_{2}},$ always coincides with ${\mathcal{W}}.$

Concerning the behavior at the intermediate zero of ${S}_{2}^{z},$

$$\begin{eqnarray}&&{k}_{c}=\mathrm{arccos}\left(-\mu \mathrm{sin}(B/2)/\mathrm{sin}(B)\right),\end{eqnarray} \tag{ A.16 }$$

one has ${S}_{2}^{x}({k}_{c})=\mathrm{cos}\left[\alpha ({k}_{c})\right]\propto {\rm{\Omega }}-J,$ therefore in the trivial region where ${{\rm{\Omega }}}^{2}\gt {{\rm{\Omega }}}_{c,2}^{2}$ the spin winding number is trivial. The other trivial region, ${{\rm{\Omega }}}^{2}\lt {{\rm{\Omega }}}_{c,1}^{2},$ is instead divided into two subregimes characterized by a spin winding number ${{\mathcal{S}}}_{2}=0$ or ${{\mathcal{S}}}_{2}=2$ separated in $J={\rm{\Omega }}.$ This is shown in the last column of figure 4 where the blue lines either does not wind around the origin (first row, corresponding to ${\rm{\Omega }}\lt J\lt {{\rm{\Omega }}}_{c,1}$) or winds twice around it (second row, ${{\rm{\Omega }}}_{c,1}\gt {\rm{\Omega }}\gt J$).

Similar conclusions hold for negative values of μ and it can be shown that the behavior of the spin winding number is not affected by open or periodic boundary conditions.