Rice–Mele model with topological solitons in an optical lattice

Conducting polymers [1] are particularly interesting one-dimensional systems due to their unusual topological properties characterized by a nontrivial Zak phase [2, 3], degenerate ground states, topological solitons [4] and a fractional charge [5]. Such polymers can be described in a simple manner by fermions moving in a lattice with dimerized tunneling amplitudes [4, 6]. The corresponding models can be constructed with ultracold atoms in optical lattices which give unprecedented tunability and control over the system parameters. Recently, ultracold bosons in optical superlattices were used to prepare experimentally [7] a model of conducting polymers (namely the Rice–Mele (RM) model [6], one of the simplest 1D models of nontrivial topology) and the corresponding Zak phase was measured. In a parallel work topological edge states in a similar potential were studied theoretically [8].

Models based on superlattice potentials are relatively easy to realize in experiment, however, they have also some limitations. An optical lattice potential is typically defect-free due to its origin, which creates a difficulty in realizing topological solitons. Such solitons typically emerge on defects that are the domain walls between topologically distinct phases. Forming the signature of nontrivial topology, they are the essence of the RM model.

In the present paper, we show how to realize the RM model with controlled defects using a system of attractive ultracold fermions [9–11] in a simple shaken one-dimensional optical lattice. Shaking, i.e. the periodic driving of system parameters (e.g. the optical potential depth or position), has been successfully implemented in cold atomic systems in order to induce various effects [12–17] following the seminal proposition [18]. In [19] we have shown that such a shaking combined with attractive interactions in a two-dimensional triangular lattice can result in an emergent Dice structure with topological properties. Here, we show that in the case of a one-dimensional system there exist a regime of parameters where atoms self-organize into a dimerized structure. The ground state is then two-fold degenerate. The corresponding states represent two topologically distinct dimerized configurations. Due to the emergent nature of the dimerized state, both configurations, separated by domain walls, may be simultaneously present in the lattice. Moreover, by controlling the filling fraction, impurities may be added to the configurations. Such defects—domain walls and impurities—naturally give rise to topologically protected solitons or bound states with a fractionalized particle number.

Our system consists of a two-species (denoted as $\downarrow ,\uparrow$) fermionic mixture trapped in an optical lattice potential ${{V}_{{\rm latt}}}={{V}_{\parallel }}{{{\rm sin} }^{2}}(\pi x/a)+{{V}_{\bot }}({{{\rm sin} }^{2}}(\pi y/a)+{{{\rm sin} }^{2}}(\pi z/a))$, where a is the lattice constant. For ${{V}_{\bot }}\gg {{V}_{\parallel }}$ the system is effectively one-dimensional. To control the system we use a familiar lateral (horizontal) lattice shaking [18]. Importantly, however, we introduce also periodic changes in the potential depth which we call here vertical shaking: ${{V}_{\parallel }}={{V}_{0}}+\delta {{V}_{0}}{\rm cos} \;\omega t$. δV₀ is an amplitude of the lattice depth shaking and ω the frequency, common to the lateral and vertical shaking. We assume fermionic species of equal mass, M, with different fillings: ${{n}_{\downarrow }}\approx 1$ and ${{n}_{\uparrow }}\approx 1/2$. The interaction between atoms of different species is assumed to be attractive. In effect fermions of different spin tend to pair creating composites [20] with the density given by the minority $\uparrow$-fermion density ${{n}^{\uparrow }}$. We include p-bands in the model and effectively have the composites that occupy s-bands and excess $\downarrow$-fermions that may occupy both p and s-bands.

The Rice–Mele model [6] contains two essential ingredients: two types of sites and asymmetric couplings between them. The former is realized in our model by a density-wave self-arranged configuration of composites. Such a configuration is energetically favorable when intra-band tunnelings are switched off by appropriate adjustment of shaking amplitude while making density-dependent inter-band tunneling resonant by adjusting the shaking frequency. To obtain asymmetric coupling with nearest neighbours the additional vertical shaking is necessary with an appropriate phase shift with respect to standard lateral shaking. This phase difference breaks the left–right symmetry of the problem.

To write the effective Hamiltonian of the model, we construct the time dependent Hamiltonian, H(t) and average it in time [18]. The minimal Hamiltonian of our system contains tunnelings, density induced tunnelings, renormalized interactions and shaking: $\hat{H}={{\hat{H}}_{{\rm tun}}}+{{\hat{H}}_{{\rm dit}}}+{{\hat{H}}_{\operatorname{int}}}+{{\hat{H}}_{{\rm sh}}}(t),$ where:

$$\begin{eqnarray} {{{\hat{H}}}_{{\rm tun}}} & = & {{J}_{0}}\mathop{\sum }\limits_{\langle ij\rangle }\left[ \hat{s}_{i}^{\dagger }{{{\hat{s}}}_{j}}+\widehat{{{s}_{\uparrow }}}_{i}^{\dagger }{{\widehat{{{s}_{\uparrow }}}}_{j}} \right]+{{J}_{1}}\mathop{\sum }\limits_{\langle ij\rangle }\hat{p}_{i}^{\dagger }{{{\hat{p}}}_{j}}, \\ {{{\hat{H}}}_{{\rm dit}}} & = & \mathop{\sum }\limits_{\langle ij\rangle }\left[ {{T}_{0}}\widehat{{{s}_{\uparrow }}}_{i}^{\dagger }(\hat{n}_{i}^{{}}+\hat{n}_{j}^{{}}){{\widehat{{{s}_{\uparrow }}}}_{j}}+{{T}_{1}}\hat{p}_{i}^{\dagger }(\hat{n}_{i}^{\uparrow }+\hat{n}_{j}^{\uparrow }){{{\hat{p}}}_{j}} \right. \\ {} & {} & +\left. {{T}_{01}}((j-i)\hat{p}_{i}^{\dagger }\hat{n}_{i}^{\uparrow }\hat{s}_{j}^{{}}+h.c) \right], \\ {{{\hat{H}}}_{\operatorname{int}}} & = & U_{0}^{{}}\mathop{\sum }\limits_{i}\hat{n}_{i}^{\uparrow }{{{\hat{n}}}_{i}}+U_{1}^{{}}\mathop{\sum }\limits_{i}\hat{p}_{i}^{\dagger }{{{\hat{p}}}_{i}}\hat{n}_{i}^{\uparrow }+{{E}_{1}}\mathop{\sum }\limits_{i}\hat{p}_{i}^{\dagger }{{{\hat{p}}}_{i}}, \\ {{{\hat{H}}}_{{\rm sh}}}(t) & = & K{\rm cos} \omega t\mathop{\sum }\limits_{j}j(\hat{n}_{j}^{\uparrow }+\hat{s}_{j}^{\dagger }{{{\hat{s}}}_{j}}+\hat{p}_{j}^{\dagger }{{{\hat{p}}}_{j}})+\delta {{E}_{1}}{\rm cos} (\omega t+\varphi )\mathop{\sum }\limits_{i}\hat{p}_{i}^{\dagger }\hat{p}_{i}^{{}}. \\ \end{eqnarray} \tag{ 1 }$$

