Quantum transport in graphene nanoribbon networks: complexity reduction by a network decimation algorithm

We study electronic quantum transport in graphene nanoribbon (GNR) networks on mesoscopic length scales. We focus on zigzag GNRs and investigate the conductance properties of statistical networks. To this end we use a density-functional-based tight-binding model to determine the electronic structure and quantum transport theory to calculate electronic transport properties. We then introduce a new efficient network decimation algorithm that reduces the complexity in generic three-diemnsional GNR networks. We compare our results to semi-classical calculations based on the nodal analysis approach and discuss the dependence of the conductance on network density and network size. We show that a nodal analysis model cannot reproduce the quantum transport results nor their dependence on model parameters well. Thus, solving the quantum network by our efficient approach is mandatory for accurate modelling the electron transport through GNR networks.

In this publication, we study GNR networks of different sizes and densities. We tackle the question how transferable quantum and (semi-)classical results are in the mesoscopic range where quantum effects may be relevant. For this, we perform quantum transport (QT) calculations as well as nodal analysis (NA) in comparison. This work is organized as follows: In section section 2 we first introduce the model system of the GNR network. Second, we give a brief overview about the general basic QT equations. Third, we present our new recursive network decimation scheme algorithm to efficiently calculate the key quantities. We then extend the linear complexity scaling of the standard recursive algorithm for linear chains to higher dimensional sparse networks. Finally, we present results of GNR networks for varying network density and size (network base area). In section 3 we first introduce the NA model adjusted for our GNR quantum network purpose. Then we present results of GNR networks with identical geometric parameters as for the QT calculations and compare the two approaches. We conclude that a single NA model cannot replicate all QT features, summarize our results and end with an outlook in section 4.
2 Quantum transport in nanoribbon networks

Model System
The model system we study are GNR networks, see figure 1, which are constructed as follows. A given number N of GNRs are randomly distributed within the network base area A by randomly selecting a point and a direction. Here, we focus on one type of nanoribbon: The networks consist of zigzag GNRs with a width of 3 unit cells (also called 6-zGNR), which equals 12 carbon atoms (corresponding to a width of 1.3 nm). We use this type of GNR because it is large enough to describe realistic devices yet sufficiently small to obtain fast results. Each GNR has a length of 20 unit cells, which equals 5 nm. They are passivated to avoid unphysical dangling bonds, i.e. a hydrogen passivation is The colour of the ribbons corresponds to their z-coordinate, ranging from red (lowest) to yellow (highest). Blue cells mark the semi-infinite electrodes where the alignment of the carbon atoms (black) and the attached hydrogen atoms is indicated.
used as this occurs during the fabrication process [26]. Thick strips in figure 1 represent a single GNR (color coded for its spatial position), each of its small sub-stripes marks one unit cell. Each new GNR is placed in the lowest possible layer. This means that if it would intersect with an GNR already placed, it is placed at the next higher layer with a vertical distance of 3.35 (red / orange / yellow in figure 1) as this corresponds to the distance between two graphene layers [27]. That means, the GNRs are assumed to be flat, and curvature effects are neglected. In reality, overhanging ends of stacked GNRs could bend. The resulting curvature would lead to a reduced in-plane transmission and thus a reduced conductance. This effect is in the range of up to 40% reduction, depending on the curvature [28]. Furthermore, additional connections between GNRs could occur, which would effectively lead to a larger network density ρ, thus counteracting the bending effect. Curvature effects are not included here, not least to keep the model simple.
At the edges of the network base area, protruding GNR unit cells are shifted to the opposite side of the network imposing periodic boundary conditions. For this system, the dimensionless network density is defined as with the area covered by a single nanoribbon Ω = 6.39 nm 2 (for 6-zGNRs) and network area A. For example, = 1.5 means that the GNRs could fill one and a half layers if aligned properly. Furthermore, one ribbon at the left and right network boundary is elongated infinitely to act as an electrode within the transport formalism (indicated in blue in figure 1).

