Modelling electron-phonon interactions in graphene with curved space hydrodynamics

TL;DR

This paper presents a novel approach to modeling electron-phonon interactions in graphene using curved space hydrodynamics, capturing lattice vibrations as curvature effects influencing electronic flow.

Contribution

It introduces a new geometric framework for electron-phonon interactions in graphene, linking lattice vibrations to curvature in a hydrodynamic model.

Findings

Successfully reproduces linear resistivity scaling at high temperatures

Provides a temperature-independent coupling constant model

Offers a new perspective applicable to other Fermi liquid materials

Abstract

We introduce a different perspective describing electron-phonon interactions in graphene based on curved space hydrodynamics. Interactions of phonons with charge carriers increase the electrical resistivity of the material. Our approach captures the lattice vibrations as curvature changes in the space through which electrons move following hydrodynamic equations. In this picture, inertial corrections to the electronic flow arise naturally effectively producing electron-phonon interactions. The strength of the interaction is controlled by a coupling constant, which is temperature independent. We apply this model to graphene and recover satisfactorily the linear scaling law for the resistivity that is expected at high temperatures. Our findings open up a new perspective of treating electron-phonon interactions in graphene, and also in other materials where electrons can be described by the Fermi liquid theory.

Figures

Equations

$${\Gamma }_{\mu \nu }^{\sigma }=\frac{1}{2}{g}^{\sigma \rho }\left(\frac{\partial g_{\rho \mu }}{\partial x^{\nu }}+\frac{\partial g_{\rho \nu }}{\partial x^{\mu }}-\frac{\partial g_{\mu \nu }}{\partial x^{\rho }}\right).$$

$${\partial }_{\mu }{N}^{\mu }=\alpha F_N,$$

$${\partial }_{\mu }{T}_0^{\mu \nu }=\alpha F_T^{\nu },$$

$$F_N=-{\Gamma }_{\mu \lambda }^{\mu }{N}^{\lambda }$$

$$F_T^{\nu }=-{\Gamma }_{\mu \lambda }^{\mu }{T}_0^{\lambda \nu }-{\Gamma }_{\mu \lambda }^{\nu }{T}_0^{\mu \lambda }-({\eta }^{\mu \nu }-g^{\mu \nu }){\partial }_{\mu }p.$$

$$g_{\mu \nu }=\frac{\partial {\Phi }^{\alpha }(ct,x,y)}{\partial x^{\mu }}\frac{\partial {\Phi }^{\beta }(ct,x,y)}{\partial x^{\nu }}{\eta }_{\alpha \beta }.$$

$$\kappa (T)=c_{\nu }n(\frac{\mu }{k_{B}T}{)}^{\frac{3}{2}},$$

$$I(t)={\int }_0^{L_y}en(ct,x=L_x,y)U^x(ct,x=L_x,y)dy.$$

$$\rho =\gamma T\equiv \frac{\pi D^2{k}_B}{4e^2\hslash {\rho }_m{v}_{ph}^2{v}_F^2}T,$$

$${\partial }_{\mu }{N}^{\mu }=0,$$

$${\partial }_{\mu }{T}_0^{\mu \nu }=0,$$

$$T_0^{\mu \nu }=(ϵ+p)\frac{U^{\mu }{U}^{\nu }}{v_F^2}-p{\eta }^{\mu \nu }+{\pi }^{\mu \nu }.$$

$${\nabla }_{\mu }{N}^{\mu }=0,$$

$${\nabla }_{\mu }{T}^{\mu \nu }=0,$$

$${\nabla }_{\mu }{N}^{\mu }={\partial }_{\mu }{N}^{\mu }+{\Gamma }_{\mu \lambda }^{\mu }{N}^{\lambda },$$

$${\nabla }_{\mu }{T}^{\mu \nu }={\partial }_{\mu }{T}^{\mu \nu }+{\Gamma }_{\mu \lambda }^{\mu }{T}^{\lambda \nu }+{\Gamma }_{\mu \lambda }^{\nu }{T}^{\mu \lambda },$$

