# Fluctuation-induced currents in suspended graphene nanoribbons: Adiabatic
quantum pumping approach

Federico D. Ribetto  Instituto de Física Enrique Gaviola (CONICET) and FaMAF,
Universidad Nacional de Córdoba, Argentina  Departamento de Física,
Universidad Nacional de Río Cuarto, Río Cuarto, Argentina  Silvina A. Elaskar
Facultad de Ciencias Químicas, Universidad Nacional de Córdoba, Argentina
Hernán L. Calvo  Instituto de Física Enrique Gaviola (CONICET) and FaMAF,
Universidad Nacional de Córdoba, Argentina  Raúl A. Bustos-Marún
Corresponding author<EMAIL_ADDRESS>Instituto de Física Enrique
Gaviola (CONICET) and FaMAF, Universidad Nacional de Córdoba, Argentina
Facultad de Ciencias Químicas, Universidad Nacional de Córdoba, Argentina

###### Abstract

Graphene nanoribbons (GNRs) are thin strips of graphene with unique properties
due to their structure and nanometric dimensions. They stand out as basic
components for the construction of different types of nanoelectromechanical
systems (NEMS), including some very promising sensors and pumps. However,
various phenomena, such as unintended mechanical vibrations, can induce
undesired electrical currents in these devices. Here, we take a quantum
mechanical approach to analyze how currents induced by fluctuations (either
thermal or of some other kind) in suspended GNRs contribute to the electric
current. In particular, we study the pumping current induced by the adiabatic
variation of the Hamiltonian of the system when a transverse vibration
(flexural mode) of a GNR suspended over a gate is excited. Our theoretical
approach and results provide useful tools and rules of thumb to understand and
control the charge current induced by fluctuations in GNR-based NEMS, which is
important for their applications in nanoscale sensors, pumps, and energy
harvesting devices.

##  I Introduction

Figure 1: (a) Scheme of the studied system: A suspended GNR between two
metallic contacts freely oscillating, at a distance $d_{0}$, over a gate with
gate voltage $V_{g}$. We only consider the lowest frequency (acoustic) modes
for the transversal vibration of the GNRs in a classical approximation and
assume a rectangular membrane model to describe them. The shown mode
corresponds to $(n_{x},n_{y})=(1,0)$. (b) Other examples of normal modes with
(from left to right): $(n_{x},n_{y})=(1,1)$, $(2,0)$, and $(2,1)$. (c) and (d)
Logarithm of the scaled emissivity $\tilde{I}_{k}$ in the low temperature
limit as a function of $n_{x}$ and $n_{y}$ for square nanoribbons with
armchair (aGNR) and zigzag (zGNR) edges. $\tilde{I}_{k}$ has units of 1/eV and
is proportional to the maximum value of the pumping current during a period.
The size of the ribbons is $L=449\,a$ (with $a=0.246$ nm) and the Fermi energy
is $\varepsilon_{\text{F}}=0.1$ eV.

Being one of the strongest materials ever tested Lee  _et al._ (2008) with
large thermal and electrical conductivities, Chen  _et al._ (2008, 2011)
graphene and, in particular, the more tunable graphene nanoribbons (GNRs) have
become the conductive building blocks of innumerable nanoelectromechanical
systems (NEMS). GNR-based NEMS have been studied as electromechanical
resonators; Bunch  _et al._ (2007); Garcia-Sanchez  _et al._ (2008); Chen  _et
al._ (2009); Croy  _et al._ (2012) electron pumps; Zhu and Chen (2009); Prada
_et al._ (2009); Foa Torres  _et al._ (2011); San-Jose  _et al._ (2011); Low
_et al._ (2012); Ingaramo and Foa Torres (2013); Connolly  _et al._ (2013)
sensors of mass, pressure, strain, and temperature, Zang  _et al._ (2015);
Davaji  _et al._ (2017); Khan  _et al._ (2017) as well as detectors of
vibrations Moreno-Garcia  _et al._ (2022) and gases Gupta  _et al._ (2021);
switches; Liao and Koide (2011) ultrasmall accelerometers; Fan  _et al._
(2019) and even viral detectors for COVID-19. Ben-Shimon  _et al._ (2022)

In all of the above examples, it is essential to study the various phenomena
that can limit the use of the proposed devices, such as the noise in
electrical currents. Besides the thermal (or Nyquist-Johnson) noise, Johnson
(1928); Nyquist (1928) other phenomena can also interfere with the detected
electrical currents. For example, suspended GNRs are in constant mechanical
motion at room temperature, which can induce measurable oscillations in the
electrical current, as recently shown. Thibado  _et al._ (2020) Of course,
oscillations are not only due to thermal excitations and other phenomena, such
as the propagation of vibrational waves through the material, can also cause
them.

Typically, the calculations of instantaneous electric currents induced by
vibrations are based on classical models of time-dependent capacitances.
Thibado  _et al._ (2020) Here, instead, we adopt a quantum mechanical approach
to this problem and analyze the contributions to the current due to the
adiabatic quantum pumping Brouwer (1998); Bode  _et al._ (2011); Bustos-Marún
(2018) arising from the stochastic oscillation of GNRs. Quantum pumping
currents are a consequence of the delayed response of electronic wave
functions to a time-dependent Hamiltonian which, in the present case,
originates from the movement of the nuclei. In this regard, a comment is in
order to avoid confusions. Since the oscillations of individual vibrational
modes are independent, the generation of pumping currents produced by them is
expected to average to zero over a period. This is because, to have a finite
pumped charge per cycle, more than one time-dependent parameter needs to be
moved with a fixed phase. Brouwer (1998) However, even if the pumping currents
have a null mean value, this does not prevent them from contributing
instantaneously to the total current and, thus, to the current noise, when the
Hamiltonian parameters move stochastically. This crucial point marks a clear
difference with some previous works, e.g., Refs. [Low  _et al._ , 2012; Foa
Torres  _et al._ , 2011; Ingaramo and Foa Torres, 2013; Connolly  _et al._ ,
2013]. While there the authors study the pumped charge in GNRs (meaning the
average DC current per cycle), here we are interested in the maximum value of
the instantaneous adiabatic pumping currents induced by stochastic
fluctuations of GNRs.

In particular, we study a common configuration of GNR-based NEMS consisting of
a GNR suspended over a controllable gate, see Fig. 1(a). There, we evaluate
the contribution to the electric current due to quantum pumping induced by the
movement of the lowest frequency (acoustic) modes for the transversal
vibration of suspended GNRs, also known as flexural modes, Low  _et al._
(2012) see Figs. 1(b) and (c). To this end, we adapted the theoretical
description of adiabatic quantum pumping Brouwer (1998) to the generic case of
vibrational normal modes. Bustos-Marún (2018) Our methodology and results are
not only relevant for different types of sensors and pumps based on GNRs, but
they can also be extended to other applications such as energy harvesting at
the nanoscale Kim  _et al._ (2017); Bustos-Marún (2018); Thibado  _et al._
(2020) or to other systems such as carbon nanotubes.

The paper is organized as follows. In Sec. II, we present the theoretical
framework, including the different models and approximations used in our
calculations. In Sec. III, we explore the role of different system’s
parameters on the fluctuation-induced currents, including the length, type of
border, and Fermi energy of GNRs. At the end of this section we also present a
realistic estimation of the expected value of currents. In Sec. IV we develop
a semiclassical theory to evaluate the zero-frequency noise of thermally
induced pumping current and use it to compare with Nyquist-Johnson zero-
frequency noise. Finally, in Sec. V we summarize and discuss the main results.

##  II Theoretical framework

Due to the vast number of variables that can potentially influence the current
induced by the oscillations of GNRs, throughout the manuscript we focus on
certain limits, approximations, and simplified models that allow us to obtain
simple expressions that nonetheless can be used to understand the general
features of the studied phenomenon.

###  II.1 GNR: electronic modelling

Hamiltonian. Considering energies close to the Fermi energy, the description
of the electrons in graphene can be carried out by means of a tight-binding
model. Foa Torres  _et al._ (2020) Given the energy difference between the
$\sigma$ and $\pi$ molecular orbitals, and the fact that the band of $\pi$
orbitals is the one found around the Fermi energy, we will only consider
$p_{z}$ atomic orbitals. Thus, we have the following Hamiltonian:

$\hat{H}=\sum_{i}\varepsilon_{i}\hat{c}_{i}^{\dagger}\hat{c}_{i}-\sum_{\braket{i,j}}t_{i,j}\hat{c}_{i}^{\dagger}\hat{c}_{j},$
(1)

where $\varepsilon_{i}$ is the energy of site $i$, $\hat{c}_{i}^{\dagger}$ and
$\hat{c}_{i}$ are the creation and annihilation operators in the $p_{z}$
orbital of site $i$, and $t_{i,j}$ represents the hopping amplitude between
sites $i$ and $j$. Furthermore, the sum $\braket{i,j}$ is restricted to
nearest neighbor sites. In the absence of defects and external disturbances,
$\varepsilon_{i}$ is set to zero for simplicity, while the bare hopping
amplitude takes the value $t_{i,j}\equiv t_{0}=2.66$ eV Foa Torres  _et al._
(2020).

Electronic properties. Graphene nanoribbons are usually cuts in a certain
direction. Based on the direction of the cut, typically two edge types are
described: zigzag edge [Fig. 2(a)] and armchair edge [Fig. 2(b)]. To classify
the ribbons, the following convention will be used: GNRs with armchair
(zigzag) edges are classified by the number of dimer lines (zigzag lines)
across the width of the ribbon. In addition, the notation $N$-aGNR ($N$-zGNR)
will be used for armchair (zigzag) GNRs, where $N$ is the number of dimer
lines (zigzag lines). Foa Torres  _et al._ (2020)

$N$-aGNR and $N$-zGNR have very different electronic properties that arise
from their contrasting boundary conditions. Some aGNR exhibit semiconductor
behavior, while others are metallic. An analytical calculation of the
eigenvalues of the tight-binding Hamiltonian allows us to show that the energy
gap $\Delta_{N}$ of $N$-aGNR oscillates with the width of the ribbon, Cresti
_et al._ (2008) besides the obvious decaying limit $\Delta_{N}\rightarrow 0$
for $N\rightarrow\infty$. In particular, we have that $\Delta_{N}=0$ for
$N=3\ell+2$, where $\ell$ is an integer, making it metallic or semiconducting
otherwise. Dutta and Pati (2010); Foa Torres  _et al._ (2020)

Unlike the previous case, the zGNRs retain the semimetallic character of
graphene, regardless of their width. Another interesting feature is the
formation of a pronounced peak in the density of states for $\varepsilon=0$,
which results from the formation of partially flat and degenerate bands with
zero energy. Foa Torres  _et al._ (2020) The presence of these highly confined
electronic edge states has been confirmed by scanning tunneling microscopy and
spectroscopy. Foa Torres  _et al._ (2020) They have a “topological protection”
that makes them robust against different types of disturbances such as
vacancies or Anderson-type noise Anderson (1978); Foa Torres  _et al._ (2020).

System and contacts. In our calculations we will assume the configurations
shown in Fig. 2 of the system depicted in Fig. 1(a). There, we define a finite
central region (blue and orange dots in Fig. 2), that is able to oscillate and
is coupled to two semi-infinite and identical contacts, left ($L$) and right
($R$) (red dots), that are also made of GNRs of the same type than that of the
central region.

Figure 2: Examples of GNRs with zigzag (zGNR) and armchair (aGNR) edges (blue
and orange sites) connected to two semi-infinite contacts on the left and
right (red sites). Note that the contacts are also made of GNRs of the same
width and type. The green rectangles show how a row is defined for each type
of GNR (to be used in Sec. III.2).

###  II.2 GNR: Oscillating membrane model

We are going to assume that the GNRs are large enough to be able to
approximate its transversal normal modes of vibration by those of a
rectangular membrane (limit of the continuum for low-frequency acoustic
modes). We also assume that this membrane is placed in the $xy$ plane, and
$z(x,y,t)$ is a function that describes its displacement with respect to said
plane, see Fig. 1(a). Consequently, the equation that governs the transverse
movement of the membrane is (see Appendix A for more details)

