First Principles Study of Electronic Structure and Transport in Graphene Grain Boundaries

Grain boundaries play a major role for electron transport in graphene sheets grown by chemical vapor deposition. Here we investigate the electronic structure and transport properties of idealized graphene grain boundaries (GBs) in bi-crystals using first principles density functional theory (DFT) and non-equilibrium Greens functions (NEGF). We generated 150 different grain boundaries using an automated workflow where their geometry is relaxed with DFT. We find that the GBs generally show a quasi-1D bandstructure along the GB. We group the GBs in four classes based on their conductive properties: transparent, opaque, insulating, and spin-polarizing and show how this is related to angular mismatch, quantum mechanical interference, and out-of-plane buckling. Especially, we find that spin-polarization in the GB correlates with out-of-plane buck-ling. We further investigate the characteristics of these classes in simulated scanning tunnelling spectroscopy and diffusive transport along the GB which demonstrate how current can be guided along the GB.


Introduction
Grain boundaries (GBs) are line defects that occur in monolayer graphene grown by chemical vapor deposition (CVD) on common catalytic surfaces such as Cu where non-epitaxial growth from multiple nucleation centers causes orientation mismatch at the boundary between individual grains or domains [1,2,3].This is illustrated in Fig. 1.In these one-dimensional regions where the grains interface, Figure 1: Cartoon illustration of grain formation in CVD-grown graphene, inspired by Ref. [3].The four different colors indicate different angles of orientation of the graphene lattice in each grain.Part B) is not to scale in A) the two graphene lattices deform and reconstruct into disordered structures which exhibit varying degrees of order along the boundary [1,4,5,6,7,3,8].As CVD growth of polycrystalline graphene is the most common means of producing large quantities of graphene for use in electronics or optoelectronics [9,10,11], the impact of GBs on the overall and local transport properties has been an area of research for more than a decade [12,13,14,15,16,17,18,19,20,21,22,23,24].In addition, the one-dimensional interface between two misoriented graphene grains has been shown to have several surprising and useful properties, including Yu-Shiba-Rusinov states and spontaneous time-reversal symmetry breaking [25,26], localized voltage drops [27,28], optical response [29,30] and they could be used as electron wave guides [31].Previous studies have considered the geometric effects in two lattices twisted relative to each other, which gives rise to two different types of GBs; conducting and non-conducting, depending on if the GB is momentum mismatched or not [32].On the other hand, quantum interference effects related to how the carbon atoms bond and to the local electronic structure at the GB in general can also have significant influence on the conduction of electrons [33,34].
There are many different possible ways in which two graphene sheets with a different orientation can interface and reconstruct, and the details of how the atomic structure influences the sheet conductivity of polycrystalline graphene is still not well understood.Here we take an approach inspired by "big data" and generate a set of 150 GBs to get a deeper insight into the effects of the atomic arrangements, their classifications and trends.

From pristine graphene to a grain boundary
The mathematics of interfacing two graphene grains with a GB is closely related to how moiré cells are constructed [35], with the notable simplification that only one lattice vector needs to be common between the two sides of the GB.Having this condition fulfilled allows the construction of a GB as sketched in Fig. 1B and now formalise how to satisfy this condition in a simple way for two graphene lattices with a relative rotation and strain.Let the unit vectors of the hexagonal cell of the primitive graphene unit cell be a 1 and a 2 , let A = [a 1 , a 2 ], let M, N ∈ {S ∈ Z 2×2 |det(S) ̸ = 0}, let R θ be a rotation in 2D, ϵ a scalar strain, and search for solutions to M = m 00 m 01 m 10 m 11 N = n 00 n 01 n 10 n 11 µ = 0 or 1, ν = 0 or 1 The action of the rotation matrix R θ and the supercell matrices M and N is seen in Fig. 2. Once two matching lattices with a small tolerable strain have been found, we construct an initial guess for the atomic coordinates.When the solution is found it is trivial to order the coordinate system so that the ydirection is along the GB, while the two lattice vectors describing the graphene periodicity on the left and right hand side of the GB are freely chosen from the rotated primitive lattice vectors.Methods for constructing GBs have also been considered previously and have then been relaxed using various potentials [36,37,4].However, here we furthermore perform a relaxation using DFT to achieve more realistic GB models.The steps of the workflow devised for this work are listed in Table 1: Table 1 0. Find Solutions using random initialisation (2D) • Loop over possible M, N, and by trial find a good initial guess for θ. • Solve eqs.(1) by a numerical solver to obtain (M, N, θ, ϵ).

