
\section{The tight-binding Hamiltonian of the model system}
\label{S-1}

\begin{figure}
\centering
\includegraphics[width=0.7\textwidth]{Figs1.pdf}
\caption{
  (a)  The spatial distribution of the confining potential $V^{}_{\rm DQD-in}(x)$ that determines the intra-cell tunneling amplitude within the double quantum-dot chain in the main text (see Fig.~1). (b) Confining potential $V^{}_{\rm DQD-ex}(x)$ that determines the inter-cell tunneling amplitude. (c) The intra/inter-cell hopping amplitude of an InAs double quantum dot chain, $t^{}_{\rm in/ex,0}$, as a function  of the potential difference $\Delta V=V^{}_{1}-V^{}_{2}$, with the  periodic length $2L=120$~nm and  $d=10$~nm. Here, $t^{}_{\rm in,0}$ ($t^{}_{\rm ex,0}$) is denoted by the solid (dashed) line and, meanwhile, the variations of the  chemical potential $V_{\rm c}$ are also displayed by the dash-doted curve.  (d)   The intra- and inter-cell spin-flip tunneling amplitudes, $t^{\prime}_{\rm in}$ and  $t^{\prime}_{\rm ex}$, of the double-quantum-dot chain as functions of the tilting angle $\theta$, with the potential difference $\Delta V^{}_{}=-10$~mev,  the  magnetic field strength $B=1.0$~T, and the spin-orbit coupling scaled by $x^{}_{\rm so}\equiv \hbar/(m_{e}\alpha)=180$~nm.}
\label{Figs1}
\end{figure}





 Here is outlined the  derivation of the  nearest-neighbor tight-binding (TB) Hamiltonian of our model.
 The derivation is based on the original Hamiltonian, Eq.~(5) in the main text. The periodic gate potential there, $U(x)$,  is viewed as composed of two alternating quantum dots, QD1 and QD2, as depicted in Fig.~\ref{Figs1}(a) and (b). Each quantum dot obeys the Hamiltonian
 \begin{align}
 H^{}_{\rm 1/2}=\frac{p^{2}}{2m_{e}}+ V^{}_{\rm qd,1/2}(x)+\alpha p\sigma_{y}+\frac{\Delta_{\rm z}}{2}\sigma^{}_{n},
 \label{H12}
 \end{align}
 where  $\alpha$ is  the strength of the Rashba spin-orbit coupling (SOC), $\Delta_{\rm z}$ is the Zeeman splitting, $\sigma^{}_{n}= \hat{\boldsymbol{\sigma}}\cdot \hat{\mathbf{n}}$ with  $\hat{\mathbf{n}}=\{0, \cos\theta, \sin\theta\}$ being the direction vector of the external magnetic field  and  $\hat{\boldsymbol{\sigma}}$ denotes  the Pauli matrices,  $\hat{\boldsymbol{\sigma}}=\{\sigma^{}_{x},\sigma^{}_{y},\sigma^{}_{z}\}$.  The  confining  potentials that define the quantum dots are
 \begin{align}
 V^{}_{\rm qd,1}(x)=\begin{cases}
 V_{1} ~~~~~~~~~~~~~x\geq -d\\
 0~~~~~~~ -L+d\leq x<-d,\\
 V_{2}~~~~~~~~~~~x<-L+d
 \end{cases}~~~~   V^{}_{\rm qd,2}(x)=\begin{cases}
 V_{1} ~~~~~~~~~~~~~x< d\\
 0~~~~~~~ d\leq x<L-d.\\
 V_{2}~~~~~~~~x\geq L-d
 \end{cases}
 \label{lcp}
 \end{align}

 We next solve for the energy eigenstates of the two dots using Eq. (\ref{H12}). Moreover,  in order to preserve the locality of the wave functions, we assume  the electron energy eigenvalue   $E<V^{}_{1},V^{}_{2}$. It
is expedient to consider the Schr\"{o}dinger equation pertaining to all three regions in Eq. (\ref{lcp}) by introducing
 \begin{align}
 \widetilde{H}^{}_{\zeta=0,1,2}\psi^{}_{\zeta }(x)=E\psi_{\zeta}(x).
 \label{WQD}
 \end{align}
 Here,   $\psi^{}_{\zeta}$ is the corresponding eigenfunction on each portion  and
\begin{align}
\widetilde{H}^{}_{\zeta}=\frac{p^{2}}{2m_{e}}+\alpha p\sigma^{}_{y}+\frac{\Delta^{}_{\rm z}}{2}\sigma^{}_{n}+V^{}_{\zeta},
\label{LQD}
 \end{align}
  where
 $V_{0}=0$, and $V_{1,2}$ are given in Eqs. (\ref{lcp}).    
The eigenfunctions of $\widetilde{H}^{}_{\zeta}$ are spinors,
 \begin{align}
 \psi^{}_{\zeta}(x)=\exp\left(i\gamma^{}_{\zeta}x\right)\begin{pmatrix}
 w_{\zeta}\\
 d_{\zeta}
 \end{pmatrix}.
 \label{BEFs}
 \end{align}
For a fixed value of $E$ there exist four independent solutions to Eq.~(\ref{WQD}), with the corresponding wave vectors  $\gamma^{}_{\zeta, \lambda=1,2,3,4}$  determined by the  quartic equation
  \begin{align}
  \left(\frac{\hbar^{2}{\gamma^{2}_{\zeta}}}{2m_{e}}+V^{}_{\zeta}-E\right)^{2}-\hbar^{2}\gamma^{2}_{\zeta}\alpha^{2}-\Delta_{\rm z}\hbar \gamma^{}_{\zeta}\alpha\cos(\theta)-\frac{\Delta^{2}_{\rm z}}{4}=0.
  \label{wk}
  \end{align}
  The respective spinor components are
 \begin{align}
w^{}_{\zeta,\lambda}= i\hbar\alpha\gamma^{}_{\zeta, \lambda}+i\frac{\Delta^{}_{\rm z}}{2}\cos(\theta),~~~~~~d^{}_{\zeta,\lambda}=\frac{\hbar^{2}\gamma^{2}_{\zeta, \lambda}}{2m_{e}}+\frac{\Delta_{\rm z}}{2}\sin(\theta)+V^{}_{\zeta}-E.
  \label{wkc}
 \end{align}
 In particular, when $\theta=\pi/2$, the expressions for the wave vectors can be derived analytically
 \begin{align}
 \gamma^{}_{\zeta,1/2}=&\pm\sqrt{\frac{2m_{e}}{\hbar^{2}}}\sqrt{E-V_{\chi}+m_{e}\alpha^{2}+\sqrt{m^{2}_{e}\alpha^{4}+2m_{e}\alpha^{2}(E-V_{\chi})+\Delta^{2}_{\rm z}/4}} \nonumber\\
 \gamma^{}_{\zeta,3/4}=&\pm\sqrt{\frac{2m_{e}}{\hbar^{2}}}\sqrt{E-V_{\chi}+m_{e}\alpha^{2}-\sqrt{m^{2}_{e}\alpha^{4}+2m_{e}\alpha^{2}(E-V_{\chi})+\Delta^{2}_{\rm z}/4}}.

 \end{align}
