Universal spin diffusion length in polycrystalline graphene

TL;DR

This study shows that in polycrystalline graphene, the spin diffusion length is unaffected by grain size and is solely determined by substrate-induced spin-orbit coupling, which is promising for large-scale spintronic applications.

Contribution

The paper introduces a first-principles derived tight-binding model to evaluate spin transport in polycrystalline graphene with grain boundaries, revealing the independence of spin diffusion length from grain size.

Findings

Spin diffusion length is independent of grain size.

Spin relaxation is governed by substrate-induced spin-orbit coupling.

Results support the D'yakonov-Perel' mechanism in the diffusive regime.

Abstract

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 spintronics; despite being polycrystalline, spin transport in CVD graphene has been measured over lengths up to 30 $\mu$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 impact 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.

Tables

Table 1. Table 1: Tight-binding parameters for pristine graphene obtained from fits to first-principles calculations, as illustrated in Figs.  1 (a-c).

$t_1^o$	${\beta }_1$	$t_2^o$	${\beta }_2$	${\lambda }_{\text{I}}$	${\lambda }_{\text{R}}$
-2.414 eV	1.847	-0.168 eV	3.077	13.437 $\mu$eV	3.383 $\frac{\mu \text{eV}}{\text{V/nm}}$

Table 2. Table 2: Tight-binding parameters describing the renormalization of charge doping and energy dispersion at graphene GBs, obtained from fits to first-principles simulations of carbon haeckelite structures, as illustrated in Figs.  1 (d-j).

$\Delta {ϵ}^5$	$\Delta {ϵ}^7$	$\Delta t_1^5$	$\Delta t_1^7$	$\Delta t_2^5$	$\Delta t_2^7$
0.4988 eV	-0.497 eV	0.0005	0.0414	-0.04	-0.4095

Equations

$$\hat{\mathcal{H}}$$

$${\hat{\mathcal{H}}}_{\text{kin}}$$

$${\hat{\mathcal{H}}}_{\text{SOC}}$$

$$+\frac{2\text{i}}{3}{\lambda }_{\text{R}}\sum _{⟨i,j⟩}{\hat{c}}_{i\sigma }^{†}{\left[\mathbf{E}\cdot \left(\hat{\mathbf{s}}\times {\mathbf{d}}_{ij}\right)\right]}_{\sigma {\sigma }^{\prime }}{\hat{c}}_{j{\sigma }^{\prime }},$$

$$t_1=t_1^o{e}^{-{\beta }_1(r-a)}\text{and}{t}_2=t_2^o{e}^{-{\beta }_2(r-b)}.$$

$${\tilde{\mathcal{H}}}_{\text{loc}}=\sum _i{\hat{c}}_{i\sigma }^{†}{ϵ}_i{\hat{c}}_{i\sigma }\text{ with }{ϵ}_i=\sum _{m\in \{{\mathcal{L}}_i\}}\Delta {ϵ}^m.$$

$${\tilde{\mathcal{H}}}_{\text{kin}}=\sum _{⟨i,j⟩}{\hat{c}}_{i\sigma }^{†}{t}_{1ij}{\hat{c}}_{j\sigma }+\sum _{⟨⟨i,j⟩⟩}{\hat{c}}_{i\sigma }^{†}{t}_{2ij}{\hat{c}}_{j\sigma },$$

$$t_{1ij}$$

$$t_{2ij}$$

$$\tilde{\mathcal{H}}={\tilde{\mathcal{H}}}_{\text{loc}}+{\tilde{\mathcal{H}}}_{\text{kin}}+{\hat{\mathcal{H}}}_{\text{SOC}}.$$

$$|\psi (0)⟩=\frac{1}{\sqrt{N}}\sum _{n=1}^N{\text{e}}^{\text{i}{\xi }_n}|n⟩,$$

$$|{\psi }_{\uparrow }(0)⟩=\left[\begin{array}{c}{\mathcal{I}}_N\\ 0\end{array}\right]|\psi (0)⟩,$$

$$\Delta X^2(E_{\text{F}},t)=\frac{⟨{\psi }_{\uparrow }^X(t)|\delta (E_{\text{F}}-\tilde{\mathcal{H}})|{\psi }_{\uparrow }^X(t)⟩}{⟨{\psi }_{\uparrow }(0)|\delta (E_{\text{F}}-\tilde{\mathcal{H}})|{\psi }_{\uparrow }(0)⟩},$$