Geometry Creation (2D)
• Repeat one cell to the left, repeat the other cell to the right.Remove atoms that are within 0.85 times the carbon-carbon bond length d CC to other atoms.This length has been chosen ad-hoc.
• Place atoms at reasonable distances from each other, removing atoms that are too close and adding atoms where there are holes to be filled.A computationally cheap cost function V (see SM.9.1) depending only on distances and angles to the surrounding neighbor atoms is used.An atom can be inserted if V ({r i } ∪ {r new }) < V ({r i }).Atoms further than ∆ = 8 Å away from the GB where the two flakes meet are considered as frozen and are not considered further in the calculation.

Force field relaxation (3D)
• Geometry refinement with force-field relaxation using GULP [38] and the Brenner potential [39] for the atoms within ∆ of the GB.• If more than 2 counts of bond-angles are outside the interval of [θ min , θ max ] = [100 • , 160 • ] the calculation is discontinued, and a new calculation is instead started from step 0. This interval is introduced ad hoc so that most bonds are closer to the 120 • bond angle of graphene.

Geometry relaxation using DFT (3D)
• Spin-degenerate DFT relaxation using the SIESTA code [40] of the atoms within ∆ of the GB.Atoms outside this region are moved as a single, internally fixed structure on each side.Singlezeta basis is used and forces are required to be below 0.01eV/ Å.

Final geometry relaxation using DFT-NEGF (3D)
• Spin-degenerate DFT relaxation using TranSIESTA [41] which employs open boundary conditions and avoids crosstalk between neighboring GBs in a calculation with periodic boundary conditions.A single-zeta polarized basis is used and forces are required to be below 0.01eV/ Å.This workflow is also visualized in Fig. 3 which shows the intermediate steps.In Fig. 4 an example of an initial structure is shown together with the intermediate and final steps.It can also be seen that the atoms close to the GB move substantially, while the atoms further away move in a rigid way.
Figure 3: Flowchart of the workflow outlined in Table 1 for GB construction.
In the end of the workflow, the resulting GB may have relatively large outof-plane buckling.A GB buckling in the z-direction can be a possibly sound configuration, but since the relaxing geometry has been constrained to be flat when |x| > ∆, any in-plane and out-of-plane distortions are restricted to the region with |x| < ∆.
A selection of the relaxed structures coming out of the workflow can be seen in Fig. 5, where the out-of-plane buckling also has been defined as ∆z = max i {z i − z avg }.Furthermore, the initial structures have been chosen such that the length of the period of the GB is smaller than a maximal value L max GB in order to restrict the required computational resources.

