Topological domain wall states in a non-symmorphic chiral chain
Wojciech Brzezicki Timo Hyart
cond-mat.str-el






International Research Centre MagTop, Institute of Physics, Polish
Academy of Sciences, 

 Aleja Lotnikow 32/46, PL-02668 Warsaw, Poland



International Research Centre MagTop, Institute of Physics, Polish
Academy of Sciences, 

 Aleja Lotnikow 32/46, PL-02668 Warsaw, Poland

The Su-Schrieffer-Heeger (SSH) model, containing dimerized hopping
and a constant onsite energy, has become a paradigmatic model for one-dimensional topological phases, soliton excitations and fractionalized charge  in the presence of chiral symmetry. Motivated by the recent developments in engineering artificial lattices, we study an alternative model
where hopping is constant but the onsite energy is dimerized. We find that it has a non-symmorphic
chiral symmetry and supports topologically distinct phases described by a ℤ_2 invariant ν. In the case of multimode ribbon we also find topological phases protected by hidden
symmetries and  we uncover the corresponding ℤ_2 invariants  ν_n. We show that, in contrast to the SSH case, zero-energy states do not necessarily appear at the boundary between topologically distinct phases, but instead these systems support a new kind of bulk-boundary correspondence: The energy of the topological domain wall
states typically scales to zero as 1/w, where  w is the width of the domain wall separating phases with different topology. Moreover, under specific circumstances we also find a faster scaling e^-w/ξ, where ξ is an intrinsic length scale.
We show that the spectral flow of these states and  the charge of the domain walls are different than in the case of the SSH model.



Topological domain wall states in a non-symmorphic chiral chain
    Timo Hyart
    Received: date / Accepted: date
===============================================================




The Su-Schrieffer-Heeger (SSH) model was originally introduced to describe the properties of conducting polymers, where the spontaneous symmetry breaking leads to dimerization of the sites along the chain <cit.>.  Due to two-fold degeneracy of the ground state  a new type of excitation, a domain wall (DW) between different bonding structures, can exist. For the conducting polymers the width of the DW excitations is large and  they can propagate along the chain. Thus, they can be considered as solitons in analogy to the shape-preserving propagating solutions of the nonlinear differential equations <cit.>. Moreover, the solitons in the SSH model have a remarkable effect on the electronic spectrum leading to an appearance of a bound state in the middle of the energy gap. This midgap state is understood as a topologically protected boundary mode and the SSH model serves as a paradigmatic example of chiral symmetric topological insulator <cit.>. Namely, the chiral symmetry allows to block-off-diagonalize the Hamiltonian and the winding of the determinant z_k of the off-diagonal block around the origin as a function k determines a topological invariant (see Fig. <ref>). Because this invariant is different on two sides of the DW, each DW carries zero-energy bound state. The DWs come in pairs so that the spectral flow is symmetric around zero-energy  and each DW carries a charge q=±1/2 in analogy to the fractionally charged excitations studied in the quantum field theory  <cit.>.
This can be also understood in terms of modern notions of bulk obstructions and filling anomalies <cit.>.




The idea of soliton excitations reappears in the context
of 1D diatomic polymers in a form of the Rice-Mele model <cit.>, which has been studied also in contexts of ferroelectricity
<cit.> and  organic salts <cit.>. In this model
 not only bond length alternates but also the onsite energy (mass) takes
opposite sign for the even/odd lattice sites (Fig. <ref>).
A very interesting
feature of such a model is that its solitons can carry irrational
charge q=1/2(1± f) <cit.>, where f describes the breaking of the  chiral symmetry <cit.>. The interest for these models has revived because they can engineered in photonic systems <cit.>, optical lattices <cit.> and nanostructures
<cit.>  in a controlled way, and in these systems also their emergence from spontaneous symmetry breaking  <cit.> and the properties of the solitons can be tuned using external parameters <cit.>. Motivated by the new possibilities opened by these recent
developments we focus  on a special case of the Rice-Mele
model where all the hopping amplitudes are equal and only the mass
term alternates. We show that in this case the model has an interesting non-symmorphic (NS) chiral symmetry and it supports a topologically nontrivial phase described by a non-symmorphic chiral ℤ_2 invariant ν. This invariant was found by Shiozaki, Sato and Gomi in their pioneering work on non-symmorphic topological insulators <cit.> and therefore we name the special case of the Rice-Mele model as the Shiozaki-Sato-Gomi (SSG) model.  The peculiar property of the SSG model is that the bulk topological
invariant does not guarantee the existence of the end states in an open system, because the boundary always breaks the NS chiral symmetry <cit.>.



In this paper we analytically derive
exact phase diagrams of SSG nanoribbons of arbitrary width  and uncover hidden symmetries
relying on interchange of transverse and longitudinal modes. In addition to the NS chiral ℤ_2 invariant ν the multimode ribbons 
 support ℤ_2 invariants  ν_n protected by the hidden
symmetries.
These invariants lead to a new kind of bulk-boundary correspondence: The energy of the topological domain wall
states typically scales to zero as 1/w, where  w is the width of the domain wall separating phases with different topology. Moreover, under specific circumstances we also find a faster scaling e^-w/ξ, where ξ is an intrinsic length scale (Figs. <ref>, <ref> and <ref>). 
The NS chiral symmetry in SSG model leads to several important differences in comparison to the SSH model (Fig. <ref>): (i) In SSH model the topological zero-energy end or DW states come in pairs and have zero energy for any DW width, whereas the SSG model supports unpaired DW states approaching zero energy with increasing w.  Note that here we consider finite chains with open boundary conditions and the chains are made out of integer number of unit cells. An inifinite SSH chain could also have a unpaired DW state but in finite system the end states always guarantee the the zero-energy states come in pairs. The SSG model is different because even in this kind of situation we can obtain a single zero-energy state at the smooth domain wall but there is no low-energy states at the sharp interface with vacuum. (ii) In SSH model the charge of the DWs is q=±1/2, whereas for the SSG model we get  irrational charges q=±1-δ/2, ±1+δ/2 for solitons and antisolitons depending on whether the zero-energy  state is occupied or empty. The DWs in SSG model separate
regions with different onsite energies ± m_1 and ± m_2 (mass terms) in the two sublattices (see Figs. <ref> and <ref>), and δ=1/2(ζ_2-ζ_1), where ζ_i
is the difference of the bulk filling factors of the two sublattices in the region with mass m_i.