$$s_z(E_{\text{F}},t)=\mathrm{Re}\left\{\frac{⟨{\psi }_{\uparrow }(t)|\delta (E_{\text{F}}-\tilde{\mathcal{H}}){\sigma }_z|{\psi }_{\uparrow }(t)⟩}{⟨{\psi }_{\uparrow }(0)|\delta (E_{\text{F}}-\tilde{\mathcal{H}})|{\psi }_{\uparrow }(0)⟩}\right\},$$

$$D(E_{\text{F}},t)$$

$$\sigma (E_{\text{F}},t)$$

$${\tau }_{\text{p}}(E_{\text{F}})$$

$$s_z\left(L\right)=\exp \left(-\Delta X^2\cdot {\left(\frac{{\Omega }_{\text{R}}}{v_{\text{F}}}\right)}^2\right).$$

Full text

graphene; CVD; polycrystalline; grain boundaries; spin relaxation; spintronics

The growth of graphene via chemical vapor deposition (CVD) is the most promising approach for realizing industrial-scale applications of this material Zhang et al. (2013, ). One drawback of CVD-grown graphene is that it tends to be polycrystalline, with misoriented single-crystal domains separated by grain boundaries (GBs) consisting of arrays of five-, seven-, and occasionally eight-member carbon rings. In some specific cases the GBs can be characterized by a given periodicity, but more generally they tend to be complex meandering arrangements of these nonhexagonal rings Huang et al. (2011); Kim et al. (2011). Charge transport and scanning tunneling measurements have revealed that graphene GBs serve as a significant source of charge scattering, leading to enhanced resistance Yu et al. (2011); Tsen et al. (2012) and localization effects Yu et al. (2011); Tapasztó et al. (2012); Koepke et al. (2013). Only when the graphene grains become larger than 1-10 $\upmu$m do the GBs cease to dominate the charge transport in CVD graphene Isacsson et al. (2016).

Graphene also has clear advantages for spintronic applications, owing to its low intrinsic spin-orbit coupling (SOC) Han et al. (2014); Roche et al. (2015); Lin et al. (2019). Combined with its high electron mobility, this can lead to spin diffusion lengths as long as 30 $\upmu$m in clean exfoliated graphene-based devices (with mobility up to several $10,000$ cm2/V-s) Drögeler et al. (2016). Intriguingly, spin diffusion lengths as long as 10 $\upmu$m and spin currents measurable over channel lengths of 30 $\upmu$m have also been reported in disordered CVD graphene with much lower charge mobility Kamalakar et al. (2015); Gebeyehu et al. (2019). These measurements can be explained in two different ways: either none or very few GBs were present in the measured devices; or GBs in CVD graphene, while impeding charge transport, have little effect on spin transport and relaxation. Although the former hypothesis remains plausible, the lack of a theoretical foundation concerning spin transport in the presence of disordered networks of GBs leaves the second hypothesis open as a possibility, which could have profound consequences for the optimization of graphene-based spintronic devices.

Here we use numerical simulations to study the impact of GBs on spin transport in polycrystalline graphene. We first develop an effective tight-binding (TB) model of polycrystalline graphene in the presence of intrinsic and substrate-induced SOC, which is derived from extensive first-principles calculations. The model is based on simulations of a variety of carbon-based haeckelites and is found to be general and transferable to the complex morphologies of graphene GBs. We then use this model to perform spin transport simulations in realistic models of polycrystalline graphene, using an efficient linear-scaling methodology that gives direct access to spin relaxation and propagation. Our simulations reveal that the spin diffusion length in polycrystalline graphene is independent of grain size and depends only on the strength of the substrate-induced SOC. This result is fully consistent with the D’yakonov-Perel’ mechanism of spin relaxation in the diffusive regime, but here it is also shown to be robust to the contributions of quantum interference induced by disorder. These findings indicate that in the presence of SOC the graphene GBs serve as scatterers of charge, but they do not play a direct role in spin relaxation. In other words, our results suggest that grain size is not a limiting factor for spin transport in CVD-grown graphene, a result of genuine relevance for the future optimization of graphene-based spin devices and architectures in the context of memory or spin logic technologies.