Quantum transport theory
Electronic transport is described quantum mechanically by using quantum transport theory [29] in combination with an underlying electron structure theory. The conductance in the limit of a small bias can be computed from the transmission based on the Landauer-Büttiker formula [30] Here, G 0 = 2e 2 /h is the conductance quantum, T (E) the transmission function, and f (E) the Fermi distribution function. In order to calculate the transmission function, the system is divided into a finite channel (C) and two semi-periodic electrodes (L -left and R -right) as depicted in figure 2. Each part is described by a Hamiltonian matrix H L,C,R and coupling matrices τ LC,CR . The resulting channel Green's function is where S C is the channel overlap matrix, S LC,CR are the coupling matrices of the overlap, and Σ L = (τ CL − ES CL )G L (τ LC − ES LC ) and Σ R = (τ CR − ES CR )G R (τ RC − ES RC ) are the self energy corrections due to the coupling of C to the electrodes L and R. Finally, η is a small numerical value for improving convergence, which shifts the singularities at eigenenergies into the complex plane.  The electrode surface Green's functions G L/R are calculated efficiently using the renormalization decimation algorithm (RDA), which is a fast recursive algorithm [31,32]. We use a value of η = 10 −5 for calculating G C and for G L,R via the RDA.
Finally, the transmission function is where describe the broadening of the original channel states due to the coupling to the electrodes.
The Hamiltonian and coupling matrices themselves are computed in advance by an electronic structure calculation. To treat the large network systems studied here, we use the so-called non-self-consistent DFTB model [33,34], which readily provides the distance-and orbital-dependent Hamiltonian and overlap matrix elements. For describing the organic molecules present here, the parameter set 3ob [35,36] is used. This set uses a sp 3 -basis appropriate for our purpose and has been verified against other organic molecules, especially sp 2 -hybridised ones, with a similar structure. The cutoff distance was chosen to be a CO = 3.7 Å. This value includes the 3rd nearest-neighbor interaction in plane but excludes the 4th one, implying rather small GNR unit cells and reducing the decimation calculation time to an optimum. In addition, it ensures the inclusion of 2nd nearest-neighbor interaction out of plane, i.e. between two GNRs in different layers.

Network Transport: Reducing the Complexity through the Network Decimation Scheme
Using the aforementioned methods, the conductance for an arbitrary system can in principle be obtained for a given, adequate parameter set. However, due to the matrix inversion in (3) the number of necessary operations scales cubic with the number of atoms in the system. This is alleviated by dividing the channel into smaller cells to calculate an effective channel Green's functionG C that is equivalent to G C in (3) but of lower dimension.
Proficient algorithms benefit from this subdivision and result in a linear scaling with the number of atoms. For example, the recursive Green's function formalism (RGF) [10,11] has been widely used for similar calculations of linear chains [16][17][18][19][20][21][22]. Two slightly different versions of the RGF, the renormalization decimation scheme (RDS) and forward iteration scheme (FIS) are utilized [37] and adopted to the networks studied here, as shown in figure 3.
These methods are combined in the network decimation scheme introduced here to reduce the complexity of the network incrementally by virtually decreasing the number of nodes that need to be considered. The starting point is a system where each GNR unit cell, cf. figure 1, corresponds to one node.  The RDS is given by the equations The indices i, j, k refer to the respective nodes. τ i,j is the coupling matrix from node i to node j. The presence of a node causes a shift of all the neighbouring energy levels. Using RDS, this shift is applied to the Hamiltonian and coupling matrices of the surrounding nodes. Thus, after applying (5) node i is effectively eliminated as all its information is now stored in its neighbours and the connections between them. The FIS takes advantage of the fact that the overall Hamiltonian is sparse and that the electrodes only couple to a small part of the system. It is defined by which represents the influence of the electrodes and the neighbourhoods of the cells propagating through a linear chain of cells. In the endG C is obtained, which can be understood as an effective Green's function of the system.
As a network of cells as shown in figure 1 is far from constituting a linear chain, (5) and (6) are combined to obtaiñ G C . In order to minimise the number of computer operations, they are applied as follows: We start the decimation scheme at the outermost edges, where the nodes have only one neighbour (see circled cells in figure 3a). These cells are successively decimated using RDS, until only clusters and connections between them remain (see figure 3b). Then, connections between clusters are treated, where nodes only have two edges each. Subsequently, the RDS is applied to the clusters of many interconnected cells, where the nodes with the fewest edges are decimated one at a time (see figure 3c). In the end, only the cells connected to the electrodes remain present (see figure 3d), as they can't be decimated using RDS. In this final step, FIS is used and G C is obtained. We implemented this QT algorithm in python using numpy and scipy. For further details we would like to refer the interested reader to [38].
This hierarchic decimation procedure works especially well for weakly connected networks, i.e. networks with few overlapping areas compared to the total network. The standard recursive procedure, which is often used for linear chains, would provide a subdivision in transport direction only. For linear chains this results in an O(N ) scaling with the number of atoms N . For 2D and 3D networks, it yields a O(N 2 ) and O(N 7/3 ), respectively, scaling, because only one dimension with N 1/d atoms scales linear and the other dimensions with N (d−1)/d atoms scale cubical due to the matrix inversion. The decimation scheme presented here treats this more efficiently as all directions are subject to decimation. In the limit of weakly connected networks, this implies a vanishing fraction of overlapping areas (node clusters in figure 3) compared to the original network, and the algorithm scales like O(N ) independent of the network dimension. The complexity scaling is summarised in table 1.
However, such a O(N ) scaling is not achievable for dense, realistic systems in general. Clusters with a high degree of interconnections scale with O(N 2 ) and contain typically between 50 to 90% of node sites. Nevertheless, in comparison to the standard procedure with a N 2 scaling for 2D systems, a scaling between N and N 2 can be achieved using the novel network decimation scheme, which allows one to treat larger systems by quantum transport methods, as the one we discuss here. Of course, in the limit of denser networks, the linear scaling breaks down. In the extreme limit where the whole network is one cluster, a quadratic scaling O(N 2 ) is obtained, which is the same as for 2D systems treated with the standard algorithm of linear chains. For 3D systems, a slightly better scaling may be achievable, but this is out of the scope of our work as we focused on quasi 2D networks with few layers.

