Introduction

Topological phases of matter are among the most exciting developments of modern condensed matter physics,1,2,3,4 owing to their rich phenomenology and wide-ranging potential applications from metrology5 to quantum computation.6 Many experimental platforms have been used to realize these exotic phases of matter such as cold atoms,7 photonic lattices,8,9 and engineered solid-state systems including graphene nanoribbons,10,11 arrays of carbon monoxide molecules,12,13 and chlorine monolayers14 on a copper surface. The band theory of topological insulators (TIs)15 based on the independent-electron approximation is well developed and has had many successes. However, in many of the possible experimental platforms for quantum simulation of TIs using electrons in solids, such as dopant atoms and gate-defined quantum dots in semiconductors,16,17 the electron–electron interaction is much stronger than the hopping amplitude of the electrons18,19 and therefore the independent-electron approximation is poor. Topological phases of strongly correlated models form a topic of ongoing active research with intense theoretical20 and experimental effort, including recent implementations in cold atoms21 and two-dimensional materials.22 There have been various proposals for the equivalent of the single-particle Berry phase (or Zak phase in one dimension) for the characterization of interacting topological phases, from the magnetic-flux-induced Berry phase23,24 to Green’s functions25 and entanglement.26,27,28

Here we discuss one of the simplest one-dimensional (1D) models of strongly correlated TIs, the Su–Schrieffer–Heeger–Hubbard (SSHH) model, whose topological properties in various contexts have been investigated using the entanglement entropy,26,27 the entanglement spectrum,28 correlation functions,29 quench dynamics,30 and Berry phase.31 The SSHH model describes electrons hopping on a 1D superlattice with staggered hopping amplitudes but uniform local interaction. In this model, there exists a charge excitation gap at half-filling due to the on-site interaction (the Mott gap) and another gap at quarter-filling due to dimerization. This opens the possibility of realizing these fillings in experiments, for example by measuring transport while varying the chemical potential and looking for vanishing conductance when the chemical potential lies in the gaps.18 For this reason, we focus on these two fillings.

We introduce the concept of the reduced many-body Zak phase based on the reduced density matrix of a subsystem and show that this phase, rather than the normal many-body Zak phase of the full system, should be used for classifying the topological phases at half-filling. This phase jumps from 0 to π as the hopping amplitude difference between the even and odd sites changes sign. At half-filling, the usual bulk–edge correspondence is manifested in the topological phase transition: the closing and reopening of the eigenenergy gap at the transition point accompanies the appearance of uncorrelated edge states. This is evident in the triplon-excitation spectrum of the dimer chain. In contrast, at quarter-filling the edges remain correlated to the bulk for both signs of the hopping amplitude difference, because of the presence of a long-range antiferromagnetic (AFM) order. There is also no gap in the eigenenergy spectrum due to the presence of gapless spin excitations. So the quarter-filled state does not show the characteristics of a TI. However, we show that applying an external magnetic field leads to a transition to the topological ground state of the non-interacting Su–Schrieffer–Heeger (SSH) model. Importantly, the strong on-site interaction significantly reduces the critical field strength required for the transition. Thus our analysis paves the way for the observation of electronic 1D topological insulator states in nanofabricated semiconductor devices. We propose a device architecture for observing this transition in a 1D chain of dopant atoms or quantum dots. The transition can be probed by measuring the tunneling current through the edges of the chain, which we estimate using a many-body formulation for the conductance of coupled quantum dots.32,33,34

Results

The SSHH model

The SSHH Hamiltonian is

$$H={H}_{0}+V,$$
(1)

where

$${H}_{0}=\sum _{j,\sigma =\uparrow \downarrow }-\left[1+{(-1)}^{j}\Delta t\right]{c}_{j+1,\sigma }^{\dagger }{c}_{j,\sigma }+{\rm{h.c.}},$$

and

$$V=U\sum _{j}{n}_{j,\uparrow }{n}_{j,\downarrow }.$$

H0 is the well-known non-interacting SSH model35 of a particle hopping along a chain with staggered hopping amplitudes, t± = 1 ± Δt, as shown in Fig. 1a, and V is the on-site interaction. Here \({c}_{j,\sigma }^{\dagger }\) denotes the creation operator for the particle at site j and spin σ. All energies in this paper are scaled by the mean value of the two hopping amplitudes.

Fig. 1: Energy bands and edge states of the SSH model.
figure 1

a The non-interacting SSH model of particles hopping along a 1D chain of potential wells with alternating hopping amplitudes. b Particle distribution of a typical bulk state (diamond) and an edge state (circle) for Δt = 0.5. Inset: Eigenenergies calculated for a chain with N = 20 sites; the bulk energy gap closes at Δt = 0 and opens again with the appearance of the mid-gap edge states.

Full size image

We first describe briefly the topological phases of the SSH model given by H0.36,37 For 1D periodic systems of independent particles, the Berry phase picked up during an adiabatic process when the particle moves across the Bloch states in the Brillouin zone, first discussed by Zak,38 is

$$\phi =i{\int _{-\pi /d}^{\pi /d}}dk\left\langle {u}_{k}\right|{\partial }_{k}\left|{u}_{k}\right\rangle ,$$
(2)

where uk is the periodic part of the Bloch wavefunction, k the crystal momentum, and d the length of the unit cell.

In a chain with open boundary conditions (OBC) the single-particle eigenstates of this Hamiltonian consist of two distinct types: zero-energy edge states that are localized at the left and right edges, and extended bulk states that avoid the edges. The energy spectrum of the bulk states under periodic boundary conditions (PBC) splits into two bands, \({E_\pm}(k)={\pm {[2(1+\Delta {t}^{2})+2(1-\Delta {t}^{2})\cos (kd)]}^{1/2}}\). The bulk state wavefunction in PBC has the Bloch form