We first develop an effective TB model for the description of itinerant electrons in polycrystalline graphene, which is derived from fitting to first-principles simulations. The model includes the impact of electrostatic barriers and atomic disorder present at GBs, and explicitly accounts for both intrinsic and extrinsic SOC. The first-principles simulations were carried out with the all-electron FP-LAPW method as implemented in the 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 Bohr*-1*. A $33\times 33$ k-point mesh was used to sample the first Brillouin zone of pristine graphene, and equivalent k-point densities were used for the supercell calculations.

The starting point of our one-orbital TB model is the now well-known functional form of the second-nearest-neighbor hopping Hamiltonian of graphene. The Hamiltonian can be conveniently decomposed into a kinetic operator and a SOC operator,

[TABLE]

where $\hat{c}_{i\sigma}^{\dagger}$ ($\hat{c}_{i\sigma}$) is the creation (annihilation) operator for the $p_{z}$ orbital with spin $\sigma$ at lattice site $i$, $\mathbf{s}$ is the spin Pauli matrix, $\mathbf{E}$ is the external electric field, $\mathbf{d}_{ij}$ is the unit vector pointing from site $j$ to $i$ and $\nu_{ij}$ is $+1$ ($-1$) for clockwise (counter-clockwise) hopping paths from site $j$ to $i$. Single (double) brackets stand for summation over the nearest (second-nearest) neighboring lattice sites. In our model we assume an exponential dependence of the $t_{1}$ and $t_{2}$ hopping parameters with respect to the intersite distance,

[TABLE]

Here, $a$ ($b$) is the equilibrium distance between the first (second) nearest atomic sites in graphene at equilibrium ($a=1.42$ Å, $b=\sqrt{3}a=2.46$ Å). The third term is the intrinsic SOC operator which connects second nearest neighbor sites, and the last term is the extrinsic (Rashba) SOC operator that arises from an external perpendicular electric field $\mathbf{E}$ Kane and Mele (2005). This electric field can be directly applied via a gate voltage, or it can be an effective field that arises from placing graphene on a substrate. The Rashba term connects nearest neighbor sites and linearly depends on the strength of the perpendicular electric field. The hopping integrals $t_{1}$ and $t_{2}$, as well as the parameters $\beta_{1}$ and $\beta_{2}$ ruling their distance dependence, were fitted together with the SOC parameters $\lambda_{\text{I}}$ and $\lambda_{\text{R}}$ with respect to first-principles band structure calculations of graphene at equilibrium and under isotropic strains up to 20$\%$. As our focus is on the transport properties of low-energy itinerant electrons, the fitting has been limited to the $[-1,1]$ eV energy window around the Fermi level. The optimized parameters are given in Table 1. The radial dependence of the hopping integral is shown in Fig. 1(a). The perfect agreement of the model with reference band structure calculations is illustrated in Figs. 1(b,c).

To account for the crystallographic disorder introduced by dislocations at GBs, the TB model is enriched by introducing topology-dependent renormalization factors. Within a simple ball-and-stick model, any $sp^{2}$ lattice can be mapped to a tiling of the plane by polygons, where each site is the shared apex of three polygons. In our case, it can be three hexagons as in graphene, or any combination of pentagons, hexagons and heptagons due to the presence of disclinations in the lattice. For the purpose of our model, we denote $\mathcal{L}_{i}$ as the geometric environment of a given site $i$. $\mathcal{L}_{i}$ can be any set of three elements out of the ensemble of considered carbon polygons $\mathcal{E}=\{5,6,7\}$. Similarly, $\mathcal{L}_{ij}$ are defined as the local environment of a pair of neighboring sites $i$ and $j$. If the sites are first nearest neighbors, the vector joining $i$ and $j$ is the shared edge of two polygons and $\mathcal{L}_{ij}$ can be any combination of two elements of the ensemble $\mathcal{E}$. If the sites are second nearest neighbors, the vector joining the sites is inscribed within a polygon and $\mathcal{L}_{ij}$ is reduced to a single element of $\mathcal{E}$. Finally, an extra term is added to the on-site energy of the Hamiltonian, which describes the local redistribution of charge around pentagon and heptagon carbon rings,

