Spin-polarizing electron beam splitter from crossed graphene nanoribbons

Junctions composed of two crossed graphene nanoribbons (GNRs) have been theoretically proposed as electron beam splitters where incoming electron waves in one GNR can be split coherently into propagating waves in \emph{two} outgoing terminals with nearly equal amplitude and zero back-scattering. Here we scrutinize this effect for devices composed of narrow zigzag GNRs taking explicitly into account the role of Coulomb repulsion that leads to spin-polarized edge states within mean-field theory. We show that the beam-splitting effect survives the opening of the well-known correlation gap and, more strikingly, that a \emph{spin-dependent} scattering potential emerges which spin-polarizes the transmitted electrons in the two outputs. A near-perfect polarization can be achieved by joining several junctions in series. Our findings suggest that GNRs are interesting building blocks in spintronics and quantum technologies with applications for interferometry and entanglement.

Graphene is an exceptional material with attractive properties to explore fundamental physics and for use in technological applications [1].While ideal graphene is non-magnetic, custom-shaped graphene nanostructures can be designed to exhibit complex magnetic phenomenology with promising possibilities for a new generation of nanoscale spintronics devices [2,3].In fact, graphene π-magnetism is more delocalized and isotropic than conventional magnetism arising from d or f orbitals, which makes it electrically accessible [4] and stable even at room temperature [5].The intrinsically weak spinorbit and hyperfine couplings in graphene lead to long spin coherence and relaxation times [6] as well as a long spin-diffusion length that is expected to reach ∼ 10µm even at room temperature [7].This makes graphene an interesting platform for designing functionalities such as spin filters [8][9][10][11], spin qubits [12,13] and electron quantum optics setups [14].
Graphene nanoribbons (GNRs) have emerged as particularly attractive building blocks for molecular-scale electronic devices because they inherit some of the exceptional properties from graphene while having tunable electronic properties, such as the band gap dependency on their width and edge topology [8].With the advent of bottom-up fabrication techniques, long defect-free samples of narrow GNRs can now be chemically produced via on-surface synthesis as demonstrated in the seminal works for armchair (AGNR) [15] and zigzag (ZGNR) [16] ribbons.Furthermore, manipulation of GNRs with scanning tunneling probes [17,18] opens the possibility to build two-dimensional multi-terminal graphene-based electronic circuits [19], where their spin properties can be addressed by using spin polarized tips [20] and probed by shot noise measurements [21].
Indeed, electron transport in GNR networks has been theoretically explored with the Landauer-Büttiker formalism [22] for a rich variety of multi-terminal device configurations [23][24][25][26].Most recently, crossed GNR junctions have been proposed as electron beam splitters for electron quantum optics [27][28][29].In these works it was found that by placing one GNR on top of another with a relative angle of 60 • the electron transfer process between the ribbons is strongly enhanced.This enables to split incoming low-energy electron waves between two outgoing ports with a tunable ratio and negligible reflection probability, an effect with roots in valley (chirality) preservation in the low-energy bands of ZGNRs [30,31].However, since ZGNRs develop spin-polarized edge states, as theoretically [32] and experimentally [5,33] demonstrated, one may expect that Coulomb repulsion could give rise to additional interesting features for the charge and spin transport in crossed ZGNRs.For instance, it has been shown that the introduction of one rough zigzag edge can be used to boost spin injection [34].
In this manuscript we analyze the electronic structure and quantum transport properties of junctions composed of two infinite ZGNRs crossed with a relative angle of 60 • using the mean-field Hubbard (MFH) model in combination with nonequilibrium Green's functions (NEGF) to describe the open quantum systems [35].We show how the Coulomb repulsion opens a transport band gap and generates a spin-dependent scattering potential in the junction, which enables the devices to be operated as a spin-polarizing beam splitter.
For a transparent analysis and efficient numerics we use the Hubbard Hamiltonian [36] within the mean field ap- proximation, well suited to describe sp 2 carbon systems [2], for both semi-infinite electrodes and device region as shown in Fig. 1, i.e., Here c iσ is the annihilation operator of an electron at site i with spin σ = {↑, ↓} and n iσ = c † iσ c iσ the corresponding number operator.The matrix element t ij is computed by a two-center integral based on a Slater-Koster parametrization as explained in Ref. [29] and U accounts for the Coulomb interaction between two electrons occupying the same p z orbital.We fix U = 3 eV which is in the typical range that yields a good agreement with ab initio calculations [2,9,11,37,38].The open system described by Eq. ( 1) is solved self-consistently using the NEGF method [35,39,40] as detailed in the supplemental material (SM) [41].The corresponding many-electron state thus takes the form of a single Slater determinant of the occupied single-particle states from the MFH-NEGF equations.
Figure 1(a) shows the device structure for two ABstacked ZGNRs, each with a width of 8 carbon atoms (8-ZGNRs).In principle, away from the crossing (but within the spin correlation length), each of the four electrodes can be imposed one of the two possible symmetrybroken spin configurations at the edges, leading to 2 4 /2 = 8 unique boundary conditions for the device region.The self-consistent solutions to this problem are shown in Figs.S6 and S7 [41] for AB-and AA-stacked junctions, respectively, along with the electronic energy differences.The spin polarization for the two lowest-energy states with AB-stacking are shown in Fig. 1(a-b).In the following we label these as ↑↓ and ↑↑ , where the first (second) arrow refers to the spin orientation of the lower edge of the horizontal (inclined) GNR.Although the electronic energy of ↑↑ is found to be 82 meV above that of ↑↓ with AB-stacking, it is interesting to consider both configurations as this (constant) energy penalty may be compensated by a (length-dependent) energy preference for a certain polarization on the extended GNRs through interactions with their environment.
The spin-and energy-resolved transmission probability between any pair of electrodes can be computed from , where G σ is the device Green's function and Γ ασ = i(Σ ασ − Σ † ασ ) the broadening matrix related to the self-energy Σ ασ from electrode α and for spin-orientation σ [22,39].Similarly, the siteresolved density of scattering states can be computed as Figures 1(c-d) show the spatial distribution of the scattering states incoming from electrode 1 in the conduction band.At each lattice site the disk size is proportional to the density of states (summed over spin) while its color indicates the local majority spin.The electron energy is chosen at E = 0.5 eV above the Fermi energy E F = 0, i.e., slightly away from the window with edge states.This implies mode propagation involving only a single GNR subband, cf.Figs.S3 and S4 [41], as well as robustness against edge disorder [42].Figures 1(c-d) also illustrate how the transmitted wave-for both spin configurations ↑↓ and ↑↑ -is split into electrodes 2 and 3 with negligible reflection and amplitude in electrode 4, as expected for the beam splitter.Conceptually, this is expressed with the representation in Fig. 1(e-f) along with the computed transmission probabilities.
Remarkably, ↑↓ and ↑↑ differ substantially when one considers the spin-resolved transmissions.Whereas ↑↓ does not polarize the current, since the transmission probabilities for both spin channels are equal, the ↑↑ configuration leads to a ratio T ↓ 12 /T ↑ 12 = 0.4, i.e., a spin filtering effect.
For further quantitative analysis, Fig. 2 reports the spin-and energy-resolved transmission and reflection probabilities for an electron injected from terminal 1 into the ↑↓ (panels a-d) and ↑↑ (panels e-h) configu-rations.For comparison, each panel includes the corresponding results for the unpolarized device (U = 0, dashed gray lines) reported previously [29].The introduction of Coulomb repulsion has two direct consequences: (i) it opens a transport gap near zero energy due to polarization of the edge bands, and (ii) it shifts the states at the Brillouin zone boundary (Figs.S3 and S4 [41]) resulting in the formation of two transverse modes at very low energy.While the beam splitting effect in the two-mode energy range is hampered by substantial scattering and reflection (Fig. 2d,h), it is completely restored in the energy range with only a single mode, i.e., 0.4 eV < |E| < 1.3 eV, a condition already identified for unpolarized devices [29].In fact, the transmission properties for ↑↓ coincides there with those of the unpolarized device (Fig. 2a-d).On the other hand, for the ↑↑ configuration the probabilities T 12 and T 13 show a strong spin splitting (Fig. 2e-h), revealing that the spin-filtering effect emphasized in Fig. 1(d,f) exists for the whole band.
This qualitative difference between ↑↓ and ↑↑ can be understood by considering the different symmetries that apply to these two configurations.Geometrically, the considered AB-stacked structure possesses one mirror-symmetry plane as shown by the dashed lines in Fig. 1(a,b) [29].The difference emerges when one considers symmetry-lowering by the spin polarization: For ↑↓ the spin index maps into the opposite through the mirror operation (red axis) while for ↑↑ the spin index is conserved.More specifically for ↑↓ , these spatial symmetries impose constraints in the transmission probabilities between the spin channels, e.g., that T σ 12 = T σ 43 , T σ 13 = T σ 42 , etc. Further, considering probability conservation for injection from electrodes 1 or 2 one has the relations ). Together with time-reversal symmetry (T σ ij = T σ ji ) it follows that T σ 12 = T σ 12 in the case of ↑↓ , i.e., that the transmissions are spin independent.
For ↑↑ no such condition applies and the spin channels are decoupled and the transmission probabilities may be very different.Indeed, this is directly seen in our calculations.
If we consider junction imperfections the aforementioned symmetry constraint would be absent and the spin-polarizing effect no longer symmetry forbidden.To examine the relationship between geometry and transport properties we use as a measure the spin polarization in the transmission between a pair of electrodes: Figure 3 shows P 12 at E = 0.5 eV as a function of inplane translations of one ribbon with respect to the other for both ↑↑ and ↑↓ configurations.The AB-and AAstacked geometries are indicated with symbols in the density plots.Evidently, away from these high-symmetry situations the spin-polarizing effect is generally present.The same conclusion holds true also for a range of twist angles (Sec.S11 [41]).
At this point it should be noted that it may be difficult to prepare the device in one specific spin configuration, such as the low-energy states ↑↓ and ↑↑ discussed up to now.For instance, it is not possible to tune which one is the energetically lower (and thus at low temperatures thermally stable) state by a homogeneous magnetic field as the Zeeman energy is the same for both solutions.On the other hand, transverse electric fields across the individual electrodes [8] or injection of spinpolarized currents at the edges from the tip of an STM [43] could potentially be strategies to control their magnetization.Nevertheless, our fundamental assumption is that the different collective spin states of the device are sufficiently long lived and robust to be probed by a transient current pulse.This assumption is supported by the fact that our calculations predict that the electronic energy is increased by about 0.20 eV when a magnetic domain wall is inserted into a 8-ZGNR (Figs.S6 and S7 [41]), an indication of a very large barrier even compared to room temperature.
The spin-polarizing effect of a single junction discussed above can be enhanced by placing several consecutive crossings to form an array of scatterers as displayed in Fig. 4(a).Because back-scattering is negligible in the single-mode energy region, we can approximate the overall transmission probability across an array of N crossings as where T is the transmission of the ith junction.This approximation was tested for the case of N = 3 and shows an excellent agreement compared to a calculation of the full device (Sec.S9 [41]).This idea is exemplified in Fig. 4(b) for two different scenarios: an ideal arrangement of identical ↑↑ AB-stacked configurations (blue circles) as well as a more realistic situation corresponding to random sampling (orange squares, Sec.S12 [41]) over different spin-, intersection angle (within 55-65 • ), and translation configurations.This shows that with four crossings the total current polarization can reach P 12 ∼ 95% with a transmission of T ↑ 12 ∼ 32% in the ideal case.Even in the pessimistic case with random junctions, where partial cancellation can occur due to sign changes in the individual P 12 , the spin polarization |P 12 | of the array approaches 1 exponentially in √ N (black curve, see [41]).The best 1st percentile (top gray curve) of the sampled arrays still reaches P 12 ∼ 80% for N = 8.Although this statistical analysis-described in detail in Sec.S12 [41]-is based on the simplifying assumption of equal weights of the configurations, it serves to illustrate that arrays can be interesting even if one does not have precise control over the individual junctions.
In conclusion, we have analyzed the spin-dependent transport properties of crossed ZGNRs using MFH and NEGF theory, and found that the beam-splitting effect reported previously survives in presence of Coulomb repulsions with two distinct modifications: a transport gap opens at low energies and a spin-dependent scattering potential emerges.Except for specific high-symmetry configurations, this class of electronic devices are generally predicted to behave as spin-polarizing beam splitters with interesting possibilities for electron quantum optics [44].Such spin-dependent scattering potentials are also obtained with other edge-polarized nanoribbons (Sec.S13 [41]).By constructing arrays of junctions the spin-polarizing effect can be enhanced.
Although the proposed devices are ahead of current experiments, a rapid progress in bottom-up fabrication and scanning probe techniques makes it conceivable to assemble nearly defect-free junctions on insulating thin films [45], to drive coherent electron dynamics [46,47], and to characterize electron transport by multi-probe setups [48] or through single-photon emission [49].Our results add to the vision of using GNR-based devices for spintronics and quantum technologies.For instance, two spinpolarizing beam splitters in combination with a charge detector can be used to deterministically entangle a moving spin qubit [50].Conversely, a spin-polarizing beam splitter can also be used to determine the entanglement of injected pairs of spins [51].As an additional application, a high-fidelity spin filter allows "spin-to-charge" conversion and thus a charge-measurement-based spin determination.