Here, $\hat{s}_{i}^{\dagger },\hat{s}_{i}$, $\hat{p}_{i}^{\dagger },\hat{p}_{i}^{{}}$ are creation and annihilation operators of $\downarrow$-fermions in the s- and p-bands respectively, while $\hat{s}_{\uparrow i}^{\dagger },\hat{s}_{\uparrow i}^{{}}$ are s-band creation and annihilation operators for the $\uparrow$-fermion. Accordingly, ${{\hat{n}}_{i}},\hat{n}_{i}^{p}$ and $\hat{n}_{i}^{\uparrow }$ are the corresponding number operators. In the on-site interaction Hamiltonian, ${{\hat{H}}_{\operatorname{int}}}$, the renormalized self-energy of the composites is denoted by U₀, the on-site renormalized interaction between the composite and an excess $\downarrow$-fermion is characterized by U₁, and E₁ is the energy of the p-band. The modulus of a negative U₀ value is the largest interaction energy scale assumed in the model with $\mid {{U}_{1}}\mid \lt \mid {{U}_{0}}\mid$ (this follows from properties of Wannier functions, see appendix A). Attractive interactions between two species lead then to the creation of composites in the s band. All the tunneling amplitudes are assumed to be much smaller than ω. The range of ω values will be set by the required resonance condition (2).

${{\hat{H}}_{{\rm tun}}}$ corresponds to standard tunnelings with amplitudes J₀ and J₁ for the s and p bands while ${{\hat{H}}_{{\rm dit}}}$ describes often neglected density induced tunnelings ([21–24]) with amplitudes T_i. Observe that $\uparrow$-tunneling from site i to site j may happen only when there is a composite on site i and a free $\downarrow$-fermion in the s-band on site j. In effect, this tunneling creates a composite on site j. In the case of $\downarrow$-fermions, the analogical situation does not take place because there are no free $\uparrow$-fermions. The last—most important for the mechanism discussed later—term of ${{H}_{{\rm dit}}}$ couples s and p levels and describes process occurring when a composite-empty site adjoins a composite-occupied one (for details see appendix A).

${{\hat{H}}_{{\rm sh}}}(t)$ is a time-periodic Hamiltonian with K denoting the amplitude of the lateral lattice shaking while δE₁ denotes the strength of the time-variation of single-particle energy in the p-band which is induced by periodic driving of the lattice depth while φ is a relative phase between the lateral and vertical drivings. Additional effects due to the vertical shaking that are negligible for moderate δV₀ are discussed in appendices.

Next, we describe the averaging process (see appendix B for more details). First, we apply the unitary transformation, $\hat{U}={\rm exp} [-{\rm i}{{\hat{H}}_{\operatorname{int}}}t-{\rm i}\int _{0}^{t}{{\hat{H}}_{{\rm sh}}}(t^{\prime} ){\rm d}t^{\prime} ]$ that gives us a new Hamiltonian $\hat{H}^{\prime} ={{\hat{U}}^{\dagger }}\hat{H}\hat{U}-{\rm i}{{\hat{U}}^{\dagger }}[{{d}_{t}}\hat{U}]$. In comparison with $\hat{H}$ the shaking and the on-site interaction parts are removed by $\hat{U}$ while the tunneling part is dressed in oscillating terms. We aim at the situation when inter-band density dependent sp tunneling makes the dominant tunneling contribution. Accordingly, we assume the resonant condition

$$\begin{eqnarray}&&{{E}_{1}}+{{U}_{1}}=N\omega +2\Delta ,\end{eqnarray} \tag{ 2 }$$

where N is integer and $\Delta \ll \omega$ is the detuning and we time-average the Hamiltonian. In the process fast oscillating terms ∼1/ω are neglected. Necessarily the shaking frequency ω is chosen large compared to all the tunneling amplitudes. We obtain the effective Hamiltonian ${{H}_{{\rm eff}}}={{\tilde{H}}_{{\rm tun}}}+{{\tilde{H}}_{{\rm dit}}}+{{\tilde{H}}_{{\rm ons}}}$, with