Non-equilibrium Greens Functions and Electronic Transport
The theoretical framework for the DFT-based structural relaxation and electronic structure is the non-equilibrium Green's function framework implemented in the TranSIESTA code [41].The central object in this framework is the retarded Green's function, expressed in the device orbital basis, is given as [41,42,43] The Green's function is at z = E +iη (for a numerically small η) from the device overlap matrix S ky , device Hamiltonian H ky , electrode self-energies Σ r ky (z) = α Σ r α,ky (z) corresponding to the open boundary conditions in the non-periodic directions.A general matrix in k y -space C ky can be written as a Bloch sum Figure 4: The individual steps of refinement of the workflow from Table 1 performed after solving the problem outlined in eq. ( 1).Black lines indicate the unit-cell of the pristine graphene sheet on either side of the GB.where C n is the collection of hopping matrix elements from the zero'th to the n'th unit cell.In the LCAO basis, this is the general prescription for how to calculate a matrix in k y -space for a structure periodic in the y-direction.It is also useful, cf.Eq. ( 3) to define the unitless ky = k y R y ∈ (−0. Equation ( 4) looks like the quasi-particle equation in many-body theory, but there are no many-body contributions to the self-energy in this case.The lifetime instead comes from the finite size of the simulation region.We will call the states obtained by solving eq. ( 4) localized states (LS) from now on.Furthermore, the matrix element determines the life-time τ .The state ψ i ky is an approximate eigenstate of the system Hamiltonian with a LS half life τ i ky .The self-energies Σ r α,ky are calculated from a bulk electrode Hamiltonian using a recursive algorithm for the calculation of the surface Green's functions of the electrodes [44].This algorithm is implemented in e.g TBtrans or sisl codes [45,41].Given that the GBs considered here are periodic in the transversal ydirection, the real-space self-energy in a single cell Σ R can furthermore be found by a k y -integral of the device Green's function G r ky over the transverse Brillouin zone (TBZ) [46] as, Further computational details can be found in the SM.subsection 9.2.
In the general case when having calculated the electrode self-energies, the electrode broadening matrices Γ α,ky (z) = i Σ r α,ky (z) − Σ a α,ky (z) can be evaluated.This makes it possible to calculate the electronic transmission function from left to right as [47] T T (E, ky ) d ky .
Equation ( 9) gives the elastic transmission probability accounting for the electron wave interfering with itself as it moves through the structure.The transmission function calculated with TBtrans will in following plots be normalized with respect to the associated GB period, L GB .
Lastly, the spectral function A α,ky can be calculated from the Green's function and the electrode broadening function, given as [41] A It is useful for quantifying bound states in the system, as A α,ky only contains the states which couple into electrode α, and thus does not contain the contribution from the LS at the GB.

Density Functional Theory
The final step of relaxing the structures and the transport calculations is performed using the PBE exchange-correlation functional [48], an energy cutoff of 150Ry and a SZP basis set in TranSIESTA [41].A simulation of a system being doped e.g. by a back-gate or dopant atoms, is done using mixed pseudopotentials, also known as the virtual crystal approximation [49].Here a weighted sum of carbon and nitrogen or boron is used to add or remove charge from the system.

Atomic Geometries
The relative angle of rotation θ of the two lattices forming the periodic GB, together with the out-of-plane buckling are basic parameters that influence the conductive properties.A histogram of the angles contained in the data set is shown in Fig. 6A.This data set spans a wide number of angles, but exhibits gaps, for example in the vicinity of θ = 30 • where no GB structures were identified.A scatter-plot of the maximal out-of-plane buckling of the atomic z-displacements can be seen in Fig. 6B.From these figures it appears that a majority of the 150 GBs produced shows a significant buckling within the structure, while only 48 GBs are found below the line at ∆z S = 0.5 Å in Fig. 6B.The significance of out-of-plane buckling is to enable hybridization between the π-and σ-subsystems that are decoupled in bulk graphene.This might also be the reason for the lack of structures around the line ∆z S = 0.5 Å in Fig. 6B, where the carbon atoms either tend to hybridize in a sp 3 configuration with buckling or by sp 2 in-plane bonding.Notably, the GB structures below this line cluster around relatively few angles.The angles available in this data set are limited by the maximum GB length L max GB since it is certainly possible to generate GBs at smaller angles, they just have a very long periodicity, thus requiring more computational resources.The whole data set presented in Fig. 6 is available in Ref. [50].