The k-space SSG Hamiltonian for a multimode wire is 

    H_k    =     -mσ_zτ_z+2t_xcosk2(cosk2σ_x-sink2σ_y)
        +     t_yτ_x+2t_dcosk2(sink2σ_x+cosk2σ_y)τ_y,

where m is the mass, t_x
and t_y are hopping amplitudes in x
and y directions and t_d is the diagonal hopping amplitude with opposite signs in the two sublattices [Fig. <ref>(a)]. Here  σ_α are
Pauli matrices describing  the unit cell in the x direction, τ_α are L_y× L_y matrices describing
transverse hopping (τ_x)_pq=δ_1,|p-q|, (τ_y)_pq=iδ_1,|p-q|(-1)^p,
(τ_z)_pq=δ_pq(-1)^p and L_y is the width in the y direction.  For L_y=2
operators τ_α are equivalent to σ_α and
in the special case L_y=1 we set τ_z=1
and τ_x=τ_y=0. Because the SSG model belongs to symmetry class AI <cit.> (for list of symmetries see
Appendix <ref>) the only known topological invariant is the NS chiral ℤ_2 invariant ν <cit.>.





 
To calculate the topological invariant ν we rewrite the Hamiltonian in a block off-diagonal form in
the eigenbasis of NS chiral operator S_k=sink2σ_xτ_z+cosk2σ_yτ_z
that anticommutes with H_k. The determinant z_k of
the off-diagonal block is a complex number and its trajectory in the
complex plane as a function k determines the topological invariant <cit.>. Namely, due to non-symmorphicity
of S_k the period of z_k is 4π and in a properly chosen basis it satisfies  constraints ℑ𝔪z_k=-ℑ𝔪z_k+2π
and ℜ𝔢z_k=ℜ𝔢z_k+2π, so that the trajectory z_k starts at k=0 and ends at k=2π with same real part but with opposite imaginary part. Thus the parity of the number of times
the trajectory z_k crosses the positive real semiaxis for k∈[0,2π]
is a ℤ_2 topological invariant because it cannot be
changed without closing the gap or breaking the NS chiral symmetry (see Fig. <ref>). 
In our case the mirror symmetry M_x=cosk2-isink2σ_z
becomes identity in the eigenbasis of S_k (see Appendix  <ref>) so that 
z_k=z_-k. For this reason ℑ𝔪z_π=0 and the
formula for the ν gets simplified to 

    ν= sign ℜ𝔢z_π

in analogy to the simplification of the invariant for topological insulators in the presence of inversion symmetry <cit.>. The band-inversion
corresponding to a change of ν happens at k=π and m=0. 
We find that

    ν=12[1+ sign(m)]     if  L_y=2n-1 
    
    			0,     if   L_y=2n   (n∈ℕ_+)




In Fig. <ref> we show the topological phase diagrams of the
SSG model at L_y=6,7 as functions of m/t_x and t_d/t_x,
setting t_y=t_x. Surprisingly we find more phases than predicted
by the ν invariant. The gap closes not only for m=0  at k=π
when L_y is odd but also for any L_y along lines m=m_n at k=k_n, where 

    m_n=t_dt_y/t_xε_n^2, k_n=±2arccos[t_y/2t_xε_n],  ε_n=2cosnπL_y+1

 and n=1,2,…,⌊ L_y/2⌋
provided that 

    |t_y/2t_xε_n|≤1.

This means that in the limit of very wide ribbon (L_y→∞)
the phase diagram consists of a quasi-continuous set of lines t_d=γ_nm
with slopes γ_n ranging between t_x4t_y
and ∞. The natural question to ask now is what is the origin
of these gap-closing lines? The answer are the hidden symmetries,
that can be found at the magical k_n points, yielding to new 
ℤ_2 invariants ν_n. 




To see the hidden symmetries we rotate σ_α
matrices by angle k2 around the z-axis and use the eigenbasis of τ_x to transform the operators τ_α in a block-diagonal form, where the blocks are given by (see Appendix
<ref>) 

    τ_x,n=ε_nσ'_x,  τ_y,n=ε_nσ'_y,  τ_z,n=σ'_z.

and σ'_α is a new set of Pauli matrices. For odd L_y the
blocks n=0 is given by τ_x,0=τ_y,0=0
and τ_z,0=1. After this transformation the Hamiltonian (<ref>)
also has a block-diagonal form 

    H'_k,n    =     -mσ_zσ'_z+2t_xcosk2σ_x
        +     2t_dε_ncosk2σ_yσ'_y+t_yε_nσ'_x.

Now we notice that H'_k,n is invariant under interchange
of σ and σ' operators if 2t_xcosk2=t_yε_n
which provides the condition for gap closing points in the k-space [Eq. (<ref>)].
The spin-interchange X_12σ⃗X_12=σ⃗' and vice-versa is realized by  operator  X_12=12(1+σ⃗·σ⃗') 
 <cit.>.
The spectrum of X_12
consists of single -1 (singlet state) and three +1 (triplet states) eigenvalues. Thus H'_k,n
in the eigenbasis of X_12 becomes block-diagonal with one block
being 1×1 and the other being 3×3. Therefore we can define a topological
ℤ_2 invariant based on the sign of the matrix element of the 1×1 block. It changes at the gap closing
lines defined by Eq. (<ref>) so that it takes the form 

    ν_n=12{1+ sign [m-m_n]}.

