Quantum Transport in Graphene in Presence of Strain-Induced Pseudo-Landau Levels

We report on mesoscopic transport fingerprints in disordered graphene caused by strain-field induced pseudomagnetic Landau levels (pLLs). Efficient numerical real space calculations of the Kubo formula are performed for an ordered network of nanobubbles in graphene, creating pseudomagnetic fields up to several hundreds of Tesla, values inaccessible by real magnetic fields. Strain-induced pLLs yield enhanced scattering effects across the energy spectrum resulting in lower mean free path and enhanced localization effects. In the vicinity of the zeroth order pLL, we demonstrate an anomalous transport regime, where the mean free paths increases with disorder. We attribute this puzzling behavior to the low-energy sub-lattice polarization induced by the zeroth order pLL, which is unique to pseudomagnetic fields preserving time-reversal symmetry. These results, combined with the experimental feasibility of reversible deformation fields, open the way to tailor a metal-insulator transition driven by pseudomagnetic fields.

Inhomogeneous lattice deformations in graphene generate an effective gauge field modulating the electronic spectrum 1-3 . However, compared to a real magnetic field, the formation of a pseudomagnetic field preserves time reversal symmetry, having an opposite sign in the two inequivalent K and K' valleys [4][5][6] . This leads to different behavior than for real magnetic field especially when introducing disorder. Experimentally, scanning-tunneling measurements on graphene nanobubbles have revealed an electronic spectrum consisting of pseudo-Landau levels (pLL), including a zero-energy peak, showing that moderate spatial deformations can introduce pseudomagnetic field values reaching hundreds of Tesla [7][8][9] . Pseudomagnetic fields (PMF) have also been analyzed in deformed crystals by an atomically controlled arrangement of CO molecules on a gold surface 10 , graphene on Ir with intercalated Pb monolayer islands 11 , or have been harnessed for designing innovative electronic and photonic graphene devices [12][13][14][15][16][17] .
However, the intrinsic quantum transport fingerprints of graphene in presence of pseudomagnetic fields still needs investigation, especially in disordered systems. In particular, random strain fluctuations are believed to be the dominating disorder source in high-quality onsubstrate graphene devices 18,19 . A signature of the pseudomagnetic n = 0 pLL state has been predicted for stretched graphene ribbons in the form of a quadruplet low-energy conductance resonance split by edge-induced valley mixing 16 , but its experimental confirmation remains challenging, and the variable range of transport features in deformed graphene still require further indepth exploration.
Here we study the effect of pLLs, generated by an ordered network of graphene nanobubbles (or pseudomag-netic dots), using an efficient real space Kubo quantum transport methodology. We consider samples containing electron-hole puddles caused by substrate interactions where the presence of pLLs leads to several anomalous transport features resulting from an intertwined contribution of pseudomagnetic field and disorder effects. The formation of the zero-energy pLL reduces the mean free path by orders of magnitudes in the limit of low defect density (c, modeled by a random distribution of Gaussian impurities), but scales as e ∼ c up to c ∼ 1%, in contradiction with the usual Fermi golden rule argument 20 predicting e ∼ 1/c which we obtain by considering the unstrained structure. Additionally, our simulations show that e ∼ 1/ √ B s which evidences a superimposed scaling with the pseudomagnetic length. Finally, we demonstrate further unconventional behavior of the conductivity close to zero-energy where it becomes largely independent of energy for large enough disorder (above 1%). Although in the limit of sample length L → ∞, one expects an insulating behavior at zero temperature. These features manifest the strong influence that strain induced pseudomagnetic fields have on the quantum transport properties and also suggest possibilities for metalinsulator transition driven by strain fields.