Conductive Properties and Electronic Structure
The transmission of current across GBs is important for the sheet resistance of large-area polycrystalline graphene films and the performance (and variability of performance) of any devices made thereof.We consider here the transmission function T (E) per unit length of GB.The Landauer-Büttiker formula determines the current density running across the GB at zero temperatures by [41] where s is spin and T (E) beyond linear response also depends on the applied voltage.G 0 = 2e 2 h and h is Planck's constant.The quantitative behavior of the current density relative to the twist angle is shown in Fig. 7A and it displays a large variance of how well the GBs conduct current.Some GBs conduct current as well as graphene ("transparent regime"), while in other GBs, the current is numerically negligible ("blocking regime").Yet more GBs are somewhere in between these limits ("opaque regime").The categorisation in Fig. 7A is of course subjective and the transition between opaque and transparent is smooth.The cause of the blocking behavior for some GBs is momentum mismatch between the available electrode states on each side as identified in ref. [32].The intermediate, "opaque regime", is interesting because a combination of available scattering states on either side and wave interference effects cause the transmission function to be significantly suppressed.
The current density in Fig. 7A is summed over spin, meaning spin-dependent effects are not visible.In Fig. 7B the spin filtering of the GBs in the data set is quantified though the Transmission Spin Filtering(TSF) coefficient, defined as A plot of the maximal value of TSF(E) with −0.5eV < E < 0.5eV for the various structures is seen Fig. 7B, showing that out-of-plane buckling makes some GBs spin-polarize and can work as a spin-filter.In the most extreme case the TSF comes out at close to ∼ 100% in Fig. 7B, indicating a GB that completely polarizes the current.We can identify four overall types of GBs in this dataset, which are summarized in Table 2. 3. Blocking GBs J < 0.05J pristine with transport gaps arising from mismatch in crystal momentum [32].
4. Spin-filtering GBs that spin-polarize and transmit a preferred spin.This case can fall into any of the previous three categories.
In Fig. 8A/B/C/(D+E) the electronic structure of the GBs in Figs.5A-D is characterised by the DOS obtained from the imaginary part of the Green's function.The DOS is resolved in both energy and transverse k y -point in the TBZ, and the plot also shows the LS lifetime τ i ky from eq. ( 4).The Dirac cone of graphene is clear in all five panels of Fig. 8.There are, however, also other band-like features.These come from LS around the GB.These states act as a separate subsystem with its own electronic structure and they couple very weakly to the bulk graphene electrode states.The structures shown in Figs.5A-D fall into cases 1-4 of Table 2, respectively.Ayuela et al. [51] used tight-binding calculations to show how the number of LS bands depends on the edge-states of the two semi-infinite graphene sheets that connect via the GB, and thus classified these in terms of the edges of either side.In this work we instead choose to describe the GBs in terms of their conductive properties cf.Fig. 8 and Table 2.We may characterize the localized GB states using the quasi-particle formalism and associate a life-time cf.eqs.( 4) and (5).The intrinsic LS lifetime in graphene has been found in experiments to be above ∼ 0.5 ps corresponding to η ≈ 10meV and a mean free path on the order of 100nm [52].In our calculation the LS lifetime is an artefact due to the truncation of the LS in the finite device region between the electrodes (Always >40 Å until the electrodes start on either side of the GB).The very long lifetime we obtain in the calculation (above ∼ 0.5ns) compared to intrinsic graphene means that we should consider these as localized.On the other hand we may use the relative change in the computed lifetimes as a measure of their decay and overlap with the electrode regions.With this in mind, we see how the lifetime (circles in Fig. 6) grow with distance to the electrode states in the Dirac cone.In Fig. 8A the bands merge with the Dirac cone as they get close to the projection of the graphene K-point on the GB direction.This is, however, not the case in Fig. 8B where signatures of the bands associated with the LS can also be observed inside the region of the Dirac cone.A similar situation appears in Fig. 8D and Fig. 8E.
The effects of the LS bands intersecting the Dirac cone are clear from T (E, k y ), shown in Fig. 9.The presence of the LS bands facilitates transmission at the intersection point, but also gives rise to a sharp decline when moving above/below the LS bands, see Fig. 9B-D.Lastly, the LS of the first conduction band of the GB in Fig. 8A is plotted in Fig. 10A.The GB is symmetric around x = 0 and the state shown in this figure has negative parity i.e. it changes sign across the GB.In Fig. 10B the LS of the flat band in Fig. 8C is shown.This state is instead skewed towards one of the electrodes with a longer decay length and thus stronger truncation, which explains the shorter lifetime of the flat-band in Fig. 10B compared to the state in Fig. 10A.Fig. 8C furthermore shows a peculiar situation with a different number of bands on either side of the k y = 1/3 point.In all DOS-plots, except Fig. 8C, the Dirac cones of either side of the GB are on top of each other, while in Fig. 8C the geometry is such that the Dirac cone is centered at different k y -points, which shows how this case is momentum-mismatched yielding a transport gap in Fig. 9C.The bands seen in Fig. 8 give rise to a DOS that is significantly modified relative to pristine graphene, with the notable introduction of van Hove singularities in the DOS where the LS bands are flat.This behavior has for instance been observed in ref. [53,18] by voltage-dependent differential conductance ( dI dV ) measurements by STM-spectroscopy on top of the GB.