We conclude that the full topological description  the SSG model is given by a vector {ν,ν_1,ν_2,…,ν_⌊ L_y/2⌋}
because changes of these invariants coincide with all the gap closing
lines in the phase diagram. Each hidden symmetry n=1,2,…,⌊ L_y/2⌋ exists only if condition (<ref>) is satisfied. Therefore, it is possible to connect phases with different ν_n without closing the bulk gap by choosing a path in the parameter space which goes outside the region where the hidden symmetry exists [see Fig. <ref>(c)], distinguishing the hidden-symmetric topological phases from the ones protected by structural symmetries.
Note that in the two-dimensional limit where t_x=t_y and the system is periodic in both directions the hidden symmetries become a mirror symmetry with respect to the x̂+ŷ line, but the hidden symmetries can exist even if the mirror symmetry is absent (see Appendix <ref>).


After establishing topological properties of the SSG ribbon we now turn our attention to the bulk-boundary correspondence. In Ref. <cit.> it was argued that a NS chiral ℤ_2 invariant does not generically support
end states in one dimension because the boundary necessarily breaks
the S_k symmetry. This however does not exclude a special type
of smooth DWs from having zero energy bound states. To obtain analytical insights we can develop a continuum model for odd L_y by
expanding H_k around gap closing point at k=π and
get 

    H_ eff=vδ kσ_x-mσ_y (v>0).

Now we create a DW of width w in the real space 

    m(x)=m_2+m_12+m_2-m_12tanhxw

between regions with positive mass m_1=m_0 and
negative mass m_2=-m_0  (m_0>0), separating phases with ν=1 and ν=0 [Figs. <ref>(b) and <ref>(a)], and we find that a zero-energy eigenstate of H_ eff exists in a form

    ψ(x)=([                0 (coshxw)^-m_0w/v ])^T/ N,

where N is a normalization factor.
This however does not take into account the fact that the chiral symmetry
of H_ eff becomes NS if one goes beyond linear order
in δ k. Therefore, we implemented numerically such a DW, 
for which |m| is constant within a unit cell as follows from Eq. (<ref>),
and calculated the energy of the DW state as a 
function of w [Figs. <ref>(c)(d) and <ref>(a)]. This way we find that the energy of the topological DW state approaches zero as 1/w whereas the energies of the other bound states  scale as 1/√(w).
This  means that the topological DW state can be distinguished from other bound states based on the scaling behavior because its energy approaches zero faster than the energies of the other states. 
We emphasize that the scaling of the energy to zero as 1/w or faster is  a robust property of these topological DW states.  In the continuum model this state would have zero-energy  and the lattice effects can give maximally a correction proportional to 1/w.  The robusness against perturbations preserving the NS chiral symmetry is demonstrated in Appendix <ref>.
On the other hand, in the Appendix <ref> we show that the DW state is exponentially localized around the center of the DW, i.e. such position x=x_0 that m(x_0)=0, and the localization length decreases with increasing w, as one could expect.




The 1/√(w) scaling of the bulk gap can be
understood as well. By inserting an expansion m(x)≃ m_0 xw around x=0   to Hamiltonian  H_ eff and eliminating ψ_1 we obtain a Harmonic oscillator equation for ψ_2

    -v^2 ψ_2”(x)+m_0^2/w^2x^2ψ_2(x)=(E^2+m_0 vw)ψ_2(x).

The energies of this problem are given by

    E_n=√(2n m_0 v/w),  n=0,1,2,...

The n=0 solution gives the topological DW state and the energies of the other states scale as 1/√(w). 
We have just shown that the in-gap state at the DW between two topologically distinct
domains follows from the continuum-limit model and adiabatic evolution of the states when DW is sufficiently smooth.
This is quite a different case to the one discussed in <cit.> where a surface state also appears in a presence of a NS 
symmetry but it also requires a surface potential. 







We find even more striking bulk-boundary correspondence for phases
described by ν_n invariants. In Figs. <ref>(d) and <ref>(c),(d)
we show the scaling of energies of the topological DW states and the non-topological states as a function of w when the masses 
 m_1 and m_2 are chosen so that 
the DW separates two different ν_n phases. The energies of the non-topological states behave in the same way as before, scaling
as 1/√(w), but the energies of the hidden-symmetry protected topological DW states scale as e^-w/ξ. By expanding Hamiltonian around the gap closing points ± k_n  we get two
similar bound state solutions as for k=π. However, these solutions do not give zero-energy states even in the continuum model  because the gap closes at two different momenta so that these bound states hybridize leading to non-zero energy. Nevertheless, in Appendix <ref>
we show using the properties of the Schwartz functions that in the continuum model their
overlap vanishes exponentially fast with w which allows the possibility of the exponential scaling. Nevertheless, the lattice effects could lead to 1/w corrections also in this case. Therefore, we have studied more carefully how the scaling of the energies depends on the details of the DW. We find that  for a specific type of DW   it is possible to obtain an exponential decay of energy e^-w/ξ also in the case of the DWs separating different ν phases,  and conversely it is possible to obtain 1/w scaling for a DW separating different ν_n phases (see Fig. <ref>). 
In this way, by manipulating the details of the DWs, we are able to see that the behavior of both types of DW states are equivalent.
This leads to a robust conclusion that the energies of the DW states scale as 1/w or faster. While the 1/w scaling is expected to be generic based on the analytical arguments given above, the faster exponential scaling, which we numerically find in some specific circumstances, means that the 1/w corrections are not always present. The analytical model-independent theoretical understanding of the conditions for the exponential scaling is an interesting direction for future research. By setting m_1 and
m_2 is such a way that all the gap closing lines are crossed
on the way from m_1 to m_2 we can always obtain extensive
number of DW states L_y  both for even and odd L_y.


An interesting proterty of the SSG model is that when the width of the DW w increases a single state separates
from the bulk spectrum and tends to zero from above or below [Fig. <ref>(c)]. This asymmetric spectral flow needs to be taken into account when calculating the charges for solitons and antisolitons  (see Appendix <ref>). For L_y=1 we obtain q_0=±1-δ/2, ±1+δ/2 for solitons and antisolitons depending on whether the zero-energy  state is occupied or empty. Here δ=ζ_2-ζ_1/2,

    ζ_i=2/πm_i/√(m_i^2+4t_x^2)K(4t_x^2/m_i^2+4t_x^2)

