Crossed graphene nanoribbons as beam splitters and mirrors for electron quantum optics

We analyze theoretically 4-terminal electronic devices composed of two crossed graphene nanoribbons (GNRs) and show that they can function as beam splitters or mirrors. These features are identified for electrons in the low-energy region where a single valence or conduction band is present. Our modeling is based on $p_z$ orbital tight-binding with Slater--Koster type matrix elements fitted to accurately reproduce the low-energy bands from density functional theory calculations. We analyze systematically all devices that can be constructed with either zigzag or armchair GNRs in AA and AB stackings. From Green's function theory the elastic electron transport properties are quantified as a function of the ribbon width. We find that devices composed of relatively narrow zigzag GNRs and AA-stacked armchair GNRs are the most interesting candidates to realize electron beam splitters with a close to 50-50 ratio in the two outgoing terminals. Structures with wider ribbons instead provide electron mirrors, where the electron wave is mostly transferred into the outgoing terminal of the other ribbon, or electron filters where the scattering depends sensitively on the wavelength of the propagating electron. We also test the robustness of these transport properties against variations in intersection angle, stacking pattern, lattice deformation (uniaxial strain), inter-GNR separation, and electrostatic potential differences between the layers. These generic features show that GNRs are interesting basic components to construct electronic quantum optical setups.


I. INTRODUCTION
The similarities between the wave nature of electrons propagating coherently in ballistic conductors with photon propagation in optical wave guides has spawned the field of electron quantum optics [1,2]. In this way several electronic analogues of optical setups-such as the Mach-Zehnder [3,4] and Fabry-Pérot [5][6][7] interferometers, as well as the Hanbury Brown-Twiss [8][9][10][11] geometry to study the Fermion antibunching and the two-particle Aharonov-Bohm [12] effects-have been implemented. Fundamental components for these setups include mirrors (M), beam splitters (BS, i.e., partially transparent mirrors), and wavelength filters. Such control elements for electron beams are important in the fields of quantum information and solid-state quantum computation: By sending a single electron through a BS one can generate a mode-entangled state that can be used to violate a Bell inequality [13] or for quantum teleportation [14,15]. A BS is the central building block of the Hong-Ou-Mandel setup to test the indistinguishability [16] or the entanglement [17] of electrons entering in the two input ports. With two BSs and two oriented Ms the Mach-Zehnder interferometer can be fully implemented, which has been demonstrated to work as a quantum logic processor [18].
A platform with remarkable prospects for electron quantum optics are graphene-based systems, in which several pioneering experiments on electron beam splitters * sofia.sanz@dipc.org † thomas frederiksen@ehu.eus and related devices have been performed [19,20]. More recently graphene nanoribbons (GNRs) [21,22] have emerged as attractive candidates for the construction of molecular-scale electronic devices [23] because they inherit some of the exceptional properties from graphene while having tunable electronic properties, such as the opening of a band gap depending on their width and edge topology [24][25][26][27][28]. The electron coherence length in GNRs can be long, with values of the order ∼ 100 nm being reported [29][30][31]. Furthermore, ballistic transport can be rather insensitive to edge defects because of the presence of localized edge states (e.g., in ZGNRs) and the dominating Dirac-like physics, that make the current flow maximally through the center of the ribbon [32]. With the advent of bottom-up fabrication techniques, long defect-free samples can be chemically synthesized with both armchair (AGNR) [33] and zigzag (ZGNR) [34] edge topologies via on-surface synthesis. Manipulation of GNRs with scanning tunneling probes has been also demonstrated [35,36], opening the possibility to build two-dimensional multi-terminal graphenebased electronic circuits [37][38][39][40][41].
Recently, it has been shown theoretically that two crossed GNRs with a relative angle of 60 • can behave as a BS [42,43] for valence-and conduction-band electrons, since such four-terminal devices were found to divide the electron beam into two out of the four arms with an equal transmission probability of 50%. In this paper we analyze this possibility more generally and show that all the mentioned beam-control elements (BS, M, filters) can be realized with a suitable choice of two crossed GNRs. More specifically, we compute the electron transport properties of these devices in terms of the edge topology and width of the GNRs, and the precise alignment and stacking of the ribbons.
The problem is theoretically approached by means of tight-binding (TB) modeling, which is known to reproduce graphite-like systems with sufficient accuracy [44][45][46][47][48], while allowing to explore a large number of systems of considerable sizes in a fast and transparent way. For instance, the geometry of a crossing between two 50-atom wide GNRs readily comprises around 10 000 atoms. The main complexity of the modeling lies in the description of the interlayer couplings, for which we use a Slater-Koster parametrization [49] that has proven successful for describing the band structure and velocity renormalization of Dirac electrons in twisted bilayer graphene [50,51]. The employed technique can describe arbitrary device geometries and therefore allows us to also study of the robustness of the predicted transport properties against variations in intersection angle, stacking pattern, lattice deformation (uniaxial strain), inter-GNR separation, and electrostatic potential differences between the layers. With this, we give a complete analysis of the transport properties of crossed GNRs, highlight their tunability, and provide quantitative data that can serve as a guide for design optimization.
This paper is organized as follows: in Sec. II we introduce the general TB Hamiltonian used to describe the kinetics of the electrons travelling through the different devices as well as the scattering formalism used to compute transmission and reflection probabilities of incoming electron waves from the different leads. In Sec. III we present a complete analysis of the transport properties based on the key combinations of stacking pattern, edge topology and width of the GNRs. Finally, the conclusions and remarks are provided in Sec. IV.