$$\begin{eqnarray}\begin{array}{rcl} {{\tilde{H}}_{{\rm tun}}} & = & {{\tilde{J}}_{0}}\mathop{\sum }\limits_{\langle ij\rangle }\left[ \hat{s}_{i}^{\dagger }{{{\hat{s}}}_{j}}+\widehat{{{s}_{\uparrow }}}_{i}^{\dagger }{{\widehat{{{s}_{\uparrow }}}}_{j}} \right]+{{\tilde{J}}_{1}}\mathop{\sum }\limits_{\langle ij\rangle }\hat{p}_{i}^{\dagger }{{{\hat{p}}}_{j}} \\ {{\tilde{H}}_{{\rm dit}}} & = & \mathop{\sum }\limits_{\langle ij\rangle }\left[ {{\tilde{T}}_{0}}\widehat{{{s}_{\uparrow }}}_{i}^{\dagger }(\hat{n}_{i}^{{}}+\hat{n}_{j}^{{}}){{\widehat{{{s}_{\uparrow }}}}_{j}}+{{\tilde{T}}_{1}}\hat{p}_{i}^{\dagger }(\hat{n}_{i}^{\uparrow }+\hat{n}_{j}^{\uparrow }){{{\hat{p}}}_{j}} \right] \\ {} & {} & +\tilde{T}_{01}^{-}\mathop{\sum }\limits_{i}\left( \hat{s}_{2i}^{\dagger }\hat{p}_{2i-1}^{{}}+h.c. \right)-\tilde{T}_{01}^{+}\mathop{\sum }\limits_{i}\left( \hat{s}_{2i}^{\dagger }\hat{p}_{2i+1}^{{}}+h.c. \right) \\ {{\tilde{H}}_{{\rm ons}}} & = & \Delta \mathop{\sum }\limits_{\langle ij\rangle }\hat{p}_{i}^{\dagger }p_{i}^{{}}-\Delta \mathop{\sum }\limits_{\langle ij\rangle }\hat{s}_{i}^{\dagger }{{s}_{i}}. \\ \end{array}\end{eqnarray} \tag{ 3 }$$

The intra-band tunneling parts above are modified in the standard manner [18]: ${{\tilde{J}}_{l}}={{\mathcal{J}}_{0}}\left( K/\omega \right){{J}_{l}}$ (as well as ${{\tilde{T}}_{l}}={{\mathcal{J}}_{0}}\left( K/\omega \right){{T}_{l}}$) for $l\in \{0,1\}$, where ${{\mathcal{J}}_{0}}\left( K/\omega \right)$ is the ordinary Bessel function of order zero. In the case of the inter-band part, time averaging brings us, however, a new effect. The inter-band hopping is modified by the Bessel function of order N with different amplitudes depending on the direction of this process (+ or −): $\tilde{T}_{01}^{\pm }={{\mathcal{J}}_{N}}\left( {{A}^{\pm }}/\omega \right){{T}_{01}}$ where ${{A}^{\pm }}=\sqrt{{{(K\pm \delta {{E}_{1}}\;{\rm cos} \;\varphi )}^{2}}+{{K}^{2}}{{{\rm sin} }^{2}}\varphi }$. These amplitudes depend on the relative phase of the drivings, φ, which, we believe, can be controlled in real experiments with a good precision. The detuning, Δ, leads to residual on-site potential ${{\tilde{H}}_{{\rm ons}}}$.

Now we can tune the hopping parameters. The intra-band amplitudes may be made very small by choosing $K/\omega$ such that ${{\mathcal{J}}_{0}}(K/\omega )\approx 0$. For slightly different $K/\omega$ ss hopping remains negligible (so the composites may be still considered as immobile) while the typically much larger pp hopping start to plays a role and has to be taken into account. At the same time the inter-band hopping is large since it depends on Bessel functions of order $N\ne 0$.

From now on we set the recoil energy, ${{E}_{{\rm R}}}={{h}^{2}}/(8M{{a}^{2}})$, as an energy unit and consider the ground state structure of ${{H}_{{\rm eff}}}$ on an exemplary case of lattice depths ${{V}_{0}}=8$, ${{V}_{\bot }}=25$ and interaction strength $\alpha ={{a}_{s}}/a=-0.1$ (with a_s being the (negative) scattering length). We choose the vertical shaking to be in phase with the lateral one ($\varphi =0$), giving $\tilde{{\bf T}}_{{\bf 01}}^{-}\gt \tilde{{\bf T}}_{{\bf 01}}^{+}$. In the region with dominant inter-orbital tunneling, we expect that the ground state is given by the density wave configuration (DW) with every second side occupied by composites. The reason for that is quite simple. In such a configuration neighboring sites may contain excess fermions only and the sp tunneling lowers the energy for such a situation. On the other hand if two composites reside in the consecutive sites then the inter-band sp tunneling would break one of the composites (which costs the energy) while intraband pp tunneling for excess fermions is assumed small so it cannot lower the energy. To confirm that prediction, we first assume the composites to be immobilized due to the negligible value of ${{\tilde{{\bf J}}}_{{\bf 0}}}$ and ${{\tilde{{\bf T}}}_{{\bf 0}}}$ when compared to other hopping amplitudes in the whole regime under consideration. Then finding the ground state configuration boils down to solving the single particle Hamiltonian for a group of all possible configurations of $\hat{n}_{i}^{\uparrow }$ (that determines composite distribution). Within the approximation of immobilized $\uparrow$-fermions, we can replace their number operators by c-numbers $n_{i}^{\uparrow }=1,0$ depending on the presence or absence of the fermion on site i. Since the search space grows exponentially with the number of sites, it quickly becomes too large for the exact diagonalization. Thus we apply simulated annealing [25] to find the dependence of the lowest energy configuration on $K/\omega$. Our calculations have been performed for the lattice of 24, 40 and 60 sites and the obtained results do not depend on the number of sites.

The configurations obtained for $\Delta =0$ are shown in figure 1. We find indeed two possible configurations of composites: 1) clustered phase (CL) where the composites cluster together with the rest of the lattice being empty or 2) DW phase (the shaded region) where we have alternation of occupied and empty sites. We see that DW structures occur for the shaking parameter, $K/\omega$, for which $\mid \tilde{T}_{01}^{-}\mid +\mid \tilde{T}_{01}^{+}\mid \geqslant \mid {{J}_{1}}\mid {{\mathcal{J}}_{0}}(K/\omega )$. Simulations performed for other values of Δ gave similar dependence of the ground state configuration on $K/\Omega$ with the region of DW configuration broadening slightly for $\Delta \gt 0$. The obtained phase diagram is stable under small fluctuations in the fermion densities for both $\uparrow$- and $\downarrow$-fermions—the DW structure is preserved with some sparse defects appearing.

