Current flow paths in deformed graphene: from quantum transport to classical trajectories in curved space

In this work we compare two fundamentally different approaches to the electronic transport in deformed graphene: a) the condensed matter approach in which current flow paths are obtained by applying the non-equilibrium Green's function (NEGF) method to the tight-binding model with local strain, b) the general relativistic approach in which classical trajectories of relativistic point particles moving in a curved surface with a pseudo-magnetic field are calculated. The connection between the two is established in the long-wave limit via an effective Dirac Hamiltonian in curved space. Geometrical optics approximation, applied to focused current beams, allows us to directly compare the wave and the particle pictures. We obtain very good numerical agreement between the quantum and the classical approaches for a fairly wide set of parameters, improving with the increasing size of the system. The presented method offers an enormous reduction of complexity from irregular tight-binding Hamiltonians defined on large lattices to geometric language for curved continuous surfaces. It facilitates a comfortable and efficient tool for predicting electronic transport properties in graphene nanostructures with complicated geometries. Combination of the curvature and the pseudo-magnetic field paves the way to new interesting transport phenomena such as bending or focusing (lensing) of currents depending on the shape of the deformation. It can be applied in designing ultrasensitive sensors or in nanoelectronics.


I. INTRODUCTION
It is a well known fact in relativistic particle physics that electromagnetic and gravitational fields couple differently to particles and antiparticles -the first distinguishes their opposite charges, the second treats equivalently their identical masses. In graphene, both these phenomena can be simulated by the action of magnetic and pseudo-magnetic fields as well as the influence of curvature on the dynamics of particle-like excitations and holes. Continuous models of elastically deformed graphene, describing the excitations effectively by a twodimensional Dirac equation for massless fermions coupled to an artificial magnetic field, have been extensively studied in the literature in recent years, see e.g. [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16] for a representative selection or [17] for the most recent overview of the topic. Further geometric aspects, related to the Dirac equation in curved space, such as coupling to the effective metric leading to position and direction dependent Fermi velocity, have not yet attracted as much attention [18][19][20][21][22][23][24][25]. Especially, the electronic transport in the presence of curvature centers [24] still lacks detailed qualitative and quantitative results. The key question elaborated here, on the relation between the electric currents and classical trajectories, has been addressed before in the presence of pure pseudo-magnetic [26,27], real magnetic [28], combined pseudo-magnetic and real magnetic [29] or pure electric fields [30] without taking into account further influence from the curvature. How- * e-mail:stegmann@icf.unam.mx † e-mail:nikodem.szpak@uni-due.de ever, as we will show below, the metric effect can be as relevant as that of the pseudo-magnetic field.
In this work, we study the effect of curvature on the current flow lines in elastically deformed graphene sheets. Using a tight-binding model, we first apply the nonequilibrium Green's function (NEGF) method to obtain the current flow paths in the presence of a localized deformation. Then, utilizing the specific dispersion relation of graphene at low energies, we turn to the continuous approximation for long wavelengths, which leads us to the two-dimensional Dirac equation coupled to an effective pseudo-magnetic field and an attractive curvature. In order to visualize the action of both factors, we consider coherent current beams and apply the eikonal approximation. This simplifies the approach to the semiclassical particle picture in which the current lines turn out to be geodesics for relativistic charged massless particles moving in a curved two-dimensional surface in the presence of a magnetic field. We compare the numerically obtained current flows from the NEGF calculations with the classical geodesic lines. We find very good agreement, improving with the increasing size of the system, for a wide range of parameters of experimental interest. The combination of both emergent forces -curvature and the pseudo-magnetic field -paves the way to interesting transport phenomena such as bending or focusing (lensing) of currents depending on the shape of the deformation. It can find technological applications in designing ultrasensitive sensors or in nanoelectronics. The presented analogy facilitates a comfortable and efficient tool for calculating electronic transport properties in newly designed graphene nanostructures with complicated geometries. It offers an enormous reduction of complexity from hopping Hamiltonians defined on large deformed lattices to a semiclassical geometric language for the description of curvature effects in continuous surfaces. We conclude with the proposition of a geometrical lens for currents flowing through the deformed graphene. Similar electron optics phenomena in graphene have been proposed in [30][31][32][33] but the focusing was caused by p-n junctions and by the electric potential.
Here, we do not take into account either the electronelectron interaction [34], the recombination of electronic density due to deformation [35] or the relaxation of the lattice structure to its minimal elastic energy [36], leaving these topics open to further investigation.