$z(x,y,t)=\sum_{n_{x},n_{y}=0}^{\infty}\sin\left(\frac{n_{x}\pi}{L_{x}}x\right)\cos\left(\frac{n_{y}\pi}{L_{y}}y\right)q_{n_{x},n_{y}},$
(2)

where $L_{x}$ and $L_{y}$ are the length and width of the central region,
respectively, and the amplitude of the mode $q_{n_{x},n_{y}}$ (defined by the
pair of integers $n_{x}=1,2,...$ and $n_{y}=0,1,...$) is

$q_{n_{x},n_{y}}\equiv
A_{n_{x},n_{y}}\cos(\omega_{n_{x},n_{y}}t+\phi_{n_{x},n_{y}}).$ (3)

Fig. 1(b) shows some examples of modes for different pairs $(n_{x},n_{y})$. In
the above equation, the phase $\phi_{n_{x},n_{y}}$ of the oscillation is
arbitrary, the maximum amplitude of $q$ will be discussed afterwards, and the
frequencies of the normal modes are given by

$\omega_{n_{x},n_{y}}=v\sqrt{\left(\frac{n_{x}\pi}{L_{x}}\right)^{2}+\left(\frac{n_{y}\pi}{L_{y}}\right)^{2}},$
(4)

where $v$ is the speed of sound in graphene (for low frequency transversal
acoustic waves). For simplicity, from now on we will collapse the pair of
indices $n_{x}$ and $n_{y}$ into a single one $k$. A central quantity that
will be needed in our calculations is the maximum velocity of the normal
modes, given by

$\dot{q}_{k,\text{max}}=A_{k}\omega_{k}.$ (5)

The energy associated with the normal modes of graphene, see Appendix B, can
be written as

$E_{k}=m_{\text{c}}\left(\frac{1}{2}\dot{q}_{k}^{2}+\frac{1}{2}\omega_{k}^{2}q_{k}^{2}\right)=\frac{1}{2}m_{\text{c}}A_{k}^{2}\omega_{k}^{2},$
(6)

where $m_{\text{c}}$ is the mass of a carbon atom. Now we can use
equipartition theorem together with the above equation to derive an estimation
for $A_{k}$, which yields

$A_{k}=\frac{1}{\omega_{k}}\sqrt{\frac{2k_{\text{B}}T}{m_{\text{c}}}}.$ (7)

As we will see in Eq. (11), the above expression implies that the maximum
value of the pumping current is independent of the frequency of the mode.

###  II.3 Quantum pumping current

The objective of this section is to derive an expression that allows us to
calculate the maximum value of the pumping current on the reservoir $r$,
$I_{r,k,\text{max}}^{(1)}\equiv\max[I_{r,k}^{(1)}(t)]$, over a period of
oscillation of a given normal mode $k$ of a suspended GNR. This pumping
current, typically denoted with a superscript (1), is associated with the
first order term of an adiabatic expansion of the observable current, and
arises from the delayed electronic response of the system to the time-
variation of some classical parameters, see Appendix C. It has been thoroughly
studied by using different formalisms Brouwer (1998); Splettstoesser  _et al._
(2005); Arrachea and Moskalets (2006) and in a variety contexts (even outside
quantum transport Li  _et al._ (2019)) but only as a way of generating DC
currents controlled by the movement of at least two system’s parameters.

Since our interest here lies in stochastic variations of the system’s
parameters (the vibrational normal modes $q_{k}$), it is necessary to
highlight that the associated pumping currents average to zero in one period
of $q_{k}(t)$. However, this does not mean that $I_{r,k}^{(1)}(t)$ is always
zero. On the contrary, vibration-induced currents can instantaneously (at a
given time) contribute to the total current which, as discussed in the
introduction, can potentially interfere with current measurements in sensors
based on GNRs, for example.

We will work under a perturbative limit of the modes, i.e., small amplitudes
of the oscillation, and thus $\bm{q}(t)\approx\bm{q}_{0}$, where $\bm{q}$ is
the vector containing the amplitudes $q_{k}$ of all normal modes, and
$\bm{q}_{0}\equiv\bm{0}$ is the equilibrium position of the system. 111In our
model, the equilibrium position (no excitation of vibrational modes for the
oscillation of GNRs) corresponds to a flat still membrane, equivalent to
setting $\bm{q}_{0}=\bm{0}$ Therefore, we can approximate $I_{r}^{(1)}(t)$ as

$I_{r}^{(1)}(t)=e\sum_{k}\frac{\text{d}n_{r}(\bm{q})}{\text{d}q_{k}}\dot{q}_{k}(t)\approx
e\sum_{k}\left(\frac{\text{d}n_{r}}{\text{d}q_{k}}\right)_{\bm{q}_{0}}\dot{q}_{k}(t),$
(8)

where $(\text{d}n_{r}/\text{d}q_{k})_{\bm{q}_{0}}$ is the emissivity due to
the $k$-th normal mode evaluated at $\bm{q}_{0}$. This emissivity can be
written in terms of the scattering matrix $S$ of the system as Brouwer (1998)

$\left(\frac{\text{d}n_{r}}{\text{d}q_{k}}\right)_{\bm{q}_{0}}=\int\frac{\text{d}\epsilon}{2\pi}\left(-\frac{\partial
f}{\partial\epsilon}\right)\sum_{\alpha\in
r,\beta}\text{Im}\left[\frac{\partial S_{\alpha\beta}}{\partial
q_{k}}S_{\alpha\beta}^{*}\right]_{\bm{q}_{0}},$ (9)

where $f$ is the equilibrium Fermi function (we are assuming a zero bias
voltage configuration), $\alpha$ and $\beta$ are channel indices of the
contacts, and we require $\alpha\in r$, since we want to calculate the
emissivity at contact $r$ to determine the associated quantum pumping current
$I_{r}^{(1)}$. If we additionally take the limit of low temperatures, Eq. (9)
reduces to

$\left(\frac{\text{d}n_{r}}{\text{d}q_{k}}\right)_{\varepsilon_{\text{F}},\bm{q}_{0}}=\sum_{\alpha\in
r,\beta}\frac{1}{2\pi}\text{Im}\left[\frac{\partial S_{\alpha\beta}}{\partial
q_{k}}S_{\alpha\beta}^{*}\right]_{\varepsilon_{\text{F}},\bm{q}_{0}}.$ (10)

It is interesting to note that the above equation holds even at finite
temperatures if the right-hand side of Eq. (9) depends linearly on the energy
around $\varepsilon_{\text{F}}$. Due to the approximation used, emissivities
do not depend on time. Then, the maximum value of the $k$-th contribution to
the pumping current $I_{r,k}^{(1)}$ (defined through
$I_{r}^{(1)}(t)=\sum_{k}I_{r,k}^{(1)}$) is simply

$I_{r,k,\text{max}}^{(1)}=\left|e\left(\frac{\text{d}n_{r}}{\text{d}q_{k}}\right)_{\varepsilon_{\text{F}},\bm{q}_{0}}\dot{q}_{k,\text{max}}\right|,$
(11)

where $\dot{q}_{k,\text{max}}$ is the maximum value of $\dot{q}_{k}$.

###  II.4 Emissivities

The explicit expression of the emissivity depends on the Hamiltonian model
used. Here, we assume small displacements in the direction perpendicular to
the GNR plane, $z$. Thus, we can take a linear regime of the diagonal elements
of the Hamiltonian,

$H_{ii}=E_{0}+\left(\frac{\partial E}{\partial z}\right)\delta z_{i},$ (12)

where $\delta
z_{i}=[z_{i}(x_{i},y_{i},\bm{q})-z_{i}(x_{i},y_{i},\bm{q}_{0})]$,
$E_{0}=\varepsilon_{i}=0$, and the factor $(\partial E/\partial z)$ is the
same for all sites. Taking the limit of large sizes for the GNR, within which
our model should behave like a classical parallel plate capacitor, it can be
shown (see Appendix D) that this factor should satisfy the relation

$\left(\frac{\partial E}{\partial
z}\right)=-\epsilon_{0}\frac{A_{\text{site}}V_{g}^{2}}{d_{0}^{2}},$ (13)

where $A_{\text{site}}$ is the area associated with each individual site,
$V_{g}$ is the potential difference between the membrane and the gate (gate
voltage), $d_{0}$ is the distance between the GNR and the gate, and
$\epsilon_{0}$ is the vacuum permittivity.

In GNRs, the hopping amplitudes in the Hamiltonian usually take an exponential
dependence on the distance between first neighbors Pereira  _et al._ (2009)

$t_{i,j}=t_{0}e^{-b\left(\frac{|\bm{r}_{i}-\bm{r}_{j}|}{a_{\text{cc}}}-1\right)},$
(14)

where $a_{\text{cc}}$ is the equilibrium distance between neighboring atoms
and $b$ is a constant that sets how much the hoppings change with the distance
$|\bm{r}_{i}-\bm{r}_{j}|$. Here, the distance $|\bm{r}_{i}-\bm{r}_{j}|$ is a
quantity that depends on the modes that are being excited, and can be
expressed as

$|\bm{r}_{i}-\bm{r}_{j}|=\sqrt{a_{\text{cc}}^{2}+\left(\delta z_{i}-\delta
z_{j}\right)^{2}}.$ (15)

Given that we are considering that the amplitude $A_{k}$ of the modes is
small, we can approximate

$t_{i,j}\approx t_{0}\left[1-\frac{b}{2}\left(\frac{\delta z_{i}-\delta
z_{j}}{a_{\text{cc}}}\right)^{2}\right].$ (16)

The Hamiltonian of the system can then be divided into two parts,
$\bm{H}=\bm{H}^{(E)}+\bm{H}^{(V)}$, where $\bm{H}^{(E)}$ contains
$\bm{H}\left(\bm{q_{0}}\right)$ and the site energies dependence on $\bm{q}$,

$H_{ij}^{(E)}=\begin{cases}E_{0}+\left(\dfrac{\partial E}{\partial
z}\right)\delta z_{i}&i=j\\\ -t_{0}&i=j\pm 1,\end{cases}$ (17)

and $\bm{H}^{(V)}$ contains only the hoppings dependence on $\bm{q}$

$H_{ij}^{(V)}=\begin{cases}0&i=j\\\ \dfrac{t_{0}b}{2}\left(\dfrac{\delta
z_{i}-\delta z_{j}}{a_{\text{cc}}}\right)^{2}&i=j\pm 1.\end{cases}$ (18)

Then, using the Fisher-Lee formula one can prove (see Appendix E) that the
derivatives with respect to $q_{k}$ of the scattering matrix can be decomposed
into

$\frac{\partial\bm{S}}{\partial q_{k}}=\frac{\partial\bm{S}^{(E)}}{\partial
q_{k}}+\frac{\partial\bm{S}^{(V)}}{\partial q_{k}},$ (19)

where $\bm{S}^{(E)}$ and $\bm{S}^{(V)}$ are the scattering matrices obtained
from Hamitonians $\bm{H}^{(E)}$ and $\bm{H}^{(V)}$, respectively. It can be
readily shown that, for small oscillations, $\partial\bm{S}^{(V)}/\partial
q_{k}$ can be neglected, see Appendix E. Using this, and noticing that
$\bm{S}(\bm{q}_{0})=\bm{S}^{(E)}(\bm{q}_{0})$, the emissivity takes the form

$\left(\frac{\text{d}n_{r}}{\text{d}q_{k}}\right)_{\varepsilon_{\text{F}},\bm{q}_{0}}\approx\left(\frac{\text{d}n_{r}^{(E)}}{\text{d}q_{k}}\right)_{\varepsilon_{\text{F}},\bm{q}_{0}},$
(20)

where $\text{d}n_{r}^{(E)}/\text{d}q_{k}$ is the emissivity calculated with
Hamiltonian $\bm{H}^{(E)}$.

To simplify the analysis of the next sections, we define the scaled
emissivity, which removes from $\text{d}n_{r}^{(E)}/\text{d}q_{k}$ the
parameters related to particularities of the studied system, i.e.,