Zoom In Zoom Out Reset image size

When the DW-configuration of the composites minimizes the energy of the system, the dominant hopping process is the inter-band sp one and the effective Hamiltonian for the excess $\downarrow$-fermions corresponds to the Rice–Mele model [6],

$$\begin{eqnarray}\begin{array}{rcl} {{H}_{{\rm DW}}} & = & \tilde{T}_{01}^{-}\mathop{\sum }\limits_{i}\left( \hat{s}_{2i}^{\dagger }\hat{p}_{2i-1}^{{}}+h.c. \right)-\tilde{T}_{01}^{+}\mathop{\sum }\limits_{i}\left( \hat{s}_{2i}^{\dagger }\hat{p}_{2i+1}^{{}}+h.c. \right) \\ {} & {} & +\Delta \mathop{\sum }\limits_{\langle i\rangle }\hat{p}_{2i+1}^{\dagger }p_{2i+1}^{{}}-\Delta \mathop{\sum }\limits_{\langle i\rangle }\hat{s}_{2i}^{\dagger }{{s}_{2i}}. \\ \end{array}\end{eqnarray} \tag{ 4 }$$

From now we shall drop the $\widetilde{\;}$ sign over tunneling amplitudes as we shall consider effective tunnelings only restricting to (4). The above Hamiltonian describes a perfect lattice without defects. However, if the defects in the lattice are sparse with comparison to the edge mode length (see section 4) and we tune the shaking to make intra-band hopping small, then each of the domains may be separately described by the Rice–Mele Hamiltonian (4). Otherwise, a proper description of the system requires including also intra-band tunnelings.

To write the the Hamiltonian in the momentum space we specify a unit cell to contain two neighboring sites of which only one is always occupied by a composite. Such a unit cell can be chosen in two different ways depending whether the composite resides in the first (we call such a configuration 'PS') or the second site ('SP') of the open chain (see figure 2). As expected for the RM model [6] these two choices of the unit cell give rise to topologically distinct states. When written in the momentum space, the Hamiltonian (4) reads:

$$\begin{eqnarray}&&H_{{\rm DW}}^{\pm }(k)=\pm (T_{01}^{\mp }-T_{01}^{\pm }{\rm cos} (2ka)){{\sigma }_{x}}+T_{01}^{\pm }{\rm sin} (2ka){{\sigma }_{y}}-\Delta {{\sigma }_{z}},\end{eqnarray} \tag{ 5 }$$

where the $H_{{\rm DW}}^{+}(k)$ corresponds to the 'PS' configuration, $H_{{\rm DW}}^{-}(k)$ corresponds to the 'SP' one, and ${{\sigma }_{x,y,x}}$ are Pauli matrices. The dispersion relations are the same for both configurations:

$$\begin{eqnarray}&&{{\epsilon }_{\pm }}(k)=\pm \sqrt{{{\Delta }^{2}}+{{(T_{01}^{+}+T_{01}^{-})}^{2}}+4T_{01}^{+}T_{01}^{-}{{{\rm cos} }^{2}}(ka)}.\end{eqnarray} \tag{ 6 }$$

Zoom In Zoom Out Reset image size

Topologically distinct configurations are characterized by different Zak phases [2] (i.e. Berry phases acquired across the Brillouin zone). Zak phases of particular states depend on the choice of the unit cell, their difference forms an invariant of the system. The Zak phase is given by [2, 3]:

$$\begin{eqnarray}&&{{\phi }_{{\rm Zak}}}=i\int _{-\pi /2}^{\pi /2}\langle {{u}_{k}}\mid {{\partial }_{k}}\mid {{u}_{k}}\rangle {\rm d}k,\end{eqnarray} \tag{ 7 }$$

where $\mid {{u}_{k}}\rangle$ are Bloch functions of the system i.e. eigenfunctions of the Hamiltonian (5). For $\Delta =0$, when the Hamiltonian is equivalent to the SSH model [1], we obtain $\phi _{{\rm Zak}}^{{\rm SP}}-\phi _{{\rm Zak}}^{{\rm PS}}=\pi$. This indicates that 'SP' and 'PS' phases are topologically distinct—one of them must be nontrivial. For nonzero Δ we obtain fractional (in units of π) Zak phase differences changing from π to $1.67\pi$ for $\Delta \in [0,0.002\omega ]$. To determine which configuration has a nontrivial topology, we investigate the existence of edge modes.

In our model, defects arise naturally due to the emergent nature of the DW structure. As discussed in [20], the timescale required to reach a particular DW lattice configuration is set by the minority component tunneling rate. Subsequently the timescale to form the entire DW configuration is governed by the corresponding Lieb–Robinson bound [26]. When the time of creation is not sufficiently long, smaller regions of different DW configurations, separated by domain walls, will be created. Moreover, due to number fluctuations present for trapped atoms, the composites will not be exactly at half-filling. Any deviation from this filling will result in a defect in the form of a vacancy or a filled site.