SUMMARY OF CONTENT
In this supplementary material we describe our methodology and present additional calculations that may be interesting for a deeper understanding of the reported effects.In Sec.S1 we explain the details of the calculation with the mean-field Hubbard model (MFH).
In Sec.S2 we show the effect of the Coulomb repulsion on the bands of bilayer and monolayer graphene nanoribbons, while in Sec.S3 we show the effect that the size of the scattering region has on the local magnetization of the device.Next in Sec.S4 we plot all the possible spin configurations resulting from the possible combinations of the spin densities of the four electrodes.In Secs.S6 and S7 we show the figure of merit, where we analyze the quality of the device as a spin-beam splitter or mirror, and the spin polarization at a different electronic energy from the one shown in the main text, respectively.In Sec.S8 we compare the transport properties for different ribbon widths.In Sec.S9 we show the quality of the independent scatterers approximation compared to the exact result for a device with three crossings.In Sec.S10 we compare the spin-averaged transmission probabilities with the unpolarized case.In Sec.S11 we explore how possible distortions, such as a small twist angle or a lateral translation of the on-top ribbon with respect to the bottom one, affect the transport properties of individual crossings.In Sec.S12 we describe the used statistical sampling and analytical expression for the averaged spin-polarization shown in Fig. 4(b) of the main text.Finally, in Sec.S13 we analyze the spin-polarizing beam-splitting effect that can be found in junctions formed with other edge-polarized GNRs, such as those built with crossed bearded GNRs.* sofia.sanz@dipc.org† thomas frederiksen@ehu.eus