$\left(\frac{\text{d}\tilde{n}_{r}^{(E)}}{\text{d}q_{k}}\right)\equiv\left(-\epsilon_{0}\frac{A_{\text{site}}V_{g}^{2}}{d_{0}^{2}}\right)^{-1}\left(\frac{\text{d}n_{r}^{(E)}}{\text{d}q_{k}}\right).$
(21)

Finally, the maximum value of the quantum pumping current, given by the
oscillation of a mode $q_{k}$ with energy $k_{\text{B}}T$, yields

$I_{r,k,\text{max}}^{(1)}=\left|\frac{e\epsilon_{0}A_{\text{site}}V_{g}^{2}}{d_{0}^{2}}\right|\sqrt{\frac{2k_{\text{B}}T}{m_{\text{c}}}}\left|\frac{\text{d}\tilde{n}_{r}^{(E)}}{\text{d}q_{k}}\right|_{\varepsilon_{\text{F}},\bm{q}_{0}}.$
(22)

This expression can be divided into three contributions: the first one from
left to right is the scaling factor and it accounts for the gate voltage and
the distance between the GNR and the gate; the second one is the amplitude
factor and it accounts for the amplitude of the $k$ mode oscillation (in this
case given only by the temperature), and the last one is the scaled emissivity
which is independent of the other quantities. To simplify the notation, from
now on the scaled emissivity will be denoted as $\tilde{I}_{k}$, i.e.,

$\tilde{I}_{k}\equiv\left(\frac{\text{d}\tilde{n}_{r}^{(E)}}{\text{d}q_{k}}\right)_{\varepsilon_{\text{F}},\bm{q}_{0}},$
(23)

where we omitted the $r$ subindex since it is irrelevant for the present case
with equal contacts and zero bias voltage. This quantity, which has units of
one over energy, can be transformed into a true current simply by multiplying
it by the scaling and amplitude factors, see Eq. (22).

##  III Pumping currents induced by oscillations of the GNR

The main quantity to be discussed in Secs. III.1, III.2, and III.3, and
plotted in Figs. 1, 3, 4, and 5, is the scaled emissivity $\tilde{I}_{k}$ in
the low temperature limit, where $\tilde{I}_{k}$ is proportional to the
maximum value of the pumping current during a cycle of the $k$-mode. In
sections III.1, III.2, and III.3 we mainly study the dependence of
$\tilde{I}_{k}$ on different characteristics of the GNRs, while in section
III.4 we discuss a realistic estimation of the maximum value of the vibration-
induced pumping current in typical GNRs. To carry out these tasks, we use the
theory and models previously presented, combining them with the numerical
tools provided by the KWANT package Groth  _et al._ (2014) for the calculation
of the scattering matrices. In the plots all energies are in eV while
distances are in units of $a\equiv\sqrt{3}a_{\text{cc}}=0.246$ nm, which is
the length of the primitive vectors of the associated Bravais lattice. Foa
Torres  _et al._ (2020)

###  III.1 Mode dependence

Here we study the effect of the mode oscillation on the quantum pumping
current through the calculation of $\tilde{I}_{k}$ for different values of
$n_{x}$ (in the range [1,10]) and $n_{y}$ (in the range [0,10]), taking into
account both zigzag and armchair edges. For the calculations we consider a
square membrane with length $L_{x}$ and width $L_{y}$ equal to $449\,a\sim
0.1$ $\mu$m. This value of $L_{y}$ ensures that the aGNR is metallic (see Sec.
II.1) and that the system is large enough to validate the rectangular membrane
model (Sec. II.2). Regarding the Fermi energy, the value
$\varepsilon_{\text{F}}=0.1$ eV was taken, which ensures that we are working
far from edge states for the zigzag case, which can cause the pumping current
to diverge (see Sec. III.2).

The dependence of $\tilde{I}_{k}$ on the excited mode is illustrated in Fig.
1(c) for armchair edges and in Fig. 1(d) for zigzag edges. It is clear, in
both cases, a relationship between the symmetry of the normal modes and the
pumping current: If $n_{x}$ is even or if $n_{y}$ is odd, the pumped current
is negligible. This behavior can be explained by analyzing the contribution of
each atom $l$ to the total emissivity of a mode $k$, that is

$\left(\frac{\text{d}\tilde{n}_{r}}{\text{d}q_{k}}\right)=\sum_{l}\left(\frac{\text{d}\tilde{n}_{r}}{\text{d}z_{l}}\right)\frac{\partial
z_{l}}{\partial q_{k}},$ (24)

where $(\text{d}\tilde{n}_{r}/\text{d}z_{l})$ is the emissivity given by
moving only atom $l$. If we now assume that the system is infinite and that
all atoms are equivalent, in the sense that their variation equally affects
the scattering matrix, we can write

$\left(\frac{\text{d}\tilde{n}_{r}}{\text{d}q_{k}}\right)\approx\left(\frac{\text{d}\tilde{n}_{r}}{\text{d}z_{l}}\right)\sum_{l}\frac{\partial
z_{l}}{\partial q_{k}}.$ (25)

Then, based on this assumption, we see that it is only the position of each
atom [$z_{l}\equiv z_{l}(x_{l},y_{l},\bm{q})$] and the normal mode involved
$q_{k}$ which determines its contribution to the total current. Thus, for
those vibrations where $n_{x}$ is even or $n_{y}$ is odd, any contribution to
the pumping current is canceled by a contribution of the opposite sign. It is
important to note that, even if the previous approximation is not completely
valid (for example, due to the finite size of the system), it can still be
expected in certain cases that the maximum value of the pumping current is
zero. This occurs when the system presents a reflection symmetry with respect
to a line that cuts the nanoribbon in half. In that case, for any pair of
sites $i$ and $j$ in mirror positions with respect to a line parallel (or
perpendicular) to the direction of transport, the emissivities per atom must
be equivalent, i.e.,

$\left(\frac{\text{d}\tilde{n}_{r}}{\text{d}z_{i}}\right)=\left(\frac{\text{d}\tilde{n}_{r}}{\text{d}z_{j}}\right).$
(26)

Given that for $n_{y}$ odd ($n_{x}$ even), the following relation holds

$\frac{\partial z_{i}}{\partial q_{k}}=-\frac{\partial z_{j}}{\partial
q_{k}},$ (27)

it is clear that $\tilde{I}_{k}=0$.

As can be seen in Fig. 1, the above arguments hold quite well for even $n_{x}$
on both types of ribbons. This is reasonable considering that the system is
truly infinite in the $x$ direction. Let us recall that we are using as
contacts semi-infinite nanoribbons of the same width and type as that of the
suspended system. In this way, the pumping current is strictly zero for
$n_{x}$ even, when the positions of the equivalent atoms coincide exactly with
a reflection with respect to a line perpendicular to the direction of
transport that cuts the system in half. The deviations from this are very
small for the values of $L$ used and therefore the pumping currents are
negligible, as can be seen. With respect to the $y$-direction, the system is
finite and, furthermore, for GNRs with zigzag edges there is no strict
reflection symmetry along this direction, see Fig. 2(a). For this reason,
deviations from our previous reasoning are expected, and this explains that
some current (although orders of magnitude smaller) can be appreciated for odd
$n_{y}$. GNRs with armchair edges, on the other hand, do present a reflection
symmetry along the $y$ axis, see Fig. 2(b), so it is expected that the pumping
currents are exactly zero for odd $n_{y}$.

We verified that the same type of behavior shown in Figs. 1(c) and (d) is also
repeated for different Fermi energies and other sizes of the GNR, as long as
the ribbon is not too small. Another general characteristic of the studied
system is that the lower frequency modes yield, by far, the largest
contribution to the pumping current. In particular, the mode that generates
the highest value of $\tilde{I}_{k}$ is the one with $n_{x}=1$ and $n_{y}=0$.
This is reasonable since, for this mode, the movement of all the sites
contributes with the same sign to the pumping current and, thus, there are no
cancellations. In Fig. 1(a) we depict the shape of this mode, while in Fig.
1(b) we show some other modes with negligible contributions to
$\tilde{I}_{k}$. Based on the previous results, from now on we will only work
with the mode with $n_{x}=1$ and $n_{y}=0$, denoting the scaled emissivity as
$\tilde{I}_{(1,0)}$.

###  III.2 Fermi energy dependence

In this section we will study how the maximum pumping current behaves for
different values of the Fermi energy, $\varepsilon_{\text{F}}$. Again, we will
consider both edges for the graphene nanoribbon. In Fig. 3 we show
$\tilde{I}_{(1,0)}$ as a function of the Fermi energy for square nanoribbons
of sides $L_{x}=L_{y}=449\,a$ and for the mode with $n_{x}=1$ and $n_{y}=0$.
In this figure it can be seen that for large Fermi energies the type of edge
does not seem to have an effect on the maximum value of the pumping current.
In both cases, the pumping current increases roughly linearly with the number
of conduction channels, which increases with the energy. For small Fermi
energies, however, there is a clear difference in the behavior of the curves.
In the case of aGNRs, $\tilde{I}_{(1,0)}$ converges to a fixed value as
$\varepsilon_{\text{F}}$ goes to zero. Let us recall that because of the
values of $L_{y}$ used, there is always one conduction channel at
$\varepsilon_{\text{F}}\approx 0$ for this type of nanoribbon. On the other
hand, for zGNRs, the maximum value of the current increases rapidly as
$\varepsilon_{\text{F}}$ approaches zero. In fact, for
$\varepsilon_{\text{F}}=0$, the pumping current diverges, causing numerical
problems in our calculations. This effect can be explained by noticing that
zGNRs possess edge states at this energy with zero group velocity, Foa Torres
_et al._ (2020) causing the divergence of the associated density of states.
The relation between the density of states and the pumped currents have been
established in many different contexts.Ingaramo and Foa Torres (2013); Deghi
_et al._ (2021)

In order to confirm that edge states are causing the abrupt increase of the
pumping currents, we evaluate the contribution of each row of atoms to the
total current by means of Eq. (24). In Fig. 2 it can be seen how these rows
are defined for each nanoribbon type (see the green area). Fig. 4 shows the
contribution per row to $\tilde{I}_{(1,0)}$ for $\varepsilon_{\text{F}}=0.001$
eV. As can be seen, for zGNR rows close to the edges (located at $y=0$ and
$y=449$) are the ones that contribute the most to $\tilde{I}_{(1,0)}$. On the
contrary, for aGNR rows near the edges do not seem to play any particular role
in the currents. It should be mentioned that, in real systems, a true
divergence of the pumping current is not expected since the coupling with the
environment (causing for example decoherence Pastawski and Medina (2001))
should regularize it. However, at least a sharp peak of the pumping current
should appear in potential experiments for $\varepsilon_{\text{F}}=0$. In our
case, without a specific model for the interaction with the environment, we
cannot estimate the value of the pumping current precisely at this value.
Finally, it is interesting to discuss the effect of finite temperatures in
Fig. 3. Since the result of this is the convolution of the emissivity with the
derivative of the Fermi-Dirac function [compare Eqs. (9) and (10)], its
expected effect is simply to smooth (or average) the curves over a range of
the order $k_{\text{B}}T$. 222 Take into account that for Fermi energies
$\varepsilon_{\text{F}}$ too close to the Dirac point
($|\varepsilon_{\text{F}}|\lesssim k_{\text{B}}T$) the calculation of pumping
currents at finite temperatures would present numerical issues due to the
divergence of emissivities there, see Eq. (9)

According to the above discussion, there are two ways of increasing the pumped
currents for zGNRs: either by adjusting the Fermi energy close to zero (and
thus taking advantage of edge states and the divergence of the density of
states that they cause), or by setting the Fermi energy to a large value
(which naturally includes more conduction channels through which the current
can be pumped). For aGNRs, the only available strategy is the second one.

Figure 3: Scaled emissivity $\tilde{I}_{(1,0)}$ in units of 1/eV, as a
function of the Fermi energy $\varepsilon_{\text{F}}$ for square GNRs with
$L=449\,a$. Figure 4: Contribution of each row of atoms [calculated using Eq.
(24)] to the scaled emissivity $\tilde{I}_{(1,0)}$, for square GNRs of side
$L=449\,a$, and $\varepsilon_{\text{F}}=0.001$ eV. The value of $y$ is the
position of the row, in units of $a$. See Fig. 2 for examples of rows for each
type of GNR.