$${\Gamma }_{\mu \nu }^{\sigma }=\frac{1}{2}{g}^{\sigma \rho }\left(\frac{\partial g_{\rho \mu }}{\partial x^{\nu }}+\frac{\partial g_{\rho \nu }}{\partial x^{\mu }}-\frac{\partial g_{\mu \nu }}{\partial x^{\rho }}\right).$$

$$T^{\mu \nu }=(ϵ+p)\frac{U^{\mu }{U}^{\nu }}{v_F^2}-p{g}^{\mu \nu }+{\pi }^{\mu \nu },$$

$${\nabla }_{\mu }{T}^{\mu \nu }={\nabla }_{\mu }{T}_0^{\mu \nu }+({\eta }^{\mu \nu }-g^{\mu \nu }){\nabla }_{\mu }p.$$

$$0$$

$$={\partial }_{\mu }{T}_0^{\mu \nu }+{\Gamma }_{\mu \lambda }^{\mu }{T}_0^{\lambda \nu }+{\Gamma }_{\mu \lambda }^{\nu }{T}_0^{\mu \lambda }+({\eta }^{\mu \nu }-g^{\mu \nu }){\partial }_{\mu }p,$$

$${\partial }_{\mu }{N}^{\mu }=\alpha F_N,$$

$${\partial }_{\mu }{T}_0^{\mu \nu }=\alpha F_T^{\nu },$$

$$F_N=-{\Gamma }_{\mu \lambda }^{\mu }{N}^{\lambda }$$

$$F_T^{\nu }=-{\Gamma }_{\mu \lambda }^{\mu }{T}_0^{\lambda \nu }-{\Gamma }_{\mu \lambda }^{\nu }{T}_0^{\mu \lambda }-({\eta }^{\mu \nu }-g^{\mu \nu }){\partial }_{\mu }p.$$

$$\frac{1}{x_0^3}{\partial }_{\mu }^{\prime }{N}^{\prime \mu }=\frac{1}{x_0^3}{F}_N^{\prime },$$

$$\frac{{ϵ}_0}{x_0}{\partial }_{\mu }^{\prime }{T}^{\prime \mu \nu (eq)}$$

$$=\frac{{ϵ}_0}{x_0}{F}_T^{\prime \nu }+\frac{e_0{E}_0}{x_0^2}{n}^{\prime }{e}^{\prime }{E}^{\prime \nu }.$$

$$A_1=\frac{{\kappa }_0{u}_0}{x_0{ϵ}_0},$$

$$A_2=\frac{e_0{E}_0}{x_0{ϵ}_0}.$$

Full text

At finite temperatures, phonons interact with charge carriers, and therefore, contribute to the electrical resistivity of the respective material Rössler (2009). In suspended graphene, the temperature dependence of the resistivity follows a linear increase due to electron-phonon interactions above 100 K.

111 In graphene, the Debye temperature (corresponding to the theoretical highest phonon frequency in the material) is very high, $T\approx 2300$ K and $T\approx 1300$ K for planar and out-of-plane phonons, respectively Tewary and Yang (2009). However, for low electron densities, the Fermi energy can be substantially smaller than the Debye energy, and only phonons with energy smaller than two times the Fermi energy, corresponding to a full backscattering of electrons, can scatter with the electrons. This defines the Bloch-Grun̈eisen temperature $T_{BG}$ Bloch (1930), which for doped suspended graphene is around $100$ K Efetov and Kim (2010); Hwang and Das Sarma (2008). The linear relation also holds when external strain is applied to the graphene samples Castro et al. (2010).