[TABLE]

This term is crucial for the description of the resonance peaks introduced by GBs in the low-energy electronic spectrum. It also enables a geometry-dependent electrostatic alignment of the GBs with respect to the graphene grains. Note that without loss of generality we can impose the constraint $\Delta\epsilon^{6}=0$, which conveniently defines the reference energy as the charge neutrality point of pristine graphene.

To account for the local variation of aromaticity and energy-momentum dispersion relation in the presence of GBs, the kinetic operator is also renormalized,

[TABLE]

with

[TABLE]

We impose $\Delta t_{1}^{6}=\Delta t_{2}^{6}=0$ in order to preserve the energy-momentum dispersion relation of pristine graphene. The full Hamiltonian of our model now reads

[TABLE]

The additional parameters of our topological model, $\{\Delta\epsilon^{5},\ \Delta\epsilon^{7},\ \Delta t_{1}^{5},\ \Delta t_{1}^{7},\ \Delta t_{2}^{5},\ \Delta t_{2}^{7}\}$, have been fitted against a set of first-principles band structures corresponding to periodic carbon $sp^{2}$ systems containing five-, six-, and seven-membered rings. The structures used for the fitting of our model, known as haeckelites Terrones et al. (2000), are illustrated in Figs. 1(f-j). These haeckelite structures have been created via the incorporation of disclinations into $2\times 2$ and $3\times 3$ graphene supercells. The considered geometries have been fully relaxed with respect to both atomic and cell degrees of freedom. The topological parameters obtained by fitting the low-energy band structures are given in Table 2.

We note that the renormalization of the nearest neighbor hopping integral is very weak. This directly reflects the $sp^{2}$ character of the fitting set, and additional benchmarks revealed that these renormalization factors contribute only marginally to the improved description of $sp^{2}$ carbon structures. The main contribution to the success of the model comes from the geometry-dependent renormalization of the onsite energies. Figures 1(d,e) illustrate the improved description of itinerant electrons due to the renormalization. While $\Delta\epsilon^{5}$ and $\Delta\epsilon^{7}$ improve the description of the low-energy bands of the haeckelite structures, as shown in Fig. 1(d), they are also instrumental 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 simulations of haeckelite structures, it must also be transferrable and applicable to the variety of GB structures that can be present in polycrystalline graphene. We therefore test our model on two prototypical GBs, associated 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 by Yazyev Yazyev and Louie (2010). 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 computational cells are larger in these cases, we turned to numerical atomic orbitals (NAOs) to compute the reference band structures. The self-consistent calculations have been performed with the OpenMX package Ozaki (2003); Ozaki and Kino (2004, 2005), relying on a double-$\zeta$ polarized basis set for the expansion of the eigenstates. The first Brillouin zones of the (2,1)$|$(1,2) and (5,0)$|$(3,3) GBs were sampled with k-point grids of $16\times 3$ and $8\times 3$, respectively. As shown in Fig. 2, the topological model enables an accurate description of the electronic structure of the two prototypical GBs. Additionally, the renormalization is crucial for obtaining a qualitative and a quantitative description of the SOC-induced spin splitting of the low-energy bands, as the TB model of pristine graphene (Eq. (1)) underestimates this splitting by orders of magnitude (see the dashed green curves in Figs. 2(c,f)). Finally, we note that our model does not capture the effect of vacancies or imperfect bonding configurations, both of which may be present in CVD-grown graphene. However, these account for only \sim$$0.1\% of the atoms in our MD-generated polycrystalline samples, and thus we do not expect them to have a significant impact on our results. Furthermore, progress in fabrication techniques now allows for very clean large-area polycrystalline graphene that is free from cracks and other strong lattice imperfections Zhang et al. .

With an accurate TB model in hand, we now turn to simulations of charge and spin transport in polycrystalline graphene. The polycrystalline samples have been created via molecular dynamics simulations that mimic the CVD growth process Kotakoski and Meyer (2012); Van Tuan et al. (2013). We use three samples that were previously used in Ref. 23, with average grain diameters of 14.8, 20.9, and 29.7 nm. For this study we have also created two additional samples with average grain diameters of 14.9 and 21 nm. The 20.9-nm sample is shown in Fig. 3(a), where the magenta regions are the pristine graphene grains and the black areas depict the nonhexagonal rings in the structure. The inset is a zoom-in of one graphene GB, showing a disordered array of pentagons and heptagons between two misoriented grains.