S1. SOLVING THE MFH WITH OPEN BOUNDARY CONDITIONS
The Hubbard Hamiltonian in the mean-field approximation model has proven to be in remarkable agreement with quantum Monte Carlo simulations for ZGNR for moderate Coulomb interactions [1].The complete geometry is divided into the relevant parts here involved: the semi-infinite electrodes and the scattering area (device), where the latter contains the crossing between the two infinite ribbons, , where H d is the device Hamiltonian and H α , H αd are the αth electrode Hamiltonian and its coupling to the device region.The occupation of the electronic states is defined by Fermi-Dirac statistics with a temperature set to T = 300 K.For our system the band gap opening is of the order 500 meV, i.e., ∼ 25 × kT at this temperature (and k the Boltzmann constant), which implies that the distribution function is basically a step function.Therefore, practically there is no difference with the results one would obtain at lower temperatures.We also assume an orthogonal basis of localized atomic orbitals.Further, the qualitative picture presented in the main text is not affected by the numerical choice for the Coulomb repulsion parameter U = 3 eV, but only the quantitative results.To show this last statement we have calculated the transmission probabilities for the ↑↑ device obtained with U = 1.5 eV in Fig. S1.The computational setup begins by solving the electrodes, where the spin densities n i,σ are found by diagonalization of the infinite ZGNR's unit-cell at each k-point over the Brillouin zone.The spin-densities will be used to update the Hamiltonian of Eq. ( 1) in each iteration, until the convergence criterion is achieved, as implemented in our Python package Hubbard [2].To properly account for the effect of the semi-infinite leads in the device Hamiltonian, the spin densities in equilibrium are computed by an integration of the Green's function along a predefined energy contour in the complex plane [3] that we obtained from TranSiesta [4], where n F (ε − µ) is the Fermi distribution with µ the electrochemical potential, For equilibrium calculations the electronic contribution to the total energy can be calculated as where the left term of Eq. (S2) is the integration of the occupied states while the right term comes from the interaction term of the Hamiltonian.To perform the numerical transport calculations we have used the free and open-source code TBtrans [4].