I. SUPERLATTICE OF PSEUDOMAGNETIC DOTS
To describe the electronic properties of graphene, we use the common tight-binding (TB) model of graphenê DOS for a strain array (blue) with L = 200a, R = 40a, σ = 10a and local strength corresponding to Bs = 450 T (black dashed curve is the unstrained graphene DOS). Both calculations include a 0.05% concentration of impurities. The red curve indicate the LDOS in the center of the pseudomagnetic dot averaged over both sublattices for r < 1 nm (note that the zeroth level resides completely on the Bsublattice 21,22 ). The full DOS (blue) becomes the sum of the unstrained region (dashed, black), the inner part of the strained region (red) and the outer part (not shown). The symbols refer to Figure 3. Left inset shows the superlattice of pseudomagnetic dots with a lattice constant L. The central dot has been magnified for visibility illustrating the central region (r < R) with a constant pseudomagnetic field surrounded by a region with a pseudomagnetic field of opposite sign (R < r < R + 3σ). Right inset shows the energies of the DOS peaks as a function of the peak number, √ n, confirming their pLL nature as En = sign(n) 2e v 2 F Bs|n|.
where i is the onsite energy and the sum over i, j runs over nearest neighbor sites with γ 0 = −2.7 eV. To generate pseudomagnetic fields, we consider a dot with radius R subjected to a planar triaxial displacement, which in polar coordinates are given by 22,23 u = (u r , u θ ) = (u 0 r 2 sin(3θ), u 0 r 2 cos(θ)), where r and θ are the polar coordinates and u 0 determines the strength of the strain field. This displacement field gives rise to a constant PMF 23 given by B s . Here we ignore out-of-plane and curvature 21,24,25 components of the strain field in order to get a simple connection between the strain and the magnitude of the PMF. This is justified when in-plane strain dominates (to induce pLLs) and no sharp bends are present [26][27][28] . The PMF is inherently related to the first order expansion of the TB model, but we emphasize that the calculations reported below do not rely on such approximations, instead we use the modified atomic positions to modify the TB parameters. Changing the atomic positions according to the displacement field u alters the bond lengths d ij , and thereby leads to renormalized TB hopping parameters, where β = ∂log(γ)/∂log(a)| a=a0 ≈ 3.37 1 . Here a ij = a 2 0 + xx x 2 ij + yy y 2 ij +2 xy x ij y ij /a 0 denotes the modified bond length caused by the deformation, where a 0 = 0.142 nm is the equilibrium bond length and the strain tensor is νµ = ∂ µ u ν + ∂ ν u µ /2 with ν, µ = x, y 29 . Using these definitions the deformation field in Equation (2) gives rise to a magnetic field of B s = 8u 0 β/2ea 0 . At last, we note that Equation (3) can easily incorporate out-of-plane components of the strain as it only depends on the change in bond lengths. The connection between strain tensor and bond length deformation is only approximate 30 but sufficient for the present analysis of generic quantum transport fingerprints of strain-induced pLLs. Especially, we note that the effect of the PMF on the LDOS is qualitatively unchanged by relaxation using molecular dynamics methods 22,27,[31][32][33][34][35] and that second nearest neighbor terms only contribute with a scalar potential without generating pseudomagnetic effects 22,26 .
Outside the central dot region, we apply a smoothing to the strain tensor to assure a soft transition to the strain-free pristine regions. This is accomplished by applying a Gaussian transformation where R is the radius of the PMF region as shown in the left inset of Figure 1. The triaxial strain for r < R gives rise to a constant PMF, whereas for r > R an r-dependent PMF of opposite sign develops 21 . The field of opposite sign within the smoothing region arises because we apply the smoothing to the physical strain, and not to the derived PMF. We repeat this pseudomagnetic dot deformation in a periodic array with lattice constant L (see left inset of Figure 1). This system is an idealized representation of the array of self-formed bubble deformations giving rise to the pLL structure envisioned experimentally by J. Lu et al. 8 . We note that similar calculations have been performed for other array symmetries yielding qualitatively the same conclusions.

II. KUBO TRANSPORT METHODOLOGY
We study the quantum transport using an order-N , real space implementation of the Kubo approach for the conductivity σ xx (E, t) 20,36,37 . The scaling properties of σ xx are followed through the dynamics of electronic wavepackets using is the density of states, D x (E, t) = ∆X 2 (E, t)/t is the diffusion coefficient and the time-and energy-dependent mean square displacement of the wavepacket is (4) whereX(t) is the position operator in Heisenberg representation. The semi-classical conductivity can be calculated using σ SC xx (E) = e 2 ρ(E)max t (D x (E, t))/2. Calculations are performed on systems containing approximately 5 × 10 6 atoms. The state propagation is followed through a total simulation time of 14 ps and split into three intervals with different time steps ∆t 1 = 1 fs, ∆t 2 = 15 fs and ∆t 3 = 60 fs. The expansion of the time evolution operator uses Chebyshev polynomials corresponding to coefficients larger than 10 −12 . Finally, we approximate the traces using random phase states and the Lanczos method with 1500 iterations and a broadening of η = 5 meV 38 .
To mimick the effects of electron-hole puddles induced by the substrate 39 , the Hamiltonian incorporates long range impurities defined by onsite potential on the n'th is the maximum onsite energy and the range is ξ = 3a 0 . A low concentration of impurities (c = 0.05%, c = 0.1% and c = 0.2%) is distributed randomly throughout the sample allowing us to reach the diffusive regime 20,37,39,40 .

