\documentclass[prl,twocolumn,showpacs,preprintnumbers,amsmath,amssymb,superscriptaddress]{revtex4-2}

\usepackage{graphicx}
\usepackage{dcolumn}
\usepackage{bm}
\usepackage{bbm}
\usepackage{color}
\usepackage{bbm}
\usepackage{amsfonts}
\usepackage{soul}
%\usepackage{ulem}
\usepackage[colorlinks=true,citecolor=blue,linkcolor=blue,urlcolor=blue]{hyperref}
%\renewcommand\refname{}

\newcommand{\zlx}[1]{\textcolor{blue}{#1}}
\newcommand{\xll}[1]{\textcolor{red}{#1}}

\setcounter{secnumdepth}{3}

\hyphenpenalty=2000
\tolerance=800


\begin{document}

\cleardoublepage
\onecolumngrid
\clearpage
\setcounter{equation}{0}
\setcounter{figure}{0}
\setcounter{table}{0}
\makeatletter
\renewcommand{\thesection}{\Roman{section}}
\renewcommand{\theequation}{S\arabic{equation}}
\renewcommand{\thefigure}{S\arabic{figure}}
\renewcommand{\thetable}{S\Roman{table}}
\renewcommand{\bibnumfmt}[1]{[S#1]}
\renewcommand{\citenumfont}[1]{S#1}





\begin{center}
        \textbf{\large Supplementary material for: \\ Controllable Operations of Edge States in Cross-One-dimensional Topological Chains}
\end{center}



\section{Topological edge states in a cross pair of SSH chains}
The model we consider in the main text is described by two identical SSH chains and a central four-node junction, where four edge states emerge when the the system enters topological phase. The Hamiltonian of SSH model is described by
\begin{equation}
H_\mathrm{SSH}=\sum_{i=1}^{N}J_1(\sigma_{i,A}^{+}\sigma_{i,B}^{-}+\mathrm{H.c.})+\sum_{i=1}^{N-1}J_2(\sigma_{i+1,A}^{+}\sigma_{i,B}^{-}+\mathrm{H.c.}),
\end{equation}
where $\sigma_{i,A}^{-}/\sigma_{i,B}^{-}$ ($\sigma_{i,A}^{+}/\sigma_{i,B}^{+}$) are the lowering (raising) operators of atoms in sublattice $A/B$ at $i$-th unit cell; $J_1$ and $J_2$ are the inter- and intra- cell hopping strength. We can write the bulk Hamiltonian with periodic boundary condition (PBC) in momentum space as $H_\mathrm{SSH}=\sum_k S_k^\dagger \mathcal{H}(k) S_k$, $S_k=(\sigma_{A,k}^{-},\sigma_{B,k}^{-})^T$, with the single-particle bulk Hamiltonian given by
\begin{equation}
\mathcal{H}(k)=
\left(\begin{array}{cc}
0 & \rho(k) \\
\rho^*(k) & 0
\end{array}\right)=\mathrm{Re}(\rho)\sigma_x+\mathrm{Im}(\rho)\sigma_y,\quad  \rho(k)=J_1+J_2 e^{-i k}=|\rho(k)|e^{i\phi(k)},
\end{equation}
where $\sigma_{A,k}^{-}=\sum_n e^{-ikn}\sigma_{n,A}/\sqrt{N}$, $\sigma_{B,k}^{-}=\sum_n e^{-ikn}\sigma_{n,B}/\sqrt{N}$ are the lowering operators in momentum space after Fourier transformation, and the phase is given by $\phi(k)=\mathrm{arccot}\left(\frac{J_1}{J_2 \sin{k}}+\cot{k}\right)$.
\begin{figure}[ht]
\centering
\includegraphics[width=0.9\linewidth]{FigS1}
\caption{SSH model. Energy spectrum with (a) periodic boundary condition (b) open boundary condition. (c). The winding number $\nu=1$ in topological non-trivial phase ($J_1<J_2$).}
\label{FigS1}
\end{figure}
By diagonalizing $\mathcal{H}(k)$ we obtain the bulk energy band $\epsilon_{\pm}(k)=\pm |\rho(k)|=\sqrt{J_1^2+J_2^2+2J_1 J_2 \cos{k}}$, as shown in Fig.~\ref{FigS1}(a), and eigenfunctions
\begin{equation}
|\psi_\pm(k)\rangle=\frac{1}{\sqrt{2}}
\left(\begin{array}{c}
\pm e^{-i \phi(k)} \\ 1
\end{array}\right).
\end{equation} 
As shown in Figs. \ref{FigS1}(b) and \ref{FigS1}(c), when $J_1<J_2$, the system enters topological phase with non-trivial winding number
\begin{equation}
\nu=\frac{1}{2 \pi i} \int_{-\pi}^\pi d k \frac{d}{d k} \ln \rho(k)=\frac{1}{2\pi} \oint d k \frac{d \phi}{d k}=1.
\end{equation}
According to the bulk-boundary correspondence, the zero-energy edge states appear in the bandgap and localize at the left and right end of the chain for open boundary condition (OBC). Defining the ﬁrst-quantized Hamiltonian $\mathcal{H}_\mathrm{SSH}$ represented by a $2N\times 2N$ matrix in real space with $H_\mathrm{SSH}=\Psi^\dagger \mathcal{H}_\mathrm{SSH}\Psi$ and $\Psi=(\sigma_{1,A}^{-},\sigma_{1,B}^{-},\cdots,\sigma_{N,B}^{-})^{T}$, we introduce the chiral symmetry and inversion symmetry of the system both in the OBC
\begin{equation}\label{symmetryinOBC}
\begin{gathered}
\Gamma^{-1}\mathcal{H}_\mathrm{SSH} \Gamma=-\mathcal{H}_\mathrm{SSH}, \quad \Gamma =\sum_n|n,A\rangle\langle n,A|-|n,B\rangle\langle n,B|=\bigoplus_{n=1}^N \sigma_z, \\
I^{-1}\mathcal{H}_\mathrm{SSH}I=\mathcal{H}_\mathrm{SSH},\quad I=\sum_n|N+1-n,B\rangle\langle n,A|+|N+1-n,A\rangle\langle n,B|,
\end{gathered}
\end{equation}
and PBC
\begin{equation}
\begin{gathered}
\mathcal{C}^{-1}\mathcal{H}(k)\mathcal{C}=-\mathcal{H}(k),\quad \mathcal{C}=\sigma_z,\\
\mathcal{I}^{-1}\mathcal{H}(k)\mathcal{I}=\mathcal{H}(-k), \quad \mathcal{I}=\sigma_x.
\end{gathered}
\end{equation}
For eigenstate with energy $E_i$, the chiral symmetry predicts another eigenstate with energy $-E_i$:
\begin{equation}
H \Gamma |\psi_i\rangle=-\Gamma H|\psi_i\rangle=- E_i \Gamma|\psi_i\rangle.
\end{equation}
Thus, the spectrum of bulk and edge states are symmetric, and the wave functions solved with two different boundary conditions have odd or even parity. Additionally, the chiral symmetry operator follows anti-commutation relations $\{I,\Gamma\}=\{\mathcal{I},\mathcal{C}\}=0$. For an eigenstate with certain parity $I|\psi_i\rangle=p_i |\psi_i\rangle$, $p_i=\pm 1$, we have
\begin{equation}
I\Gamma|\psi_i\rangle=-\Gamma I|\psi_i\rangle=-p_i \Gamma\psi_i\rangle.
\end{equation}
Therefore, $\Gamma|\psi_i\rangle$ has opposite parity and eigenenergy compared with $|\psi_i\rangle$. The conclusions above are the same when eigenstates are solved in the PBC.

In the absence of translation symmetry, the initial Hamiltonian $H_\mathrm{tot}$ we introduced in the main text can be effectively described by applying the perturbation theory, where the total Hamiltonian can be decomposed as 
\begin{equation}\label{Htotpert}
H_\mathrm{tot}\equiv H^0+H'=H_\mathrm{SSH,1}+H_\mathrm{SSH,2}+H_I.
\end{equation}
Since we only concern with the dynamics between edge states, the \emph{model space} can be restricted to the four edge states while the complementary space contains all bulk states. Assume that the inter-chain coupling $K_i$ is sufficiently weak that the energy perturbation coefficient $K_i/\Delta_\mathrm{gap}\ll 1$, where $\Delta_\mathrm{gap}=|J_1-J_2|$ denotes the large gap between edge and bulk states. To partition the total Hamiltonian, we should construct diagonal unperturbed Hamiltonian $H^0$ first. We use $|v_{k_j},\pm\rangle$ to denote the wave functions solved with OBC, which can be reconstructed by the set of bulk-states basis in PBC~\cite{PhysRevB.84.195452,PhysRevResearch.2.012076},
\begin{equation}\label{bulkstates}
|v_{k_j},\pm\rangle = \frac{1}{\sqrt{N}}\sum_n \sin{[nk_j-\phi(k_j)]}|n,A\rangle \pm \sin{(n k_j)}|n,B\rangle,
\end{equation}
where $k_j$ fulfills the quantization condition $k_j(N+1)-\phi(k_j)=j \pi$, with $j=1,\cdots,N-1$.
As shown in Fig.~\ref{FigS1}(b), each $k_j$ corresponds to eigenstates $|\psi_j\rangle\equiv|v_{k_j},-\rangle$ and $\Gamma|\psi_j\rangle=|v_{k_j},+\rangle=|\psi_{2N+1-j}\rangle$ from lower band and upper band respectively, where $|\psi_j\rangle$ is in energy ascending order. On the other hand, due to the finite size of the topological chain, the edge states hybridize into symmetric and anti-symmetric states,
\begin{equation}
|\psi_{S/A}\rangle=\frac{1}{\sqrt{\mathcal{N}_h}}\sum_{n=1}^{N}(-1)^{n-1}\sinh{(N+1-n)\lambda}|n,A\rangle \pm (-1)^{N-n} \sinh{(n\lambda)}|n,B\rangle,
\end{equation}
where $\mathcal{N}_h$ is the normalization constant of the \emph{hybridized edge states}. The wave functions vanish at the boundary, which leads to $\sinh{(N \lambda)}/\sinh{\left[(N+1) \lambda\right]}=J_1/J_2$, and with chiral symmetry their energy level splits into
\begin{equation}\label{E_SE_A}
\epsilon_S=-\epsilon_A=\langle\psi_R|H_\mathrm{SSH}|\psi_L\rangle=(-1)^{N+1} \frac{J_2\sinh \lambda}{\sinh(N+1) \lambda}.
\end{equation}
The SSH Hamiltonian $H_\mathrm{SSH}$ can be expressed with two unhybridized edge states and $2N-2$ bulk states solved with OBC, i.e., $H_\mathrm{SSH}=H_\mathrm{SSH}^0+V_e$. In the single-particle representation, the SSH Hamiltonian is explicitly expressed as
\begin{equation}\label{HSSHpert}
\mathcal{H}_\mathrm{SSH}^0=\sum_j \epsilon_{\pm}(k_j)|v_{k_j}\rangle \langle v_{k_j}|,\quad \mathcal{V}_e=\epsilon_S |\psi_L\rangle\langle \psi_R| + \mathrm{H.c.},
\end{equation}
where the unhybridized edge states are the eigenstates of $H_\mathrm{SSH}^0$ and given by
\begin{equation}
|\psi_L\rangle=\frac{|\psi_S\rangle+|\psi_A\rangle}{\sqrt{2}}, \quad |\psi_R\rangle=\frac{|\psi_S\rangle-|\psi_A\rangle}{\sqrt{2}},
\end{equation}
with eigenenergies $\epsilon_L=\epsilon_R=0$. When the system is far from the phase transition critical point ($J_1/J_2=1$), the edge states are highly localized and thus we have $N\lambda\gg1$, $e^{-\lambda}\simeq J_1/J_2$. The unhybridized edge states can be replaced by approximate solutions~\cite{asboth2016short},
\begin{equation}\label{apsolutions}
|\psi_L\rangle=\sum_n \frac{1}{\sqrt{\mathcal{N}}}\left(-\frac{J_1}{J_2}\right)^{n}|n,A\rangle,\quad |\psi_R\rangle=\sum_n \frac{1}{\sqrt\mathcal{N}}\left(-\frac{J_1}{J_2}\right)^{N+1-n}|n,B\rangle,
\end{equation}
where $\mathcal{N}$ is the normalization constant of \emph{unhybridized edge states}.

With Eqs. \eqref{Htotpert} and \eqref{HSSHpert} we denote the zeroth-order Hamiltonian $H^0$ and perturbation $H'$ as
\begin{equation}
H^0=H_\mathrm{SSH,1}^0+H_\mathrm{SSH,2}^0,\quad H'=V_e^1+V_e^2+H_I,
\end{equation}
where the junction interaction $H_I$ in the main text reads
\begin{equation}
H_I=\left[(K_1\sigma_{1,s}^{+}+K_2\sigma_{1,s+1}^{+})\sigma_{2,s}^{-}+(K_3\sigma_{1,s}^{+}+K_4\sigma_{1,s+1}^{+})\sigma_{2,s+1}^{-}+\mathrm{H.c.}\right],
\end{equation}
and the index $s$ represents the $N$-th atom of the chain. For an even or odd number of cells, it's given by
\begin{equation}
\begin{cases}
{s=\{\lfloor(N+1)/2\rfloor,B\},~s+1=\{\lfloor N/2\rfloor+1,A\}, \quad \text{ for } N\in \text{even},}\\ 
{s=\{\lfloor N/2\rfloor+1,A\},~s+1=\{\lfloor(N+1)/2\rfloor,B\}, \quad \text{ for } N\in \text{odd},}
\end{cases}
\end{equation}
where $\lfloor \cdot \rfloor$ means floor function. The model space $\mathcal{M}$ is spanned by $\{|\psi_{i, L/R}\}$ with projector $\mathbb{P}=\sum_{i=1}^2 |\psi_{iL}\rangle\langle \psi_{iL}|+|\psi_{iR}\rangle\langle \psi_{iR}|$, whereas the complementary space is spanned by bulk states $\{|v_{i,k_j},\pm\rangle\}$ with projector $\mathbb{Q}=1-\mathbb{P}$. 

we use the bulk states and unhybridized edge states of $H^0$ to expand the single-excitation state $|i,A_n\rangle$ and $|i,B_n\rangle$. we denote $\xi$ as the orthogonal matrix that transform the real-space basis $\{\sigma_{i,A}^+,\sigma_{i,B}^+\}$ to eigenstates basis $\{\Psi_j^\dagger\}$ with energy ascending order,
\begin{equation}
(\sigma_{1,A}^+,\sigma_{1,B}^+,\cdots,\sigma_{N,B}^+)=(\Psi_1^\dagger,\cdots,\Psi_L^\dagger,\Psi_R^\dagger,\cdots,\Psi_{2N}^\dagger)\left(\begin{array}{cccc}
\xi_{1,1} & \xi_{1,2} & \cdots & \xi_{1, 2N} \\
\xi_{2,1} & \xi_{2,2} & \cdots & \xi_{2, 2N} \\
\vdots & \vdots & \ddots & \vdots \\
\xi_{2N, 1} & \xi_{2N, 2} & \cdots & \xi_{2N, 2N}
\end{array}\right),
\end{equation}
and thus we have
\begin{equation}
\sigma_s^{+}|G\rangle=\sum_j\xi_{j,N}\Psi_j^\dagger|G\rangle=\sum_j\xi_{j,N}|\psi_j\rangle, \quad \sigma_{s+1}^{+}|G\rangle=\sum_j\xi_{j,N+1}\Psi_j^\dagger|G\rangle=\sum_j\xi_{j,N+1}|\psi_j\rangle,
\end{equation}
where $|G\rangle=|g_1 g_2\cdots g_{2N}\rangle$ is the ground state. The inverse transformation $\xi^{-1}=\xi^{T}$ gives $|\psi_j\rangle=\sum_i \xi_{j,2i-1}|A_i\rangle+\xi_{j,2i}|B_i\rangle$. Therefore, the matrix element $\xi_{j,N}$ and $\xi_{j,N+1}$ are determined by the expression of $|\psi_j\rangle$.

The junction interaction can be decomposed into three parts:
\begin{equation}
H_I=H_{ee}+H_{eb}+H_{bb},
\end{equation}
corresponding to the edge-edge coupling, edge-bulk coupling, and bulk-bulk coupling, respectively. By projecting onto model space and complementary space, one can write the secular equation as $(H^0+H')(\mathbb{P}+\mathbb{Q})|\psi\rangle = E|\psi\rangle$, which gives the following two equations:
\begin{align}
\label{PHP}(\mathbb{P}H^0 \mathbb{P}+\mathbb{P}H' \mathbb{P})(\mathbb{P}|\psi\rangle)+\mathbb{P}H' \mathbb{Q}(\mathbb{Q}|\psi\rangle)&=E \mathbb{P}|\psi\rangle, \\
\label{QHQ}(\mathbb{Q}H^0 \mathbb{Q}+\mathbb{Q}H' \mathbb{Q})(\mathbb{Q}|\psi\rangle)+\mathbb{Q}H' \mathbb{P}(\mathbb{P}|\psi\rangle)&=E \mathbb{Q}|\psi\rangle.
\end{align}
We solve $\mathbb{Q}|\psi\rangle$ from Eq. \eqref{QHQ} and substitute it into Eq. \eqref{PHP}, which leads to
\begin{equation}
H_\mathrm{eff} \mathbb{P}|\psi\rangle=E \mathbb{P}|\psi\rangle,
\end{equation}
where the effective Hamiltonian reads
\begin{equation}
H_\mathrm{eff}=\mathbb{P} H^0 \mathbb{P}+ \mathbb{P}H'\mathbb{P}+\mathbb{P}H'\mathbb{Q}\frac{1}{E-\mathbb{Q}H^0 \mathbb{Q}-\mathbb{Q}H' \mathbb{Q}}\mathbb{Q}H' \mathbb{P},
\end{equation}
and $\mathbb{P} H^0 \mathbb{P}=E_P$ is the degenerate energy in model space. The effective Hamiltonian to second-order perturbation is given by
\begin{equation}\label{expanofHeff}
H_\mathrm{eff}=E_P+H^{(1)}+H^{(2)},
\end{equation}
where
\begin{equation}\label{pertofHeff}
H^{(1)}=\mathbb{P} H' \mathbb{P}=V_e^1+V_e^2+H_{ee},\quad H^{(2)}=PH_{eb} \frac{\mathbb{Q}}{E_P-H^0}H_{eb}\mathbb{P} .
\end{equation}







\subsection{\texorpdfstring{$C_4$}{} symmetry}
In the simplest case, all the inter-chain couplings $K_i$ are identical, i.e., $K_i=K$. The system has $C_4$ symmetry in the quasi-two-dimensional plane, which in turn results in the degenerate subspaces of odd-parity states. The junction Hamiltonian reads (in the following content we only consider the single-excitation case)
\begin{equation}
H_I^\mathrm{id}=K\sum_j (\xi_{j,N}+\xi_{j,N+1})|1,\psi_j\rangle \sum_{j'}(\xi_{j',N}+\xi_{j',N+1})\langle 2,\psi_{j'}|+\mathrm{H.c.},
\end{equation}
where the index $i=1,2$ in $|i,\psi_j\rangle$ denotes the $i$-th chain. Due to the inversion symmetry of the system, the coefficient $\xi_{j,N+1}=(-1)^{p_j}\xi_{j,N}$ are associated with odd or even parity $p_j$ of $|\psi_j\rangle$. Thus, there are interaction terms between symmetric eigenstates, and anti-symmetric states are dark in degenerate subspaces, as shown in Fig.~\ref{FigS2}(a). The parity of bulk states is given by
\begin{equation}\label{parityofBS}
\begin{aligned}
\hat{\mathcal{P}}|v_{k_j},-\rangle&=\frac{1}{\sqrt{N}}\sum_n \sin{[nk_j-\phi(k_j)]}|N+1-n,B\rangle- \sin{(nk_j)}|N+1-n,A\rangle\\
&=\frac{1}{\sqrt{N}}\sum_n \sin{[(N+1-n)k_j-\phi(k_j)]}|n,B\rangle- \sin{[(N+1-n) k_j]}|n,A\rangle\\
&=\frac{1}{\sqrt{N}}\sum_n \sin{[j\pi -nk_j]}|n,B\rangle- \sin{[j\pi+\phi(k_j) -nk_j]}|n,A\rangle\\
&=(-1)^j |v_{k_j},-\rangle, \quad \text{Similarly }\hat{\mathcal{P}}|v_{k_j},+\rangle=(-1)^{j+1} |v_{k_j},+\rangle.
\end{aligned}
\end{equation}
According to Eq. \eqref{parityofBS}, for bulk states $|\psi_j\rangle$ in complementary space we have $\hat{\mathcal{P}}|\psi_j\rangle=(-1)^j|\psi_j\rangle, j\neq L,R$, and $\xi_{j,i}=(-1)^{j}\xi_{j,2N+1-i}$.

Therefore we can expand the junction Hamiltonian as
\begin{equation}
\begin{aligned}
H_I^\mathrm{id}&=H_{ee}+H_{eb}+H_{bb}\\
&=g |\psi_{1S}\rangle\langle\psi_{2S}|+K \eta \sum_{n} \xi_{2n,N}\left(|\psi_{1S}\rangle\langle 2,\psi_{2n}|+|\psi_{2S}\rangle\langle 1,\psi_{2n}|\right)+K\sum_{n,n'}\xi_{2n,N}\xi_{2n',N}|1,\psi_{2n}\rangle\langle 2,\psi_{2n'}|+\mathrm{H.c.}
\end{aligned}
\end{equation}
where
\begin{equation}\label{geta}
g=2K\eta^2,\quad \eta=\langle A_{\lfloor N/2\rfloor+1}|\psi_L\rangle=\langle B_{\lfloor(N+1)\rfloor/2}|\psi_R\rangle=\frac{1}{\sqrt{\mathcal{N}_h}}\sinh \lfloor(N+1)/2\rfloor\lambda.
\end{equation}
The $C_4$ symmetry configuration results in the couplings for symmetric edge states $|\psi_{1S}\rangle, |\psi_{2S}\rangle$ and bulk states $|\psi_{2n}\rangle$, while the anti-symmetric edge states $|\psi_{1A}\rangle, |\psi_{2A}\rangle$ are degenerate dark states, as mentioned above. In other words, only symmetric edge states exchange between the two topological chains. For convenience, we omit the anti-symmetric part and denote the projector of model space $\mathcal{M}$ as 
\begin{equation}
\mathbb{P}=|\psi_{1S}\rangle\langle\psi_{1S}|+|\psi_{2S}\rangle\langle\psi_{2S}|.
\end{equation}
\begin{figure}[htbp]
\centering
\includegraphics[width=0.9\linewidth]{FigS2}
\caption{(a) The eigenenergies of XSSH with $C_4$ symmetry. The squares are the degenerate states with odd parity, and the circles are the coupled even-parity states. (b) Second-order perturbation of $\psi_{iS}$ from bulk states. The inset shows perturbation energy with increasing $K$ for $N=5$. Here, $J_1=8$, $J_2=20$, $K_i=K=4$, and the total number of cells in (a) is $2N=10$.}
\label{FigS2}
\end{figure}
With Eqs. \eqref{expanofHeff} and \eqref{pertofHeff}, the effective Hamiltonian is now given by
\begin{equation}
H_\mathrm{eff}^\mathrm{id}=\epsilon_S\cdot \mathbbm{1}_{2\times 2}+\left(\begin{array}{cc}
0 & g \\
g & 0
\end{array}\right)
+H^{(2)},
\end{equation}
where
\begin{equation}\label{secondpert}
\begin{gathered}
H^{(2)}_{11}=H^{(2)}_{22}=K^2\eta^2 \left(\sum_{2n\in \text{lower}} \frac{\xi_{2n,N}^2}{\epsilon_S-\epsilon_{2n}}+ \sum_{2n\in \text{upper}} \frac{\xi_{2n,N}^2}{\epsilon_S-\epsilon_{2n}}\right)=(K\eta)^2 \sum_{n=1}^{N-1} \frac{\xi_{2n,N}^2}{\epsilon_S-(-1)^n\epsilon_{n}},\\
H_{12}^{(2)}=H_{21}^{(2)}=0.
\end{gathered}
\end{equation}
The summation of the lower and upper band can be transformed into a single summation of the lower band because from Eq. \eqref{bulkstates} we can learn that $\xi_{2n,N}^2=\xi_{2N+1-2n,N}^2$. We note that the denominator of the perturbation energy in Eq. \eqref{secondpert} is the order of the large bandgap $\Delta_{\mathrm{gap}}$. As shown in Fig.~\ref{FigS2}(b), the second-order perturbation from bulk states is considerably weak when the number of cells $N$ is odd and can be reasonably ignored. Therefore, together with the anti-symmetric part, the effective Hamiltonian reduces to
\begin{equation}
H_\mathrm{eff}^\mathrm{id}=\widetilde{H}_0+g\left(|\Psi_{1S}\rangle\langle\Psi_{2S}|+\mathrm{H.c.}\right),
\end{equation}
where $\widetilde{H}_0=\sum_{l=1,2}\epsilon_S \left(|\Psi_{lS}\rangle\langle\Psi_{lS}|-|\Psi_{lA}\rangle\langle\Psi_{lA}|\right)$. If we take $|\psi_{1S}\rangle$ as $\mid \uparrow\rangle$, and $|\psi_{2S}\rangle$ as $\mid\downarrow\rangle$, then we obtain the single-qubit Hamiltonian
\begin{equation}
H_\mathrm{eff}^\mathrm{id}=g\sigma_x.
\end{equation}
If we prepare the initial state as $|\psi\rangle_i=|\psi(t=0)\rangle=|\psi_{1S}\rangle$, after the dynamical evolution with time $T=\pi/2g$ the final state will be $|\psi\rangle_f=|\psi(T)\rangle=-i|\psi_{2S}\rangle$. Hence, we have implemented a topological state transfer between two chains in a non-adiabatic process and without global modulations, as shown in the main text. Moreover, by adjusting the strength of $K_i$ we can also control the transfer rates.


\subsection{\texorpdfstring{$\mathbb{Z}_2$}{} symmetry}
In the rotationally symmetric configuration, i.e., $K_1=K_3=K_+$, $K_2=K_4=K_-$, the system has $\mathbb{Z}_2$ instead of $C_4$ symmetry. The junction interaction is
\begin{equation}
H_I^\mathrm{rot}=\sum_{j,j'}\left[K_+(\xi_{j,N}\xi_{j',N}+\xi_{j,N+1}\xi_{j',N+1})+K_-(\xi_{j,N}\xi_{j',N+1}+\xi_{j,N+1}\xi_{j',N})\right]|1,\psi_j\rangle\langle 2,\psi_j'\rangle+\mathrm{H.c.}.
\end{equation}
In this case, the edge-edge coupling is given by
\begin{equation}
\begin{aligned}
H_{ee}&=\eta^2 K_+(|\psi_{1L}\rangle\langle\psi_{2L}|+|\psi_{1R}\rangle\langle\psi_{2R}|)+\eta^2 K_-(|\psi_{1L}\rangle\langle\psi_{2R}|+|\psi_{1R}\rangle\langle\psi_{2L}|)+\mathrm{H.c.}\\
&=(g_+ + g_-)|\psi_{1S}\rangle\langle\psi_{2S}|+(g_+ - g_-)|\psi_{1A}\rangle\langle\psi_{2A}| + \mathrm{H.c.},
\end{aligned}
\end{equation}
where $g_+=\eta^2 K_+$, $g_-=\eta^2 K_-$. Following the same derivation steps, we obtain the effective Hamiltonian
\begin{equation}
H_\mathrm{eff}^\mathrm{rot}= \widetilde{H}_0 + (g_+ + g_-)|\psi_{1S}\rangle\langle\psi_{2S}|+(g_+ - g_-)|\psi_{1A}\rangle\langle\psi_{2A}| + \mathrm{H.c.},
\end{equation}
which also shows the coupling between anti-symmetric edge states $|\psi_{iA}\rangle$, leading to the transition from dark states to bright states. The $\mathbb{Z}_2$ configuration makes the (anti-)symmetric states interact with each other in their subspaces with different strengths. We note that the $\mathbb{Z}_2$ symmetry opens another interacting channel for the qubit-exchange process, which enables us to construct a two-qubit gate, as described in the main text.

Indeed, symmetry plays a crucial role as a constraint for the system Hamiltonian. In the general case with no symmetry, there are four interacting channels for edge states:
\begin{equation}
H_\mathrm{eff} = \widetilde{H}_0 + g_{SS} |\psi_{1S}\rangle\langle\psi_{2S}| + g_{SA}|\psi_{1S}\rangle\langle\psi_{2A}| + g_{AS} |\psi_{1A}\rangle\langle\psi_{2S}| + g_{AA} |\psi_{1A}\rangle\langle\psi_{2A}| + \mathrm{H.c.},
\end{equation}
where the coupling strengths
\begin{equation}
\begin{aligned}
g_{SS} = (K_1+K_2+K_3+K_4)/2, \quad g_{SA} = (K_1+K_2-K_3-K_4)/2,\\
g_{AS} = (K_1-K_2+K_3-K_4)/2, \quad g_{AA} = (K_1-K_2-K_3+K_4)/2.
\end{aligned}
\end{equation}

\section{Effective spin-exchange model with tunable interactions}\label{SWAPgate}
By omitting the second-order perturbation (the error caused by bulk states is about $1\% $ when $N=5$), we can also solely use the symmetry to analyze the effective Hamiltonian in both cases. Here we restrict the total Hamiltonian $H_\mathrm{tot}$ in the zero-energy manifold $\mathcal{M}$,
\begin{equation}
H_{\mathrm{eff},\mathcal{M}}=\left(\begin{array}{cccc}
0 & \langle\psi_{1L}|H_\mathrm{tot}|\psi_{1R}\rangle & \langle\psi_{1L}|H_\mathrm{tot}|\psi_{2L}\rangle & \langle\psi_{1L}|H_\mathrm{tot}|\psi_{2R}\rangle\\
\langle\psi_{1R}|H_\mathrm{tot}|\psi_{1L}\rangle & 0 & \langle\psi_{1R}|H_\mathrm{tot}|\psi_{2L}\rangle & \langle\psi_{1R}|H_\mathrm{tot}|\psi_{2R}\rangle\\
\langle\psi_{2L}|H_\mathrm{tot}|\psi_{1L}\rangle & \langle\psi_{2L}|H_\mathrm{tot}|\psi_{1R}\rangle & 0 & \langle\psi_{2L}|H_\mathrm{tot}|\psi_{2R}\rangle\\
\langle\psi_{2R}|H_\mathrm{tot}|\psi_{1L}\rangle & \langle\psi_{2R}|H_\mathrm{tot}|\psi_{1R}\rangle & \langle\psi_{2R}|H_\mathrm{tot}|\psi_{2L}\rangle & 0
\end{array}\right).
\end{equation}
The edges in $C_4$ symmetric case gives rise to $\langle\psi_{1L}|H_\mathrm{tot}|\psi_{2L}\rangle=\langle\psi_{1L}|H_\mathrm{tot}|\psi_{2R}\rangle=\langle\psi_{1R}|H_\mathrm{tot}|\psi_{2L}\rangle=\langle\psi_{1R}|H_\mathrm{tot}|\psi_{2R}\rangle$. Therefore, the Hamiltonian with such symmetry can be written as
\begin{equation}
H_{\mathrm{eff},\mathcal{M}}^\mathrm{id}=\left(\begin{array}{cccc}
0 & \epsilon_S & r & r\\
\epsilon_S & 0 & r & r\\
r & r & 0 & \epsilon_S\\
r & r & \epsilon_S & 0
\end{array}\right).
\end{equation}
The eigenstates and eigenenergies are given by $|\psi_{1A}\rangle,|\psi_{2A}\rangle$, $E_A=-\epsilon_S$ as well as $|\psi_S^\pm\rangle=(|\psi_{1S}\rangle\pm \psi_{2S}\rangle)/\sqrt{2}$, $E_{S}^\pm=\epsilon_S\pm 2r$, which is consistent with the result analyzed in the previous section. Similarly, the Hamiltonian in $\mathbb{Z}_2$ symmetry shows $I_1^{-1} H_\mathrm{tot} I_2 = I_2^{-1} H_\mathrm{tot} I_1$, with $I_i|\psi_{iL}\rangle = |\psi_{iR}\rangle$. Thus, we obtain $\langle\psi_{1L}|H_\mathrm{tot}|\psi_{2L}\rangle=\langle\psi_{1R}|H_\mathrm{tot}|\psi_{2R}\rangle$, $\langle\psi_{1L}|H_\mathrm{tot}|\psi_{2R}\rangle=\langle\psi_{1R}|H_\mathrm{tot}|\psi_{2L}\rangle$, and the effective Hamiltonian is then given by
\begin{equation}
H_{\mathrm{eff},\mathcal{M}}^\mathrm{rot}=\left(\begin{array}{cccc}
    0 & \epsilon_S & u & v \\
    \epsilon_S & 0 & v & u \\
    u & v & 0 & \epsilon_S \\
    v & u & \epsilon_S & 0
    \end{array}\right).
\end{equation}
The eigenstates and eigenenergies are given by $|\psi_S^\pm\rangle=(|\psi_{1S}\rangle\pm \psi_{2S}\rangle)/\sqrt{2}$, $E_{S}^\pm=\epsilon_S\pm (u+v)$ and $|\psi_A^\pm\rangle=(|\psi_{1A}\rangle\pm|\psi_{2A}\rangle)/\sqrt{2}$, $E_A^\pm=-\epsilon_S\pm (u-v)$. The hybridized edge states $|\psi_{iS}\rangle$ and $|\psi_{iA}\rangle$ oscillate in their own subspaces.
\begin{table}[ht]
\centering
\caption{four level to two-qubit basis}
\begin{tabular}{ccccc}
\hline\hline\noalign{\smallskip}
energy $\epsilon$& $\epsilon_S$ & $-\epsilon_S$ & $-\epsilon_S$ & $\epsilon_S$ \\
\noalign{\smallskip}\hline\noalign{\smallskip}
edge state & $|\psi_{1S}\rangle$ & $|\psi_{1A}\rangle$ & $|\psi_{2A}\rangle$ & $|\psi_{2S}\rangle$\\
qubit state & $\mid\uparrow\uparrow\rangle$ & $\mid\uparrow\downarrow\rangle$ & $\mid\downarrow\uparrow\rangle$ & $\mid\downarrow\downarrow\rangle$\\
\noalign{\smallskip}
\hline\hline
\end{tabular}
\label{qubitbasis}
\end{table}
We denote them with the two-qubit basis as in Table \ref{qubitbasis}, and thus the effective Hamiltonian becomes
\begin{equation}\label{Heffrot}
H_\text{eff}^\mathrm{rot}=\epsilon_S \sigma_z \otimes \sigma_z + u\sigma_x\otimes\sigma_x-v\sigma_y\otimes\sigma_y,
\end{equation}
where $u=g_+$ and $v=g_-$. This shows an anisotropic Heisenberg model with two spins, enabling the gate operations of two qubits. 

\subsection{SWAP gate of topologically protected qubits}

For any initial state $|\psi\rangle_i=|\psi(0)\rangle$, the final state after coherent evolution can be expressed as
\begin{equation}\label{gatematrix}
|\psi\rangle_f=e^{-i H_\text{eff}^\mathrm{rot}t/\hbar}|\psi(0)\rangle=
\frac{1}{2}\left(\begin{array}{cccc}
U_S^+ + U_S^- & 0 & 0 & U_S^+ - U_S^- \\
0 & U_A^+ + U_A^- & U_A^+ - U_A^- & 0 \\
0 & U_A^+ - U_A^- & U_A^+ + U_A^- & 0 \\
U_S^+ - U_S^- & 0 & 0 & U_S^+ + U_S^-
\end{array}\right)|\psi\rangle_i,
\end{equation}
where $U_\nu^\pm(t)=e^{-i E_\nu^\pm t/\hbar}$, $\nu=S,A$. If we let $|\psi(0)\rangle=|\psi_{2A}\rangle=\mid\downarrow\rangle\otimes\mid\uparrow\rangle=(0,0,1,0)^T$, the final state is given by
\begin{equation}
|\psi(t)\rangle=\left(0,\frac{U_A^+(t)-U_A^-(t)}{2},\frac{U_A^+(t)+U_A^-(t)}{2},0\right)^T.
\end{equation}
We regard $U_\nu^\pm(t)$ as rotating vectors in the complex plane with different velocities, which results in the interference of different components of $|\psi\rangle_f$. The probabilities for spin states are 
\begin{equation}
P(\mid \downarrow\uparrow\rangle \text{ or } \mid \uparrow\downarrow\rangle)=\frac{1\mp  \cos\Delta_{A} t}{2},
\end{equation}
where $\Delta_{A}=E_A^+ - E_A^- =2(u-v)$. For each time period  $T_A=\pi/\Delta_{A}$, there is a spin-flipping between the two qubits. Similarly, the time period of double-spin-flipping in subspace $\Pi$ is $T_S=\pi/\Delta_S=\pi/2(u+v)$. Due to the different time periods in the two subspaces, we can implement a SWAP gate between the two qubits, where the rotating vectors should fulfill $U_S^+=U_S^-, U_A^+=-U_A^-, U_S^+=U_A^+$, and thus
\begin{equation}\label{TSWAP}
T_\mathrm{SWAP}=\frac{2n\pi}{2(u+v)}=\frac{m\pi}{2(u-v)}=\frac{2k\pi}{2(\epsilon_S+v)},\quad n,m,k\in \mathbb{Z}.
\end{equation}
The appropriate choice is $n=m=1$, i.e., $u=3v$ and $\epsilon_S=(4k-1)v$, which implies that, for instance if $k=2$, it's determined by
\begin{figure}[t]
\centering
\includegraphics[width=0.85\linewidth]{FigS3.pdf}
\caption{(a) Gate fidelity for different $K_+$ and $K_-$. The bright regions represent different $k$ that correspond to a perfect SWAP gate. (b) Fidelity for $K_-=K_+/3$ and varying $K_+$. The red dots are the fidelity for each $k$ and $\bar{\mathcal{F}}_k>0.999$. Here we consider $N=5$ cells for each chain, and $J_1=0.51 J_2$.}
\label{FigS3}
\end{figure}
\begin{equation}
K_+=3K_-, \quad \epsilon_S=7g_-. 
\end{equation}
This determines the sweet point $(K_+^0,K_-^0)$ in the main text for maximum fidelity in $k=2$ region. The latter condition can be approximately achieved by setting $J_1=7K_-$ when $N$ is odd, since from Eqs. \eqref{E_SE_A}, \eqref{apsolutions} and \eqref{geta} we derive
\begin{gather}
\epsilon_S\simeq J_2\frac{\frac{J_2}{J_1}-\frac{J_1}{J_2}}{\left(J_2/J_1\right)^{N+1}}=J_1\left(\frac{J_1}{J_2}\right)^{N-1}\frac{J_2^2-J_1^2}{J_2^2},\\
v=g_-\simeq K_-\left(\frac{J_2}{J_1}\right)^{N-1}(1-\frac{J_1^2}{J_2^2})=K_- \left(\frac{J_1}{J_2}\right)^{N-1}\frac{J_2^2-J_1^2}{J_2^2},
\end{gather}
and thus $\epsilon_S/g_- \simeq J_1/K_- $. Fig.~\ref{FigS3}(a) shows different regions of high-fidelity for $K_+=3K_-$, $\epsilon_S=(4k-1)g_-$, $k\in \mathbb{Z}$. The phase diagram is symmetric under the interchange $K_+\leftrightarrow K_-$, manifesting the $\mathbb{Z}_2$ symmetry of the system. In particular, when $K_+=K_-$, there is a degenerate between $\psi_{A2}$ and $\psi_{A1}$ and thus we have $U_A^+=U_A^-$ and $\bar{\mathcal{F}}=0.5$. The fidelity approaches 1 $(\bar{\mathcal{F}}_k>0.999)$ corresponding to different $k$, as shown in Fig.~\ref{FigS3}(b). We note that all the data is calculated by using the initial Hamiltonian $H_\mathrm{tot}$, which demonstrates excellent consistency with our perturbation theory and verifies the validity of the SWAP gate.
\begin{figure}[htbp]
\centering
\includegraphics[width=0.9\linewidth]{FigS4.pdf}
\caption{Gate time $T_\mathrm{SWAP}$ for different strength of $J_2$. The gate time reduces quickly ($T\sim K_+^{-N}$) with the increasing junction interactions while retaining the fidelity $\bar{\mathcal{F}}>0.99$. The optimal $J_1$ obeys a nearly linear relationship with $K_+$. Here, we fix $K_-=K_+/3$.}
\label{FigS4}
\end{figure}

The gate time $T_\mathrm{SWAP}$ is sensitive to $K_\pm$ and $J_1/J_2$, as the edge-edge coupling is proportional to the junction interactions and probability density $\eta^2$, as shown in Fig.~\ref{FigS4}. Here we fix $K_-=K_+/3$ to show how $T_\mathrm{SWAP}$ changes with $K_+$ for different values of $J_2$, where $J_1$ is chosen to maximize the fidelity. Approximately, the optimal value for $J_1$ is $J_1\simeq 7K_+/3$, and together with Eq. \eqref{apsolutions} we have
\begin{equation}
T_\mathrm{SWAP}^{-1}\sim K_+(\frac{J_1}{J_2})^{N-1}\left( 1-\frac{J_1^2}{J_2^2} \right)\sim K_+^N(1-\alpha), \quad \alpha=\left(\frac{7K_+}{3J_2}\right)^2< 1.
\end{equation}
When $\alpha \ll 1$, the gate time decreases in power law about $K_i$, i.e., $T_\mathrm{SWAP}\sim K_+^{-N}$. We find that there is a sufficient range of parameters available to achieve a short-time gate with high fidelity. Likewise, the effective spin Hamiltonian can be also employed to realize the iSWAP gate with a relative phase $U_A^+=e^{-i\pi/2}U_S^+$, corresponding to $\epsilon_S=4kv$, $k\in\mathbb{Z}$.





\subsection{Suppressing the effects from disorder with \texorpdfstring{$K_i$}{}}
One of the most exotic features of topological systems is their robustness when they confront disorder. In the SSH model, the edge states are protected by the chiral symmetry and robust to off-diagonal disorder. However, the disorder can change the energy levels of the edge states, resulting in the destruction of condition Eq. \eqref{TSWAP}. In this sense, the edge-state transfer and quantum gate are subject to the disorder~\cite{lang2017topological,longhi2019landau} due to the energy shifts.

In our XSSH model, the energy levels of the edge states are also related to the junction interactions. Therefore, by tuning $K_i$ we can reconstruct the SWAP gate for the system disorder. To model the off-diagonal disorder in the Hamiltonian, we add $\delta_i$ sampled from a uniform distribution over the range $[-\delta, \delta]$ to $J_i$ such that $J_i\rightarrow J_i+\delta_i$. For a single SSH chain with the off-diagonal disorder, we denote the Hamiltonian as $H_\mathrm{SSH}'$, where the chiral symmetry is still preserved. Therefore, we denote the shifted energy as $\epsilon_{S'}=-\epsilon_{A'}$ for edge states $|\psi_{S'}\rangle$ and $|\psi_{A'}\rangle$. Also, the zero-energy unhybridized edge states $|\psi_{L'}\rangle$ and $|\psi_{R'}\rangle$ of $H_\mathrm{SSH}'^0$ can still be used as the unperturbed basis. Their wave functions are protected by chiral symmetry, i.e., $\langle \psi_L|\psi_{L'}\rangle \approx \langle \psi_R|\psi_{R'}\rangle \approx 1$, and exclusively occupy one sublattice:
\begin{equation}
\langle B_i | \psi_{L'}\rangle = \langle A_i | \psi_{R'}\rangle =0.
\end{equation}
The interaction Hamiltonian (we consider $\mathbb{Z}_2$ symmetry) in the presence of disorder becomes
\begin{equation}
H_I'=H_I^\mathrm{rot}=\sum_{i,j}\left[K_+(\xi_{i,N}'\xi_{j,N}'+\xi_{i,N+1}'\xi_{j,N+1}')+K_-(\xi_{i,N}'\xi_{j,N+1}'+\xi_{i,N+1}'\xi_{j,N})'\right]|1,\psi_i'\rangle\langle 2,\psi_j'\rangle+\mathrm{H.c.},
\end{equation}
and the edge-edge coupling arising from $H_I'$ is then given by
\begin{equation}\label{eeCoupDis}
\begin{aligned}
H_\mathrm{ee}'&=\sum_{\mu} (K_+\eta^2 +\Delta_{\mu,\mu})|\psi_{1\mu}\rangle\langle\psi_{2\mu}|+\sum_{\mu\neq\nu}(K_-\eta^2+\Delta_{\mu,\nu})|\psi_{1\mu}\rangle\langle\psi_{2\nu}|+\mathrm{H.c.}, \quad \mu,\nu \in \{L',R'\}\\
&=(g_+ + g_-)|\psi_{1S'}\rangle\langle\psi_{2S'}|+(g_+ - g_-)|\psi_{1A'}\rangle\langle\psi_{2A'}| + \mathrm{H.c.}+H_\mathrm{disorder}',
\end{aligned}
\end{equation}
where $\Delta_{\mu,\nu}\ll K_i \eta^2$ is caused by the small deviation of the edge-state wave functions at the central lattice. The impact of disorder is denoted by $H_\mathrm{disorder}'$, which contains couplings of hybridized edge states with strength $\sim \Delta_{\mu,\nu}$. The effective Hamiltonian in the two-qubit basis is governed by (the bulk-edge energy gap remains sufficiently large to guarantee the validity of perturbation theory)
\begin{equation}\label{HeffrotDis}
H_\mathrm{eff}^\mathrm{\prime rot}\simeq\left(\begin{array}{cccc}
t_1 & 0 & 0 & g_+ + g_- \\
0 & -t_1 & g_+ - g_- & g_+ \\
0 & g_+ - g_- & -t_2 & 0 \\
g_+ + g_- & g_+ & 0 & t_2
\end{array}\right),
\end{equation}
where $t_1$ and $t_2$ are the energies of hybridized edge states without junction interactions. In other words, the inversion symmetry is broken due to the disorder such that the $\mathbb{Z}_2$ symmetry is broken even if we have $K_1=K_3=K_+$, $K_2=K_4=K_-$, and thus $\langle \psi_{1S'} |H_\mathrm{ee}'| \psi_{2A'}\rangle \neq 0$, $\langle \psi_{1A'} |H_\mathrm{ee}'| \psi_{2S'}\rangle \neq 0$. However, the chiral-symmetry-protected wave functions ensure that the unexpected couplings are negligible in Eq. \eqref{eeCoupDis}:
\begin{equation}
\langle \psi_{1S'} |H_\mathrm{ee}'| \psi_{2A'}\rangle, \langle \psi_{2S'} |H_\mathrm{ee}'| \psi_{1A'}\rangle \ll \langle \psi_{1S'} |H_\mathrm{ee}'| \psi_{2S'}\rangle, \langle \psi_{1A'} |H_\mathrm{ee}'| \psi_{2A'}\rangle.
\end{equation}




\begin{figure}[t]
\centering
\includegraphics[width=0.8\linewidth]{FigS5.pdf}
\caption{Gate performance in the presence of disorder. (a) Gate fidelity with disorder strength $\delta=0.1J_2$. The maximum fidelity is no longer at $K_+=3K_-$ and $\epsilon_S=7g_-$. However, we can still find the sweet point by adjusting $K_\pm$. (b) The fidelity $\bar{\mathcal{F}}$ of the SWAP gate in $k=2$ region. (c) The gate fidelity without modification. (d) The maximal fidelity $\bar{\mathcal{F}}_m$ by adjusting $K_\pm$. Each point corresponds to the mean value of $5000$ instances in (c) and $100$ instances in (d). The gray and pink regions in (c) and (d) denote the value of a standard deviation.}
\label{FigS5}
\end{figure}

The only difference between Eqs. \eqref{Heffrot} and \eqref{HeffrotDis} is the unperturbed energies of hybridized edge states, while the perfect SWAP gate condition is no longer given by Eq. \eqref{TSWAP} since there is no well-defined $\epsilon_S$. However, the symmetric and anti-symmetric edge states still oscillate independently in the two subspaces. The four edge states are
\begin{equation}
\begin{aligned}|\psi_{S^{\prime}}^+\rangle&\sim|\psi_{1S^{\prime}}\rangle+\lambda_{S}^+|\psi_{2S^{\prime}}\rangle,&|\psi_{S^{\prime}}^-\rangle&\sim|\psi_{1S^{\prime}}\rangle+\lambda_{S}^-|\psi_{2S^{\prime}}\rangle\\|\psi_{A^{\prime}}^+\rangle&\sim|\psi_{1A^{\prime}}\rangle+\lambda_{A}^+|\psi_{2A^{\prime}}\rangle,&|\psi_{A^{\prime}}^-\rangle&\sim|\psi_{1A^{\prime}}\rangle+\lambda_{A}^-|\psi_{2A^{\prime}}\rangle,
\end{aligned}
\end{equation}
where $\lambda_{A'/S'}^\pm$ can be analytically expressed from \eqref{HeffrotDis}. When $|t_1-t_2|\ll g_+-g_-$ we obtain $\lambda_{A'/S'}^\pm\rightarrow 1$, which reduces to the previous results. Here we denote the new basis for two-qubit state in Table \ref{qubitbasis} as
\begin{equation}
|\widetilde{\psi}_{1S'/2S'}\rangle=\frac{|\psi_{S'}^+\rangle\pm|\psi_{S'}^-\rangle}{\sqrt{2}}, \quad 
|\widetilde{\psi}_{1A'/2A'}\rangle=\frac{|\psi_{A'}^+\rangle\pm|\psi_{A'}^+\rangle}{\sqrt{2}}.
\end{equation}
In this case, the evolution matrix with rotating vectors $U_{\nu'}^{\pm}(t)$ has the same form as in Eq. \eqref{gatematrix}. Therefore, the condition $U_{S'}^+= U_{S'}^-$, $U_{A'}^+=- U_{A'}^-$, and $U_{S'}^+=U_{A'}^+$ can be simultaneously satisfied by continuously adjusting the junction interactions $K_\pm$, enabling us to achieve a perfect SWAP gate in the presence of disorder. As shown in Figs.~\ref{FigS5}(a) and \ref{FigS5}(b), we find that the gate fidelity can still reach 0.99 although the maximum is no longer at $K_-=K_+/3$. The diagram slightly deviates from Fig.~\ref{FigS3}(a) with smaller high-fidelity regions together with the broken $\mathbb{Z}_2$ symmetry due to the disorder. In Figs. \ref{FigS5}(c) and \ref{FigS5}(d) we show the fidelity as a function of the disorder strength up to $0.1J_2$, where $K_\pm^0$ denotes the sweet point without disorder as in the main text. When the disorder strength $\delta$ is greater than $K_-^0$, the mean value of fidelity decays rapidly if the junction interactions $K_\pm$ remain unchanged, as shown in Fig.~\ref{FigS5}(c). In contrast, Fig.~\ref{FigS5}(d) reveals that the disorder on $J_i$ is suppressed and we can always achieve a satisfactory SWAP gate by modifying $K_\pm$.

 

\section{Topological edge states interact with 1D electromagnetic environment}
We consider that the atoms of each SSH chain are fabricated homogeneously to couple with a one-dimensional (1D) waveguide. As each chain interacts with its respective electromagnetic environment, we concentrate on the dissipative dynamics of one chain, while the coherent dynamics analyzed earlier remain unchanged. The Hamiltonian describing the SSH atomic chain coupled to a 1D waveguide is given by
\begin{equation}
H=H_\mathrm{A}+H_\mathrm{EM}+H_\mathrm{int}=H_0+H_\mathrm{SSH}+H_\mathrm{EM}+H_\mathrm{int},
\end{equation}
where $H_0=\sum_j \hbar\omega_a \sigma_j^+\sigma_j^-$ is the bare Hamiltonian of two-level atoms with transition frequency $\omega_a$, and the Hamiltonian for electromagnetic field is
\begin{equation}
H_{\mathrm{EM}}=\int_{-\infty}^\infty dk~ \hbar \omega(k) a_k^\dagger a_k.
\end{equation}
The general form of atom-filed (including artificial atom) interaction is described by the dipole coupling Hamiltonian
\begin{equation}
H_\mathrm{int}=\hbar \bar{g} \sum_j \int_{-\infty}^\infty d k ~i\sqrt{\omega(k)}(a_k e^{ikx_j} - a_k^\dagger e^{-ikx_j})(\sigma_j^+ + \sigma_j^-),
\end{equation}
where the coupling strength $\bar g$ has the dimension of square root of velocity, i.e., $[\bar{g}^2]=[v]$. Additionally, we consider the propagating field has linear dispersion $\omega=|k| v_g$ here ($v_g$ is the group velocity)~\cite{chang2012cavity}. We take the electromagnetic field as the large bath and transform the system into the interaction picture. Thus, the dynamics of the SSH chain coupled to the 1D waveguide is given by the Nakajima-Zwanzig equation~\cite{gardiner2004quantum}
\begin{equation}\label{mastereq}
\dot{\tilde{\rho}}_A(t)\simeq -\frac{1}{\hbar^2} \int_0^{\infty} d \tau~ \mathrm{Tr_B}\left\{\left[H_{\mathrm{int}}(t),\left[H_{\mathrm{int}}(t-\tau), \tilde{\rho}_{\mathrm{A}}(t) \otimes \rho_{\mathrm{B}}(0)\right]\right]\right\}
\end{equation}
where
\begin{gather}
\tilde{\rho}_A(t)=e^{i(H_\mathrm{A}+H_\mathrm{EM})t/\hbar}\rho_A(t) e^{-i(H_\mathrm{A}+H_\mathrm{EM})t/\hbar},\\
H_\mathrm{int}(t)=\hbar \bar{g} \sum_j \int_{-\infty}^\infty d k ~i\sqrt{\omega_{|k|}}\left[a_k e^{i(k x_j-\omega t)} - a_k^\dagger e^{-i(kx_j-\omega t)}\right](\sigma_j^+ e^{i\omega_a t} + \sigma_j^- e^{-i\omega_a t}).
\end{gather}
Here, we have assumed that the atomic hopping $J_i$ in SSH chains are much smaller than the frequency $\omega_a$, and thus the energy spectrum of $H_A$ is nearly localized at $\omega_a$ in the single-excitation subspace. We substitute the interaction Hamiltonian $H_\mathrm{int}(t)$ into \eqref{mastereq} and note that for vacuum electromagnetic environment, only four terms containing non-trivial bath correlation function $\mathrm{Tr_B}[a_k a^\dagger_{k^\prime}\rho_B(0)]=\langle a_k a_{k'}^\dagger\rangle_B=\delta(k-k')$ contribute to the integral. One of the terms reads
\begin{equation}\label{integralofoneterm}
\begin{aligned}
I_1=&-\frac{1}{\hbar^2}\int_0^\infty d\tau~ \mathrm{Tr_B}\{H_\mathrm{int}(t)H_\mathrm{int}(t-\tau)\tilde{\rho}_A(t)\otimes\rho_B(0)\}\\
=&-\bar{g}^2 \sum_{jj'}\int_{-\infty}^\infty dk~ \omega_{|k|} e^{ik(x_j-x_{j'})}\int_0^\infty d\tau~ \left[e^{i(\omega_a-\omega)}\sigma_j^+\sigma_{j'}^- +e^{-i(\omega_a+\omega)\tau}\sigma_{j}^- \sigma_{j'}^+ \right]\tilde{\rho}_A(t)\\
=&-\bar{g}^2\sum_{jj'}\int_{0}^\infty  \frac{d\omega}{v_g}~2\omega\cos{\omega(x_j-x_{j'})/v}\times \left[\zeta(\omega_a-\omega)\sigma_{j}^+\sigma_{j'}^- +\zeta(-(\omega_a+\omega))\sigma_{j'}^-\sigma_{j}^+ \right]\tilde{\rho}_A(t)\\
=&-\bar{g}^2\sum_{jj'} \frac{2\pi\omega_a}{v_g}\cos{\omega_a|x_j-x_{j'}|/v_g}\sigma_j^+ \sigma_{j'}^-\tilde{\rho}_A(t)-i\bar{g}^2\sum_{jj'} \hat F(x_j,x_j',\omega_a)\tilde{\rho}_A(t),
\end{aligned}
\end{equation}
with
\begin{equation}
\hat{F}(x_j,x_j',\omega_a)=\mathrm{P}\int_0^\infty \frac{d\omega}{v_g}~2\omega\cos{\omega|x_j-x_{j'}|/v_g}\left(\frac{1}{\omega_a-\omega}\sigma_j^+\sigma_{j'}^- -\frac{1}{\omega_a+\omega}\sigma_{j'}^- \sigma_{j}^+\right),
\end{equation}
where we have made the replacement
\begin{equation}
\zeta(\pm x)=\int_0^{\infty} d \tau e^{\pm i x \tau}=\pi \delta(x)\pm i \mathrm{P}\left(\frac{1}{x}\right).
\end{equation}
To illustrate the meanings of the components in \eqref{integralofoneterm}, we discuss two different cases: $j'=j$ and $j'\neq j$, which correspond to the effects on every single atom and mutual coupling of atoms brought by the 1D vacuum electromagnetic field environment:
\begin{itemize}
    \item For $j=j'$, it reduces to
    \begin{equation}
    \begin{aligned}
    \hat{F}(x_j,x_j',\omega_a)&=\mathrm{P}\int_0^\infty \frac{d\omega}{v_g}~2\omega\left[\frac{1}{\omega_a-\omega}\sigma_j^+\sigma_{j}^- -\frac{1}{\omega_a+\omega}(\sigma_j^+\sigma_{j}^- -1)\right]\\
    &=\mathrm{P}\int_0^\infty \frac{d\omega }{v_g} \frac{4\omega_a \omega}{\omega_a^2 -\omega^2} \sigma_j^+\sigma_{j}^-,
    \end{aligned}
    \end{equation}
    where we have dropped the constant term that only results in a phase shift on $\tilde{\rho}_A$. We separate the integral $I_1$ into real and imaginary parts to show the dissipative and coherent dynamical properties:
    \begin{equation}
    I_1= \sum_j  \left(\mathfrak{Re}\{I_1\}_j+ i  \mathfrak{Im}\{I_1\}_j \right)\sigma_j^+ \sigma_j^- \tilde{\rho}_A(t).
    \end{equation}
    The real part of the integral characterizes the spontaneous decay rate of each atom,
    \begin{equation}
    \mathfrak{Re}\{I_1\}_j =- \bar{g}^2\cdot 2\pi   k_a= -\gamma_0/2,
    \end{equation}
    where $k_a=\omega_a/v_g$ is the wave vector at atom frequency. And the imaginary part represents the Lamb shift on atomic transition frequency (we can redefine the atomic transition frequency to absorb this shift),
    \begin{equation}
    \mathfrak{Im}\{I_1\}_j=-\bar{g}^2 P\int_0^{\omega_c} \frac{d\omega}{v_g} \frac{4\omega_a \omega}{\omega_a^2 -\omega^2},
    \end{equation}
    where $\omega_c$ is the cut-off frequency by taking account of the relativistic effects at high energy.
    
    \item For $j\neq j'$, applying the Kramers-Kronig relation we obtain
    \begin{equation}
    \begin{aligned}
    \hat{F}(x_j,x_j',\omega_a)&=\mathrm{P}\int_0^\infty \frac{d\omega}{v_g}~ \cos{\omega|k_j-k_{j'}|/v_g}\frac{4\omega^2}{\omega_a}\frac{\omega_a}{\omega_a^2-\omega^2}\sigma_j^+\sigma_{j'}^- \\
    &=2\pi k_a \sin{k_a|x_j-x_{j'}|} \sigma_j^+\sigma_{j'}^-.
    \end{aligned}
    \end{equation}
    The integral manifests the collective dynamics of the atoms, with the position-dependent correlated decay and non-local interaction:
    \begin{align}
    & \gamma_{jj'}= -2~\mathfrak{Re}\{I_1\}_{jj'}  = \gamma_0\cos{k_a|x_j-x_{j'}|},\\
    & g_{jj'}=-~\mathfrak{Im}\{I_1\}_{jj'}= \gamma_0\sin{k_a|x_j-x_{j'}|}/2.
    \end{align}
\end{itemize}
In the case of real atoms coupled to the 3D electromagnetic vacuum field, the interaction Hamiltonian reads
\begin{equation}
H_\mathrm{int}=- \sum_j \boldsymbol{\wp}_j \cdot \boldsymbol{E}(\boldsymbol{r})= -\sum_j  \sum_\lambda \int d^3 \boldsymbol{k}~ \sqrt{\frac{\hbar \omega_{|\boldsymbol{k}|}}{2 \varepsilon_0(2 \pi)^3}} \boldsymbol{\epsilon}_{\lambda,\boldsymbol{k}}(a_{\lambda, \boldsymbol{k}} e^{i \boldsymbol{k} \cdot \boldsymbol{r}}-a_{\lambda, \boldsymbol{k}}^{\dagger} e^{-i \boldsymbol{k} \cdot \boldsymbol{r}}) (\boldsymbol{\wp}_j^{eg}\sigma^-+\boldsymbol{\wp}_j^{ge}\sigma^+),
\end{equation}
where $\boldsymbol{\wp}_j=\boldsymbol{\wp}_j^{eg}\sigma^-+\boldsymbol{\wp}_j^{ge}\sigma^+$ is the dipole moment operator and $\boldsymbol{\epsilon}_{\lambda,\boldsymbol{k}}$ denotes the polarization vectors. The coupling strength is expressed as
\begin{equation}
\bar{g}^2= \sum_\lambda \int d\Omega~ \frac{1}{2\hbar \varepsilon_0(2 \pi)^3} |\boldsymbol{\wp}_j^{eg}|^2 \sin{\theta}^2,
\end{equation}
where $\Omega$ is the solid angle. Therefore, the spontaneous emission rate of each atom in 3D free space is then given by
\begin{equation}
\Gamma_{3D}= 2\bar{g}^2 \pi \omega_a^3/c^3=\frac{|\boldsymbol{\wp}_j^{eg}|^2 \omega_a^3}{3 \pi \varepsilon_0 \hbar c^3},
\end{equation}
which is consistent with the Wigner-Weisskopf theory. We note that both the rotating wave and counter-rotating wave terms are essential to stimulate the photon-mediated interactions between atoms, which correspond to $\omega_a-\omega$ and $\omega_a+\omega$ terms in the imaginary part of the integral, but only the rotating wave term contributes to the dissipative dynamics since the radiation of real photons depends on the resonant part of frequency. 

After gathering the remaining three terms and transforming them back into the Schrödinger picture, we obtain the master equation for the SSH chain (here we replace $\rho_A(t)$ with $\rho(t)$ for simplicity):
\begin{equation}
\begin{aligned}
\dot{\rho}(t)&=-\frac{i}{\hbar} \left[H_A+H_\mathrm{nl},\rho\right]+\mathcal{L}\rho\\
&=-\frac{i}{\hbar} \left[H_A+H_\mathrm{nl},\rho\right]+\sum_{ij}\frac{\gamma_{ij}}{2}(2\sigma_i^- \rho \sigma_{j}^+ - \sigma_i^+ \sigma_{j}^-\rho - \rho\sigma_i^+ \sigma_{j}^- ),
\end{aligned}
\end{equation}
where the non-local interaction mediated by virtual photons is
\begin{equation}
H_\mathrm{nl}=\sum_{ij} \hbar g_{ij} (\sigma_i^+ \sigma_j^- +\sigma_i^- \sigma_j^+).
\end{equation}
Alternatively, the same result can be derived by using Green's function approach~\cite{dung2002resonant,asenjo2017atom}, with the 1D Green's function for the electric field given by
\begin{equation}
G(x_i,x_j,\omega_a)=-i\frac{1}{2k_a} e^{ik_a|x_i-x_{j}|}.
\end{equation}
In particular, the correlated decay $\gamma_{ij}$ changes the decay behavior of the atomic chain, leading to the emergence of distinct dissipative eigenmodes instead of independent spontaneous emissions. We formulate the Lindblad operator in the form of the inner product,
\begin{equation}
\mathcal{L}\rho= \bm{\sigma}_-^T \bm{\Gamma}\rho \bm{\sigma}_+ -\frac{1}{2}\bm{\sigma}_+^T \bm{\Gamma} \bm{\sigma}_-\rho -\frac{1}{2}\rho \bm{\sigma}_+^T \bm{\Gamma}\bm{\sigma}_-,
\end{equation}
where $\bm{\sigma}_\pm=(\sigma_{1}^\pm,\sigma_{2}^\pm,\cdots,\sigma_{2N}^\pm)^T$, and $\bm{\Gamma}_{ij}=\gamma_{ij}$ is the dissipation matrix. We diagonalize this real symmetric matrix $T^{-1}\bm{\Gamma}T = \bm{\Lambda}$ to find the $2N$ dissipative eigenmodes with associated decay rate $\bm{\Lambda}_{ii}=\Gamma_i$, which exhibits superradiance and subradiance when $\Gamma_i>\gamma_0$ or $\Gamma_i<\gamma_0$. Hence, the Lindblad operator is decoupled to $2N$ individual dissipative channels, with
\begin{equation}
\mathcal{L}\rho=\sum_i \Gamma_i \mathcal{D}[\tilde{\sigma}_i^-]\rho= \sum_i \frac{\Gamma_i}{2}(2\tilde{\sigma}_i^- \rho \tilde{\sigma}_i^+ - \tilde{\sigma}_i^+ \tilde{\sigma}_i^- \rho- \rho \tilde{\sigma}_i^+ \tilde{\sigma}_i^-),
\end{equation}
where the dissipative eigenmode is given by $\tilde{\sigma}_i^- = \sum_{j} T^{-1}_{ij} \sigma_j^-$.

In our system, the atomic chain consists of $2N-2$ bulk atoms and two edge atoms. We set the distance between nearest bulk atoms as $d_B=|x_{j+1}-x_j|=\lambda_a/2$, where $\lambda_a=2\pi/k_a$ denotes the characteristic wavelength, while the two edge atoms are separated from their nearest bulk atoms by a distance of $d_{EB}=|x_2-x_1|=|x_N-x_{N-1}|=\lambda_a/4$. In this case, the dissipative dynamics of edge atoms are decoupled with the bulk atoms, 
\begin{equation}
\bm{\Gamma}=\left(\begin{array}{cc}
\bm\Gamma_E & 0 \\
0 & \bm\Gamma_B
\end{array}\right),\quad \text{where} \quad
\bm\Gamma_E=\gamma_0\left(\begin{array}{cc}
1 & -1 \\
-1 & 1
\end{array}\right),
\end{equation}
and the dissipative eigenmodes for edge atoms are
\begin{equation}
|\Gamma_{E,+}\rangle=\frac{|1,A\rangle+ |N,B\rangle}{\sqrt{2}},\quad |\Gamma_{E,-}\rangle=\frac{|1,A\rangle- |N,B\rangle}{\sqrt{2}},
\end{equation}
with $\Gamma_{E,-}=2\gamma_0$ and $\Gamma_{E,+}=0$, which represent the dissipative bright and dark states. Similarly, there exists one collective bright state for bulk atoms:
\begin{equation}
|\Gamma_B\rangle=\frac{1}{\sqrt{2N-2}}\sum_{i=2}^{2N-1}(-1)^i \sigma_i^+ |G\rangle,
\end{equation}
with $\Gamma_B=(2N-2)\gamma_0$, and the other bulk eigenmodes are dark states. In this sense, the ensemble of bulk atoms behaves as a superatom with a broadening spectrum. Additionally, as shown in Fig.~\ref{FigS6}(a), both of the two bright states exhibit odd parity (while the dark state $|\Gamma_{E,+}\rangle$ has even parity). To obtain describe the decay rates of the topological edge states, we write the effective non-Hermitian Hamiltonian: 
\begin{figure}[t]
\centering
\includegraphics[width=0.9\linewidth]{FigS6.pdf}
\caption{Dissipative dynamic of SSH and XSSH model. (a) The wave functions of the two dissipative eigenmodes. They characterize the decay rates of edge atoms and bulk atoms respectively. Time evolution of the population of (b) edge states and (c) bulk states in SSH chains with the waveguide, where $\psi_j$ is the $j$-th eigenstate in energy ascending order. The gray dash-dotted line shows the standard decay curve $e^{-\gamma_0 t}$. The solid lines calculated by the master equation are in correspondence with the effective decay rates. (d)-(f) show similar dynamics in our XSSH model, where $P_g$ is the ground state population. (d) Energy levels of edge states and ground state. There are exchange interactions between edge states from different chains, with dissipation occurring exclusively in the anti-symmetric states. Here the cell number of a single chain are $N=5$ and $N=3$ in (b)(c) and (e)(f). Other parameters are $\gamma_0=0.035J_2$, $J_1=0.25J_2$. }
\label{FigS6}
\end{figure}
\begin{equation}
H_\mathrm{nh} = H_\mathrm{SSH}  -i\Gamma_{E,-}|\Gamma_{E,-}\rangle\langle \Gamma_{E,-}| - i\Gamma_B|\Gamma_B\rangle\langle \Gamma_B|.
\end{equation}
Here we drop the bare Hamiltonian $H_0$, and the non-local interaction Hamiltonian $H_{nl}$ is also omitted since it can be treated as the disorder in such a configuration of distance. The strength of this interaction is the same order as $\gamma_0$, i.e., $\mathrm{max}(g_{ij})=\gamma_0 \ll J_1,J_2$. The non-Hermitian part brings diagonal and off-diagonal matrix elements for the odd-parity eigenstates of $H_\mathrm{SSH}$, which presents the direct dissipation and dissipative couplings respectively. We can expand $|\Gamma_{E,-}\rangle$ and $|\Gamma_B\rangle$, and then diagonalize $H_\mathrm{nh}$ to give the exact decay rates for any initial state. Nevertheless, due to the similarity between $|\psi_A\rangle$ and $|\Gamma_{E,-}\rangle$, we have $\langle \psi_A|\Gamma_{E,-}\rangle \approx 1, \langle \psi_A|\Gamma_{B}\rangle \approx 0$. Thus, for the edge states, we only keep the diagonal elements and ignore the off-diagonal part. The edge states should decay through both the edge and bulk dissipative channels. Then, we obtain the total decay rate ($\gamma_\mathrm{tot}$) for symmetric edge states
\begin{equation}
\gamma_S = \Gamma_{E,-}|\langle\Gamma_{E,-}|\psi_S\rangle|^2 + \Gamma_{B} |\langle\Gamma_{B}|\psi_S\rangle|^2=0.
\end{equation}
Attribute to its even parity, the symmetric edge state $|\psi_S\rangle$ becomes a perfect dark state that is immune to the decoherence, as shown in Fig.~\ref{FigS6}(b), while the odd parity of $|\psi_A\rangle$ allows it to exhibit superradiance,
\begin{equation}
\gamma_A = \Gamma_{E,-}|\langle\Gamma_{E,-}|\psi_A\rangle|^2 + \Gamma_{B} |\langle\Gamma_{B}|\psi_A\rangle|^2 \approx 2\gamma_0.
\end{equation}
Similarly, there is a superradiant bulk state $|\psi_1\rangle$ and others are subradiant states or perfect dark states as shown in Fig.~\ref{FigS6}(c). We emphasize that this feature comes from the well-defined parity of their wave functions. To characterize such property, we define the parity function as 
\begin{equation}
\mathcal{P}(\psi)=\sum_{j=1}^{N} \left|\langle\psi|\sigma_j^+|G\rangle + \langle\psi|I\sigma_j^+|G\rangle \right|^2-1=\sum_{j=1}^N \left|\langle j,A|\psi\rangle+\langle N+1-j,B|\psi\rangle\right|^2-1,
\end{equation}
where $I$ is the reflection operator defined in Eq. \eqref{symmetryinOBC}. The wave function that has even or odd parity yields $\mathcal{P}=1$ or $\mathcal{P}=-1$, respectively. In Figs. \ref{FigS6}(d)-(f) we show similar dynamics in our XSSH model. The only difference is that there are also couplings between these edge states, where their total population decays with the same rate as in the SSH case.

The disorder can definitely break the inversion symmetry of the SSH chain, resulting in the destruction of such parity. However, as shown in Figs. \ref{FigS7}(a) and \ref{FigS7}(b), the decay rates of topological edge states remain almost unchanged with the increasing disorder. The underlying cause is that the paired edge states are protected by chiral symmetry either, and their wave functions deviate slightly such that $\mathcal{P}(\psi_{S/A}') \approx 1$, which presents the dissipative robustness of the edge states, as shown in Figs. \ref{FigS7}(c) and \ref{FigS7}(d). As regards the bulk states, the subradiant state with higher energy exhibits weaker dissipation as its wave function has fewer zigzags, which therefore contributes only a small overlap with bulk bright states $|\Gamma_B\rangle$, resulting in the small decay rate and even in the presence of large disorder. But when the disorder exists, the parity of bulk states is not well-defined, and their decay rates have larger fluctuations compared with the edge states.

\begin{figure}[t]
\centering
\includegraphics[width=0.75\linewidth]{FigS7.pdf}
\caption{Decay rates and parity of topological edge states in the presence of disorder. (a) and (b) show the total decay rates $\gamma_\mathrm{tot}$ and fluctuations for the symmetric and anti-symmetric edge states, respectively. (c) and (d) exhibit their parity $\mathcal{P}(\psi)$ function. Here, the parameters are $\delta_\mathrm{max}=0.4J_1=0.1J_2$, and total cells $N=5$.}
\label{FigS7}
\end{figure}


\section{Experimental feasibility}

\begin{figure}[t]
\centering
\includegraphics[width=0.88\linewidth]{FigS9.pdf}
\caption{(a) Schematic of transmon qubit array coupled to the 1D waveguide. (b) Structure of the four-node junction. (c) Tunable junction interaction $K_i$ between the central atoms via a microwave cavity bus.}
\label{FigS9}
\end{figure}

Our implementation of this system utilizes a 2D integration structure of superconducting circuits. We construct a frequency-tunable transmon qubit array and tunable couplers to form an SSH chain, with hopping rates $J_1$ and $J_2$ controlled by the couplers~\cite{Kannan2023} (not shown in the figure). The frequency of a transmon qubit is about 5 GHz, and the qubit-coupling strength $J_i/2\pi=0\sim 35$ MHz. As shown in Fig.~\ref{FigS9}(a), the tunable transmon qubits are identically and capacitively coupled to a waveguide with a spatial distance $\lambda_a/2$ (or $\lambda_a/4$ for edge atoms)~\cite{mirhosseini2019cavity,wen2019large}, which induces the individual decay $\gamma_0\sim 1$ MHz. We can change the distance between atoms by effectively tuning their transition frequency~\cite{van2013photon}. As for the central four-node junction interactions, we can employ the tunable couplers~\cite{chen2014qubit,yan2018tunable} to bridge these central atoms, as shown in Figs.~\ref{FigS9}(b). By exchanging the virtual photons, the interaction strength of artificial atoms is about $K_i \sim 20$ MHz. The aforementioned results can be observed with the given parameters. To obtain a faster gate, the intra-chain hopping rates $K_i$ can be tuned $k=1$ regime such that $K_i\simeq J_1$ as in Fig.~\ref{FigS3}.

\section{Switch on/off of the State transfer and SWAP gate}
In our system, the interaction channels between edge states can be established or suspended by tuning the junction interactions.  In superconducting circuits, this control is achieved through control pulses. The quantum state transfer and SWAP gate of edge-state qubits require the switching of $K_i$ and involve distinct phases: switching on, waiting, and switching off. In this context, the total Hamiltonian $H_\mathrm{tot}(t)$ becomes time-dependent, allowing for the manipulation of edge states without the need for adiabatic modulations. The dynamic control of the Hamiltonian provides a more flexible and efficient method for managing quantum state transitions within the system. For example, if the system is initially prepared in a symmetric edge state, i.e., $|\psi\rangle_i=|\psi(0)\rangle=|\psi_{1S}\rangle$ with the $C_4$ configuration, one should proceed with the following steps to apply the edge-state transfer: (i) For $t<0$ with junction interactions $K_i=0$, there is no transition between these edge states, and thus the system remains stationary. (ii) At $t=0$, the junction interactions are adjusted to $K_i=K$, corresponding to $C_4$ symmetry. The system then evolves in the subspace spanned by $|\psi_{1S}\rangle$ and $|\psi_{2S}\rangle$. (iii) When $t=\pi/2g=T_t$, the edge state is transferred to another topological chain. At this point, by turning off the junction interactions, the final stationary state $|\psi\rangle_f=|\psi(T_t)\rangle=|\psi_{2S}\rangle$ is obtained. Similarly, the implementation of the SWAP gate follows the same preparation-evolution-termination process, with a gate time of $t=T_\mathrm{SWAP}$. 

\begin{figure}[t]
\centering
\includegraphics[width=0.98\linewidth]{FigSwitchonoff.pdf}
\caption{(a) Control scheme of the junction interactions. (b) The fidelity of edge-state transfer versus rise/fall time with arbitrary decay rate $\gamma_0$. (c) Gate fidelity for a series of decay rates. Here, we have used the parameters in Fig.~\ref{FigS4}(c): $J_1=0.35 J_2$, $K_+=0.15 J_2$.}
\label{FigSwitchonoff}
\end{figure}

In Fig.~\ref{FigSwitchonoff}, we show the dynamical switching of edge-state transfer and quantum gates with time-dependent junction interactions $K_i(t)$ and dissipation. The parameters satisfying $K_-=K_+$ and $K_-=K_+/3$ correspond to the cases of edge-state transfer and the SWAP gate, respectively. To precisely control the evolution time, the junction interactions $K_i(t)$ should be varied over short switching intervals, including both rise and fall times. Here, we use the master equation presented in the main text to simulate the dynamics of these operations for various rise and fall times $\tau_{r/f}$, where the time-dependent junction interactions are characterized by trapezoidal pulses, as depicted in Fig.~\ref{FigSwitchonoff}(a). The ratio $\tau_{r/f}/T_t$ and $\tau_{r/f}/T_\mathrm{SWAP}$ determine the fidelity of these operations. In Fig.~\ref{FigSwitchonoff}(b), the fidelity of edge-state transfer improves with an increasing $\tau_{r/f}/T_t$ ratio. The transfer fidelity exceeds 0.99 when $\tau_{r/f}/T_t>10$, a condition that can be easily achieved in common artificial systems. In this plot, the transfer fidelity curve is independent of the atomic decay rate $\gamma_0$ since the symmetric edge states are dissipative dark states. In contrast, the anti-symmetric edge states suffer from atomic dissipation, causing the fidelity of the SWAP gate to decrease as the decay rate $\gamma_0$ increases, as shown in Fig.~\ref{FigSwitchonoff}(c).


\bibliography{reference_SM.bib}
\end{document}