To study charge and spin transport in the polycrystalline samples, we employ a real-space wave packet propagation method that has been used to study charge and spin transport in a wide variety of disordered systems Fan et al. ; Roche (1999); Roche and Mayou (1997). At time $t=0$ we define

[TABLE]

where $\xi_{n}\in\left[0,2\pi\right)$ is a random phase associated with each atomic site $\left|n\right>$ in the polycrystalline sample, $\left|n\right>$ is an $N\times 1$ vector with a $1$ at row $n$ and zeros elsewhere, and $N$ is the total number of atoms in the sample. This random phase state is then spin-polarized along the $z$-axis according to

[TABLE]

where $\mathcal{I}_{N}$ is the $N\times N$ identity matrix. We let this wave packet evolve in time and we calculate its mean square displacement $\Delta X^{2}(E_{\text{F}},t)$ and its out-of-plane spin polarization $s_{z}(E_{\text{F}},t)$ as a function of Fermi energy and time,

[TABLE]

where $\left|\psi_{\uparrow}^{X}(t)\right>=[\hat{X},\hat{U}(t)]\left|\psi_{\uparrow}(0)\right>$, $\hat{X}$ is the position operator, $\hat{U}(t)=\exp(-\text{i}\mathcal{\tilde{H}}t/\hbar)$ is the time evolution operator, $\left|\psi_{\uparrow}(t)\right>=\hat{U}(t)\left|\psi_{\uparrow}(0)\right>$, and $\sigma_{z}$ is the Pauli matrix for spin along the $z$-axis. The time evolution operator and the energy projection operator $\delta(E_{\text{F}}-\mathcal{\tilde{H}})$ are both expanded in a numerically efficient way using Chebyshev polynomials. From the mean square displacement we calculate the diffusion coefficient $D$, the electrical conductivity $\sigma$, and the momentum relaxation time $\tau_{\text{p}}$,

[TABLE]

where $\rho(E_{\text{F}})=2\left<\psi_{\uparrow}(0)\right|\delta(E_{\text{F}}-\mathcal{\tilde{H}})\left|\psi_{\uparrow}(0)\right>$ is the density of states, $v_{\text{F}}$ is the Fermi velocity of graphene, and $D_{\text{max}}(E_{\text{F}})=\max_{t}\left\{D(E_{\text{F}},t)\right\}$.

We simulated charge and spin transport in the five different polycrystalline samples mentioned above, using both the pristine graphene TB model of Eq. (1) and the renormalized GB model of Eq. (8), for Rashba SOC strengths of $\lambda_{\text{R}}|\mathbf{E}|=$ 12.5 $\upmu$eV, 25 $\upmu$eV, 50 $\upmu$eV, 100 $\upmu$eV, and 1 meV. The values in the lower $\upmu$eV range are typical of those seen in graphene on SiO2 or hBN substrates Zollner et al. (2019), while the range $\left[100~{}\upmu\text{eV},1~{}\text{meV}\right]$ corresponds to what is seen when graphene is placed in contact with a transition metal dichalcogenide or a topological insulator Gmitra et al. (2016); Song et al. (2018).

Typical results are shown in Figs. 3(b,c) for the 20.9-nm sample with the full GB model, $\lambda_{\text{R}}|\mathbf{E}|=100$ $\upmu$eV, and $E_{\text{F}}=0.2$ eV, corresponding to a carrier density of $n\approx 5\times 10^{12}$ cm${}^{\text{-2}}$. Figure 3(b) shows the electrical conductivity $\sigma$ as a function of time. At short times $\sigma$ increases linearly, indicative of the ballistic regime of transport, and saturates to a maximum value as scattering off the GBs forces transport into the diffusive regime. At longer times, quantum interference leads to localized behavior, as indicated by the decay of $\sigma$ with time. In the inset we plot $\sigma$ as a function of the propagation length $L\equiv 2\sqrt{\Delta X^{2}}$. The blue solid line shows the numerical results and the red dashed line is the decay expected from weak localization (WL) theory Lee and Ramakrishnan (1985), $\sigma(L)=\sigma_{\text{sc}}-(G_{0}/\pi)\ln(L/l_{\text{e}})$, where $\sigma_{\text{sc}}=\max_{t}\left\{\sigma(t)\right\}$ is the semiclassical conductivity, $G_{0}=2e^{2}/h$ is the quantum of conductance, and $l_{\text{e}}=\tau_{\text{p}}v_{\text{F}}$ is the mean free path. The reasonable agreement of WL theory to the numerical results, without any fitting parameters, indicates that at longer times the charge transport in the polycrystalline samples is in the weakly localized regime. This trend appears to break down for $L>200$ nm, which could be due to periodicity effects, as the size of the periodic polycrystalline sample is $180\times 180$ nm.