ZGNRS
In this section we show the band structures for periodic mono-and bilayer w-ZGNRs for widths of w = 8, 16 carbon atoms across.The geometries of the AA-and AB-stacked ZGNRs are shown in Fig. S3(a,b).In Fig. S3(c-h) we show the band structures for the unpolarized (black dashed lines) and polarized (color lines) Hamiltonians with U = 3 eV.The main effect of the interaction term on the electronic structures of these systems is the opening of the correlation gap around E F .While the larger hybridization in the AA-stacking pushes the edge states further in energy for the unpolarized case (panels (e,f)), which competes with the Coulomb repulsion parameter, in the AB-stacking case the presence of the flat bands (edge states) give rise to the opening of a correlation gap around E F .The spin distribution of the bilayer ZGNRs (of lowest energy) is composed of two polarized monolayer ZGNRs (antiferromagnetic alignment between the edges of the ribbons) while the atoms that are vertically aligned are also antiferromagnetically aligned.These results are in line with DFT calculations [7].
In Fig. S4 we plot the polarization represented by the center of mass of the wavefunctions ψ ↑ nk , calculated as y|ψ ↑ nk | 2 dy, for the case of the monolayer 8-ZGNR.For an unpolarized wave the center of mass of the wave coincides with the geometrical center of the unit cell, given the inversion symmetry of the ZGNR unit cell.However, for the polarized ZGNR, this is not necessarily the case given the symmetry breaking between the two spin components.
We observe a large polarization of the wave especially in the valence and conduction bands.We also find an interesting behavior of the spin-wave distribution, where the center of mass of the wave transitions from the lower half of the unit cell to the upper half as k goes from Γ to X.This explains the spin majority on the different sublattices depending on the electron energy.