II. THE EFFECTIVE DIRAC EQUATION IN CURVED SPACE
A. Tight-binding model for deformed graphene We begin with the tight-binding Hamiltonian describing the hopping of electronic excitations in the lowest band (in units in which = e = 1) The creation and annihilation operatorsâ † n ,â n and b † n+d l ,b n+d l belong to the sublattices A and B, respectively. Vector n runs over all points of the sublattice A and the three vectors d l point along the links to the nearest neighbors in the sublattice B. We choose The hopping parameters in pristine graphene are all equal, t n,l = t 0 , and become position (n) and direction (l) dependent in deformed graphene which reflects modified tunneling probabilities between neighboring atoms. In the following, energies are measured in multiples of t 0 = 2.8 eV and distances in multiples of d 0 = 0.142 nm, if not given explicitly.
Since we do not take into account any interaction effects between the electrons, the second quantized Hamiltonian (1) effectively reduces to a one-particle sector and can be written in a basis of localized states 1 |ψ n for n ∈ A ∪ B as In regular graphene, where all t n,l = t 0 , the spectrum of the Hamiltonian exhibits conical Dirac points K (s) at which the dispersion relation is approximately linear and K (s) are obtained by rotation of K 0 by an angle π 3 s with s = 0, ..., 5, as shown in Fig 1 (cf. [37] for more fundamentals of graphene theory). Since only two of them are physically inequivalent we will choose the pair K (0) , K (3) and denote it K ± = 0, ±4π/(3 √ 3d 0 ) . In the long wave limit, the discrete tight-binding Hamiltonian can be approximated by the twodimensional Dirac Hamiltonian [2] separately for each of the Dirac points s. Here, v F = 3 2 t 0 d 0 is the Fermi velocity of the excited electrons and σ l are Pauli matrices. As long as regular planar graphene is considered, the constant K (s) vectors can be gauged away by multiplication of the wavefunction with an appropriate phase.
In this work we concentrate on small out-of-plane perturbations described by the height function h(x, y) (cf. Fig 2, top). We assume for simplicity that the carbon atoms are just lifted by the deformed surface in the perpendicular direction z n = h(x n , y n ) while keeping their original (x n , y n ) coordinates in the plane. This approximation is acceptable as long as the perturbations of the distances stay small 2 , i.e. δd n,l = ε n d n,l d 0 , where ε is the strain tensor describing deformation of the graphene lattice [38]. Knowing the positions of the atoms the perturbations of their distances can be directly calculated via the Euclidean formula and for small deformations approximated by δ|d n,l | ≈ (δz n /d 0 ) 2 /2 ≈ (∂ l h) 2 /2 where ∂ l h is a directional derivative along the link d l .
Accordingly, the hopping parameters change slightly to become t n,l = t 0 + δt n,l and vary slowly over the lattice. Applying an empirical relation between the distances d n,l with β = 3.37 [38][39][40]. For small length variations δ|d n,l | d 0 the latter can be approximated by the linear relation δt n,l /t 0 = −β δ|d n,l |/d 0 which will allow us to relate δt n,l linearly to the strain tensor ε n .
It can be shown that the Hamiltonian (2) with such modified hopping parameters 4 leads in the long wave-length limit to the continuous Hamiltonian [22,23,25,41] resembling the one for the Dirac fermions in curved space. The low energy electronic excitations would then satisfy an effective 5 two-dimensional Dirac equation i∂ t Ψ = H D Ψ with the evolution generator given by Here, e a (x) (a = 1, 2) plays the role of a local frame (zweibein) and gives rise to an effective (inverse) metric describing an emergent curved geometry in which the electronic excitations propagate. K l (x) is a vector potential whose curl gives rise to an effective pseudomagnetic field which in two dimensions has only one component. B (s) acts oppositely on the excitations in different valleys s (hence the prefix "pseudo") and has therefore to be distinguished from a true magnetic field which would break the time-reversal symmetry of the system. The latter situation is excluded here because the perturbations of the hopping parameters in (2) are purely real. Γ l (x) corresponds to the spin-connection and guarantees hermiticity of the above Hamiltonian when the frame e i a (x) is position-dependent [22,23]. In contrast to the vector potential and the frame perturbations, both proportional to the strain, spin-connection is proportional to its gradient. It has no classical counterpart but in quantum scattering on steep surface deformations containing closed geodesics it can contribute to transmission resonances at wavelengths comparable with the size of the deformation [42,43]. For shallow deformations varying smoothly over the surface, such as considered here, the spin-connection becomes subdominant and can be skipped.
Both remaining fields are continuous extrapolations of their discrete counterparts on the lattice, given in Appendix A, and their values can be related to an abstract strain tensor ε which transforms the local frame 6 5 All effective fields will be distinguished by a tilde. 6 There is a small discrepancy between the frame transformation given in (9) and e i a interpreted as local Fermi velocity in [22,23] by a strain-trace term ε i i . It is related to different normalizations of the wavefunction Ψ(x). Since on the lattice nψ n ψn = 1 in curved space Ψ Ψ √ g d 2 x = 1 which forces us to choose the discretization rule ψn = g 1/4 (xn)Ψ(xn) here. defines the metric of the curved continuous space and the pseudo-magnetic vector potential [2, 18,38]. It is not a priori known what the effective geometry will look like since, for the given curved surface of graphene, the effective metric is obtained by the expansion of the Hamiltonian (2) around the Dirac points shifted by K (s) [44]. However, due to the assumed proportionality of δt n,l ∼ β δ|d n,l |, the effective strain tensor ε is proportional to the real strain applied to graphene, ε = β ε, and thus the effective geometry for the electronic excitations is nearly identical to the real geometry of the deformed graphene sheet but the deformation is magnified by the factor β > 1.
When the graphene lattice is perturbed the Dirac cones do not have a global meaning any more. Instead, everywhere on the lattice local cones can be associated to the Dirac operator in (6) along which wavefronts of local perturbations would propagate (cf. microlocal analysis of operators [45]). While the pseudo-magnetic vector potential K (s) (x) locally shifts the Dirac cones, the effective inverse metric locally deforms them in the momentum space (see Fig. 3). In consequence, the pseudomomentum vectors k belonging to the deformed cones located at the shifted Dirac points satisfy locally In the continuous approximation the fields can be written in terms of the height function h(x) and its derivatives: in the lowest order the effective strain reads and the pseudo-magnetic field B(x) = rot K(x) is In the examples discussed below we will consider rotationally symmetric elastic deformations of graphene for which the height function h(x, y) = h(r) uniquely parameterizes the lattice geometry. Such geometries seem to be well under experimental control as they appear when a localized perpendicular force is applied to the surface, e.g. from an STM-head [46], and forms a tip-like deformation. Then, the last formula simplifies (in polar coordinates) to The angular prefactor cos(3ϕ) changes sign six times on the interval (0, 2π) and hence creates 6 zones of alternating sign of B as can be seen in Fig. 2.
The definitions of both emergent fields, the metric and the pseudo-magnetic field, directly in terms of the change of the hopping parameters δt n,l are given in Appendix A. While here, δt n,l arise in response to an elastic deformation of the graphene structure, the developed formalism can be applied to any lattice system described by the Hamiltonian (2) in which δt n,l can be of a different origin (cf. hopping of atoms in optical lattices [47][48][49]).