$${\psi }_{k}=\sum _{j}{e}^{ijkd/2}{u}_{k}(j)=\sum _{j}{e}^{ijkd/2}{e}^{i{\theta }_{j}(k)}{c}_{j}^{\dagger }\left|\varnothing \right\rangle ,$$
(3)

where \(\left|\varnothing \right\rangle\) is the vacuum and the phase shift θj(k) = 0 for the odd sub-lattice and \(| {E_\pm} (k)| \exp [i{\theta }_{j}(k)]=2\cos (kd/2)+2i\Delta t\sin (kd/2)\) for the even sub-lattice. For the SSH model, the Zak phase is quantized, more specifically it can only be 0 or π depending on the sign of Δt.7,36,37 There is a topological phase transition from the trivial phase (Δt < 0), where there is no edge state, to the non-trivial phase (Δt > 0) where the edge states appear. The energy levels of the bulk states form two bands separated by a gap in both phases, and the energy levels of the edge states appear in the middle of the band gap in the non-trivial phase (see Fig. 1b). The trivial and the non-trivial phases are characterized by the Zak phase of 0 and π, respectively.36

Charge excitation gap

With interaction the single-particle picture is no longer valid, but insight into the topological phases can be gained from looking at the addition energy spectrum, also known as the charge excitation spectrum, Ead(n) = E0(n) − E0(n − 1) where E0(n) is the many-body ground energy for filling n. We use exact diagonalization based on the Lanczos algorithm for an open chain with N = 12 sites to compute the addition energy spectrum and show it in Fig. 2. In this paper, we set N = 12 in all the numerical computations for the correlated case. Features that survive in the thermodynamic limit are either obvious or stated explicitly. At U = 0, the addition energy spectrum reduces to the single particle spectrum of the SSH model with the zero-energy edge state in the middle of the gap at half-filling. At large interaction, the Mott gap forms at half-filling separating the lower and upper Hubbard bands as expected of the Hubbard model. The spectrum has reflection symmetry through the middle point of the Mott gap due to the particle–hole symmetry. Interestingly, there are further gaps at one-quarter- and three-quarter-fillings, and in the non-trivial phase the edge states of the charge excitation cross to lie in these gaps (see Fig. 2a, d). A more detailed description of these edge states is given below in the discussion of the quarter-filled system. The formation of the quarter-filling gap is due to the combination of the on-site repulsion U and the hopping amplitude difference Δt and has been studied previously.39 Our numerical analysis shows that this gap is approximately Δt in the large-U limit.

Fig. 2: Addition energy spectrum and edge population of the SSHH model.
figure 2

The addition energy spectrum of the Su–Schrieffer–Heeger–Hubbard model as a function of interaction strengths in the a non-trivial phase and b trivial phase and c as a function of hopping amplitude difference in the strongly correlated limit. To keep track of the edge states, the energy levels are colored according to the population at the two ends of the chain (nedge); changes of color thus correspond to occupation of edge states. At strong interaction, there are three gaps formed, at half-filling and at quarter-/three-quarter-filling. The edge population shows that in the non-trivial phase the mid-gap edge states shift from half-filling at U = 0 to quarter-filling and three-quarter-fillings at large U [marked by the horizontal arrow in a]. The edge population as a function of filling in the strongly correlated limit is shown in d: In the trivial phase, the edge population increases gradually, while in the non-trivial phase the edge population increases sharply at quarter- and three-quarter-fillings, indicating the existence of an edge state for the charge excitation at these fillings.

Full size image

The charge gap and the mid-gap edge state in the addition energy spectrum at quarter-filling can be explained analytically in the full dimerization limit where t− = 0 and t+ > 0. The ground state energy level of each Hubbard dimer (coupled by t+) when there is a single particle is the bonding state E1 = −t+, while the ground state energy level of the same dimer with two particles is

$${E}_{2}=(1/2)\left[U-\sqrt{{U}^{2}+16{t}_{+}^{2}}\right].$$
(4)

The energy of a particle at the isolated edges is zero. As the particles are added into the chain, they first fill the dimer bonding states as they are lower in energy. When each dimer is filled with one particle, we reach the point of quarter-filling. Now if another particle is added to a dimer, the energy cost is ΔE = E2 − E1, while if a particle is added to the edges the energy cost is zero. Hence, if ΔE > 0 the edges are filled first, otherwise the dimers get filled with two particles until the point of half-filling and only then the edges are filled. Thus the transition of the addition energy level of the edge states from half-filling to quarter-filling happens at the critical interaction Uc such that ΔE = 0, or Uc = 3t+. It is obvious from the above discussion that the addition energy gap at quarter-filling is ΔE ≈ t+ in the strongly correlated limit where U ≫ t+. We note that in the general case where neither hopping amplitude is zero the gaps at half-filling and quarter-filling remain in the thermodynamic limit.39 We provide further evidence by a finite-size scaling analysis of the gap in Supplementary Material.

The reduced many-body Zak phase

For correlated systems one can define a many-body Zak phase, first introduced as a measure of macroscopic polarization,40,41,42 from the ground state of H satisfying a twisted boundary condition

$${\Psi }_{\kappa }({x}_{1},...,{x}_{j}+L,...,{x}_{N})={e}^{i\kappa L}{\Psi }_{\kappa }({x}_{1},...,{x}_{j},...,{x}_{N}),$$
(5)