S3. DISTANCE FROM THE SCATTERING CENTER TO THE ELECTRODES
Here we show results for the crossed 8-ZGNRs with different GNR lengths, i.e., different distances from the scattering center to the leads.We performed calculations for the same crossing with the same spin configuration ( ↑↑ ) and compare the spin densities to that corresponding to the perfect (periodic) ribbon solution, i.e., uncoupled ribbons.To do so we plot in real space the difference between the local magnetization, defined as for the self-consistent solution of the full device calculation and the local magnetization corresponding to the periodic ribbons, m 0 i .In order to preserve the bulk nature for the electrodes, there must be a sufficient distance separating the scattering center from the leads.However, the highly localized Coulomb repulsion term seems to prevent the spin density of the device from being affected by this size effect as Fig. S5 shows that the spin polarization of the coupled ribbons is essentially unaffected by the position of the electrodes.In other words, the inter-GNR coupling induces changes in the electronic structure only in the near vicinity to the crossing.

S4. ALL INEQUIVALENT SPIN CONFIGURATIONS AND TRANSMISSION CURVES FOR AB-AND AA-STACKINGS
Here we show all the unique spin configurations for AA-and AB-stacked 8-ZGNRs and the corresponding transmission and reflection probabilities for each spin component.We also list the electronic energy according to Eq. (S2) of each configuration.As mentioned in the main text, each spin configuration for this open quantum system is found by fixing spin density distribution of the 4 electrodes, leading to 2 4 /2 = 8 unique solutions to the imposed boundary conditions (excluding the trivial global inversion of the spin).To show these solutions we plot the difference between the spin densities for the ↑ and ↓ spin components, i.e., n ↑ − n ↓ , in Figs.S6 and S7 for the AB-and AA-stacking, respectively.
In Tab.S1 we compare the relative electronic energies for the AA-and AB-stacked devices with the 8 unique spin-configurations (a-h) that arise from imposing the spin-densities in the electrodes.S6) and AA-stacked (Fig. S7) crossed 8-ZGNRs devices.The energies (in eV) are compared to the configuration of minimum energy corresponding to AB-stacked ↑↓ in Fig. S6(a).
The first observation is that we find the AA-stacked devices to have larger electronic energy than the AB-stacked ones for all spin configurations (see Tab. S1).This result is in line with other results of this kind [8].The second observation is that the ground state is the configuration with antiferromagnetic (AFM) alignment between layers, i.e., the atoms that lie one on top of the other have opposite spin index, in line with results published in Ref. [9].On the other hand, some of the spin configurations involve the presence of domain walls (grain boundaries), if the spin densities of the electrodes belonging to the same ribbon are inverted.
We observe that for configurations that have a single grain boundary, the domain wall moves along the ribbons to leave the scattering area (crossing) with the spin density distri-bution that corresponds to the ground state configuration [cf.Figs.S6 and S7].For instance, for the AB-stacked device, the grain boundary in Fig. S6(e-h) moves along the ribbons to leave the crossing with the spin density distribution of Fig. S6(a).Similarly, for the AAstacked devices, the grain boundary moves along the ribbons to leave the crossing with the spin density distribution of Fig. S7(b).We find that, for the two cases without grain boundaries, the spin configuration of each ribbon deviates very little from that corresponding to the perfect (translationally invariant) 8-ZGNR, showing the weak coupling between them.
On the other hand, Fig. S8 and Fig. S9 shows the transmission and reflection probabilities for the AB-and AA-stacked devices with the different possible spin distributions.As it can be seen in these figures, and following the argument and discussion from the main text, the only systems that have T σ ij = T σ ij are those that only display black symmetry axes (break the symmetry between the spin indices).The transmission probability is computed as mentioned in the main text, while the reflection probability can be computed by subtracting the total outgoing transmission probability to the number of channels/modes existing at a certain energy, M α , We observe that both R s 1 and T s 14 are zero for energies lying in the single-band energy region for both stackings and all spin configurations, according to the unpolarized case [10].