Both kinds of impurities—domain walls and lattice defects—give rise to topologically protected localized modes [27] (visualized in figure 3). If we tune the shaking to make both ss and pp tunnelings negligibly small, then on domain walls we effectively create open boundary conditions. This will result in an appearance of edge modes in SP configuration when $T_{01}^{-}\gt T_{01}^{+}$ and in PS configuration otherwise. These modes vanish sharply on the edges (compare figures 3(a) and (c)) indicating that the configuration has a nontrivial topology. Let us here focus on the case when $T_{01}^{-}\gt T_{01}^{+}$. The edge modes have energies $\pm \Delta$ and in the continuum limit, their eigenvectors are given by the spinor: $({{\psi }_{s}}(x),{{\psi }_{p}}(x))$ Depending on which edge we are, setting x = 0 on the boundary, we get: the edge mode on the left end (the one ending with an S site) with energy $-\Delta$ where ${{\psi }_{s}}(x)=A({{e}^{-{{\lambda }^{+}}x}}-{{e}^{-{{\lambda }^{-}}x}})$, ${{\psi }_{p}}(x)=0$ and the edge mode on the right end (the one ending with a P site) with energy Δ, where ${{\psi }_{s}}(x)=0$, ${{\psi }_{p}}(x)=A({{e}^{{{\lambda }^{+}}x}}-{{e}^{{{\lambda }^{-}}x}})$ with

$$\begin{eqnarray}&&{{\lambda }^{\pm }}=\frac{T_{01}^{-}\pm \sqrt{T_{01}^{-}(2T_{01}^{+}-T_{01}^{-})}}{2T_{01}^{-}},\end{eqnarray} \tag{ 8 }$$

and A being the normalization constant. A detailed and tutorial discussion of edge modes in the dimer model is given in [28].

Zoom In Zoom Out Reset image size

When tuning $K/\omega$ further from the zero point of the Bessel function, the ss hopping is still negligible, but the pp hopping becomes significant. Therefore, on those boundaries that are separated by P-sites (figures 3(b) and (d)), particles can tunnel through the boundary and the mode vanishes exponentially on both sites giving topological solitons with energy $\Delta .$ Defects occurring inside the 'SP' configuration (figure 3(e)) give rise to two localized modes on both sides of the impurity. Depending on ss and pp tunneling rates they may end sharply on the boundary or smoothly vanish inside the defect. The width of the edge states depends on the hopping amplitudes and can be changed by tuning the value of $K/\omega$. For $K/\omega =2.3$ the edge state is about 15 lattice sites long and it becomes narrower with higher values of $K/\omega$.

Defects present in the lattice are associated with local changes of fermion number by fraction, ${{N}_{{\rm frac}}}={{f}_{0}}/\pi$ where ${{f}_{0}}={{{\rm tan} }^{-1}}\left[ \mid T_{01}^{+}-T_{01}^{-}\mid /2\Delta \right]$ [6] at zero temperature. At finite temperatures, $T\ll \Delta$, the corresponding fractional fermion number for the localized mode is given by the thermal expectation value [29, 30], $N_{{\rm frac}}^{T}={{f}_{0}}/\pi -{\rm sgn} ({{f}_{0}}){\rm exp} [-\Delta /T]$, where ${\rm sgn} ()$ is the sign function. For a typical value of $\Delta =0.01{{E}_{{\rm R}}}$ the relevant temperature is in the nano-Kelvin regime for $^{40}K$.

Finally let us briefly comment on the dynamics of solitons. The solitonic localized modes are pinned to the defects. Their dynamics is affected by tunnelings as well as thermal excitations at finite temperature. A single density dependent tunneling event [given by the T₀ term in the original Hamiltonian (1)] will change the 'PS' pair to 'SP'. This corresponds to a motion of the defect and thus the motion of the soliton localized on the edge (e.g. compare figures 3(a) or (b)). While by adjusting the frequency we minimize the influence of such processes they will be still partially present due to, e.g., a frequency mismatch with respect to the exact zero of the appropriate ${{\mathcal{J}}_{0}}(K/\omega )$ Bessel function or higher-order terms, discussed in the next section. One may envision also that once the system is formed, the frequency/amplitude of shaking is changed a little to stimulate the motion of defects.

Consider, however, the situation represented in figure 3(e). The single 'PS' to 'SP' tunneling will remove two defects, that would correspond to collisional annihilation of two localized modes. The stable 'solitonic' solution corresponds thus to situations with well-separated defects that cannot easily be removed by tunnelings.

Let us first discuss the timescales needed to realize the described system experimentally. The time for the formation of the crystal is determined by the tunneling rate of $\uparrow$-fermions and bounded from above by loss rates. With an example of $^{40}K$ we estimate the tunneling time to be of the order of 10 ms (see figure 1) except at the close vicinity of $K/\omega =2.4$. Close to $K/\omega =2.4$ when the tunneling in the averaged Hamiltonian vanishes there will be still residual higher-order tunneling. Its non-resonant effect can be estimated in an analogous way to Bloch–Siegert shifts in quantum optics [31] being proportional to $J_{0}^{2}/4\omega$ which again for $^{40}K$ is about 10 ms. The formation of the crystal will take a few tunneling times. The precise estimate would require full dynamical calculation of the crystal formation which is beyond the scope of the present paper.

The shorter time will lead to numerous defects. The number of defects is dependent also on the temperature. For the nano-Kelvin regime and assuming a sufficient preparation time, the number of defects (which can be calculated [32] comparing the energy of creating a defect with temperature assuming a canonical ensemble) will be of the order of 1%—which allows one to obtain well-separated edge states. With the system prepared, we may probe its topological properties. The Zak phase can be measured experimentally in a way that has been proposed in [7], with application of coherent Bloch oscillations combined with Ramsey interferometry. At exact half filling of excess fermions the topologically nontrivial band will be filled and the standard time-of-flight method cannot detect the Zak phase⁵ . Slightly lower filling of excess fermions does not affect the DW structure of composites (so the Rice–Mele model is applicable). Then such a measurement of the Zak phase should be possible. For this measurement also the number of defects should not be too large. However, to detect the localized states, defects can be helpful. Localized states can be observed with photo-emission spectroscopy [33] with time-of-flight where the number of edge modes in our system results in an increased peak intensity near zero momentum making the signal less susceptible to noise. At half filling, the only localized states that can be occupied are those of the negative energy. Fermion number fractionalization can be probed on defects with application of single site imaging.