Non-equilibrium effects
However, in the transport gapped case, shown in Fig. 11C-E, the LS bands are not seen to be broadened, but are rather modified and moving in position with bias.In particular, a bending of the LS band running from ky = 0.0 to ky = 1/3 appears as a consequence of the finite bias.This behavior is qualitatively different from Fig. 11A and Fig. 11B and stems from the different locations of the Dirac point in the TBZ.Thus, the bias over the GB changes the band-structure and group velocity of the transverse, localized one-dimensional GB state.In a local transport measurement along the GB, this band-bending could manifest itself as an enhanced transversal conductivity when a bias is applied from the left to right.This in essence means an electric field in the xdirection will modify the local conductivity close to the GB in the y-direction.
The momentum mismatch between the electrodes surrounding the GB, leading to a gap in the transmission, may be lifted by phonon scattering supplying the missing momentum.This inelastic channel may, at finite bias, result in local Joule heating around the GB.This effect has been experimentally observed near GBs in graphene [54].We here qualitatively treat this phenomenon using the Special Thermal Displacement (STD) method [55].We use QuantumATK [56,57] and the Brenner potential [39] to generate the STD atomic displacements for various temperatures within a region of width ±1nm around the GB.In order to include the right phonon momentum for a system with a gap corresponding to ky = 1/3, we consider phonons in a 3 times wider supercell.
Using the STD method we can calculate the transmission function for the wider supercell, with and without the presence of phonons at a given temperature which, in turn, determines the statistical atomic displacements along the various phonons.We use this as a simple way to quantify the impact of the phonons on the conductive properties.In Fig. 12A,B,C the transmission function is shown at various temperatures with the STD.As the temperature rises, the pristine transmission function gets modified slightly by the presence of phonons.However, in the transport gapped case, the phonons furthermore have the effect of allowing current to run through the GB.This originates from the fact that the STD has a longer period than the primitive cell initially considered for the GB in Fig. 5C, which allows mixing between the two Dirac cones from Fig. 8C that previously were momentum mismatched.The E and k y resolved transmission function and DOS can furthermore be seen in Fig. 12D,E,F,G, and shows in detail how this longer periodicity induced by the STD results in a modified transmission function which has a visible imprint of the GB DOS.The mechanism for altering the transmission function is the same also for the transparent and opaque cases, where the transmission function also gets modified when the DOS of the GB is large.For example, the dip from the STD in Fig. 12A around E = 0.3eV is from the edge of the first conduction band shown in Fig. 8A.

Effects of Electrostatic Gating
Graphene can be electrostatically gated using a back-gate (under or above the graphene sheet) to induce a change of charge in the graphene, shifting the Fermi level in the Dirac cone.Because the DOS of the bulk states in graphene will be the determining factor for where the Fermi-level will be located given a gate voltage, the Fermi-level in the GB will follow accordingly.We do not account for any re-relaxation of the GB from gating.Fig. 13 shows the effects of electrostatic gating on transmission function the GBs in Fig. 5 where the Fermi-level is moving approximately ±0.3eV from the charge neutrality point.The "transparent" GB from Fig. 5A whose k y -resolved DOS is seen in Fig. 8A is not subject to significant changes to its transmission function under introduction of additional holes or electrons by gating.The blocking GB in Fig. 13C displays only minor changes with gating.However, for both the "opaque" and "spin-polarized" GBs the transmission functions change significantly with gate voltage.The reason of the dips and peaks in the transmission function close to the charge neutrality point has previously in Fig. 8 and Fig. 9 been attributed to the presence of the localized states and their bands in the GB.When charge carriers are introduced into the GB, the bands bend significantly, as is evident from Fig. 14B and Fig. 14E.This band bending makes the bands intersect the Dirac cone in different places, giving rise to an altered transmission function as shown in Fig. 14C and Fig. 14F.The results demonstrate how the transport for the large class of opaque GBs especially sensitive to gating.