are the differences of the bulk filling factors of the two sublattices in the region with mass m_1>0 and m_2<0
and K(x) is the complete elliptic integral of the first kind. The DWs between different hidden-symmetric topological phases carry charges q_n=0, ± 1, so that for a general DW the charge is the sum of   q_0  and the charges q_n contributed by the transverse modes supporting transitions between different hidden-symmetric topological phases.

To summarize, we have analytically described the topological properties of the SSG model and propose it as a paradigmatic model for NS chiral-symmetric topological phases. We have shown that a smooth
DW supports zero energy state(s) if the DW separates regions
with different NS chiral invariant ν or different hidden-symmetry
invariants ν_n.  In addition to engineered artificial lattices <cit.>  our findings are also relevant  in the context of low-dimensional binary compounds
supporting surface atomic steps. In these systems the surface steps lead to one-dimensional topological modes and the system obeys NS chiral symmetry, so that DWs between topologically distinct phases can support  DW states <cit.>, providing a possible explanation for the zero-bias conductance peak observed in the recent experiment <cit.>.


The work is supported by the Foundation for Polish Science through
the IRA Programme co-financed by EU within SG OP Programme.
W.B. also acknowledges support by Narodowe Centrum Nauki 
(NCN, National Science Centre, Poland) Project No. 2016/23/B/ST3/00839.


10

SSH79 W. P. Su, J. R. Schrieffer, and A. J. Heeger,
“Solitons in polyacetylene”, Phys. Rev. Lett. 42, 1698-1701
(1979). 

SSH80 W. P. Su, J. R. Schrieffer, and A. J. Heeger,
“Soliton excitations in polyacetylene”, Phys. Rev. B 22,
2099-2111 (1980). 

SSH88 A. J. Heeger, S. Kivelson, J. R. Schrieffer,
and W. P. Su, “Solitons in conducting polymers”, Rev. Mod. Phys.
60, 781-850 (1988). 


Kane2010
M. Z. Hasan and C. L. Kane, "Colloquium: Topological insulators"
Rev. Mod. Phys. 82, 3045 (2010).






Chiu2016
C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, "Classification of topological quantum matter with symmetries",
Rev. Mod. Phys. 88, 035005 (2016).

Jakiw76 R. Jackiw and C. Rebbi, "Solitons with fermion
number 1/2”, Phys. Rev. D 13, 3398 (1976). 

Kha19 E. Khalaf, W. A. Benalcazar, T. L. Hughes, and R. Queiroz,  "Boundary-obstructed topological phases", 
arxiv: 1908.00011; D. V. Else, H. C. Po, and H. Watanabe, "Fragile topological phases in interacting systems", arxiv: 1809.02128.

Rice82 M. J. Rice and E. J. Mele, "Elementary Excitations
of a Linearly Conjugated Diatomic Polymer”, Phys. Rev. Lett. 49,
1455 (1982).

Brazowskii80 S. A. Brazovskii, "Self-localized excitations
in the Peierls-Fröhlich state”, Zh. Eksp. Teor. Fiz. 78,
677-699 (1980).

Jackiw83 R. Jackiw and G. Semenoff, "Continuum Quantum
Field Theory for a Linearly Conjugated Diatomic Polymer with Fermion
Fractionization”, Phys. Rev. Lett. 50, 439 (1983).

Vaderbilt93 D. Vanderbilt and R. D. King-Smith, "Electric
polarization as a bulk quantity and its relation to surface charge”,
Phys. Rev. B 48, 4442 (1993). 

Onoda04 S. Onoda, S.Murakami, and N. Nagaosa, "Topological
Nature of Polarization and Charge Pumping in Ferroelectrics”, Phys.
Rev. Lett. 93, 167602 (2004). 

Soos07Z. G. Soos and A. Painelli, "Metastable domains
and potential energy surfaces in organic charge-transfer salts with
neutral-ionic phase transitions”, Phys. Rev. B 75, 155119
(2007). 

Tsuchiizu16 M. Tsuchiizu, H. Yoshioka, and H. Seo,
"Phase Competition, Solitons, and Domain Walls in Neutral–Ionic
Transition Systems”, J. Phys. Soc. Jpn. 85, 104705 (2016).



Wilczek81 J. Goldstone and F. Wilczek, “Fractional
Quantum Numbers on Solitons”, Phys. Rev. Lett. 47, 986
(1981). 


Kitagawa12
T. Kitagawa, M. A. Broome, A. Fedrizzi, M. S. Rudner, E. Berg, I. Kassal, A. Aspuru-Guzik, E. Demler, and A. G. White, "Observation of topologically protected bound states in photonic quantum walks", Nat. Commun. 3, 882 (2012). 

Schomerus
J. Arkinstall, M. H. Teimourpour, L. Feng, R. El-Ganainy, H. Schomerus, "Topological tight-binding models from non-trivial square roots", Phys. Rev. B 95, 165109 (2017).

Atala13
M. Atala, M. Aidelsburger, Julio T. Barreiro, D. Abanin, T. Kitagawa, E. Demler, and I. Bloch, "Direct measurement of the Zak phase in topological Bloch bands", Nature Physics 9, 795 (2013). 

Wang13 
L. Wang, M. Troyer, and X. Dai, “Topological
Charge Pumping in a One-Dimensional Optical Lattice”, Phys. Rev.
Lett. 111, 026802 (2013). 

Lohse16
M. Lohse, C. Schweizer, O. Zilberberg, M. Aidelsburger, I. Bloch,  "A Thouless quantum pump with ultracold bosonic atoms in an optical superlattice", Nature Physics 12, 350 (2016).

Nakajima16
S. Nakajima, T. Tomita, S. Taie, T. Ichinose, H. Ozawa, L. Wang, M. Troyer, Y. Takahashi,
"Topological Thouless pumping of ultracold fermions",  Nature Physics 12, 296 (2016).