We have shown that a combination of shaking and attractive interactions in a one-dimensional optical lattice can give rise to a topologically nontrivial system. We have used standard lateral shaking but also introduced an additional vertical shaking. Together, they result in a dimerized tunneling structure. Moreover, by tuning the onsite energy slightly out of the resonance we can induce the staggered potential. By controlling the filling, we have shown further the presence of topologically protected localized modes. We believe that such modes can be experimentally verified at accessible temperatures.

Numerical calculations were performed at the Academic Computer Center in Gdańsk. This work was realized under National Science Center (Poland) project no. DEC-2012/04/A/ST2/00088. AP is also supported by the International PhD Project 'Physics of future quantum-based information technologies', grant MPD/2009-3/4 from Foundation for Polish Science and by the University of Gdansk grant BW 538-5400-B524-14.

We derive the minimal model in a standard manner starting from a many body Hamiltonian of a dilute gas of atoms in a second quantization representation [24, 34]. We consider two species (denoted by $\uparrow$-fermions and $\downarrow$-fermions) of equal masses which can occupy the lowest band. The $\downarrow$-fermions have an occupation close to unity, for them we consider also the excited, p orbital. Different species undergo contact interactions. The parameters in the Hamiltonian (1) in the main text are given by integrals of Wannier functions $\mathcal{W}_{i}^{0(1)}(x,y)$ on s(p)-bands, where i is a site index.

Specifically the single particle ss and pp hoppings do not depend on the type of species and read

$$\begin{eqnarray}\begin{array}{rcl} {{J}_{0}} & = & \int {{\left[ \mathcal{W}_{i}^{0}(x) \right]}^{*}}{{H}_{{\rm latt}}}\mathcal{W}_{i+1}^{0}(x){\rm d}x, \\ {{J}_{1}} & = & \int {{\left[ \mathcal{W}_{i}^{1}(x) \right]}^{*}}{{H}_{{\rm latt}}}\mathcal{W}_{i+1}^{1}(x){\rm d}x, \\ \end{array}\end{eqnarray} \tag{ A.1 }$$

where ${{H}_{{\rm latt}}}=-\frac{{{\partial }^{2}}}{\partial {{x}^{2}}}+{{V}_{0}}{{{\rm sin} }^{2}}(\pi x/a)$ is a single particle Hamiltonian for a static lattice. Observe the lack of two in the kinetic energy as we work in recoil units. The contact interactions between different species lead to density induced tunnelings [21–24]. The corresponding part of the Hamiltonian may be expressed as

$$\begin{eqnarray}\begin{array}{rcl} {{{\hat{H}}}_{{\rm dit}}} & = & \mathop{\sum }\limits_{\langle ij\rangle }\left[ {{T}_{01}}((j-i)\hat{p}_{i}^{\dagger }\hat{n}_{i}^{\uparrow }\hat{s}_{j}^{{}}+h.c)+{{T}_{1}}\hat{p}_{i}^{\dagger }(\hat{n}_{i}^{\uparrow }+\hat{n}_{j}^{\uparrow }){{{\hat{p}}}_{j}} \right. \\ {} & {} & +\left. T_{1}^{{\prime} }\widehat{{{s}_{\uparrow }}}_{i}^{\dagger }(\hat{n}_{i}^{p}+\hat{n}_{j}^{p}){{\widehat{{{s}_{\uparrow }}}}_{j}}+{{T}_{0}}\widehat{{{s}_{\uparrow }}}_{i}^{\dagger }(\hat{n}_{i}^{{}}+\hat{n}_{j}^{{}}){{\widehat{{{s}_{\uparrow }}}}_{j}}+T_{0}^{{\prime} }\hat{s}_{i}^{\dagger }(\hat{n}_{i}^{\uparrow }+\hat{n}_{j}^{\uparrow }){{{\hat{s}}}_{j}} \right], \\ \end{array}\end{eqnarray} \tag{ A.2 }$$

where, let us recall, $\hat{s}_{i}^{\dagger },\hat{s}_{i}^{{}}$, $\hat{p}_{i}^{\dagger },\hat{p}_{i}^{{}}$ are the creation and annihilation operators of the $\downarrow$-fermions in the s- and p-bands respectively, while $\hat{s}_{\uparrow i}^{\dagger },\hat{s}_{\uparrow i}^{{}}$ are s-band creation and annihilation operators for the $\uparrow$-fermion. ${{\hat{n}}_{i}},\hat{n}_{i}^{p}$ and $\hat{n}_{i}^{\uparrow }$ denote the corresponding number operators. Throughout the paper we assume that minority $\uparrow$-fermions appear in pairs only due to strong attractive interactions. Thus some of the processes included above vanish. In particular, the term proportional to $T_{0}^{{\prime} }$ should be excluded as occupation of the i-site by a $\uparrow$-fermion means that there is a $\downarrow$-fermion occupying this site already, so the Pauli principle inhibits tunneling into this site. For a different reason $T_{1}^{{\prime} }$ may also be omitted as for the ground state we focus on the occupation of the i-site by a p fermion is possible energetically only if there is a composite there. The presence of a composite prohibits tunneling into this site of an s-fermion. There is a possibility that a site occupied by a composite and a p fermion neigbors a site with two fermions: one in the s and the other in the p band. Such sites can exchange s-type fermions by the $T_{1}^{{\prime} }$ process. One should keep in mind that p-fermions appear in the system only due to resonant shaking (otherwise they cost E₁)—their presence in neighboring sites should be a rare event. The remaining terms form ${{\hat{H}}_{{\rm dit}}}$ included in the Hamiltonian (1) of the paper.

The amplitudes, T, are given by integrals over four Wannier functions and take the form