Results
For each set of geometric parameters (GNR type, area, number of GNRs) 500 random and percolating networks were generated, for which the QT calculations are performed and statistically analysed. Percolating or non-percolating networks could arise due to the finite interaction distance given by the cutoff distance a CO .
We first state that the usual percolation behaviour is observed: Figure 4 shows the percolation probability as a function of the network density for different network sizes (base areas A denoted by colour). There is a minimal network density below which the percolation probability is exactly zero because the cumulative length of the GNRs is smaller than the side length of the network base area. At the percolation threshold network density the percolation probability increases markedly and tends towards one for large densities. The overall trend can be approximately described by a logistic function fit (red line in figure 4) [39]. A finite-size effect can be seen as larger network base areas result in steeper curves, that is, sharper transitions. However, percolation aspects are not the focus of this work, and we refer the reader to [40,41], while we focus on percolating network ensembles in the following. From an engineering point of view, non-percolating networks would be "broken devices" and not relevant for applications. Figure 5 shows the average conductance of percolating networks as a function of the network density for various network sizes (base area A denoted by colour as before). For each data set (i.e., for a fixed base area A), two different regimes can cleraly be distinguished: (A) The "linear chain regime" for low densities < 1. Here, the networks are weakly connected and possess few, often only one transport path. They thus behave like an effectively linear chain with few perturbations. In the limit of a network consisting only of one GNR covering the whole system (not included in figure 5), this behaviour tends towards the ballistic transport regime with constant conductance G = 2G 0 . In a macroscopic picture for larger network sizes, this can be interpreted as a series circuit, where the conductance decreases roughly like G ∼ 1/N with N being the number of tunnelling regions between any two GNRs contribution to the chain. N increases with increasing effective chain length eff , which in turn increases with density as long as only one chain exists: N ∼ eff ∼ . Thus, the conductance decreases with increasing density G ∼ 1/ , as can be seen in figure 5, especially for small base areas A.
(B) The "interference dominated diffusive regime" for high densities > 1. Here, the networks are highly connected, i.e. a significant amount of the networks contributes to percolation and many transport paths are present. An increase of the network density results in further increasing the number of connections and transport paths, thereby increasing the conductance. In a macroscopic picture this can be interpreted as a parallel circuit, where the conductance increases roughly like G ∼ M with M being the number of parallel GNR paths. M increases with increasing number of layers or system height H, which in turn increases with increasing network density , leading to M ∼ H ∼ . Thus, the conductance increases with increasing density G ∼ , as can be seen in figure 5 for all base areas A. A linear regression for ≥ 2.5, i.e. the range where virtually all networks are percolating, has been done (solid lines in figure 5), yielding a slope ∂G ∂ = 2.55 · 10 −3 G 0 .  Both trends (A) and (B) can be seen best for the smaller network base areas, e.g. the system with a base area of 14 × 6 nm 2 . In between both regimes (A) and (B), a transition region must be present and is realised as a plateau-like minimum around ≈ 1.
In addition, the conductance for a fixed density depends on the network size (base area A indicated by colours as in figure 5). To this end we increased the area A = W L in width W and length L, while fixing the aspect ratio W/L = const. In the "linear chain regime" (A) increasing A leads also to an increased effective chain length eff and thus a decreased conductance, as can be seen in figure 5. In the "interference dominated diffusion regime" (B) increasing A leads on the one hand to an increased number of parallel paths due to the increased width W , which will increase the conductance G ∼ W . On the other hand the increased length requires more tunnelling events between two GNRs. But when the system gets very large the scattering events in the diffusive regime of 1D systems lead to enhanced interference and strong localisation effects, which will decrease the conductance exponentially with length L. In the mesoscopic range of the networks discussed here, the strong localisation regime may not be fully reached yet (as shown for similar purely 1D systems [20][21][22]). In summary, we get G ∼ H W/e L , which explains the increasing conductance with increasing ∼ H and the decreasing conductance with equally increasing W and L (i.e. increasing A).