where L is the total length of the chain. Writing \({\Psi }_{\kappa }={e}^{i\kappa \sum _{j=1}^{n}{x}_{j}}{\Phi }_{\kappa }\) with n the number of particles, then Φκ is the ground state of \(H(\kappa )={e}^{-i\kappa \sum _{j=1}^{n}{x}_{j}}H{e}^{i\kappa \sum _{j=1}^{n}{x}_{j}}\) that satisfies the periodic boundary condition in all coordinates xj. It can be seen as the many-body analog of uk in the non-interacting case. H(κ) is the Hamiltonian of a ring threaded with the magnetic flux κL,43 and can be obtained from H by the simple replacement \({t}_{j}\to {t}_{j}{e}^{-i\kappa ({x}_{j+1}-{x}_{j})}\). The many-body Zak phase is then defined as the adiabatic phase picked up by the many-body ground state when this magnetic flux is changed by one flux quantum

$$\phi =i{{\int _{-\pi /L}^{\pi /L}}}d\kappa \left\langle {\Phi }_{\kappa }\right|{\partial }_{\kappa }\left|{\Phi }_{\kappa }\right\rangle .$$
(6)

From Eq. (5), we see that Ψ−π∕L = Ψπ∕L since they satisfy the same anti-periodic boundary condition, hence the function Φκ at the initial and the end points are related by Φπ∕L = WΦ−π∕L where \(W={e}^{-i(2\pi /L)\sum _{j=1}^{n}{x}_{j}}\). For numerical computation, κ is discretized in a grid of M points κl from −π∕L to π∕L, and it can be shown that ϕ is simply the phase of a complex number23

$$\phi =\arg (Z)=\arg \left(\prod _{l=1}^{M-1}\left\langle {\Phi }_{{\kappa }_{l}}| {\Phi }_{{\kappa }_{l+1}}\right\rangle \right).$$
(7)

This phase can be rewritten in terms of the density matrix \(\rho (\kappa )=\left|{\Phi }_{\kappa }\right\rangle \left\langle {\Phi }_{\kappa }\right|\) as

$$\phi =\arg \left[{\rm{tr}}\left(W\prod _{l=1}^{M-1}\rho ({\kappa }_{l})\right)\right].$$
(8)

The SSHH Hamiltonian H has inversion symmetry, thus changing xj to −xj maps H(κ) to H(−κ). This implies that \(H(\kappa )={\mathcal{U}}H(-\kappa ){{\mathcal{U}}}^{\dagger }\) where \({\mathcal{U}}\) is the unitary operator of inversion. As a result E(κ) = E(−κ), and \({\Phi }_{\kappa }={e}^{i{\alpha }_{\kappa }}{\mathcal{U}}{\Phi }_{-\kappa }\), where ακ is an arbitrary phase. It follows that \({\rho }_{\kappa }={\mathcal{U}}{\rho }_{-\kappa }{{\mathcal{U}}}^{\dagger }\), and with a grid centered around κ = 0 such that κl = −κM−l, we have \({Z}^{* }={\rm{tr}}\left(\prod _{l=M-1}^{1}\rho ({\kappa }_{l}){W}^{\dagger }\right)={\rm{tr}}\left(\prod _{l=1}^{M-1}\rho ({\kappa }_{l}){{\mathcal{U}}}^{\dagger }{W}^{\dagger }{\mathcal{U}}\right)\). As the inversion operation transforms xj to −xj, we have \({{\mathcal{U}}}^{\dagger }{W}^{\dagger }{\mathcal{U}}=W\) and thus Z* = Z, or Z is real, meaning the many-body Zak phase must be either 0 or π depending on whether Z is positive or negative.

The total particle number n is a good quantum number for the eigenstates of the SSHH Hamiltonian. We first study the half-filled spinful case (n = N). The ground state Φ(κ) of H(κ) is computed with PBC and the Zak phase is obtained using the discrete formula of Eq. (7). We carry out the computation for −0.5 ≤ Δt ≤ 0.5 to see whether there is a topological phase change when Δt changes sign as in the non-interacting SSH model and for 0 ≤ U ≤ 10 to study the effect of interaction. The phase is found to be 0 for all values of U and Δt, so it does not reveal any phase transition for either weak or strong interaction. This result is expected at zero interaction, as we then have two copies of the SSH model, one with spin up and the other with spin down. For each copy the Zak phase changes from 0 to π, and it is straightforward that the Zak phase of the joint state (given by a Slater determinant) is the sum of the individual phases, hence its values are 0 and 2π, which are equal since the phase is defined only modulo 2π. Our result shows that adding interaction does not alter the Zak phase.

For revealing the topological phase transition, we introduce the reduced Zak phase of a subsystem within a larger correlated system. It can be obtained from the reduced density matrix of that subsystem. The definition is a generalization of the discrete formula given in Eq. (8), with the density matrix of the total system replaced by the reduced density matrix of the subsystem, which we denote by A:

$${\phi }_{A}=\arg ({Z}_{A})=\arg \left[{\rm{tr}}\left({W}_{A}\prod _{l=1}^{M-1}{\rho }_{A}({\kappa }_{l})\right)\right],$$
(9)

where ρA is the reduced density matrix obtained by the partial trace over the complement, \({\rho }_{A}(\kappa )={{\rm{tr}}}_{\bar{A}}\left({\rho }_{\kappa }\right)\), and \({W}_{A}={e}^{-i(2\pi /L)\sum _{j\in A}{x}_{j}}\). For the specific case of the half-filled SSHH model, we look at the reduced Zak phases of the spin-up and spin-down subsystem. By symmetry, the two phases are equal. We find that both change from 0 to π when Δt changes sign for both weak and strong interaction, i.e., for all interaction strength in the interval [0, 10].

Bulk-edge correspondence and triplon excitations