$$\begin{eqnarray}\begin{array}{rcl} {{T}_{0}} & = & {{g}_{1D}}(\alpha )\int {{\left[ \mathcal{W}_{i}^{0}(x) \right]}^{*}}{\left|\mathcal{W}_{i}^{0}{(x)}\right|}^{2}\mathcal{W}_{i+1}^{0}(x){\rm d}x, \\ {{T}_{1}} & = & {{g}_{1D}}(\alpha )\int {{\left[ \mathcal{W}_{i}^{1}(x) \right]}^{*}}{\left|\mathcal{W}_{i}^{0}{(x)}\right|}^{2}\mathcal{W}_{i+1}^{1}(x){\rm d}x, \\ {{T}_{01}} & = & {{g}_{1D}}(\alpha )\int {{\left[ \mathcal{W}_{i}^{1}(x) \right]}^{*}}{\left|\mathcal{W}_{i}^{0}{(x)}\right|}^{2}\mathcal{W}_{i+1}^{0}(x){\rm d}x, \\ \end{array}\end{eqnarray} \tag{ A.3 }$$

where ${{g}_{1D}}(\alpha )$ is a renormalized 1D coupling constant [35] and $\alpha ={{a}_{s}}/a$ is the ratio of the interaction strength to the lattice spacing.

A pictorial representation of different tunneling processes is shown in figure A1 . Note that, since $\uparrow$-fermions are minority fermions, they are always paired and the probability of their tunneling to the p-band is negligibly small. That is why in the case of these fermions we consider only ss tunneling. On the other hand, the presence of $\uparrow$-fermions (and therefore composites) stimulates pp tunneling of $\downarrow$-fermions. Observe that the first term in (A.2) is the inter-band sp hopping which has a staggered nature (reflected by the $(j-i)$ sign) and may happen when a composite-empty site adjoins a composite-occupied one.

Zoom In Zoom Out Reset image size

Let us now discuss the on-site energies present in the Hamiltonian. The corresponding term reads:

$$\begin{eqnarray}&&{{\hat{H}}_{\operatorname{int}}}=U_{0}^{{}}\mathop{\sum }\limits_{i}\hat{n}_{i}^{\uparrow }{{\hat{n}}_{i}}+U_{1}^{{}}\mathop{\sum }\limits_{i}\hat{p}_{i}^{\dagger }{{\hat{p}}_{i}}\hat{n}_{i}^{\uparrow }+{{E}_{1}}\mathop{\sum }\limits_{i}\hat{p}_{i}^{\dagger }{{\hat{p}}_{i}}.\end{eqnarray} \tag{ A.4 }$$

U₀, U₁ are given by:

$$\begin{eqnarray}&&{{U}_{0}}={{g}_{1D}}(\alpha )\int \mathcal{W}_{i}^{0}{{(x)}^{4}}{\rm d}x,\end{eqnarray} \tag{ A.5 }$$

$$\begin{eqnarray}&&{{U}_{1}}={{g}_{1D}}(\alpha )\int \mathcal{W}_{i}^{0}{{(x)}^{2}}\mathcal{W}_{i}^{1}{{(x)}^{2}}{\rm d}x.\end{eqnarray} \tag{ A.6 }$$

U₀ is by far the biggest (on the modulus) energy scale and is responsible for pairing. We assume that composites are formed in the s band only. The composites could form also in the p band with (negative) energy U₀^p given by (A.5) with $\left|\mathcal{W}_{i}^{1}(x)\right|$ instead of $\left|\mathcal{W}_{i}^{0}(x)\right|$. Since p-orbitals are extended by comparison with s-functions, $U_{0}^{p}\approx 0.6{{U}_{0}}$ for typical lattice depths (e.g. at ${{V}_{0}}=8{{E}_{{\rm R}}}$ as assumed in numerical calculations). Similarly ${{U}_{1}}\approx 0.4{{U}_{0}}$. Single particle energy of occupying p-band E₁ reads:

$$\begin{eqnarray}&&{{E}_{1}}=\frac{1}{2}\int {{\left[ \mathcal{W}_{i}^{1}(x) \right]}^{*}}{{H}_{{\rm latt}}}\mathcal{W}_{i}^{1}(x){\rm d}x-\frac{1}{2}\int {{\left[ \mathcal{W}_{i}^{0}(x) \right]}^{*}}{{H}_{{\rm latt}}}\mathcal{W}_{i}^{0}(x){\rm d}x,\end{eqnarray} \tag{ A.7 }$$

with the origin of the energy axis corresponding to the s-fermion single particle energy. E₁ may be larger than $\mid {{U}_{0}}\mid$.

Consider now the effects due to lateral and vertical shaking. The former is quite standard [18] and leads to a familiar term

$$\begin{eqnarray*}&&K{\rm cos} \omega t\mathop{\sum }\limits_{j}j(\hat{n}_{j}^{\uparrow }+\hat{s}_{j}^{\dagger }{{\hat{s}}_{j}}+\hat{p}_{j}^{\dagger }{{\hat{p}}_{j}}),\end{eqnarray*}$$

where K is the shaking amplitude. The vertical shaking of the lattice depth (assumed to be not too large) causes periodic changes of single particle hoppings ${{J}_{z}}(t)={{J}_{z}}+\delta {{J}_{z}}{\rm cos} \omega t$ for z = 0,1, with amplitudes:

$$\begin{eqnarray}&&\delta {{J}_{z}}=\int {{\left[ \mathcal{W}_{i}^{z}(x) \right]}^{*}}\left( \delta {{V}_{0}}{{{\rm sin} }^{2}}\frac{\pi x}{a} \right)\mathcal{W}_{i+1}^{z}(x){\rm d}x.\end{eqnarray} \tag{ A.8 }$$

On time averaging we will see that these periodic changes have negligibly small influence on the system and can be omitted. That is the reason why they do not appear in the Hamiltonian (1) of the main text.