B. Geodesic lines in curved continuous geometry
In order to develop an efficient way of predicting the current flow paths in deformed graphene we will make use of the geometrical optics (eikonal approximation) in which the propagation of waves is replaced by the tracing of rays. The latter satisfy "equations of motion" typical for point-like particles. Eikonal approximation applied to the Dirac equation in curved space leads to the Mathisson-Papapetrou equation of motion describing a spinning particle in curved space [50].
Since the issue of spin in graphene is quite delicate we postpone its treatment to future work and ignore the spin degree of freedom here 7 . Then, by squaring the Dirac Hamiltonian and applying the eikonal approximation, we arrive at the geodesic equation for null geodesics (due to massless fermions) in a given curved 2D-surface. Note that in a static, i.e. time-independent, 2+1D geometry the full spacetime 2+1D geodesics match with the shortest (extreme) lines of the spatial 2D geometry. Therefore, it is sufficient to solve the geodesic equation for the 2D curved surface where the "velocity" Here, with spin we mean the pseudo-spin appearing effectively due to the symmetries of the honeycomb lattice [51]. The real spin of electrons is usually ignored in graphene anyway.
bols for the metric g ij and ij is the completely antisymmetric Levi-Civita symbol ( 12 = − 21 = 1). We are interested in a family of geodesics starting parallel to each other and representing a current injected at the contact with the initial momentum p i = k i − K j (s) , with k being the wave vector satisfying v F |k| = E. Therefore, the initial velocity v i (x(τ 0 )) for a geodesic starting at the point x(τ 0 ) should be chosen as i.e. different for each geodesic line. If B (s) is negligibly small the solutions do not depend on the value of the initial velocity (i.e. depend only on its direction) nor on the initial energy (also not on its sign) and hence particles and antiparticles follow the same lines. The presence of the magnetic field B (s) breaks this symmetry and deflects particles and antiparticles in opposite directions. Also, the initial value of the velocity (energy) starts to play a role -the lower it is the more the magnetic field has influence on the trajectory.
The curvature of a typical bump is clearly positive in the middle (and only slightly negative outside, cf. Fig.  2, top) and hence it bends the geodesic lines on the surface towards the center -as attractive gravitational forces do. The pseudo-magnetic field, due to its six-fold symmetry (cf. Fig. 2, bottom), bends the geodesics inwards and outwards in an alternating way. The two forces, the geometric Γ i kl and pseudo-magnetic B, are both proportional to the gradient of the strain tensor ε ij and hence to the product of first and second derivatives of the height function, symbolically ∼ ∂h · ∂ 2 h. Hence, both contributions are of the same order of magnitude so the correct prediction of electric currents in deformed graphene will require taking both of them into account. These two contributions to the bending of trajectories are compared in Fig. 4 for a typical bump geometry.