One of the most interesting aspect of topological band insulators is bulk–edge correspondence: The change of the Zak phase, which is a property of the eigenstates of the bulk, is accompanied by the appearance of localized edge states and signaled by the closing and reopening of the bulk band gap in the excitation energy spectrum, as demonstrated for the non-interacting SSH model in Fig. 1. To investigate bulk–edge correspondence in the strongly correlated case at half-filling, we compute the excitation energy spectrum of an open chain for U = 10 in Fig. 3. In our calculations, we classify states according to the spin projection, \({S}_{z}=(1/2)\sum _{j}({n}_{j,\uparrow }-{n}_{j,\downarrow })\), which commutes with the SSHH Hamiltonian. Owing to the symmetry between spin up and spin down, the states with Sz and −Sz have the same energy; the total spin S is also a good quantum number. In the trivial phase (Δt < 0), the ground state is non-degenerate and has Sz = 0 and S = 0. The energy gap closes at Δt = 0, at which point the lowest triplet excitation (S = 1, comprising the second lowest energy level with Sz = 0 and the two lowest energy levels with Sz = ±1) come down and become degenerate ground states in the non-trivial phase (Δt > 0). Thus the ground-state degeneracy χ changes to four across the topological phase transition.

Fig. 3: Properties of a half-filled SSHH model with OBC.
figure 3

a Eigenenergy spectrum in the strongly correlated regime (U = 10), colored according to the spin ∣Sz∣ of the eigenstates (states with ±Sz have the same energy). For each value of Sz up to 60 lowest energies are shown. The ground-state degeneracy χ changes from one to four at the topological phase transition due to the four degenerate edge states. b Eigenenergy spectrum in the non-trivial phase (Δt = 0.5) as a function of the on-site interaction. The ground-state degeneracy reduces from χ = 6 in the non-interacting regime to χ = 4 in the interacting regime due to the raising in energy of two ionic states where the edge is double occupied (see text). c The entanglement entropy between the first subsystem consisting of the two ends and the second one consisting of the remaining sites in the middle of the chain. d Spin correlation \(\left\langle {S}_{z,j}{S}_{z,k}\right\rangle\) in the non-trivial phase. The abscissa and ordinate are the sites j and k. The bulk dimers have perfect AFM correlation but the edges are uncorrelated with the rest; e same as d for Δt = 0 showing long-range AFM correlation; and f same as e for Δt = −0.5.

Full size image

When U ≫ t±, the half-filled SSHH model is well approximated by the Heisenberg model of localized spin-1/2 particles with alternating exchange interaction \({J_{\pm }}=4{t_{\pm }^{2}/U}\).44 For Δt far from zero, we find that, in the ground state, each dimer consisting of two sites coupled by the stronger exchange interaction, J>, is in a singlet state and is uncorrelated from the other dimers. The excitation gap between the adjacent bands of states in Fig. 3a is due to the excitation to the triplet states of these dimers. This triplon excitation gap is hence approximately J>, and the nth excited band has n triplons.

In the non-trivial phase (Δt > 0), the two edge spins at the two ends of the chain are decoupled from the dimers, and they can be in either a singlet state \(\left|{S}_{0}\right\rangle =\left(\left|\uparrow \downarrow \right\rangle -\left|\downarrow \uparrow \right\rangle \right)/\sqrt{2}\) or one of the three triplet states \(\left|{T}_{0}\right\rangle =\left(\left|\uparrow \downarrow \right\rangle +\left|\downarrow \uparrow \right\rangle \right)/\sqrt{2}\), \(\left|{T}_{+}\right\rangle =\left|\uparrow \uparrow \right\rangle\), \(\left|{T}_{-}\right\rangle =\left|\downarrow \downarrow \right\rangle\), all having zero energy. This explains the fourfold degeneracy of the ground state in the non-trivial phase. More specifically, at Δ = 0.5 and U = 10, our numerical calculation shows that the reduced density matrices of the two edges in the four ground states have around 99% overlap with the singlet and the triplet states, that is, \(\left\langle {\phi }_{j}\right|{\rho }_{j}^{{\rm{edges}}}\left|{\phi }_{j}\right\rangle \approx 0.99\) for j = 1, …, 4 where \(\left|{\phi }_{j}\right\rangle\) are \(\left|{S}_{0}\right\rangle ,\left|{T}_{0}\right\rangle\), and \(\left|{T_{\pm }}\right\rangle\).

The reduced density matrix of each strongly coupled dimer in the bulk at both Δt = ±0.5 satisfies \(\left\langle {S}_{0}\right|{\rho }^{{\rm{dimer}}}\left|{S}_{0}\right\rangle \approx 0.92\). Despite the significant nonvanishing value of the coupling J<, the dimers are almost totally uncorrelated with each other and the edges. This is due to monogamy of entanglement: Since each spin in a dimer is maximally entangled with its partner in the singlet state, it must be unentangled from any other spin.45

The decoupled dimer picture is further confirmed by the spin correlation \(\left\langle {S}_{z,j}{S}_{z,k}\right\rangle\) in the ground state with Sz = 0, as shown in Fig. 3d. In the non-trivial phase, the two spins in the dimer are perfectly anti-correlated with each other but uncorrelated with other spins, and the edge spins at the two ends are free. The spin correlation in the trivial phase is the same but without the two free ends (see Fig. 3f), since there is no weakly coupled edges in this case. A long-range AFM order develops near Δt = 0 as expected for a 1D Hubbard model46 (see Fig. 3e).

In the non-interacting regime, the edge states are identified by the localization of a single-particle state at the edges, while in the strongly correlated regime we have an effective spin model and the edge states are identified as uncorrelated spins. In order to identify the edge states at arbitrary interaction, we calculate the von Neumann entropy of the entanglement47 between the two ends of the chain and the rest in Fig. 3c. The formation of the edge states is indicated by the sharp drop in entanglement since these states are uncorrelated from the rest. Figure 3c shows a clear transition line where the edge states are formed, and the transition to the edge phase is more abrupt at large U. At U = 0, the transition point is not at Δt = 0 as expected for the SSH model owing to finite-size effects; we expect the transition point would become closer to Δt = 0 with increasing system size (see Supplementary Material).

