\section{S1.\quad Symmetry constraints on the S matrix}

First, we consider $C_3$ rotation symmetry which preserves the valley:
\begin{equation}
\mathcal S = C_3 \mathcal S C_3^{-1},
\end{equation}
where $C_3$ is a cyclic permutation $(a_1,a_1') \rightarrow (a_2,a_2') \rightarrow (a_3,a_3') \rightarrow (a_1,a_1')$ of the incoming modes which are defined in Fig.\ {\color{red} 1}(b) of the main text, and similar for outgoing modes. Next, we discuss the effect of $C_2$ rotation symmetry and time-reversal symmetry $T$. As these symmetries do not conserve the valley, we need to consider both valleys:
\begin{equation}
\begin{pmatrix} b_K \\ b_{K'} \end{pmatrix} = \begin{pmatrix} \mathcal S_K & 0 \\ 0 & \mathcal S_{K'} \end{pmatrix} \begin{pmatrix} a_K \\ a_{K'} \end{pmatrix}.
\end{equation}
In the basis shown in Fig.\ {\color{red} 1}(b) of the main text, $C_2$ rotation symmetry gives
\begin{equation}
\begin{pmatrix} b_{K'} \\ b_K \end{pmatrix} = \begin{pmatrix} \mathcal S_K & 0 \\ 0 & \mathcal S_{K'} \end{pmatrix} \begin{pmatrix} a_{K'} \\ a_K \end{pmatrix},
\end{equation}
such that $\mathcal S_{K'} = \mathcal S_K$. On the other hand, under time-reversal symmetry we have
\begin{equation}
\begin{pmatrix} a_{K'}^* \\ a_K^* \end{pmatrix} = \begin{pmatrix} \mathcal S_K & 0 \\ 0 & \mathcal S_{K'} \end{pmatrix} \begin{pmatrix} b_{K'}^* \\ b_K^* \end{pmatrix},
\end{equation}
such that $\mathcal S_{K'} = (\mathcal S_K)^t$. Hence, the combination $C_2T$ enforces $\mathcal S_K = (\mathcal S_K)^t$.

\section{S2.\quad S matrix without forward scattering}

The $S$ matrix relates valley Hall states that propagate along AB/BA domain walls at the scattering nodes (AA regions) such that $(b,b')^t = \mathcal S (a, a')^t$ with $a,a'$ six incoming modes and $b,b'$ six outgoing modes where the prime distinguishes the two valley Hall states as illustrated in Fig.\ {\color{red} 1}(b) of the main text. In the absence of forward scattering, we find that the most general $S$ matrix consistent with unitarity and $C_3$ and $C_2T$ symmetry is given by
\begin{equation} \label{eq:Smat1}
\mathcal S(\phi,P_{d1}) = e^{i\varphi} \begin{pmatrix}
0 & e^{i \phi} \sqrt{P_{d1}} & e^{i \phi} \sqrt{P_{d1}} & 0 & \sqrt{P_{d2R}} & -\sqrt{P_{d2L}} \\
e^{i \phi} \sqrt{P_{d1}} & 0 & e^{i \phi} \sqrt{P_{d1}} & -\sqrt{P_{d2L}} & 0 & \sqrt{P_{d2R}} \\
e^{i \phi} \sqrt{P_{d1}} & e^{i \phi} \sqrt{P_{d1}} & 0 & \sqrt{P_{d2R}} & -\sqrt{P_{d2L}} & 0 \\
0 & -\sqrt{P_{d2L}} & \sqrt{P_{d2R}} & 0 & -e^{-i\phi} \sqrt{P_{d1}} & -e^{-i\phi} \sqrt{P_{d1}} \\
\sqrt{P_{d2R}} & 0 & -\sqrt{P_{d2L}} & -e^{-i\phi} \sqrt{P_{d1}} & 0 & -e^{-i\phi} \sqrt{P_{d1}} \\
-\sqrt{P_{d2L}} & \sqrt{P_{d2R}} & 0 & -e^{-i\phi} \sqrt{P_{d1}} &-e^{-i\phi} \sqrt{P_{d1}} & 0
\end{pmatrix},
\end{equation}
with $\varphi$ and $\phi$ real, $P_{d1}$ the intrachannel deflection probability, and $P_{d2R}$ ($P_{d2L}$) the probability for interchannel deflections to the right (left). Unitarity gives the conditions $2P_{d1}+P_{d2R}+P_{d2L}=1$ and $P_{d1}=\sqrt{P_{d2R} P_{d2L}}$. This has two solutions: either all probabilities are nonzero with $0<P_{d1} \leq 1/4$ the only independent parameter or $P_{d1}=0$ and either $P_{d2R}$ or $P_{d2L}$ zero, which is equivalent to what we call the zigzag regime below. Hence, we consider the former solution. In this case, the secular equations yields
\begin{equation} \label{eq:secular1}
1 - \lambda^6 + \lambda^2 \left[ \lambda^2 h(\bm k) - h(\bm k)^* \right] \left( 1 - 4 P_{d1} \sin^2 \phi \right) + 2i \lambda^3 ( 2 \sqrt{P_{d1}} \sin \phi )^3=0,
\end{equation}
with $\lambda = e^{i(El/\hbar v+\varphi)}$ and $h(\bm k)=e^{ik_1}+e^{ik_2}+e^{-i(k_1+k_2)}$. If we define $\sin \phi' = 2\sqrt{P_{d1}} \sin \phi$, which always has a solution for $\phi'$ since $0<P_{d1} \leq 1/4$, we can write the secular equation as 
\begin{equation} \label{eq:secular2}
1 - \lambda^6 + \lambda^2 \left[ \lambda^2 h(\bm k) - h(\bm k)^* \right] \cos^2 \phi' + 2i \lambda^3 \sin^3 \phi'=0,
\end{equation}
which is equivalent to the case $P_{d1}=P_{d2R}=P_{d2L}=1/4$ and $\phi \rightarrow \phi'$. Hence, it is reasonable to assume that $\mathcal S(\phi,P_{d1})$ is unitary equivalent to $\mathcal S(\phi',1/4)$. The latter $S$ matrix is given in Eq.\ ({\color{red} 1}) of the main text with $\varphi=0$ as it only shifts the overal energy, and where we drop the prime on $\phi$ from now on. For general $\phi$, Eq.\ \eqref{eq:secular2} has analytical solutions only at $\bar \Gamma$, in which case $h=3$ and we find
\begin{equation}
\lambda_{1\pm} = \pm e^{\mp i \phi}, \qquad \lambda_{2\pm} = \pm \exp \left[ \pm i\arctan \left( \frac{\sin \phi}{\sqrt{4-\sin^2\phi}} \right) \right],
\end{equation}
where the $\lambda_{2\pm}$ are doubly degenerate. Hence, for each network energy period $2\pi \hbar v/l$, there are always two protected nodes at the $\bar \Gamma$ point. This is shown in Fig.\ \ref{fig:sfig1} where we show the network spectrum along high-symmetry lines of the moir\'e Brillouin zone (MBZ) for several values of $\phi$. The same statement also holds at $\bar K$ and $\bar K'$. Furthermore, when $\phi\neq n\pi$ ($n \in \mathbb Z$), we find that the triple degeneracy at the $\bar \Gamma$, $\bar K$, and $\bar K'$ points of the MBZ is lifted, while the bands remain doubly degenerate at these points for all $\phi$, even after including forward scattering, as these crossings are protected by $C_3$ and $C_2 T$ symmetry.
\begin{figure}
\centering
\includegraphics[width=\linewidth]{sfig1}
\caption{Network spectrum in the absence of forward scattering over one energy period $2\pi \hbar v/l$ along high-symmetry lines of the moir\'e Brillouin zone, as shown in the inset of the rightmost panel. Here, we put $\varphi=-\pi \hbar v/6l$ which centers the flatbands for $\phi=\pi/2$ around zero energy,  and solid (dashed) curves correspond to the $K$ ($K'$) valley.}
\label{fig:sfig1}
\end{figure}

We have shown that up to a unitary transformation, the most general $S$ matrix in the absence of forward scattering is given by $\mathcal S(\phi,1/4)$ for which the left and right interchannel deflection amplitudes are equal. Hence, we consider this case from now on and perform the unitary transformation
\begin{equation}
U \mathcal S U^\dag = e^{i\varphi} \begin{pmatrix}
i \sin \phi \, s_0 & \cos \phi \, s_0^t \\
\cos \phi \, s_0 & i \sin \phi \, s_0^t
\end{pmatrix}, \qquad \textrm{with} \qquad s_0 = \begin{pmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{pmatrix},
\end{equation}
where $U = e^{-i\pi \sigma_y/4} e^{i\phi\sigma_z/2} \otimes \mathds 1_3$ transforms $(a,a') \rightarrow (a_+,a_-)$ with $a_\pm = ( a e^{i \phi / 2} \mp a' e^{-i \phi / 2})/\sqrt{2}$ and similar for outgoing modes. We see that for $\phi=n \pi$ ($n \in \mathbb Z$),
\begin{equation}
U \mathcal S U^\dag =  (-1)^n e^{i\varphi} \begin{pmatrix} 0 & s_0^t \\ s_0 & 0 \end{pmatrix},
\end{equation}
such that scattering modes form three independent chiral zigzag channels. Furthermore, in this case we see from Eqs.\ \eqref{eq:secular1} and \eqref{eq:secular2} that $\phi'=\phi$ such that the network supports chiral zigzag modes for any allowed values of the deflection probabilities. On the other hand, for $\phi = (n+1/2) \pi$ ($n \in \mathbb Z$),
\begin{equation}
U \mathcal S U^\dag = (-1)^n i e^{i\varphi} \begin{pmatrix} s_0 & 0 \\ 0 & s_0^t \end{pmatrix},
\end{equation}
such that scattering modes perform closed orbits around AB and BA domains.

\section{S3.\quad S matrix with forward scattering}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{sfig2}
\caption{Network spectrum for the case where zigzag modes along different directions are decoupled ($\phi=0$ and $P_{d1}=P_{d2}$) over one energy period $2\pi \hbar v/l$ along high-symmetry lines of the moir\'e Brillouin zone with $\varphi =-\pi \hbar v/6l$. The bands are shown for $K$ (solid) and $K'$ (dashed) for (a) no forward scattering, (b) $P_{f1}=0.02$ and $P_{f2}=0$, and (c) $P_{f1}=0$ and $P_{f2}=0.02$.}
\label{fig:sfig2}
\end{figure}

When we allow for forward scattering, the $S$ matrix can be written as
\begin{equation} \label{eq:Smat2}
\mathcal S =
e^{i\varphi} \begin{pmatrix}
e^{i(\phi+\chi)} \sqrt{P_{f1}} & e^{i\phi} \sqrt{P_{d1}} & e^{i\phi} \sqrt{P_{d1}} & -\sqrt{P_{f2}} & \sqrt{P_{d2}} & - \sqrt{P_{d2}} \\
e^{i\phi} \sqrt{P_{d1}} & e^{i(\phi+\chi)} \sqrt{P_{f1}} & e^{i\phi} \sqrt{P_{d1}} & -\sqrt{P_{d2}} & -\sqrt{P_{f2}} & \sqrt{P_{d2}} \\
e^{i\phi} \sqrt{P_{d1}} & e^{i\phi} \sqrt{P_{d1}} & e^{i(\phi+\chi)} \sqrt{P_{f1}} & \sqrt{P_{d2}} & - \sqrt{P_{d2}} & -\sqrt{P_{f2}} \\
-\sqrt{P_{f2}} & -\sqrt{P_{d2}} & \sqrt{P_{d2}} & -e^{-i(\phi+\chi)} \sqrt{P_{f1}} & -e^{-i\phi} \sqrt{P_{d1}} & -e^{-i\phi} \sqrt{P_{d1}} \\
\sqrt{P_{d2}} & -\sqrt{P_{f2}} & -\sqrt{P_{d2}} & -e^{-i\phi} \sqrt{P_{d1}} & -e^{-i(\phi+\chi)} \sqrt{P_{f1}} & -e^{-i\phi} \sqrt{P_{d1}} \\
-\sqrt{P_{d2}} & \sqrt{P_{d2}} & -\sqrt{P_{f2}} & -e^{-i\phi} \sqrt{P_{d1}} & -e^{-i\phi} \sqrt{P_{d1}} & -e^{-i(\phi+\chi)} \sqrt{P_{f1}}
\end{pmatrix},
\end{equation}
with $\cos\chi = \left( P_{d2} - P_{d1} \right) / 2 \sqrt{P_{f1} P_{d1}}$ such that $\chi$ is real, so that $2 \sqrt{P_{f1} P_{d1}} \geq | P_{d2}-P_{d1} |$. Note that when $P_{d2}=0$, this condition gives a lower bound on forward scattering $P_{f1} \geq P_{d1}/2$. Here, we assumed that the probability for intrachannel processes is the same for the two valley Hall states. Current conservation then requires $2(P_{d1} + P_{d2}) + P_{f1} + P_{f2} = 1$, where $P_{f1}$ ($P_{d1}$) and $P_{f2}$ ($P_{d2}$) are the probabilities for intra- and interchannel forward scattering (deflections), respectively. To investigate the couplings between zigzag modes, it is preferable to work in the new basis: $(b_+,b_-)^t = U\mathcal SU^{-1}(a_+,a_-)^t$. For $P_d=P_{d1}=P_{d2}$, we find
\begin{equation}
U \mathcal S U^\dag = e^{i\varphi}
\begin{pmatrix}
f \cos \phi & 0 & 2i \sqrt{P_d} \sin \phi & -if^* \sin \phi & 2 \sqrt{P_d} \cos \phi & 0 \\
2i \sqrt{P_d} \sin \phi & f \cos \phi & 0 & 0 & -if^* \sin \phi & 2 \sqrt{P_d} \cos \phi \\
0 & 2i \sqrt{P_d} \sin \phi & f \cos \phi & 2 \sqrt{P_d} \cos \phi & 0 & -if^* \sin \phi \\
if \sin \phi & 0 & 2 \sqrt{P_d} \cos \phi  & -f^* \cos \phi & 2i \sqrt{P_d} \sin \phi & 0 \\
2 \sqrt{P_d} \cos \phi & if \sin \phi  & 0 & 0 & -f^* \cos \phi & 2i \sqrt{P_d} \sin \phi \\
0 & 2 \sqrt{P_d} \cos \phi & if \sin \phi  & 2i \sqrt{P_d} \sin \phi & 0 & -f^* \cos \phi
\end{pmatrix},
\end{equation}
where $f=\sqrt{P_{f2}}+i\sqrt{P_{f1}}$.

In the absence of coupling between zigzag modes propagating along different direction ($\phi=0$ and $P_{d1}=P_{d2}$), the network spectrum is given by $(j=1,2,3)$
\begin{equation} \label{eq:spectrum2}
E_{j,\pm,n}(\bm k) = \frac{\hbar v}{l} \left[ 2 \pi n - \varphi - i \log \frac{1}{2} \left( F_j(\bm k) \pm \sqrt{4e^{i k_j} + F_j(\bm k)^2} \right) \right],
\end{equation}
where $F_j(\bm k) = f^* e^{-ik_{j+1}} - f e^{-ik_{j+2}}$ with $j$ defined cyclically and which is shown in Fig.\ \ref{fig:sfig2}. As can be seen in Fig.\ {\color{red} 2}(c) of the main text, coupling between parallel zigzag channels warps the Fermi surface, in a manner that depends on the type of forward scattering. For $P_{f2}=0$, the bands are symmetric about $k_y$ as in this case $F_1(k_x,-k_y)=F_2(k_x,k_y)$ and $F_3(k_x,-k_y)=F_3(k_x,k_y)$. Moreover, states at the $\bar \Gamma$ and $\pm\bar K$ points in the moir\'e Brillouin zone (MBZ) remain triply degenerate, but are shifted as $F_j(\bar \Gamma) = -2i\sqrt{P_{f1}}$ and $F_j(\pm \bar K) = \pm 2 e^{\pm i\pi/6} \sqrt{P_{f1}}$ for all $j$.

\section{S4.\quad Transport calculation}

We consider a network strip as shown in Fig.\ \ref{fig:sfig3} with length $L = Nl$ or $(N-1/2)l$ ($N=1,2,\ldots$) and width $W \gg L$, where $l \approx 14 \left( \theta^\circ \right)^{-1} \,$nm is the moir\'e lattice constant. To investigate transport, we calculate the transfer matrix that relates amplitudes between two sides, e.g.\ $\psi_{\textrm{right}} = \mathcal M \psi_{\textrm{left}}$. For the transfer matrix, unitarity is expressed as $\mathcal M^\dag J \mathcal M = J$, with $J = \diag \left( 1, 1,-1,-1,1,1\right)$ the current operator for a bulk node, so that $|\textrm{det} \, \mathcal M| =1$.

\subsection{A.\quad Transfer-matrix method}

The total transfer matrix $T$ links one lead to the other lead, e.g.\ $\psi_{N+1} = T \psi_0$ with
\begin{equation} \label{eq:transfer}
T = 
\begin{cases}
B T_N \cdots T_1 & \qquad L = Nl, \\
T_N \cdots T_1 & \qquad L = (N-\frac{1}{2})l,
\end{cases}
\end{equation}
where $T_n = C_n D A_n B$ with $A_n$, $B$, $C_n$, and $D$ defined in Fig.\ \ref{fig:sfig3}. We consider a network strip with $W\gg L$ and introduce the transverse momentum $0\leq k<2\pi/\sqrt{3}l$. If we label the modes from top to bottom, the transfer matrices of the unit cell can be written as $B = \diag \left(\mathds 1_2, \mathcal M \right)$ and
\begin{equation}
D(k) = 
\begin{pmatrix} \mathcal M_{22} & 0 & \mathcal M_{21} e^{ik\sqrt{3}l} \\ 0 & \mathds 1_2 & 0 \\ \mathcal M_{12} e^{-ik\sqrt{3}l} & 0 & \mathcal M_{11} \end{pmatrix}, \qquad \mathcal M = \begin{pmatrix} \mathcal M_{11} & \mathcal M_{12} \\ \mathcal M_{21} & \mathcal M_{22} \end{pmatrix},
\end{equation}
where we defined the submatrices $\mathcal M_{ij}$ of the transfer matrix $\mathcal M$ for a single node. Here, $\mathcal M_{11}$ and $\mathcal M_{22}$ are square matrices of dimension $2$ and $4$, respectively. In the absence of disorder and magnetic fields, $A_n = C_n = \diag \left( \lambda^{-1/2},  \lambda , \lambda^{-1/2}, \lambda \right) \otimes \mathds 1_2$ with $\lambda = e^{i El/\hbar v}$ the dynamical phase, independent of $n$. 

The transmission probability is then calculated by imposing boundary conditions on the leads. For transport in the positive $x$ direction, we have $\psi_0^{(\alpha)} = \left( r_{1\alpha}, a_{1\alpha} , r_{2\alpha}, a_{2\alpha} \right)^t$ and $\psi_{N+1}^{(\alpha)} = \left( 0, t_{1\alpha}, 0, t_{2\alpha} \right)^t$ where $r_i$ ($t_i$) are two-component reflection (transmission) amplitudes, and $\alpha$ labels the incoming modes with amplitudes $a_{11}=a_{23}=(1,0)$ and $a_{12}=a_{24}=(0,1)$ and $(0,0)$ otherwise. The transmission probability for a transverse mode with momentum $k$ is then given by $\mathcal T_k = 1 - (1/4) \sum_{\alpha=1}^4 r_\alpha^\dag r_\alpha$ with $r_\alpha = \left( r_{1\alpha} , r_{2\alpha} \right)^t$. The total transmission probability becomes 
\begin{equation}
\mathcal T = (\sqrt{3}l/2\pi) \int_0^{2\pi/\sqrt{3}l} dk \, \mathcal T_k.
\end{equation}

Up to now, we considered a single valley. For the other valley, the propagation direction of valley Hall states shown in Fig.\ \ref{fig:sfig3} is reversed, as the valleys are related by time reversal. Hence, the transmission for $K'$ in the positive $x$ direction is given by the transmission of $K$ in the negative $x$ direction. It follows that $\mathcal T_K = \mathcal T_{K'}$ in the limit $W \gg L$.
\begin{figure}
\centering
\includegraphics[width=.75\linewidth]{sfig3}
\caption{Unit cell of a network strip with length $L=Nl$ where the dashed vertical lines outline the slice corresponding to the transfer matrix $T_1$ and dotted lines mark subslices with transfer matrices $A$, $B$, $C$, and $D$ as indicated [see Eq.\ \eqref{eq:transfer}].}
\label{fig:sfig3}
\end{figure}

\subsection{B.\quad Magnetotransport}

The magnetic field introduces an additional Peierls phase accumulated by the valley Hall states during propagation along links. In the Landau gauge $\bm A = Bx \bm e_y$, the Peierls phase along a link starting at $x=ml/2$ ($m \in \mathbb Z$) is zero for horizontal links (see Fig.\ \ref{fig:sfig3}) and for an upward or downward diagonal link,
\begin{equation}
\Phi_{\pm}(m) = \mp \pi \left( m + \frac{1}{2} \right) \frac{\Phi}{\Phi_0} ,
\end{equation}
respectively, where $\Phi_0 = h/e$ is the flux quantum and $\Phi = B \mathcal A$ is the flux through the moir\'e cell, comprising an AB and BA triangle, with $B$ the magnetic field and $\mathcal A = \sqrt{3} \, l^2/2$ the moir\'e cell area. If we assume that the magnetic field is sufficiently small such that the structure of the $S$ matrix is unchanged, the transfer matrices $B$ and $D$ remain the same, while $A_n(\Phi) = Q_{2n-2}(\Phi)$ and $C_n(\Phi) = Q_{2n-1}(-\Phi)$ with $Q_m = \diag \left( \lambda^{-1/2},  \lambda e^{i\Phi_+(m)} , \lambda^{-1/2}, \lambda e^{i\Phi_-(m)} \right) \otimes \mathds 1_2$.

\subsubsection*{Quasi-1D regime}

In the quasi-1D regime ($P_d=P_{d1}=P_{d2}$ and $\phi =0$), there are three separate contributions to the conductance from zigzag (ZZ) modes along the $\bm l_{1,2,3}$ directions for a given valley and spin. In the setup shown in Fig.\ \ref{fig:sfig3}, the contribution from $\bm l_3$ is always given by $\mathcal T_{\bm l_3} = 2$, and the transmission for $\bm l_{1,2}$ is identical. The latter vanish in the absence of forward scattering and can be calculated straightforwardly for an infinitely wide system of length $L/l=1/2,1,3/2$ by summing Feynman paths. The results are
\begin{figure}
\centering
\includegraphics[width=.7\linewidth]{sfig4}
\caption{Feynman paths in the quasi-1D regime of the $\bm l_2$ ZZ modes for a network strip of length $L=3l/2$, where $f=\sqrt{P_{f2}} + i \sqrt{P_{f1}}$ is the forward scattering amplitude in the ZZ basis between solid arrows and $\Phi$ is the flux through a moir\'e cell.}
\label{fig:sfig4}
\end{figure}
\begin{alignat}{3} \label{eq:Ga}
& \left. G \right|_{L=\frac{l}{2}} && = \left( 2 + 2 P_f^2 \right) \frac{G_0W}{\sqrt{3}l}, \\
& \left. G \right|_{L=l} && = \left( 2 + \frac{2 P_f^3}{1-P_f (4P_d)^2} \right) \frac{G_0W}{\sqrt{3}l}, \\
& \left. G \right|_{L=\frac{3l}{2}} && = \left( 2 + \frac{2 P_f^4}{1-P_f (4P_d)^2\left[ 2 \cos ( \pi \Phi/\Phi_0 ) \right]^2} \right) \frac{G_0W}{\sqrt{3}l},
\end{alignat}
where $P_f = P_{f1} + P_{f2}$, which is illustrated in Fig.\ \ref{fig:sfig4} for $L=3l/2$. Note that these expressions only depend on the total forward scattering probability $P_f$, which is always the case in the quasi-1D regime. The conductance is shown in Fig.\ \ref{fig:sfig5} at zero magnetic field as a function of $P_f$ for several lengths. For larger systems, the Feynman paths are more involved due to contributions from paths that have segments that cut horizontal along the sample, opposite to the transport direction. If we calculate $\mathcal T_{\bm l_{1,2}}$ up to second order in the path length, we find
\begin{equation} \label{eq:sGa2} 
G \simeq \left( 2 + 2 P_f^{2L/l+1} \left[ 1 + P_f (4P_d)^2 \left( \frac{\sin \left[ \pi (2L/l-1) \Phi/\Phi_0 \right]}{\sin (\pi \Phi/\Phi_0)} \right)^2 \right] \right) \frac{G_0W}{\sqrt{3}l},
\end{equation}
which includes $\mathcal T_{\bm l_3} = 2$ and contributions to $\mathcal T_{\bm l_{1,2}}$ from paths with lengths $2L$ and $2L+3l$, and corresponds to Eq.\ ({\color{red} 6}) of the main text.
\begin{figure}
\centering
\includegraphics[width=.65\linewidth]{sfig5}
\caption{Zero-field conductance for a network strip with width $W \gg L$ in the quasi-1D regime ($\phi=0$) as a function of the total forward scattering probability $P_f = P_{f1}+P_{f2}$ for different lengths $L = Nl$.}
\label{fig:sfig5}
\end{figure}

Remarkably, there exist formal analytical expressions for $L>3l/2$ which can be obtained by explicit calculation of the total transfer matrix and integration over transverse momentum. We find that the momentum-dependent transmission can be written as
\begin{equation}
\mathcal T_{\bm l_{1,2}}(k) = \frac{P_f^{2L/l+1}}{\sum_{n=0}^{m_L} a_n \cos^n k},
\end{equation}
where $m_L= \textrm{floor}(L/l)$ and the coefficients $a_n$ are functions of $P_f$ and the flux $\Phi$. We have observed this explicitly up to $L=2l$ for finite $\Phi$ and $L=7l/2$ for zero $\Phi$ and we expect it holds for all $L$. For example, for $L=2l$ we find that
\begin{align}
a_0 & = 1 - 4P_d P_f (1 + P_f^2) + 8 (4P_d)^2 P_f \cos^2 \frac{\pi \Phi}{\Phi_0} \cos \frac{2\pi\Phi}{\Phi_0}, \\
a_1 & = 2 \sqrt{P_f} \left[ 1 + (P_f - 2) P_f^2 \right] \left( 1 + 2 \cos \frac{2\pi\Phi}{\Phi_0} \right), \\
a_2 & = 16 P_d P_f, 
\end{align}
where $4P_d=1-P_f$. The total transmission function is given by
\begin{equation}
\mathcal T_{\bm l_{1,2}} = \frac{P_f^{2L/l+1}}{2\pi} \int_0^{2\pi} \frac{dk}{\sum_{n=0}^{m_L} a_n \cos^n k} = \frac{P_f^{2L/l+1}}{2\pi i} \oint \frac{(1/z) dz}{\sum_{n=0}^{m_L} \frac{a_n}{2^n} \left( z + 1/z \right)^n},
\end{equation}
where $z=e^{ik}$ and the contour goes along the unit circle. The integral can be worked out and becomes
\begin{equation}
\frac{1}{2\pi i} \oint \frac{z^{m_L-1} dz}{\sum_{n=0}^{2m_L} b_n z^n} = \frac{2^{m_L}}{a_{2m_L}} \left( \sum_{|z_i|<1} + \frac{1}{2} \sum_{|z_i|=1} \right) \frac{z_i^{m_L-1}}{\prod_{j\neq i} \left( z_i - z_j \right)},
\end{equation}
where $z_i$ ($i=1,\ldots,2m_L$) are the roots of the polynomial of degree $2m_L$ in the denominator on the left-hand side of the equation and we assumed that the poles are simple. If this is not the case, the results can be obtained with the general residue theorem. Moreover, up to $L=5l/2$, there exist closed-form expressions for the roots. For $L=2$, we find there are only simple poles, given by
\begin{equation}
z_{\zeta,\xi} = - \frac{a_1 + \zeta \sqrt{a_1^2-4a_0a_2} + \xi \sqrt{2a_1 \left( a_1 + \zeta \sqrt{a_1^2 - 4a_0a_2} \right) - 4a_2 (a_0 + a_2)}}{2a_2}
\end{equation}
where $\zeta,\xi=\pm$. Furthermore, we find there are always two poles inside the unit circle, labeled $z_1$ and $z_2$ (for $\Phi=0$, they are given by $z_{\pm,-}$), and two outside the unit circle, labeled $z_3$ and $z_4$, so that
\begin{equation}
\left. \mathcal T_{\bm l_{1,2}} \right|_{L=2l} = \frac{4P_f^4}{a_2} \left[ \frac{z_1}{(z_1-z_2)(z_1-z_3)(z_1-z_4)} + \frac{z_2}{(z_2-z_1)(z_2-z_3)(z_2-z_4)} \right].
\end{equation}

\subsection{C.\quad Disorder}

We distinguish three types of disorder: (1) disorder sharp on the interatomic scale that couples the valleys, (2) disorder on a larger length scale that couples the two valley Hall channels (for one valley and spin) during propagation along a link, and (3) smooth disorder that adds random phases to the dynamical phase accumulated during propagation along a link. The first type of disorder could, for example, arise from stacking faults where the twist angle changes abruptly. However, the phenomena that we want to address in this work occur in samples with an almost homogeneous twist angle. Therefore, the inclusion of stacking faults in the transport calculation is left for future research. The length scale of the second type of disorder is given by $(\hbar v /\gamma_\perp) \left( \gamma_\perp/U \right)^{1/3}$, which follows from the four-band continuum model for bilayer graphene for the case with two semi-infinite regions of AB ($x<0$) and BA ($x>0$) stacking with $\gamma_\perp \approx 0.3$~eV the interlayer coupling and $U$ the interlayer bias. This type of disorder leads to scattering between the two valley Hall states during propagation along a link. It can be included by letting
\begin{equation}
\lambda \mathds 1_6 \rightarrow \lambda \left[ \cos \tau \, \mathds 1_6 + i \sin \tau (\sigma_y \otimes \mathds 1_3) \right],
\end{equation}
where $\tau$ is the disorder strength. This couples $a_\pm$ modes along a link with probability $\cos^2 \phi \, \sin^2\tau$. Hence, it leads to couplings between ZZ modes similar to the processes shown in Fig.\ {\color{red} 2}(d) of the main text. Finally, disorder due to smooth fluctuations of the Fermi energy or the twist angle, for example, does not couple the valley Hall states along a link but introduces randomness in the scattering parameters at the nodes as well as random phases after propagation along links.

We can include the third type of disorder in the transport calculations by letting $\lambda \rightarrow \lambda e^{i2\pi \gamma(\bm r)}$, where $\gamma(\bm r)$ is a normal-distributed number with zero mean and standard deviation $\sigma$, where we keep the $S$ matrix constant everywhere for simplicity. To this end, we consider a network strip with a finite width $W=\sqrt{3} Ml$, as shown in Fig.\ \ref{fig:sfig7}. We model the edges by introducing an edge $S$ matrix $\mathcal S_b$. For concreteness, we take $U\mathcal S_b U^{-1} = \sigma_x$ so that the edge does not affect $\bm l_3$ zigzag modes, but mixes the $\bm l_1$ and $\bm l_2$ zigzag modes. In any case, in the limit $W \gg L$, the boundary effect becomes negligible. In Fig.\ \ref{fig:sfig8}, we show the disorder-averaged magnetoconductance for a network strip of width $W=60\sqrt{3}\,l$ in the quasi-1D regime for the same parameters as Fig.\ {\color{red} 3}(b) of the main text with $\sigma=0.025$, $0.05$, $0.075$, and $0.1$. Note that except for the case $P_f=0.8$, the background is smaller than its value $8e^2/h$ for the infinitely-wide system, due to scattering at the edges, but which is not affected by disorder given its chiral origin. Nevertheless, we see that the A-B resonances are broadened and that the peak heights are reduced when the disorder strength $\sigma$ is increased. This can be understood from the fact that different contributions with the same path length are individually shifted in $\Phi$ by the disorder. In Fig.\ \ref{fig:sfig9}, we show the height of the main A-B peak versus $\sigma$. 
\begin{figure}
\centering
\includegraphics[width=.8\linewidth]{sfig7}
\caption{Finite network strip of length $L=Nl$ and width $W=2\sqrt{3}l$, where black (gray) dots are bulk (edge) scattering nodes.}
\label{fig:sfig7}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{sfig8}
\caption{Disorder-averaged magnetoconductance for a network strip of length $L = 10 \, l$ and width $W=60\sqrt{3} \,l$ in the quasi-1D regime ($\phi=0$) for several $P_f$ as a function of $\Phi=B\mathcal A$. The disorder strength which is given by the standard deviation $\sigma$ of the random phases is shown on the figures, and the average is over $10$ disorder configurations.}
\label{fig:sfig8}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=.75\linewidth]{sfig9}
\caption{Peak height of the conductance at $\Phi=0$ in the quasi-1D regime ($\phi=0$) as a function of the standard deviation $\sigma$ of the Gaussian-distributed random link phases with zero mean for a network strip of length $L=10\, l$ and width $W=60\sqrt{3} \,l$, averaged over $40$ disorder configurations.}
\label{fig:sfig9}
\end{figure}

\end{document} 