STM Simulation
As mentioned above STM spectroscopy has been used to investigate the local density of states around GBs based on dI dV -measurements [53,5,58,6,18,59].The resulting dI dV curve is a measure of the local DOS right under the tip [60].Here we calculate the transmission function for a STM tip positioned a distance z tip over a GB.Here we use the real space self-energy from eq. ( 7) to model the infinite extent of the graphene sheet, and a regular gold tip [61] to model a STM tip in contact with the graphene GB at a height of z tip = 3.5 Å.This makes it possible to simulate a STM-measurement where there is a tunnel contact between the tip and the GB.The geometry is shown in Fig. 15 and for the DFT calculation a SZ basis is used.When choosing how many times to tile the minimal cell for the Bloch folding of the Green's function in eq. ( 7), the GB length is taken to be > 3.5 nm in all cases.The resulting transmission function, pristine GB DOS and GB DOS with tip is shown in Fig. 16.The correspondence between the pristine DOS and peaks in the transmission functions is clear from these figures, even though the LS states broaden a bit when the tip is included in the calculation.Furthermore, the transmission function has sharp peaks in the unoccupied part of the spectrum for structures Fig. 5A and Fig. 5B indicating that the LS bands above the Fermi-level couple strongly to the gold tip.The states in the unoccupied part of the spectrum are however much less pronounced in the transmission function, even through we know from the DOS they are present.
Large-scale tight-binding simulation (using parameters as in ref. [62]) is shown in Fig. 17.These demonstrate conduction along the GB for electrons injected on a single atomic site 7nm away from the GB (see SM. subsection 9.4 for computational details).

Kubo-Greenwood approach to transport along the GB
The conductivity of the GB in the y-direction related to the one-dimensional GB states is troublesome to evaluate by the Landauer formula, so instead a response function, σ, can be calculated using the Kubo-Greenwood formula, given as [63] with A cell being the unit cell area, ⟨ • ⟩ ky denoting the TBZ average, and v y ky the transverse component of the velocity operator in the tight-binding representation, which can be taken as the k y -derivative of the tight-binding Hamiltonian [64,63].σ yy describes the response to a static electric field as J y = σ yy E y where J y is the current running per unit cell of the simulation domain.The Green's function contains both bulk and GB contributions, where the bulk contribution will increase with size of the cell chosen in the x-direction, but the GB part will not.This approach enables a quantification of the conductive properties of localized states with dispersion transverse to the standard transport direction (computational details are presented in SM. subsection 9.3).
Using the Kubo-Greenwood formula eq. ( 13), a diffusive conductivity for the y-direction can be calculated for the same three GBs as in the STM simulations Fig. 16.The Kubo-Greenwood formula in a periodic solid evaluates to a sum over k-derivatives of the available bands at the considered energy [65].In the present case, there are additional LS band(s) available, which lead to notable peaks in the transverse conductivity in Fig. 18.Here the peaks are instead located in the occupied part of the spectrum.It is also noteworthy that the GB with the transport gap a nonzero conductivity at E = 0 coming from the flat-band in Fig. 8C.