Figure 3b shows the eigenenergy spectrum in the non-trivial phase as the on-site interaction changes from weak to strong. At U = 0, the ground state degeneracy is χ = 6 since in addition to the one singlet state and three triplet states the edges have two more degenerate ionic states where two electrons occupy the same edge, i.e. \(\left|\uparrow \downarrow ,\varnothing \right\rangle\) and \(\left|\varnothing ,\uparrow \downarrow \right\rangle\). For U > 0, the energy of these ionic states is lifted owing to the on-site repulsion, thus reducing the ground-state degeneracy to χ = 4.

Persistent long-range AFM order at quarter-filling

Unlike the half-filled case, the quarter-filled system does not exhibit conventional bulk–edge correspondence; this is clear from the lack of a gap in the eigenenergy spectrum (see Fig 4a). Recall that while there exists a gap in the charge excitation spectrum at quarter-filling discussed earlier, the absence of an eigenenergy gap is due to gapless spin excitations.44 The picture of decoupled dimers each in a singlet state no longer applies: deep in both phases (Δt far from 0) each dimer is occupied by roughly a single electron, and in the non-trivial phase the two ends are shared by one electron (see Fig. 4b for illustration). In contrast to the isolated edge spins at half-filling, we find that the edges are strongly correlated with the dimers in the bulk through the spin degree of freedom, and there is long-range AFM order in both trivial and non-trivial phases. To show this, we again plot the spin correlation \(\left\langle {S}_{z,j}{S}_{z,k}\right\rangle\), however, now with a different definition for the “sites” j and k: in the non-trivial phase we denote the first value j, k = 1 and the last value j, k = N∕2 + 1, for the left and right edges, and each value in between, j, k = 2, . . ., N∕2, is assigned to a dimer in the bulk. A similar definition of effective sites is used in the trivial phase but without the edges. The spin correlation between the edges and the bulk dimers in the non-trivial phase is shown in Fig. 4c; long-range AFM order is visible. The spin correlation in the trivial phase is exactly the same but with the two edges removed as evident in Fig. 4e. Long-range AFM order persists in both phases and also at the point of the phase transition (see Fig. 4d). This persistent long-range AFM order is the reason why the entanglement between the two edges and the bulk shows no clear drop when Δt changes sign (see Fig. 5b and the discussion below).

Fig. 4: Properties of a quarter-filled SSHH model with OBC.
figure 4

a Eigenenergy spectrum in the strongly correlated regime (U = 10), colored according to the spin ∣Sz∣ of the eigenstates. For each value of Sz the 40 lowest energies are shown. b An illustration of the particle occupation at each site and the long-range AFM correlation between the electron’s spins in the non-trivial phase. Each dimer in the middle of the chain is occupied by an electron, and the two edges are shared by one electron. c Spin correlation \(\left\langle {S}_{z,j}{S}_{z,k}\right\rangle\) in the non-trivial phase (Δt = 0.5). The abscissa and ordinate are the “effective sites” j and k defined as illustrated in b: The first and the last value refers to the two edges, while each value in between refers to each bulk dimer which is occupied by a single electron; d same as c for Δt = 0; and e same as d for Δt = −0.5. The long-range AFM correlation persists for all values of Δt.

Full size image
Fig. 5: Properties of a quarter-filled SSHH model with OBC under an external magnetic field.
figure 5

a Eigenenergy spectrum as a function of the field strength in the non-trivial phase. The levels are colored according to the spin projection Sz of the eigenstates. Only states with positive Sz are considered as the ones with negative Sz rise in energy in a magnetic field, and for each value of Sz the 10 lowest energies are shown. At the critical value \({E}_{B}^{(c)}\approx 0.1\) (indicated by the arrow), the ground state becomes maximally ferromagnetic (with all spins aligned along the field axis) and reduces to the non-interacting limit of the SSH model. b Entanglement entropy between the edges and the bulk for various field strength and Δt. Each contour separates regions where the ground state has different Sz (Sz = 3 indicates the maximally ferromagnetic ground state as there are 6 particles at quarter-filling). At zero field, the entanglement does not drop as Δt changes sign due to the long-range AFM correlation between the edges and the bulk at quarter-filling discussed above. At high field, the system reduces to the SSH model, and the entanglement drops for Δt > 0 due to the formation of the localized edge states, signaling the transition to the non-trivial phase. c Critical field strength at various values of the on-site interaction and hopping amplitude difference. d 1D slices of c for U = 0 (left panel) and Δt = 0.5 (right panel) showing the sharp decrease in the critical field strength with increasing on-site interaction.

Full size image

Discussion

Magnetic-field-induced transition to SSH ground state

The absence of an eigenenergy gap and uncorrelated edge states suggests that the SSHH model at quarter-filling is not a TI. However, the ground state can undergo a transition to the topological ground state of the non-interacting SSH model if a magnetic field is applied, resulting in the total Hamiltonian \({H}_{1}=H-({E}_{B}/2)\sum _{j}({n}_{j,\uparrow }-{n}_{j,\downarrow })\) where EB = gμBB and μB is the Bohr magneton and g the g-factor of electrons in the material. At a critical magnetic field, the ground state becomes the maximally ferromagnetic state with all of the electron spins aligned along the field axis (see Fig. 5a), and an energy gap is opened. Since there is no on-site interaction between particles with the same spin, this ground state must be the ground state of the non-interacting SSH model with N/2 electrons, and due to Pauli exclusion principle the highest energetic electron must occupy the mid-gap edge state shown in Fig. 1b. We find that the ground-state degeneracy at field strength larger than the critical value is χ = 2, agreeing with the fact that there are two degenerate edge states (left and right) in the SSH model. The transition to the SSH ground state is further confirmed in the entanglement entropy between the edges and the bulk in Fig. 5b. Without the field, the entanglement does not drop as Δt changes sign since the edges are correlated with the bulk through the persistent long-range AFM order at quarter-filling. At field strengths beyond the critical value, the entanglement drops sharply for Δt > 0 owing to the formation of the localized edge states in the non-trivial phase of the SSH model.