Recently it has been shown that electronic flow in graphene can be modelled using relativistic hydrodynamic equations Müller et al. (2009); Principi et al. (2016); Levitov and Falkovich (2016); Mendoza et al. (2013, 2011); Furtmaier et al. (2015); Bistritzer and MacDonald (2009); Svintsov et al. (2012). Usually, in this hydrodynamic formalism, electron-phonon interactions are included into the equation for the electrical current density $\vec{J}$ as a damping term, $-\vec{J}/\tau_{p}$, with a characteristic relaxation time $\tau_{p}$ Bistritzer and MacDonald (2009); Svintsov et al. (2012). However, a proper inclusion of lattice-fluid interactions are missing. Here, we present a novel approach where we account for the electron-phonon interactions by including inertial corrections due to the deformations of the graphene sheet. We address the question whether these inertial corrections can recover the linear temperature relation of the electrical resistivity in the high temperature regime.

For this purpose, we use molecular dynamics simulations to simulate different graphene membranes at different temperatures and imposed strains. From the position of the atoms, we build the coordinate system in which the electrons flow, and extract the inertial corrections to include them into the two-dimensional hydrodynamic equations. In Fig.1, we can observe how the presence of thermal fluctuations produces curved streamlines of the electron velocity field.

To simulate the graphene sheets with molecular dynamics we use the adaptive intermolecular reactive bond-order (AIREBO) potential Stuart et al. (2000). This many-body potential has been developed to simulate molecules made of carbon and hydrogen atoms, and has been very successful to reproduce the phonon dispersion curves of graphene Koukaras et al. (2015). In our simulations we consider suspended graphene sheets of $L_{x}=99.2$ Å, (see Fig. 2), but also two stretched cases with $100.6$ Å and $102.1$ Å. As shown in Fig. 2 the zigzag edge is located along the $x$-direction (top and bottom) and accordingly the armchair edge along the $y$-direction (left and right). The atoms can freely move in all directions, except the ones at the left and right boundary, representing the electrical contacts (grey in Fig. 2). The distance between contacts will be denoted by $L_{x}$. Note that since the left and right boundaries are fixed, we are effectively applying strain to the graphene samples. We simulate several samples at different temperatures within the range of $100$-$600$ K controlling the temperature of the system with a Nosé-Hoover thermostat. All molecular dynamics simulations run (during the thermalisation procedure) with a time step $\Delta t=1$ fs, which is small enough to capture the dynamics of the carbon-carbon interactions accurately. Randomised velocities are attributed to the atoms at the beginning of each simulation. After reaching thermal equilibrium we are ready to couple the graphene sheet to the fluid solver and introduce electrical currents.

The electronic flow in flat graphene can be modelled using relativistic hydrodynamic equations, which are responsible for the conservation of particles, $\partial_{\mu}N^{\mu}=0$, and energy and momentum, $\partial_{\mu}T_{0}^{\mu\nu}=0$, where $N^{\mu}=nU^{\mu}$ and $T_{0}$ are the $3$-particle flow and the energy-momentum tensor, respectively. Here, $n$ is the number of particle density and $U^{\mu}$ the $3$-velocity of the electronic flow. The energy momentum tensor reads $T_{0}=(\epsilon+p)U^{\mu}U^{\nu}/v_{F}^{2}-p\eta^{\mu\nu}+\pi^{\mu\nu}$, where $\eta_{\mu\nu}$ are the components of the Minkowski metric and $\pi^{\mu\nu}$ is the shear-stress tensor, which can be approximated by the equation $\pi^{\mu\nu}\approx\kappa(\eta^{\mu\lambda}\partial_{\lambda}U^{\nu}+\eta^{\nu\lambda}\partial_{\lambda}U^{\mu})$, with $\kappa$ being the shear viscosity. $\epsilon$ and $p$ are energy density and pressure, respectively. The equation of state completes this set of equations, which for the case of graphene is given by $\epsilon=2P$. Here, the Fermi velocity is denoted by $v_{F}$. To take into account the curvature of suspended graphene sheets, we consider inertial corrections in these equations. For this purpose, we first use the covariant formulation of the hydrodynamic equations in curved space to determine the terms responsible for the inertial corrections, and afterwards, we introduce them as a forcing term into the respective equations for flat graphene. The conservation equations in curved manifolds are derived by replacing the partial derivative by the covariant derivative and extending the energy-momentum tensor to curved space, i.e. $T^{\mu\nu}=T_{0}^{\mu\nu}+p(\eta^{\mu\nu}-g^{\mu\nu})$. The covariant derivative $\nabla_{\mu}$ in the curvilinear coordinate system can be expressed through the respective partial derivative by $\nabla_{\mu}N^{\mu}=\partial_{\mu}N^{\mu}+\Gamma^{\mu}_{\mu\lambda}N^{\lambda}$, and $\nabla_{\mu}T^{\mu\nu}=\partial_{\mu}T^{\mu\nu}+\Gamma^{\mu}_{\mu\lambda}T^{\lambda\nu}+\Gamma^{\nu}_{\mu\lambda}T^{\mu\lambda}$, with $\Gamma^{\mu}_{\mu\lambda}$ being the Christoffel symbols computed as follows:

[TABLE]

At this stage we include a temperature independent coupling constant $\alpha$ in front of the inertial corrections which accounts for the strength of the electron-phonon interaction. Thus, the set of equations can be written as

[TABLE]

with

[TABLE]

To solve the hydrodynamic equations, we use the lattice Boltzmann solver described in Ref. Oettinger et al. (2013); Furtmaier et al. (2015), and add the inertial contributions as external forces (see Supplementary Information). The model discretises space as regular triangular lattice, which we couple to graphene’s hexagonal lattice. As shown in Fig. 2, the atoms form a two-dimensional manifold which can be described by a discrete mapping $\Phi(ct,x,y)$ from the curved space to the three-dimensional flat space (reference frame of the laboratory, where the metric is given by the Minkowski-metric). The metric tensor can be computed by

[TABLE]

Similarly, the Christoffel symbols are obtained from Eq. (1). Further details on this calculations can be found in the Supplementary Information. The graphene sheet possesses zigzag boundaries at the left and right end. For the fluid solver we impose periodic boundary conditions at these boundaries, which correspond to the in- and outlet. The free sides of the graphene sheet possesses armchair geometry and, therefore, we impose free slip boundary conditions.

In order to obtain general results and allow for more insight, we make the hydrodynamic equations dimensionless (See Supplementary Information for the detailed analysis). Considering that the shear-stress tensor depends linearly on the shear viscosity (Cercignani and Kremer, 2002, p. 109) we find that the physics of the electronic flow only depends on two dimensionless numbers: $A_{1}=\kappa_{0}v_{F}^{\prime 2}u_{0}/x_{0}\epsilon_{0}$, and $A_{2}=e_{0}E_{0}/x_{0}\epsilon_{0}v_{F}^{\prime 2}$, where $\epsilon_{0}$, $u_{0}$, $x_{0}$, $\kappa_{0}$, $e_{0}$, and $E_{0}$ are characteristic values for the energy density, velocity, length, shear viscosity, electric charge, and electric field, respectively. Note that $1/A_{1}$ is related to the Reynolds number Pozrikidis (2009). The parameter range of our simulations is $0.14\leq A_{1}\leq 2.05$ and $10^{-8}\leq A_{2}\leq 10^{-6}$.

To produce an electrical current, we apply an external electric field in $x$ direction, $E^{x}=1.43\times 10^{9}-10^{11}$ V/m. In our simulations, we set $n=1.36\times 10^{13}$ /cm2, chemical potential $\mu=\hbar v_{F}\sqrt{\pi n}\approx 6.89\times 10^{-20}$ kgm2/s${}^{2}\approx 0.43$ eV, and $p=\frac{1}{3}\mu n\approx 3.12\times 10^{-3}$ kg/s2. We use the shear viscosity for doped graphene given by Principi et al. (2016)

[TABLE]

