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

We begin with the tight-binding Hamiltonian describing the hopping of electronic excitations in the lowest band (in units in which ${\rm{\hslash }}=e=1$ )

$$\begin{eqnarray}&&\hat{H}=\displaystyle \sum _{{\boldsymbol{n}}\in {\mathsf{A}}}\displaystyle \sum _{l=1}^{3}({t}_{{\boldsymbol{n}},l}\ \hat{b}{}_{{\boldsymbol{n}}+{{\boldsymbol{d}}}_{l}}^{\dagger }\ {\hat{a}}_{{\boldsymbol{n}}}+{t}_{{\boldsymbol{n}},l}^{\ast }\ \hat{a}{}_{{\boldsymbol{n}}}^{\dagger }\ {\hat{b}}_{{\boldsymbol{n}}+{{\boldsymbol{d}}}_{l}}).\end{eqnarray} \tag{ 1 }$$

The creation and annihilation operators $\hat{a}{}_{{\boldsymbol{n}}}^{\dagger },{\hat{a}}_{{\boldsymbol{n}}}$ and $\hat{b}{}_{{\boldsymbol{n}}+{{\boldsymbol{d}}}_{l}}^{\dagger },{\hat{b}}_{{\boldsymbol{n}}+{{\boldsymbol{d}}}_{l}}$ belong to the sublattices ${\mathsf{A}}$ and ${\mathsf{B}}$, respectively. Vector ${\boldsymbol{n}}$ runs over all points of the sublattice ${\mathsf{A}}$ and the three vectors ${{\boldsymbol{d}}}_{l}$ point along the links to the nearest neighbors in the sublattice ${\mathsf{B}}$. We choose ${{\boldsymbol{d}}}_{1}=\left(\tfrac{1}{2},\tfrac{\sqrt{3}}{2}\right){d}_{0}$, ${{\boldsymbol{d}}}_{2}=(-1,0){d}_{0}$, ${{\boldsymbol{d}}}_{3}=\left(\tfrac{1}{2},-\tfrac{\sqrt{3}}{2}\right){d}_{0}$, where d₀ is the interatomic distance in pristine graphene. The hopping parameters in pristine graphene are all equal, ${t}_{{\boldsymbol{n}},l}={t}_{0}$, and become position (${\boldsymbol{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\;\mathrm{eV}$ and distances in multiples of ${d}_{0}=0.142\;\mathrm{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⁴ $| {\psi }_{{\boldsymbol{n}}}\rangle$ for ${\boldsymbol{n}}\in {\mathsf{A}}\cup {\mathsf{B}}$ as

$$\begin{eqnarray}&&{\bf{H}}=\displaystyle \sum _{{\boldsymbol{n}}\in {\mathsf{A}}}\;\displaystyle \sum _{l=1}^{3}{t}_{{\boldsymbol{n}},l}\ | {\psi }_{{\boldsymbol{n}}+{{\boldsymbol{d}}}_{l}}^{{\mathsf{B}}}\rangle \langle {\psi }_{{\boldsymbol{n}}}^{{\mathsf{A}}}| +{\rm{h.c.}}\end{eqnarray} \tag{ 2 }$$

In regular graphene, where all ${t}_{{\boldsymbol{n}},l}={t}_{0}$, the spectrum of the Hamiltonian exhibits conical Dirac points ${{\boldsymbol{K}}}^{(s)}$ at which the dispersion relation is approximately linear $E({\boldsymbol{k}})\sim | {\boldsymbol{k}}-{{\boldsymbol{K}}}^{(s)}|$, where ${{\boldsymbol{K}}}^{(0)}=(0,4\pi /(3\sqrt{3}{d}_{0}))$ and ${K}^{(s)}$ are obtained by rotation of ${{\boldsymbol{K}}}_{0}$ by an angle $\tfrac{\pi }{3}s$ with $s=0,\ldots ,5$, as shown in figure 1 (see [37] for more fundamentals of graphene theory). Since only two of them are physically inequivalent we will choose the pair ${{\boldsymbol{K}}}^{(0)},{{\boldsymbol{K}}}^{(3)}$ and denote it ${{\boldsymbol{K}}}^{\pm }=(0,\pm 4\pi /(3\sqrt{3}{d}_{0}))$.

Zoom In Zoom Out Reset image size

In the long wave limit, the discrete tight-binding Hamiltonian can be approximated by the two-dimensional Dirac Hamiltonian [2]

$$\begin{eqnarray}&&{H}_{D}={\rm{i}}{v}_{F}{\sigma }^{l}({\partial }_{l}-{\rm{i}}{K}_{l}^{(s)})\end{eqnarray} \tag{ 3 }$$

separately for each of the Dirac points s. Here, ${v}_{F}=\tfrac{3}{2}{t}_{0}{d}_{0}$ is the Fermi velocity of the excited electrons and ${\sigma }^{l}$ are Pauli matrices. As long as regular planar graphene is considered, the constant ${{\boldsymbol{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)$ (see figure 2, top). We assume for simplicity that the carbon atoms are just lifted by the deformed surface in the perpendicular direction ${z}_{{\boldsymbol{n}}}=h({x}_{{\boldsymbol{n}}},{y}_{{\boldsymbol{n}}})$ while keeping their original $({x}_{{\boldsymbol{n}}},{y}_{{\boldsymbol{n}}})$ coordinates in the plane. This approximation is acceptable as long as the perturbations of the distances stay small⁵ , i.e. $\delta {{\boldsymbol{d}}}_{{\boldsymbol{n}},l}={\underline{{\boldsymbol{\varepsilon }}}}_{{\boldsymbol{n}}}{{\boldsymbol{d}}}_{{\boldsymbol{n}},l}\ll {d}_{0}$, where $\underline{{\boldsymbol{\varepsilon }}}$ 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 $\delta | {{\boldsymbol{d}}}_{{\boldsymbol{n}},l}| \approx {(\delta {z}_{{\boldsymbol{n}}}/{d}_{0})}^{2}/2\approx {({\partial }_{l}h)}^{2}/2$ where ${\partial }_{l}h$ is a directional derivative along the link ${{\boldsymbol{d}}}_{l}$.

Zoom In Zoom Out Reset image size

Accordingly, the hopping parameters change slightly to become ${\tilde{t}}_{{\boldsymbol{n}},l}={t}_{0}+\delta {t}_{{\boldsymbol{n}},l}$ and vary slowly over the lattice. Applying an empirical relation between the distances ${{\boldsymbol{d}}}_{{\boldsymbol{n}},l}$ and the hopping terms ${\tilde{t}}_{{\boldsymbol{n}},l}$ we get⁶

$$\begin{eqnarray}&&{\tilde{t}}_{{\boldsymbol{n}},l}={t}_{0}\mathrm{exp}(-\beta \;\delta | {{\boldsymbol{d}}}_{{\boldsymbol{n}},l}| /{d}_{0})\end{eqnarray} \tag{ 4 }$$

with $\beta =3.37$ [38–40]. For small length variations $\delta | {{\boldsymbol{d}}}_{{\boldsymbol{n}},l}| \ll {d}_{0}$ the latter can be approximated by the linear relation $\delta {t}_{{\boldsymbol{n}},l}/{t}_{0}=-\beta \delta | {{\boldsymbol{d}}}_{{\boldsymbol{n}},l}| /{d}_{0}$ which will allow us to relate $\delta {t}_{{\boldsymbol{n}},l}$ linearly to the strain tensor ${\underline{{\boldsymbol{\varepsilon }}}}_{{\boldsymbol{n}}}$.

It can be shown that the Hamiltonian (2) with such modified hopping parameters⁷ leads in the long wavelength limit to the continuous Hamiltonian [22, 23, 25, 41]

$$\begin{eqnarray}&&\tilde{{\bf{H}}}=\int {d}^{2}{\boldsymbol{x}}\sqrt{| \tilde{g}({\boldsymbol{x}})| }\;{\Psi }^{\dagger }({\boldsymbol{x}}){{iv}}_{F}{\sigma }^{a}{\tilde{e}}_{a}^{\ l}({\boldsymbol{x}})({\partial }_{l}-{\rm{i}}{\tilde{K}}_{l}^{(s)}({\boldsymbol{x}})+{\tilde{{\rm{\Gamma }}}}_{l}({\boldsymbol{x}}))\Psi ({\boldsymbol{x}})\end{eqnarray} \tag{ 5 }$$

resembling the one for the Dirac fermions in curved space.

The low energy electronic excitations would then satisfy an effective⁸ two-dimensional Dirac equation $i{\partial }_{t}{\rm{\Psi }}={\tilde{H}}_{D}{\rm{\Psi }}$ with the evolution generator given by

$$\begin{eqnarray}&&{\tilde{H}}_{D}={\rm{i}}{v}_{F}{\sigma }^{a}{\tilde{e}}_{a}^{\ l}({\boldsymbol{x}})({\partial }_{l}-{\rm{i}}{\tilde{K}}_{l}^{(s)}({\boldsymbol{x}})+{\tilde{{\rm{\Gamma }}}}_{l}({\boldsymbol{x}})).\end{eqnarray} \tag{ 6 }$$

Here, ${\tilde{{\boldsymbol{e}}}}_{a}({\boldsymbol{x}})$ (a = 1, 2) plays the role of a local frame (zweibein) and gives rise to an effective (inverse) metric

$$\begin{eqnarray}&&{\tilde{g}}^{{ij}}({\boldsymbol{x}})={\tilde{e}}_{a}^{i}({\boldsymbol{x}})\ {\tilde{e}}_{b}^{j}({\boldsymbol{x}})\ {\delta }^{{ab}}\end{eqnarray} \tag{ 7 }$$

describing an emergent curved geometry in which the electronic excitations propagate. ${\tilde{K}}_{l}^{(s)}({\boldsymbol{x}})$ is a vector potential whose curl gives rise to an effective pseudo-magnetic field

$$\begin{eqnarray}&&{\tilde{B}}^{(s)}({\boldsymbol{x}})={\rm{rot}}\tilde{{\boldsymbol{K}}}{}^{(s)}({\boldsymbol{x}})={(-1)}^{s}\tilde{B}({\boldsymbol{x}})\end{eqnarray} \tag{ 8 }$$

which in two dimensions has only one component. ${\tilde{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. ${\tilde{{\rm{\Gamma }}}}_{l}({\boldsymbol{x}})$ corresponds to the spin-connection and guarantees hermiticity of the above Hamiltonian when the frame ${\tilde{e}}_{a}^{i}({\boldsymbol{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 $\underline{\tilde{{\boldsymbol{\varepsilon }}}}$ which transforms the local frame⁹

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{e}}}}_{a}({\boldsymbol{x}})=({\bf{1}}-\underline{\tilde{{\boldsymbol{\varepsilon }}}}({\boldsymbol{x}})){{\boldsymbol{e}}}_{a},\end{eqnarray} \tag{ 9 }$$

defines the metric of the curved continuous space

$$\begin{eqnarray}&&{\tilde{g}}^{{ij}}({\boldsymbol{x}})={\delta }_{{ij}}-2\;{\tilde{{\boldsymbol{\varepsilon }}}}_{{ij}}({\boldsymbol{x}})\end{eqnarray} \tag{ 10 }$$

and the pseudo-magnetic vector potential

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{K}}}}^{(s)}({\boldsymbol{x}})={{\boldsymbol{K}}}^{(s)}+({-1)}^{s}\displaystyle \frac{\beta }{2}\left(-2{\tilde{{\boldsymbol{\varepsilon }}}}_{{xy}},{\tilde{{\boldsymbol{\varepsilon }}}}_{{yy}}-{\tilde{{\boldsymbol{\varepsilon }}}}_{{xx}}\right)\end{eqnarray} \tag{ 11 }$$

[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 ${\tilde{{\boldsymbol{K}}}}^{(s)}$ [44]. However, due to the assumed proportionality of $\delta {t}_{{\boldsymbol{n}},l}\sim \beta \delta | {{\boldsymbol{d}}}_{{\boldsymbol{n}},l}|$, the effective strain tensor $\underline{\tilde{{\boldsymbol{\varepsilon }}}}$ is proportional to the real strain applied to graphene, $\underline{\tilde{{\boldsymbol{\varepsilon }}}}=\beta \;\underline{{\boldsymbol{\varepsilon }}}$, 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 $\beta \gt 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 (see microlocal analysis of operators [45]). While the pseudo-magnetic vector potential $\tilde{{\boldsymbol{K}}}{}^{(s)}({\boldsymbol{x}})$ locally shifts the Dirac cones, the effective inverse metric locally deforms them in the momentum space (see figure 3). In consequence, the pseudo-momentum vectors ${\boldsymbol{k}}$ belonging to the deformed cones located at the shifted Dirac points satisfy locally

$$\begin{eqnarray}&&{\tilde{g}}^{{ij}}({\boldsymbol{x}})({k}_{i}-{\tilde{K}}_{i}{}^{(s)}({\boldsymbol{x}}))({k}_{j}-{\tilde{K}}_{j}{}^{(s)}({\boldsymbol{x}}))=0.\end{eqnarray} \tag{ 12 }$$

Zoom In Zoom Out Reset image size

In the continuous approximation the fields can be written in terms of the height function $h({\boldsymbol{x}})$ and its derivatives: in the lowest order the effective strain reads

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{\varepsilon }}}}_{{ij}}({\boldsymbol{x}})=\displaystyle \frac{\beta }{2}\;{\partial }_{i}h({\boldsymbol{x}}){\partial }_{j}h({\boldsymbol{x}})\end{eqnarray} \tag{ 13 }$$

and the pseudo-magnetic field $\tilde{B}({\boldsymbol{x}})={\rm{rot}}\tilde{{\boldsymbol{K}}}({\boldsymbol{x}})$ is

$$\begin{eqnarray}&&\tilde{B}({\boldsymbol{x}})={\partial }_{x}{\tilde{{\boldsymbol{\varepsilon }}}}_{{yy}}({\boldsymbol{x}})-{\partial }_{x}{\tilde{{\boldsymbol{\varepsilon }}}}_{{xx}}({\boldsymbol{x}})+2\;{\partial }_{y}{\tilde{{\boldsymbol{\varepsilon }}}}_{{xy}}({\boldsymbol{x}}).\end{eqnarray} \tag{ 14 }$$

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

$$\begin{eqnarray}&&\tilde{B}(r,\varphi )=\displaystyle \frac{\beta \;\mathrm{cos}(3\varphi )}{2}\;{h}^{\prime }(r)\left[h^{\prime \prime} (r)-\displaystyle \frac{1}{r}{h}^{\prime }(r)\right].\end{eqnarray} \tag{ 15 }$$

The angular prefactor $\mathrm{cos}(3\varphi )$ changes sign six times on the interval $(0,2\pi )$ and hence creates 6 zones of alternating sign of $\tilde{B}$ as can be seen in figure 2.

The definitions of both emergent fields, the metric and the pseudo-magnetic field, directly in terms of the change of the hopping parameters $\delta {t}_{{\boldsymbol{n}},l}$ are given in appendix A. While here, $\delta {t}_{{\boldsymbol{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 $\delta {t}_{{\boldsymbol{n}},l}$ can be of a different origin (see hopping of atoms in optical lattices [47–49]).

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¹⁰ . 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

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{v}^{i}}{{\rm{d}}\tau }+{\tilde{{\rm{\Gamma }}}}_{{kl}}^{i}{v}^{k}{v}^{l}=\sqrt{\tilde{g}}\ {\tilde{g}}^{{ij}}{\epsilon }_{{jk}}\;{v}^{k}{\tilde{B}}^{(s)}\end{eqnarray} \tag{ 16 }$$

where the 'velocity' ${v}^{i}(\tau )={\rm{d}}{x}^{i}(\tau )/{\rm{d}}\tau$. Here ${\tilde{{\rm{\Gamma }}}}_{{kl}}^{i}=\tfrac{1}{2}\;{\tilde{g}}^{{ij}}({\partial }_{k}{\tilde{g}}_{{jl}}+{\partial }_{l}{\tilde{g}}_{{kj}}-{\partial }_{j}{\tilde{g}}_{{kl}})$ denote the Christoffel symbols for the metric ${\tilde{g}}_{{ij}}$ and ${\epsilon }_{{ij}}$ is the completely anti-symmetric Levi-Civita symbol (${\epsilon }_{12}=-{\epsilon }_{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}-{\tilde{K}}_{j}{}^{(s)}$, with ${\boldsymbol{k}}$ being the wave vector satisfying ${v}_{F}| {\boldsymbol{k}}| =E$. Therefore, the initial velocity ${v}^{i}({\boldsymbol{x}}({\tau }_{0}))$ for a geodesic starting at the point ${\boldsymbol{x}}({\tau }_{0})$ should be chosen as

$$\begin{eqnarray}&&{v}^{i}({\boldsymbol{x}}({\tau }_{0}))={\tilde{g}}^{{ij}}({\boldsymbol{x}}({\tau }_{0}))\left[{k}_{j}-{\tilde{K}}_{j}{}^{(s)}({\boldsymbol{x}}({\tau }_{0}))\right],\end{eqnarray} \tag{ 17 }$$

i.e. different for each geodesic line.

If ${\tilde{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 ${\tilde{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, see figure 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 (see figure 2, bottom), bends the geodesics inwards and outwards in an alternating way. The two forces, the geometric ${\tilde{{\rm{\Gamma }}}}_{{kl}}^{i}$ and pseudo-magnetic $\tilde{B}$, are both proportional to the gradient of the strain tensor ${\varepsilon }_{{ij}}$ and hence to the product of first and second derivatives of the height function, symbolically $\sim \partial h\cdot {\partial }^{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 figure 4 for a typical bump geometry.

Zoom In Zoom Out Reset image size

In order to determine the current flow paths in deformed graphene we calculate numerically the local current density ${I}_{{\boldsymbol{n}},{\boldsymbol{m}}}$ between the neighboring lattice sites ${\boldsymbol{n}},{\boldsymbol{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–54], here we only summarize briefly the required equations.

Starting from the tight-binding Hamiltonian ${\bf{H}}$ of a deformed graphene ribbon (2), the matrix elements of the Green's function ${\bf{G}}$ are given by

$$\begin{eqnarray}&&{G}_{{\boldsymbol{n}},{\boldsymbol{m}}}=\langle {\psi }_{{\boldsymbol{n}}}| {({\bf{H}}-E{\bf{1}}-{\boldsymbol{\Sigma }})}^{-1}| {\psi }_{{\boldsymbol{m}}}\rangle .\end{eqnarray} \tag{ 18 }$$

E is the energy of the injected electrons and

$$\begin{eqnarray}&&{\boldsymbol{\Sigma }}=-{\rm{i}}\eta \displaystyle \sum _{{\boldsymbol{n}}\in {\mathsf{C}}}| {\psi }_{{\boldsymbol{n}}}\rangle \langle {\psi }_{{\boldsymbol{n}}}| ,\end{eqnarray} \tag{ 19 }$$

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 ${\mathsf{C}}$, see green sites in figure 2, and $\eta \gt 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 ${\mathsf{P}}$ (blue sites in figure 2, top), which we model by the inscattering function

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{{\mathsf{P}}}^{\mathrm{in}}=\nu \displaystyle \sum _{{\boldsymbol{n}},{\boldsymbol{m}}\in {\mathsf{P}}}{\chi }_{{\boldsymbol{n}}}^{\ast }\;{\chi }_{{\boldsymbol{m}}}| {\psi }_{{\boldsymbol{n}}}\rangle \langle {\psi }_{{\boldsymbol{m}}}| ,\end{eqnarray} \tag{ 20 }$$

where $\nu \gt 0$ is a constant and ${\chi }_{{\boldsymbol{n}}}$ are amplitudes encoding the initial energy and momentum of the injected plane wave at the contact ${\mathsf{P}}$. They are given by separate expressions on both sublattices ${\mathsf{A}},{\mathsf{B}}$, namely

$$\begin{eqnarray}{\chi }_{{\boldsymbol{n}}}=\left\{\begin{array}{ll}{C}_{-}{{\rm{e}}}^{{\rm{i}}({\boldsymbol{k}}+{{\boldsymbol{K}}}^{-}){\boldsymbol{n}}}+{C}_{+}{{\rm{e}}}^{{\rm{i}}({\boldsymbol{k}}+{{\boldsymbol{K}}}^{+}){\boldsymbol{n}}}, & {\boldsymbol{n}}\in {\mathsf{A}}\\ \sigma \;{C}_{-}{{\rm{e}}}^{{\rm{i}}({\boldsymbol{k}}+{{\boldsymbol{K}}}^{-}){\boldsymbol{n}}+{\rm{i}}\phi }-\sigma \;{C}_{+}{{\rm{e}}}^{{\rm{i}}({\boldsymbol{k}}+{{\boldsymbol{K}}}^{+}){\boldsymbol{n}}-{\rm{i}}\phi }, & {\boldsymbol{n}}\in {\mathsf{B}},\end{array}\right.\end{eqnarray} \tag{ 21 }$$

where $\phi =\mathrm{arctan}({k}_{x}/{k}_{y})$, $\sigma =\mathrm{sign}(E)$ and C_± are amplitudes of the excitations around the ${{\boldsymbol{K}}}^{\pm }$ valleys. This inscattering function corresponds to the injection of plane waves with momentum ${\boldsymbol{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 ${{\boldsymbol{\Sigma }}}_{{\mathsf{P}}}=0$ for the self-energy of the contact ${\mathsf{P}}$.

The local current of electrons, which originate from the contact ${\mathsf{P}}$ with energy E and which flow from atom ${\boldsymbol{n}}$ to the neighboring atom ${\boldsymbol{m}}$, is given by [56, 57]

$$\begin{eqnarray}&&{I}_{{\boldsymbol{n}},{\boldsymbol{m}}}^{{\mathsf{P}}}=\mathrm{Im}({t}_{{\boldsymbol{n}},{\boldsymbol{m}}}^{\ast }\;{{ \mathcal G }}_{{\boldsymbol{n}},{\boldsymbol{m}}}^{{\mathsf{P}}}),\end{eqnarray} \tag{ 22 }$$

(in the natural units $e\;{t}_{0}/h$) where

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal G }}}^{{\mathsf{P}}}={\bf{G}}\;{{\boldsymbol{\Sigma }}}_{{\mathsf{P}}}^{\mathrm{in}}\;{{\bf{G}}}^{\dagger }.\end{eqnarray} \tag{ 23 }$$

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 ${{\boldsymbol{K}}}^{(\pm )}$, 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.

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}\times {L}_{y}$ between 120 × 85 and 360 × 310 (multiples of d₀). The bump-like deformation (see figure 2, top)

$$\begin{eqnarray}&&h(r)=\displaystyle \frac{{h}_{0}}{1+{(r/{r}_{0})}^{2}}\end{eqnarray} \tag{ 24 }$$

is placed in the center of the ribbon with a spatial extension r₀ varying between 100 and 200 and the height h₀ varying between $0.5{r}_{0}$ and $1.5{r}_{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 $\lambda$ of the injected quantum current (see 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 (see figure 5). A simple estimation of the length scales leads to the gap¹¹ ${\rm{\Delta }}E=3\pi /L$ where $L=\mathrm{min}({L}_{x},{L}_{y})$. It reflects the absence of long waves with wavelengths $\lambda$ larger than the system size L and limits from below the range of admissible current energies E. In our systems ${\rm{\Delta }}E$ is between 0.03 and 0.11 and is kept well separated from the considered current energies.

Zoom In Zoom Out Reset image size

On the one hand, for small energies E the wavelengths $\lambda \approx \tfrac{3\pi }{E}$ must satisfy $\lambda \ll {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 $\lambda \gg {d}_{0}$ for the continuous approximation to hold. Further on, for the geometrical optics approximation to hold, the wavelengths $\lambda$ must be shorter than the scale of the geometric structures on which the wave is supposed to scatter which here means $\lambda \ll {r}_{0}$. This gives a hierarchy of length scales ${d}_{0}\ll \lambda \ll {r}_{0}\lt {L}_{x},{L}_{y}$ which must be satisfied in all numerical computations.

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 $\mathrm{sin}(\theta )=\lambda /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}^{\prime }\approx L\mathrm{sin}(\theta )=L\lambda /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\approx {D}^{\prime }\approx \sqrt{L\lambda }$ holds.

In order to obtain possibly narrow current flows we additionally give the injected 'beam' a Gaussian form ${| {\psi }^{\mathrm{in}}(x)| }^{2}\sim \mathrm{exp}(-4\mathrm{log}(2){(x-{x}_{0})}^{2}/{D}^{2})$ which is known for low spreading¹² and has half-width $D/2$ (see figure 6).

Zoom In Zoom Out Reset image size

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\gt 0$) behave as Dirac particles while the holes ($E\lt 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 $\pm {\boldsymbol{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¹³ . The contact formula (21) enables a simple filtering of valleys when the contact is placed parallel to the armchair edge and the injection momentum ${\boldsymbol{k}}$ is chosen in the zigzag direction. In such a case, (21) reduces to ${\chi }_{{\boldsymbol{n}}}={C}_{-}\pm {C}_{+}$ for ${\boldsymbol{n}}\in {\mathsf{A}},{\mathsf{B}}$, respectively (with $\sigma =+1$). By choosing ${\chi }_{{\boldsymbol{n}}\in {\mathsf{A}}}={\chi }_{{\boldsymbol{n}}\in {\mathsf{B}}}$ we ensure injection only at valley ${{\boldsymbol{K}}}^{-}$. The current flows then through the whole lattice mainly through that one channel—projections onto the corresponding subspaces lead to amplitude ratios ${{\boldsymbol{K}}}^{-}\;:{{\boldsymbol{K}}}^{+}\gt$ 10:1 at all studied energies.

We observed that the form of the current flowing via valley ${{\boldsymbol{K}}}^{-}$ is more compact while the one via valley ${{\boldsymbol{K}}}^{+}$ is clearly wider. At energies $E\gt 0.3$, the Dirac-cones become anisotropic, slightly triangular (see figure 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\lt 0.3$, the dispersion relation is almost round and both currents, via ${{\boldsymbol{K}}}^{\pm }$, flow (without curvature) almost identically.

A typical current flow around a bump, from a contact at the bottom edge towards an opposite boundary is visualized in figure 7. The maximal value of the pseudo-magnetic field, with spatial distribution as shown in figure 2 (bottom), is given by ${B}_{\mathrm{max}}=0.377\;{B}_{0}\;\beta \;{h}_{0}^{2}\;{r}_{0}^{-3}$ with ${B}_{0}=h/(e\;{d}_{0}^{2})=205\;100\;{\rm{T}}$. In the examples discussed below, this value varies in the range $100-1000\;{\rm{T}}$ which is consistent with $300\;{\rm{T}}$ estimated from the Landau levels for nanobubbles of similar size ($10\;\mathrm{nm}\approx 70\;{d}_{0}$) [8].

Zoom In Zoom Out Reset image size

In the next section, we discuss the influence of geometry on the current flows for various choices of parameters.

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 figure 8 the current is always injected at energy E = 0.2 while the amplitude h₀ 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.

Zoom In Zoom Out Reset image size

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. Figure 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.

Zoom In Zoom Out Reset image size

The last example, presented in figure 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\geqslant 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.

Zoom In Zoom Out Reset image size

The transformation ${\psi }_{{\mathsf{B}}}\to -{\psi }_{{\mathsf{B}}}$, which changes the sign of the wavefunction ψ on the whole sublattice ${\mathsf{B}}$ (while preserving the values of ${\psi }_{{\mathsf{A}}}$), is a symmetry of the hopping Hamiltonian (2). Under its action, ${\bf{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\lt 0$ from solutions with positive $E\gt 0$ in a simple way.

Performing numerical computations, we observed that the antiparticle (hole) current with $E=-{E}_{0}\lt 0$ follows the same path as the corresponding particle current with $E={E}_{0}\gt 0$ and all other conditions remain unchanged. The only difference is that the current (22) and (B5) changes its overall sign, i.e. flows in the opposite direction (see figure 11). This is a direct consequence of the transformation ${\psi }_{{\mathsf{B}}}\to -{\psi }_{{\mathsf{B}}}$ or, in other words, of the time-reversal symmetry.

Zoom In Zoom Out Reset image size

The time-reversal symmetry enforces also the transformation ${\boldsymbol{k}}\to -{\boldsymbol{k}}$ (while keeping the same propagation direction) which exchanges the valleys ${{\boldsymbol{K}}}^{-}\to {{\boldsymbol{K}}}^{+}$ and, in consequence, the sign of the pseudo-magnetic potential ${\boldsymbol{A}}\to -{\boldsymbol{A}}$. Therefore, the antiparticles (figure 11, right) 'feel' the opposite pseudo-magnetic field $-\tilde{B}$ and their current injected at the same contact and in the same direction follows the same path as that of particles (figure 11, left).

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 position dependent metric and pseudo-magnetic field would guide the electron waves through graphene (see figure 12) as has been discussed in [67, 68] in the case of pure metric lenses. One such setup is presented in figure 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size