Discussion
The bands that lie outside the graphene Dirac cones are a recurring feature of all the GBs and the different types of band-structures from Fig. 8 should therefore be possible to observe experimentally.The LS bands may also show up as polarization dependence in optical measurements using scanning probes, since the electric field below the tip can point along the GB or perpendicular to it.In the first case the localized states can conduct current along it, while being just a scatterer in the second case.This type of experiment has already been carried out in ref. [29] and might also show additional features related to polarization direction and gating as well as depend on the class of GB, c.f. Fig. 8 and Fig. 14.
The phenomenon of the GBs conducting current well along the GB shown in Fig. 18 could also be detectable by scanning magnetometry to map the local current density as has been done in e.g.ref. [66].This type of experiment should show a large current density at the GB if the Fermi-level can be adjusted to match one of the peaks in Fig. 18.Furthermore, as described in refs.[67,68], quantum capacitance measurements are directly related to the electronic density of states.This in turn means the LS bands might be detectable in quantum capacitance measurements of samples which are densely populated by a certain type of GBs.In particular the transport gapped GB from Fig. 8C and the band running from one Dirac cone to the other might show a significant peak in its quantum capacitance at V = 0, given that there are enough of these GBs to contribute significantly to the total sample DOS.This has also been calculated in relation to GBs in ref. [69].Evidence of Peierls instability [70,71] has also been found in TMD GBs [72] but there is (the the authors' knowledge) very limited evidence in graphene GBs [73].It is however a possibility because a periodic lattice distortion could in some cases lower the total energy.We have however not made any considerations of longer-range relaxation when creating the workflow from Fig. 3 and relaxing the GBs.
The band-bending induced by gating that was predicted in Fig. 14 should also be experimentally detectable with an STM measurement, since the van Hove singularities of the band edges have already been detected and during gating these peaks should move.The fact that the band structure of the GB states in some cases can be manipulated could make GBs a test-bed for exotic physics, such as a flat band that is tuneable by bias, see Fig. 11C and Fig. 11D or a gate tunable spin-filter, see Fig. 7B and Fig. 13.
The structures generated in this work are created without the inclusion of contaminants such as water and other specimens present during the CVD growth.Molecules can however subsequently adsorb onto the graphene surface during storage, device processing or measurement [74,75].In some particular cases, such as hydrogen and some other species, it is favorable to make bonds to carbon atoms inside the GB [76,77].The inclusion of such adsorbates, while interesting, would many-double the number of calculations needed to quantify the effects as there are different outcomes depending on where the absorbate bonds in the GB [78].However, identifying which GBs belong to which of the four classes listed in Table 2 is useful for knowing which GBs to select for a given property or function of a graphene device, even through it may be modified .Followingly growth optimizing conditions for generating a higher prevalence of this type of GB could then be carried out.
Another unknown factor that could impact validity of the modeling in this work is the fact that the GB phonons couples to some degree to the electrons around the GB as was demonstrated in Fig. 12.It is therefore also plausible that inelastic effects could significantly modify the ability of the GBs to conduct current.However doing the inelastic calculations have been outside of computational feasibility for the structures considered in this paper, but as previously stated, Joule-heating have been observed at GBs [54], which is an inelastic phenomenon.

Conclusions
An algorithm for automated GB generation has been constructed and used to build a dataset containing 150 GB structures (Fig. 6).These GBs have been grouped into four categories according to their ability to conduct (spinpolarized) current (Fig. 7 and Table 2).The conductive properties are shown to be linked to the local electronic structure at the GB, namely the LS states existing there, which in some cases are a major inhibitor of the GBs ability to conduct current.These LS states can in some cases be significantly affected by charge doping, and in the case of a transport-gapped GB, be manipulated by applying a bias to the device.The conductive properties have been characterised in several ways, using both a regular transport setup for conduction directly across the GB (Fig. 9) in conjunction with electrostatic gating (Fig. 13 and Fig. 14), in addition to with STM-simulations (Fig. 16) and diffusive conductivity calculations (Fig. 18).The electronic structure has furthermore been characterised in equilibrium (Fig. 8, Fig. 10 and Fig. 14) and out-of-equilibrium (Fig. 11).
This work is our attempt at providing a detailed classification of graphene GBs.This should be useful for understanding and reducing detrimental effects and for creating unique transport properties in nanoscale systems with potential for device applications.

Data availability Statement
The structures used for this article and code used to generate these can be found in ref. [50] and further data can be retrieved by request at abalo@dtu.dkby 9 Supplementary Materials This cost function is used to determine whenever to remove or add a carbon atom at various places, and has been constructed empirically for the purpose of yielding reasonable initial guesses for the GB structures.The condition for placing an atom is if V ({r i } ∪ {r new }) < V ({r i }).The first part of the potential takes care of the distance dependence of the energy when the distance between carbon atoms are changed, and the second part takes makes the result tend to form bonds are closer to 120 • as is the case in the hexagonal lattice of graphene.

Computational Details for Real Space Self Energy
The k y -integral in eq. ( 7) is calculated using the adaptive integration routine quad vec of the scipy integrate module, NumPy arrays and the Numba package [79,80,81].An upfolding of the minimal cell greens function can furthermore be employed in this calculation [46].This calculation has been automated in the siesta python code (ref.[82]) in an energy-parallel fashion and with the possibility for left and right electrodes that are different in number of orbitals and possibly electronic structure.The siesta python code mainly consists of an ASE-inspired [83,84] python class for handling the writing, processing and reading input-files for electronic structure and transport codes, primarily the SIESTA and TBtrans codes.

Computational Details for Diffusive Transverse Conductivity
The Greens function going into eq.( 13) can furthermore be calculated using the efficient algorithm for calculating the inverse of a block-tridiagonal matrix [41,85].This algorithm is available in the Block matrices code (See ref. [86]).An optimized pivoting scheme for making the matrix block-tridiagonal is obtained using the reverse Cuthill-McKee algorithm and the TBZ average is calculated using an adaptive integrator, both from SciPy [87,79].Wrapper functions are available in the siesta python code (ref.[82]).

Large-scale Simulations
For the bond-current calculations in Fig. 17 a simple tight-binding model using a nearest-neighbor model with hopping t = −2.7eVwith CAPs having been placed at < d CAP and y > L y − d CAP [62].The tip self-energy is modelled as a single site chain coupling only to one carbon site.The hopping element is t chain = −2.7eVmaking it effectively a wide-band self-energy in the energy range considered.

Figure 5 :
Figure 5: Examples of GBs with varying period and twist-angle.The structures are obtained using the workflow in Table 1.Pristine graphene semi-infinite electrodes are attached to the left and right.A) The commonly found 5-7 GB with θ = 21.79 • .B) A variation of the 5-7 GB with θ = 26.00• .C) Long GB with out-of-plane buckling and with θ = 25.28 • .D) Long GB containing two buckling carbon atoms and with θ = 17.89 • .Different scales in each subfigure.