A. The NEGF method
In order to determine the current flow paths in deformed graphene we calculate numerically the local current density I n,m between the neighboring lattice sites n, m and plot its integral lines. To find the current in the presence of a stationary source, injecting electrons into the ribbon at a constant pace, we apply the non-equilibrium Green's function method (NEGF). As this method has been described in various textbooks, see e.g. [52][53][54], here we only summarize briefly the required equations.
Starting from the tight-binding Hamiltonian H of a deformed graphene ribbon (2), the matrix elements of the Green's function G are given by E is the energy of the injected electrons and is an imaginary self-energy by means of which we introduce absorbing boundaries [55], mimicking infinite dimensions of the graphene sheet. (The sum runs over the edge atoms C, cf. green sites in Fig. 2, and η > 0 is a constant.) These boundaries absorb impinging particles and suppress finite system size effects, such as standing waves between the boundaries of the system. The electrons are injected in the graphene ribbon by a contact P (blue sites in Fig. 2, top), which we model by the inscattering function where ν > 0 is a constant and χ n are amplitudes encoding the initial energy and momentum of the injected plane wave at the contact P. They are given by separate expressions on both sublattices A, B, namely where φ = arctan(k x /k y ), σ = sign(E) and C ± are amplitudes of the excitations around the K ± valleys. This inscattering function corresponds to the injection of plane waves with momentum k and energy E. Since we consider idealized contacts the injection of electrons does not affect their propagation in the graphene ribbon and thus we take Σ P = 0 for the self-energy of the contact P.
The local current of electrons, which originate from the contact P with energy E and which flow from atom n to the neighboring atom m, is given by [56,57] I P n,m = Im t * n,m G P n,m , (in the natural units e t 0 /h) where is the correlation function of the injected electrons. Interestingly, the lattice current formula (22) can also be derived as discretization of the Dirac current in curved space, see Appendix B.
Since, in general situations, we deal in graphene with contributions from both inequivalent Dirac points K (±) , interference effects in the current (due to its quadratic dependence on the wavefunctions) are expected. In order to eliminate the highly oscillatory behavior on the lattice scale from the plots, we average the numerical values of the current over the hexagons. We leave the examination of the interference patterns to a future work.