with $c_{\nu}\approx 1.33\times 10^{-34}\text{kgm}^{2}\text{/s}\approx\frac{5}{4}\hbar$. We couple the electronic fluid to the atomistic simulation using the same length-scale $x_{0}=\Delta x_{MD}=1$ Å  and time scale $t_{0}=\Delta t_{MD}=10^{-5}$ fs. After each iteration we compute the metric tensor and the Christoffel symbols and simulate the electronic flow until we obtain the electrical current in steady state. The momentum change in the fluid is imposed on the atoms to ensure momentum conservation of the system.

We perform simulations for different strains and different values of the coupling constant and measure the resistivity of the graphene sheets using Ohm’s law, $\langle\rho\rangle=E^{x}L_{y}/I(t\to\infty)$, where the electrical current is given by

[TABLE]

In Fig. 3, we observe that all our calculations of the electrical resistivity of graphene exhibit a linear dependence with temperature, $\rho=\gamma_{s}(\alpha)T+\rho_{0}(\alpha,L_{x})$. This linear dependence is in agreement with experimental measurements and theoretical predictions Efetov and Kim (2010); Hwang and Das Sarma (2008). Additionally, the slope $\gamma_{s}(\alpha)$ only depends on the strength of the coupling constant $\alpha$ but not on the applied strain, which is also in agreement with previous works Castro et al. (2010). We can also observe from the inset of Fig. 3 that $\gamma_{s}\propto\alpha^{2.15\pm 0.03}$. The factor $\gamma_{s}$ can be compared with the values from the Boltzmann transport theory Hwang and Das Sarma (2008), where

[TABLE]

with $\rho_{m}=7.6\times 10^{-7}$ kg/m2 being the mass density of graphene, $v_{ph}$ the phonon velocity and $D$ the deformation-potential coupling constant. Experimental and theoretical results suggest that the deformation-potential varies within the range of $4.5\leq D\leq 50$ eV Borysenko et al. (2010); Chen et al. (2008); Bolotin et al. (2008), and phonon velocities within $2\times 10^{4}\leq v_{ph}\leq 3\times 10^{4}$ m/s. Thus, one expects $\gamma$ to be in the range $0.012\leq\gamma\leq 2.8~{}\Omega$/K. From our results we deduce that $0.01\leq\alpha\leq 0.14$ to obtain results compatible with $\gamma$. The residual resistivity $\rho_{0}$ is $704\pm 33~{}\Omega$ for $\alpha=0.1$ and is a consequence of the static corrugations (ripples) present at any temperature. These ripples stabilize the two-dimensional crystal circumventing the Mermin Wagner theorem, which states that crystalline order cannot exist in two dimensions Mermin (1968).

The existence of a resistivity points to the presence of dissipation. Recently, it has been shown that energy dissipation in flows through curved manifolds arises from the curvature Debus et al. (2015). The Ricci scalar (or curvature scalar) is a measure for the curvature and is given by $R^{\mu}_{\mu}=g^{\nu\sigma}(\Gamma^{\mu}_{\nu\sigma,\mu}-\Gamma^{\mu}_{\nu\mu,\sigma}+\Gamma^{\delta}_{\nu\sigma}\Gamma^{\mu}_{\mu\delta}+\Gamma^{\delta}_{\nu\mu}\Gamma^{\mu}_{\sigma\delta})$. We have analysed the Ricci scalar for different strains and temperatures and found that its maximal absolute value remains constant during time evolution and depends only on the temperature of the system (see inset of Fig. 4). As function of temperature the Ricci scalar follows a power-law with an exponent that increases with the applied strain (see main panel in Fig. 4).

We have also studied the standard deviation of the heights of our graphene sheets, and found that they are virtually not influenced by temperature, within the range of $100-600$ K (see Fig. 5 of Supplementary Information). Therefore, we conclude that, an increase of temperature mainly induces in-plane atomistic motion, increasing the average effective bond length (see Figs. 6 a-c in the Supplementary Information), and consequently, increasing the local curvature. Since curvature induces energy dissipation Debus et al. (2015, 2016), this explains why larger temperatures result in larger electrical resistivities.