S5. ELECTRONIC TOTAL ENERGY AND MAGNETIZATION
Fig. S10 shows the electronic part of the total energy of each device compared to the energy of the uncoupled system, E 0 , as a function of the translation of the on top ribbon with respect to the bottom one starting with the geometry of the AB-stacked device.This energy can be understood as a binding energy between the ribbons.We note that this energy lacks some contributions that are not taken into account in this approximation, e.g., the change in the Van der Waals forces as the two ribbons are translated with respect to each other which determines the precise inter-GNR separation, etc.Here, the geometrical distortion is only encoded through the Slater-Koster parametrization [10].
As mentioned above, the closest stackings to the AB pattern lie in a global minimum.The AA-stacking is a local minimum but more energetic than the AB-stacking.One interesting result is that the plots Fig. S10(a,b) are very similar, showing that the largest contribution to the electronic energy comes from the geometry and that the relative spin density distribution of the ribbons plays a minor role in this physical quantity.
Fig. S10(c,d) show the sum of local magnetization changes induced by the inter-GNR interaction, defined as i |m i | − |m 0 i |.This shows that the magnetization of the device is always lower than that corresponding to the uncoupled ribbons.This last statement makes sense since the effect of the hopping amplitude between the GNRs goes "against" the localization of the electrons, and thus the local magnetization of the system.Therefore the local magnetization of the coupled ribbons will be lower than the magnetization of the uncoupled ribbons, especially in the coupled area.Fig. S10(c,d) also shows that for the configurations with FM inter-layer coupling (AB in Fig. S6b and AA in Fig. S7a, i.e., the atoms that are vertically aligned have equal spin indices), the magnetization decreases more with respect to the perfect 8-ZGNR than configurations with AFM inter-layer coupling (AB in Fig. S6a and AA in Fig. S7b, i.e., the atoms that are vertically aligned have opposite spin indices).
In panels Fig. S10(e,f) we show the maximum backscattering for configurations ↑↓ and ↑↑ as a function of the translation of the on-top ribbon with respect to the other one.We hereby see that the low reflection is general for the crossed ZGNRs.

S8. ROLE OF RIBBON WIDTH
In this section we compute the spin polarization distribution and transmission and reflection probabilities as a function of the electron energy for different ribbon widths for the same crossing and spin configuration (AB-↑↑ ).We see that the different transport behavior for the two different spin channels is general for this crossing.On the other hand, we see that the inter-transmission probability (T s 13 ) grows with the ribbon width, while the opposite behavior is found for the intra-transmission probability (T s 12 ).Furthermore, losses (R s 1 +T s 14 ) remain absent independently of the width of the ribbon.These results are in line with the transport properties found for the unpolarized case [cf.Ref [10]].