B. System parameters
We consider rectangular graphene ribbons of size varying from 80×50 to 240×180 carbon rings in the armchair and zigzag direction, respectively. This corresponds to sizes L x × L y between 120 × 85 and 360 × 310 (multiples of d 0 ). The bump-like deformation (cf. Fig 2, top) h is placed in the center of the ribbon with a spatial extension r 0 varying between 100 and 200 and the height h 0 varying between 0.5r 0 and 1.5r 0 . The ratio h 0 /r 0 controls the steepness of the bump and thus the maximal strain which should not exceed 10% for the continuous model to hold [6].
We choose the electric current to flow along the zigzag direction, from a contact placed on the armchair edge. Its width D is adjusted to the wavelength λ of the injected quantum current (cf. discussion below) and takes values between 10% and 30% of the edge length L x .
In our numerical simulations we can consider only graphene sheets of finite size (nanoribbons). Consequently, the available energy spectrum is discrete. However, for large ribbons the separation of the discrete energies shrinks to zero and reaches the continuous band structure in the limit. To test the granularity, it is instructive to plot the density of states and check the behavior near E = 0 which is crucial for our approximations (cf. Fig 5). A simple estimation of the length scales leads to the gap 8 ∆E = 3π/L where L = min(L x , L y ). It reflects the absence of long waves with wavelengths λ larger 8 The isolated states at E = 0 visible in Fig. 5 correspond to dispersionless zz-edge states present in graphene nanoribbons which are known to not contribute to the electronic transport [2,58].
than the system size L and limits from below the range of admissible current energies E. In our systems ∆E is between 0.03 and 0.11 and is kept well separated from the considered current energies. On the one hand, for small energies E the wavelengths λ ≈ 3π E must satisfy λ L x , L y in order to be compatible with the size of the ribbon. On the other hand, they must be much larger than the interactomic distance λ d 0 for the continuous approximation to hold. Further on, for the geometrical optics approximation to hold, the wavelengths λ must be shorter than the scale of the geometric structures on which the wave is supposed to scatter which here means λ r 0 . This gives a hierarchy of length scales d 0 λ r 0 < L x , L y which must be satisfied in all numerical computations.

C. Propagation of plane waves
In order to see the electric current flowing along the geodesics in the curved geometry of graphene it is necessary to stay close to the regime of plane waves propagating through the ribbon. While exact plane waves would require infinitely wide regions, in finite graphene ribbons we will have to deal with disturbing boundary effects. As a solution, we choose a finite contact, in the middle of one boundary and well separated from others, at which the wave will be injected locally as plane. In a sense, it corresponds to a single-slit diffraction: a hypothetical plane wave comes from outside the ribbon and entering the ribbon gets diffracted at the slit (here contact). For narrow contacts the wave will naturally spread across the ribbon. Contacts which are wider than the wavelength λ produce interference effects at angles sin(θ) = λ/D where D is the contact width. Hence, at the opposite end of the ribbon the width of the "beam" (measured between the two first interference minima) is D ≈ L sin(θ) = Lλ/D, where L is the length of the ribbon. Since we want the current flow to have approximately constant width during its propagation we choose the parameters so that D ≈ D ≈ √ Lλ holds. In order to obtain possibly narrow current flows we additionally give the injected "beam" a Gaussian form |ψ in (x)| 2 ∼ exp(−4 log(2) (x − x 0 ) 2 /D 2 ) which is known for low spreading 9 and has half-width D/2 (cf. Fig. 6).