Figure 6 :
Figure 6: Statistics of the dataset.A) Histogram of GBs at each angle.B) Maximal difference between highest and lowest (in z-direction) atoms in the GB structure.

Figure 7 :
Figure 7: A) Current density across a GB at a bias V = 0.2V without including non-equilibrium effects.The current density J = J ↑ + J ↓ is calculated using the equilibrium transmission function and a symmetric voltage drop, V R = −V L = 0.1V.B) Scatter-plot of the TSF against maximal buckling of the structure.Maximum taken over E ∈ [−1/2, 1/2]eV

Figure 8 :
Figure 8: A-D) GB DOS as a function of energy (E) and transverse wavevector (k y ) for the GBs shown in Fig.3A-D, respectively.The size of the black circles scale with the lifetime τ i ky of the LS (largest bubble corresponds to ≈ 0.5ns).In D) the two panels are for the separate spins of the spin-polarized GB.Note the different color scheme used for spin down in the fifth panel from the left.The color bar associated with the second color-scheme is the second from the left.

Figure 9 :Figure 10 :Figure 11 :
Figure 9: A-D) Transmissions as a function of energy (E) and transverse wavevector (k y ) for the GBs shown in Fig. 3 A-D, respectively, representing the transparent, opaque, gapped, and spin-polarized cases.A broadening η = 10 −3 eV has been used.

Figure 12 :
Figure 12: A,B,C) The k y -integrated transmission function with and without the STD at various temperatures.D,E) The k y -resolved DOS with and without the presence of a STD at T = 600K using the structure from Fig. 5C.F,G) The k y -resolved transmission function with and without the presence of a STD at T = 600K using the structure from Fig. 5C.

Figure 13 :
Figure 13: The transmission functions of the four example structures, calculated at different dopings indicated in the legend.Same labelling as in Fig. 5.The percentage is relative to the number of carbon-nuclei.The doping percentage corresponds to a majority carrier concentration of ≈ 10 13 cm −2

Figure 14 :Figure 15
Figure 14: DOS (A+C) and transmission function (B+D) for the "opaque" (A+B) and "spin-polarized" case(C+D) as a function of k y , showing bandbending effects.Both cases gated such that there are +0.25%electrons per atom.A and D are repeated DOS from Fig. 8.

Figure 16 :
Figure 16: Plots of STM setup DOS minus gold electrode spectral DOS and transmission function from graphene to gold electrode.A) Results of GB in Fig. 5A.B) Results of GB in Fig. 5B.C) Results of GB in Fig. 5C.ADOS is the spectral DOS of the gold electrode.

Figure 17 :
Figure 17: Bond-transmission for electron injection on a single site.GB is located at x = 0.The GB structure is the same as in Fig. 5A.The bondtransmissions are evaluated at E = 0.05eV.

Figure 18 :
Figure 18: Transversal diffusive conductivity of three GBs compared to the conductivity of pristine graphene calculated in units of G 0 /A U C = 2e 2 hA U C .