The convergence of the wave function at  $|x|\rightarrow \infty$ requires the wave vectors to have ${\rm Im}[\gamma^{}_{1}]>0$  for dot 1 and ${\rm Im}[\gamma^{}_{2}]>0$ for dot 2.  Meanwhile, the convergence condition when $x\rightarrow -\infty$ demands  ${\rm Im}[\gamma^{}_{2}]<0$  for dot 1 and ${\rm Im}[\gamma^{}_{1}]<0$ for dot 2.
Then, the energy eigenfunctions of QD1 and QD2 can be written as
\begin{align}
\Psi_{1}(x)=&\sum^{4}_{\lambda=1}c^{}_{0,\lambda}[1-\Theta(x+d)]\Theta(x+L-d)\psi^{}_{0,\lambda}(x)+\Theta(x+d)[c^{}_{1,1}\psi^{}_{1,1}(x)+c^{}_{1,4}\psi^{}_{1,4}(x)]\nonumber\\
&~~~~~~~+[1-\Theta(x+L-d)][c^{}_{2,2}\psi^{}_{2,2}(x)+c^{}_{2,3}\psi^{}_{2,3}(x)]\\
\Psi_{2}(x)=&\sum^{4}_{\lambda=1}v^{}_{0,\lambda}[1-\Theta(x-L+d)]\Theta(x-d)\psi^{}_{0,\lambda}(x)+[1-\Theta(x-d)][v^{}_{1,2}\psi^{}_{1,2}(x)+v^{}_{1,3}\psi^{}_{1,3}(x)]\nonumber\\
&~~~~~~~+\Theta(x-L+d)[v^{}_{2,1}\psi^{}_{2,1}(x)+v^{}_{2,4}\psi^{}_{2,4}(x)],
\end{align}
where $\Theta(x)$ is the Heaviside step function  $\Theta(x)=\begin{cases}
0  ~~x\leq0\\
1 ~~x>0\end{cases}$, the coefficients $c^{}_{\zeta,\lambda}$ and  $v^{}_{\zeta,\lambda}$  are determined by the boundary conditions, i.e., from the continuity  equations at $x^{}_{1/3}=\mp d$ and $x^{}_{2/4}=\mp L\pm d$.
 Using the explicit form for $V^{}_{qd,1}(x)$ in Eq.~(\ref{lcp}), the boundary conditions for the localized wave function on QD1 are \cite{com}
\begin{align}
\lim_{\epsilon \rightarrow 0 }[\Psi_{1}(x_{1}+\epsilon)-\Psi_{1}(x_{1}-\epsilon)]=0 ,~~~~\lim_{\epsilon\rightarrow 0} [\Psi_{1}(x_{2}+\epsilon)-\Psi_{1}(x_{2}-\epsilon)]=0 , \nonumber\\
\lim_{\epsilon \rightarrow 0 }[\Psi'_{1}(x_{1}+\epsilon)-\Psi'_{1}(x_{1}-\epsilon)]=0 ,~~~~\lim_{\epsilon\rightarrow 0} [\Psi'_{1}(x_{2}+\epsilon)-\Psi'_{1}(x_{2}-\epsilon)]=0.
\label{CD-1}
\end{align}
By regarding $\hat{z}_{}^{T}= \{c^{}_{0,1},c^{}_{0,2},c^{}_{0,3},c^{}_{0,4},c^{}_{1,1},c^{}_{1,4},c^{}_{2,2},c^{}_{2,3}\}$ as the variable vector, Eq. (\ref{CD-1}) can be written in a matrix form
\begin{align}
\mathbf{M} \cdot \hat{z}=0,
\label{MAE}
\end{align}
where
\begin{footnotesize}
\begin{align}
\mathbf{M} =\begin{pmatrix}
\psi^{}_{0,1,1}(x_{1})&\psi^{}_{0,2,1}(x_{1})&\psi^{}_{0,3,1}(x_{1})&\psi^{}_{0,4,1}(x_{1})&-\psi_{1,1,1}(x_{1})&-\psi_{1,4,1}(x_{1})&0&0\\
\psi^{}_{0,1,2}(x_{1})&\psi^{}_{0,2,2}(x_{1})&\psi^{}_{0,3,2}(x_{1})&\psi^{}_{0,4,2}(x_{1})&-\psi_{1,1,2}(x_{1})&-\psi_{1,4,2}(x_{1})&0&0\\
\psi^{}_{0,1,1}(x_{2})&\psi^{}_{0,2,1}(x_{2})&\psi^{}_{0,3,1}(x_{2})&\psi^{}_{0,4,1}(x_{2})&0&0&-\psi_{2,2,1}(x_{2})&-\psi_{2,3,1}(x_{2})\\
\psi^{}_{0,1,2}(x_{2})&\psi^{}_{0,2,2}(x_{2})&\psi^{}_{0,3,2}(x_{2})&\psi^{}_{0,4,2}(x_{2})&0&0&-\psi_{2,2,2}(x_{2})&-\psi_{2,3,2}(x_{2})\\
\gamma^{}_{0,1}\psi^{}_{0,1,1}(x_{1})&\gamma^{}_{0,2}\psi^{}_{0,2,1}(x_{1})&\gamma^{}_{0,3}\psi^{}_{0,3,1}(x_{1})&\gamma^{}_{0,4}\psi^{}_{0,4,1}(x_{1})&-\gamma^{}_{1,1}\psi_{1,1,1}(x_{1})&-\gamma^{}_{1,4}\psi_{1,4,1}(x_{1})&0&0\\
\gamma^{}_{0,1}\psi^{}_{0,1,2}(x_{1})&\gamma^{}_{0,2}\psi^{}_{0,2,2}(x_{1})&\gamma^{}_{0,3}\psi^{}_{0,3,2}(x_{1})&\gamma^{}_{0,4}\psi^{}_{0,4,2}(x_{1})&-\gamma^{}_{1,1}\psi_{1,1,2}(x_{1})&-\gamma^{}_{1,4}\psi_{1,4,2}(x_{1})&0&0\\
\gamma^{}_{0,1}\psi^{}_{0,1,1}(x_{2})&\gamma^{}_{0,2}\psi^{}_{0,2,1}(x_{2})&\gamma^{}_{0,3}\psi^{}_{0,3,1}(x_{2})&\gamma^{}_{0,4}\psi^{}_{0,4,1}(x_{2})&0&0&-\gamma^{}_{2,2}\psi_{2,2,1}(x_{2})&-\gamma^{}_{2,3}\psi_{2,3,1}(x_{2})\\
\gamma^{}_{0,1}\psi^{}_{0,1,2}(x_{2})&\gamma^{}_{0,2}\psi^{}_{0,2,2}(x_{2})&\gamma^{}_{0,3}\psi^{}_{0,3,2}(x_{2})&\gamma^{}_{0,4}\psi^{}_{0,4,2}(x_{2})&0&0&-\gamma^{}_{2,2}\psi_{2,2,2}(x_{2})&-\gamma^{}_{2,3}\psi_{2,3,2}(x_{2})
\end{pmatrix},
\label{MES}
\end{align}
\end{footnotesize}
with the two spin components $\psi^{}_{\zeta,\lambda,1}(x)=w^{}_{\zeta,\lambda}\exp(i\gamma^{}_{\zeta,\lambda} x)$ and $\psi^{}_{\zeta,\lambda,2}(x)= d^{}_{\zeta,\lambda}\exp(i\gamma^{}_{\zeta,\lambda} x)$. Obviously, the nontrivial solution to Eq.~(\ref{MES}) requires
\begin{align}
\det[ \mathbf{M}(E)]=0,
\label{DEM}
\end{align}
which yields an equation for the energy eigenvalues $E$, see Eqs.~(\ref{wk})  and (\ref{wkc}).
 Let $E_{ \Uparrow}$ and $E_{ \Downarrow}$  denote the two lowest energy solutions to Eq.~(\ref{DEM}), with $E_{ \Uparrow}\geq E_{ \Downarrow}$. The explicit expressions for the two corresponding spinors can also be completely determined by the boundary conditions   Eqs.~ (\ref{CD-1}) in conjunction with the normalization condition $\int \Psi^{\dagger}_{1,\tau=\Uparrow,\Downarrow}(x)\Psi^{}_{1,\tau}(x)dx=1$.
