Nonlocal Spin Dynamics in the Crossover from Diffusive to Ballistic Transport

Improved fabrication techniques have enabled the possibility of ballistic transport and unprecedented spin manipulation in ultraclean graphene devices. Spin transport in graphene is typically probed in a nonlocal spin valve and is analyzed using spin diffusion theory, but this theory is not necessarily applicable when charge transport becomes ballistic or when the spin diffusion length is exceptionally long. Here, we study these regimes by performing quantum simulations of graphene nonlocal spin valves. We find that conventional spin diffusion theory fails to capture the crossover to the ballistic regime as well as the limit of long spin diffusion length. We show that the latter can be described by an extension of the current theoretical framework. Finally, by covering the whole range of spin dynamics, our study opens a new perspective to predict and scrutinize spin transport in graphene and other two-dimensional material-based ultraclean devices.

Since the seminal work of Tombros and coworkers [1], who first measured long spin diffusion length in graphene nonlocal spin devices, a considerable number of studies have explored how to improve the material quality and the efficiency of spin injection and detection so as to reach the upper limit of spin transport [2][3][4][5][6][7][8][9][10][11]. After fifteen years of progress, the fabrication of ultraclean (ballistic) graphene devices is now a reality with mean free paths reaching hundreds of nanometers and as long as tens of µm at lower temperatures [12,13]. Theoretical analysis of experimental data is usually based on the spin diffusion equations [14][15][16][17][18][19][20], but their validity in new regimes of spin transport, especially the ballistic regime, is under question [21,22]. A solution is to use quantum transport simulations to describe the spin dynamics in a realistic device geometry [23,24]. However, theoretical effort in this direction is currently lacking, which not only limits the understanding of spin transport in ultraclean devices but also restrains further improvement for spintronic applications based on two-dimensional materials and van der Waals heterostructures [25].
In this Letter we use quantum transport simulations to explore the physics of spin dynamics in graphene nonlocal spin valves (NSVs). By calculating exactly the nonlocal resistance R nl , we are able to capture regimes that conventional analysis fails to describe. In the diffusive regime we show that the typically overlooked drain and reference electrodes (see Fig. 1 below) play a fundamental role in limiting R nl when spin relaxation is weak. When approaching the quasiballistic regime of spin transport, simulated Hanle precession curves reveal the failure of the diffusive equations. By extending the theory of spin diffusion following Refs. [26,27], we obtain more general formulas for properly obtaining the spin diffusion length λ s in the former case, and for understanding the evolution of Hanle curves in the crossover from diffusive to ballistic transport. Our findings demonstrate that brute force quantum simulation of nonlocal transport is fundamental to properly analyze spin dynamics in unconventional regimes and for allowing a direct comparison with experimental data.
Lateral NSVs are widely used to probe spin transport in disordered materials because they decouple electrical currents from spin currents, allowing for better device sensitivity [14,26,[28][29][30][31]. In such a device configuration, shown in Fig. 1, a ferromagnetic (FM) electrode (labeled "2") drives a spin-polarized current I 0 to a drain electrode ("1"), a spin accumulation develops below the FM, and this spin diffuses to the right along the channel. This spin is detected by another FM contact ("3") as a nonlocal voltage V nl , which is normalized to a nonlocal resistance R nl ≡ V nl /I 0 . Owing to inherent spin relaxation, R nl usually decays exponentionally with channel length d, R nl ∝ exp(−d/λ s ). Extracting λ s from the length dependence of R nl is quite difficult experimentally, but this can be avoided by applying a perpendicular magnetic field which provokes precession and additional dephasing of the spins. According to traditional spin diffusion theory, R nl then takes the form [14,26,[30][31][32][33] is the polarization of the injector (detector) FM contact, σ is the electrical conductivity, w is the channel width, D is the diffusion coefficient, and ω = gµ B B/ is the Larmor spin precession frequency induced by the magnetic field B with g the g-factor, µ B the Bohr magneton, and the Planck constant. By fitting this expression to a measurement of R nl vs. B, one can extract the spin diffusion length λ s . This approach, known as a Hanle measurement, has been a cornerstone of the exploration of spin dynamics in a large variety of materials including metals [26,34], semiconductors [14,[35][36][37][38], and graphene [7,21]. However, Eq. (1) is based on several assumptions which may be violated in ultraclean devices. Specifically, Eq. (1) assumes that transport is fully diffusive, that relaxation is fast enough so that no spin signal reaches the reference electrode (lead 4), and neglects diffusion along the left-hand part of the NSV (between leads 2 and 1). It is therefore important to revisit and extend the current theoretical framework to cope with new spin transport regimes which are emerging in today's ultraclean nonlocal spin devices.
To study the behavior of graphene NSVs, we employ the Landauer-Büttiker formalism, as implemented in Kwant [39], to the device setup in Fig. 1. The graphene layer is described in a single-π-orbital tight-binding basis, with a Hamiltonian given bŷ where c † iη (c iη ) is the creation (annihilation) operator with spin η on site i. The first term (t = −2.6 eV) denotes nearest-neighbor hopping in the graphene honeycomb lattice. The second term is Anderson disorder defined by a random potential at each site The third term is magnetic disorder mainly affecting the spin dynamics. It is defined as a magnetic exchange coupling with strength J and random orientation at each site i, , with θ and φ spherical angles and s the spin Pauli matrices. The last term is the Zeeman exchange induced by an external magnetic field B (note that orbital effects of the magnetic field are neglected). In general, U is taken to be much larger than J, such that U dictates the charge transport regime, whereas the spin relaxation is driven by J. The modeling of the leads is described in Supplemental Material [40].
The calculation of R nl is performed by evaluating all transmission probabilities between different leads. We then construct the conductance matrix G [41] and solve the linear system I = GV , where I and V are vectors including the current and voltage conditions at each electrode. We fix a current I 0 from lead 2 (injector) to lead 1 (drain) while enforcing that no current flows in leads 3 (detector) and 4 (reference). This ensures zero charge current in the channel since any current going to the right from the injector will be compensated with an oppositely spin-polarized current injected by the reference lead. We also ground the drain (V 1 = 0) and solve the system to obtain the other voltages. The nonlocal resistance is then calculated as We first investigate spin dynamics in the diffusive regime of charge transport, which is identified from the scaling of the two-terminal conductance G 2T with the channel length x = L + l by removing leads 2 and 3. We evaluate the mean free path l e by fitting the numerical result to G 2T = 2e 2 /h × M l e /x, with M the number of propagating modes per spin. Also, for x l e we calculate the localization length l loc using ln(G 2T ) ∝ −x/l loc . By choosing w = 20.1 nm (164-aGNR), Fermi energy E F = 0.4 eV, M = 9 and U = 1.04 eV, we obtain l e = 117 nm and l loc = 880 nm [40]. We take L = 250 nm and l = 1000 nm so that most of the transport occurs between l e and l loc , and we compute R nl vs. the channel length d at B = 0 for different magnetic disorder strengths J. The results are plotted in Fig. 2.
For large values of J, R nl decays exponentially with channel length, as predicted by Eq. (1). However, as spin relaxation slows with decreasing J, the decay of R nl becomes linear instead of exponential. Even for J = 0, corresponding to λ s → ∞, there is a loss of spin signal with channel length [40], a result not captured by Eq. (1). Conventional spin diffusion theory assumes the spin accumulation vanishes at x → +∞, or at least at x = l [14,31]. However, this condition is violated for the lowest values of J in our simulations, and may also be the case in recent experiments for which λ s reaches tens of µm [42].
To describe the proper length dependence of R nl , we solve the spin diffusion equations taking the full device geometry into account [40]; not only are spins injected from lead 2, but leads 1 and 4 are explicitly included (lead 3 does not perturb the system). From this, R nl becomes where β = R c wσα and R c is the contact resistance between leads 1 and 4 and the graphene. In the case of perfectly transparent contacts, the interface resistance is not zero but dictated by the Sharvin resistance R S = h/(2e 2 M ) [43,44]. If one takes the limits λ s L, l, Eq. (1) is recovered. Importantly, Eq. (3) becomes linear when λ s → ∞, where R L = L/wσ and R l = l/wσ are the sheet resistance of the left and right device regions, respectively. The black dot-dashed lines in Fig. 2 show the fits of the numerical results to Eq. (3), indicating that this expression is able to capture the scaling of R nl for any value of J. Equation (4) shows that when λ s ≥ L, l the nonlocal spin signal still decays with length. This decay is no longer related solely to spin relaxation but also to charge diffusion and the presence of the leads. Recall that R nl depends on the conductance matrix G, which consists of the transmission between all leads and the imposed current/voltage conditions. In the limit of long λ s , the drain and reference electrodes act as spin sinks, fixing the value of R nl in order to meet the conditions I = I 0 and I = 0 at leads 1 and 4, respectively. We note that this spin sinking effect occurs despite the absence of spin relaxation in the leads. Rather it is the result of these leads absorbing and reinjecting spin current under the imposed boundary conditions. This contrasts with the contact-induced spin dephasing discussed in some experiments [15,19,45]. Equation (4) shows that at the reference electrode (d = l) R nl is proportional to R c to leading order. Thus, in the limit of weak spin relaxation a small R c will suppress the nonlocal spin signal. Another consequence of long λ s is that the transmission between the drain and reference electrodes becomes crucial. The condition of zero charge current in the channel forces the injection of spin-down current from lead 4 to 1 so that lead 2 can inject up-spins that diffuse towards lead 3. If lead 4 (1) is unable to inject (absorb) down-spins to (from) the system, up-spins will not be able to diffuse along the channel and R nl will be suppressed. In Fig. S3 [40], such effect is evidenced further by changing lead 1 from nonmagnetic to FM, which reduces R nl by more than three orders of magnitude. This suggests not employing FM materials for leads 1 and 4 in experiments. Another important consequence is that Eq. (1) will underestimate the spin diffusion length because it neglects an extra source of spin signal decay. This is shown in Fig. 2 (inset), where λ s is plotted vs. spin relaxation strength. The gray squares are extracted from fits to Eq. (1), while the black circles are from Eq. (3). The spin diffusion length is the same when λ s < L, l (large J), but for small J Eq. (1) significantly underestimates the value of λ s . According to the theory of spin relaxation arising from exchange fluctuations, λ s should scale as 1/J [14,40], which is captured by the fits to Eq. (3).
We now extend the analysis to Hanle precession and plot R nl vs. B in Fig. 3 using a channel length of d = 500 nm. We note that large magnetic fields are required for computational convenience (our device is smaller than in experiments) but no spurious effect is introduced since the Zeeman splitting remains much smaller than the subband energy separation and orbital effects are excluded. The simulation data are fitted with Eqs. (1) and (3) using D = v F l e and the resulting λ s are compared in Fig. 3 (inset). Similar to Fig. 2, in the limit of weak spin relaxation λ s is underestimated when using the con-ventional equation.
To estimate how strong the underestimation of λ s may be in state-of-the-art devices, we calculate a Hanle curve with Eq. (3) using realistic parameters (L = 5 µm, l = 20 µm, d = 15 µm, λ s = 10 µm, D = 0.05 m 2 /s) and fit it with Eq. (1). We obtain λ s = 7.26 µm, about 25% less than the real value. More importantly, the spin lifetime (τ s = λ 2 s /D) is underestimated by nearly 100%; the real value is 2 ns while the fit gives 1.05 ns. These results thus call for a revised analysis of Hanle spin precession measurements taking into account the device geometry and our more general formula (Eq. (3)).
Finally, we examine the quasiballistic limit. When l e ∼ d, only a few scattering events occur during transport through the channel. This situation has been discussed for spin relaxation in ultraclean graphene [46], but little is known about its impact on Hanle measurements. We keep all simulation parameters the same as before and reduce the Anderson disorder to U = 0.52 eV, giving a mean free path l e = 500 − 600 nm. The solid lines in Fig. 4 show the Hanle curves for this quasiballistic regime, and the dashed lines show fits using Eq. (3). The behavior of R nl is now substantially different, especially with respect to the dependence on B. The unit B 0 corresponds to the magnetic field needed for spin to precess 2π radians upon reaching the detector (B 0 = 2πvF γd , with γ = gµ B / the gyromagnetic ratio and v F the Fermi velocity). The first rotation of the spins occurs at B 0 = 1, but is followed by a dispersion of frequencies for larger B. This can be understood if we examine the origin of such oscillations. By performing a simulation in a purely ballistic regime, U = J = 0, we observe in Fig. 4 (inset) that the main oscillation has the same period, but is superimposed with other frequencies. This arises because the nonlocal signal is the sum of each propagating mode moving at a different velocity. To verify this, we follow Ref. 26 and sum the contributions of spin over all transport times to get R nl ∝ M i cos (γdB/v F,i ). By taking only the Fermi velocities of each mode of the system, the simulations are very well reproduced.
From these results we can conclude that in the quasiballistic regime the scattering is weak enough for the precession to follow that of ballistic transport, but is also strong enough to average the beating pattern to one main frequency. This explains why neither Eq. (3) nor the sum of cosines is able to fit the quasiballistic Hanle curves. This difficulty in capturing the crossover from diffusive to ballistic transport in a single expression highlights the importance of brute force quantum simulations to understand such a regime. Finally, in the limit of a 2D graphene flake, with most electrons moving at the same Fermi velocity, one would expect the signal to be determined by a single frequency. To observe this effect at low magnetic fields (B ≤ 0.5 T), the channel length needs to be d ≥ 2πvF γB ≈ 50 µm. In conclusion, we have developed a simulation tool that provides a more global picture of nonlocal spin transport when the material quality drives the system towards the quasiballistic regime, as well as an extended theoretical frame to analyze systems with long spin diffusion lengths. In this limit, the drain and reference electrodes become the limiting factors, and one should aim for these to be nonmagnetic and optimize their contact resistance to reach the upper limit for spin information transfer. Beyond guiding future nonlocal spin transport measurements in graphene devices, the developed methods and findings should be also relevant for other types of two-dimensional materials and van der Waals heterostructures.
We thank Ivan Vera-Marun for highly fruitful comments. M.V. acknowledges support from "La Caixa" Foundation. S.R.P. acknowledges funding from the Irish Research Council under the Laureate awards programme. X.W. acknowledges the ANR GRANSPORT funding. All authors were supported by the European Union Horizon 2020 research and innovation programme under Grant Agreement No. 785219 (Graphene Flagship). ICN2 is funded by the CERCA Programme/Generalitat de Catalunya, and is supported by the Severo Ochoa program from Spanish MINECO (Grant No. SEV-2017-0706).