The solid line in Fig. 3(c) shows that the out-of-plane spin polarization $s_{z}$ decays with time, indicating spin relaxation induced by scattering off the graphene GBs. The dashed line shows the expected decay assuming the D’yakonov-Perel’ (DP) mechanism Dyakonov and Perel (1972), which is typically the dominant mechanism of spin relaxation in two-dimensional electron systems with Rashba SOC Žutić et al. (2004). This decay is given by $s_{z}(t)=\exp(-t/\tau_{\text{s}})$, where $\tau_{\text{s}}=1/(\Omega_{\text{R}}^{2}\tau_{\text{p}})$ is the spin lifetime, $\Omega_{\text{R}}=2\lambda_{\text{R}}|\mathbf{E}|/\hbar$ is the spin precession frequency, and $\tau_{\text{p}}$ is determined numerically from Eq. (15). At short times the decay of $s_{z}$ follows the DP relation, but at longer times there is a clear slowing down of spin relaxation. This coincides with the onset of localization, the result being that one cannot define a unique spin lifetime in this system.

To obtain a global picture of spin transport and relaxation in polycrystalline graphene, in Fig. 4 we plot all our simulation results. Figure 4(a) shows the conductivity as a function of time, where the dashed lines are for the pristine graphene TB model of Eq. (1) and the solid lines are for the renormalized GB model of Eq. (8). Comparing the solid and dashed lines reveals that the redistribution of the charge around the GBs leads to enhanced scattering, reducing $\tau_{\text{p}}$ by an average of \sim$$10\%. Despite these differences in scattering strength, localization is clearly present in all cases. For each sample the six curves corresponding to different values of $\lambda_{\text{R}}|\mathbf{E}|$ lie nearly on top of each other, as the presence of SOC does not significantly affect the charge transport.

Figure 4(b) shows the out-of-plane spin polarization $s_{z}$ as a function of time. In this figure a few trends can be observed. First, the rate of spin relaxation increases with Rashba SOC strength and with grain size (and thus with $\tau_{\text{p}}$), which is qualitatively similar to the behavior expected from DP spin relaxation. Second, the spin relaxation is faster for the pristine graphene model compared to the full GB model. This is also consistent with DP theory, as weaker scattering leads to faster spin relaxation. However, as discussed above, the rate of spin 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 transport in the diffusive regime, and in this regime the mean square displacement of charge carriers grows linearly in time, $\Delta X^{2}=2Dt$. Using the definition of $\tau_{\text{p}}$ in Eq. (15) and plugging both into the theory of DP spin relaxation gives the spin polarization versus the propagation length

[TABLE]

This expression indicates that the spin propagation is independent of disorder. This is a well-known consequence of DP spin relaxation and can be rationalized as follows Žutić et al. (2004). In the diffusive regime of transport the spin diffusion length is defined as $L_{\text{s}}\equiv\sqrt{D\tau_{\text{s}}}$. Because $D\propto\tau_{\text{p}}$ and $\tau_{\text{s}}\propto 1/\tau_{\text{p}}$, we have a disorder-independent spin diffusion length $L_{\text{s}}=v_{\text{F}}/\Omega_{\text{R}}$.