To summarise, we have shown that inertial corrections lead to a linear dependence of the electrical resistivity with temperature in suspended graphene, which is in agreement with experimental measurements of the electrical resistivity due to electron-phonon interactions at high temperatures. To characterize the strength of the electron-phonon interaction in our model, we have introduced a coupling constant $\alpha$, and for values of $\alpha\approx 0.1$ one recovers the same values of $\gamma$ known from theoretical predictions and experimental measurements. Additionally, the resistivity does not depend on the strain imposed to the sample, which is also in agreement with previous work Castro et al. (2010). We also analysed the maximal absolute value of the curvature as function of temperature and found a power-law behaviour with an exponent that changes slightly with the applied strain. Finally, by studying the average height fluctuations of the graphene sheets, we also discovered that by increasing the temperature, the height fluctuations remain almost unchanged, and consequently, including larger local values of curvature. In our model, larger curvature implies more energy dissipation, and consequently, larger electrical resistivity.

Our finding opens up many interesting questions, as for instance, if one can discern the influence of different types of vibrational phonons, e.g. flexural, acoustic, and optical, by studying how they introduce curvature into the hydrodynamic system. By extending our molecular dynamics simulations, one could also explore the electrical resistivity due to the electron-phonon interactions when the graphene sample is placed on a substrate Chen et al. (2008).

Finally, our approach can also be applied to other two- and three-dimensional crystals, where electrons are well described as Fermi liquid. In those structures, ion displacements due to phonons can induce intrinsic curvature, and consequently, energy dissipation in the electronic flow. Applying our model to other materials will be an interesting subject for future research.

Appendix A Conservation equations for curved spaces

The lattice Boltzmann model for flat space reported in Ref. Oettinger et al. (2013) needs to be extended in order to simulate relativistic hydrodynamics in curved spaces. This numerical model solves the following conservation equations:

[TABLE]

where $T_{0}^{\mu\nu}$ is the energy momentum tensor for flat space,

[TABLE]

In curved spaces, the energy momentum tensor has to be conserved as well. However, conservation of quantities in a curved manifold needs to be expressed through covariant derivatives, instead of partial derivatives. The conservation equations for the number of particles and energy momentum tensor are given by

[TABLE]

where the covariant derivatives $\nabla_{\mu}$ of the 3-particle flow and energy-momentum tensor relate to the partial derivatives $\partial_{\mu}$ as follows:

[TABLE]

where $\Gamma^{\mu}_{\mu\lambda}$ denote the Christoffel symbols which are computed with the metric tensor using

[TABLE]

Additionally, the energy momentum tensor in curved manifolds reads

[TABLE]

where $g^{\mu\nu}$ is the metric tensor and $\pi^{\mu\nu}$ the shear-stress tensor, which can be approximated by the equation $\pi^{\mu\nu}\approx\kappa(g^{\mu\lambda}\nabla_{\lambda}U^{\nu}+g^{\nu\lambda}\nabla_{\lambda}U^{\mu})$. Using the property that $\nabla_{\mu}g^{\mu\nu}=0$, we can rewrite the covariant derivative of the energy-momentum tensor as

[TABLE]

Thus the conservation equations (10) expressed with the energy-momentum tensor for the flat space $T_{0}^{\mu\nu}$ reads

[TABLE]

where we have used the fact that partial derivatives equal covariant derivatives when acting on scalars. At this stage we also include a temperature independent coupling constant $\alpha$ which accounts for the strength of the inertial corrections. Thus, this set of equations can be written as

[TABLE]

with

[TABLE]

Here $F_{T}^{\nu}$ and $F_{N}$ can be introduced in the numerical model as external forces using the technique described in Ref. Furtmaier et al. (2015).

Appendix B Dimensionless quantities

In order to obtain general results we transform the hydrodynamic equations to dimensionless form and determine the dimensionless quantities that characterise our systems. We also include a force density $\vec{F}^{ext}=ne\vec{E}$ describing an external electric field $\vec{E}$ to the energy-momentum conservation equation (16).

