## 4 Visualizing Electron Flow across the Interface  
To visualize the electron flow across the interface, we computed the local current corresponding to a number of specific energies, as indicated in Figure 2. The local current is computed using the current density operator, expressed in the Green's function framework and projected on the atomic basis [30,39]:  
\[I_{ij}=\frac{e}{h}\int_{-\infty}^{\infty}dE\Big{(}t_{ij}G_{ji}^{<}(E)-t_{ji}G_ {ij}^{<}(E)\Big{)} \tag{1}\]  
where \(I_{ij}\) is the current between any two atomic sites, \(i\) and \(j\), \(t_{ij}=t_{ji}^{*}\) is the hopping between these sites and \(G^{<}\) is the lesser GF. Neglecting electronic correlations, the retarded GF (\(G_{R}\)) can be used to obtain \(G^{<}\) by [30]:  
\[G^{<}(E)=G_{R}\Sigma^{<}G_{R}^{\dagger};\qquad\Sigma^{<}=\sum_{l}f_{l}(\Sigma_ {l}^{\dagger}-\Sigma_{l}) \tag{2}\]  
where the sum runs over all the terminals, \(\Sigma_{l}\) is the \(l\)-th terminal self-energy and \(f_{l}\) is the corresponding Fermi distribution that describes the filling of the terminal states up to the chemical potential \(\mu_{l}\). Here, we consider a Heaviside-like distribution, in the zero temperature limit. The real space visualization is obtained by representing the local current data as:  
\[I(x,y)=\sum_{i\
eq j}\alpha_{ij}\langle\mathbf{r}|i\rangle\langle i|I|j\rangle \langle j|\mathbf{r}\rangle \tag{3}\]  
where \(\langle i|I|j\rangle=I_{ij}\) and \(\langle\mathbf{r}|i\rangle=e^{-r^{2}/a_{CC}^{2}}\) with \(r\) being the distance from the point \((x,y)\) to the position of site \(i\) and \(a_{CC}\) is the bond length in graphene. As \(I_{ij}=-I_{ji}\) (see Equation (1)), we select only the positive values of \(I_{ij}\) by setting \(\alpha_{ij}=1\) (\(0\)) when \(I_{ij}\geq 0\) (\(I_{ij}<0\)). Even though the equations for local current are only formally correct for full non-equilibrium, self-consistent calculations, they remain applicable for a qualitative analysis at low-bias potential. While the study of non-equilibrium transport, including inelastic effects, is out of the scope of this work, such a methodology can be expanded to include accurately the effect of finite bias [40,41]. In addition, it is expected that the application of a small bias will not disrupt the spin ordering, as is the case with GNRs conducting high currents [42,43]. For comparison purposes, we also plot the current flow distribution at the same energy for the corresponding homojunction. The plots for the homojunctions confirm that the electrons with energy lying in Region (II) proceed quasi-exclusively along the edges. We now discuss the zGNR-6-based junction in more detail: we note that the current density along the edge is significantly more intense for energies farther away from the Fermi energy. As we will show now, this result is due to the symmetry properties of the eigenstates in the pristine lead. To clarify this point, we examine the local charge density for the highest occupied band at different values of the \(\mathbf{k}\) vector for the zGNR-6 calculated with the TBU approach (Figure 3). The selected \(\mathbf{k}\) values marked from 1 to 3 (Region (I)) correspond to states spatially delocalized over the entire ribbon width. Conversely, the states become more and more localized at the edge as we move from Point 4 to Point 11 (Point 6 is located at the top of the band, and Point 11 marks the edge of the Brillouin zone), in Region (II). Clearly, the electronic transport for energies in Region (II) is supported by two states with distinct spatial localization: one state is localized on the ribbon's edge (between \(\mathbf{k}\) Points 6 and 11), while the other is delocalized over the width of the system. Note also that the relative difference in the degree of localization is higher for energies on the left-hand side of the two-channel plateau (energies close to that from Point 11) than for energies on the right (close to the energy corresponding to Point 6).