Next we have periodic changes in the onsite energy with amplitudes:

$$\begin{eqnarray}\begin{array}{rcl} \delta {{E}_{1}} & = & \frac{1}{2}\int {{\left[ \mathcal{W}_{i}^{1}(x) \right]}^{*}}\left( \delta {{V}_{0}}{{{\rm sin} }^{2}}\frac{\pi x}{a} \right)\mathcal{W}_{i}^{1}(x){\rm d}x \\ {} & {} & -\frac{1}{2}\int {{\left[ \mathcal{W}_{i}^{0}(x) \right]}^{*}}\left( \delta {{V}_{0}}{{{\rm sin} }^{2}}\frac{\pi x}{a} \right)\mathcal{W}_{i}^{0}(x){\rm d}x. \\ \end{array}\end{eqnarray} \tag{ A.9 }$$

Contrary to changes in the tunneling, this onsite effect is very important for the model and allows us to realize tunneling dimerization of the RM model.

The standard time-averaging procedure can be obtained applying Floquet theorem following [36] and deriving the effective Hamiltonian via repeated commutation of the time independent Hamiltonian with operator $F(t)=-\int _{0}^{t}{\rm d}t^{\prime} {{H}_{{\rm sh}}}(t)$. Using this approach one can verify that, already in the second commutator, terms containing periodic changes in hopping parameters become negligibly small; which allows us to omit them from our considerations.

As mentioned in the main text we invoke instead a time-dependent unitary transformation $\hat{U}={\rm exp} [-{\rm i}{{\hat{H}}_{\operatorname{int}}}t-{\rm i}\int _{0}^{t}{{\hat{H}}_{{\rm sh}}}(t^{\prime} ){\rm d}t^{\prime} ]$. We include in the transformation also the on-site terms with the aim of locating resonant coupling between bands. We obtain the transformed Hamiltonian, $\hat{H}^{\prime} ={{\hat{U}}^{\dagger }}\hat{H}\hat{U}-{\rm i}{{\hat{U}}^{\dagger }}[{{d}_{t}}\hat{U}]$ in the form:

$$\begin{eqnarray} \hat{H}^{\prime} & = & \left[ {\rm exp} \left( {\rm i}\frac{K}{\omega }{\rm sin} \omega t \right)\mathop{\sum }\limits_{j}{\rm exp} \left( {\rm i}{{U}_{0}}(n_{j}^{\uparrow }-n_{j+1}^{\uparrow })t \right)\left( {{J}_{0}}\left[ \hat{s}_{j+1}^{\dagger }{{{\hat{s}}}_{j}}+\widehat{{{s}_{\uparrow }}}_{j+1}^{\dagger }{{\widehat{{{s}_{\uparrow }}}}_{j}} \right] \right. \right. \\ {} & {} & +\left. \left. {{T}_{0}}\widehat{{{s}_{\uparrow }}}_{j+1}^{\dagger }(\hat{n}_{j+1}^{{}}+\hat{n}_{j}^{{}}){{\widehat{{{s}_{\uparrow }}}}_{j}} \right)+h.c. \right] \\ {} & {} & +\left[ {\rm exp} \left( {\rm i}\frac{K}{\omega }{\rm sin} \omega t \right)\mathop{\sum }\limits_{j}{\rm exp} \left( {\rm i}{{U}_{1}}(n_{j}^{\uparrow }-n_{j+1}^{\uparrow })t \right)\left( {{J}_{1}}\hat{p}_{j+1}^{\dagger }{{{\hat{p}}}_{j}}+{{T}_{1}}\hat{p}_{j+1}^{\dagger }(\hat{n}_{j+1}^{\uparrow }+\hat{n}_{j}^{\uparrow }){{{\hat{p}}}_{j}} \right)+h.c. \right] \\ {} & {} & +{{T}_{01}}\mathop{\sum }\limits_{\langle ij\rangle }(j-i){\rm exp} \left[ {\rm i}({{E}_{1}}+{{U}_{1}}\hat{n}_{i}^{\uparrow })t+{\rm i}\left( (i-j)\frac{K}{\omega }+\frac{\delta {{E}_{1}}}{\omega } \right){\rm sin} \omega t \right]\hat{p}_{i}^{\dagger }\hat{n}_{i}^{\uparrow }\hat{s}_{j}^{{}} \\ {} & {} & +{{T}_{01}}\mathop{\sum }\limits_{\langle ij\rangle }(j-i){\rm exp} \left[ -{\rm i}({{E}_{1}}+{{U}_{1}}\hat{n}_{i}^{\uparrow })t-{\rm i}\left( (i-j)\frac{K}{\omega }+\frac{\delta {{E}_{1}}}{\omega } \right){\rm sin} \omega t \right]\hat{s}_{j}^{\dagger }\hat{n}_{i}^{\uparrow }\hat{p}_{i}^{{}}. \\ \end{eqnarray} \tag{ B.1 }$$

Now we assume the resonant condition

$${{E}_{1}}+{{U}_{1}}=N\omega +2\Delta , \tag{ B.2 }$$

where N is integer and $\Delta \ll \omega$ is the detuning. The shaking frequency ω chosen is large compared to all the tunneling amplitudes. Before standard time averaging one more simplification is made. We consider a low energy Hilbert subspace, where due to strong attractive interactions all $\uparrow$-fermions are paired. The resonant condition (B.2) may be fulfilled only for sites occupied by the composites, i.e with $n_{i}^{\uparrow }=\langle \hat{n}_{i}^{\uparrow }\rangle =1$. We average the Hamiltonian over the oscillation period and neglect terms (∼1/ω) obtaining ${{H}_{{\rm eff}}}$ (3). Let us note also that on-site direct excitation of the p-band due to periodic shaking (see e.g. [37]) is negligible in our model due to the resonance condition (B.2) involving composite binding energy.