S10. AVERAGED TRANSMISSION PROBABILITIES
In this section we compute the averaged transmission probabilities, calculated as Tij = (T ↑ ij + T ↓ ij )/2, for both spin configurations ↑↑ and ↑↓ compared to the unpolarized case in Fig. S16.By comparing these results and the ones shown in the main text, it can be seen that while the spin-averaged transmission of the device reproduces our earlier results [10] (i.e., it is affected very little by inclusion of the mean-field Coulomb interaction), we find a strong spin-dependence of the transmission, something that is completely absent in the non-interacting case.

S13. A SPIN-POLARIZING BEAM SPLITTER WITH BEARDED GNRS
In this section we analyze the transport properties of intersecting bearded GNRs.The band structure and geometry of this structure is shown in Fig. S21.Although, to our knowledge, these systems have not yet been synthesized, we discuss these devices to complete our analysis on spin polarizing beam splitters formed by two crossed general GNRs.This type of GNRs are oriented along the zigzag direction but display different physical edges than ZGNRs.However, they present similar spin-polarized edge states as seen in Fig. S21.
In Fig. S22 we plot both the self-consistent solutions ↑↑ and ↑↓ and the scattering states for this device, while in Fig. S23 we show the more detailed transmission probabilities for an incoming electron into terminal 1 as a function of the electronic energy.

FIG. 1 .
FIG.1.Transport setup and spin-dependent properties for AB-stacked 8-ZGNRs devices.(a,b) Two different self-consistent solutions for the spin-density distribution in the device region, labeled ↑↓ and ↑↑ , respectively, defined by the spin orientation of the lower edge of each GNR.The up (down) spin density is shown in red (blue).The lower, horizontal ribbon is plotted in black, while the upper, intersecting at an angle of 60 • , is depicted in gray.Electrodes 1-4 are indicated.The ribbons are separated by a distance d = 3.34 Å along the z-axis, as displayed in the side view (lower part of b).The dashed lines in each configuration indicate a symmetry axis that maps the device geometry to itself through mirror operations, where the red (black) color of the axis further indicates that the spin index is inverted (conserved) by the symmetry operation.(c,d) Spin-resolved density of states of scattering states incoming from electrode 1 for the ↑↓ and ↑↑ spin configuration, respectively, computed at E − EF = 0.5 eV.The dominant spin on each site at this energy is shown in red for up-spins and in blue for down-spins.(e,f) Sketch of incoming and outgoing waves through the scattering center (represented by the circled cross) and the corresponding transmission probabilities from calculations.

14 − 1 1 − 1 FIG. 2 .
FIG.2.Spin-and energy-resolved transmission probabilities T12, T13, T14 and reflection R1 for (a-d) the ↑↓ and (e-h) ↑↑ configurations of Fig.1.Electrons are injected from electrode 1.The red (blue) curves correspond to the up (down) spin components with U = 3 eV.For comparison, the corresponding calculations for the unpolarized case is indicated by dashed gray lines.

FIG. 3 .
FIG. 3. Spin polarization P12 of the current from electrode 1 to 2 as function of in-plane translations of one ribbon with respect to the other for (a) ↑↑ and (b) ↑↓ configurations introduced in Fig. 1.The electron energy is in the conduction band at E = 0.5 eV.The in-plane unit cell (dashed lines) has lattice vectors a1 and a2, where a0 = 2.46 Å is the graphene lattice constant.The red crosses (green disks) indicate the high-symmetry configurations with AB (AA) stacking.