Leder16
M. Leder, C. Grossert, L. Sitta, M. Genske, A. Rosch, M. Weitz,
"Real-space imaging of a topologically protected edge state with ultracold atoms in an amplitude-chirped optical lattice",
Nature Communications 7, 13112 (2016). 


Meier16
E. J. Meier, F. A. An and B. Gadway, "Observation of the topological soliton state in the Su-Schrieffer-Heeger model", Nature Communications 7, 13986 (2016). 



Yan19 L. Yan and P. Liljeroth, "Engineered electronic
states in atomically precise nanostructures: artificial lattices and
graphene nanoribbons”, arxiv: 1905.03328. 

Park16 J.-H. Park, G. Yang, J. Klinovaja, P. Stano, and D. Loss
Phys. Rev. B 94, 075416  (2016).

Drost17
R. Drost, T. Ojanen, A. Harju, P. Liljeroth, "Topological states in engineered atomic lattices", Nature Physics 13, 668 (2017).

Huda18
M. Nurul Huda, S. Kezilebieke, T. Ojanen, R. Drost, P. Liljeroth, "Tuneable topological domain wall states in engineered atomic chains", arXiv:1806.08614.

Fischer18
D. J. Rizzo, G. Veber, T. Cao, C. Bronner, T. Chen, F. Zhao, H. Rodriguez, S. G. Louie, M. F. Crommie, F. R. Fischer,  "Topological band engineering of graphene nanoribbons", Nature 560, 204 (2018). 

Fasel18
O. Gröning, S. Wang, X. Yao, C. A. Pignedoli, G. Borin Barin, C. Daniels, A. Cupo, V. Meunier, X. Feng, A. Narita, K. Müllen, P. Ruffieux, R. Fasel,  "Engineering of robust topological quantum phases in graphene nanoribbons", Nature 560, 209 (2018). 










Weber 19 M. Weber, F. P. Toldin, and M. Hohenadler.
“Competing Orders and Unconventional Criticality in the Su-Schrieffer-Heeger
Model”, arxiv: 1905.05218. 

Pupillo08
G. Pupillo, A. Griessner, A. Micheli, M. Ortner, D.-W. Wang, and P. Zoller, "Cold Atoms and Molecules in Self-Assembled Dipolar Lattices",
Phys. Rev. Lett. 100, 050402  (2008).


Herrera2011
 F. Herrera, R. V. Krems, "Tunable Holstein model with cold polar molecules", Phys. Rev. A 84, 051401 (2011).

Hague12
 P. Hague and C. MacCormick,  "Quantum simulation of electron-phonon interactions in strongly deformable materials", New J. Phys. 14,  033019 (2012).

Mezzacapo12
A. Mezzacapo, J. Casanova, L. Lamata, E. Solano, "Digital Quantum Simulation of the Holstein Model in Trapped Ions", Phys. Rev. Lett. 109, 200501  (2012).

Bissbort13
U. Bissbort, D. Cocks, A. Negretti, Z. Idziaszek, T. Calarco, F. Schmidt-Kaler, W. Hofstetter, and R. Gerritsma, "Emulating Solid-State Physics with a Hybrid System of Ultracold Ions and Atoms", Phys. Rev. Lett. 111, 080501 (2013).

Marcin18
M. Płodzień, T. Sowiński, S. Kokkelmans, "Simulating polaron biophysics with Rydberg atoms", Sci. Rep. 8, 9247 (2018).



Sato15 K. Shiozaki, M. Sato, and K. Gomi, “Z_2
topology in nonsymmorphic crystalline insulators: Möbius twist
in surface states”, Phys. Rev. B 91, 155120
(2015).

Fu07 Liang Fu and C. L. Kane, “Topological insulators
with inversion symmetry”, Phys. Rev. B 76, 045302 (2007).

Brz14W. Brzezicki, J. Dziarmaga, and A. M. Oleś,
“Topological Order in an Entangled SU(2)⊗XY
Spin-Orbital Ring”, Phys. Rev. Lett. 112,
117204 (2014).

Top17A. Topp, R. Queiroz, A. Grüneis, L. Müchler, A. W. Rost, A. Varykhalov, D. Marchenko, 
M. Krivenkov, F. Rodolakis, J. L. McChesney, B. V. Lotsch, L. M. Schoop, and C. R. Ast, 
“Surface Floating 2D Bands in Layered Nonsymmorphic Semimetals: ZrSiS and Related Compounds 
”, Phys. Rev. X 7, 041073 (2017).


Brz18 W. Brzezicki, M. Wysokiński, and T. Hyart,
“Topological properties of multilayers and surface
steps in the SnTe material class”, Phys. Rev. B 100, 121107 (2019).

Mazur
G. P. Mazur, K. Dybko, A. Szczerbakow, A. Kazakov, M. Zgirski, E. Lusakowska, S. Kret, J. Korczak, T. Story, M. Sawicki, T. Dietl, 
“Experimental search for the origin of zero-energy modes in topological materials”,Phys. Rev. B 100, 041408(R)
(2019).

Hou07
C.-Y. Hou, C. Chamon, and C. Mudry, “Electron Fractionalization in Two-Dimensional Graphenelike Structures”, Phys. Rev. Lett. 98, 186809 (2007).

















§ HAMILTONIAN AND ITS SYMMETRIES


The SSG Hamiltonian on a square lattice has a form

    H    =    t_x∑_j⃗ (c_j⃗^†c_j⃗+x̂+h.c)+  t_y∑_j⃗ (c_j⃗^†c_j⃗+ŷ+h.c.) 
          +t_d∑_s=±1∑_j⃗(-1)^j_x+j_y+1(c_j⃗^†c_j⃗+x̂+sŷ+h.c.)
          + m∑_j⃗(-1)^j_x+j_y+1c_j⃗^†c_j⃗

where m is the mass term, t_x
and t_y are hopping amplitudes along x
and y directions and t_d is the diagonal hopping amplitude.
The k-space form is given by

    H_k    =     -mσ_zτ_z+2t_xcosk2(cosk2σ_x-sink2σ_y)
        +     t_yτ_x+2t_dcosk2(sink2σ_x+cosk2σ_y)τ_y,