Our simulations are well-described by this scaling behavior, as shown in Fig. 4(c), where the spin polarization $s_{z}$ is plotted as a function of the mean square displacement $\Delta 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 calculations 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 and heptagons that lie at the interface between misoriented graphene grains. We then used this model to study charge and spin transport in realistic models of polycrystalline graphene. In the presence of Rashba and intrinsic SOC, we found that spin relaxation is determined entirely by the D’yakonov-Perel’ mechanism. Within this mechanism the spin diffusion length is independent of disorder and depends only on the strength of the Rashba SOC. Our numerical simulations confirm this behavior, both in the diffusive and in the weakly localized regimes of transport.

One interesting result of our simulations is that the disorder-independence of the spin diffusion length extends to the weakly localized regime of transport. This effect has been studied analytically, and it was found that in the presence of weak localization the diffusion coefficient and the spin relaxation rate are both renormalized by the same factor, meaning that $L_{\text{s}}\equiv\sqrt{D\tau_{\text{s}}}$ remains constant Mal‘shukov et al. (1996). It has also been suggested that the independence of $L_{\text{s}}$ on disorder should extend to the strongly localized or insulating regime of transport Shklovskii (2006). Our simulations are a direct confirmation of the universality of $L_{\text{s}}$ in the diffusive and weakly localized regimes, but further numerical simulations with stronger disorder would be needed to confirm this behavior in the insulating regime.

In addition to the DP mechanism, it has been predicted that in the presence of Rashba SOC the out-of-plane spin can also be relaxed by the Elliott-Yafet (EY) mechanism Ochoa et al. (2012). In this mechanism each scattering event has a finite probability to flip the spin, and the spin lifetime is given by $\tau_{\text{s}}=(E_{\text{F}}/(2\lambda_{\text{R}}|\mathbf{E}|))^{2}\tau_{\text{p}}$. The dashed lines in Fig. 4(c) show the expected spin relaxation in the presence of both DP and EY mechanisms assuming $\tau_{\text{p}}=6.5$ fs, corresponding to the average $\tau_{\text{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 $\tau_{\text{s}}$ with $E_{\text{F}}$ was a consequence of the Dirac-like linear dispersion of the electrons. However, this linear dispersion relation does not necessarily describe the transport of electrons in graphene GBs, and thus the above expression 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 diffusion length in polycrystalline graphene is independent of grain size, $L_{\text{s}}=v_{\text{F}}/\Omega_{\text{R}}$. This means that for a Rashba SOC strength of $50$ $\upmu$eV ($5$ $\upmu$eV), the spin diffusion length will be $L_{\text{s}}\approx 5$ $\upmu$m ($50$ $\upmu$m) whether the graphene grains are 10 nm or 10 $\upmu$m in diameter. Thus, for spintronics applications, single-domain graphene may not be a necessity and focus can be placed on eliminating other sources of spin relaxation, such as magnetic impurities Kochan et al. (2014). 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 relax spin. Theoretical studies of single-crystal graphene concluded that the impact of corrugations is weak, with $L_{\text{s}}>100$ $\upmu$m Vicent et al. (2017). Meanwhile, a variety of AFM studies were unable to identify any height variation arising from graphene GBs Červenka et al. (2009); Yu et al. (2011); Fei et al. (2013); Ogawa et al. (2014), and if present it is limited to <$$0.3 Å Červenka and Flipse (2009). It is unclear whether localized corrugations of this height would contribute to spin relaxation in CVD graphene, and this would therefore be an intriguing future direction of research. Measurements of CVD graphene suggest that spin relaxation at GBs is not very strong Kamalakar et al. (2015); Gebeyehu et al. (2019), but in the worst-case scenario we would expect corrugations at the GBs to limit $L_{\text{s}}$ to the graphene grain size, $L_{\text{s}}\sim 10$ $\upmu$m. Beyond corrugations, historically it is believed that spin transport in graphene has been limited by paramagnetic impurities and contact-induced dephasing. In early measurements these effects limited spin lifetimes to hundreds of picoseconds and $L_{\text{s}}\sim 1$ $\upmu$m Tombros et al. (2007). Improvements in these areas now permit the measurement of $L_{\text{s}}>10$ $\upmu$m, and as high as $30$ $\upmu$m in the best nonlocal measurements to date Drögeler et al. (2016). Altogether, the spin diffusion length arising from DP spin relaxation is therefore on par with these other sources of spin relaxation, and may even be the dominant mechanism in the highest-quality CVD samples currently being produced.