FIG. 4 .
FIG. 4. (a) Sketch of an array of three consecutive AB-stacked 8-ZGNRs crossings to enhance the spin-polarized current at the output electrode 2. (b) Spin polarization P12(N ) (filled symbols) and majority-spin transmission T maj 12 (N ) (open symbols) as a function of the number of crossings N between terminals 1 and 2 in the conduction band at E = 0.5 eV.Two different scenarios are considered: an ideal arrangement of identical ↑↑ AB-crossings (blue circles) as well as a random sampling (orange squares) over 10 7 different spin-, intersection angle (within 55-65 • ), and translation configurations drawn from the data in Sec.S11 [41] assuming equal weights.The average polarization |P12(N )| follows an analytic expression (black line, see Sec.S12 [41]) approaching 1 exponentially in √ N .The gray lines indicate the best (1, 5, 10, 25, 50)th percentiles of the random distribution.
FIG. S1.Transmission probability for an incoming electron from terminal 1 of the ↑↑ device obtained with (a) U = 3 eV and (b) U = 1.5 eV.
σ is the retarded Green's function for each spin component σ.Σ α,σ is the self-energy matrix that accounts for the coupling between the αth semi-infinite lead with spin component σ to the scattering region.The self-energy matrices are converged with the Lopez-Sancho recursive method[5]  as implemented in the open source, Python-based SISL package[4, 6]  using a small broadening of η = 1 meV.For clarification, the flow diagram of the self-consistent cycle is plotted in Fig.S2.
FIG. S3.Band structures of mono-and bilayer w-ZGNRs of different widths.(a,b) Geometry of the AA-and AB-stacked bilayer 8-ZGNRs.Band structure along the Γ-X path for (c,d) monolayer and (e,f) bilayer AA-stacked and (g,h) bilayer AB-stacked 8-ZGNRs.Black dashed lines show the band structures for U = 0 and colored lines those obtained after convergence with U = 3 eV.
FIG. S5.Convergence study with respect to the size of the device region for AB-stacked 8-ZGNRs in the ↑↑ configuration.(a-c) Difference between the calculated spin polarization for the 4-terminal device (m i ) and the spin polarization for decoupled ribbons (m 0 i ) plotted in real space for different device sizes.We compare three cases, (a) one with a repetition of 16 (∼ 38 Å), (b) another with 20 (∼ 48 Å) and the last one with (c) 28 (∼ 68 Å) ZGNR unit cells.The latter corresponds to the size of the device shown in the main text.(d-f) Profile of m i − m 0 i along the ribbon axis for different positions across the perpendicular direction (confinement direction).The legend on the left side indicates the transverse position of each profile.
FIG. S6.Spin configurations for the AB-stacked 8-ZGNRs.The electronic energy E tot (in eV) of each configuration (relative to AB-stacked ↑↓ in panel a) is noted in the bottom left corner in each panel.Red and black dashed lines indicate the symmetry axes.
FIG. S10.(a,b) Electronic part of the total energy E tot relative to the energy of the uncoupled system E 0 , (c,d) sum of the absolute value of the magnetization i |m i | compared to the uncoupled system i |m 0 i |, and (e,f) maximum backscattering into terminal 1 as function of the translation of the top ribbon with respect to the bottom one for the two possible spin configurations ↑↑ and ↑↓ , respectively.The dashed parallelogram corresponds to the primitive cell of the problem.Red crosses and green circles indicate the AB-and AA-stacking, respectively.
FIG. S13.(a,b) Polarization for an electron incoming from terminal 1 and outgoing in terminal 2 and (c,d) [(e,f)] figure of merit (FM) for σ =↑ [σ =↓] as function of the translation of the top ribbon along the graphene periodicity vectors a 1 , a 2 with respect to the bottom one for the two possible spin configurations ( ↑↑ and ↑↓ , respectively) obtained at E = −0.5 eV.Red crosses and green circles indicate the AB-and AA-stacking, respectively.
FIG. S15.(a) Geometry of the device with three consecutive crossings.All 8 terminals are indicated in red squares.(b) Transmission probabilities between the different pairs of terminals for an incoming electron in terminal 1 as a function of the electron energy.Solid lines represent the transmission probability for the full device, while open circles represent the obtained transmission probability using the independent-scatterers approximation.

8 FIG
FIG. S16.Black solid lines stand for Tij = (T ↑ij + T ↓ ij )/2 while the gray dashed lines stand for the unpolarized transmission probabilities.

FIG
FIG. S21.(a) Band structure for the bearded GNR calculated with U = 0 (dashed black lines) and U = 3 eV (solid red lines).(b) Spin density for the periodic structure (U = 3 eV).
FIG.S23.Spin-and energy-resolved transmission probabilities T 12 , T 13 , T 14 , and reflection R 1 for (a-d) the ↑↓ and (e-h) ↑↑ configurations.Electrons are injected from electrode 1.The red (blue) curves correspond to the up (down) spin components with U = 3 eV.For comparison, the corresponding calculations for the unpolarized case is indicated by dashed gray lines.

TABLE S1
. Electronic energies E tot according to Eq. (S2) for different spin configurations of the AB-(Fig.