where $\tau=\Uparrow,~\Downarrow$ indicate the two quasi-spin eigenstates on the dot.


The two lowest-energy spin states of the other quantum dot, QD2,
i.e., $\Psi^{}_{2, \Uparrow}(x)$ and $\Psi^{}_{2,\Downarrow }(x)$,
 can be derived similarly to the discussion above. Due to the inversion symmetry  which implies  $V^{}_{\rm qd,1}(x)=V^{}_{\rm qd,2}(-x)$, the corresponding energy eigenvalues are identical to  those of QD1.


 Exploiting the localized wave functions on each quantum dot, we next derive the tunneling amplitudes of the double-quantum-dot (DQD) model, utilizing linear combinations of the atomic orbitals.
The intra-cell amplitude is expediently found by considering  the interdot tunneling in the double
 quantum-dot Hamiltonian [see Fig.~\ref{Figs1}(a)]
 \begin{align}
 H^{}_{\rm DQD-in}(x)=\frac{p^{2}}{2m_{e}}+\alpha p\sigma_{y}+\frac{\Delta^{}_{\rm z}}{2}\sigma_{n}+V^{}_{\rm DQD-in}(x)-V^{}_{\rm c},
 \label{ODQD1}
 \end{align}
with the double-well potential
\begin{align}
 V^{}_{\rm DQD-in}(x)=V^{}_{1}[1-\Theta(|x|-d)]+V^{}_{2}\Theta(|x|-L+d)
\end{align}
and the constant chemical potential $V^{}_{\rm c}=\frac{1}{2}(E^{}_{\Uparrow}+E^{}_{\Downarrow})$ to offset the on-site orbital energy (which is conventionally set at zero for the "standard" SSH model).
To this end, one first derives the orthonormal basis $\left\{\Phi^{}_{1,\Uparrow}(x),\Phi^{}_{1,\Downarrow}(x),\Phi^{}_{2,\Uparrow}(x),
\Phi^{}_{2,\Downarrow}(x)\right\}$ using the four localized states $\Psi^{}_{1/2,\Uparrow}(x)$ and  $\Psi^{}_{1/2,\Downarrow}(x)$. Those are then used to present the electron field operators as
$
\Psi_{\rm e }(x)= \sum^{}_{\sigma=\Uparrow,\Downarrow}[a^{}_{\sigma}\Phi^{}_{1,\sigma}(x) + b^{}_{\sigma}\Phi^{}_{2,\sigma}(x)]$,
 where $a^{}_{\sigma}$ ($a^{\dagger}_{\sigma}$) and $b^{}_{\sigma}$ ($b^{\dagger}_{\sigma}$) are the annihilation (creation) operators for the quasi-spin states on QD1 and QD2,
 respectively. The second quantized form of $H^{}_{\rm DQD-in}$ is then
 \begin{align}
 H ^{}_{\rm DQD-in}=& \int \Psi^{\ast}_{\rm e}(x) H^{}_{\rm DQD-in}(x) \Psi^{}_{\rm e }(x)dx \nonumber\\
=&\frac{\Delta^{}_{\rm z}}{2} \left(a^{\dagger}_{ }\sigma^{}_{z}a^{}_{ }+b^{\dagger}_{ }\sigma^{}_{z}b^{}_{ }\right)+\left(a^{\dagger}t^{}_{\rm in}b+{\rm H.c.}\right)\ ,
\label{DQD-in}
\end{align}
with  the spinors $a^{}_{}=\{a^{}_{\Uparrow},a^{}_{\Downarrow}\}$ and $b^{}_{}= \{b^{}_{\Uparrow},b^{}_{\Downarrow}\}$. In fact, the on-site energies of different quasi-spin states in each dot are modified by the spin-orbit interaction. However, as compared to the large spin-orbit length of the InAs nanowire, i.e., $x^{}_{\rm so}=180$~nm, the small characteristic length of the dots,  i.e., $L-2d=40$~nm, renders the SOC-induced modification negligible, and then it is a good approximation to write  the  quasi-spin up (down) on-site energy as the positive (negative) half Zeeman splitting energy as shown in Eq.~(\ref{DQD-in}). The quasi-spin tunneling amplitudes can be written as a matrix
\begin{align}
t^{}_{\rm in}=t^{}_{{\rm in}, 0}\exp\left[i\varphi^{}_{{\rm in} } (\hat{\mathbf{v}}^{}_{\rm in}\cdot\hat{\boldsymbol{\sigma}})\right]\ ,
\label{T-IN}
\end{align}
with
\begin{align}
[t^{}_{\rm in}]^{}_{\sigma \sigma'}=\int \Phi^{\dagger} _{1,\sigma}(x)H^{}_{\rm DQD-in}\Phi_{2,\sigma' }(x)dx\ .
\label{s19}
\end{align}



As seen, the effect of the  spin-orbit coupling appears as a (matrix) phase factor multiplying the intra-cell spin tunneling matrix.
 Here, $t^{}_{\rm in,0}$ is the real tunneling amplitude dominated by the potential difference $\Delta V=V_{1}-V_{2}$ as shown in Fig.~\ref{Figs1}(c), $\varphi^{}_{\rm in }$ is the spin rotation phase induced by the spin-orbit coupling, and $\hat{\mathbf{v}}^{}_{\rm in}=\{\sin(\vartheta^{}_{\rm in}),0,\cos(\vartheta^{}_{\rm in})\}$ is the unit vector determining the direction of the axis of rotation. Numerically, $\varphi^{}_{\rm in }$ is proportional to the ratio between the interdot distance and the spin-orbit length $L/x^{}_{\rm so}$ with $x^{}_{\rm so}= \hbar/(m^{}_{e}\alpha)$,  and that is consistent with
the analytical analysis in Ref.~\cite{Liu2018}


It is important to note that the coordinate system for defining $\hat{\mathbf{v}}^{}_{{\rm in}}$ spin state
 does not coincide with that defining the spin operators in Eq.~(\ref{ODQD1}).
The $x$-component of $\hat{\mathbf{v}}^{}_{{\rm in}}$ corresponds to the interdot spin-flip tunneling between the localized quasi-spin states, and vanishes in the absence of spin mixing.
Therefore, for zero magnetic field, i.e., when $\Delta^{}_{\rm z}=0$, or when $\theta=0$, the vector $\hat{\mathbf{v}}^{}_{\rm in}$ is purely along the $z$ direction since then $[H^{}_{\rm DQD-in},\sigma^{}_{y}]=0$. In other cases,
the specific value of $\vartheta^{}_{\rm in }$ depends on the degree of spin mixing on each localized quasi-spin state, which in turn is determined by the angle between the external magnetic field and the one induced by the spin-orbit  interaction. The amplitude of the spin-flip tunneling in Eq.~(\ref{T-IN}), $t^{\prime}_{\rm in }=t^{}_{\rm in, 0}\sin(\varphi^{}_{\rm in }) \sin(\vartheta^{}_{\rm in })$, as a function of the tilting angle $\theta$ for a fixed magnitude of the magnetic field, is portrayed in Fig.~\ref{Figs1}(d), showing that it vanishes for $\theta=0$.