###  III.3 Size dependence

In experiments with GNRs, these systems can sometimes have sizes on the order
of a micrometer Bolotin  _et al._ (2008a, b); Singh  _et al._ (2010); Ma  _et
al._ (2013). Since this involves a large number of atoms, performing numerical
calculations for systems with those sizes involves an enormous computational
cost. One way of estimating the value of the current $\tilde{I}_{k}$ for long
nanoribbons is to study the behavior of this quantity as a function of size to
make an extrapolation.

We calculate the maximum value of the pumping current for square nanoribbons
of different sizes $L=L_{x}=L_{y}$ and for both types of edges. The particular
widths used ($L_{y}$ values) were chosen such as to guaranty that the aGNRs
are metallic. Based on the previous analysis, we took $n_{x}=1$ and $n_{y}=0$,
and chose $\varepsilon_{\text{F}}=0.1$ eV for the calculations. This value of
$\varepsilon_{\text{F}}$ is far enough from the Dirac point (thus avoiding
numerical problems) and is, in principle, experimentally accessible. Zhang
_et al._ (2020)

In Fig. 5 we see, for both edges, that the log-log plot of $|\tilde{I}_{k}|$
with the size of the nanoribbons follows an approximate linear behavior for
large values of $L$. There are also some dispersion around this linear
behavior. However, this is not surprising since we are using the zero
temperature limit of $\tilde{I}_{k}$, which is very sensitive to the density
of states at the Fermi energy. Therefore, every time a new conduction channel
appears at this energy (consequence of an increase in $L$), a jump in the
density of states is expected (due to the emergence of van-Hove singularities)
causing also a jump in $\tilde{I}_{k}$. The fittings shows that both edges
follow an approximate quadratic dependence with $L$, for both edges. That is,
the pumping currents scale with the surface of the system ($L^{2}$), and not
with the width ($L$), as it is the case of the zero-order current $I^{(0)}$
(proportional to the number of conduction channels at the studied energy).

According to the above results, we can estimate the value of $\tilde{I}_{k}$
for a square nanoribbon with sides close to one micrometer, $L\approx
4065\,a$. This gives an extrapolated scaled emissivity of approximately
$\tilde{I}_{k}\sim 4\times 10^{4}\,\text{eV}^{-1}$. This value will be useful
in the following section.

Figure 5: Scaled emissivity $\tilde{I}_{(1,0)}$ as a function of the system’s
size $L$ for square GNRs, with $\varepsilon_{\text{F}}=0.1$ eV. The linear
fits in dotted lines, show the approximate quadratic dependence between
$\tilde{I}_{(1,0)}$ and $L$. $\tilde{I}_{(1,0)}=10^{-2.01}L^{2.11}$ for aGNR
and $\tilde{I}_{(1,0)}=10^{-2.97}L^{2.09}$ for zGNRs.

###  III.4 Realistic estimation of pumping currents

With the information collected so far, we are able to estimate the maximum
value of the pumping current, $I_{L,(1,0),\text{max}}^{(1)}$, in a realistic
context. Just for the sake of comparison, we will take the amplitude factor as
given by temperature only, see Eq. (22).

For the reasons discussed in the previous sections we choose $n_{x}=1$,
$n_{y}=0$, $\varepsilon_{\text{F}}=0.1$ eV and $L\approx 4065\,a$. With these
parameters, the extrapolated value of the scaled emissivity in the low
temperature limit yields between $4.03\times 10^{5}\,\text{eV}^{-1}$
(extrapolating for aGNRs) and $3.74\times 10^{4}\,\text{eV}^{-1}$
(extrapolating for zGNRs). According to Eq. (22), we see that we still need to
specify $d_{0}$, $T$, $V_{g}$ and $A_{\text{site}}$ to obtain a concrete value
for the current. The factor $A_{\text{site}}$ is simply the area of the
graphene unit cell divided by two, $A_{\text{site}}=0.026$ nm2. For the rest
of the parameters we took Ref. [Bolotin  _et al._ , 2008a] as a guide , fixing
the rest of the quantities at: $d_{0}=150$ nm, $T=5$ K, and $V_{g}=1$ V.

By substituting all the above values in Eq. (22), we obtain a value of
$I_{L,(1,0),\text{max}}^{(1)}$ between $1$ and $10$ pA. These contributions
can be increased by taking advantage of the dependence of
$I_{L,(1,0),\text{max}}^{(1)}$ on $d_{0}$. For example, reducing the
separation to $d_{0}=15$ nm, it follows that $I_{L,(1,0),\text{max}}^{(1)}$
lies between $0.1$ and $1$ nA.

Just to put into context the above values of pumping currents we can compare
them with some measured values of currents in nanodevices. In Ref. [Thibado
_et al._ , 2020], fluctuation-induced current of the order of 10 pA are
observed working with micrometric sheets of suspended graphene at temperatures
two orders of magnitude higher than that used in our estimate. Another example
is that of irradiated graphene experiments such as the one in Ref. [Higuchi
_et al._ , 2017], in which electric currents of the order of 10 pA at room
temperature have also been measured. Working with mechanical resonators based
on monolayer graphene, the authors of Ref. [Chen  _et al._ , 2009] measured
currents of the order of pA at temperatures of 5 K. In conclusion, the
estimated values of pumping currents induced by thermal vibrations seem
plausible of being measured and, in principle, their effect should not be
neglected since they can potentially interfere in some experiments.

Before continuing, here a comment is in order. The formalism used to describe
the pumping current assumes that the vibrational degree of freedom behaves
classically. This is valid as long as the vibrational energy is far from the
zero point energy, which for thermal excitation implies
$\hbar\omega_{k}/2k_{\text{B}}\ll T$. Considering a GNR with $L_{x}\sim
10^{-6}$ m, $n_{x}=1$, $n_{y}=0$, and taking $v\approx 12.9$ km/s, Cong  _et
al._ (2019) we estimate $\hbar\omega_{(1,0)}/2k_{\text{B}}\sim 0.05$ K. This
very small temperature value also shows that the application of the
equipartition theorem for the mechanical degrees of freedom used in Sec. II.2
does not contradict the low temperature limit for the energy dependence of the
electronic scattering matrix used in Eq. (10).

Finally, it is important to highlight that the presented estimations
correspond to pumping currents induced by thermal excitations of the
vibrational modes. Other forms of ambient vibrations, such as sound waves
traversing the sample, could give rise to considerably larger amplitudes of
the normal modes [increasing the amplitude factor in Eq. (22)] and thus much
larger values of $I_{L,(1,0),\text{max}}^{(1)}$ are possible. Therefore, the
above estimations are the minimal expected values of the pumping currents in a
real case scenario.

##  IV Current-noise induced by thermal oscillations of GNRs

In previous sections, we show that fluctuation-induced pumping currents can be
strong enough to be measurable, even for thermal excitations of the
nanoribbons, affecting the observed currents of different nanodevices.
Considering that stochastic fluctuations of the nanoribbons should induce what
is essentially noise in the current determinations, it is fair to wonder how
strong is this new form of noise, which we dubbed current-noise induced by
thermal vibrations (CNITV), as compared with more standard forms of current
noise.

In this section, we develop a semiclassical approach to evaluate CNITV. This
form of current-noise is associated with the first-order terms of the
adiabatic expansion of the current $\hat{I}^{(1)}$ (or pumping current), and
it should not be confused with Nyquist-Johnson or shot noise which are due to
zero order terms of the current operator $\hat{I}^{(0)}$. In particular, we
are wondering what is the correlation function of the pumping currents between
contacts $\alpha$ and $\beta$, $S_{\alpha\beta}^{(1)}(\tau)$, induced by
stochastic thermal variations of the parameters of a system’s Hamiltonian.
Here, we evaluate this quantity by means of

$S_{\alpha\beta}^{(1,\text{sc})}(\tau)=\Braket{\left\\{\Delta
I_{\alpha}^{(1)}\left(t\right)\Delta
I_{\beta}^{(1)}\left(t^{\prime}\right)\right\\}_{\text{s}}}_{\text{c}},$ (28)

where $\tau=t-t^{\prime}$,
$\\{A_{\alpha}A_{\beta}\\}_{\text{s}}=(A_{\alpha}A_{\beta}+A_{\beta}A_{\alpha})/2$,
$\braket{...}_{\text{c}}$ stands for the average over stochastic trajectories
of the classical parameters of the Hamiltonian, 333The electronic Hamiltonian
$H$ varies when the GNR oscillates with a given normal mode (as explained in
Sec. II.4). Let us recall that the movement of the GNR is considered
classically here, so $\bm{q}$ (the vector containing the amplitude of the
normal modes) is a parameter of $H$. In our work, we are considering small
amplitudes of the normal modes, so $\bm{q}(t)\approx\bm{q}_{0}$. However,
$\bm{\dot{q}}$ will still change stochastically due to thermal fluctuations.
Therefore, when calculating the current noise one has to average over the
stochastic variation of $\bm{\dot{q}}$, this is $\Braket{...}_{\text{c}}$,
which is different from $\Braket{...}$, the expectation value of some quantity
and $\Delta
I_{\alpha}^{(1)}=I_{\alpha}^{(1)}-\braket{I_{\alpha}^{(1)}}_{\text{c}}$. The
pumping current $I_{\alpha}^{(1)}$ is the quantum expectation value
$\braket{...}$ of the first order current operator $\hat{I}_{\alpha}^{(1)}$,
i.e., $I_{\alpha}^{(1)}=\braket{\hat{I}_{\alpha}^{(1)}}$. This expectation
value if zero once averaged over stochastic trajectories, i.e.,
$\braket{I_{\alpha}^{(1)}}_{\text{c}}=\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\hat{I}_{\alpha}^{(1)}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\text{c}}=0$.
Therefore, $\Delta I_{\alpha}^{(1)}$ is simply $I_{\alpha}^{(1)}$.

Note that Eq. (28) corresponds to a semiclassical approximation of the fully
quantum current-noise Polianski  _et al._ (2002); Riwar  _et al._ (2013)

$S_{\alpha\beta}^{(1)}(\tau)=\mathopen{\hbox{\set@color${\Big{\langle}}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\Big{\langle}}$}}\left\\{\Delta\hat{I}_{\alpha}^{(1)}(t)\Delta\hat{I}_{\beta}^{(1)}(t^{\prime})\right\\}_{\text{s}}\mathclose{\hbox{\set@color${\Big{\rangle}}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\Big{\rangle}}$}}_{\text{c}},$
(29)

where we replaced the current operators by their quantum expectation values.
In this sense, at least for $\alpha=\beta$, Eq. (28) can be thought of as a
lower limit to $S^{(1)}(\tau)$. First, note that their difference $\Delta
S=S_{\alpha\alpha}^{(1)}(\tau)-S_{\alpha\alpha}^{(1,\text{sc})}(\tau)$ gives

$\Delta
S=\left\\{\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\hat{I}_{\alpha}^{(1)}(t)\hat{I}_{\alpha}^{(1)}(t^{\prime})\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\text{c}}-\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\hat{I}_{\alpha}^{(1)}(t)\rangle\langle\hat{I}_{\alpha}^{(1)}(t^{\prime})\mathclose{\hbox{\set<EMAIL_ADDRESS>(30)

Now, for $\tau=0$, it holds

$\Braket{\left[\hat{I}_{\alpha}^{(1)}(t)\right]^{2}}-\Braket{\hat{I}_{\alpha}^{(1)}(t)}^{2}\geq
0,$ (31)

due to the Cauchy–Schwarz inequality, while for $\tau\rightarrow\infty$,

$\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\hat{I}_{\alpha}^{(1)}(t)\hat{I}_{\alpha}^{(1)}(t^{\prime})\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\text{c}}=\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\hat{I}_{\alpha}^{(1)}(t)\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\text{c}}=0,$
(32)