The magnetic-field-induced transition enables experimental realization of the SSH model in systems with local interactions. It is shown in Fig. 5c and more clearly in Fig. 5d, right panel, that the critical magnetic field reduces dramatically with increasing on-site interaction, meaning it is easier to realize the SSH model if the local interaction is stronger. For reaching the maximally ferromagnetic state, the last spin-down particle needs to be pumped to the next unoccupied single-particle energy level shown in Fig. 1b. When there is strong on-site interaction, this spin-down particle interacts strongly with the other particles with opposite spins, raising the energy, hence it costs less energy to pump this spin-down particle to a higher energy level, leading to a smaller required magnetic field.

Typical parameter values for dopant atoms in silicon and quantum dots in GaAs are \(\bar{t} \sim 1\) meV, \(U/\bar{t} \sim 10\), and assuming a g-factor of 2 our calculation gives a critical magnetic field of Bc ~ 2 T deep in the non-trivial phase at \(\Delta t/\bar{t}=0.5\), which is feasible. For electrons bounded to impurities or a quantum dot in a semiconductor host, the g-factor can deviate from the free electron value of 2,48 and the calculated critical field need to be rescaled. As long as the deviation is not too large, the critical magnetic field remains in the realistic range. One sees from Fig. 5d that, without the strong on-site interaction, at U = 0 for example, the required magnetic field is around ten times larger and therefore may not be realistically attainable for larger hopping amplitude. This emphasizes the importance of interaction.

We note that the jump in the critical field as Δt changes sign in the U = 0 limit (see Fig. 5d, left panel) is due to the formation of the mid-gap edge state in the non-trivial phase of the SSH model. In the trivial phase, the last spin-down particle needs to be pumped from the lowest energy to the highest energy in the lower bulk band of Fig. 1b, while in the non-trivial phase it needs to be pumped to the mid-gap edge state, which is higher in energy; hence, a larger magnetic field is necessary.

Experimental realization with nanofabricated semiconductor devices

In the last section of the paper, we propose a device architecture for realizing the transition to the topological phase of the SSH model described above. With state-of-the-art fabrication technology, it is possible to fabricate 1D chains of dopant atoms in silicon with scanning tunneling microscope (STM)49,50 or gate-defined quantum dots in GaAs,17 illustrated in Fig. 6a. Two leads, source and drain, are positioned close to one edge of the chain. Naturally there are potential barriers between these leads and the chain. Electrons from the source can tunnel through the barrier into the many-body state of the chain and then out to the drain. The A side gates are for tuning the on-site energy, and thus the chemical potential, by applying a voltage. The B side gates between the sites are for controlling the hopping amplitude (also commonly referred to as the “tunnel coupling” in these systems). A similar device without the source and the drain was fabricated for a chain of three quantum dots in ref. 17.

Fig. 6: Experimental proposal.
figure 6

a A device architecture for realizing and probing the topological phase transition of the SSH model. A chain of STM-fabricated dopant atoms, or gate defined quantum dots, is positioned in the center of a collection of electrodes and gates. The source and the drain leads are used for measuring the current of the electrons tunneling from the source, through the many-body state of the chain, to the drain. Since the source and the drain are close to the edges, the conductance is proportional to the charge excitation density at the edge of the many-body state (see text). The A side gates are used for tuning the on-site energy, and thus the chemical potential. The B side gates between the sites are used to control the hopping amplitudes. b The conductance spectrum at kT = 0.02 as a function of the chemical potential for Δt = 0.5 (solid red) and Δt = −0.5 (dotted blue) at zero field (top panel) and at a field strength just above the critical value (bottom panel). The spectrum for both values of Δt shows the lower Hubbard band (LHB) and upper Hubbard band (UHB) separated by the Mott gap, and each band is further separated by the charge excitation gap at quarter-filling. For Δt > 0, there is a conductance peak in the middle of this gap due to the tunneling through the edge states of the charge excitation. The formation of the strongly localized edge state in the SSH model above the critical field strength leads to a sharp increase in the conductance peak at quarter-filling in the bottom panel. c Density plot of the conductance around quarter-filling, showing the sharp increase of the peak at the critical field strength. All energies are scaled by the average value of the two hopping amplitudes.

Full size image

We now show how a transport measurement of the proposed device can probe the transition to the SSH ground state at quarter-filling and also the transition between the trivial and non-trivial SSH topological phases. When the tunneling rate Γ of the electron from the source/drain to the nearest site is much smaller than the hopping amplitudes between the sites and kT, which is typical for dopant atoms and quantum dots, we are in the sequential tunneling regime.18,51 As the chemical potential is varied, each time it matches an addition energy (shown in Fig. 2), the electron in the lead has enough energy to tunnel into the many-body state of the chain and out to the other lead, resulting in a peak in the tunneling current (see Fig. 6b). Thus the set of peaks in the conductance spectrum maps the addition energy spectrum. The conductance in the linear response regime, applicable when the bias between the source and the drain is much smaller than the hopping amplitudes and kT, is computed with the Beenakker’s formula.18,32,33,34