where τ_α operators are given by matrices

    τ_x=([ 0 1 0 0 ⋯; 1 0 1 0 ⋯; 0 1 0 1; 0 0 1 0; ⋮ ⋮     ⋱ ]), τ_y=([  0 -i  0  0  ⋯;  i  0  i  0  ⋯;  0 -i  0 -i;  0  0  i  0;  ⋮  ⋮        ⋱ ]),

and

    τ_z=([  1  0  0  0  ⋯;  0 -1  0  0  ⋯;  0  0  1  0;  0  0  0 -1;  ⋮  ⋮        ⋱ ]), τ_a=([ ⋯ 0 0 0 1; ⋯ 0 0 1 0;   0 1 0 0;   1 0 0 0; ⋰     ⋮ ⋮ ]).

Here we defined additional matrix τ_a which is needed to construct
some of the symmetry operators.

Depending on system width L_y being even or odd the system have
different symmetry properties. However, some symmetries are common for
both cases. The one which is most relevant here is the non-symmorphic
(NS) chiral symmetry defined as S_k=sink2σ_xτ_z+cosk2σ_yτ_z
which satisfies S_k H_kS_k^-1=- H_k.
The k-dependence in S_k is intrinsic and follows from the
half lattice translation that is needed to go from one sublattice
to the other. We also have a time-reversal
symmetry for spinless particles T H_k T^-1= H_-k, where  T= K is complex
conjugation. Finally,
for any L_y we have a symmetry with respect to a mirror line
perpendicular to x direction, passing through a lattice site, taking
a form of M_x=cosk2-isink2σ_z and
acting as M_x H_kM_x^-1= H_-k. Despite the
k-dependence this is a symmorphic symmetry. By shifting a mirror
line to cut a middle of a bond we can also get a chiral mirror symmetry
M_x=σ_yτ_z yielding relation M_x H_kM_x^-1=- H_-k.

For odd L_y we have another mirror symmetry with respect to line perpendicular to y direction M_y H_kM_y^-1= H_k, where
M_y=τ_a.
For even L_y mirror M_y does not exist but we have a particle-hole
symmetry C=i Kσ_zτ_aτ_z and inversion
symmetry I=σ_xτ_a, yielding relations C H_k C^-1=- H_-k
and I H_kI^-1= H_-k.

It is important to notice that in the eigenbasis of  S_k the mirror symmetry operator M_x 
operator transforms to identity. This is possible because if we put eigenvectors of  S_k in the columns
of unitary matrix U_k then  M_x is transformed as M̃_x= U_-k^†M_x U_k.
Note that this is not similarity transformation so the spectrum of M_x is not left invariant. The form of 
transformation is dictated by the mirror-symmetry relation with the Hamiltonian in the new basis. Denoting
H̃_k= U_k^†M_x U_k we get that M̃_xH̃_kM̃_x^-1=H̃_-k.



§ HIDDEN SYMMETRIES


To see the hidden symmetries we first transform H_k to H'_k =  R_k/2^† H_kR_k/2 using
R_k/2=exp(ik4σ_z) to get

    H'_k = -mσ_zτ_z+2t_xcosk2σ_x
     +  2t_dcosk2σ_yτ_y+t_yτ_x.

The eigenfunctions of τ_x corresponding to eigenvalues ε_n are ψ_n(j)=√(2L_y+1)sinnjπL_y+1,
where j labels sites in the y direction and n labels modes
(j,n=1,2,…,L_y). We can use these transverse modes to construct a new basis 
 |ϕ_2n-1⟩ =(|ψ_n⟩ +|ψ_L_y+1-n⟩)/√(2)
and |ϕ_2n⟩ =(|ψ_n⟩ -|ψ_L_y+1-n⟩)/√(2)
for n=1,2,…,⌊ L_y/2⌋ and if L_y
is odd |ϕ_0⟩ =|ψ_⌈ L_y/2⌉⟩. In this basis τ_x, τ_y and τ_z have block-diagonal
forms with ⌊ L_y/2⌋ diagonal blocks given
by

    τ_x,n=ε_nσ'_x,  τ_y,n=ε_nσ'_y,  τ_z,n=σ'_z,

where σ'_α is a new set of Pauli matrices, and for odd
L_y the block n=0 is given
by τ_x,0=τ_y,0=0 and τ_z,0=1. Thus
the Hamiltonian  has a block-diagonal form

    H'_k,n    =     -mσ_zσ'_z+2t_xcosk2σ_x
        +     2t_dε_ncosk2σ_yσ'_y+t_yε_nσ'_x

supporting the hidden symmetries discussed in the main text.



§ HIDDEN SYMMETRIES IN A TWO-DIMENSIONAL LIMIT 


In case when the system is periodic both in x and y direction it can be described by a 
k-space Hamiltonian written in a form of


    H_k⃗    =     -mσ_zτ_z+2t_xcosk_x2(cosk_x2σ_x-sink_x2σ_y)
        +     2t_ycosk_y2(cosk_y2τ_x-sink_y2τ_y)
        +     4t_dcosk_x2cosk_y2(sink_x2σ_x+cosk_x2σ_y)
          (sink_y2τ_x+cosk_y2τ_y),


where σ_α and τ_α are the Pauli matrices describing the unit cell 
along x and y directions. Similarly as before, to see the hidden symmetries we transform
H_k⃗ to H'_k⃗ =  R_k⃗^† H_k⃗R_k⃗ using
R_k⃗=exp(i4[k_xσ_z+k_yτ_z]) to get


    H'_k⃗   =    -mσ_zτ_z+2t_xcosk_x2σ_x+2t_ycosk_y2τ_x
       +    4t_dcosk_x2cosk_y2σ_yτ_y.


Obviously H'_k⃗ (and already  H_k⃗) is invariant under interchange
of σ and τ operators if 


    t_xcosk_x2=t_ycosk_y2.