II. METHODOLOGY
The general setup of this study, illustrated in Fig. 1a, is comprised of two infinite GNRs crossed with a relative angle θ = 60 • . The scattering region (intersection) breaks the translational invariance of the infinite ribbons, for which we will use the Green's function formalism to solve the Schrödinger equation with open boundary conditions.
The system is divided into the device (scattering) region that contains the intersection between the two ribbons, and the four semi-infinite GNRs (periodic electrodes), represented as red rectangles in Fig. 1a. The total Hamiltonian is correspondingly split into the different parts where H d is the device Hamiltonian, H α the α-electrode Hamiltonian, and H αd the coupling between these two subsystems.
(a) Definition of the width W of (b) ZGNRs and (c) AGNR in terms of the number of carbon atoms N across the ribbon. The interatomic distance is denoted by a.

A. Tight-binding Hamiltonian
The use of a local basis in combination with Green's function techniques provides an efficient way for obtaining the transport properties in terms of microscopic parameters. We write the single-particle TB Hamiltonian in an orthogonal basis as where c † i (c i ) creates (annihilates) an electron on site i with energy i . We further define the Fermi level as E F = i , corresponding to half-filled carbon p z orbitals. The matrix element t ij between orbitals i and j is described by Slater-Koster type two-centre π and σ bond integrals between two p z atomic orbitals [49] where l is the cosine of the angle formed between the distance vectorr ij for the ij atom pair and the unit vector that defines the z-direction (cf. Fig. 1a), i.e., l =r ij ·ê z /|r ij |. The two-centre integrals are expressed as where t (t ⊥ ) is the intra-GNR (inter -GNR) hopping parameter between atoms separated by the interatomic (interlayer) distance fixed to a = 1.42Å (d = 3.34Å) in our model, see Fig. 1. The decay rates of the bond integrals with the atomic separation, q σ and q π , are isotropic and therefore related by q σ /d = q π /a. This model, characterized by t , t ⊥ and the decay rate (which can be determined by fixing the second nearest neighbors coupling), successfully describes π electrons in twisted bilayer graphene [51]. However, it does not capture manybody effects like, e.g., the difference in nearest-neighbor hopping parameter for different lattice sites as in the Slonczewski-Weiss-McClure (SWM) model for graphite [45,[52][53][54].
In this work we use t = 2.682 eV and t ⊥ = 0.371 eV. For the third model parameter we refer to the in-plane next-nearest neighbor matrix element t = 0.0027 eV. These parameters were obtained by fitting to the lowenergy band structure of AB-stacked bilayer graphene simulated with Siesta [55] as explained in Appendix A. The satisfactory agreement between TB and DFT ( Fig. 14) further justifies that, at least for our purposes, many-body effects like in the SWM model can be neglected.

B. Transport calculations
In order to perform transport calculations we use the nonequilibrium Green's function (NEGF) method [56][57][58]. In particular, to obtain the transmission probabilities (T αβ ) between the different pairs of electrodes (α = β), we use the Landauer-Büttiker formula [59], is the broadening matrix, related to the nonhermitian part of the retarded electrode self-energy Σ α , due to the coupling of the αth semiinfinite lead to the scattering center and α, β = 1, . . . , 4, cf. Fig. 1. Further, is the retarded Green's function of the device region and I the identity matrix (orthogonal basis). The dependency on the electron energy E of these key quantities is implicit.
The reflection probability (T αα = R α ) can be conveniently written as a difference between the bulk electrode transmission M α (i.e., the number of open channels/modes in electrode α at a given energy) and the scattered part into the other electrodes ( β T αβ ) as From Eq. (7) we can also obtain the spectral function A α for states coupled to electrode α The diagonal elements A α (i, i)/2π correspond to the local density of states (DOS) at sites i of the scattering states originating from electrode α. Computationally, we constructed the Hamiltonian matrix with the sisl package [60,61] and computed transmission probabilities and spectral DOS with TBTrans [61].

III. RESULTS
In this section we present results for the electron transport properties through an extensive set of four-terminal devices formed of two crossed ribbons. We analyze the role of the precise stacking and alignment of the crossing area for both ZGNR-and AGNR-based devices in all their possible configurations.