since each trajectory of the classical parameters over which we are averaging
is stochastic. Therefore, $S_{\alpha\alpha}^{(1)}\geq
S_{\alpha\alpha}^{(1,\text{sc})}$ for $\tau=0$ and
$S_{\alpha\alpha}^{(1)}=S_{\alpha\alpha}^{(1,\text{sc})}=0$ for
$\tau\rightarrow\infty$.

In our case, the time-dependent variation of the acoustic transversal modes of
graphene $q_{i}(t)$ is the origin of the current noise, and thus

$\Delta
I_{\alpha}^{(1)}(t)=e\sum_{k}\left(\frac{\text{d}n_{\alpha}}{\text{d}q_{k}}\right)\left(\dot{q}_{k}(t)-\braket{\dot{q}_{k}(t)}_{\text{c}}\right).$
(33)

Replacing this in Eq. (28), taking into account that velocities of different
modes are not correlated, i.e.,
$\braket{\dot{q}_{k}\dot{q}_{k^{\prime}}}_{\text{c}}=\braket{\dot{q}_{k}}_{\text{c}}\braket{\dot{q}_{k^{\prime}}}_{\text{c}}$,
and that their stochastic average is zero,
$\braket{\dot{q}_{k}(t)}_{\text{c}}=0$, we obtain

$S_{\alpha\beta}^{(1,\text{sc})}(\tau)=e^{2}\sum_{k}\left(\frac{\text{d}n_{\alpha}}{\text{d}q_{k}}\right)\\!\left(\frac{\text{d}n_{\beta}}{\text{d}q_{k}}\right)\braket{\dot{q}_{k}(t)\dot{q}_{k}(t^{\prime})}_{\text{c}}.$
(34)

Now we will assume a generic model for the autocorrelation function of
velocities

$\braket{\dot{q}_{k}(t)\dot{q}_{k}(t^{\prime})}_{\text{c}}\approx\frac{k_{\text{B}}T}{m_{\text{c}}}e^{-|t-t^{\prime}|/\tau_{k}},$
(35)

which fulfills equipartition theorem for $t=t^{\prime}$, while including the
effects of a finite correlation time $\tau_{k}$ for each normal mode $k$. This
kind of exponentially decaying correlation function can indeed be derived
analytically for simple models, see for example Ref. [Balakrishnan, 2021].
Inserting Eq. (35) into Eq. (34) and taking the limit of small $\tau_{k}$
gives 444This limit means that correlations decay instantaneously as compared
with the measurement time. This would lead $S(\tau)$ to be proportional to the
Dirac delta function (with a proportionality factor equal to the so called
zero-frequency noise), which is an usual assumption in noise treatment.

$S_{\alpha\beta}^{(1,\text{sc})}(\tau)\approx\frac{2e^{2}k_{\text{B}}T}{m_{\text{c}}}\left[\sum_{k}\left(\frac{\text{d}n_{\alpha}}{\text{d}q_{k}}\right)\\!\left(\frac{\text{d}n_{\beta}}{\text{d}q_{k}}\right)\tau_{k}\right]\delta(\tau),$
(36)

where the term in brackets is the zero-frequency noise of the CNITV,
$S_{\alpha\beta}^{(1,\text{sc})}(\omega=0)$. It is interesting to compare this
term with the Nyquist-Johnson zero-frequency noise,
$S_{\alpha\beta}^{(\text{NJ})}(\omega=0)$, typically used to assess current-
noise. In the limit of small temperature, the latter is given as Blanter and
Büttiker (2000)

$S_{\alpha\beta}^{(\text{NJ})}(\omega=0)=-4k_{\text{B}}T\frac{e^{2}}{h}\sum_{\eta\in\alpha,\eta^{\prime}\in\beta}T_{\eta\eta^{\prime}},$
(37)

where $T_{\eta\eta^{\prime}}$ is the transmission coefficient between
conduction channels $\eta$ and $\eta^{\prime}$ belonging to contacts $\alpha$
and $\beta$, respectively, and
$T_{\eta\eta}=-\sum_{\eta\neq\eta^{\prime}}T_{\eta\eta^{\prime}}$. If we
wonder when $S_{\alpha\alpha}^{(1,\text{sc})}(\omega=0)$ will be the dominant
source of current noise, i.e.,
$S_{\alpha\alpha}^{(1,\text{sc})}>S_{\alpha\alpha}^{(\text{NJ})}$, this
implies the condition, in the low temperature limit,

$\frac{h}{2m_{\text{c}}}\left(\frac{\partial E}{\partial
z}\right)^{2}\tilde{I}_{(1,0)}^{2}\tau_{(1,0)}\gtrsim-\sum_{\eta,\eta^{\prime}\in\alpha}T_{\eta\eta^{\prime}},$
(38)

where the right-hand side of the equation is the total transmittance and, in
agreement with the previous sections, we assume that the mode
$(n_{x},n_{y})=(1,0)$ gives the dominant contribution to the pumped current.

Typically, defects in quasi-one-dimensional systems like GNRs reduces
dramatically the transmittance while pumping currents are not necessarily
affected in the same way. Indeed, the peaks in the density of states caused by
defects can favorably affect the pumping currents. Ingaramo and Foa Torres
(2013) The scaling with the length of the sample also favors pumping currents.
Note that the emissivity is approximately proportional to the length of the
system [as Eq. (24) suggests], while the total transmittance in the coherent
limit of the current is, at best, independent of the length of the sample. Let
us recall that, in quasi-one-dimensional systems with defects, the
transmittance is exponentially suppressed with the length. Cattena  _et al._
(2010); Fernández-Alcázar  _et al._ (2019)

The above reasoning sets some general conditions where one would expect the
CNITV to be the dominant contribution to the current noise: elongated samples
with defects. For example, we consider a GNR with $L_{x}=800\,a$,
$L_{y}=200\,a$ (which gives $\omega_{(1,0)}=1.9\times 10^{11}$ rad/s), 0.5% of
vacancies randomly distributed, $\varepsilon_{\text{F}}=0.1$ eV, $d_{0}=15$
nm, and $V_{g}=1$ V. With these parameters, we obtained, after 70 independent
runs giving different pairs of scaled emissivities and transmittances, that in
83% of the samples $\tau_{(1,0)}=1.05\times 10^{-8}$ s is enough to make the
CNITV the dominant contribution to current noise. This means that a quality
factor (taking $Q=\tau_{k}\omega_{k}/2$) of only $10^{3}$ is required, which
is a modest value compared with what is possible experimentally for graphene
resonators at low temperatures. Takamura  _et al._ (2016) Increasing vacancies
up to 1% makes that a quality factor of $10^{3}$ suffices for 92% of the
samples.

##  V Conclusions

Motivated by the implementation of GNRs as pumps, sensors, and recent
experiments on thermal fluctuation-induced currents, here we developed a
theoretical and numerical model to study the instantaneous value of quantum
pumping currents in these systems. In particular, we focus on calculating the
maximum value of the pumping current that arises from the excitation of the
transverse acoustic modes of suspended GNRs. Although our calculations are
strictly valid under certain conditions (namely classical limit of vibrations,
small oscillations, low frequency transversal acoustic modes, large ribbon
sizes, and the low temperature limit), our results can be helpful in
estimating the magnitude of this type of currents. In this context, the
developed theoretical tools and our results allow one to understand the role
of the different parameters in the fluctuation-induced currents and assess
under which conditions these currents can play an important role in actual
experiments.

As particular results of the geometry studied, we found that: Pumping currents
increase with the size of the system, being approximately proportional to the
area of the ribbon and independent of the edge geometry for large sizes. The
lowest frequency modes contribute the most to the pumping currents. The
contribution of modes with $n_{x}$ even or $n_{y}$ odd is negligible, being
strictly zero when the size of the system tends to infinity. While for large
Fermi energies there are no significant differences between zigzag and
armchair GNRs, for energies close to the Dirac point marked divergences
appear. The behavior of aGNRs depends on their width, which makes them
metallic or semiconductor. Metallic aGNRs present a relatively small and
energy-independent pumping current, while semiconductor aGNRs have zero
pumping current for Fermi energy close to the Dirac point. On the contrary,
zGNRs present a strong increase of pumping currents close to the Dirac point
due to the presence of their edge states.

Our estimations of GNRs’ pumping currents in the studied geometries indicate
that, in principle, even thermal fluctuations of the ribbons should produce
measurable currents and that their effects should not be negligible for
several experimental setups. It is worth noticing that other forms of
excitation of the GNRs, such as propagating sound waves traversing the system,
are potentially even more favorable for inducing pumping currents.

Finally, we developed a semiclassical theory to evaluate the current noise
induced by thermal vibrations (CNITV). The comparison with Nyquist-Johnson
noise allowed us to set some general conditions where one would expect the
CNITV to be the dominant contribution to the current noise: elongated samples
with defects.

Our results, which include quantum effects, could provide a more accurate
calculation of fluctuation-induced currents in experiments like the one in
Ref. [Thibado  _et al._ , 2020]. Given that sensing is ultimately limited by
the signal-to-noise ratio, our results can also help to improve the design of
sensors based on suspended GNRs. This could be done, e.g., in optimal
configurations that minimize current noise and/or maximize the signal of the
fluctuation-induced current coming from a specific oscillation of the ribbon.
If fluctuating currents arising from the oscillation of GNRs could be
rectified (that would depend on the frequency and the magnitude of the
current), then our results could be used to design energy harvesters that turn
ambient mechanical excitations (above the thermal level) into electric power.

There are several further aspects that could be studied, including other
vibrational modes, and different systems like carbon nanotubes or
heterogeneous multilayer 2D materials. To assess the role of anharmonicities
in the pumping currents, large oscillation amplitudes of the modes can be
explored. Also, the effect of finite bias voltages on the induced current, or
even the consequences of local voltage thermal fluctuations, give rise to
interesting directions to follow. In this regard, the present work paves the
way to explore these and other exciting directions.

##  VI Acknowledgments

This project was supported by Consejo Nacional de Investigaciones Científicas
y Técnicas (CONICET, PIP-2022-59241); Secretaría de Ciencia y Tecnología de la
Universidad Nacional de Córdoba (SECyT-UNC, Proyecto Formar 2020); and Agencia
Nacional de Promoción Científica y Tecnológica (ANPCyT, PICT-2018-03587). We
acknowledge the useful comments from Lucas Fernández-Alcázar and Sebastián
Deghi. FDR, HLC, and RABM thank Luis E. F. Foà Torres for fruitful discussions
and hospitality during the 2023 ICTP conference at Universidad de Chile,
Santiago de Chile.

##  Appendix A Transverse normal modes of a rectangular elastic membrane

As mentioned in Sec. II.2, to study the vibrations in a GNR we will consider
that the suspended system is large enough to be approximated by a rectangular
elastic membrane. Assuming that the equilibrium position of the membrane
defines the $xy$ plane, we will call $z(x,y,t)$ the displacement of the
membrane with respect to said plane. This variable $z$ satisfies the wave
equation

$\nabla^{2}z=\frac{1}{v^{2}}\frac{\partial^{2}z}{\partial t^{2}},$ (39)

where $v$ is the speed of propagation of the wave. Making the replacement
$z=F(x,y)q(t)$, we separate the above equation into a spatial equation
(Helmholtz equation) and a temporal equation:

$(\nabla^{2}+k^{2})F=0,\quad\ddot{q}+k^{2}v^{2}q=0.$ (40)

The solutions to the temporal equation are of the form

$q(t)=A\cos(kvt)+B\sin(kvt),$ (41)

while for the spatial part it is necessary to specify the boundary conditions.
Once the corresponding eigenfunctions have been obtained, the general time-
dependent solution, $z(x,y,t)$, can be constructed by superposition. Going
back to the Helmholtz equation, it has the form

$\frac{\partial^{2}F}{\partial x^{2}}+\frac{\partial^{2}F}{\partial
y^{2}}+k^{2}F=0.$ (42)

If the boundary conditions are separated in the $x$ and $y$ directions, we can
apply separation of variables in the form $F(x,y)=X(x)Y(y)$. Substituting this
into Eq. (42) and dividing by $F$ gives

$k^{2}+\frac{1}{X}\frac{\text{d}X^{2}}{\text{d}x^{2}}=-\frac{1}{Y}\frac{\text{d}Y^{2}}{\text{d}y^{2}}.$
(43)