D. Pseudo-magnetic field
The action of geometry is fundamentally different from the action of a magnetic field when applied to particles and antiparticles. Here, the electronic excitations above the Fermi level (E > 0) behave as Dirac particles while the holes (E < 0) behave as Dirac antiparticles with opposite charge. Both should react identically to curvature and oppositely to the pseudo-magnetic field.
There is, however, one problem in graphene: the number of Dirac excitations is doubled due to the existence of two inequivalent Dirac cones in the dispersion relation on the hexagonal lattice. In deformed graphene, each of these two kinds "feels" the opposite sign of the pseudo-magnetic potential ±A, which globally reflects the time-reversal symmetry present in the system.
The way out of this symmetric catch, implemented in our calculations, is an asymmetric injection of electrons at both valleys 10 . The contact formula (21) enables a simple filtering of valleys when the contact is placed parallel to the armchair edge and the injection momentum k is chosen in the zigzag direction. In such a case, (21) reduces to χ n = C − ±C + for n ∈ A, B, respectively (with σ = +1). By choosing χ n∈A = χ n∈B we ensure injection only at valley K − . The current flows then through the whole lattice mainly through that one channel -projections onto the corresponding subspaces lead to amplitude ratios K − : K + > 10 : 1 at all studied energies.
We observed that the form of the current flowing via valley K − is more compact while the one via valley K + is clearly wider. At energies E > 0.3, the Dirac-cones become anisotropic, slightly triangular (cf. Fig. 3, left). The flattened part helps the current to flow in one direction while the triangle-edge disperses the current away from the main direction. At lower energies, E < 0.3, the dispersion relation is almost round and both currents, via K ± , flow (without curvature) almost identically. FIG. 7: Electric current flow in curved graphene. Current streamlines (yellow arrows) and current density (red color shading) calculated by the NEGF method compared to classical geodesics (blue solid lines) for massless charged particles moving in a magnetic field on the continuous curved surface.
A typical current flow around a bump, from a contact at the bottom edge towards an opposite boundary is visualized in Fig. 7. The maximal value of the pseudomagnetic field, with spatial distribution as shown in is consistent with 300 T estimated from the Landau levels for nanobubbles of similar size (10 nm ≈ 70 d 0 ) [8].
In the next Section, we discuss the influence of geometry on the current flows for various choices of parameters.

A. Bending of current lines
In the following examples we demonstrate the influence of different geometric factors on the shape of the current flow and its agreement with classically predicted geodesic lines. In Fig. 8 the current is always injected at energy E = 0.2 while the amplitude h 0 of the central deformation varies between 0.5 and 1.0 times the size of the deformation given by r 0 = 150. The shape of the contact is chosen to be optimal according to the diffraction formula given above to keep the width of the current narrow along its whole path.
However, the effect of crossing of geodesics helps to focus the current and leads to a narrower flow than what can be expected from the standard diffraction. Fig. 9 shows two such examples at E = 0.2 and E = 0.3 where the contact width has been reduced by half and the current keeps its narrow width along the whole path while without the geodesic focusing a widening by factor 4 would be expected at the upper boundary. The enhanced focusing takes place at the cost of slight disagreement appearing between the flow and the geodesics. Its origin lays in the breakdown of the eikonal approximation relating waves to the current lines at caustics, i.e. where geodesics cross.
The last example, presented in Fig. 10, compares flows of current injected at different energies E = 0.3, 0.5 and 0.7 in the same geometry. In each case, optimal contact width has been chosen in an energy dependent way. With increasing energy the influence of the pseudo-magnetic field becomes less relevant and the flow becomes visibly straighter. The anisotropy of the Dirac cones at E ≥ 0.3 leads for each valley to propagation in three dominant directions [66]. Only the one pointing straight up contributes to the flow from the source located at the bottom edge (the other two point downwards). This explains the mismatch with geodesics for which an isotropic propagation is assumed. The transformation ψ B → −ψ B , which changes the sign of the wavefunction ψ on the whole sublattice B (while preserving the values of ψ A ), is a symmetry of the hopping Hamiltonian (2). Under its action, H changes its overall sign, independently of whether the hopping terms are all equal or perturbed. This enables us to generate solutions with negative energies E < 0 from solutions with positive E > 0 in a simple way.
Performing numerical computations, we observed that the antiparticle (hole) current with E = −E 0 < 0 follows the same path as the corresponding particle current with E = E 0 > 0 and all other conditions remain unchanged. The only difference is that the current (22), (B5) changes its overall sign, i.e. flows in the opposite direction (cf. Fig. 11). This is a direct consequence of the transformation ψ B → −ψ B or, in other words, of the time-reversal symmetry.
FIG. 11: Currents of particles with E = E 0 > 0 (left) and antiparticles with E = −E 0 < 0 (right) follow identical paths as they react to the pseudo-magnetic fields B of opposite signs and equal curvature. Since for antiparticles the current I is anti-parallel to the pseudo-momentum k the flow can be also interpreted as a flow of particles moving in the opposite direction.
The time-reversal symmetry enforces also the transformation k → −k (while keeping the same propagation direction) which exchanges the valleys K − → K + and, in consequence, the sign of the pseudo-magnetic potential A → −A. Therefore, the antiparticles (Fig. 11, right) "feel" the opposite pseudo-magnetic field − B and their current injected at the same contact and in the same direction follows the same path as that of particles (Fig.  11, left).
FIG. 12: Geometric lensing: geodesic lines for curved geometry and pseudo-magnetic field (blue solid lines) compared to pure geometry without pseudo-magnetic field (red dotted lines, background shading as in Fig 4). The pseudo-magnetic field shifts the focus away from behind the bump, creating two symmetric foci.