Comparison with semi-classical transport and nodal analysis
In this chapter we will address the question to what extent our results for GNR networks could be reproduced based on a classical transport approach. Quantum-classical correspondence has been established as a useful tool in mesoscopic physics. The classical approach neglects all interference effects, and comparison with the quantum transport results allows one thus to classify their importance and a deeper characterisation of the system properties.

Model System
In order to approach the GNR networks in section 2 using classical transport based on nodal analysis (NA), the networks are represented by layers of two-dimensional (2D) polygons (see figure 6). Each polygon has a width and length corresponding to their atomic (quantum) counterpart, which is roughly 1.3 × 5 nm 2 in the 6-zGNR case. These stripes are now randomly placed on the base area A = LW , starting in the lowest layer. Analogous to section 2.1, overhanging stripes are shifted to the opposing side of the base area.
Two stripes are deemed interacting, if they are in adjacent layers and overlap in the x, y plane. The centre of the overlapping region determines the x, y coordinates of the two resulting nodes, one on each ribbon taking part in the interaction. These coordinates are then used to order the emerging nodes of a ribbon, so that each node is connected to its nearest neighbours. This way the network consisting of ribbons is represented as nodes and edges, which are then investigated using classical NA.  Using Kirchhoff's current law, an arbitrary circuit consisting of ohmic resistances can be expressed aŝ where the conductance matrixĜ represents the connections between the nodes, the vector ϕ indicates the potential at each node and I contains the net currents in each node. The currents are all zero due to Kirchhoff's first law, with the exception of two contacts where the current is inserted into or extracted from the system (see blue markers in figure 6). G must be chosen to best reflect the results obtained for the quantum system in section 2. Thus, the following definition is used: nodes not connected .
The choice for the j = k and the unconnected node case is motivated by the NA definition. Additionally, one needs to distinguish between intra-ribbon connections and inter-ribbon connections. The conductance between two nodes on the same ribbon is assumed to be G intra = 2G 0 , which corresponds to the conductance of an ideal GNR with ballistic transport. For two nodes on separate ribbonŝ is chosen, effectively describing electron tunnelling. This conductance is assumed proportional to the overlapping area A ij = A ji of the ribbons between two nodes i and j (see figure 6). The factor g t is to be optimized as to yield the best agreement with the quantum transport results obtained in the previous section 2. Figure 7 depicts the percolation probability of the NA approach for different network base areas. The qualitative behaviour is the same as described in section 2.4. Quantitatively there are some small differences in the threshold density and slope in comparison to the QT results that can be related to the differences between the models used. Figure 8 shows the conductance as a function of the network density for different network base areas (denoted by colour), calculated with the NA approach. Analogous to the QT results in section 2.4, two regimes can be found for each data set (fixed base area).