The spin-interchange X_12σ⃗X_12=σ⃗' and vice-versa is realized 
by  operator  X_12=12(1+σ⃗·τ⃗). Coming back to the original
basis of Eq.(<ref>) we find that X_12 operator gives rise to two different symmetry operations,


    M_xy   ≡    R_k_y,k_xX_12R_k_x,k_y^†=12(1+σ⃗·τ⃗), 
     X_k⃗   ≡    R_k_x,kyX_12R_k_x,k_y^†=12(1+cosk_x-k_y2[σ_xτ_x+σ_yτ_y]  
       +   sink_x-k_y2[σ_xτ_y-σ_yτ_x]+σ_zτ_z).


The first one is the mirror symmetry with respect to the x̂+ŷ line. One can verify that 
if t_x=t_y then 


    M_xy H_k_x,k_yM_xy^-1= H_k_y,k_x.


Therefore M_xy commutes with H_k⃗ for k_x=k_y. Note that this also holds 
for a finite periodic system where L_x≠L_y,  when the discretized quasimomentum satisfies k_x=k_y.
On the other hand, the hidden symmetry operator X_k⃗  satisfies


    X_k_x,f(k_x) H_k_x,f(k_x) X_k_x,f(k_x)^-1= H_k_x,f(k_x),


as long as condition (<ref>) holds which means that 


    k_y=f(k_x)=2arccos(t_x/t_ycosk_x/2).


Note that this does not require that t_x=t_y. Hence, the hidden symmetry is a different operator than
the mirror symmetry M_xy eventhough it originates from the same X_12 operator in the
R_k⃗–transformed basis. It can be regarded as an interchange of longitudinal and transverse modes
in the system. Note that condition (<ref>) trivially generalizes to the condition found for a multimode wire. 




§ ROBUSTNESS OF THE DW STATES




To check whether the DW states protected by the NS chiral symmetry are indeed topologically protected we 
consider a single-mode SSG wire. An arbitrary perturbation that preserve the NS chiral symmetry S_k 
has a form of


    V_k =  a_k σ_z + b_k(cosk2σ_x-sink2σ_y),


where k-dependent coefficients should have such form that V_k is 2π-periodic in k. Hence we have


    a_k    =    p_0 + ∑_i=1^r(p_2i-1sin(ik)+p_2icos(ik)),
    
    b_k    =   ∑_i=1^r(q_2i-1sin(2i-1)k2+q_2icos(2i-1)k2),


where r is  the range of hoppings involved in V_k. 

In Fig. <ref> we show the scaling of the energy of the DW state between two regions of 
opposite masses for a single-mode wire in presence of perturbation V_k. We choose
r=10 and we pick random p_n and q_n coefficients from the uniform distribution over
the interval [-0.1t_x,0.1t_x]. The masses are chosen as m_1,2=± t_x  so that such perturbation
does not close the bulk gap. As we can see in Fig. <ref>(a) the energy E_0 of the DW state still scales
to zero as 1/w. On the other hand, if we add a perturbation that breaks the NS chiral symmetry, for
instance V_k = δσ_x, the DW state should converge to a non-zero energy. In
Fig. <ref>(b) we see that this is indeed the case for our single-mode wire E_0 scales like 1/w to
δ.



§ LOCALIZATION OF THE DW STATES




In Fig. <ref> we show the spectral flow and the local density of states for a DW state in a representative multimode chain with L_y=7 for increasing
width w of a domain wall between two domains with different NS ℤ_2 invariant. We note that the DW state is exponentially localized around
the center of the domain wall at x=x_0 (defined by m(x_0)=0) and the localization length decreases with increasing w .




§ FERMI VELOCITY AT DIRAC POINTS K=± K_N


We notice that the Hamiltonian H'_k,n given by Eq. (<ref>)
commutes with P=σ_xσ'_x. The bands that cross
at k=± k_n can be found in the block P=-1 of the
Hamiltonian that takes the form,


    H'_k,n,-    =     -4sink+k_n/4sink-k_n/4(t_xσ_x+t_dε_nσ_z)
        -    (m-m_n)σ_z.

The remaining P=+1 block is related by a shift of 2π
in the k-space and unitary transformation, namely H'_k,n,+=σ_z H'_k+2π,n,-σ_z.
Consequently, the Fermi velocity at the Dirac points is given by

    v_n=± t_x√((1-t_y^2/4t_x^2ε_n^2)(1+t_d^2/t_x^2ε_n^2)).




§ OVERLAP OF TWO DOMAIN-WALL FUNCTIONS