The left hand side is a function of $x$ only, while the right hand side is a
function of $y$ only. Therefore, this relation is only valid if both sides are
equal to a constant. With this in mind, we write

$\frac{\text{d}X^{2}}{\text{d}x^{2}}=-k_{x}^{2}X,$ (44)

and

$\frac{\text{d}Y^{2}}{\text{d}y^{2}}=-k_{y}^{2}Y,$ (45)

with $k_{x}$ and $k_{y}$ satisfying the relation $k_{x}^{2}+k_{y}^{2}=k^{2}$.
Thus, the two-dimensional problem has been reduced to two one-dimensional
problems, whose solutions are

$\displaystyle X(x)$ $\displaystyle=C_{x}\cos(k_{x}x)+D_{x}\sin(k_{x}x),$
$\displaystyle Y(y)$ $\displaystyle=C_{y}\cos(k_{y}y)+D_{y}\sin(k_{y}y).$ (46)

We are now in position to specify the boundary conditions of the problem. We
will think of a rectangular membrane with side $L_{x}$ on the $x$-axis and
side $L_{y}$ on the $y$-axis, with the two edges in the $x$-direction fixed
and the two edges in the $y$-direction free. This implies that the separate
solutions of Eqs. (44) and (45) must comply with $X(0)=X(L_{x})=0$ and
$Y^{\prime}(0)=Y^{\prime}(L_{y})=0$, to ensure that the membrane does not move
at the fixed edges and that there is zero slope at the free edges. Applying
these conditions, we arrive at

$C_{x}=0,\quad k_{x}=\frac{n_{x}\pi}{L_{x}},\quad n_{x}=0,1,2,\ldots,$ (47)

and

$D_{y}=0,\quad k_{y}=\frac{n_{y}\pi}{L_{y}},\quad n_{y}=0,1,2,\ldots,$ (48)

and in consequence:

$X(x)=\sin(k_{x}x),\quad Y(y)=\cos(k_{y}y).$ (49)

In conclusion, the solution of the spatial problem is of the form

$F_{n_{x},n_{y}}(x,y)=\sin\left(\frac{n_{x}\pi}{L_{x}}x\right)\cos\left(\frac{n_{y}\pi}{L_{y}}y\right),$
(50)

with

$k_{n_{x},n_{y}}^{2}=\frac{\omega_{n_{x},n_{y}}^{2}}{v^{2}}=\left(\frac{n_{x}\pi}{L_{x}}\right)^{2}+\left(\frac{n_{y}\pi}{L_{y}}\right)^{2}.$
(51)

The general solution to the two-dimensional wave equation can then be obtained
by superposition of the normal modes, resulting in this case in Eq. (2) of the
main text:

$z(x,y,t)=\sum_{n_{x}=0}^{\infty}\sum_{n_{y}=0}^{\infty}\sin\left(\frac{n_{x}\pi}{L_{x}}x\right)\cos\left(\frac{n_{y}\pi}{L_{y}}y\right)q(t).$
(52)

##  Appendix B Energy of normal modes

The total energy of a system composed of atoms labeled by $i$, and considering
only the $z$ component of their respective positions for simplicity, is

$E_{\text{tot}}=\sum_{i}m_{i}\frac{\dot{z}_{i}^{2}}{2}+U_{\text{tot}}(\bm{z}),$
(53)

where $m_{i}$ is the mass of the atom $i$ and $U_{\text{tot}}$ is the total
potential energy of the system. As usual in any harmonic approximation we
expand the potential up to second order in the displacement around the
equilibrium position $\bm{z}_{0}$ (here sets equal to zero for simplicity).
Using this, we obtain

$E_{\text{tot}}-U_{\text{tot}}(\bm{z}_{0})\approx\frac{m_{\text{c}}}{2}\left(\dot{\bm{z}}^{\text{T}}\dot{\bm{z}}+\frac{1}{m_{\text{c}}}\bm{z}^{\text{T}}\bm{\mathcal{H}}\bm{z}\right),$
(54)

where $\bm{\mathcal{H}}$ is the Hessian matrix. Now, assuming all masses are
equal ($m_{i}=m_{\text{c}}$) and inserting $\bm{U}^{\text{T}}\bm{U}=\bm{I}$,
where $\bm{U}\bm{\mathcal{H}}\bm{U}^{T}/m_{\text{c}}=\bm{\Omega}_{q}$ with
$\bm{\Omega}_{q}$ the diagonal matrix containing the square frequency of
normal modes, the energy results in

$E_{\text{tot}}-U_{\text{tot}}(\bm{z}_{0})\approx\frac{m_{\text{c}}}{2}\left(\dot{\bm{z}}^{\text{T}}\bm{U}^{\text{T}}\bm{U}\dot{\bm{z}}+\bm{z}^{\text{T}}\bm{U}^{\text{T}}\bm{\Omega}_{q}\bm{U}\bm{z}\right).$
(55)

Finally, we recognize $\bm{q}=\bm{U}\bm{z}$ as the coordinate vector of the
normal modes [which in our case we approximate as the $q_{k}$ of Eq. (2)].
Then, the energy can be written as

$E_{\text{tot}}-U_{\text{tot}}(\bm{z}_{0})\approx\frac{m_{\text{c}}}{2}\left(\dot{\bm{q}}^{\text{T}}\dot{\bm{q}}+\bm{q}^{\text{T}}\bm{\Omega}_{q}\bm{q}\right)=\sum_{k}E_{k},$
(56)

with $E_{k}=m_{\text{c}}(\dot{q}_{k}^{2}+\omega_{k}^{2}q_{k}^{2})/2$ being the
energy of the normal mode $k$.

##  Appendix C Validity of the adiabatic approximation

Adiabatic quantum pumping arises from the first-order correction (in a
frequency expansion) to the adiabatic approximation of an observable, the
current in our case. Technical details of its derivation can be found among
different contexts. Splettstoesser  _et al._ (2006); Arrachea and Moskalets
(2006); Bode  _et al._ (2011); Deghi  _et al._ (2021); Ribetto  _et al._
(2021) In the present case, the adiabatic approximation and its first-order
correction will be valid when the time it takes for the electrons to move
along the lattice sites is much shorter than the time it takes for the
Hamiltonian to change. To assess this, we take a mode with some value $n_{x}$
of a GNR of length $L_{x}$ ($n_{y}=0$ for simplicity). The time
$\tau_{\text{e}}$ electrons take to move (in a fluctuating GNR) between
regions with the maximum and minimum values of the onsite energies, is
$\tau_{\text{e}}=L_{x}/(2n_{x}v_{\text{F}})$, where $v_{\text{F}}$ is the
Fermi velocity. The frequency of the normal modes is given by Eq. (4), which
sets the variation frequency of the electronic Hamiltonian. Taking
$v_{\text{F}}=10^{6}$ m/s and $v_{\text{T}}=12.9\times 10^{3}$ m/s (the
velocity of sound for GNR transversal modes), yields a ratio between
$\tau_{\text{e}}$ and the period of the vibration
$\tau_{\text{mech}}=2\pi/\omega_{n_{x},n_{y}}$ of $3.23\times 10^{-3}$, which
is a reasonably small value for the expansion. Remarkably, for $n_{y}=0$ this
ratio is independent of the vibrational mode $n_{x}$, and given the expression
for $\tau_{\text{e}}/\tau_{\text{mech}}$ we do not expect large variations
from this ratio for small values of $n_{y}$.

##  Appendix D Parallel plate model

The goal of this appendix is to give an expression for the factor
$\partial_{z}E$. The key idea is using the principle of correspondence between
the quantum system described in Sec. II.1 and its classical analog taken as
the capacitor formed by a continuous conductive membrane (the GNR) separated
by a distance $d_{0}$ from the gate with a potential difference $V_{g}$. To
simplify the derivation, we will assume that the GNR can only move as a whole
in the $z$-direction, such that the classical analog is simply a parallel
plate capacitor with separation $d_{0}$ between the plates. In this way,
$z_{i}$ will be the same for all sites and the change in the onsite energies
will also be the same, i.e. $\partial_{z}\bm{H}^{(E)}=(\partial_{z}E)\bm{I}$,
with $\bm{I}$ the identity matrix. To simplify the notation, from now on we
will omit the superscript ($E$) since we will be always working with the part
of the Hamiltonian that represents the variations in the site energies of the
system. Now, the question that we want to answer is the following: What is the
value of $\partial_{z}E$ in the tight-binding Hamiltonian such that the system
behaves as a parallel plate capacitor in the limit of large membranes?

Consider the energy of a capacitor with parallel plates of area $A$, a
potential $V_{g}$ and separated by a distance $d$. The energy stored in this
“classical” capacitor is

$E^{(\text{classic})}=\frac{1}{2}CV_{g}^{2}=\frac{\epsilon_{0}A}{2d}V_{g}^{2}.$
(57)

If we now think of small $z$-displacements from a given equilibrium position
$d_{0}$, we can approximate the new energy as

$\displaystyle E^{(\text{classic})}$ $\displaystyle\approx
E_{0}^{(\text{classic})}+\left.\frac{\partial E^{(\text{classic})}}{\partial
z}\right|_{d_{0}}z,$
$\displaystyle=E_{0}^{(\text{classic})}+\left(-\frac{\epsilon_{0}A}{2d_{0}^{2}}V_{g}^{2}\right)z.$
(58)

Now considering the Hamiltonian model of the system, the total energy (near
equilibrium) is

$\displaystyle\braket{\hat{H}}$ $\displaystyle=$
$\displaystyle\text{Tr}[\hat{H}\hat{\rho}^{(\text{eq})}]$ (59)
$\displaystyle\approx$
$\displaystyle\braket{\hat{H}_{0}}+\text{Tr}[(\partial_{z}\hat{H})\hat{\rho}^{(\text{eq})}]z+\text{Tr}[\hat{H}\partial_{z}\hat{\rho}^{(\text{eq})}]z$
$\displaystyle\approx$
$\displaystyle\braket{\hat{H}_{0}}+\left(\sum_{i}(\partial_{z}H_{ii})\rho_{ii}^{(\text{eq})}\right)z$
$\displaystyle+\left(\sum_{ij}H_{ij}\partial_{z}\rho_{ji}^{(\text{eq})}\right)z,$

where we perform an expansion in $z$ around the equilibrium position $z_{0}$.
Here $\langle\hat{H}_{0}\rangle$ is the average energy of the system at
equilibrium ($z=0$). Note that $\hat{\rho}^{(\text{eq})}$ is diagonal in the
energy basis, but not necessarily in the $i$ position basis.

To find $\partial_{z}E$ we are going to demand that, in the limit of large
systems, the change in the total energy of the system corresponds to that
expected for the parallel plate capacitor described above. This implies
working Eq. (59) and then performing a comparison with Eq. (57) in the limit
of macroscopic systems.

First, starting from Eq. (59), we will show that the term with
$\partial_{z}\rho_{ji}^{(\text{eq})}$ is zero. In order to easily calculate
this amount, we are going to assume that

$\hat{\rho}_{\text{tot}}^{(\text{eq})}=\hat{\rho}_{\text{leads}}^{(\text{eq})}\otimes\hat{\rho}_{\text{sys}}^{(\text{eq})},\quad\hat{\rho}_{\text{sys}}^{(\text{eq})}=\frac{\exp[-\beta(\hat{H}_{\text{sys}}-\mu\hat{N})]}{\Phi_{\text{sys}}},$
(60)

where $\Phi_{\text{sys}}=\text{Tr}[\exp(-\beta(\hat{H}-\mu\hat{N}))]$ and
$\beta=1/k_{\text{B}}T$. Then, taking the Hamiltonian of the system [see Eq.
(12)] and setting $H_{\text{sys}}\equiv H$, we see that its derivatives are of
the form

$\partial_{z}\hat{H}=\left(\frac{\partial E}{\partial
z}\right)\hat{1}\quad\Rightarrow\quad\partial_{z}\hat{H}^{n}=n\left(\frac{\partial
E}{\partial z}\right)\hat{H}^{n-1}.$ (61)

Similarly, it holds

