## 2 Systems and Methods  
The systems under scrutiny consist of AFM heterospin junctions made up of a single defect-free zGNR (Figure 1a). To design such junctions, we connect two semi-infinite pristine zGNRs prepared in their AFM ground states, with the essential difference between them being that the degenerate states are chosen to correspond to an AFM and inverted AFM (\(\overline{\mathrm{AFM}}\)) distributions far from the connecting region. The spin distribution is then obtained self-consistently in the rest of the system. This constraint, applied far from the interface, leads to the development of an effective AFM\(|\overline{\mathrm{AFM}}\) heterospin junction. Examples of an AFM-homojunction and the corresponding heterojunction are depicted in Figure 1b,c. We will consider a series of heterojunctions composed of different zGNR-\(n\), where the number \(n\) of zigzag strips across the ribbon width takes on a set of values ranging from 4 to 40. The electronic properties are computed using a self-consistent \(\pi\)-orbital tight-binding description (TBU) that includes interactions up to third nearest-neighbors and a Hubbard term to account for spin-spin interactions, as detailed in [8]). The numerical values for the tight-binding hopping integrals are \(\gamma_{1}=3.2\) eV, \(\gamma_{2}=0\) and \(\gamma_{3}=0.3\) eV for first-, second- and third-neighbor sites (overlap integrals are taken as zero). The first-neighbor hopping integral is reduced by a \(\Delta\gamma_{1}=0.2\) eV correction on the edge atoms to account for the different chemical environment at those sites [28]. The Hubbard on-site parameter is given by \(U=2.944\) eV, as determined by fitting with density functional theory calculations [8]. The converged Hamiltonian matrices are recast into an open-system setup appropriate for a Landauer-Green's function framework to compute the quantum conductance in extended systems [29]. [PERSON]'s function (GF) for the central scattering region is obtained by the application of a set of efficient recursive algorithms recently proposed [30, 31]. To validate results obtained with TBU, spin-polarized density functional theory (DFT) calculations were performed using the VASP package [32, 33] within the generalized gradient approximation (GGA) using the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [34]. In DFT, we performed non-collinear spin calculations to make it possible to impose proper boundary conditions corresponding to the heterospin junction. Comparison with TBU was performed for the system described by a supercell composed of 18 unit cells of a zGNR-6 (Figure 1d). We constrained the spin polarization on the edges for the first three and the other three intermediate unit cells to an AFM and an \(\overline{\mathrm{AFM}}\) configuration, respectively. For DFT, only the direction of magnetic moments was constrained while their magnitude was allowed to relax. While imposing the spin constraints in those cells, the spin distribution on the other cells was unconstrained and obtained in a self-consistent manner. The resulting graph in Figure 1d shows atomically-resolved magnetic moments along one edge for both DFT and TBU and indicates that all of the features of the spin distribution are in quantitative agreement. Note that, for the TBU-based transport calculations, we performed a self-consistent calculation to obtain the converged Hamiltonians using a larger number of unit cells in the central region while constraining the spin to the bulk AFM distribution for the cells corresponding to the leads, to ensure proper convergence with the system size. Figure 1: (**a**) Atomic structure of a periodic zigzag-edged graphene nanoribbon (zGNR)-\(n\), where \(n\) is the number of zigzag strips along its width (\(n=6\) in this example). \(M=98\) unit cells are used as the spin-relaxation region sandwiched between two semi-infinite leads with spin distributions matching those of the bulk AFM state (the leads unit cells are highlighted in red).