The inter-cell hopping is addressed by considering a neighboring DQD [as shown in Fig.~\ref{Figs1}(b)] which obeys the Hamiltonian
\begin{align}
H^{}_{\rm DQD-ex}(x)=\frac{p^{2}}{2m_{e}}+\alpha p \sigma_{y}+\frac{\Delta_{\rm z}}{2}\sigma^{}_{n}+V^{}_{\rm DQD-ex}(x)-V^{}_{\rm c},
\end{align}
where the double-well potential takes the form
\begin{align}
V^{}_{\rm DQD-ex}(x)=V^{}_{2}[1-\Theta(|x|-d)]+V^{}_{1}\Theta(|x|-L+d).
\end{align}
Similar to the intra-cell spin tunneling discussed above, the spin transfer matrix for the inter-cell tunneling  takes the form
\begin{align}
t^{}_{\rm ex}= t^{}_{{\rm ex},0}\exp\left[i\varphi^{}_{{\rm ex} } (\hat{\mathbf{v}}^{}_{\rm ex}\cdot\hat{\boldsymbol{\sigma}})\right]\ ,
\label{T-EX}
\end{align}
with $t^{}_{{\rm ex},0}$ and $\varphi^{}_{{\rm ex}}$ corresponding to the real amplitude and the spin rotation phase of the interdot tunneling, and $\hat{\mathbf{v}}^{}_{{\rm ex}}=\{\sin(\vartheta^{}_{{\rm ex}}),0,\cos(\vartheta^{}_{{\rm ex}})\}$  describing the direction of the axis of rotation.
 Analogously, the corresponding spin-flip tunneling amplitude  $t^{\prime}_{\rm ex}=t^{}_{{\rm ex},0}\sin(\varphi^{}_{\rm ex})\sin(\vartheta^{}_{\rm ex})$ can be regulated by manipulating the tilting angle $\theta$, as shown in Fig.~\ref{Figs1}(d).
Therefore, by including the intra-cell spin tunneling, the second quantized Hamiltonian of the double-quantum-dot chain is
\begin{align}
H^{}_{\rm T}= \sum^{}_{n}\frac{\Delta^{}_{\rm z}}{2}\left(a^{\dagger}_{n} \sigma^{}_{z}a^{}_{n}+b^{\dagger}_{n}\sigma^{}_{z}b_{n} \right)+\sum^{}_{n}\left(a^{\dagger}_{n}t^{}_{\rm in }b^{}_{n}+a^{\dagger}_{n+1}t^{}_{\rm ex}b^{}_{n}+{\rm h.c.}\right)\ ,
\label{TBT}
\end{align}
where $n$ is the index of the unit cell,  and   $a^{}_{n}=\{a^{}_{n\Uparrow},a^{}_{n\Downarrow}\}$ and $b^{}_{n}=\{b^{}_{n\Uparrow},b^{}_{n\Downarrow}\}$  are spinors.




For a spinless DQD chain, the tunneling matrices in Eqs.~(\ref{T-IN}) and (\ref{T-EX}) become real numbers, $t^{}_{{\rm in},0}$ and $t^{}_{{\rm ex},0}$, and the effective tight-binding Hamiltonian  is identical to that of the single-level Su-Schrieffer-Heeger model~\cite{Su1979}
\begin{align}
H^{}_{{\rm T},0}=\sum^{}_{n}\left(t^{}_{{\rm in},0 } a^{\dagger}_{n}b^{}_{n}+t^{}_{{\rm ex},0 } a^{\dagger}_{n+1}b^{}_{n}+{\rm h.c.}\right),
\label{TB0}
\end{align}
with $a^{}_{n}$ and $b^{}_{n}$ corresponding to the operators  for the lowest localized orbitals on the QD1 and QD2, respectively.


\section{The Bloch spectrum of the double-quantum-dot chain ~\label{S-2}}



\begin{figure}
\centering
\includegraphics[width=0.4\textwidth]{Figs2.pdf}
\caption{  The spatial distribution of the confinement potential $U(x)$ in a unit cell; separate the unit cell into four regions (denoted by $I$, $II$, $III$ and $IV$)  in which the Bloch functions are derived (see text).}
\label{Figs2}
\end{figure}


Here the detailed calculation of the Bloch bands of our system is given.
For a one-dimensional lattice obeying periodic boundary conditions, the eigenstates are  the Bloch wave functions
\begin{align}
\Psi^{}_{k}(x)=\exp\left(ikx\right)u^{}_{k}(x),
\end{align}
 where $k$ is the wave vector and $u^{}_{k}(x)$ is the cell-periodic Bloch wave function. Obviously, the corresponding energy spectrum $\varepsilon(k)$ is given by
 \begin{align}
 H \Psi^{}_{k}(x)=\varepsilon(k)\Psi^{}_{k}(x),
 \label{B-F}
 \end{align}
where
 \begin{align}
  H =\frac{p^{2}}{2m_{e}}-V^{}_{\rm c}+U(x)+\alpha p \sigma_{y}+\frac{\Delta^{}_{\rm z}}{2}\sigma^{}_{n}.
\end{align}

 Similarly to the derivation in Sec.~\ref{S-1}, the solutions of Eq.~(\ref{B-F}) are obtained from the eigenfunctions of the Schr\"{o}dinger equations with different gate potentials
\begin{align}
\widetilde{H}^{}_{\zeta=0,1,2}\Psi^{}_{k,\zeta }(x)=E \Psi_{k,\zeta}(x).
\label{BSZ}
\end{align}
where $\Psi^{}_{k,\zeta}(x)=\exp\left(ikx\right)u^{}_{k,\zeta}(x)$  and $\varepsilon(k)=E-V^{}_{\rm c}$ (measured from the chemical potential $V_{\rm c}$).
 Then, analogously to the eigenstates in Eq.~(\ref{LQD}),  there exist four-fold degenerate solutions to Eq.~(\ref{BSZ}) for a fixed $k$, of the form
 \begin{align}