$\displaystyle\partial_{z}(\hat{H}-\mu\hat{N})$ $\displaystyle=$
$\displaystyle\left(\frac{\partial E}{\partial z}\right)\hat{1},$ (62)

and then,

$\displaystyle\partial_{z}(\hat{H}-\hat{N})^{n}$ $\displaystyle=$
$\displaystyle n\left(\frac{\partial E}{\partial
z}\right)(\hat{H}-\mu\hat{N})^{n-1}.$ (63)

Using these results, it is possible to calculate the derivative of the
equilibrium density matrix

$\displaystyle\partial_{z}\rho_{ji}^{(\text{eq})}$ $\displaystyle=$
$\displaystyle\langle
j|\frac{\partial_{z}e^{-\beta(\hat{H}-\mu\hat{N})}}{\Phi}|i\rangle$ (64)
$\displaystyle+\langle
j|e^{-\beta(\hat{H}-\mu\hat{N})}\partial_{z}\left(\frac{1}{\Phi}\right)|i\rangle.$

Working with the first term, we get

$\displaystyle\langle
j|\partial_{z}\left(\frac{e^{-\beta(\hat{H}-\mu\hat{N})}}{\Phi}\right)|i\rangle$
$\displaystyle=$ $\displaystyle(-\beta)\left(\frac{\partial E}{\partial
z}\right)\rho_{ji}^{(\text{eq})}.$ (65)

On the other hand, working the second term we arrive at

$\displaystyle\langle
j|e^{-\beta(\hat{H}-\mu\hat{N})}\partial_{z}\left(\frac{1}{\Phi}\right)|i\rangle$
$\displaystyle=$ $\displaystyle\beta\left(\frac{\partial E}{\partial
z}\right)\rho_{ji}^{(\text{eq})}.$ (66)

Thus, joining the two previous results yields
$\partial_{z}\rho_{ji}^{(\text{eq})}=0$.

The fact that $\partial_{z}\rho_{ji}^{(\text{eq})}=0$ and the assumption
$\partial_{z}H=(\partial E/\partial z)\bm{I}$ allows us to write $\langle
H\rangle$ as

$\langle H\rangle\approx\langle H_{0}\rangle+\left(\frac{\partial E}{\partial
z}\right)n_{\text{sys}}^{(\text{eq})}z,$ (67)

where
$n_{\text{sys}}^{(\text{eq})}=\left(\sum_{i}\rho_{ii}^{(\text{eq})}\right)$ is
the mean number of particles within the system at equilibrium. Therefore, for
our Hamiltonian model to be consistent with macroscopic theory, we must do

$\displaystyle\underset{A\rightarrow\infty}{\lim}\left(-\frac{1}{2}\epsilon_{0}\frac{A}{d_{0}^{2}}V_{g}^{2}\right)z$
$\displaystyle=$
$\displaystyle\underset{A\rightarrow\infty}{\lim}\left(\frac{\partial
E}{\partial z}\right)n_{\text{sys}}^{(\text{eq})}z,$ (68)

which gives

$\displaystyle\left(\frac{\partial E}{\partial z}\right)$ $\displaystyle=$
$\displaystyle\dfrac{-\dfrac{1}{2}\dfrac{\epsilon_{0}V_{g}^{2}}{d_{0}^{2}}}{\underset{A\rightarrow\infty}{\lim}\dfrac{n_{\text{sys}}^{(\text{eq})}}{A}}.$
(69)

Note that $n_{\text{sys}}^{(\text{eq})}$ is the number of sites multiplied by
$1/2$ (since two electrons enter per site according to the Pauli exclusion
principle), and that the total area $A$ can be thought of as the number of
sites multiplied by the area occupied by each site, $A_{\text{site}}$. Then we
can write

$\underset{A\rightarrow\infty}{\lim}\frac{n_{\text{sys}}^{(\text{eq})}}{A}=\frac{1}{2A_{\text{site}}},$
(70)

and in consequence,

$\left(\frac{\partial E}{\partial
z}\right)=-\epsilon_{0}\frac{A_{\text{site}}V_{g}^{2}}{d_{0}^{2}}.$ (71)

##  Appendix E Pumping currents due to coupling variations

In this appendix we will show that, in the limit of small oscillations of
flexural modes, the contribution to the pumping current due to hopping
variations between neighboring atoms can be neglected. More specifically, this
involves deriving Eq. (19) and proving that $\partial\bm{S}^{(V)}/\partial
q_{k}$ tends to zero in the perturbative limit. To carry out this task, we
will use the matrix form of the Fisher-Lee formula Fisher and Lee (1981); Bode
_et al._ (2012); Bustos-Marún (2018)

$\bm{S}=\bm{I}-2i\bm{W}^{\dagger}\bm{G}^{R}\bm{W}.$ (72)

Here, $\bm{G}^{R}$ is the retarded Green function given by

$\bm{G}^{R}=\underset{\eta\rightarrow
0^{+}}{\lim}\left[(\varepsilon+i\eta)\bm{I}-\bm{H}-\bm{\Sigma}^{R}\right]^{-1},$
(73)

where $\bm{H}$ is the Hamiltonian of the system without the contacts,
$\bm{\Sigma}^{R}$ is the retarded self-energy due to the contacts, and
$\varepsilon$ is the energy of the electrons. The matrix $\bm{W}$ comes from
the relation

$\bm{\Gamma}_{\alpha}=\bm{W}^{\dagger}\bm{\Pi}_{\alpha}\bm{W},$ (74)

where $\bm{\Pi}_{\alpha}$ is the projection operator towards the channel
$\alpha$ of some reservoir $r$, and $\bm{\Gamma}_{\alpha}$ is the
contribution, due to channel $\alpha$, of the imaginary part of the self-
energy $\bm{\Sigma}^{R}$, i.e., $\bm{\Gamma}=-\text{Im}(\bm{\Sigma}^{R})$ and
$\bm{\Gamma}=\sum_{\alpha}\bm{\Gamma}_{\alpha}$.

Using Eq. (72) and
$\partial\bm{G}^{r}=-\bm{G}^{r}\partial\left[\bm{G}^{r}\right]^{-1}\bm{G}^{r}$
we find

$\displaystyle\frac{\partial\bm{S}}{\partial q_{k}}$ $\displaystyle=$
$\displaystyle(-2i)\bm{W}^{\dagger}\bm{G}^{r}\frac{\partial\bm{H}^{(E)}}{\partial
q_{k}}\bm{G}^{r}\bm{W}$ (75)
$\displaystyle+(-2i)\bm{W}^{\dagger}\bm{G}^{r}\frac{\partial\bm{H}^{(V)}}{\partial
q_{k}}\bm{G}^{r}\bm{W}$ $\displaystyle=$
$\displaystyle\frac{\partial\bm{S}^{(E)}}{\partial
q_{k}}+\frac{\partial\bm{S}^{(V)}}{\partial q_{k}},$

proving Eq. (19).

Now, using the chain rule

$\displaystyle\frac{\partial\bm{H}^{(V)}}{\partial q_{k}}$ $\displaystyle=$
$\displaystyle\sum_{\ell}\frac{\partial\bm{H}^{(V)}}{\partial
z_{\ell}}\frac{\partial z_{\ell}}{\partial q_{k}},$ (76)

where $z_{\ell}$ is the displacement of atom $\ell$ in the $z$ direction, and
studying the matrix elements of $\partial\bm{H}^{(V)}/\partial z_{\ell}$ when
the system is in equilibrium we arrive at that

$\left(\frac{\partial H_{ij}^{(V)}}{\partial
z_{\ell}}\right)_{\bm{q}_{0}}=\left[(z_{i}-z_{j})\frac{t_{0}b}{a_{\text{cc}}}\frac{\partial}{\partial
z_{\ell}}(z_{i}-z_{j})\right]_{\bm{q}_{0}}.$ (77)

Since evaluating at $\bm{q}_{0}$ implies taking $z_{i}=0$, $\forall i$, we
have that

$\left(\frac{\partial H_{ij}^{(V)}}{\partial
z_{\ell}}\right)_{\bm{q}_{0}}\sim\left(z_{i}-z_{j}\right)_{\bm{q}_{0}}=0.$
(78)

In short, for small oscillations around the equilibrium position, the hopping
variation does not contribute to the pumping current induced by flexural
modes.