For the measurement using the source and drain on the left of Fig. 6a, the conductance peak at filling n is proportional to G0Dn where G0 = e2Γ∕(\(\hbar\)kT) and \({D}_{n}=| \left\langle {\Psi }_{0}^{(n)}| {c}_{1\uparrow }^{\dagger }+{c}_{1\downarrow }^{\dagger }| {\Psi }_{0}^{(n-1)}\right\rangle {| }^{2}\) where \({\Psi }_{0}^{(n)}\) is the many-body ground state of the chain at filling n. Dn can be interpreted as the charge excitation density at the left edge at filling n. Figure 6b, top panel, shows the conductance spectrum in the strongly correlated case (U = 10) for both signs of Δt, revealing the lower and upper Hubbard bands, separated by the Mott gap, of the addition energy spectrum in Fig. 2. There is also evidence of the charge gap at quarter-filling in the lower band and three quarter-filling in the upper band. For Δt > 0, there is a sharp conductance peak in the middle of the quarter-filling gap due to the edge state of the charge excitation at this filling, which has a high density of charge excitation at the edges. The lower panel shows the conductance spectrum at an applied field strength just above the critical value of the transition to the SSH model. The mid-gap conductance peaks are much higher since in the SSH model the edge states are much more localized. Observing a sharp increase in the conductance peak at quarter-filling, as shown by our calculation in Fig. 5c, can serve as experimental evidence of the transition to the SSH model of topological insulators. And the appearance of this peak as Δt changes sign from negative to positive can be a probe of the topological phase transition.

As an example, we consider parameter values that are typical of phosphorous donors in silicon: Γ = 0.001 meV, \(\bar{t}=4\) meV, and U = 40 meV. The conductance peak at quarter-filling is then of the order of 10−6 S at 1 K, leading to a tunneling current of 0.1 nA at 0.1 mV bias, which is large enough for detection.

In the above, we discuss a relatively long chain of 12 sites so that the lower and upper Hubbard bands in the conductance spectrum appear dense, but finite-size signature of the edge state can be observed with a much smaller number of sites in experiments, as low as N = 4. For N = 4, the “bulk” consists of a single dimer in the non-trivial phase and each half-band either side of the quarter-filling gap in Fig. 6b has a single peak. It is better to have N = 6 so that each half-band has two close peaks and hence can be easily distinguished from the edge-state peak in the middle of the gap. These system sizes are feasible with current technology.

Robustness against disorders

Disorders in the on-site energy and hopping amplitude are unavoidable in real experiments. We investigate the effect of disorders by adding to the SSHH Hamiltonian the on-site energy term \({H}_{{\rm{on-site}}}=\sum _{j,\sigma =\uparrow \downarrow }{\epsilon }_{j}{n}_{j,\sigma }\), where each ϵj is chosen uniformly at random from a range [−δE, δE], and we add to each hopping amplitude, tj = 1 + (−1)jΔt, a variation chosen uniformly at random from a range [−δt, δt]. To study the robustness of the magnetic-field-induced topological phases above the critical field strength, we look at the distribution of the addition energy gap at quarter filling, the distribution of the critical magnetic field required, the signature of the edge state in the conductance spectrum, and the distribution of the many-body Zak phase.

We find that, in the presence of both on-site energy and hopping amplitude disorder, the topological phases are robust as long as δt + δE < ∣Δt∣. This is expected since the gap at quarter-filling in the fully spin-polarized regime is Δt, and the above condition makes sure that the gap is not closed by the disorder. Also, δt < ∣Δt∣ means that the weak–strong order between the odd and even hopping amplitudes is preserved, that is, the couplings within the bulk dimers are always stronger than those between them and with the edges. We show in Fig. 7 the numerical evidence for robustness when Δt = 0.5 and δt = δE = 0.5Δt. Note that only in this figure we choose to use absolute unit in meV for a more direct connection with experiments.

Fig. 7: Robustness of the magnetic-field-induced topological phases against disorders for a chain with N = 6 sites.
figure 7

The mean values of the weak and strong hopping amplitudes are 2 and 6 meV, and hence ∣Δt∣, which is half of the hopping amplitude difference, is 2 meV. The on-site interaction is U = 40 meV. The maximum variations in the on-site energy and hopping amplitude disorder are both ∣Δt∣∕2 (see text). The distributions are generated from a sample of 5000 random instances. a Probability distributions of the addition energy gap at quarter filling for B > Bc (dark green) and a finite-size addition energy separation within the lower Hubbard band (dark blue), showing that the gap is much larger and distinct from other energy differences arising from finite-size effects. b Probability distribution of the critical magnetic field Bc assuming a g-factor of 2. c A typical conductance spectrum of the disordered system for B > Bc in the trivial phase (dashed blue) and non-trivial phase (solid red). The appearance of the high rising edge-state peak in the non-trivial phase is clearly visible. d Probability distribution of the many-body Zak phase for B > Bc, which shows only a small deviation from the ideal value of π.

Full size image

Figure 7a shows the distribution of the addition energy gap at quarter filling, which is Ead(N∕2 + 1) − Ead(N∕2) in the trivial phase, for B > Bc, compared with the distribution of an addition energy separation due to finite-size effect at a lower filling within the lower Hubbard band, Ead(N∕2) − Ead(N∕2 − 1). One sees that it is highly probable that the gap is much larger than the finite-size energy separations and thus can be identified in experiments. The critical magnetic field required for the transition to the non-interacting SSH limit varies within a realistic range, as demonstrated in Fig. 7b. The critical magnetic field can be reduced if one chooses a smaller hopping amplitude, but this will leave less room for disorder, hence there is a trade-off. A typical conductance spectrum (for B > Bc) of the disordered system in Fig. 7c shows clearly the high-rising edge-state peak in the non-trivial phase in the middle of the split Hubbard bands of the trivial phase, which can be used for identifying the topological phase transition. We see similar clear signatures in all the 20 random instances we generated for the conductance spectrum. The many-body Zak phase in the non-trivial phase for B > Bc also deviates very little from the ideal value of π.