A. Possible device configurations
For a systematic analysis we begin by considering all the possible devices that can be built with two crossed AA-or AB-stacked GNRs with a relative angle of 60 • . These are sketched Fig. 2. In the case of crossed ZGNRs there exist two configurations, the AB-stacking [labeled AB, Fig. 2(a)] and the AA-stacking [labeled AA, Fig. 2(b)]. These two geometries have different symmetries, indicated by the reflection planes (dashed lines) in Fig. 2. While AB has only one reflection symmetry, AA has two. Here, and in the following, we refer only to symmetries in the xy-plane. The additional operation of reflection in the z-direction to interchange top and bottom ribbon is physically not important and therefore implicit.
In the case of AGNRs there are two different AA-stacked configurations [labeled AA-1 and AA-2, lower one. Similarly, AB-1 (AB-2) can be obtained from AA-1 by translating the upper (red) ribbon by +ax (−ax) with respect to the lower one. Again, these four generic configurations have different symmetries as indicated in Fig. 2(c-f).
The reflection planes imply that there are operations which leave the scattering potential (created by the intersection of the two ribbons) unchanged. This is, if we apply one or more reflections across the indicated axes, the Hamiltonian of the new device does not change. Consequently, the Green's function and all the transport properties derived from it, will also remain unchanged under some particular electrode permutations.
Let us begin by discussing the properties of these six different configurations with particular examples constructed from 8-ZGNRs and 11-AGNRs. In Fig. 3 we show the spectral DOS of scattering electrons that come in from electrode α = 1 as obtained from Eq. (9) for each configuration at specific energies. In this real-space representation it is easy to see where the scattered electron wave propagates after being injected into the device. The large DOS that appears in the input electrode region does not correspond to the backscattered electrons, but rather to the DOS of the incoming electrons (as we will show later on). This is also illustrated in Fig. S1 [62], where we complement the results shown in Fig. 3 by plotting the bond currents between nearest neighbor atoms, where the arrows indicate the direction of the flowing electrons.
For the ZGNR devices, Fig. 3(a,b) and Fig. S1(a,b) show that an electron injected from α = 1 in both cases only escapes from the scattering center into electrodes β = 2, 3, i.e., terminals 1 and 4 are suppressed. Similarly, for the two AA-stacked AGNR devices Fig. 3(c,d) and Fig. S1(c,d) show that the outgoing terminals β = 1 and β = 3 (β = 4) for AA-1 (AA-2) are suppressed. These two cases are interesting since their relative displacement of only √ 3aŷ leads to very different electron transport: for AA-1 the split electron turns by 60 • , while for AA-2 the bend is 120 • . In the case of the two AB-stacked ribbons, Fig. 3(e,f) show that an electron wave in these devices is scattered qualitatively (yet not quantitatively) similarly and into all outgoing electrodes.

B. Symmetry considerations
Since we deal with 4-terminal devices, the matrix of transmission and reflection probabilities, Eq. (6) and Eq. (8), has the general form However, due to symmetries there are not 16 independent quantities in this matrix. First of all, in absence of a magnetic field, time reversal symmetry forces T αβ = T βα . This reduces the matrix to 10 independent elements, e.g., those without the dark gray background (α > β) in Eq. (10). Secondly, the symmetries indicated in Fig. 2 reduce the number of independent elements of the matrix further. The reflection plane y = sin(−60 • )x maps the electrode labels (1, 2, 3, 4) ↔ (4, 3, 2, 1) with unchanged transmissions, e.g., which allows to consider R 3 , R 4 , T 24 and T 34 as dependent variables [4 of the light gray elements in Eq. (10)]. Similarly, the reflection plane y = sin(30 • )x implies (1, 2, 3, 4) ↔ (3, 4, 1, 2) and R 3 , R 4 , T 23 , and T 34 as possible dependent variables (4 of the light gray elements). The combination of both reflection planes further implies (1, 2, 3, 4) ↔ (2, 1, 4, 3) and R 2 and T 23 as further dependent variables (i.e., all gray elements in this case). In summary, depending on the number of symmetries, the transmission probabilities of any given device will be fully characterized by either 4, 6 or 10 independent matrix elements. Figure 4 shows the full, energy-resolved transmission matrix [Eq. (10)] obtained numerically for devices formed of two crossed ZGNRs in the AB configuration for a range of different ribbon widths W . As ZGNR AB displays only one reflection plane, the transmission probabilities for these systems are, in principle, characterized by 6 independent quantities.
However, qualitatively only 4 independent ones are readily identified in Fig. 4. Only upon close inspection of the data all the expected differences emerge. The reason for the seemingly higher symmetry (corresponding to two reflection planes) is the fact that the scattering potential created by the crossings between the GNRs, depends exponentially on the atomic distances between the GNRs and, therefore, is dominated by the closest atom pairs. These atom pairs, shown in Fig. S2(a) [62], are in fact symmetric with respect to both reflection planes. More generally, for all the configurations in Fig. 2 we find that the scattering potentials are dominated by terms with at least one reflection plane (Fig. S2). For all practical purposes, the effective symmetry appears higher and it suffices to describe the transmission probabilities with only 4 or 6 independent quantities.
In the following we will thus only consider it sufficient to discuss electrons incoming from terminal α = 1. However, for completeness we show the full transmission matrices for all the considered systems in Figs. S3-S15 [62].

C. Beam splitters and mirrors
Looking again at Fig. 4 and focusing on the first row (electron beam injected from terminal α = 1), we observe distinct regimes where the devices would present particular electron quantum optical characteristics. We are especially interested in geometries for which the transmission matrix allows to designate two input and two output electrodes in the sense that any electron sent in through one of the input ports is scattered predominantly into the two output ports with very little reflection or transmission into the other input. For instance, the green areas in the plots show where the device behaves as a BS, since they show that the electron beam is scattered only into two out of the four possible arms with a transmission probability that lies around T 12 ∼ T 13 ∼ 0.5. One can also identify regimes in which the device can work as a M where T 13 ∼ 1. This situation corresponds to the red areas in Fig. 4, since the electron would enter from terminal α = 1 and turn 120 • to go out exclusively into terminal β = 3 with low reflection. The energy-dependence of the transmission functions is very symmetric with respect to the Fermi level, reflecting the approximate particle-hole symmetry characteristic of a half-filled bipartite lattice. Nevertheless, the presence of next-nearest neighbor couplings in our TB model in principle breaks this symmetry.
On one hand, we note that for energies close to the Fermi level (|E − E F | < 0.07 eV) in Fig. 4, the electron is scattered into all the four output ports, which makes this small energy window not very interesting for electron quantum optical purposes. These features probably arise due to the hybridization of states from the flat bands of the individual ribbons in the overlapping area. The band structures for both monolayer and bilayer ZGNRs are shown in Fig. S16 in the SM [62]. On the other hand, we note here that outside the low-energy region (where there is more than one electronic band) we find for all systems that reflection and interband scattering play a larger role in the electron transport through these devices, as the number of open channels (modes) grows with energy. In other words, it was not possible to identify conditions for realizing BS or M at energies with multiple subbands in the GNRs. Therefore the following discussion is focused on the energy window corresponding to a single (conduction or valence) band, since the most interesting physics related to the electron quantum optical features were identified here.

D. Quality of the realized mirrors and beam splitters
To obtain a qualitative picture across all the possible systems of the most suitable candidates for BSs or Ms, we construct in the following a figure of merit (FM). On the one hand, we look for candidate systems where a significant part of the scattered electron wave can be transferred to the other ribbon, i.e., that T 13 or T 14 is large. We encode this property in the quantity τ ≡ max(T 13 , T 14 ). On the other hand, for a suitable BS or M it is important that the reflection and transmission to a third electrode should be small. This property is encoded as a "loss" function λ ≡ R 1 + min(T 12 , T 13 , T 14 ).
Our FM is then defined as FM = e −20λ tanh 1 20   We use a linear color scale where BSs (FM = −1) as black, Ms (FM = 1) as red, and uninteresting systems (FM = 0) appear as white. We set FM equal to zero whenever there is more then one band per GNR at the energy considered (as it happens, e.g., for large values of |E − E F |). In that case the sum of all transmission probabilities is equal to the number of bands and thus larger than 1. This case is not useful for the devices we have in mind, though a more careful analysis may show how to also use the systems in this energy range. In other words, λ determines the intensity of the plots while τ sets the color. The FM is chosen to be highly selective: it decays to about 1/2 of the maximum value for losses (=transmission probability into the undesired output ports) of about 3%. Similarly, the FM of a lossfree, but unbalanced BS is reduced to FM = −1/2 at an imbalance of about 57:43. Figures 5 and 6 show the FM for ZGNRs and AGNRs from the metallic 3p + 2 family, respectively, as a function of electron energy and ribbon width W . Overall, these figures show that the most interesting systems are those composed by ZGNRs or AA-stacked AGNRs. For both types of GNRs one can find devices that behave as BS or as M, respectively. For instance, Fig. 5 reflects that the 8-ZGNR devices shown in Fig. 2(a,b) are good candidates for BS, consistent with the qualitative picture of Fig. 3(a,b) and Fig. 4. For both AA and AB ZGNR devices the transmitted electron wave to the other ribbon is always bent 120 • into electrode 3 (see also the full transmission matrices in Fig. 5 and Fig. S3 [62]). To obtain a M, where an electron incoming from electrode 1 is almost entirely transferred to electrode 3, one should choose wider ZGNRs.
For the AGNRs the situation is a little more complex. As discussed in Fig. 2 it is possible to form four differ-ent stackings (AA-1, AA-2, AB-1 and AB-2). Further, the band gap of AGNRs is determined by the overall ribbon width W , which classifies them into three distinct families 3p, 3p + 1, or 3p + 2 for integer p [24,25,27,28]. This leaves us with 12 different situations, considered in terms of the full transmission matrices in Figs. S3-S15 [62]. We find that the most interesting devices are those built with (3p + 2)-AGNRs in the AA-stacked configurations. However, compared with the ZGNRs, the parameter space for desirable devices is more restricted and the losses are generally larger. Independent of width, the AB-stacked configurations lead to scattering of the electron wave into all terminals.
We also note here that the qualitative difference mentioned in Sec. III A between the 60 • turn of the transferred electron wave for AA-1 configuration versus the 120 • turn for AA-2 is a robust feature across the different families (Figs. S3-S15 [62]). Additionally, we also find very thin white regions, that do not correspond to high losses but to T 12 ∼ 1, immersed in red -e.g., seen for W = 10-15 in Fig. 5(b) and for W > 20 in Fig. 6(b). This suggests that crossed GNRs can also work as energy filters. These T 12 (T 13 ) peaks (dips), also plotted in Fig. S17 for clarification, become narrower as the width of the ribbons grows, which enhances the energy selection.

E. Robustness of the discussed properties
So far we have discussed the different transport properties that can be found in the ideal case, that is commensurate GNRs (AA-or AB-stacking) with a relative angle of θ = 60 • . However, precise control of the device geometry is likely a significant experimental challenge. In this section we therefore proceed to test the robustness of the transport properties against some perturbations of this ideal situation. More specifically, we explore now the exact roles of the intersection angle, deviations from the idealized stacking pattern, lattice deformations via uniaxial strain, variation of the inter-GNR separation, and electrostatic potential differences between the two ribbons.
Since we concluded above that ZGNR devices may be the best candidates for building electron quantum optical setups, we will focus the following discussion around them. We take as the reference device the crossing of two AB-stacked 8-ZGNRs (Fig. 2a) and compute the transmission probabilities from terminals α = 1 to β = 1, 2, 3, 4 for each of the above mentioned perturbations. The AA-stacked 8-ZGNRs were found to display qualitatively similar trends as can be seen from the Figs. S18-S22 [62]. We will see that the low-loss property of these devices is thus preserved for the applied variations and in some cases the FM is even significantly enhanced, indicating that almost perfect BS or M could be obtained by tuning the above mentioned parameters.

Intersection angle
We first discuss the effect of small rotations of the ontop ribbon starting from the ideal configuration where θ = 60 • . To isolate the effect of the intersection angle from that of the precise stacking pattern (translation), we apply the rotation around the center of the scattering region (crossing) indicated with a black dot in Fig. 7(a). This ensures that the center of the junction remains unchanged and the effect of the rotation angle perturbs mostly the edge zones of the crossing. Figure 8 shows the reflection and transmission probabilities for varying angles δθ = ±2 • . The results for the reference case of θ = 60 • is shown as black lines in all panels. We first note that the reflection probability R 1 does not vary much from its initial value ∼ 0. The same holds for the (unwanted) transmission T 14 . The main effect is the precise distribution between the transmissions T 12 and T 13 .
This shows that the angle can be a physical knob to tune the transmission ratio between the two outgoing terminals of a BS. On the other hand, the approximate particle-hole symmetry found for the ideal AB-or AAstacking goes away as the lattice mismatch grows. The reflection plane shown in Fig. 2a is also lost for δθ = 0 (and other geometrical distortions), however we still identify only 4 qualitative independent elements in T for all cases.

Lateral translations
To study the precise lattice matching in the crossing area, we performed a series of calculations where the top GNR is translated by ∆ x along the x-axis with respect to the bottom GNR, see Fig. 7(b). Due to periodicity it is sufficient to consider translation vectors with modulus ∆ x ≤ 2a sin(60 • ) ≈ 2.46Å. Figure 9 shows the reflection and transmission proba- bilities as a function of such translations. Again, the results for the ideal AB-stacking is shown as black lines. As for small variations in the intersection angle, even though this geometrical distortion also intensifies the particlehole asymmetry as the system goes away from the ideal stacking, R 1 and T 14 remain rather unaffected by translation. In other words, the low-loss situation of these devices is robust with respect to translations. On the other hand, the inter -ribbon transfer process of electrons be- comes mostly less effective. We interpret this as due to an overall elongation of interlayer atom distances. For this reason T 13 slightly decreases with the translating of the on-top ribbon, while T 12 slightly increases with respect to the reference curves (black lines) for most of the cases.

Uniaxial strain
For experimentally grown GNRs it is relevant to consider the strain-induced deformations, e.g., a lattice mismatch with the supporting substrate [63]. But strain can also be applied in a controlled way [64] for example using a piezoelectric substrate for shrinking or elongating samples by applying a bias voltage [65]. In these directions we study here a simplified scenario of applying the same uniaxial strain ε to both GNRs in the device as defined in Fig. 7(c). As explained in the case of variation of the intersecting angle, to isolate the effect of strain on the transport properties of the device, we apply the strain with respect to the center of the crossing area (as depicted in Fig. 7c). Otherwise arbitrary lattice mismatches could further modify the transmission curves. The main effect of uniaxial strain is that it induces an anisotropy between the atomic bonds and therefore in the electronic structure of the individual GNRs. Addi-tionally, a strain induces some mismatch of the lattices in the crossing region, and therefore changes the scattering potential. The transport properties of the devices are therefore expected to be sensitive to strain. Figure  10 explores uniaxial strain in the range from −1% (compression) to 1% (stretching). Again, both R 1 and T 14 are not affected by the lattice deformation, and remain very close to zero in the single-channel energy region. Looking at the intra-and inter -transmissions T 12 and T 13 the curves vary smoothly around the reference values (black lines). The effects of compression and stretching of the GNRs is quite different: GNR compression causes T 12 (T 13 ) to increase (decrease), while stretching has the opposite effect. Again, strain can be seen as a physical knob to engineer the device properties. For instance, a strain of ε ∼ 1% brings the system closer to the ideal BS with T 12 = T 13 = 50%, while keeping both R 1 ∼ T 14 ∼ 0. In fact, our FM graph of panel Fig. 10e shows a significant enhancement of the performance of the device as a BS when stretching the device.

Interlayer separation
The exponential distance-dependence of electron transport in the tunneling regime suggests that the separation between ribbons may considerably affect the transport properties. Figure 11 shows the reflection and transmission probabilities as a function of the GNR separation d within an interval determined by ±0.15Å around a typical van der Waals distance d = 3.34Å [43,66,67] (black lines in all panels). Apart from the flat-band energy region very close to E = E F , the loss channels characterized by R 1 and T 14 are largely unaffected.
The main effect of varying d is to control the ratio between the intra-and inter -transmissions T 12 and T 13 , which varies smoothly to almost 0:1 as the ribbon separation d is decreased. In the other direction, the ratio goes (unsurprisingly) to 1:0 as the ribbon separation is increased and therefore eventually become decoupled.
The strong variation with the inter -GNR separation suggests that this is a key parameter to tune the transport properties. An ideal 50-50 BS may thus be obtained by applying an external force to the junction for d ∼ 3.30 A. While a perfect M is found for d < 3.20Å, as seen in Fig. 11e, where the plateaus at FM = 1 show this behavior. The possibility to use such electromechanical switching has been also proposed to be used for suspended multilayer graphene [68], crossed AGNRs [43] and crossed carbon nanotubes [69].

Electrostatic potential differences
Here we discuss the effect of an electrostatic potential difference between the two ribbons. This could for instance correspond to an experimental situation where a bias voltage is applied to the GNR electrodes. We consider a potential difference V that modifies uniformly the onsite energies to i − E F = −V /2 (and consequently the chemical potentials of the electrodes) of the top (red) ribbon and i − E F = V /2 of the bottom (blue) ribbon (see Fig. 1). Figure 12 shows the reflection and transmission probabilities for the range |V | ≤ 0.5 V. Drastic changes are observed in the energy range between the electrode chemical potentials, [−V /2, V /2], where valence bands (VB) and conduction bands (CB) of the two GNRs now overlap. In fact, the mixing of VB and CB leads to an interchange of the propagation direction: A transferred electron in the bias window turns 60 • instead of 120 • . In fact, our FM (Fig. 12e) shows that the performance of the device is enhanced in the energy window |E − E F | ≤ V /2, compared to the unbiased case (black curve). In contrast, the single-channel energy region slightly shrinks, as the chemical potential shifting produces the edge of the single-mode part of the CB (VB) of the bottom (top) ribbon to overlap with more than one mode in the top (bottom) ribbon. The presence of multiple bands in any of the incoming or outgoing electrodes is responsible for the larger reflection and transmission into the other output, e.g., as it happens for energies |E − E F | > 1.0 eV in panels Fig. 12(a,d).
Outside the bias window the curves are hardly changed, reflecting a low variability of the transport properties even when the elastically transferred electron wave to the other ribbon is now propagating with a different momentum due to the energy offsets of their band structures.

IV. CONCLUSIONS AND OUTLOOK
In this paper we studied the electronic transport properties of 4-terminal devices formed of two intersecting GNRs with a nominal crossing angle of θ = 60 • . We presented a complete classification and characterization of the different functionalities that can be found in these type of junctions by varying the edge topology of the GNRs (zigzag or armchair), stacking sequence (AA or AB), width of the ribbons, and energy for the propagating electrons in the valence or conduction bands.
We determined the number of independent transmission probability matrix elements in Eq. (10) that fully characterize their transport behavior: 10, 6 or 4 depending on the degree of symmetry that a given device displays. In practice, however, we found that for low-energy electrons it suffices qualitatively to describe the transmission probabilities with only 4 independent elements. The reason for this is the fact that the dominant part of the scattering potential contains more symmetries than that of the device geometry as a whole. Implicitly, this result also means that the strict geometrical symmetry behind the systems is not critical for the GNR crossings to function as beam splitters.
Besides the BS property, we also identified other interesting electron quantum optical functionalities of these devices. For instance, depending on the GNR width and electron energy the device can also behave as a mirror or an energy filter.
Interestingly, for AA-stacked AGNRs we discovered that there exist two different configurations (AA-1 and AA-2) that show little geometrical difference but behave very differently from each other in terms of the electron transport for low-energy electrons. In the particular case of 3p + 2-AGNR crossings, the electron beam is only allowed to turn 60 • for the AA-1 configuration, as opposed to to 120 • for the AA-2 configuration. On the other hand, AB-stacked AGNR devices do not show good electron quantum optical features due to the comparatively larger losses and low inter-GNR transmission. Unfortunately, AA-stacked configurations are probably harder to realize in practice (not the most stable energetically) compared to the AB-stacked one [70]. Combined with a generally larger variability of the AGNR transport behavior, these facts indicate that ZGNRs are more interesting objects for the considered device applications than AGNRs.
We further tested the robustness of the predicted transport properties by studying small variations on the intersecting angle between the ribbons, lattice matching in the crossing area, uniaxial strain, interlayer separation, and finite potential differences for devices composed of 8-ZGNRs. While the overall qualitative behavior was found to be robust under these modifications, a strong quantitative response can be obtained -indicating the need to control these effects as well as there potential for tuning the crossed-GNR devices. On the other hand, in this work we considered the situation of a spindegenerate electronic structure. However, ZGNRs have been predicted to develop spin-polarized states localized along the edges of the ribbons close to the Fermi level [24]. This suggests that additional spin-dependent effects could emerge in these devices. The interplay with the physics discussed here could become an interesting topic for future research.
For electron quantum optics applications, the central feature of the considered devices is that they coherently distribute incoming electrons in the intended output ports. In our model, with a precisely given unitary scattering matrix and without considering environmental degrees of freedom, all the considered devices process the input coherently. The analysis of the operative decoherence processes in GNR-based devices is an important task for future work. In particular, a single pure-state electron injected into one arm of a BS device discussed here is mapped to an (mode-)entangled state of the output ports. Such entanglement could be verified experimentally, for example by measuring the state's Bell correlations as discussed in [13]. A second basic application of the BS device is the Hanbury Brown-Twiss setup [8][9][10][11], which can be used to study the indistinguishability of electrons prepared in different input ports by the observation of anti-bunching in the output ports of the BS. A theoretical analysis of these experiments would include the investigation of the influence of environmental degrees of freedom (phonons, electrons in the substrate, or fluctuating perturbations as the ones discussed in Sec. III E), and, in the case of the Hanbury Brown-Twiss setup, also the effect of the interaction between electrons in the BS. An important prerequisite for all such experiments are methods to inject single electrons in a well-defined mode and to reliably detect them.
Finally, we envision that the functionalities presented here may be interesting as fundamental building blocks in larger electronic networks based on GNRs. For instance, with four GNRs one could construct the electronic analogue of the Mach-Zehnder interferometer, consisting of two beam splitters and two oriented mirrors at the intersection of pairwise parallel ribbons. Such a versatile setup, sensitive to the relative phase shift between two spatially separated pathways, have a wide range of quantum technology applications, e.g., metrology, entanglement, cryptography, and information processing [18]. In this appendix we compare the results presented in the main text with DFT, another popular theoretical approach used in the field of solid state physics. In particular, we choose to compute the electronic structure of AB-stacked bilayer graphene as a model system to establish suitable parameters for the general TB Hamiltonian introduced in Sec. II. We further simulate the electron transport characteristics of the specific device geometries shown in Fig. 2 for detailed bench-marking.
We employ the self-consistent DFT and NEGF methods as implemented in the Siesta/TranSiesta [55,61,71] packages. All calculations of this kind used the vdW density functional [72] with the modified exchange by Klimeš et al. [73]. The core electrons were described with Troullier-Martins pseudopotentials [74] and a double-ζ basis set defined with a 30 meV energy shift was used to expand the valence-electron wave functions. The fineness of the real-space integration mesh was defined using a 350 Ry energy cutoff. All carbon atoms were saturated at the edges with hydrogen atoms. Figure 13 shows the calculated electronic bands along the Γ-K-M-Γ path of the Brillouin zone of AB-stacked bilayer graphene obtained with Siesta [55]. Given the usage of a double-ζ basis set, the orthogonal σ and π bands have simple representations in terms of the {s, p x , p y } and {p z } basis orbitals, respectively. To map the DFT electronic structure onto the effective TB model in Eqs. (2)-(5), it is thus sufficient to consider only the p z part of the DFT Hamiltonian. Since we are interested in the low-energy physics, we fitted the TB bands inside an energy window of |E − E F | ≤ 2 eV using non-linear least squares and obtained the following optimal hopping parameters used in the main text: t = 2.682 eV, t = 2.7 meV and t ⊥ = 0.371 eV. The corresponding TB bands with these parameters are plotted in red dashed lines in Fig. 13, showing a very good agreement in the energy range of relevance in this work. Albeit unnecessary for the purposes here, we note that the significant deviations at the π band edges can readily be improved with a non-orthogonal TB model by introduction of additional parameters for the overlap matrix.
Having fixed the parameters for the TB model we proceed to compare it against the derived transport properties from DFT-NEGF for the six characteristic devices shown Fig. 2. Figure 14 shows the computed reflection and transmission probabilities from TB (solid lines) and DFT (dotted lines) within an energy window of |E − E F | 1.5 eV. Apart from different magnitudes of the AGNR band gap (known to be related to edge effects ignored in this TB modeling [25]) the two models only show minor numerical differences. Overall the two models provide very similar shapes and quantitative results for the transmission functions. From Fig. 13 and Fig. 14 we therefore conclude that the TB method used in the main text provides an accurate description of the essential physics in the energy range we are interested here.

S1. BOND CURRENTS
In this section we analyze the transport properties of multi-terminal devices in real space by computing the bond currents [1], defined as where H ij denotes the matrix element of the Hamiltonian of Eq. (2), and A α (i, j) is the matrix element of the spectral density of scattering states [Eq. (9)] for electrons incoming from lead α = 1, between nearest neighbor atoms i, j.
There is an implicit energy dependency on J ij and A α . Similar results are shown in Fig. 3, where the spectral density of the scattering states are plotted. However, the spectral density of states, defined in Eq. (9), also contains the contribution to the DOS of non-propagating (localized) states, that do not contribute to the electron transport. For this reason, we complement those results by plotting the current flowing between the different pairs of atoms, as defined in Eq. (S1), where we see the real space distribution of the propagating scattering states.
The bond currents in Fig. S1 were obtained with TBTrans

S2. SCATTERING POTENTIALS
In this section we analyse in more detail the scattering potentials created by the inter-GNRs coupling between the crossed ribbons for each of the devices of Fig. 2. The black dots in Fig. S2 indicate the atoms of the two ribbons that lie one on top of the other (stacked atoms), i.e., that possess the same xy-coordinates. Similarly to Fig. 2, we show the reflection symmetry planes (red dashed lines) that leave the geometries of Fig. S2 unchanged. One thing that is worth mentioning, is that not only the geometry generated by the overlapping atoms (Fig. S2) determines the symmetry, but in principle also their local environment, especially for those located at the borders of the intersection. However, the atoms that lie one on top of the other will give the main contribution to the scattering potential, as they contribute with the strongest interatomic coupling elements. For this reason the symmetries indicated in Fig. S2 apply approximately to the full problem.

S3. TRANSMISSION MATRICES
In order to complement the results presented in the main text, we compute here the transmission matrices as described in Eq. (10) for all edge terminations and stacking configurations.In the main text we showed the example of AB-stacked ZGNRs (Fig. 4). Here, Figs. S3-S15 provide the analogous results for all the other cases. In addition, there are three families according to the width of the ribbon for the cases of AGNR devices (W = 3p, 3p + 1, and 3p + 2) [3][4][5][6]. Therefore, the transmission probability matrix is plotted for each configuration (AA-1, AA-2, AB-1 and AB-2) and family separately. We only show results for energies where there is only one band, i.e., white regions in Figs. S3 -S13 correspond to energies where the number of bands is zero (gap) or larger than one (multiple electronic bands).

S5. TRANSMISSION PEAKS AS A FUNCTION OF THE RIBBON WIDTH
In this section we show the reflection (R 1 ) and transmission (T 12 , T 13 , T 14 ) probabilities as a function of the ribbon widths for AA-stacked ZGNRs and AA-2-stacked (3n + 2)-AGNRs. In Fig. S17 we plot these probabilities for ZGNRs of W ∈ [8, 16] C atoms (panel (a)), and for AGNRs of W ∈ [8, 32] C atoms (panel (b)) as a function of the incoming electron energy.

S6. ROBUSTNESS OF TRANSPORT PROPERTIES FOR AA-STACKED ZGNRS
In the main text we presented and discussed the variability of the transport properties of AB-stacked 8-ZGNR devices in Sec. III.E against some perturbations. Here we complement those results with the same calculations performed on AA-stacked devices (Figs. S18-S22). All graphs are compared to the reference case (AA-stacked 8-ZGNRs), which is plotted in black lines.