## References

  * Lee  _et al._ (2008) C. Lee, X. Wei, J. W. Kysar, and J. Hone, “Measurement of the elastic properties and intrinsic strength of monolayer graphene,” Science 321, 385 (2008).
  * Chen  _et al._ (2008) J.-H. Chen, C. Jang, S. Xiao, M. Ishigami, and M. S. Fuhrer, “Intrinsic and extrinsic performance limits of graphene devices on SiO2,” Nature Nanotechnology 3, 206 (2008).
  * Chen  _et al._ (2011) S. Chen, A. L. Moore, W. Cai, J. W. Suk, J. An, C. Mishra, C. Amos, C. W. Magnuson, J. Kang, L. Shi, and R. S. Ruoff, “Raman measurements of thermal transport in suspended monolayer graphene of variable sizes in vacuum and gaseous environments,” ACS Nano 5, 321 (2011).
  * Bunch  _et al._ (2007) J. S. Bunch, A. M. van der Zande, S. S. Verbridge, I. W. Frank, D. M. Tanenbaum, J. M. Parpia, H. G. Craighead, and P. L. McEuen, “Electromechanical resonators from graphene sheets,” Science 315, 490 (2007).
  * Garcia-Sanchez  _et al._ (2008) D. Garcia-Sanchez, A. M. van der Zande, A. S. Paulo, B. Lassagne, P. L. McEuen, and A. Bachtold, “Imaging mechanical vibrations in suspended graphene sheets,” Nano Letters 8, 1399 (2008).
  * Chen  _et al._ (2009) C. Chen, S. Rosenblatt, K. I. Bolotin, W. Kalb, P. Kim, I. Kymissis, H. L. Stormer, T. F. Heinz, and J. Hone, “Performance of monolayer graphene nanomechanical resonators with electrical readout,” Nature Nanotech 4, 861 (2009).
  * Croy  _et al._ (2012) A. Croy, D. Midtvedt, A. Isacsson, and J. M. Kinaret, “Nonlinear damping in graphene resonators,” Phys. Rev. B 86, 235435 (2012).
  * Zhu and Chen (2009) R. Zhu and H. Chen, “Quantum pumping with adiabatically modulated barriers in graphene,” Applied Physics Letters 95, 122111 (2009).
  * Prada  _et al._ (2009) E. Prada, P. San-Jose, and H. Schomerus, “Quantum pumping in graphene,” Phys. Rev. B 80, 245414 (2009).
  * Foa Torres  _et al._ (2011) L. E. F. Foa Torres, H. L. Calvo, C. G. Rocha, and G. Cuniberti, “Enhancing single-parameter quantum charge pumping in carbon-based devices,” Applied Physics Letters 99, 092102 (2011).
  * San-Jose  _et al._ (2011) P. San-Jose, E. Prada, S. Kohler, and H. Schomerus, “Single-parameter pumping in graphene,” Phys. Rev. B 84, 155408 (2011).
  * Low  _et al._ (2012) T. Low, Y. Jiang, M. Katsnelson, and F. Guinea, “Electron pumping in graphene mechanical resonators,” Nano Letters 12, 850 (2012).
  * Ingaramo and Foa Torres (2013) L. H. Ingaramo and L. E. F. Foa Torres, “Quantum charge pumping in graphene-based devices: When lattice defects do help,” Applied Physics Letters 103, 123508 (2013).
  * Connolly  _et al._ (2013) M. Connolly, K. Chiu, and S. e. a. Giblin, “Gigahertz quantized charge pumping in graphene quantum dots.” Nature Nanotech 8, 417 (2013).
  * Zang  _et al._ (2015) X. Zang, Q. Zhou, J. Chang, Y. Liu, and L. Lin, “Graphene and carbon nanotube (cnt) in mems/nems applications,” Microelectronic Engineering 132, 192 (2015).
  * Davaji  _et al._ (2017) B. Davaji, H. Cho, M. Malakoutian, J.-K. Lee, G. Panin, T. Kang, and C. Lee, “A patterned single layer graphene resistance temperature sensor.” Sci. Rep. 7, 8811 (2017).
  * Khan  _et al._ (2017) Z. H. Khan, A. R. Kermany, A. Öchsner, and F. Iacopi, “Mechanical and electromechanical properties of graphene and their potential application in mems,” Journal of Physics D: Applied Physics 50, 053003 (2017).
  * Moreno-Garcia  _et al._ (2022) D. Moreno-Garcia, X. Fan, A. D. Smith, M. C. Lemme, V. Messina, C. Martin-Olmos, F. Niklaus, and L. G. Villanueva, “A resonant graphene nems vibrometer,” Small 18, 2201816 (2022).
  * Gupta  _et al._ (2021) R. K. Gupta, F. H. Alqahtani, O. M. Dawood, M. Carini, A. Criado, M. Prato, S. K. Garlapati, G. Jones, J. Sexton, K. C. Persaud, C. Dang, U. Monteverde, M. Missous, R. J. Young, S. Boult, N. Dixon, L. Majewski, and M. A. Migliorato, “Suspended graphene arrays for gas sensing applications,” 2D Materials 8, 025006 (2021).
  * Liao and Koide (2011) M. Liao and Y. Koide, “Carbon-based materials: Growth, properties, mems/nems technologies, and mem/nem switches,” Critical Reviews in Solid State and Materials Sciences 36, 66 (2011).
  * Fan  _et al._ (2019) X. Fan, F. Forsberg, A. Smith, S. Schröder, S. Wagner, H. Rödjegård, A. Fischer, M. Östling, M. C. Lemme, and F. Niklaus, “Graphene ribbons with suspended masses as transducers in ultra-small nanoelectromechanical accelerometers.” Nat Electron 2, 394–404 (2019).
  * Ben-Shimon  _et al._ (2022) Y. Ben-Shimon, C. P. Sharma, C. J. Arnusch, and A. Ya’akobovitz, “Freestanding laser-induced graphene ultrasensitive resonative viral sensors,” ACS Applied Materials & Interfaces 14, 44713 (2022).
  * Johnson (1928) J. B. Johnson, “Thermal agitation of electricity in conductors,” Physical Review 32, 97 (1928).
  * Nyquist (1928) H. Nyquist, “Thermal agitation of electric charge in conductors,” Physical Review 32, 110 (1928).
  * Thibado  _et al._ (2020) P. M. Thibado, P. Kumar, S. Singh, M. Ruiz-Garcia, A. Lasanta, and L. L. Bonilla, “Fluctuation-induced current from freestanding graphene,” Physical Review E 102, 042101 (2020).
  * Brouwer (1998) P. W. Brouwer, “Scattering approach to parametric pumping,” Physical Review B 58, R10135 (1998).
  * Bode  _et al._ (2011) N. Bode, S. V. Kusminskiy, R. Egger, and F. von Oppen, “Scattering theory of current-induced forces in mesoscopic systems,” Physical Review Letters 107, 036804 (2011).
  * Bustos-Marún (2018) R. A. Bustos-Marún, “Geometric rectification for nanoscale vibrational energy harvesting,” Physical Review B 97, 075412 (2018).
  * Kim  _et al._ (2017) S. H. Kim, C. S. Haines, N. Li, K. J. Kim, T. J. Mun, C. Choi, J. Di, Y. J. Oh, J. P. Oviedo, J. Bykova, S. Fang, N. Jiang, Z. Liu, R. Wang, P. Kumar, R. Qiao, S. Priya, K. Cho, M. Kim, M. S. Lucas, L. F. Drummy, B. Maruyama, D. Y. Lee, X. Lepró, E. Gao, D. Albarq, R. Ovalle-Robles, S. J. Kim, and R. H. Baughman, “Harvesting electrical energy from carbon nanotube yarn twist,” Science 357, 773 (2017).
  * Foa Torres  _et al._ (2020) L. Foa Torres, S. Roche, and J. Charlier,  _Introduction to Graphene-Based Nanomaterials: From Electronic Structure to Quantum Transport_ (Cambridge University Press, 2020).
  * Cresti  _et al._ (2008) A. Cresti, N. Nemec, B. Biel, G. Niebler, F. Triozon, G. Cuniberti, and S. Roche, “Charge transport in disordered graphene-based low dimensional materials,” Nano Research 1, 361 (2008).
  * Dutta and Pati (2010) S. Dutta and S. K. Pati, “Novel properties of graphene nanoribbons: a review,” Journal of Materials Chemistry 20, 8207 (2010).
  * Anderson (1978) P. W. Anderson, “Local moments and localized states,” Reviews of Modern Physics 50, 191 (1978).
  * Splettstoesser  _et al._ (2005) J. Splettstoesser, M. Governale, J. König, and R. Fazio, “Adiabatic pumping through interacting quantum dots,” Phys. Rev. Lett. 95, 246803 (2005).
  * Arrachea and Moskalets (2006) L. Arrachea and M. Moskalets, “Relation between scattering-matrix and Keldysh formalisms for quantum transport driven by time-periodic fields,” Phys. Rev. B 74, 245322 (2006).
  * Li  _et al._ (2019) H. Li, L. J. Fernández-Alcázar, F. Ellis, B. Shapiro, and T. Kottos, “Adiabatic thermal radiation pumps for thermal photonics,” Phys. Rev. Lett. 123, 165901 (2019).
  * Note (1) In our model, the equilibrium position (no excitation of vibrational modes for the oscillation of GNRs) corresponds to a flat still membrane, equivalent to setting $\bm{q}_{0}=\bm{0}$.
  * Pereira  _et al._ (2009) V. M. Pereira, A. H. Castro Neto, and N. M. R. Peres, “Tight-binding approach to uniaxial strain in graphene,” Physical Review B 80, 045401 (2009).
  * Groth  _et al._ (2014) C. W. Groth, M. Wimmer, A. R. Akhmerov, and X. Waintal, “Kwant: a software package for quantum transport,” New Journal of Physics 16, 063065 (2014).
  * Deghi  _et al._ (2021) S. E. Deghi, L. J. Fernández-Alcázar, H. M. Pastawski, and R. A. Bustos-Marún, “Current-induced forces in single-resonance systems,” Journal of Physics: Condensed Matter 33, 175303 (2021).
  * Pastawski and Medina (2001) H. M. Pastawski and E. Medina, “‘Tight Binding’ methods in quantum transport through molecules and small devices: From the coherent to the decoherent description,” Revista Mexicana de Física 47S1, 1 (2001).
  * Note (2) Take into account that for Fermi energies $\varepsilon_{\text{F}}$ too close to the Dirac point ($|\varepsilon_{\text{F}}|\lesssim k_{\text{B}}T$) the calculation of pumping currents at finite temperatures would present numerical issues due to the divergence of emissivities there, see Eq. (9).
  * Bolotin  _et al._ (2008a) K. Bolotin, K. Sikes, Z. Jiang, M. Klima, G. Fudenberg, J. Hone, P. Kim, and H. Stormer, “Ultrahigh electron mobility in suspended graphene,” Solid State Communications 146, 351 (2008a).
  * Bolotin  _et al._ (2008b) K. I. Bolotin, K. J. Sikes, J. Hone, H. L. Stormer, and P. Kim, “Temperature-dependent transport in suspended graphene,” Physical Review Letters 101, 096802 (2008b).
  * Singh  _et al._ (2010) V. Singh, S. Sengupta, H. S. Solanki, R. Dhall, A. Allain, S. Dhara, P. Pant, and M. M. Deshmukh, “Probing thermal expansion of graphene and modal dispersion at low-temperature using graphene nanoelectromechanical systems resonators,” Nanotechnology 21, 165204 (2010).
  * Ma  _et al._ (2013) L. Ma, J. Wang, and F. Ding, “Recent progress and challenges in graphene nanoribbon synthesis,” ChemPhysChem 14, 47 (2013).
  * Zhang  _et al._ (2020) Z. Zhang, X. Qi, J. Zhang, C. Guo, and Z. Zhu, “Graphene-enabled electrically tunability of metalens in the terahertz range,” Optics Express 28, 28101 (2020).
  * Higuchi  _et al._ (2017) T. Higuchi, C. Heide, K. Ullmann, H. B. Weber, and P. Hommelhoff, “Light-field-driven currents in graphene,” Nature 550, 224 (2017).
  * Cong  _et al._ (2019) X. Cong, Q.-Q. Li, X. Zhang, M.-L. Lin, J.-B. Wu, X.-L. Liu, P. Venezuela, and P.-H. Tan, “Probing the acoustic phonon dispersion and sound velocity of graphene by Raman spectroscopy,” Carbon 149, 19 (2019).
  * Note (3) The electronic Hamiltonian $H$ varies when the GNR oscillates with a given normal mode (as explained in Sec. II.4). Let us recall that the movement of the GNR is considered classically here, so $\bm{q}$ (the vector containing the amplitude of the normal modes) is a parameter of $H$. In our work, we are considering small amplitudes of the normal modes, so $\bm{q}(t)\approx\bm{q}_{0}$. However, $\bm{\dot{q}}$ will still change stochastically due to thermal fluctuations. Therefore, when calculating the current noise one has to average over the stochastic variation of $\bm{\dot{q}}$, this is $\Braket{...}_{\text{c}}$, which is different from $\Braket{...}$, the expectation value of some quantity.
  * Polianski  _et al._ (2002) M. L. Polianski, M. G. Vavilov, and P. W. Brouwer, “Noise through quantum pumps,” Phys. Rev. B 65, 245314 (2002).
  * Riwar  _et al._ (2013) R.-P. Riwar, J. Splettstoesser, and J. König, “Zero-frequency noise in adiabatically driven interacting quantum systems,” Physical Review B 87, 195407 (2013).
  * Balakrishnan (2021) V. Balakrishnan,  _Elements of Nonequilibrium Statistical Mechanics_ (Springer, Cham, 2021).
  * Note (4) This limit means that correlations decay instantaneously as compared with the measurement time. This would lead $S(\tau)$ to be proportional to the Dirac delta function (with a proportionality factor equal to the so called zero-frequency noise), which is an usual assumption in noise treatment.
  * Blanter and Büttiker (2000) Y. Blanter and M. Büttiker, “Shot noise in mesoscopic conductors,” Physics Reports 336, 1 (2000).
  * Cattena  _et al._ (2010) C. J. Cattena, R. A. Bustos-Marún, and H. M. Pastawski, “Crucial role of decoherence for electronic transport in molecular wires: Polyaniline as a case study,” Phys. Rev. B 82, 144201 (2010).
  * Fernández-Alcázar  _et al._ (2019) L. J. Fernández-Alcázar, H. M. Pastawski, and R. A. Bustos-Marún, “Nonequilibrium current-induced forces caused by quantum localization: Anderson adiabatic quantum motors,” Phys. Rev. B 99, 045403 (2019).
  * Takamura  _et al._ (2016) M. Takamura, H. Okamoto, K. Furukawa, H. Yamaguchi, and H. Hibino, “Energy dissipation in graphene mechanical resonators with and without free edges,” Micromachines 7 (2016), 10.3390/mi7090158.
  * Splettstoesser  _et al._ (2006) J. Splettstoesser, M. Governale, J. König, and R. Fazio, “Adiabatic pumping through a quantum dot with Coulomb interactions: A perturbation expansion in the tunnel coupling,” Physical Review B 74, 085305 (2006).
  * Ribetto  _et al._ (2021) F. D. Ribetto, R. A. Bustos-Marún, and H. L. Calvo, “Role of coherence in quantum-dot-based nanomachines within the Coulomb blockade regime,” Physical Review B 103, 155435 (2021).
  * Fisher and Lee (1981) D. S. Fisher and P. A. Lee, “Relation between conductivity and transmission matrix,” Physical Review B 23, 6851 (1981).
  * Bode  _et al._ (2012) N. Bode, S. Viola-Kusminskiy, R. Egger, and F. von Oppen, “Current-induced forces in mesoscopic systems: A scattering-matrix approach,” Beilstein Journal of Nanotechnology 3, 144 (2012).