We assume that the functional form of m(x) is such that it only depends on parameter x/w. In analogy to the domain wall solution for a gap closing point at
k=π we get two solutions for k=± k_n gap closings in a
form (we ignore the spinor structure which would change as a function momentum and the normalization factor which is not important for the statements below)

    ψ^±(x)=exp[± ik_nx-v_n^-1∫_0^x m(x'/w)dx'].

Their overlap is

    ⟨ψ^-|ψ^+⟩=∫_-∞^+∞exp[2ik_nx-2v_n^-1∫_0^x m(x'/w)dx']dx

and substituting x=yw we get

    ⟨ψ^-|ψ^+⟩=w∫_-∞^+∞exp[2ik_nwy-2wv_n^-1M(y)]dy,

where M(y)=∫_0^y m(y')dy'. Assuming that M(y)>0,  m(∞)>0 and m(-∞)<0 as we expect from a domain wall, we
notice that   f(y)=exp[-2v_n^-1 w M(y)] is a Schwartz function. This property does not depend on the details of the model (even if one goes beyond the linear order expansion in momentum) or the shape of the DW as long as the solutions ψ_± are smooth functions decaying exponentially (or faster) far away from the DW.
The Fourier transform of a Schwartz function is also a Schwartz function and therefore the overlap ⟨ψ^-|ψ^+⟩
must vanish quicker than any power of 1/(k_nw). 

To see the exponential dependence explicitly we can assume for example m(x/w)=m_0 x/w. Then we obtain 

    ⟨ψ^-|ψ^+⟩   =   w∫_-∞^+∞exp[2ik_nwy-wm_0v_n^-1y^2]dy 
       =   w √(π v_n/wm_0)exp[-w k_n^2 v_n/m_0].






§ CHARGE OF A DOMAIN WALL


To understand the charge of the DW we start by calculating the charge appearing at the end of single-mode SSG chain with open boundary conditions

    H=([ -m  t  0  0  ⋯;  t  m  t  0  ⋯;  0  t -m  t  ⋯;  0  0  t  m;  ⋮  ⋮  ⋮     ⋱ ]).

Assuming even number of sites L=2N the Hamiltonian can be Block diagonalized into 2x2 blocks similarly as in Appendix <ref>, and the eigenstates can be found by diagonalizing these blocks. This way we find that the eigenenergies are (n=1,2,....,N)

    E_n ± = ± E_n,  E_n=√(m^2+4 t^2 cos^2(n π/2N+1)).

The eigenstate corresponding to E_n- at lattice site j is 

    ψ_n, -(j)   =   sin(n j π/2N+1)/√(2(2N+1))[[1-(-1)^j] √(1+m/E_n)
          -[1+(-1)^j] √(1-m/E_n) ].

The total end charge q_end at lattice sites j=1,...,2 ξ (relative to the corresponding bulk charge ξ) is

    q_end   =   2/2N+1∑_n=1^N ∑_j=1^2ξsin^2(n j π/2N+1) [1-(-1)^j m/E_n]-ξ
       =   1/2N+1∑_n=1^N ∑_j=1^2ξ (-1)^j cos(n j 2π/2N+1)  m/E_n
       =    -ζ/4+1/21/2N+1∑_n=1^N  cos[(2ξ+1/2) n  2π/2N+1]/cos[1/2n  2π/2N+1]m/E_n
       =    sign(m)-ζ/4,

where

    ζ(m)    =   2/2N+1∑_n=1^N m/E_n =  2/2π∫_0^π dk  m/√(m^2+4 t^2 cos^2(k/2))
       =   2/πm/√(m^2+4t^2) K(4t^2/m^2+4t^2)

is the difference of the bulk filling factors of the two sublattices and

    K(x)= ∫_0^π/2 dθ1/√(1-x sin^2(θ))

is the complete elliptic integral of the first kind. Here we have used 

    1/4π∫_0^π dk  m cos[(2ξ+1/2) k ]/√(m^2+4 t^2 cos^2(k/2))cos[1/2k ]  = sign(m)/4,

which is valid up to corrections which decay exponentially with increasing ξ. At the other end of the chain there is a charge -q_end. 



Let's now consider a sharp DW between regions with mass m_i and m_j. If we first assume that we turn off the hopping connecting these two regions, we find that the region with mass m_i gives rise to the charge -q_end, i and the region with mass m_j gives rise to a charge q_end, j at the DW so that the total charge is

    q_DW= sign(m_j)- sign(m_i)/4-δ_ji/2,

where we have defined

    δ_ji=ζ(m_j)-ζ(m_i)/2.

This charge is exponentially localized at the DW and therefore when the hopping connecting the regions is turned on it causes only a local perturbation in the Hamiltonian (slightly redistributing the charge density locally) but does not influence the total charge localized at the DW.
If both m_i and m_j have the same sign the regions are in the same topological phase and the q_DW=-δ_ji/2. If the signs of the masses are different we have DW between two topologically distinct phases. We now consider two possible DWs: (i) Soliton where mass changes from m_1>0 to m_2<0 as a function of increasing x and (ii) antisoliton where mass changes from m_2<0 to m_1>0 as a function of increasing x (see Figs. <ref> and <ref> in the main text.)  We denote δ=δ_21.

(i) In the case of soliton the sharp DW carries a charge  q_DW=-(1+δ)/2. By increasing the width of the DW we find that the spectral flow is such that a state will approach zero energy from positive energies (see Fig. <ref> in the main text). This means that if this zero-energy state is unoccupied the charge of the DW is q_DW and if it is occupied the charge is q_DW+1. Thus we can summarize that the possible charges for soliton are 

    q_DW=± 1 -δ/2.


(ii) For antisoliton the charge of the sharp DW is  q_DW=(1+δ)/2. By increasing the width of the DW we find that the spectral flow is such that a state will approach zero energy from negative energies. This means that if this zero-energy state is unoccupied the charge of the DW is q_DW-1 and if it is occupied the charge is q_DW. Thus we can summarize that the possible charges for antisoliton are 

    q_DW=± 1 +δ/2.


In the case of multimode SSG system we can separate the transverse modes n as discussed in Appendix <ref>. We consider the cases (a) L_y even and (b) L_y odd separately.

(a) When L_y is even the Hamiltonian can be decomposed into 4x4 blocks given by Eq. (<ref>). Each of these blocks n=1,...L_y/2 supports symmorphic chiral symmetries
C_yz=σ_y σ_z' and C_zy=σ_z σ_y'. Therefore using the argument given in Ref. <cit.> we find that each transverse mode n carries possible charges q_n=-N_n/2,...,N_n/2, where N_n is the number of zero-energy states at the DW supported by transverse mode n and the value of charge q_n is determined by the number of occupied zero-energy states. In the case of smooth DWs each transverse mode supports N_n=2 zero-energy states if ν_n is different on the two sides of the DW  so that q_n=-1,0,1 (with two possible states corresponding to q_n=0). If ν_n is the same on both sides then N_n=0 and q_n=0. The total charge of the DW is 

    q_DW=∑_n=1^L_y/2 q_n.


(b) When L_y is odd the Hamiltonian can be decomposed into ⌊ L_y/2⌋ 4x4 blocks obeying the same symmorphic chiral symmetries. Each of these modes  carries charges q_n=-1,0,1 (n=1,...⌊ L_y/2⌋). Additionally there exists one transverse mode n=0 which is similar as the one studied in the case of single-mode SSG chain. Thus, this transverse mode carries a charge
q_0=(± 1 -δ)/2 in the case of solitons and q_0=(± 1 + δ)/2 in the case of antisolitons. The total charge of the DW is

    q_DW=∑_n=0^⌊ L_y/2⌋ q_n.