III. LOCAL DENSITY OF STATES AND MEAN FREE PATH
We first consider the density of states (DOS) of a pseudomagnetic dot array with L = 200a, R = 40a, σ = 10a (a = √ 3a 0 ), and a maximum strain of approximately 15%. The parameters lead to a PMF of 450 T (see Figure 1), corresponding to a pseudomagnetic length ( B = /eB s ∼ 1.22 nm). The pseudomagnetic length is smaller than the dot size which is necessary for the appearance of Landau quantization 6 . Similar conclusions are obtained for other geometries by varying L and R (not shown). The DOS of the PMF array in Figure 1 shows the formation of pLLs following the expected ∼ |n| behavior (Figure 1, right inset), where n is the pLL index. In particular, a strong zero-energy peak is formed, as for real magnetic fields. The peak features are superimposed with the linear dispersion characteristic of the unstrained graphene calculated (black, dashed).
To characterize the impact of pLL states on the quantum transport, we first consider the mean free path (see Figure 2a), extracted in the diffusive regime as e (E) = max t (D(E, t) We use v F = 8.6 × 10 5 m/s 20 even though it should be noted that the Fermi velocity v F is weakly strain dependent 15,[41][42][43][44][45][46] . In the absence of a strain field (dashed lines), the mean free path follows the Fermi golden rule (FGR) and scales with the impurity density as e (E) ∼ 1/c, with values and energy dependence dictated by intervalley scattering 37 . When superimposing the strain field, we observe a large difference between e (E) with (full lines) and without (dashed lines) the deformation field (see Figure 2a). For a given impurity concentration c, we find a systematic decrease of e (E) compared to the unstrained case. The enhanced scattering induced by strain is particularly significant at the pLL energies where the DOS is large, as revealed by the dips in e (E) in Figure 2a. At the centers of pLLs, the mean free path scaling is driven by the pseudomagnetic length B i.e. e (E) ∼ 1/ √ B s (see Figure 2b). Next, we focus on the low energy regime in the inset of Figure 2a. Here, we observe a surprising crossover in the mean free path near E = ±50 meV (indicated by vertical dashed lines in the inset of Figure 2a) where the scaling of e (E) with impurity density is reversed going from higher to lower energies. While the high-energy (|E| > 50 meV) behavior follows the FGR, e (E) ∼ 1/c (filled triangles in Figure 2c), the opposite behavior is observed for |E| < 50 meV (filled circles in Figure 2c). Consequently, we observe an anomalous regime (filled circles) at low energies for which e increases with disorder. This suggests that the pseudomagnetic field counteracts disorder effects. Increasing the disorder beyond a certain concentration disrupts this mechanism and we recover the more conventional decay of the mean free path with increasing disorder.

IV. LOCALIZATION EFFECTS
To further characterize the transport fingerprints caused by the strain-induced pLLs, we analyze the diffusion coefficients D(E, t) in Figure 3 at energies marked with symbols in Figures 1 and 2, allowing us to distinguish various transport regimes. A diffusive regime is observed at the energy of the first pLL (Figure 3c, square), with a constant asymptotic behavior of D(E, t) both with (full lines) and without a strain field (dashed lines). The additional scattering due to the strain field results in a lower diffusion coefficient leading to the shorter mean free path at the pLL energies as discussed above.
The low energy regime, however, is qualitatively different (Figure 3a-b). The DOS is dominated by the zeroth pLL, which induces a strong sublattice polarization 6,22,47 . Here the decrease of D(E, t) with time reveals the onset of localization effects. In Figure 3a in the absence of the strain field (dashed lines). Comparing the curves with and without strain, we conclude that localization effects are generated by the presence of the pseudomagnetic field. This is especially clear in Figure 3d showing the variation of D(E = 0.1 eV, t) while decreasing the pseudomagnetic length. The enhancement of the localization effects are evident through the reduced value and sharper decay of D(t) for higher pseudomagnetic field strengths.
The anomalous transport regime around zero-energy is also manifested in the diffusion coefficients D(E, t). Indeed, D(E, t) exhibits a qualitatively different behavior for |E| < 50 meV (Figure 3a) compared to |E| > 50 meV (Figure 3b) even though the strain-induced pLL causes localization in both regimes. In agreement with the discussion of the mean free path, the absolute value of the diffusion coefficient at low energy increases with impurity concentration.
Finally, we consider the semi-classical (σ SC (E)) and quantum (σ(E)) conductivities at low energies for different impurity concentrations (Figure 4a). As seen in the figure, σ SC (E) always remains larger than or equal to σ 0 = 4e 2 /πh (horizontal dashed line) over the whole energy spectrum, in agreement with the lower bound of the semi-classical value 48,49 . The larger value of σ SC (E ∼ 0) may be related to a low energy behavior, similar to the effect caused by zero-energy modes around lattice monovacancies [50][51][52] .
The value of σ 0 actually separates two different transport regimes, which are also identified by scrutinizing the long time behavior of D(E, t) [51][52][53][54][55][56] . As long as σ(E) > σ 0 , the system remains diffusive in a metallic (i.e. non insulating) regime. Correspondingly, the diffusion coefficient saturates at long times. On the other hand, the condition σ(E) < σ 0 implies the existence of localization effects, whose strength (weak or strong localization) is dictated by the impurity density and the length-scale (or time-scale) at which the quantum conductivity is evaluated.
The anomalous scaling of the mean free path with disorder is further translated to similar intriguing scaling for the quantum conductivity. In presence of strain, the zero-energy quantum conductivity (taken at the maximum calculation time, t max ) remains very close to the fundamental semi-classical limit σ 0 (see Figure 4b-inset). Similarly to the semi-classical conductivity (Figure 4a), the dynamical conductivity, σ(E, t max ), shows a complex energy-dependent profile evolution with defect concentration. According to Figure 4b, the conductivity for |E| ≥ 0.1 eV decays with defect density, i.e. σ(|E| ≥ 0.1 eV, t max ) ∼ 1/c. However, at lower energies such behavior is interrupted and σ(|E| ≤ 0.1 eV, t max ) instead increases with defect density. A special situation is observed in Figure 4b-inset where σ(|E| ∼ 0 eV, t max ) is almost constant for a wide range of disorder.

V. ENERGY-DEPENDENT RANDOM DISORDER EFFECT
In this section, we provide a tentative interpretation for the opposite effect of the Gaussian disorder on the low and high-energy behavior, in the classical transport regime (dictated by the mean free path). To this end, we scrutinize the role of real-space distribution of electronic states subjected to a pseudomagnetic field. It has been shown that low-energy states are preferentially located inside the deformed region, while high-energy states correspond to states outside the bubbles 68 . Furthermore, the zeroth pLL induces sublattice polarization 22 , where the preferred sublattice depends on the deformation direction 6 . This sublattice polarization of the first Landau level is a unique feature of pseudomagnetic fields. The strain direction considered here causes the n = 0 pLL states in the central region (red color in left inset of  On the other hand, the states in the outer ring (blue color) with an opposite pseudomagnetic field sign, is A-sublattice polarized. Consequently, the low-energy states inside the bubbles tend to be confined by this opposite sublattice polarization from the normal conducting states outside the bubbles in the undeformed regions. Now, because we introduce long-range disorder randomly on both A and B sublattices everywhere in the sample, the increased disorder concentration disturbs the perfect low-energy sublattice polarization. Thus, the spatial sublattice-confinement of the states is broken, leading to an increased mean free path and conductivity. We note that this scenario is compatible with the fact that the quantum localization effects increase with disorder, as expected (see for instance Fig. 3a).
To further confirm the sublattice-polarization scenario, we introduce a sublattice-selective (vacancy) disorder, which has been shown to provoke asymmetric transport phenomena 69,70 . We model such unreconstructed vacancies by removing carbon atoms randomly everywhere in the sample. In Figure 5, the diffusion coefficients for sublattice-selective vacancy disorder, located either on the A-or B-sublattice, are shown, both with and without strain. The inequivalent effect of A-and B-vacancies exists only when the pseudomagnetic field is present. Such vacancies increase the number of states in the opposite sublattice from where they reside. As such, the A-vacancies increase the number of states in the Bsublattice. Because the low-energy states are mostly polarized on the B-sublattice, they are strongly affected by short-range valley-mixing effects induced by A-vacancies, leading to stronger localization effects (blue curve on Figure 5). Indeed, unlike the zeroth order Landau level caused by a real magnetic field where valleys are located on opposite sublattices, the valleys are concentrated on the same sublattice in the presence of a pseudomagnetic field 6 . Finally, even if the A-vacancies yield stronger localization effects, they concurrently lead to longer meanfree paths, when compared to B-vacancies. This observation further supports our interpretation since A-vacancies (increased B-sublattice polarization) increase the number of states to percolate through the blocking A-polarized ring around the bubbles.

VI. CONCLUSION
The presence of strain-induced pseudomagnetic fields forming pseudo-Landau levels has been shown to trigger unconventional quantum transport features compared to other types of disordered graphene systems. The transport physics close to half-filling is particularly remarkable, with a huge drop of the transport mean free path upon switching on a pseudomagnetic field, indicating the possibility of a mechanically induced metal-insulator transition.
We note, however, that despite of the weak localization caused by the deformation field, the strength of quantum interferences, even in the large disorder limit, remains anomalously weak for the sublattice polarized ze-roth pseudo Landau level (compared to other types of disorder). This anomalous scaling with disorder gives a peculiar transport fingerprint of the combination of pseudomagnetic field and disorder effects, especially for the sublattice-polarization unique to the pseudomagnetic field compared to real magnetic fields.
Superlattice networks of deformation fields have been experimentally demonstrated 71 , and could be tunable under pressure 27,72,73 or by temperature. In this way, varying the mean free path upon deformation, could reversibly switch the system between metallic and insulating states (for ξ/L 1, with ξ the localization length and L the sample size). Finally, the role of pseudomagnetic field in spin transport deserves further consideration, in the context of spin manipulation and graphene spin-based devices 74 .