Results
(A) The "linear chain regime" for low densities < 1. The general qualitative trends are similar to the ones of the QT results in figure 5: Weakly connected networks in the limit of only one chain yield a conductance G ∼ 1/L ∼ 1/ that decreases with increasing length and density, which in a macroscopic interpretation corresponds to a series circuit. The decrease with density and base area A can be seen in the diagram. The effect is less pronounced than for the quantum results. Nevertheless, the limit of only one lead spanning the whole system G = 2G 0 is the same (not shown).
(B) The "ohmic regime" for high densities > 1. Also here, the trend for fixed base area is similar to the quantum transport calculations. Highly connected networks with many transport paths yield an increasing conductance G ∼  H ∼ with increasing height and density, as can be seen in the diagram. In a macroscopic interpretation this corresponds to a parallel circuit. In contrast to the quantum treatment before, however, the base area A dependence is now different: Although the conductance increases with width W as before, the length dependence differs because interference effects and thus an exponential conductance reduction with increasing length cannot be captured in NA calculations. Thus, a (only) roughly inverse length dependence G ∼ 1/L is found. In total, we find the Ohmic behaviour G ∼ HW/L. As we kept the aspect ratio W/L constant for all A, the conductance G depends only on height H and not on area A. This can nicely be seen in figure 8 as all curves roughly fall together. The solid line represents a linear regression in this ohmic regime for ≥ 2.5 where all networks are percolating. The corresponding derivative yields the slope for the specific parameter g t = 0.1 G 0 /nm 2 . As the derivative depends on g t , it could be adjusted such that the NA value equals the QT value, however, only for this specific value g t . To this end one could perform multiple NA calculations with many different g t and choosing the best fitting g t by minimising the difference between the QT and NA conductances. Nevertheless, this g t will be specific for the set of geometric parameters chosen. A generic value of g t , that brings the NA calculations in agreement with QT calculations for arbitrary geometries, cannot be found. A comparison of the density dependence of the conductance G( ) in the quantum transport (QT) and nodal analysis (NA) model reveals a semi-quantitative agreement in the behaviour for low (linear chain regime). However, for high characteristic deviations are visible, the most striking being that the conductance G does not depend on the network size (or base area A) in the NA model due to the ohmic behavior. In contrast, quantum interference effects depending on the network size A change this behavior in the QT model and result in an A dependence that cannot be reproduced  with the NA model. However, the NA parameter g t can be tuned such that a best fit to the QT case is obtained for a (single) given network of base area A. Finally, we briefly discuss the distribution of the conductance within a statistical analysis. This is shown in figure 9 for QT (blue) and NA (orange) within a) the low and b) the high density regime. Comparing both regimes, the distribution is much broader for low than for high . This is due to the fact that in the linear chain regime (low ) the conductance depends sensitively on the specific tunnelling geometry. For the large amount of geometries where the GNRs are barely touching each other the conductance is fully determined by the corresponding tunnelling distance. For the NA there is a hard cutoff and thus no such long tail is present.
In the high density regime, many transport paths exist within the network. Consequently the electrons can choose between different pathways to avoid tunnelling-restricted geometries, which leads to a suppression of the low-conductance tail of the distribution in figure 9b) in comparison to figure 9a). However, the QT calculations allow for strong localization effects as discussed before. Thus, specific geometries show an exponentially suppressed conductance even in the high situation. This explains the low conductance contributions for the QT model (blue) in figure 9b).

Summary and Outlook
The network decimation scheme introduced in this paper represents an efficient decimation scheme for quantum transport, which enables the treatment of comparably large quantum networks, that otherwise would not be accessible using standard recursive procedures or even direct inversion. The conductances of GNR networks with varying density were calculated using this novel algorithm and compared with results from nodal analysis as a classical transport method. Although these two approaches may be able to yield similar results for specific geometries in the largedensity regime for ≥ 2.5 by parameter tuning, single parameterization of a nodal analysis model cannot generically reproduce the QT results with all features for different geometric parameters. For low densities quantum effects such as tunnelling and interference dominate and yield an increased conductance, which is only partly accessible by semi-classical nodal analysis. The dependence on the network size (base area A) for high observed for QT cannot be reproduced by NA. In turn, the numerically and conceptually cheap nodal analysis can be expected to yield qualitatively reliable results for specific GNR networks with given, fixed geometric parameters and for ≥ 2.5. In summary, the efficient network decimation algorithm introduced here is both necessary to accurately calculate the conductivity of GNR-networks and allows one to extent full QT calculation to large systems that were previously only accessable by less correct semiclassical methods.
The work presented here will be the basis for manifold future investigations ranging from accessing extended parameter regimes (aspect ratio of the network, GNR width, or GNR type) to other material systems (carbon nanotubes or other types of nanowires), and is not restricted to carbon. The extension to other networks, for example the simulation of porous materials, is possible. Eventually, special geometries can be investigated, such as periodically repeating structures, bent/curved structures or multi-terminal devices.