u^{}_{k,\zeta,\lambda}(x)=\exp\big[i(\gamma^{}_{\zeta,\lambda}-k)x\big]\begin{pmatrix}
w^{}_{\zeta,\lambda}\\
d^{}_{\zeta,\lambda}
\end{pmatrix}
\label{BF}
\end{align}
with the wave vector $\gamma^{}_{\zeta,\lambda}$ and the two spin components ($u^{}_{\zeta,\lambda}$ and $d^{}_{\zeta,\lambda}$) determined by Eqs.~(\ref{wk}) and (\ref{wkc}). It is convenient to divide the unit cell into four parts,  as shown in Fig.~\ref{Figs2}(a). The Bloch functions in each segment are expanded as
\begin{align}
u^{}_{k,I} (x)=&\sum^{4}_{\lambda=1}c^{}_{I,\lambda}u^{}_{k,0,\lambda} (x)~~~~ u^{}_{k,II} (x)=\sum^{4}_{\lambda=1}c^{}_{II,\lambda}u^{}_{k,1,\lambda}(x) \nonumber\\
u^{}_{k,III} (x)=&\sum^{4}_{\lambda=1}c^{}_{III,\lambda}u^{}_{k,0,\lambda}(x) ~~~~ u^{}_{k,IV} (x)=\sum^{4}_{\lambda=1}c^{}_{IV,\lambda}u^{}_{k,2,\lambda}(x),
\label{4u}
\end{align}
with $c^{}_{\Lambda,\lambda}$ $(\Lambda=I,~II,~III,~IV)$ representing the, yet to be determined, sixteen coefficients. The entire Bloch wave function is then
\begin{align}
u^{}_{k}(x)=&\Theta(x+L-d)[1-\Theta(x+d)]u^{}_{k,I}(x)+[1-\Theta(|x|-d)]u^{}_{k,II}(x)+[1-\Theta(x-L+d)]\nonumber\\
&\times  \Theta(x-d) u^{}_{k,III}(x) + \Theta(|x|-L+d)[1-\Theta(|x|-L)]u^{}_{k,IV}(x).
\label{bfc}
\end{align}
Combined with the periodicity of the cell Bloch functions, i.e., $u^{}_{k}(x_{2})=u^{}_{k}(x_{5})$ and $u^{\prime}_{k}(x_{2})=u^{\prime}_{k}(x_{5})$ with $x^{}_{5}\equiv x^{}_{2}+2L= L+d$, the energy bands are determined by the continuity of the wave function at the interfaces, $\lim_{\epsilon \rightarrow 0 }[u^{}_{k}(\xi+\epsilon)-u^{}_{k}(\xi+\epsilon)]=0$ and $\lim_{\epsilon \rightarrow 0 }[u^{\prime}_{k}(\xi+\epsilon)-u^{\prime}_{k}(\xi+\epsilon)]=0$,  with $\xi=x^{}_{1},x^{}_{3},x^{}_{4}$. Regarding the coefficients of the combinations in Eqs. (\ref{4u}) as a variable vector, the prerequisite conditions can be written in a matrix form
\begin{align}
\mathbf{W}(k,E) \cdot \hat{s}=0,
\label{MAT2}
\end{align}
where the vector is
$\hat{s}^{T}=\{c^{}_{I,1},c^{}_{I,2},c^{}_{I,3},c^{}_{I,4},c^{}_{II,1},c^{}_{II,2},c^{}_{II,3},c^{}_{II,4},c^{}_{III,1},c^{}_{III,2},c^{}_{III,3},c^{}_{III,4},c^{}_{IV,1},c^{}_{IV,2},c^{}_{IV,3},c^{}_{IV,4}\}$ and the matrix is
\begin{footnotesize}
\begin{align}
\mathbf{W}=\left(\begin{array}{cccccccccccccccc}
u^{(1,\uparrow)}_{k,0,\lambda}&u^{(1,\uparrow)}_{k,0,\lambda}&u^{(1,\uparrow)}_{k,0,\lambda}&u^{(1,\uparrow)}_{k,0,\lambda}&-u^{(1,\uparrow)}_{k,1,\lambda}&-u^{(1,\uparrow)}_{k,1,\lambda}&-u^{(1,\uparrow)}_{k,1,\lambda} &-u^{(1,\uparrow)}_{k,1,\lambda} &0&0&0&0 &0&0 &0&0\\
u^{(1,\downarrow)}_{k,0,\lambda}&u^{(1,\downarrow)}_{k,0,\lambda}&u^{(1,\downarrow)}_{k,0,\lambda}&u^{(1,\downarrow)}_{k,0,\lambda}&-u^{(1,\downarrow)}_{k,1,\lambda}&-u^{(1,\downarrow)}_{k,1,\lambda}&-u^{(1,\downarrow)}_{k,1,\lambda} &-u^{(1,\downarrow)}_{k,1,\lambda} &0&0&0&0 &0&0 &0&0\\
0&0&0&0&-u^{(3,\uparrow)}_{k,1,\lambda}&-u^{(3,\uparrow)}_{k,1,\lambda}&-u^{(3,\uparrow)}_{k,1,\lambda} &-u^{(3,\uparrow)}_{k,1,\lambda} &u^{(3,\uparrow)}_{k,0,\lambda}&u^{(3,\uparrow)}_{k,0,\lambda}&u^{(3,\uparrow)}_{k,0,\lambda}&u^{(3,\uparrow)}_{k,0,\lambda} &0&0 &0&0\\
0&0&0&0&-u^{(3,\downarrow)}_{k,1,\lambda}&-u^{(3,\downarrow)}_{k,1,\lambda}&-u^{(3,\downarrow)}_{k,1,\lambda} &-u^{(3,\downarrow)}_{k,1,\lambda} &u^{(3,\downarrow)}_{k,0,\lambda}&u^{(3,\downarrow)}_{k,0,\lambda}&u^{(3,\downarrow)}_{k,0,\lambda}&u^{(3,\downarrow)}_{k,0,\lambda} &0&0 &0&0\\
0&0&0&0&0&0&0 &0&u^{(4,\uparrow)}_{k,0,\lambda}&u^{(4,\uparrow)}_{k,0,\lambda}&u^{(4,\uparrow)}_{k,0,\lambda}&u^{(4,\uparrow)}_{k,0,\lambda} &-u^{(4,\uparrow)}_{k,2,\lambda}&-u^{(4,\uparrow)}_{k,2,\lambda}&-u^{(4,\uparrow)}_{k,2,\lambda} &-u^{(4,\uparrow)}_{k,2,\lambda}\\
0&0&0&0&0&0&0 &0&u^{(4,\downarrow)}_{k,0,\lambda}&u^{(4,\downarrow)}_{k,0,\lambda}&u^{(4,\downarrow)}_{k,0,\lambda}&u^{(4,\downarrow)}_{k,0,\lambda} &-u^{(4,\downarrow)}_{k,2,\lambda}&-u^{(4,\downarrow)}_{k,2,\lambda}&-u^{(4,\downarrow)}_{k,2,\lambda} &-u^{(4,\downarrow)}_{k,2,\lambda}\\
u^{(2,\uparrow)}_{k,0,\lambda}&u^{(2,\uparrow)}_{k,0,\lambda}&u^{(2,\uparrow)}_{k,0,\lambda}&u^{(2,\uparrow)}_{k,0,\lambda}&0&0&0 &0 &0&0&0&0 &-u^{(5,\uparrow)}_{k,2,\lambda}&-u^{(5,\uparrow)}_{k,2,\lambda}&-u^{(5,\uparrow)}_{k,2,\lambda} &-u^{(5,\uparrow)}_{k,2,\lambda}\\
u^{(2,\downarrow)}_{k,0,\lambda}&u^{(2,\downarrow)}_{k,0,\lambda}&u^{(2,\downarrow)}_{k,0,\lambda}&u^{(2,\downarrow)}_{k,0,\lambda}&0&0&0 &0 &0&0&0&0 &-u^{(5,\downarrow)}_{k,2,\lambda}&-u^{(5,\downarrow)}_{k,2,\lambda}&-u^{(5,\downarrow)}_{k,2,\lambda} &-u^{(5,\downarrow)}_{k,2,\lambda}\\
\tilde{u}^{(1,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(1,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(1,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(1,\uparrow)}_{k,0,\lambda}&-\tilde{u}^{(1,\uparrow)}_{k,1,\lambda}&-\tilde{u}^{(1,\uparrow)}_{k,1,\lambda}&-\tilde{u}^{(1,\uparrow)}_{k,1,\lambda} &-\tilde{u}^{(1,\uparrow)}_{k,1,\lambda} &0&0&0&0 &0&0 &0&0\\
\tilde{u}^{(1,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(1,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(1,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(1,\downarrow)}_{k,0,\lambda}&-\tilde{u}^{(1,\downarrow)}_{k,1,\lambda}&-\tilde{u}^{(1,\downarrow)}_{k,1,\lambda}&-\tilde{u}^{(1,\downarrow)}_{k,1,\lambda} &-\tilde{u}^{(1,\downarrow)}_{k,1,\lambda} &0&0&0&0 &0&0 &0&0\\
0&0&0&0&-\tilde{u}^{(3,\uparrow)}_{k,1,\lambda}&-\tilde{u}^{(3,\uparrow)}_{k,1,\lambda}&-\tilde{u}^{(3,\uparrow)}_{k,1,\lambda} &-\tilde{u}^{(3,\uparrow)}_{k,1,\lambda} &\tilde{u}^{(3,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(3,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(3,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(3,\uparrow)}_{k,0,\lambda} &0&0 &0&0\\
0&0&0&0&-\tilde{u}^{(3,\downarrow)}_{k,1,\lambda}&-\tilde{u}^{(3,\downarrow)}_{k,1,\lambda}&-\tilde{u}^{(3,\downarrow)}_{k,1,\lambda} &-\tilde{u}^{(3,\downarrow)}_{k,1,\lambda} &\tilde{u}^{(3,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(3,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(3,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(3,\downarrow)}_{k,0,\lambda} &0&0 &0&0\\
0&0&0&0&0&0&0 &0&\tilde{u}^{(4,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(4,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(4,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(4,\uparrow)}_{k,0,\lambda} &-\tilde{u}^{(4,\uparrow)}_{k,2,\lambda}&-\tilde{u}^{(4,\uparrow)}_{k,2,\lambda}&-\tilde{u}^{(4,\uparrow)}_{k,2,\lambda} &-\tilde{u}^{(4,\uparrow)}_{k,2,\lambda}\\
0&0&0&0&0&0&0 &0&\tilde{u}^{(4,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(4,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(4,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(4,\downarrow)}_{k,0,\lambda} &-\tilde{u}^{(4,\downarrow)}_{k,2,\lambda}&-\tilde{u}^{(4,\downarrow)}_{k,2,\lambda}&-\tilde{u}^{(4,\downarrow)}_{k,2,\lambda} &-\tilde{u}^{(4,\downarrow)}_{k,2,\lambda}\\
\tilde{u}^{(2,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(2,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(2,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(2,\uparrow)}_{k,0,\lambda}&0&0&0 &0 &0&0&0&0 &-\tilde{u}^{(5,\uparrow)}_{k,2,\lambda}&-\tilde{u}^{(5,\uparrow)}_{k,2,\lambda}&-\tilde{u}^{(5,\uparrow)}_{k,2,\lambda} &-\tilde{u}^{(5,\uparrow)}_{k,2,\lambda}\\
\tilde{u}^{(2,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(2,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(2,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(2,\downarrow)}_{k,0,\lambda}&0&0&0 &0 &0&0&0&0 &-\tilde{u}^{(5,\downarrow)}_{k,2,\lambda}&-\tilde{u}^{(5,\downarrow)}_{k,2,\lambda}&-\tilde{u}^{(5,\downarrow)}_{k,2,\lambda} &-\tilde{u}^{(5,\downarrow)}_{k,2,\lambda}
\end{array}\right)
\label{CMAT}
\end{align}
\end{footnotesize}
with $u^{(j,\uparrow)}_{k,\zeta,\lambda}= \exp[i(\gamma^{}_{\zeta,\lambda}-k)x_{j}]w^{}_{\zeta,\lambda}$ ($j=1,2,3,4,5$), $u^{(j,\downarrow)}_{k,\zeta,\lambda}= \exp[i(\gamma^{}_{\zeta,\lambda}-k)x_{j}]d^{}_{\zeta,\lambda}$, $\tilde{u}^{(j,\uparrow)}_{k,\zeta,\lambda}= (\gamma^{}_{\zeta,\lambda}-k)\exp[i(\gamma^{}_{\zeta,\lambda}-k)x_{j}]w^{}_{\zeta,\lambda}$, and $\tilde{u}^{(j,\downarrow)}_{k,\zeta,\lambda}= (\gamma^{}_{\zeta,\lambda}-k)\exp[i(\gamma^{}_{\zeta,\lambda}-k)x_{j}]d^{}_{\zeta,\lambda}$. A nontrivial solution to Eq.~(\ref{MAT2}) requires
\begin{align}
\det [\mathbf{W} (k,E)]=0,
\label{GF}
\end{align}
 and yields an implicit equation for $k$ and $E$ of Eqs.~(\ref{wk}),~(\ref{wkc}), and ~(\ref{CMAT}). The exact energy spectrum of the Bloch band $\varepsilon (k)=E-V^{}_{\rm c}$ is found by solving numerically for $E$ at a fixed value of $k$. Substituting the so-derived values of $E$ in Eq.~(\ref{MAT2}) yields the values of the coefficients $c^{}_{\Lambda,\lambda}$, i.e., gives the specific form of the Bloch wave functions, with the normalization condition $\int^{L}_{-L}u^{\ast}_{k}(x)u^{ }_{k}(x)dx=1$.


In the absence of a magnetic field or when $\theta=0$, the two spin degrees of freedom are separable and the number of the corresponding coefficients  for ascertaining the Bloch band is reduced by half, which is similar to a spinless chain.  In this case, the matrix $\mathbf{W}(k,E)$ is  decomposed into two independent submatrices whose determinants yield the Bloch spectrum, as shown in Fig.3 (a) for the case of no magnetic field.











\section{Edge states of a finite double-quantum-dot chain}
\label{S-3}



\begin{figure}
\centering
\includegraphics[width=0.68\textwidth]{Figs3.pdf}
\caption{  The energy  spectrum (as a function of $k$) of a double quantum-dot chain made of  InAs, for different values of the magnetic field and in the absence of the spin-orbit coupling. Here $\Delta V=10$~meV and $\theta=\pi/2$. (a)  Small magnetic fields,    $\Delta^{}_{\rm z}<t^{}_{\rm ex,0 }-t^{}_{\rm in,0 }$. (b) Large  magnetic fields,  $\Delta^{}_{\rm z}>t^{}_{\rm ex,0 }+t^{}_{\rm in,0 }$. $k$ is measured in units of $1/(2L)$ and the energies of two pure spin edge states, $\varepsilon^{}_{\uparrow}$ and $\varepsilon^{}_{\downarrow}$, are also indicated.}
\label{Figs3}
\end{figure}



Here we consider the edge modes of an open  chain of DQDs.
 We begin with the configuration where there is no Zeeman field,
and the electron spin is polarized solely by the field  induced by the spin-orbit interaction.
In this case the tight-binding Hamiltonian reads
 \begin{align}
H^{}_{\rm T,S}=\sum^{}_{n} (t^{}_{\rm in,0}a^{\dagger}_{n}e^{i\varphi^{}_{\rm in}\sigma^{}_{z}}b^{}_{n}+t^{}_{\rm ex,0}a^{\dagger}_{n+1}e^{i\varphi^{}_{\rm ex}\sigma^{}_{z}}b^{}_{n}+{\rm h.c.} ),
\nonumber
 \end{align}
and there are no spin-flip transitions.

The energy eigenfunctions in an  open DQD chain are distinguished by the different quasi-spin states
\begin{align}
|\Psi^{}_{\tau=\Uparrow,~\Downarrow }\rangle=\sum^{N}_{n=1}\left( A^{}_{n,\tau} a^{\dagger}_{n\tau}+B^{}_{n,\tau}b^{\dagger}_{n\tau}\right)|0\rangle,
\label{F-LS}
\end{align}
where $|0\rangle$ is the vacuum state.  The coefficients
$A^{}_{n,\tau}$ and $B^{}_{n,\tau}$
 are determined by the spin-dependent Schr\"{o}dinger equation, $H^{}_{\rm T,S}|\Psi^{}_{\tau}\rangle=E|\Psi^{}_{\tau}\rangle$, where $E^{}_{}$ is  the energy eigenvalue. Specifically, the Schr\"{o}dinger equation for the quasi-spin up state can be written as
\begin{align}
\sum^{}_{n}\left( t^{}_{\rm in, 0}e^{i\varphi^{}_{\rm in}}a^{\dagger}_{n\Uparrow} b^{}_{n\Uparrow}+ t^{}_{\rm ex,0}e^{i\varphi^{}_{\rm ex}}a^{\dagger}_{n+1\Uparrow}b^{}_{n\Uparrow}+{\rm h.c.}\right)|\Psi^{}_{ \Uparrow }\rangle=E^{}_{ }|\Psi^{}_{ \Uparrow}\rangle.
\label{EGS}
\end{align}
Then, the explicit Shr\"{o}dinger equations in  the $n$-th cell are
\begin{align}
 t^{}_{\rm in, 0}e^{i\varphi^{}_{\rm in}}B^{}_{n,\Uparrow}+  t^{}_{\rm ex, 0}e^{i\varphi^{}_{\rm ex}}B^{}_{n-1,\Uparrow}=E^{}_{} A^{}_{n,\Uparrow}\ ,~~~~~~~~~ t^{}_{\rm in, 0}e^{-i\varphi^{}_{\rm in}}A^{}_{n,\Uparrow}+ t^{}_{\rm ex, 0}e^{-i\varphi^{}_{\rm in}} A^{}_{n+1,\Uparrow} =E^{}_{}B^{}_{n,\Uparrow}\ .
 \label{CQES}
 \end{align}
These can be reduced to a recursion formula when $E^{}_{ } \neq 0$
\begin{align}
X^{}_{n,\Uparrow}\left(t^{2}_{\rm in,0}+t^{2}_{\rm ex, 0}-E_{}^{2}\right)+ t^{}_{\rm in,0}t^{}_{\rm ex,0}e^{i(\varphi^{}_{\rm in}-\varphi^{}_{\rm ex})}X^{}_{n+1,\Uparrow}+t^{}_{\rm in,0}t^{}_{\rm ex,0}e^{i(\varphi^{}_{\rm ex}-\varphi^{}_{\rm in})}X^{}_{n-1,\Uparrow}=0\ ,
\end{align}
with $X=A,~B$. This recursion relation is solved by  the transfer-matrix scheme
\begin{align}
\begin{pmatrix}
X_{n+1,\Uparrow}\\
X_{n,\Uparrow}
\end{pmatrix}=\mathbf{T}\begin{pmatrix}
X_{n,\Uparrow}\\
X_{n-1,\Uparrow}
\end{pmatrix}\ ,
\end{align}
with
\begin{align}
\mathbf{T}=\begin{pmatrix}
 m^{}_{11}& m^{}_{12}\\
1&0
\end{pmatrix}\ ,
\end{align}
in which $m^{}_{11} =e^{i(\varphi^{}_{\rm ex}-\varphi^{}_{\rm in})}(E^{2}_{}-t^{2}_{\rm in,0}-t^{2}_{\rm ex, 0})/(t^{}_{\rm in,0}t^{}_{\rm ex,0})$ and $m^{}_{12}=e^{2i(\varphi^{}_{\rm ex}-\varphi^{}_{\rm in})}$. The formation of edge modes hinges on the eigenvalues $\lambda^{}_{\pm}= \frac{1}{2} \left(m^{}_{11}\pm\sqrt{m^{2}_{11}+4m^{}_{12}}\right)$ of the transfer matrix $\mathbf{T}$. An eigenvalue with an absolute value larger than 1, $|\lambda^{}_{+}|>1$, excludes the possibility to form edge modes with nonzero energy~\cite{Li2015}.
Thus, the edge modes only occur for $E^{}_{ }=0$. The recursion relations for this case are
\begin{align}
t^{}_{\rm in, 0}e^{-i\varphi^{}_{\rm in }} A_{n,\Uparrow}+t^{}_{\rm ex,0} e^{-i\varphi^{}_{\rm ex }} A_{n+1,\Uparrow}=0\ ,~~~~~~~~  t_{\rm in, 0} e^{i\varphi^{}_{\rm in }}B_{n,\Uparrow}+t_{\rm ex, 0}e^{i\varphi^{}_{\rm ex }} B_{n-1,\Uparrow}=0 .
\end{align}
As occurs for the Su-Schrieffer-Heeger model~\cite{Su1979}, the zero-energy edge modes appear only when $ t^{}_{\rm in,0}/t^{}_{\rm ex, 0}<1$. By introducing the notation  $ r^{}_{ }=(t^{}_{\rm in,0}/t^{}_{\rm ex, 0})e^{i(\varphi^{}_{\rm in }-\varphi^{}_{\rm ex })}$, the zero-energy edge modes are given by
\begin{align}
|E^{}_{a,\Uparrow}\rangle=A^{}_{1,\Uparrow}\sum^{N}_{n=1}(r_{}^{\ast})^{n-1}a^{\dagger}_{n\Uparrow}|0\rangle\ ,  ~~~~~~~~ |E^{}_{b,\Uparrow}\rangle=B^{}_{N,\Uparrow}\sum^{N}_{n=1}r_{}^{ N-n}b^{\dagger}_{n\Uparrow}|0\rangle,
\label{SOE}
\end{align}
with the coefficients $A^{}_{1,\Uparrow}$ and $B^{}_{N,\Uparrow}$  determined by the normalization conditions  $\langle E^{}_{a/b,\Uparrow} | E^{}_{a/b ,\Uparrow}\rangle=1$. Meanwhile, in the case of a small $N$, the overlap between $|E^{}_{a,\Uparrow}\rangle$ and $|E^{}_{b,\Uparrow}\rangle$    results in the hybridization between the two zero-energy edge modes, as indicated in Fig.~3(b) for an open spinless chain with $N=20$.
Analogously, the quasi-spin down edge states can be derived by the same method as above, and all the zero-energy edge modes are protected by  the chiral symmetry as explained in the main text.



 The situation is changed when a Zeeman interaction is applied. We first consider the Hamiltonian in the absence of spin-orbit coupling, i.e., for
 $\varphi^{}_{\rm in/ex}=0$,
\begin{align}
 H^{ }_{  \rm  T,Z}= \sum^{}_{n}\frac{\Delta^{}_{\rm z}}{2}\left(a^{\dagger}_{n}\sigma^{}_{z}a^{}_{n}+b^{\dagger}_{n}\sigma^{}_{z}b^{}_{n}\right)+\sum^{}_{n}\sum^{}_{\tau=\uparrow,\downarrow}\left(t^{}_{\rm in,0} a^{\dagger}_{n\tau}b^{}_{n\tau}+t^{}_{\rm ex, 0}a^{\dagger}_{n+1\tau}b^{}_{n\tau}+ \rm{h.c.}\right).
 \label{TBP}
  \end{align}
  Notably, because there is no spin-orbit interaction, in this case $\tau=\uparrow,\downarrow$ represent the two pure spin states  polarized by the Zeeman field.
 In Fourier space, this   Hamiltonian reads
\begin{align}
H^{}_{\rm T,Z}(k)= \begin{pmatrix}
\Delta^{}_{\rm z}/2 &t^{}_{\rm  in, 0 }+e^{-ik}t^{}_{\rm  ex,0 }&0&0\\
t^{ }_{\rm  in, 0}+e^{ik}t^{ }_{\rm ex,0} &\Delta^{}_{\rm z}/2&0&0\\
0&0&-\Delta^{}_{\rm z}/2& t^{}_{\rm  in,0 }+e^{-ik}t^{}_{\rm  ex,0}\\
0&0&t^{}_{\rm  in,0 }+e^{ik}t^{\ast }_{\rm  ex,0} &-\Delta^{}_{\rm z}/2
\end{pmatrix},
\end{align}
and the energies of the spin-dependent Bloch bands are
\begin{align}
E^{ }_{\uparrow \pm} (k)=& \left(\Delta^{}_{\rm z}/2\right)\pm \sqrt{t^{2}_{\rm in,0}+ t^{2}_{\rm ex,0} +2t^{ }_{\rm in,0}t^{ }_{\rm ex,0}\cos k}\nonumber\\
E^{ }_{\downarrow \pm} (k)=& -\left(\Delta^{}_{\rm z}/2\right) \pm \sqrt{t^{2}_{\rm in,0}+ t^{2}_{\rm ex,0} +2t^{ }_{\rm in,0}t^{ }_{\rm ex,0}\cos k}.
\label{enes}
 \end{align}
 As seen, the Hamiltonian Eq.~(\ref{TBP}) separates  into two independent   models, each being described by   a Su-Schrieffer-Heeger model with a finite on-site energy. Consequently, one expects two pairs of  edge modes in the topological regime in which $t^{}_{\rm ex,0}>t^{}_{\rm in,0}$. However, when
\begin{align}
t^{}_{\rm ex,0}- t^{}_{\rm in,0}\leq\Delta^{} _{\rm z} \leq t^{}_{\rm ex,0}+ t^{}_{\rm in,0},
\label{BBC}
\end{align}
the overlap of an edge mode belonging to
a certain spin projection with the bulk of states of the opposite projection
masks its topological character.
It follows that in the absence of spin-orbit coupling the edge states reside within the energy gap provided that
 $ \Delta^{}_{\rm z}<t^{}_{\rm ex,0}-t^{}_{\rm in,0} $ at small magnetic fields, and  $ \Delta^{}_{\rm z}>t^{}_{\rm ex,0}+t^{}_{\rm in,0} $ at large  fields, as shown in Fig.~ \ref{Figs3}.


When both the spin-orbit coupling and the Zeeman magnetic field are active, the energy eigenfunctions of an open dimerized chain are  the quasi-spin mixed states
\begin{align}
|\Psi^{}_{}\rangle =\sum^{N}_{n=1} \left(A^{}_{n,\Uparrow}a^{\dagger}_{n\Uparrow}+ A^{}_{n,\Downarrow}a^{\dagger}_{n\Downarrow}+B^{}_{n,\Uparrow}b^{\dagger}_{n\Uparrow}+B^{}_{n,\Downarrow}b^{\dagger}_{n\Downarrow}\right)|0\rangle.
\label{PSF0}
\end{align}
The coefficients are determined
 by the Schr\"{o}dinger equation
$H^{}_{\rm T}|\Psi^{}_{}\rangle =E_{ }|\Psi^{}_{}\rangle$, with the discrete Hamiltonian $H^{}_{\rm T}$ given in Eq.~(\ref{TBT}).
The edge modes are distinguished by the spatial density distributions of the state, i.e., they are residing on either side of the chain. We quantify this property by choosing
 \begin{align}
 \eta^{}_{\rm L} =\sum^{3}_{n=1}P^{}_{n}+\sum^{N}_{n=N-2}P^{}_{n}
\end{align}
with $P^{}_{n}= |A^{}_{n,\Uparrow}|^{2}+|B^{}_{n,\Uparrow}|^{2}+ |A^{}_{n,\Downarrow}|^{2}+|B^{}_{n,\Downarrow}|^{2}$, $n=1,2,3$ indicating the leftmost three sites of the chain,  and $n=N-2,N-1,N$ corresponding to rightmost three sides.  For an open InAs chain of $N=20$ sites and $\Delta V=10$~meV, the specifical values of $\eta^{}_{\rm L}$ for the edge states under different magnetic fields are scaled by the color bar shown in Fig.~\ref{Figs4}. In order to count the number of  edge modes as more as possible, here  we assume the presence of an edge state requires $\eta^{}_{\rm L}\geq 0.45$, with the  threshold value determined by the probability density distribution of a predestined edge state in  the prototype SSH model with a small degree of dimerization, $\Delta V=1.0$~meV.


Being consistent with the results shown in Fig.~3(d), the numerical values of  $\eta^{}_{\rm L}$  in Fig~\ref{Figs4}(a) confirm  the appearance of  edge states  near the zero degrees of dimerization ( i.e., $\Delta V\rightarrow 0^{+}$), but characterized with the weak localization.   Under the interplay between the spin-orbit and Zeeman interactions, the large locality  parameter $\eta^{}_{\rm L}$  in Fig.~\ref{Figs4}(b) demonstrates the presence of additional edge modes living in the gap induced by the anti-crossing of the energy levels. However, in contrast to the Zeeman-split edge states, these  extra edge modes are not stable and can be merged into the bulk states by decreasing the tilting angle $\theta$ for $0\leq\theta<\pi/2$, and that is  manifested by the reduction of $\eta^{}_{\rm L}$ shown in Fig.~\ref{Figs4}(c). Besides,  the insensitivity of the energy values to $\Delta^{}_{\rm z}$, seen in Fig~\ref{Figs4}(b), can also be illustrated by  the high degree of spin mixing intrinsic to the additional edge states, as reflected by the small spin  density distribution, i.e., $S_{n}=\left(|A^{}_{n\Uparrow}|^{2}+|B^{}_{n\Uparrow}|^{2}-|A^{}_{n\Downarrow}|^{2}-|B^{}_{n\Downarrow}|^{2}\right)/2$, shown in Fig.~\ref{Figs4}(d).

\begin{figure}
\centering
\includegraphics[width=0.7\textwidth]{Figs4.pdf}
\caption{ (a) The energy spectrum of an open InAs  chain comprising 20 sites as a function of the potential difference $\Delta V$, with the magnetic field $B=0.2$~T, the tilting angle $\theta=0.5\pi$, and the spin-orbit length $x^{}_{\rm so}=180$~nm. Moreover, for each edge state, the locality parameter $\eta^{}_{\rm L}$ is scaled by the color bar shown in the right side. (b) The energy spectrum of the open InAs chain as a function of the Zeeman splitting $\Delta_{\rm z}$, with $\Delta V=10$~nm, $\theta=0.5\pi$, and  $x^{}_{\rm so}=180$~nm. (c) For the case of $B=1.0$~T, the energy spectrum of the open chain as a function of $\theta$, when $\Delta V=10$~nm and $x^{}_{\rm so}=180$~nm.  Moreover, the spatial spin density distributions  of the two  edge modes, as indicated by a square and a diamond in panel (c), are displayed in (d).  }
\label{Figs4}
\end{figure}


 \footnotetext[1]{asb}