Similar robustness behavior is observed for the cases of pure on-site energy disorder when δE < ∣Δt∣ and pure hopping amplitude disorder when δt < ∣Δt∣. The signatures of the topological phases are of course erased for very strong disorders, for example, when δt = δE = 2∣Δt∣ (see Supplementary Material). In a donor chain in silicon, the hopping amplitude oscillates rapidly with the donor separation due to the inter-valley interference in the wavefunction. Even with a positional variation within a silicon unit cell, the hopping amplitude can drop to close to zero according to effective mass theory,18 thus it might be challenging to limit the hopping amplitude disorder to within the range [−∣Δt∣, ∣Δt∣]. This can be mitigated by fabricating inter-donor side gates depicted as B in Fig. 6 for controlling the hopping amplitude. If acceptors are used instead, there is no hopping amplitude oscillation owing to the absence of intervalley interference,52 thus the B side gates are not needed, but the A side gates are still required for varying the chemical potential.

We also investigate the robustness against disorder of the topological phases at half-filling in zero magnetic field, as shown in Fig. 3. Recall that in the large-U limit the properties of the SSHH model can be understood from an effective Heisenberg model of local spins interacting with staggered exchanges, \({J}_{\pm }=4{t}_{\pm }^{2}/U\). The characteristic spin correlation in Fig. 3d, where each dimer in the bulk is strongly correlated while the correlations between the dimers and with the edges are negligible, is preserved if all the odd exchanges are smaller than all the even exchanges. This means the disorder in the hopping amplitude, δt, should be smaller than ∣Δt∣, similar to the quarter-filling case. One interesting difference from the quarter-filling case is that the on-site energy disorder can now be much larger. In the picture of the staggered Heisenberg model, one electron is localized at each site, thus changing the on-site energy is akin to changing an energy constant in the Heisenberg Hamiltonian, which does not affect the spin excitation spectrum showing the triplon bands in Fig. 3a. This is true as long as the on-site energy disorder is smaller than the on-site interaction U. For stronger variations, where the energy of one site is lower than that of another by more than U, double occupancy will be favored, leading to the breakdown of the Heisenberg picture. When both hopping amplitude and on-site energy disorders are present, we find that the topological phases at half filling are robust when δt < ∣Δt∣ and δt + δE < U. More specifically, the gap above the ground state in Fig. 3a and the characteristic spin correlation of Fig. 3d remain intact, and the reduced Zak phase deviates very little from its ideal value. We refer the reader to Supplementary Material for the numerical results.

Finally, we comment briefly on other possible imperfections. The SSHH model assumes the electrons are phase coherent throughout the length of the chain. A finite chain of 6 sites of dopant atoms in silicon can be made <50 nm, while the phase coherence length in STM-fabricated samples at low temperature can be well >100 nm, as inferred from weak-localization experiments,53 and in GaAs-based samples the phase coherence length can be as large as a μm.54 Spin–orbit coupling is not taken into account in our model, but it can be neglected if its energy scale is much smaller than the hopping amplitude. For Si:P, even with the enhancement due to external fields the energy scale of the spin–orbit coupling is of the order of 10−6 μeV, as inferred from the spin-flip rate in the region of ms,55 which is negligible compared with a hopping amplitude in the meV or sub-meV range. In GaAs, spin–orbit coupling may lead to unwanted effect such as spin flip56 or spin-flip tunneling between the dots57; however, the energy scale of these effects are just as small. In our calculation of the conductance spectrum, we assume an energy-independent tunneling rate, Γ, between the system and the leads, but in reality electrons tunneling to the upper Hubbard bands have energies much closer to the top of the barrier between the leads and the system, resulting in a larger tunneling rate. The conductance peaks in the upper Hubbard band therefore should be much higher than those in the lower band.

In summary, we have investigated the topological phases of a 1D Fermi–Hubbard model in the strongly correlated regime. We introduce the concept of the reduced Zak phase, defined based on the reduced density matrix of a subsystem, and show that the topological phases at half-filling can be characterized by this phase. This reduced phase might be useful for studying the topological phases of a subsystem in a larger interacting system or an open system interacting with the external environment. From a study of entanglement and spin correlation, we demonstrate the bulk–edge correspondence in the half-filled system. At quarter-filling, the model does not exhibit properties of a topological insulator, but it can be transformed to the topological ground state of the non-interacting SSH model by applying a magnetic field. Finally, we propose a promising experimental realization with dopant atoms in silicon or quantum dots in GaAs. The scheme is robust against significant disorder in the hopping amplitude and on-site energy.

Methods

We use the Lanczos algorithm to diagonalize the Hamiltonian in the occupation basis58

$${\left|1,0,1,\ldots ,1\right\rangle }_{\uparrow }\otimes {\left|0,1,0,\ldots ,1\right\rangle }_{\downarrow },{\rm{etc.}},$$
(10)

where the first part is for the spin-up particles and the second for the spin-down ones. A 1 at position j means the jth site is occupied and a 0 means the reverse. These basis states are then numbered according to the decimal value of its binary string. From this representation, the reduced density matrix of the spin-up or spin-down subsystem can be computed in a straightforward manner. For computing the entanglement entropy between the edges and the rest of the chain, we define the local Hilbert space of each site as \(\left|\varnothing \right\rangle ,\left|\uparrow \right\rangle ,\left|\downarrow \right\rangle ,\left|\uparrow \downarrow \right\rangle\), that is, a qudit with dimension 4, the entanglement entropy can then be calculated for the resulting system of qudits.47