This document is the Accepted Manuscript version of a Published Work that appeared in final form in Nano Letters, copyright © American Chemical Society after peer review and technical editing by the publisher. To access the final edited and published work see: https://dx.doi.org/10.1021/acs.nanolett.9b03112. Universal spin diffusion length in polycrystalline graphene Aron W. Cummings,1, ∗ Simon M.-M. Dubois,2 Jean-Christophe Charlier,2 and Stephan Roche1, 3 1 Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and BIST, Campus UAB, Bellaterra, 08193 Barcelona, Spain 2 Institute of Condensed Matter and Nanosciences, Nanoscopic Physics group, UCLouvain, Louvain-La-Neuve, Belgium 3 ICREA–Institució Catalana de Recerca i Estudis Avançats, 08010 Barcelona, Spain (Dated: September 9, 2019) Graphene grown by chemical vapor deposition (CVD) is the most promising material for industrial-scale applications based on graphene monolayers. It also holds promise for spintron- ics; despite being polycrystalline, spin transport in CVD graphene has been measured over lengths up to 30 µm, which is on par with the best measurements made in single-crystal graphene. These results suggest that grain boundaries (GBs) in CVD graphene, while impeding charge transport, may have little effect on spin transport. However, to date very little is known about the true im- pact of disordered networks of GBs on spin relaxation. Here, by using first-principles simulations, we derive an effective tight-binding model of graphene GBs in the presence of spin-orbit coupling (SOC), which we then use to evaluate spin transport in realistic morphologies of polycrystalline graphene. The spin diffusion length is found to be independent of the grain size, and is determined only by the strength of the substrate-induced SOC. This result is consistent with the D’yakonov- Perel’ mechanism of spin relaxation in the diffusive regime, but we find that it also holds in the presence of quantum interference. These results clarify the role played by GBs and demonstrate that the average grain size does not dictate the upper limit for spin transport in CVD-grown graphene, a result of fundamental importance for optimizing large-scale graphene-based spintronic devices. The growth of graphene via chemical vapor deposi- presence of disordered networks of GBs leaves the sec- tion (CVD) is the most promising approach for realiz- ond hypothesis open as a possibility, which could have ing industrial-scale applications of this material1 . One profound consequences for the optimization of graphene- drawback of CVD-grown graphene is that it tends to be based spintronic devices. polycrystalline, with misoriented single-crystal domains Here we use numerical simulations to study the impact separated by grain boundaries (GBs) consisting of arrays of GBs on spin transport in polycrystalline graphene. of five-, seven-, and occasionally eight-member carbon We first develop an effective tight-binding (TB) model of rings. In some specific cases the GBs can be character- polycrystalline graphene in the presence of intrinsic and ized by a given periodicity, but more generally they tend substrate-induced SOC, which is derived from extensive to be complex meandering arrangements of these non- first-principles calculations. The model is based on sim- hexagonal rings2,3 . Charge transport and scanning tun- ulations of a variety of carbon-based haeckelites and is neling measurements have revealed that graphene GBs found to be general and transferable to the complex mor- serve as a significant source of charge scattering, lead- phologies of graphene GBs. We then use this model to ing to enhanced resistance4,5 and localization effects4,6,7 . perform spin transport simulations in realistic models of Only when the graphene grains become larger than 1-10 polycrystalline graphene, using an efficient linear-scaling µm do the GBs cease to dominate the charge transport methodology that gives direct access to spin relaxation in CVD graphene8 . and propagation. Our simulations reveal that the spin Graphene also has clear advantages for spintronic ap- diffusion length in polycrystalline graphene is indepen- plications, owing to its low intrinsic spin-orbit coupling dent of grain size and depends only on the strength of (SOC)9–11 . Combined with its high electron mobility, the substrate-induced SOC. This result is fully consis- this can lead to spin diffusion lengths as long as 30 µm tent with the D’yakonov-Perel’ mechanism of spin relax- in clean exfoliated graphene-based devices (with mobil- ation in the diffusive regime, but here it is also shown ity up to several 10, 000 cm2 /V-s)12 . Intriguingly, spin to be robust to the contributions of quantum interfer- diffusion lengths as long as 10 µm and spin currents mea- ence induced by disorder. These findings indicate that surable over channel lengths of 30 µm have also been in the presence of SOC the graphene GBs serve as scat- reported in disordered CVD graphene with much lower terers of charge, but they do not play a direct role in charge mobility13,14 . These measurements can be ex- spin relaxation. In other words, our results suggest that plained in two different ways: either none or very few grain size is not a limiting factor for spin transport in GBs were present in the measured devices; or GBs in CVD-grown graphene, a result of genuine relevance for CVD graphene, while impeding charge transport, have the future optimization of graphene-based spin devices little effect on spin transport and relaxation. Although and architectures in the context of memory or spin logic the former hypothesis remains plausible, the lack of a technologies. theoretical foundation concerning spin transport in the We first develop an effective TB model for the descrip- 2 tion of itinerant electrons in polycrystalline graphene, TABLE I. Tight-binding parameters for pristine graphene ob- which is derived from fitting to first-principles simula- tained from fits to first-principles calculations, as illustrated tions. The model includes the impact of electrostatic in Figs. 1(a-c). barriers and atomic disorder present at GBs, and ex- plicitly accounts for both intrinsic and extrinsic SOC. to1 β1 to2 β2 λI λR The first-principles simulations were carried out with the all-electron FP-LAPW method as implemented in the -2.414 eV 1.847 -0.168 eV 3.077 13.437 µeV 3.383 µeV V/nm Elk code (http://elk.sourceforge.net/). The self- consistent calculations with SOC have been carried out within the LDA approximation with a muffin tin radius of 1.316 Bohr for carbon atoms and an APW cutoff of 5.32 our focus is on the transport properties of low-energy Bohr−1 . A 33 × 33 k-point mesh was used to sample the itinerant electrons, the fitting has been limited to the first Brillouin zone of pristine graphene, and equivalent [−1, 1] eV energy window around the Fermi level. The k-point densities were used for the supercell calculations. optimized parameters are given in Table I. The radial The starting point of our one-orbital TB model is the dependence of the hopping integral is shown in Fig. 1(a). now well-known functional form of the second-nearest- The perfect agreement of the model with reference band neighbor hopping Hamiltonian of graphene. The Hamil- structure calculations is illustrated in Figs. 1(b,c). tonian can be conveniently decomposed into a kinetic op- To account for the crystallographic disorder introduced erator and a SOC operator, by dislocations at GBs, the TB model is enriched by introducing topology-dependent renormalization factors. Ĥ = Ĥkin + ĤSOC , (1) Within a simple ball-and-stick model, any sp2 lattice can be mapped to a tiling of the plane by polygons, where X † X † Ĥkin = t1 ĉiσ ĉjσ + t2 ĉiσ ĉjσ , (2) hi,ji hhi,jii each site is the shared apex of three polygons. In our case, it can be three hexagons as in graphene, or any i νij ĉ†iσ (ŝz )σσ0 ĉjσ0 X ĤSOC = √ λI combination of pentagons, hexagons and heptagons due 3 3 hhi,jii to the presence of disclinations in the lattice. For the 2i X † purpose of our model, we denote Li as the geometric en- + λR ĉiσ [E · (ŝ × dij )]σσ0 ĉjσ0 , (3) vironment of a given site i. Li can be any set of three 3 elements out of the ensemble of considered carbon poly- hi,ji gons E = {5, 6, 7}. Similarly, Lij are defined as the local where ĉ†iσ (ĉiσ ) is the creation (annihilation) operator for environment of a pair of neighboring sites i and j. If the the pz orbital with spin σ at lattice site i, s is the spin sites are first nearest neighbors, the vector joining i and Pauli matrix, E is the external electric field, dij is the j is the shared edge of two polygons and Lij can be any unit vector pointing from site j to i and νij is +1 (−1) for combination of two elements of the ensemble E. If the clockwise (counter-clockwise) hopping paths from site j sites are second nearest neighbors, the vector joining the to i. Single (double) brackets stand for summation over sites is inscribed within a polygon and Lij is reduced to the nearest (second-nearest) neighboring lattice sites. In a single element of E. Finally, an extra term is added to our model we assume an exponential dependence of the the on-site energy of the Hamiltonian, which describes t1 and t2 hopping parameters with respect to the intersite the local redistribution of charge around pentagon and distance, heptagon carbon rings, t1 = to1 e−β1 (r−a) and t2 = to2 e−β2 (r−b) . (4) ĉ†iσ i ĉiσ with i = X X H̃loc = ∆m . (5) Here, a (b) is the equilibrium distance between the first i m∈{Li } (second) nearest atomic √ sites in graphene at equilibrium (a = 1.42 Å, b = 3a = 2.46 Å). The third term is This term is crucial for the description of the resonance the intrinsic SOC operator which connects second nearest peaks introduced by GBs in the low-energy electronic neighbor sites, and the last term is the extrinsic (Rashba) spectrum. It also enables a geometry-dependent electro- SOC operator that arises from an external perpendicular static alignment of the GBs with respect to the graphene electric field E15 . This electric field can be directly ap- grains. Note that without loss of generality we can im- plied via a gate voltage, or it can be an effective field that pose the constraint ∆6 = 0, which conveniently defines arises from placing graphene on a substrate. The Rashba the reference energy as the charge neutrality point of pris- term connects nearest neighbor sites and linearly depends tine graphene. on the strength of the perpendicular electric field. The To account for the local variation of aromaticity and hopping integrals t1 and t2 , as well as the parameters β1 energy-momentum dispersion relation in the presence of and β2 ruling their distance dependence, were fitted to- GBs, the kinetic operator is also renormalized, gether with the SOC parameters λI and λR with respect ĉ†iσ t1ij ĉjσ + ĉ†iσ t2ij ĉjσ , X X to first-principles band structure calculations of graphene H̃kin = (6) at equilibrium and under isotropic strains up to 20%. As hi,ji hhi,jii 3 (a) (b) (d) (c) (e) (f) (g) (h) (i) (j) AAAC6XicjVHLSsNAFD2NrxpfVZduBqsgCCWVinYhFHXhsoKtghVJ4rQOnSYhmQil+APu3Ilbf8Ct/oj4B/oX3hkj+EB0QpIz595zZu69XiRFohznOWcNDY+MjuXH7YnJqemZwuxcMwnT2OcNP5RhfOS5CZci4A0llORHUczdnif5odfd0fHDCx4nIgwOVD/iJz23E4i28F1F1GlhqbXLpXJZi0eJkMSssy22ypxSpVplLcabzLbt00LRKTlmsZ+gnIEislUPC09o4QwhfKTogSOAIizhIqHnGGU4iIg7wYC4mJAwcY5L2KRNKYtThktsl74d2h1nbEB77ZkYtU+nSHpjUjIskyakvJiwPo2ZeGqcNfub98B46rv16e9lXj1iFc6J/Uv3kflfna5FoY1NU4OgmiLD6Or8zCU1XdE3Z5+qUuQQEafxGcVjwr5RfvSZGU1iate9dU38xWRqVu/9LDfFq74lDbj8fZw/QXOtVCa8XynWtrNR57GARazQPDdQwx7qaJD3Fe7xgEera11bN9bte6qVyzTz+LKsuzfglJpF ✏5 = +0.499 eV AAAC5HicjVHLSsNAFD2NrxpfVZcuHCyCq5KKqBuhqAuXFewDWilJOq1Dp0lIJoKIS3fuxK0/4Fa/RfwD/QvvjCn4QHRCkjPn3nNm7r1eJEWiHOclZ42NT0xO5aftmdm5+YXC4lI9CdPY5zU/lGHc9NyESxHwmhJK8mYUc3foSd7wBgc63jjncSLC4ERdRPx06PYD0RO+q4jqFFbbh1wql7V5lAhJzDbbY06JtRmvM9u2O4WiU3LMYj9BOQNFZKsaFp7RRhchfKQYgiOAIizhIqGnhTIcRMSd4pK4mJAwcY4r2KRNKYtThkvsgL592rUyNqC99kyM2qdTJL0xKRnWSRNSXkxYn8ZMPDXOmv3N+9J46rtd0N/LvIbEKpwR+5dulPlfna5FoYddU4OgmiLD6Or8zCU1XdE3Z5+qUuQQEadxl+IxYd8oR31mRpOY2nVvXRN/NZma1Xs/y03xpm9JAy5/H+dPUN8slQkfbxUr+9mo81jBGjZonjuo4AhV1Mj7Gg94xJPVs26sW+vuI9XKZZplfFnW/Tu6A5kj ✏6 = 0. eV AAAC6HicjVHBTttAFBzc0gYDJaXHXlaNkLgQOQgp9FApKhx6pBIJkWKEbGeTLmxsy15XQqgf0FtvVa/8AFf6JVX/oPwFs1tHokVVWcv27LyZ2X27ca5VaYLg54L36PHik6eNJX95ZfXZWvP5+qDMqiKR/STTWTGMo1Jqlcq+UUbLYV7IaBZreRSf7dn60UdZlCpLD815Lo9n0TRVE5VEhtRJsxXuS20iEcq8VJpMV7wRW0F753VXhEIOhO/7VAXtwA1xH3Rq0EI9DrLmD4QYI0OCCjNIpDDEGhFKPiN0ECAnd4wLcgWRcnWJT/DpraiSVERkz/idcjaq2ZRzm1k6d8JVNN+CToENejLqCmK7mnD1yiVb9l/ZFy7T7u2c/7jOmpE1+ED2f7658qE+24vBBLuuB8WecsfY7pI6pXKnYncu7nRlmJCTs3jMekGcOOf8nIXzlK53e7aRq/9ySsvaeVJrK9zYXfKCO39f530w2G53iN/vtHpv66tu4CVeYZP32UUP73CAPrM/4wrX+O6del+8r96331Jvofa8wB/Du7wFaqeaHQ== ✏7 = 0.497 eV FIG. 1. One-orbital topological TB model. Gold curves and dots correspond to first-principles results, blue curves are obtained within our topological TB model, and dashed green curves are obtained within the pristine TB model of graphene before topological renormalization. (a) Radial dependence of the hopping integrals considered in our one-orbital TB model. The green (purple) curve corresponds to the hopping between first (second) nearest neighbor lattice sites. (b) Low-energy band structure of graphene at equilibrium (left panel) and under 20% isotropic strain (right panel). (c) Band structure of graphene around the tip of the Dirac cones in the presence of SOC. The left panel shows the case with only intrinsic SOC, while the next three panels are for increasing transversal electric field (E = 1 V/nm, 2.5 V/nm, and 4.1 V/nm from left to right). (d) Low-energy band structure of the H5,6,7 periodic haeckelite structure. (e) Spin splitting of the two metallic bands of the H5,6,7 periodic haeckelite structure. (f-j) Ball-and-stick representation of the five haeckelite structures included in our fitting set. Colors have been added to represent the tiling of the plane by pentagons (red), hexagons (orange) and heptagons (blue). According to the literature16 , the structure represented in (f) is named H5,6,7 . with the fitting of our model, known as haeckelites16 , are il- Y lustrated in Figs. 1(f-j). These haeckelite structures have t1ij = t1 · (1 + ∆tm 1 ), been created via the incorporation of disclinations into m∈{Lij } 2 × 2 and 3 × 3 graphene supercells. The considered Y geometries have been fully relaxed with respect to both t2ij = t2 · (1 + ∆tm 2 ). (7) atomic and cell degrees of freedom. The topological pa- m∈{Lij } rameters obtained by fitting the low-energy band struc- tures are given in Table II. We impose ∆t61 = ∆t62 = 0 in order to preserve We note that the renormalization of the nearest neigh- the energy-momentum dispersion relation of pristine bor hopping integral is very weak. This directly reflects graphene. The full Hamiltonian of our model now reads the sp2 character of the fitting set, and additional bench- marks revealed that these renormalization factors con- H̃ = H̃loc + H̃kin + ĤSOC . (8) tribute only marginally to the improved description of sp2 carbon structures. The main contribution to the suc- The additional parameters of our topological model, cess of the model comes from the geometry-dependent {∆5 , ∆7 , ∆t51 , ∆t71 , ∆t52 , ∆t72 }, have been fitted renormalization of the onsite energies. Figures 1(d,e) il- against a set of first-principles band structures corre- lustrate the improved description of itinerant electrons sponding to periodic carbon sp2 systems containing five-, due to the renormalization. While ∆5 and ∆7 improve six-, and seven-membered rings. The structures used for 4 TABLE II. Tight-binding parameters describing the renor- malization of charge doping and energy dispersion at graphene (a) (b) GBs, obtained from fits to first-principles simulations of car- bon haeckelite structures, as illustrated in Figs. 1(d-j). ∆5 ∆7 ∆t51 ∆t71 ∆t52 ∆t72 0.4988 eV -0.497 eV 0.0005 0.0414 -0.04 -0.4095 (c) the description of the low-energy bands of the haeckelite structures, as shown in Fig. 1(d), they are also instrumen- tal in obtaining a qualitative and quantitative description of the spin splitting induced by SOC in the presence a transverse electric field, as shown in Fig. 1(e). While this model has been developed based on sim- (d) (e) ulations of haeckelite structures, it must also be trans- ferrable and applicable to the variety of GB structures that can be present in polycrystalline graphene. We therefore test our model on two prototypical GBs, as- sociated with the (2,1)|(1,2) and (5,0)|(3,3) interfaces. These correspond to a conducting type I GB and a type II GB with a transport gap, respectively, as classified (f) by Yazyev17 . The optimized geometries as well as the computed band structures and the spin splitting of the low-energy bands are shown in Fig. 2. As the compu- tational cells are larger in these cases, we turned to nu- merical atomic orbitals (NAOs) to compute the reference band structures. The self-consistent calculations have been performed with the OpenMX package18–20 , relying on a double-ζ polarized basis set for the expansion of the FIG. 2. Electronic structure of prototypical periodic GBs. eigenstates. The first Brillouin zones of the (2,1)|(1,2) Ball-and-stick representations of the (2,1)|(1,2) and (5,0)|(3,3) and (5,0)|(3,3) GBs were sampled with k-point grids of GBs are shown in panels (a) and (d) respectively. Colors have 16×3 and 8×3, respectively. As shown in Fig. 2, the topo- been added to visualize the pentagons (red) and heptagons logical model enables an accurate description of the elec- (blue). Yellow and green are used to depict the hexagonal tronic structure of the two prototypical GBs. Addition- carbon rings of the two separate graphene domains. (b) Low- ally, the renormalization is crucial for obtaining a quali- energy band structure of the (2,1)|(1,2) GB. (c) Spin split- tative and a quantitative description of the SOC-induced ting of the two metallic bands of the (2,1)|(1,2) GB under a spin splitting of the low-energy bands, as the TB model transverse electric field of 4 V/nm, with linear (upper panel) of pristine graphene (Eq. (1)) underestimates this split- and logarithmic (lower panel) vertical axis scale. Panels (e,f) ting by orders of magnitude (see the dashed green curves are the same as (b,c) for the (5,0)|(3,3) GB. Gold curves are first-principles results, blue curves are the renormalized TB in Figs. 2(c,f)). Finally, we note that our model does not model, and dashed green curves are the TB model of pristine capture the effect of vacancies or imperfect bonding con- graphene without renormalization due to GBs. figurations, both of which may be present in CVD-grown graphene. However, these account for only ∼0.1% of the atoms in our MD-generated polycrystalline samples, and thus we do not expect them to have a significant im- ameters of 14.8, 20.9, and 29.7 nm. For this study we pact on our results. Furthermore, progress in fabrication have also created two additional samples with average techniques now allows for very clean large-area polycrys- grain diameters of 14.9 and 21 nm. The 20.9-nm sam- talline graphene that is free from cracks and other strong ple is shown in Fig. 3(a), where the magenta regions are lattice imperfections21 . the pristine graphene grains and the black areas depict With an accurate TB model in hand, we now turn the nonhexagonal rings in the structure. The inset is to simulations of charge and spin transport in polycrys- a zoom-in of one graphene GB, showing a disordered ar- talline graphene. The polycrystalline samples have been ray of pentagons and heptagons between two misoriented created via molecular dynamics simulations that mimic grains. the CVD growth process22,23 . We use three samples that To study charge and spin transport in the polycrys- were previously used in Ref. 23, with average grain di- talline samples, we employ a real-space wave packet prop- 5 agation method that has been used to study charge and seen when graphene is placed in contact with a transition spin transport in a wide variety of disordered systems24 . metal dichalcogenide or a topological insulator26,27 . At time t = 0 we define Typical results are shown in Figs. 3(b,c) for the 20.9- nm sample with the full GB model, λR |E| = 100 µeV, N 1 X iξn and EF = 0.2 eV, corresponding to a carrier density |ψ(0)i = √ e |ni , (9) N n=1 of n ≈ 5 × 1012 cm-2 . Figure 3(b) shows the electri- cal conductivity σ as a function of time. At short times where ξn ∈ [0, 2π) is a random phase associated with σ increases linearly, indicative of the ballistic regime of each atomic site |ni in the polycrystalline sample, |ni is transport, and saturates to a maximum value as scat- an N × 1 vector with a 1 at row n and zeros elsewhere, tering off the GBs forces transport into the diffusive and N is the total number of atoms in the sample. This regime. At longer times, quantum interference leads to random phase state is then spin-polarized along the z- localized behavior, as indicated by the decay of σ with axis according to time. In the inset we√plot σ as a function of the propa- gation length L ≡ 2 ∆X 2 . The blue solid line shows the numerical results and the red dashed line is the " # I |ψ↑ (0)i = N |ψ(0)i , (10) decay expected from weak localization (WL) theory28 , 0 σ(L) = σsc − (G0 /π) ln(L/le ), where σsc = maxt {σ(t)} is the semiclassical conductivity, G0 = 2e2 /h is the quan- where IN is the N × N identity matrix. We let this tum of conductance, and le = τp vF is the mean free path. wave packet evolve in time and we calculate its mean The reasonable agreement of WL theory to the numerical square displacement ∆X 2 (EF , t) and its out-of-plane spin results, without any fitting parameters, indicates that at polarization sz (EF , t) as a function of Fermi energy and longer times the charge transport in the polycrystalline time, samples is in the weakly localized regime. This trend D E appears to break down for L > 200 nm, which could ψ↑X (t) δ(EF − H̃) ψ↑X (t) be due to periodicity effects, as the size of the periodic ∆X 2 (EF , t) = , (11) hψ↑ (0)| δ(EF − H̃) |ψ↑ (0)i polycrystalline sample is 180 × 180 nm. ( ) The solid line in Fig. 3(c) shows that the out-of- hψ↑ (t)| δ(EF − H̃)σz |ψ↑ (t)i plane spin polarization sz decays with time, indicating sz (EF , t) = Re , (12) hψ↑ (0)| δ(EF − H̃) |ψ↑ (0)i spin relaxation induced by scattering off the graphene GBs. The dashed line shows the expected decay as- suming the D’yakonov-Perel’ (DP) mechanism29 , which E where ψ↑X (t) = [X̂, Û (t)] |ψ↑ (0)i, X̂ is the position op- is typically the dominant mechanism of spin relaxation erator, Û (t) = exp(−iH̃t/~) is the time evolution oper- in two-dimensional electron systems with Rashba SOC30 . ator, |ψ↑ (t)i = Û (t) |ψ↑ (0)i, and σz is the Pauli matrix This decay is given by sz (t) = exp(−t/τs ), where τs = for spin along the z-axis. The time evolution operator 1/(Ω2R τp ) is the spin lifetime, ΩR = 2λR |E|/~ is the spin and the energy projection operator δ(EF − H̃) are both precession frequency, and τp is determined numerically expanded in a numerically efficient way using Chebyshev from Eq. (15). At short times the decay of sz follows the polynomials. From the mean square displacement we cal- DP relation, but at longer times there is a clear slowing culate the diffusion coefficient D, the electrical conduc- down of spin relaxation. This coincides with the onset tivity σ, and the momentum relaxation time τp , of localization, the result being that one cannot define a unique spin lifetime in this system. 1 d To obtain a global picture of spin transport and relax- D(EF , t) = ∆X 2 (EF , t), (13) 2 dt ation in polycrystalline graphene, in Fig. 4 we plot all σ(EF , t) = e2 ρ(EF )D(EF , t), (14) our simulation results. Figure 4(a) shows the conduc- 2Dmax (EF ) tivity as a function of time, where the dashed lines are τp (EF ) = , (15) for the pristine graphene TB model of Eq. (1) and the vF2 solid lines are for the renormalized GB model of Eq. (8). Comparing the solid and dashed lines reveals that the where ρ(EF ) = 2 hψ↑ (0)| δ(EF − H̃) |ψ↑ (0)i is the den- redistribution of the charge around the GBs leads to en- sity of states, vF is the Fermi velocity of graphene, and hanced scattering, reducing τp by an average of ∼10%. Dmax (EF ) = maxt {D(EF , t)}. Despite these differences in scattering strength, localiza- We simulated charge and spin transport in the five tion is clearly present in all cases. For each sample the different polycrystalline samples mentioned above, us- six curves corresponding to different values of λR |E| lie ing both the pristine graphene TB model of Eq. (1) and nearly on top of each other, as the presence of SOC does the renormalized GB model of Eq. (8), for Rashba SOC not significantly affect the charge transport. strengths of λR |E| = 12.5 µeV, 25 µeV, 50 µeV, 100 µeV, Figure 4(b) shows the out-of-plane spin polarization and 1 meV. The values in the lower µeV range are typical sz as a function of time. In this figure a few trends of those seen in graphene on SiO2 or hBN substrates25 , can be observed. First, the rate of spin relaxation in- while the range [100 µeV, 1 meV] corresponds to what is 6 relaxation is not constant and in no case is it possible to identify a unique spin lifetime in any of the simulated systems. However, it is possible to define a unique spin transport length in these polycrystalline graphene systems. The theory of DP spin relaxation was developed for trans- port in the diffusive regime, and in this regime the mean square displacement of charge carriers grows linearly in time, ∆X 2 = 2Dt. Using the definition of τp in Eq. (15) and plugging both into the theory of DP spin relaxation gives the spin polarization versus the propagation length  2 ! Ω R sz (L) = exp −∆X 2 · . (16) vF This expression indicates that the spin propagation is independent of disorder. This is a well-known conse- quence of DP spin relaxation and can be rationalized as follows30 . In the diffusive regime of √ transport the spin dif- fusion length is defined as Ls ≡ Dτs . Because D ∝ τp and τs ∝ 1/τp , we have a disorder-independent spin dif- fusion length Ls = vF /ΩR . Our simulations are well-described by this scaling be- havior, as shown in Fig. 4(c), where the spin polarization sz is plotted as a function of the mean square displace- ment ∆X 2 . The symbols are the numerical results, and one can see that they all collapse into a single universal decay that depends only on the strength of the Rashba SOC. The solid lines correspond to Eq. (16), indicating that in our simulations the spin relaxation is dominated entirely by the DP mechanism, including for longer times when weak localization is prevalent. To summarize, we have used first-principles calcula- tions to derive a tight-binding model for graphene that includes grain boundaries in the presence of spin-orbit coupling. We found that both the electronic and the spin properties of the GBs are captured by accounting for the redistribution of charge around the pentagons FIG. 3. (a) Structure of the polycrystalline graphene sam- and heptagons that lie at the interface between misori- ple with 20.9 nm average grain diameter. The inset shows ented graphene grains. We then used this model to study a zoom-in of one of the graphene GBs. (b) Time-dependent charge and spin transport in realistic models of polycrys- conductivity for a Rashba SOC strength λR |E| = 100 µeV at talline graphene. In the presence of Rashba and intrinsic EF = 0.2 eV, corresponding to a carrier density of n ≈ 5×1012 cm-2 . The inset shows the conductivity as a function of prop- SOC, we found that spin relaxation is determined en- agation length. The blue solid line is the numerical result and tirely by the D’yakonov-Perel’ mechanism. Within this the red dashed line is obtained from weak localization theory. mechanism the spin diffusion length is independent of (c) Corresponding out-of-plane spin polarization. The solid disorder and depends only on the strength of the Rashba line is the numerical result and the dashed line is the expected SOC. Our numerical simulations confirm this behavior, decay assuming DP spin relaxation in the diffusive regime of both in the diffusive and in the weakly localized regimes transport. of transport. One interesting result of our simulations is that the disorder-independence of the spin diffusion length ex- creases with Rashba SOC strength and with grain size tends to the weakly localized regime of transport. This (and thus with τp ), which is qualitatively similar to the effect has been studied analytically, and it was found that behavior expected from DP spin relaxation. Second, the in the presence of weak localization the diffusion coeffi- spin relaxation is faster for the pristine graphene model cient and the spin relaxation rate are both√renormalized compared to the full GB model. This is also consistent by the same factor, meaning that Ls ≡ Dτs remains with DP theory, as weaker scattering leads to faster spin constant31 . It has also been suggested that the indepen- relaxation. However, as discussed above, the rate of spin dence of Ls on disorder should extend to the strongly 7 mechanism33 . In this mechanism each scattering event has a finite probability to flip the spin, and the spin life- time is given by τs = (EF /(2λR |E|))2 τp . The dashed lines in Fig. 4(c) show the expected spin relaxation in the presence of both DP and EY mechanisms assum- ing τp = 6.5 fs, corresponding to the average τp of the polycrystalline samples. This clearly gives a faster spin relaxation than the numerical results, indicating the EY spin relaxation is either nonexistent in our simulations or much weaker than predicted. In the original prediction of EY spin relaxation in graphene, the quadratic scaling of τs with EF was a consequence of the Dirac-like linear dispersion of the electrons. However, this linear disper- sion relation does not necessarily describe the transport of electrons in graphene GBs, and thus the above expres- sion for EY spin relaxation may not be appropriate for polycrystalline graphene. The most consequential result of this work is that in the presence of Rashba and intrinsic SOC, the spin dif- fusion length in polycrystalline graphene is independent of grain size, Ls = vF /ΩR . This means that for a Rashba SOC strength of 50 µeV (5 µeV), the spin diffusion length will be Ls ≈ 5 µm (50 µm) whether the graphene grains are 10 nm or 10 µm in diameter. Thus, for spintron- ics applications, single-domain graphene may not be a necessity and focus can be placed on eliminating other sources of spin relaxation, such as magnetic impurities34 . In general, these results bode well for scalable CVD- grown graphene as an efficient transporter of spins in future spintronic applications. Finally, we point out that this study only considers uniform SOC, but other sources of spin relaxation could be present in CVD graphene. For example, out-of-plane corrugation induces local variations in SOC which can FIG. 4. (a) Time-dependent conductivity for all simulations relax spin. Theoretical studies of single-crystal graphene run in this work. Each set of colors corresponds to increasing concluded that the impact of corrugations is weak, with grain sizes of 14.8, 14.9, 20.9, 21, and 29.7 nm. The dashed Ls > 100 µm35 . Meanwhile, a variety of AFM stud- lines are for the pristine graphene TB model of Eq. (1) and ies were unable to identify any height variation arising solid lines are for the renormalized GB model of Eq. (8). (b) Time-dependent spin polarization for all simulations, where from graphene GBs4,36–38 , and if present it is limited each color corresponds to increasing Rashba SOC strength of to <0.3 Å39 . It is unclear whether localized corruga- λR |E| = 12.5 µeV, 25 µeV, 50 µeV, 100 µeV, and 1 meV. The tions of this height would contribute to spin relaxation dashed/solid lines have the same meaning as in panel (a). (c) in CVD graphene, and this would therefore be an intrigu- Spin polarization as a function of mean square displacement. ing future direction of research. Measurements of CVD Each color corresponds to a different Rashba SOC strength. graphene suggest that spin relaxation at GBs is not very Open squares (crosses) are numerical results for the pristine strong13,14 , but in the worst-case scenario we would ex- (renormalized GB) model of graphene. Solid lines correspond pect corrugations at the GBs to limit Ls to the graphene to Eq. (16), and dashed lines include the effect of Elliott-Yafet grain size, Ls ∼ 10 µm. Beyond corrugations, historically spin relaxation. The inset is a zoom to better show the case it is believed that spin transport in graphene has been of λR |E| = 1 meV. limited by paramagnetic impurities and contact-induced dephasing. In early measurements these effects limited spin lifetimes to hundreds of picoseconds and Ls ∼ 1 localized or insulating regime of transport32 . Our simu- µm40 . Improvements in these areas now permit the mea- lations are a direct confirmation of the universality of Ls surement of Ls > 10 µm, and as high as 30 µm in the in the diffusive and weakly localized regimes, but further best nonlocal measurements to date12 . Altogether, the numerical simulations with stronger disorder would be spin diffusion length arising from DP spin relaxation is needed to confirm this behavior in the insulating regime. therefore on par with these other sources of spin relax- In addition to the DP mechanism, it has been pre- ation, and may even be the dominant mechanism in the dicted that in the presence of Rashba SOC the out-of- highest-quality CVD samples currently being produced. plane spin can also be relaxed by the Elliott-Yafet (EY) 8 ACKNOWLEDGMENTS were provided by the supercomputing facilities of the Université catholique de Louvain (CISM/UCL) and the ICN2 is supported by the Severo Ochoa program “Consortium des Équipements de Calcul Intensif” en from Spanish MINECO (grant no. SEV-2017-0706) and Fédération Wallonie-Bruxelles (CECI). All authors ac- funded by the CERCA Programme / Generalitat de knowledge support from the European Union Horizon Catalunya. S.M.-M.D and J.-C.C. acknowledge the Na- 2020 Programme under grant agreement nos. 696656 and tional Fund for Scientific Research of Belgium [F.R.S.- 785219 (Graphene Flagship Core 1 and Core 2). FNRS] for financial support. Computational resources ∗ aron.cummings@icn2.cat C. Stampfer, and B. Beschoten, “Spin Lifetimes Exceed- 1 Y. Zhang, L. Zhang, and C. Zhou, “Review of Chemical ing 12 ns in Graphene Nonlocal Spin Valve Devices,” Nano Vapor Deposition of Graphene and Related Applications,” Lett. 16, 3533–3539 (2016). 13 Accounts Chem. Res. 46, 2329–2339 (2013). M. V. Kamalakar, C. Groenveld, A. Dankert, and 2 P. Y. Huang, C. S. Ruiz-Vargas, A. M. van der Zande, S. P. Dash, “Long distance spin communication in chem- W. S. Whitney, M. P. Levendorf, J. W. Kevek, S. Garg, ical vapour deposited graphene,” Nat. Commun. 6, 6766 J. S. Alden, C. J. Hustedt, Y. Zhu, J. Park, P. L. McEuen, (2015). 14 and D. A. Muller, “Grains and grain boundaries in single- Z. M. Gebeyehu, S. Parui, J. F. Sierra, M. Timmermans, layer graphene atomic patchwork quilts,” Nature 469, 389– M. J. Esplandiu, S. Brems, C. Huyghebaert, K. Garello, 392 (2011). M. V. Costache, and S. O. Valenzuela, “Spin communica- 3 K. Kim, Z. Lee, W. Regan, C. Kisielowski, M. F. Crommie, tion over 30 µm long channels of chemical vapor deposited and A. Zettl, “Grain Boundary Mapping in Polycrystalline graphene on SiO2 ,” 2D Mater. 6, 034003 (2019). 15 Graphene,” ACS Nano 5, 2142–2146 (2011). C. L. Kane and E. J. Mele, “Quantum Spin Hall Effect in 4 Q. Yu, L. A. Jauregui, W. Wu, R. Colby, J. Tian, Z. Su, Graphene,” Phys. Rev. Lett. 95, 226801 (2005). 16 H. Cao, Z. Liu, D. Pandey, D. Wei, T. F. Chung, P. Peng, H. Terrones, M. Terrones, E. Hernández, N. Grobert, J.-C. N. P. Guisinger, E. A. Stach, J. Bao, S.-S. Pei, and Charlier, and P. M. Ajayan, “New Metallic Allotropes of Y. P. Chen, “Control and characterization of individual Planar and Tubular Carbon,” Phys. Rev. Lett. 84, 1716– grains and grain boundaries in graphene grown by chemi- 1719 (2000). 17 cal vapour deposition,” Nat. Mater. 10, 443–449 (2011). O. V. Yazyev and S. G. Louie, “Electronic transport in 5 A. W. Tsen, L. Brown, M. P. Levendorf, F. Ghahari, P. Y. polycrystalline graphene,” Nat. Mater. 9, 806 (2010). 18 Huang, R. W. Havener, C. S. Ruiz-Vargas, D. A. Muller, T. Ozaki, “Variationally optimized atomic orbitals for P. Kim, and J. Park, “Tailoring Electrical Transport large-scale electronic structures,” Phys. Rev. B 67, 155108 Across Grain Boundaries in Polycrystalline Graphene,” (2003). 19 Science 336, 1143–1146 (2012). T. Ozaki and H. Kino, “Numerical atomic basis orbitals 6 L. Tapasztó, P. Nemes-Incze, G. Dobrik, K. Jae Yoo, from H to Kr,” Phys. Rev. B 69, 195113 (2004). 20 C. Hwang, and L. P. Biró, “Mapping the electronic proper- T. Ozaki and H. Kino, “Efficient projector expansion for ties of individual graphene grain boundaries,” Appl. Phys. the ab initio LCAO method,” Phys. Rev. B 72, 045121 Lett. 100, 053114 (2012). (2005). 7 21 J. C. Koepke, J. D. Wood, D. Estrada, Z.-Y. Ong, K. T. J. Zhang, K. Jia, L. Lin, W. Zhao, H. T. Quang, L. Sun, He, E. Pop, and J. W. Lyding, “Atomic-Scale Evidence T. Li, Z. Li, X. Liu, L. Zheng, R. Xue, J. Gao, Z. Luo, M. H. for Potential Barriers and Strong Carrier Scattering at Rummeli, Q. Yuan, H. Peng, and Z. Liu, “Large-Area Graphene Grain Boundaries: A Scanning Tunneling Mi- Synthesis of Superclean Graphene via Selective Etching of croscopy Study,” ACS Nano 7, 75–86 (2013). Amorphous Carbon with Carbon Dioxide,” Angew. Chem. 8 A. Isacsson, A. W. Cummings, L. Colombo, L. Colombo, Int. Ed. , doi:10.1002/anie.201905672. 22 J. M. Kinaret, and S. Roche, “Scaling properties of poly- J. Kotakoski and J. C. Meyer, “Mechanical properties crystalline graphene: a review,” 2D Mater. 4, 012002 of polycrystalline graphene based on a realistic atomistic (2016). model,” Phys. Rev. B 85, 195447 (2012). 9 23 W. Han, R. K. Kawakami, M. Gmitra, and J. Fabian, D. Van Tuan, J. Kotakoski, T. Louvet, F. Ortmann, J. C. “Graphene spintronics,” Nat. Nanotechnol. 9, 79480 Meyer, and S. Roche, “Scaling Properties of Charge Trans- (2014). port in Polycrystalline Graphene,” Nano Lett. 13, 1730– 10 S. Roche, J. Åkerman, B. Beschoten, J.-C. Charlier, 1735 (2013). 24 M. Chshiev, S. P. Dash, B. Dlubak, J. Fabian, A. Fert, Z. Fan, J. H. Garcia, A. W. Cummings, J. E. Barrios, M. Guimarães, F. Guinea, I. Grigorieva, C. Schönenberger, M. Panhans, A. Harju, F. Ortmann, and S. Roche, P. Seneor, C. Stampfer, S. O. Valenzuela, X. Waintal, and “Linear Scaling Quantum Transport Methodologies,” B. van Wees, “Graphene spintronics: the European Flag- arXiv:1811.07387 [cond-mat.mes-hall]. 25 ship perspective,” 2D Mater. 2, 030202 (2015). K. Zollner, M. Gmitra, and J. Fabian, “Heterostructures 11 X. Lin, E. Yang, K. L. Wang, and W. Zhao, “Two- of graphene and hBN: Electronic, spin-orbit, and spin re- dimensional spintronics for low-power electronics,” Nat. laxation properties from first principles,” Phys. Rev. B 99, Electron. 2, 274–283 (2019). 125151 (2019). 12 26 M. Drögeler, C. Franzen, F. Volmer, T. Pohlmann, M. Gmitra, D. Kochan, P. Högl, and J. Fabian, “Trivial L. Banszerus, M. Wolter, K. Watanabe, T. Taniguchi, and inverted Dirac bands and the emergence of quantum 9 spin Hall states in graphene on transition-metal dichalco- Mechanism in Graphene: Resonant Scattering by Mag- genides,” Phys. Rev. B 93, 155104 (2016). netic Impurities,” Phys. Rev. Lett. 112, 116602 (2014). 27 35 K. Song, D. Soriano, A. W. Cummings, R. Robles, I. M. Vicent, H. Ochoa, and F. Guinea, “Spin relaxation P. Ordejón, and S. Roche, “Spin Proximity Effects in in corrugated graphene,” Phys. Rev. B 95, 195402 (2017). 36 Graphene/Topological Insulator Heterostructures,” Nano J. Červenka, M. I. Katsnelson, and C. F. J. Flipse, “Room- Lett. 18, 2033–2039 (2018). temperature ferromagnetism in graphite driven by two- 28 P. A. Lee and T. V. Ramakrishnan, “Disordered electronic dimensional networks of point defects,” Nat. Phys. 5, 840 systems,” Rev. Mod. Phys. 57, 287–337 (1985). (2009). 29 37 M. I. Dyakonov and V. I. Perel, “Spin relaxation of con- Z. Fei, A. S. Rodin, W. Gannett, S. Dai, W. Regan, duction electrons in noncentrosymmetric semiconductors,” M. Wagner, M. K. Liu, A. S. McLeod, G. Dominguez, Sov. Phys. Solid State 13, 3023–3026 (1972). M. Thiemens, A. H. Castro Neto, F. Keilmann, A. Zettl, 30 I. Žutić, J. Fabian, and S. Das Sarma, “Spintronics: Fun- R. Hillenbrand, M. M. Fogler, and D. N. Basov, “Elec- damentals and applications,” Rev. Mod. Phys. 76, 323–410 tronic and plasmonic phenomena at graphene grain bound- (2004). aries,” Nat. Nanotechnol. 8, 821 (2013). 31 38 A. G. Mal‘shukov, K. A. Chao, and M. Willander, “Quan- Y. Ogawa, K. Komatsu, K. Kawahara, M. Tsuji, K. Tsuk- tum Localization Effects on Spin Transport in Semiconduc- agoshi, and H. Ago, “Structure and transport properties tor Quantum Wells with Zinc-Blende Crystal Structure,” of the interface between CVD-grown graphene domains,” Phys. Rev. Lett. 76, 3794–3797 (1996). Nanoscale 6, 7288–7294 (2014). 32 39 B. I. Shklovskii, “Dyakonov-Perel spin relaxation near J. Červenka and C. F. J. Flipse, “Structural and elec- the metal-insulator transition and in hopping transport,” tronic properties of grain boundaries in graphite: Planes Phys. Rev. B 73, 193201 (2006). of periodically distributed point defects,” Phys. Rev. B 79, 33 H. Ochoa, A. H. Castro Neto, and F. Guinea, “Elliot-Yafet 195429 (2009). 40 Mechanism in Graphene,” Phys. Rev. Lett. 108, 206808 N. Tombros, C. Jozsa, M. Popinciuc, H. T. Jonkman, and (2012). B. J. van Wees, “Electronic spin transport and spin pre- 34 D. Kochan, M. Gmitra, and J. Fabian, “Spin Relaxation cession in single graphene layers at room temperature,” Nature 448, 571 (2007).