This is the accepted version of the article: Vila M., Garcia J.H., Cummings A.W., Power S.R., Groth C.W., Waintal X., Roche S.. Nonlocal Spin Dynamics in the Crossover from Diffusive to Ballistic Transport. Physical Review Letters, (2020). 124. 196602: - . 10.1103/PhysRevLett.124.196602. Available at: https://dx.doi.org/10.1103/PhysRevLett.124.196602 Nonlocal Spin Dynamics in the Crossover from Diffusive to Ballistic Transport Marc Vila,1, 2 Jose H. Garcia,1 Aron W. Cummings,1 Stephen R. Power,1, 3, 4 Christoph Groth,5 Xavier Waintal,5 and Stephan Roche1, 6 1 Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and BIST, Campus UAB, Bellaterra, 08193 Barcelona, Spain 2 Department of Physics, Universitat Autònoma de Barcelona, Campus UAB, Bellaterra, 08193 Barcelona, Spain 3 Universitat Autònoma de Barcelona, Campus UAB, Bellaterra, 08193 Barcelona, Spain 4 School of Physics, Trinity College Dublin, Dublin 2, Ireland 5 Univ. Grenoble Alpes, CEA, IRIG-PHELIQS, 38000 Grenoble, France 6 ICREA–Institució Catalana de Recerca i Estudis Avançats, 08010 Barcelona, Spain Improved fabrication techniques have enabled the possibility of ballistic transport and unprece- dented 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 captured 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 cowork- the typically overlooked drain and reference elec- ers [1], who first measured long spin diffusion length trodes (see Fig. 1 below) play a fundamental role in graphene nonlocal spin devices, a considerable in limiting Rnl when spin relaxation is weak. When number of studies have explored how to improve approaching the quasiballistic regime of spin trans- the material quality and the efficiency of spin injec- port, simulated Hanle precession curves reveal the tion and detection so as to reach the upper limit of failure of the diffusive equations. By extending the spin transport [2–11]. After fifteen years of progress, theory of spin diffusion following Refs. [26, 27], we the fabrication of ultraclean (ballistic) graphene de- obtain more general formulas for properly obtain- vices is now a reality with mean free paths reaching ing the spin diffusion length λs in the former case, hundreds of nanometers and as long as tens of µm and for understanding the evolution of Hanle curves at lower temperatures [12, 13]. Theoretical anal- in the crossover from diffusive to ballistic transport. ysis of experimental data is usually based on the Our findings demonstrate that brute force quantum spin diffusion equations [14–20], but their validity in simulation of nonlocal transport is fundamental to new regimes of spin transport, especially the ballistic properly analyze spin dynamics in unconventional regime, is under question [21, 22]. A solution is to regimes and for allowing a direct comparison with use quantum transport simulations to describe the experimental data. spin dynamics in a realistic device geometry [23, 24]. However, theoretical effort in this direction is cur- Lateral NSVs are widely used to probe spin trans- rently lacking, which not only limits the understand- port in disordered materials because they decouple ing of spin transport in ultraclean devices but also electrical currents from spin currents, allowing for restrains further improvement for spintronic appli- better device sensitivity [14, 26, 28–31]. In such a de- cations based on two-dimensional materials and van vice configuration, shown in Fig. 1, a ferromagnetic der Waals heterostructures [25]. (FM) electrode (labeled “2”) drives a spin-polarized current I0 to a drain electrode (“1”), a spin accu- In this Letter we use quantum transport simu- mulation develops below the FM, and this spin dif- lations to explore the physics of spin dynamics in fuses to the right along the channel. This spin is graphene nonlocal spin valves (NSVs). By calculat- detected by another FM contact (“3”) as a nonlo- ing exactly the nonlocal resistance Rnl , we are able cal voltage Vnl , which is normalized to a nonlocal to capture regimes that conventional analysis fails resistance Rnl ≡ Vnl /I0 . Owing to inherent spin to describe. In the diffusive regime we show that relaxation, Rnl usually decays exponentionally with channel length d, Rnl ∝ exp(−d/λs ). Extracting λs from the length dependence of Rnl is quite difficult experimentally, but this can be avoided by applying a perpendicular magnetic field which provokes pre- cession and additional dephasing of the spins. Ac- cording to traditional spin diffusion theory, Rnl then takes the form [14, 26, 30–33] e−dα   Pi Pd Rnl = Re , (1) 2wσ α where α = , Pi (Pd ) is the polarization of q 1 ω λ2s + iD the injector (detector) FM contact, σ is the electri- cal conductivity, w is the channel width, D is the diffusion coefficient, and ω = gµB B/~ is the Larmor FIG. 1. Sketch of the lateral nonlocal spin valve. Red (black) regions denote the contacts (sample). The injec- spin precession frequency induced by the magnetic tor and detector contacts, labeled 2 and 3 respectively, field B with g the g-factor, µB the Bohr magneton, are ferromagnetic with their magnetization indicated by and ~ the Planck constant. By fitting this expres- arrows. Contacts 1 and 4 represent the drain and refer- sion to a measurement of Rnl vs. B, one can extract ence electrodes, respectively. the spin diffusion length λs . This approach, known as a Hanle measurement, has been a cornerstone of the exploration of spin i, with δUi ∈ [−U/2, U/2]. The third term is mag- dynamics in a large variety of materials including netic disorder mainly affecting the spin dynamics. metals [26, 34], semiconductors [14, 35–38], and It is defined as a magnetic exchange coupling with graphene [7, 21]. However, Eq. (1) is based on sev- strength J and random orientation at each site i, eral assumptions which may be violated in ultraclean Ji = J [sin(θi ) cos(φi ), sin(θi ) sin(φi ), cos(θi )], with devices. Specifically, Eq. (1) assumes that transport θ and φ spherical angles and s the spin Pauli matri- is fully diffusive, that relaxation is fast enough so ces. The last term is the Zeeman exchange induced that no spin signal reaches the reference electrode by an external magnetic field B (note that orbital (lead 4), and neglects diffusion along the left-hand effects of the magnetic field are neglected). In gen- part of the NSV (between leads 2 and 1). It is eral, U is taken to be much larger than J, such that therefore important to revisit and extend the cur- U dictates the charge transport regime, whereas the rent theoretical framework to cope with new spin spin relaxation is driven by J. The modeling of the transport regimes which are emerging in today’s ul- leads is described in Supplemental Material [40]. traclean nonlocal spin devices. The calculation of Rnl is performed by evaluat- To study the behavior of graphene NSVs, we em- ing all transmission probabilities between different ploy the Landauer-Büttiker formalism, as imple- leads. We then construct the conductance matrix G mented in Kwant [39], to the device setup in Fig. 1. [41] and solve the linear system I = GV , where I The graphene layer is described in a single-π-orbital and V are vectors including the current and voltage tight-binding basis, with a Hamiltonian given by conditions at each electrode. We fix a current I0 X † from lead 2 (injector) to lead 1 (drain) while enforc- δUi c†iη ciη + X X † Ĥ = t ciη cjη + ciη [s · Ji ]ηη0 ciη0 ing that no current flows in leads 3 (detector) and 4 hi,ji i iηη 0 (reference). This ensures zero charge current in the c†iη [s · Bi ]ηη0 ciη0 . (2) channel since any current going to the right from X + µB iηη 0 the injector will be compensated with an oppositely spin-polarized current injected by the reference lead. where c†iη (ciη ) is the creation (annihilation) operator We also ground the drain (V1 = 0) and solve the with spin η on site i. The first term (t = −2.6 eV) system to obtain the other voltages. The nonlocal denotes nearest-neighbor hopping in the graphene resistance is then calculated as Rnl = (V3 − V4 )/I0 . honeycomb lattice. The second term is Anderson We first investigate spin dynamics in the diffu- disorder defined by a random potential at each site sive regime of charge transport, which is identified 2 from the scaling of the two-terminal conductance decay of Rnl becomes linear instead of exponential. G2T with the channel length x = L + l by removing Even for J = 0, corresponding to λs → ∞, there leads 2 and 3. We evaluate the mean free path le by is a loss of spin signal with channel length, a result fitting the numerical result to G2T = 2e2 /h×M le /x, not captured by Eq. (1). Conventional spin diffu- with M the number of propagating modes. Also, for sion theory assumes the spin accumulation vanishes x  le we calculate the localization length lloc us- at x → +∞, or at least at x = l [14, 31]. However, ing hln(G2T )i ∝ −x/lloc . By choosing U = 1.04 eV, this condition is violated for the lowest values of J in we obtain le = 117 nm and lloc = 880 nm [40]. We our simulations, and may also be the case in recent take L = 250 nm and l = 1000 nm so that most experiments for which λs reaches tens of µm [42]. of the transport occurs between le and lloc , and we compute Rnl vs. the channel length d at B = 0 for To describe the proper length dependence of Rnl , different magnetic disorder strengths J. The results we solve the spin diffusion equations taking the full are plotted in Fig. 2. device geometry into account [40]; not only are spins For large values of J, Rnl decays exponentially injected from lead 2, but leads 1 and 4 are explicitly with channel length, as predicted by Eq. (1). How- included (lead 3 does not perturb the system). From ever, as spin relaxation slows with decreasing J, the this, Rnl becomes   Pi Pd [β cosh(Lα) + 4 sinh(Lα)] · [β cosh(α(d − l)) − 4 sinh(α(d − l))] Rnl = Re , (3) 2wσ α[4β cosh((L + l)α) + (8 + β 2 /2) sinh((L + l)α)] where β = Rc wσα and Rc 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 RS = h/(2e2 M ) [43, 44]. If one takes the limits λs  L, l, Eq. (1) is recovered. Importantly, Eq. (3) becomes linear when λs → ∞, Pi Pd (4RL + Rc )(−4d + 4l + Rc wσ) Rnl = , (4) 2wσ 8RL + 8Rl + 4Rc where RL = L/wσ and Rl = l/wσ are the sheet re- sistance of the left and right device regions, respec- tively. The black dot-dashed lines in Fig. 2 show the FIG. 2. Rnl as a function of injector-detector distance fits of the numerical results to Eq. (3), indicating for different strengths of magnetic disorder. Error bars that this expression is able to capture the scaling of result from the averaging of several disorder configura- Rnl for any value of J. tions (> 130). Black dot-dashed lines are the fits using Equation (4) shows that when λs ≥ L, l the non- Eq. (3). Inset: comparison of λs extracted from Eq. (1) local spin signal still decays with length. This decay (gray squares) and Eq. (3) (black circles). The red line indicates 1/J scaling of λs . is no longer related solely to spin relaxation but also to charge diffusion and the presence of the leads. Recall that Rnl depends on the conductance matrix G, which consists of the transmission between all experiments [19, 45]. Equation (4) shows that at leads and the imposed current/voltage conditions. the reference electrode (d = l) Rnl is proportional In the limit of long λs , the drain and reference elec- to Rc to leading order. Thus, in the limit of weak trodes act as spin sinks, fixing the value of Rnl in spin relaxation a small Rc will suppress the nonlo- order to meet the conditions I = I0 and I = 0 at cal spin signal. Another consequence of long λs is leads 1 and 4, respectively. We note the absence that the transmission between the drain and refer- of spin relaxation in these leads, which results from ence electrodes becomes crucial. The condition of spin current absorption and injection under the im- zero charge current in the channel forces the injec- posed boundary conditions. This contrasts with the tion of spin-down current from lead 4 to 1 so that contact-induced spin dephasing discussed in some lead 2 can inject up-spins that diffuse towards lead 3. 3 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 m2 /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 = λ2s /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 mea- surements taking into account the device geometry FIG. 3. Hanle spin precession curves for different and our more general formula (Eq. (3)). strengths of magnetic disorder. Error bars result from averaging several disorder configurations (> 90). Black Finally, we examine the quasiballistic limit. When dot-dashed line is the fit using Eq. (3). Inset: comparison le ∼ d, only a few scattering events occur during of λs extracted from Eq. (1) (gray squares) and Eq. (3) transport through the channel. This situation has (black circles). been discussed for ultraclean graphene devices [46], but little is known about its impact on Hanle mea- surements. We keep all simulation parameters the If lead 1 (4) is unable to inject (absorb) down-spins same as before and reduce the Anderson disorder to to (from) the system, up-spins will not be able to dif- U = 0.52 eV, giving a mean free path le = 500 − 600 fuse along the channel and Rnl will be suppressed. nm. The solid lines in Fig. 4 show the Hanle curves In Fig. S3 [40], such effect is evidenced further by for this quasiballistic regime, and the dashed lines changing lead 1 from nonmagnetic to FM, which re- show fits using Eq. (3). The behavior of Rnl is now duces Rnl by more than three orders of magnitude. substantially different, especially with respect to the This suggests not employing FM materials for leads dependence on B. The unit B0 corresponds to the 1 and 4 in experiments. Another important conse- magnetic field needed for spin to precess 2π radi- quence is that Eq. (1) will underestimate the spin ans upon reaching the detector (B0 = 2πv γd , with F diffusion length because it neglects an extra source γ = gµB /~ the gyromagnetic ratio and vF the Fermi of spin signal decay. This is shown in Fig. 2 (inset), velocity). The first rotation of the spins occurs at where λs is plotted vs. spin relaxation strength. The B0 = 1, but is followed by a dispersion of frequencies gray squares are extracted from fits to Eq. (1), while for larger B. This can be understood if we examine the black circles are from Eq. (3). The spin diffu- the origin of such oscillations. By performing a sim- sion length is the same when λs < L, l (large J), but ulation in a purely ballistic regime, U = J = 0, for small J Eq. (1) significantly underestimates the we observe in Fig. 4 (inset) that the main oscillation value of λs . According to the theory of spin relax- has the same period, but is superimposed with other ation arising from exchange fluctuations, λs should frequencies. This arises because the nonlocal signal scale as 1/J [14, 40], which is captured by the fits to is the sum of each propagating mode moving at a Eq. (3). different velocity. To verify this, we follow Ref. 26 We now extend the analysis to Hanle precession and sum the contributions PM of spin over all transport and plot Rnl vs. B in Fig. 3 using a channel length times to get Rnl ∝ i cos (γdB/vF,i ). By taking of d = 500 nm. We note that large magnetic fields only the Fermi velocities of each mode of the system, are required for computational convenience (our de- the simulations are very well reproduced. vice is smaller than in experiments) but no spuri- From these results we can conclude that in the ous effect is introduced since the Zeeman splitting quasiballistic regime the scattering is weak enough remains much smaller than the subband energy sep- for the precession to follow that of ballistic trans- aration and orbital effects are excluded. The sim- port, but is also strong enough to average the beat- ulation data are fitted with Eqs. (1) and (3) using ing pattern to one main frequency. This explains D = vF le and the resulting λs are compared in Fig. 3 why neither Eq. (3) nor the sum of cosines is able (inset). Similar to Fig. 2, in the limit of weak spin to fit the quasiballistic Hanle curves. This difficulty 4 and is supported by the Severo Ochoa program from Spanish MINECO (Grant No. SEV-2017-0706). [1] N. Tombros, C. Jozsa, M. Popinciuc, H. T. Jonkman, and B. J. van Wees, Nature 448, 571 (2007). [2] N. Tombros, S. Tanabe, A. Veligura, C. Jozsa, M. Popinciuc, H. T. Jonkman, and B. J. van Wees, Phys. Rev. Lett. 101, 046601 (2008). [3] W. Han and R. K. Kawakami, Phys. Rev. Lett. 107, 047207 (2011). [4] B. Dlubak, M.-B. Martin, C. Deranlot, B. Servet, FIG. 4. Hanle spin precession curves in the quasiballis- S. Xavier, R. Mattana, M. Sprinkle, C. Berger, tic regime. Solid lines correspond to simulations, while W. A. De Heer, F. Petroff, A. Anane, P. Seneor, dashed lines are fits using Eq. (3). Inset: Case with and A. Fert, Nat. Phys. 8, 557 (2012). U = J = 0, solid (dashed) line shows the simulation [5] P. J. Zomer, M. H. D. Guimarães, N. Tombros, and (Rnl ∝ M B. J. van Wees, Phys. Rev. B 86, 161416 (2012). P i cos (γdB/vF,i )). [6] M. H. D. Guimarães, P. J. Zomer, J. Ingla-Aynés, J. C. Brant, N. Tombros, and B. J. van Wees, Phys. in capturing the crossover from diffusive to ballis- Rev. Lett. 113, 086602 (2014). tic transport in a single expression highlights the [7] W. Han, R. K. Kawakami, M. Gmitra, and J. Fabian, Nat. Nanotechnol. 9, 794 (2014). importance of brute force quantum simulations to [8] M. Gurram, S. Omar, S. Zihlmann, P. Makk, understand such a regime. Finally, in the limit of C. Schönenberger, and B. J. van Wees, Phys. Rev. a 2D graphene flake, with most electrons moving at B 93, 115441 (2016). the same Fermi velocity, one would expect the signal [9] M. V. Kamalakar, C. Groenveld, A. Dankert, and to be determined by a single frequency. To observe S. P. Dash, Nat. Commun. 6, 6766 (2015). this effect at low magnetic fields (B ≤ 0.5 T), the [10] I. G. Serrano, J. Panda, F. Denoel, O. Vallin, channel length needs to be d ≥ 2πv F D. Phuyal, O. Karis, and M. V. Kamalakar, Nano γB ≈ 50 µm. Lett. 19, 666 (2019). In conclusion, we have developed a simulation tool [11] Z. M. Gebeyehu, S. Parui, J. F. Sierra, M. Timmer- that provides a more global picture of nonlocal spin mans, M. J. Esplandiu, S. Brems, C. Huyghebaert, transport when the material quality drives the sys- K. Garello, M. V. Costache, and S. O. Valenzuela, tem towards the quasiballistic regime, as well as 2D Mater. 6, 034003 (2019). an extended theoretical frame to analyze systems [12] L. Wang, I. Meric, P. Y. Huang, Q. Gao, Y. Gao, with long spin diffusion lengths. In this limit, the H. Tran, T. Taniguchi, K. Watanabe, L. M. Cam- drain and reference electrodes become the limiting pos, D. A. Muller, J. Guo, P. Kim, J. Hone, K. L. Shepard, and C. R. Dean, Science 342, 614 (2013). factors, and one should aim for these to be non- [13] L. Banszerus, M. Schmitz, S. Engels, M. Goldsche, magnetic and optimize their contact resistance to K. Watanabe, T. Taniguchi, B. Beschoten, and reach the upper limit for spin information trans- C. Stampfer, Nano Lett. 16, 1387 (2015). fer. Beyond guiding future nonlocal spin transport [14] J. Fabian, A. Matos-Abiaguea, P. Ertler, measurements in graphene devices, the developed C. snd Stano, and I. Zutic, Acta Phys. Slo- methods and findings should be also relevant for vaca 57, 565 (2007). other types of two-dimensional materials and van [15] T. Maassen, I. J. Vera-Marun, M. H. D. Guimarães, and B. J. van Wees, Phys. Rev. B 86, 235408 (2012). der Waals heterostructures. [16] M. Wojtaszek, I. J. Vera-Marun, and B. J. van We thank Ivan Vera-Marun for highly fruitful Wees, Phys. Rev. B 89, 245427 (2014). comments. M.V. acknowledges support from “La [17] E. Sosenko, H. Wei, and V. Aji, Phys. Rev. B 89, Caixa” Foundation. X.W. acknowledges the ANR 245436 (2014). GRANSPORT funding. All authors were supported [18] H. Idzuchi, Y. Fukuma, S. Takahashi, S. Maekawa, by the European Union Horizon 2020 research and and Y. Otani, Phys. Rev. B 89, 081308 (2014). innovation programme under Grant Agreement No. [19] H. Idzuchi, A. Fert, and Y. Otani, Phys. Rev. B 91, 241407 (2015). 785219 (Graphene Flagship). ICN2 is funded by [20] L. O’Brien, D. Spivak, N. Krueger, T. A. Peter- the CERCA Programme/Generalitat de Catalunya, 5 son, M. J. Erickson, B. Bolon, C. C. Geppert, J. Zhang, K. S. Madhukar, K. S. Flefxner, C. Palm- C. Leighton, and P. A. Crowell, Phys. Rev. B 94, strom, and P. Crowell, Nat. Phys. 3, 197–202 094431 (2016). (2007). [21] S. Roche, J. Åkerman, B. Beschoten, J.-C. Char- [36] G. Salis, A. Fuhrer, R. R. Schlittler, L. Gross, and lier, M. Chshiev, S. P. Dash, B. Dlubak, J. Fabian, S. F. Alvarado, Phys. Rev. B 81, 205323 (2010). A. Fert, M. Guimarães, et al., 2D Mater. 2, 030202 [37] O. vant Erve, A. Friedman, C. Li, J. Robinson, (2015). J. Connell, L. Lauhon, and B. Jonker, Nat. Com- [22] M. Gurram, S. Omar, and B. J. van Wees, 2D mun. 6, 7541 (2015). Mater. 5, 032004 (2018). [38] K. Hamaya, Y. Fujita, M. Yamada, M. Kawano, [23] H. X. Tang, F. G. Monzon, R. Lifshitz, M. C. Cross, S. Yamada, and K. Sawano, J. Phys. D 51, 393001 and M. L. Roukes, Phys. Rev. B 61, 4437 (2000). (2018). [24] A. N. M. Zainuddin, S. Hong, L. Siddiqui, S. Srini- [39] C. W. Groth, M. Wimmer, A. R. Akhmerov, and vasan, and S. Datta, Phys. Rev. B 84, 165306 X. Waintal, New J. Phys. 16, 063065 (2014). (2011). [40] See Supplemental Material for two-terminal conduc- [25] X. Lin, W. Yang, K. Wang, and W. Zhao, Nat. tance simulations, the modeling of leads, the deriva- Electron. 2, 274 (2019). tion of Eq. (3), the special case J = 0, the relation [26] F. Jedema, H. Heersche, A. Filip, J. Baselmans, and between Rnl and the reference electrode, and the B. Van Wees, Nature 416, 713 (2002). theory of spin relaxation by randomly fluctuating [27] C. Petitjean, D. Luc, and X. Waintal, Phys. Rev. magnetic moments.. Lett. 109, 117204 (2012). [41] S. Datta, Electronic transport in mesoscopic systems [28] M. Johnson and R. H. Silsbee, Phys. Rev. Lett. 55, (Cambridge University Press, 1997). 1790 (1985). [42] M. Drögeler, C. Franzen, F. Volmer, T. Pohlmann, [29] M. Johnson and R. H. Silsbee, Phys. Rev. B 37, L. Banszerus, M. Wolter, K. Watanabe, 5312 (1988). T. Taniguchi, C. Stampfer, and B. Beschoten, [30] I. Žutić, J. Fabian, and S. Das Sarma, Rev. Mod. Nano Lett. 16, 3533 (2016). Phys. 76, 323 (2004). [43] V. S. Rychkov, S. Borlenghi, H. Jaffres, A. Fert, and [31] M. Wu, J. Jiang, and M. Weng, Phys. Rep. 493, X. Waintal, Phys. Rev. Lett. 103, 066602 (2009). 61 (2010). [44] S. Borlenghi, V. Rychkov, C. Petitjean, and [32] S. Takahashi and S. Maekawa, Phys. Rev. B 67, X. Waintal, Phys. Rev. B 84, 035412 (2011). 052409 (2003). [45] W. Amamou, Z. Lin, J. van Baren, S. Turkyilmaz, [33] S. Takahashi and S. Maekawa, Sci. Technol. Adv. J. Shi, and R. K. Kawakami, APL Materials 4, Mater. 9, 014105 (2008). 032503 (2016). [34] E. Villamor, L. E. Hueso, and F. Casanonva, J. [46] A. W. Cummings and S. Roche, Phys. Rev. Lett. Appl. Phys. 17, 223911 (2015). 116, 086602 (2016). [35] X. Lou, C. Adelmann, C. Crooker, E. Garlid, 6