To write the equations in a dimensionless form, we first express all relevant quantities in dimensionless form: $x=x_{0}x^{\prime},u=u_{0}u^{\prime},t=t_{0}t^{\prime},\epsilon=\epsilon_{0}\epsilon^{\prime},E=E_{0}E^{\prime}$, $\kappa=\kappa_{0}\kappa^{\prime}$ $e=e^{\prime}e_{0}$, $n=n^{\prime}/x_{0}^{2}$ and $v_{F}=u_{0}v_{F}^{\prime}$ where all primed variables are dimensionless. Furthermore, the relation $u_{0}t_{0}=x_{0}$ holds, which allows us to write the temporal derivative $\partial_{0}=(1/v_{F})\partial_{t}=(1/v_{F}^{\prime}u_{0}t_{0})\partial_{t^{\prime}}=(1/v_{F}^{\prime}x_{0})\partial_{0}^{\prime}$, and thus for $v_{F}^{\prime}=1$, the relation $\partial_{\mu}=(1/x_{0})\partial_{\mu}^{\prime}$ is also satisfied.

We divide the energy-momentum tensor in an equilibrium part and a dissipative part $T^{\mu\nu}=T^{\mu\nu(eq)}+\pi^{\mu\nu}$. The equilibrium part can be written in its dimensionless form $T^{\mu\nu(eq)}=\epsilon_{0}T^{\prime\mu\nu(eq)}$. Assuming that the shear-stress tensor is linearly dependent on the shear viscosity $\kappa$ (Cercignani and Kremer, 2002, p. 109) and taking into account that the shear stress tensor has the same units as the energy-momentum tensor ($\kappa_{0}=\frac{\epsilon_{0}x_{0}}{u_{0}}$), we get $\pi^{\mu\nu}=\kappa_{0}(\epsilon_{0}/\kappa_{0})\pi^{\prime\mu\nu}=\kappa_{0}(u_{0}/x_{0})\pi^{\prime\mu\nu}$. The force densities (17) transform as $F_{N}=(\epsilon_{0}/x_{0})F_{N}^{\prime}$ and $F_{T}^{\nu}=(\epsilon_{0}/x_{0})F_{T}^{\prime\nu}$, respectively. This can be derived from the fact that, in our case, the metric is dimensionless and thus the Christoffel symbols (as derivative of the metric) transform as $\Gamma^{\mu\nu}_{\sigma}=(1/x_{0})\Gamma^{\prime\mu\nu}_{\sigma}$.

With the considerations above, equations (16) read

[TABLE]

While the first equation remains unchanged, we observe after multiplying by $(x_{0}/\epsilon_{0})$ that the second equation only depends on two dimensionless numbers

[TABLE]

and

[TABLE]

Appendix C Strain

All simulations are performed with the same number of carbon atoms ($N=4114$). The atoms at the left and right boundary are fixed at a distance $L_{x}$ from each other and are not allowed to move nor change their effective bond lengths. The simulations are preformed for different $L_{x}=99.2-104.2$ Å  which correspond to initial effective bond lengths $a=1.38-1.45$ Å.

The distance $L_{x}$ is related to the tension imposed on the graphene membrane. One expects that increasing $L_{x}$ results in a more stretched membrane and the height fluctuations are reduced. Figure 5 confirms that the height fluctuations $\sqrt{\langle z^{2}\rangle-\langle z\rangle^{2}}$ depend on $L_{x}$ and therefore also on the strain imposed on the graphene membrane. Interestingly, the height fluctuations are independent on the temperature. In contrast, the width of the effective bond length distribution presents a strong temperature dependence as shown in Figs. 6 a-c.

This finding is related to the Ricci scalar analysis (see Fig. 4), where it has been shown that the absolute value of the Ricci scalar has a much stronger temperature dependence than a dependence on $L_{x}$. A higher temperature results in a broader width of the effective bond length distribution. This leads to a higher average Ricci scalar and consequently to more curvature in the system. The curvature introduces shear which is responsible for saturation of the current.