C. Geometric lensing of currents
One of the most intriguing applications of the presented ideas could be a purpose-built geometric lens that is able to deliberately focus or deflect electric current by an elastic deformation of the graphene surface. The principle would be similar to the gradient-index lenses where the position dependent refractive index guides electromagnetic waves through the medium. Here, the posi- tion dependent metric and pseudo-magnetic field would guide the electron waves through graphene (cf. Fig. 12) as has been discussed in [67,68] in the case of pure metric lenses. One such setup is presented in Fig. 13. Note that the presence of the pseudo-magnetic field prevents the trajectories from crossing directly behind the bump. For the same reason closed geodesics encircling the bump and the corresponding quantum scattering effects [42,43] are not present here. Such a device might play a role as an ultra-sensitive pressure sensor built of graphene which would direct the electric current from a source contact (at the bottom edge) to one of several drain contacts (at the top edge) depending on the deformation of the surface. The differences between voltages measured at different contacts would provide information on the amplitude of the deformation. We leave elaboration of this idea to a forthcoming publication.

V. DISCUSSION AND OUTLOOK
In this work we compared two fundamentally different approaches to the electronic transport in deformed graphene. Firstly, using the NEGF method we computed the quantum currents (22) of electronic excitations directly from the hopping model (2) with local modifications of the hopping parameters (4) due to strain. Secondly, we integrated classical trajectories (16) for relativistic charged massless particles moving in a curved 2-dimensional surface z = h(x, y) in the presence of an emergent pseudo-magnetic field (15). The connection between the two approaches has been established via an effective 2-dimensional Dirac Hamiltonian in curved space (6) appearing in the long-wave limit from the special dispersion relation of graphene at low energies. By applying geometrical optics approximation to focused plane wave beams we were able to switch from the wave to the particle picture in which the current lines are represented by the above mentioned classical trajectories. We obtained very good numerical agreement between these pictures for a fairly wide set of parameters. It has also been confirmed in time-dependent numerical simulations, similar to those in [27], in which wave-packets of finite size have been sent through the deformed lattice. Their propagation has been consistent with the calculated classical trajectories [69]. The presented analogue model is based on approximations which are valid for the following hierarchy of scales lattice constant wavelength deformation scale < system dimensions accompanied by the narrow-flow condition optimal contact width ∼ system size × wavelength which are satisfied for most real systems of interest 11 . Generally, the precision of the proposed approximations increases with the system size.
Clearly, the focus of the current paper has been put on the analogy between the distorted graphene and the Dirac equation in curved 2-dimensional space. This is why coherently injected plane-wave currents were of special interest. However, beside the simulation of quantum fields in curved spaces, the presented analogy has the potential of becoming an alternative, comfortable and efficient tool for calculating the electronic properties of newly designed graphene nanostructures. It offers an enormous reduction of the complexity from irregular hopping Hamiltonians defined on large hexagonal lattices to the semiclassical geometric language for the description of curvature effects in a continuous surface.
For the proposed applications, the following extensions of the presented setup seem to be most natural. Instead of injection at one contact, the current would rather flow between two or more defined contacts and take real boundaries into account. An extended model of contacts can be included, allowing for injections of electrons at various grades of coherence. A semiclassical description of the contact and boundary properties as well as of the interference effects for mixed currents flowing through both valleys, K and K , would be required. Also the significance of corrections from the electron-electron interactions could be taken into account.
Summarizing, the main results of the current work include the demonstration of the applicability of the continuous approximation in the prediction of current flow paths, verification of the relative significance of the pseudo-magnetic and geometric effects and the observation of the geometric lensing of currents in elastically deformed graphene. (A3)