# Coulomb drag by motion of a monolayer polar crystal through graphene nano-
constriction

A. L. Chudnovskiy  I. Institut für Theoretische Physik, Universität Hamburg,
Notkestraße 9, D-22607 Hamburg, Germany

###### Abstract

We predict theoretically that the motion of a polar crystalline layer between
two graphene planes exerts Coulomb drag on electrons in graphene, inducing a
DC-drag current. The physical mechanism of this drag relies on the intervalley
scattering of charge carriers in graphene by the time-dependent potential of
the moving crystalline layer. The drag current exhibits a nonlinear Hall
effect, which can be used for the experimental measurement of the Berry
curvature. In turn, the drag of electrons in graphene exerts a back-action on
the crystalline layer, which can be expressed in terms of increased dynamic
viscosity.

Recent experiments have shown that confining water in nano-constrictions leads
to an ordering phase transition, resulting in a lattice of water molecules and
forming a two-dimensional ice crystal [1, 2, 3, 4]. Due to the inherent dipole
moment of a water molecule, this single-layer ice forms a polar crystal with a
frozen pattern of dipoles. Under applied external pressure, the ice layer can
move through constrictions. A recent example of such motion is demonstrated in
experiments in Refs. [5, 6], where water flowed through a nanoconstriction of
a single graphene layer’s width between graphene planes. The motion of the
spatially ordered dipole moments of the ice layer generates a time- and space-
periodic electric field acting on the electrons in the walls of the
constriction. Notably, the spatial period of the electric field set by the
lattice spacing of the polar crystal is comparable to the nearest neighbor
distance in graphene. For instance, for the crystal of square ice, the lattice
spacing is given by $a_{0}\approx 2.8$ Å [1], which is comparable to the
nearest neighbor distance in graphene $a=1.42$ Å [7]. Consequently, the
momentum transferred from the moving crystal to the electrons in graphene is
close to the distance between the $K$ and $K^{\prime}$ points in reciprocal
space. This opens the possibility for Coulomb drag of electrons in graphene
induced by the motion of the polar crystalline layer, facilitated by
intervalley and/or umklapp scattering processes. The fundamental physical
mechanism for this drag involves intervalley transitions of electrons induced
by the external time-dependent potential created by the moving crystalline
layer.

In this paper, we investigate the DC-drag current density in graphene induced
by the motion of a crystalline layer. We uncover several distinct features of
the DC-drag current that set it apart from conventional drag phenomena: (i)
The drag current is proportional to the velocity of the crystalline layer.
(ii) The direction of the drag current is determined by the orientation of the
moving crystalline lattice with respect to the vector connecting the two
K-points in graphene, rather than the direction of the crystalline layer’s
velocity. (iii) The drag exhibits signatures of the nonlinear Hall effect,
which allows its use for experimental measurement of the Berry curvature. (iv)
The drag sets in at a finite doping level that creates Fermi surfaces around
the K points. (v) Unlike conventional Coulomb or phonon drag, the leading
contribution to the drag current is temperature-independent.

The scattering mechanism underlying the drag effect is similar to the umklapp
scattering that causes excess resistivity in graphene/hexagonal boron nitride
heterostructures [8]. However, unlike the stationary superlattice in
heterostructures, scattering by the moving crystalline layer transfers energy
to electrons in graphene, thereby inducing a finite drag current. The drag
effect considered here also bears some analogy to the acousto-electric effect,
where electric current is generated by acoustic waves [9, 10, 11, 12, 13, 14].
The main difference lies in the very small wavelength associated with the
potential of the moving polar crystal, which can induce intervalley
transitions or umklapp processes. Although motivated by recent experiments on
the ice layer sliding through the mamo-constriction, the drag phenomenon is of
a quite general nature and is applicable to any two-dimensional (2D) polar
crystal layer (see Fig. 1a).

Let us introduce the model for the Coulomb drag by a moving polar crystalline
layer. The crystalline layer exerts a periodic potential on electrons in the
immediate vicinity determined by its crystal structure, given by $U({\bf
r})\sim\sum_{i=1,2}\cos\left({\bf Q}_{i}{\bf r}\right)$, where ${\bf Q}_{1,2}$
denote the reciprocal lattice vectors of the crystal. The potential of the
moving crystal gains time dependence through a boost transformation ${\bf
r}\rightarrow{\bf r}-{\bf v}_{0}t$, where ${\bf v}_{0}$ denotes the velocity
of the crystalline layer. Therefore, we model the periodic time-dependent
potential exerted by the layer moving with velocity $v_{0}$ by the expression
$U({\bf r,t})=\sum_{i=1,2}U_{i}\cos({\bf Q}_{i}({\bf r}-{\bf v}_{0}t)).$ The
scalar products of the vectors ${\bf Q}_{1,2}$ with the velocity of the
crystalline layer determine the frequencies of the time dependent potential
$\omega_{i}=({\bf Q}_{i}\cdot{\bf v}_{0})$. As will be shown below, the
amplitude of the drag current is proportional to the frequency of the driving
potential, hence the component of the potential with the higher frequency
provides the dominant contribution to the drag. For instance, among several
possible orderings of water molecules, the simplest one is a square ice, in
which molecules are ordered in a one-layer square lattice. In this
configuration the vectors ${\bf Q}_{1}$ and ${\bf Q}_{2}$ are orthogonal. The
maximal drag is achieved when one reciprocal wave vector, ${\bf Q}_{1}$, is
aligned with the velocity of the layer ${\bf v}_{0}$. In this case, the other
wave vector ${\bf Q}_{2}$ is perpendicular to the velocity, and its component
provides no contribution to the drag. Therefore, in what follows we consider
only one component of the periodic potential and drop the subscript at the
wave vector ${\bf Q}$.

The external potential of the crystalline layer affects both the on-site
energies and the nearest-neighbor hopping amplitudes in the graphene lattice.
In this paper, we focus on the effect on the on-site energies. Modulation of
the nearest-neighbor hopping does not lead to qualitatively new results, so we
relegate its discussion to the Supplemental Material [15]. We describe the
modulation of the on-site energies by the Hamiltonian

$\displaystyle H_{U}=U_{0}\sum_{n,m}\\!\left\\{\cos({\bf Q}{\bf
r}_{nm}\\!-\\!\omega t)\bar{\psi}_{A}({\bf r}_{nm})\psi_{A}({\bf
r}_{nm})+\right.$ $\displaystyle\left.\cos({\bf Q}({\bf
r}_{nm}\\!+{\bm{\delta}}_{3}\\!-\\!\omega t)\bar{\psi}_{B}({\bf
r}_{nm}\\!+{\bm{\delta}}_{3})\psi_{B}({\bf
r}_{nm}\\!+{\bm{\delta}}_{3})\right\\}.$ (1)

Here $A$ and $B$ denote the two sublattices of the hexagonal graphene lattice
[16, 7]. The atoms of the $A$-sublattice are situated at ${\bf r}_{nm}=n{\bf
a}_{1}+m{\bf a}_{2}$, where $n,m$ are integers and ${\bf a}_{1,2}$ denote the
lattice vectors in graphene. The positions of atoms of the $B$-sublattice are
given by ${\bf r}_{nm}+{\bm{\delta}}_{3}$ (see Fig. 1b).

Consider now the action of external potential Eqs. (1) as a source of the
intervalley and umklapp scattering. An elemetary scattering process requires
conservation of energy and quasi-momentum (with $\hbar=1$ throughout this
paper)

$\epsilon({\bf p})=\epsilon({\bf p^{\prime}})+\omega,\,\,\,{\bf
p^{\prime}}+{\bf Q}={\bf p}+n_{1}{\bf b}_{1}+n_{2}{\bf b}_{2}+m({\bf K}-{\bf
K^{\prime}}).$ (2)

Here $n_{1}$, $n_{2}$, and $m$ are integers, ${\bf b}_{1}$ and ${\bf b}_{2}$
denote the reciprocal lattice vectors in graphene, and ${\bf K}$ and ${\bf
K^{\prime}}$ denote the wave vectors at $K$-points. The wave vectors ${\bf
k^{\prime}}$ and ${\bf k}$ relate to the initial and the final states by the
scattering event. The relation of wave-vectors for a scattering event with
$n_{1}=n_{2}=0$, $m=1$ is shown in Fig. 1b). Using the data on the graphene
lattice structure one can conclude that the minimal mismatch of the vectors
${\bf k}$ and ${\bf k^{\prime}}$ is reached by the intervalley scattering,
when the reciprocal lattice vector ${\bf Q}$ of the crystalline layer is
parallel to the vector ${\bf K}-{\bf K^{\prime}}$ connecting the different
$K$-points. Taking the length of the vector ${\bf Q}$ proper to the square ice
lattice, $|{\bf Q}|=2\pi/a_{0}$, using the standard values for the $K$ points

${\bf K}=\frac{2\pi}{3a}(1,1/\sqrt{3}),\,\,\,{\bf
K^{\prime}}=\frac{2\pi}{3a}(1,-1/\sqrt{3}),$ (3)

and considering the geometry ${\bf Q}\parallel{\bf K}-{\bf
K^{\prime}}\parallel{\bf\hat{y}}$, we estimate the difference of the final and
initial wave-vectors as

${\bf q}={\bf p}-{\bf p^{\prime}}={\bf Q}-({\bf K}-{\bf
K^{\prime}})\approx\frac{0.77}{a}{\bf\hat{y}}.$ (4)

Absorption of the momentum ${\bf q}$ by the electron system of graphene
results in the drag current. Furthermore, energy conservation by scattering
requires that the vector ${\bf q}$ connects two states in the graphene energy
band, that is $\epsilon({\bf p}+{\bf q})=\epsilon({\bf p})+\omega$. The
maximal possible value of the frequency $\omega$ is given by
$\omega=v_{0}|Q|$. However, since the drift velocity of the ice layer
$v_{0}\sim 1$ m/s is much smaller that the electron velocity in graphene
$v_{F}\sim 10^{6}$ m/s, the energy of the transition between the states with
${\bf p}$ parallel to ${\bf p}+{\bf q}$, given by $v_{F}|q|$, is much larger
than $\omega$. The energy conservation condition can be satisfied though, if
the wave-vectors of the initial and final states are anti-parallel, lying at
opposite points of the Fermi-surface, as shown in Fig. 1c). Then the energy
conservation reads $v_{F}(|q|-2p_{F})=v_{0}|Q|$, and it can be satisfied at
the Fermi wave vector $p_{F}=\frac{1}{2}(|q|-\frac{v_{0}}{v_{F}}|Q|)$. This
Fermi wave vector marks the threshold for the Fermi energy $\mu=v_{F}p_{F}$,
at which the drag current starts. For larger Fermi wave vectors, the
conservation of energy and quasi-momentum can also be satisfied at some angle
beween the vectors ${\bf p}$ and ${\bf q}$, as shown in Fig. 1b). Therefore,
the intervalley scattering provides the most efficient drag mechanism under
the conditions where the velocity and one crystallographic axis of the moving
crystalline layer are parallel to the ${\bf K}-{\bf K^{\prime}}$ vector in
graphene.

Reaching the threshold chemical potential for the drag requires substantial
doping, which nevertheless lies within the reach of modern experiments. Using
the values $v_{F}\sim 10^{6}$m/s, $v_{0}\sim 1$ m/s, and $|{\bf q}|\sim|{\bf
K}-{\bf K^{\prime}}-{\bf Q}|\sim 5.5\cdot 10^{9}$m-1 (see Eq. 4), we estimate
the threshold chemical potential as $\mu\sim 0.6t\sim 2$ eV. It is noteworthy
that the threshold chemical potential is still below the van Hove singularity
at the M-point of graphene $E_{\mathrm{vH}}\sim 3eV$, so considering separate
$K$ points remains valid. In light of recent progress in doping graphene up to
and even above the M-point [17, 18], reaching the threshold doping for the
drag looks quite feasible. We estimate the electron density around the
threshold as $n_{e}\sim 10^{18}$ m-2.

Now let us proceed with calculation of the induced drag current. For the
states close to the $K$ and $K^{\prime}$ points, the electron spectrum in the
graphene layer given by a well-known Hamiltonian [16, 7]

$H_{g}=-iv_{F}\\!\int d^{2}{\bf r}\left[\Psi_{1}^{+}({\bf
r}){\bm{\sigma}}\cdot\\!\nabla\Psi_{1}({\bf r})+\Psi_{2}^{+}({\bf
r}){\bm{\sigma}^{*}}\cdot\\!\nabla\Psi_{2}({\bf r})\right],$ (5)

where the indexes $1,2$ correspond to the $K$ and $K^{\prime}$ points
respectively, and $\Psi({\bf k})=(\psi_{A}({\bf k}),\psi_{B}({\bf k}))^{T}$
denotes a spinor in the pseudospin space. Because the electron spin does not
play a role in the drag physics considered here, spin indexes are suppressed.
The interaction with the moving crystalline layer is given by Eqs. (1). Upon
projecting Eq. (1) on the states close to the $K$ points, the drag potential
acquires the form

$\displaystyle H_{d}=\frac{U_{0}}{2}\int_{\bf p}\left\\{\Psi^{+}_{1}({\bf
p}+{\bf q})\left({\bf u}({\bf q},{\bf
p})\cdot{\bm{\sigma}}\right)\Psi_{2}({\bf p})e^{-i\omega t}+\right.$
$\displaystyle\left.\Psi^{+}_{2}({\bf p})\left({\bf u}^{*}({\bf q},{\bf
p})\cdot{\bm{\sigma}}\right)\Psi_{1}({\bf p}+{\bf q})e^{i\omega t}\right\\},$
(6)

where $\left({\bf u}({\bf q},{\bf
p})\cdot{\bm{\sigma}}\right)=\sum_{i=0,3}u_{i}({\bf q},{\bf p})\sigma_{i}$,
and we set $U_{0}$ to be the energy scale characterizing the interaction
strength. The interaction constants $u_{0,3}$ account for the shift of local
on-site energies. We note in passing that modulation of the hopping amplitude
introduces the terms $\sum_{i=1,2}u_{i}({\bf q},{\bf p})\sigma_{i}$ with the
constants $u_{1,2}$. Explicit form of the dimensionless coupling constants
$u_{i}({\bf q},{\bf p})$ is provided in the Supplemental Material [15]. In
writing Eq. (55) we retained only the time-dependent exponentials that can
satisfy the energy conservation condition $v|{\bf p}|\pm\omega=v|{\bf
p}\pm{\bf q}|$ according to the dispersion in graphene close to the
$K$-points.

Figure 1:  a) Scheme of motion of a crystalline layer through a graphene nano-
constriction. b) Graphene and square ice lattices. The maximal energy transfer
$\omega=({\bf v}_{0}\cdot{\bf Q})$ is achieved when the velocity of the ice
${\bf v}_{0}$ is parallel to one of the main lattice directions (see the
text). c) Reciprocal lattice of graphene. The vectors $\bf b_{1}$, $\bf b_{2}$
indicate the first Brilluin zone. Filled circles denote the Fermi sea arond
the $K$ and $K^{\prime}$ points. Arrows indicate the momentum conservation by
intervalley transitions leading to the drag current ${\bf K}-{\bf
K^{\prime}}+{\bf q}={\bf Q}$. d) Definitions of the angles $\alpha$ and $\phi$
and their relation to the directions of the initial (${\bf p}$) and the final
(${\bf p^{\prime}}$) momenta by a single scattering event. At the threshold
doping, the vector ${\bf Q}$ is parallel to ${\bf K}-{\bf K^{\prime}}$, the
transferred quasi-momentum $\hbar q=2p_{F}$, the angle $\phi=\pm\pi$.

We calculate the drag current in frame of perturbative expansion in the
interaction strength $U_{0}$ using the Keldysh formalism. The unperturbed
Green functions for the states near the K-points are determined by the
Hamiltonian Eq. (5). They can be conveniently represented in the form that
separates a valley independent pole structure and the valey dependent
pseudospin structure $\hat{G}_{\nu}^{R/A}(\epsilon,{\bf
p})=G^{R/A}(\epsilon,{\bf p})\hat{g}_{\nu}(\epsilon,{\bf p}),$ where
$G^{R/A}(\epsilon,{\bf
p})=\epsilon/\left[(\epsilon\pm\frac{i}{2\tau})^{2}-v_{F}^{2}p^{2}\right]$
denote the pseudospin-independent parts of the Green functions, and
$\hat{g}_{K/K^{\prime}}(\epsilon,{\bf
p})=\sigma_{0}+\frac{v_{F}}{\epsilon}(p_{x}\sigma_{x}\pm p_{y}\sigma_{y})$.
Here the Pauli-matices $\sigma_{i}$ act on pseudospin. The Keldysh conponent
of the Green functions is given by $\hat{G}_{\nu}^{K}(\epsilon,{\bf
p})=\tanh\left(\frac{\epsilon}{2T}\right)\left(\hat{G}_{\nu}^{R}(\epsilon,{\bf
p})-\hat{G}_{\nu}^{A}(\epsilon,{\bf p})\right)$. The lowest order
contributions to the drag current density are given by the diagrams in Fig.
(2), where solid lines denote Green functions for the states close to $K$ and
$K^{\prime}$ points. Detailed calculations of AC- and DC- drag current
densities are provided in the Supplemental Material [15].

The AC drag current at the frequency $\omega$ is given by the diagrams of the
first order in the interaction strength in Fig. 2a). Here the two diagrams
correspond to the current induced in the $K$ and $K^{\prime}$ valleys by
intervalley scattering. However, we believe that the AC-drag will be strongly
suppressed in the realistic experiment because of defects in the ice lattice
structure. Those defects would lead to phase slips in the driving potential
and eventually to a randomization of the phase, which in turn will suppress
the AC drag current.

In contrast, the DC current is largely insensitive to small random changes in
the phase of the driving potential.

Figure 2: Diagrams for the ac (panel a) and dc (panel b) drag current
densities in the lowest order in interactions.

The total DC current density is calculated as a sum of contributions from the
two valleys, ${\bf j}_{\mathrm{dc}}=\left\langle{\bf
j}_{KK}\right\rangle+\left\langle{\bf j}_{K^{\prime}K^{\prime}}\right\rangle$,
each represented by a diagram in Fig. 2b). These diagrams depict the
interaction-induced rectification, where the AC displacement current of the
moving layer is converted into a DC drag current in graphene. Since the drag
requires a finite chemical potential $\mu$ in graphene, the temperature
effects are of minor importance for $T\ll\mu$, which is valid up to room
temperature. In this case, the calculation of the diagrams in Fig. 2b) can be
performed by setting $T=0$. Furthermore, for $\mu>0$ the states in conduction
band provide the dominant contribution to the drag current. Upon projection on
the states of conduction band, we obtain the following expression for the DC-
drag current density

${\bf
j}_{\mathrm{DC}}=\frac{\pi^{2}e}{\hbar^{3}}\frac{U_{0}^{2}p_{F}\omega\tau}{2v_{F}|q|}\frac{\mathcal{\bf
F}^{\mathrm{DC}}(\alpha,\phi)}{|\sin\phi|}=\frac{\pi^{4}}{e\hbar}\frac{v_{0}}{a_{0}}\frac{U_{0}^{2}}{\rho_{g}v_{F}|q|}\frac{\mathcal{\bf
F}^{\mathrm{DC}}(\alpha,\phi)}{|\sin\phi|}.$ (7)

In the second equation, we express the mean free time $\tau$ throught the
resistivity of graphene $\rho_{g}=\frac{\pi\hbar}{e^{2}\mu\tau}$, where
$\mu=p_{F}v_{F}$ denotes the chemical potential. The DC-drag is purely
dissipative, with the drag current beingproportional to the mean-free time
$\tau$ of electron states in graphene. Additionally, we use the relation
between the frequency and the velocity of the crystalline layer $\omega=2\pi
v_{0}/a_{0}$, assuming the velocity to be parallel to the vector ${\bf K}-{\bf
K^{\prime}}$ connecting the K-points in graphene.

The dependencies of the drag current on chemical potential for two different
directions of the transfferred quasi-momentum ${\bf q}$ are shown in Fig. 3a).
There is a sharp threshold chemical potential for the drag effect. Finite
temperature leads to a smearing of the step at the threshold
$\mu=\left(\frac{v_{F}|q|}{\omega}-1\right)/2$, but it does affect the
qualitative dependence of the current on chemical potenital. Above this
threshold, the direction of the current is determined by the vector ${\bf
F}^{\mathrm{DC}}(\alpha,\phi)$, which is given by the product of matrix
elements of the current density operator and interaction vertices ${\bf
F}^{\mathrm{DC}}(\alpha,\phi)=\Re[{\bf j}_{K}({\bf p})e^{i\phi_{\bf p}}-{\bf
j}_{K^{\prime}}({\bf p})e^{i\phi_{\bf p^{\prime}}}]\left|\langle{\bf
p+q},K^{\prime}|({\bf u}\cdot\bm{\sigma})|{\bf p},K\rangle\right|^{2}$, where
$\phi_{\bf p}=\alpha\pm\phi$ denotes the direction of the vector ${\bf p}$,
and $\phi_{\bf p^{\prime}}=\alpha\mp\phi$ denotes the direction of the vector
${\bf p^{\prime}}$, as shown in Fig. 1d). The angle $\alpha$ specifies the
direction of the transferred quasi-momentum, and one might naively expect the
drag current to align with this direction. This expectation holds true only
when the reciprocal vector of the moving lattice is parallel to the vector
connecting the K points in graphene, ${\bf Q}\parallel{\bf K^{\prime}}-{\bf
K}$. In the general case, the direction of the drag current deviates from the
transferred quasi-momentum, as shown in Fig. 3b). This deviation is a
signature of the nonlinear valley Hall effect. The nonlinear valley Hall
effect bears an analogy with a nonlinear Hall effect considered previously in
Refs. [19, 20]. It also suggests using the drag as a momentun-resolved probe
of the Berry curvature $\Omega({\bf p},{\bf q})$. Specifically, the Berry
curvature can be calculated from the relation of the longitudinal (with
respect to transferred quasimomentum ${\bf q}$ ) $j_{q}$ and Hall, $j_{H}$
components of the drag current

$\Omega({\bf p},{\bf
q})\\!=\\!\frac{\pi^{2}j_{H}/j_{q}}{v_{F}q^{2}\cos\phi}\Re\left\\{j_{q,K}({\bf
p})e^{i\phi_{\bf p}}+j_{q,K^{\prime}}({\bf p+q})e^{i\phi_{\bf p+q}}\right\\}$
(8)

Here the components of the drag current $j_{H}$, $j_{q}$ are inferred from
experimental measurements, whereas the matrix elements $j_{q,K}({\bf p})$,
$j_{q,K^{\prime}}({\bf p+q})$ are known from theoretical caluclations.
Detailed derivation fo Eq (152) is provided in the Supplemental Material [15].
Fig. 3c) shows the Berry phase as a function of chemical potential. Changing
the chemical potential, one probes the Berry curvature at different states at
the Fermi.surface, as shown at Fig. 1d).

The drag current given by Eq. (7) formally diverges at the threshold
$p_{F}=\frac{1}{2}(q-\omega/v_{F})$, corresponding to the angle $\phi=\pi$, as
illustrated in the left panel in Fig. 3. That divergence is removed by
accounting for scattering processes that change the momentum of
quasiparticles. Considering the finite lifetime $\tau$ leads to the
replacement $\sin\phi\rightarrow\frac{1}{\sqrt{v_{F}q\tau}}$ in Eq. (7) in the
region $v_{F}|2p_{F}-q|<1/\tau$ close to the threshold. The vector-function
$\mathcal{\bf F}^{\mathrm{DC}}(\alpha,\phi)$ remains nonsingular at the
threshold value $\phi=\pi$.

Figure 3: a) DC drag current as function of chemical potential for two
different directions of the transferred quasimomentum. The threshold chemical
potential $\mu_{0}$ corresponds to the transferred quasi-momentum $q=2p_{F}$;
b) The angle between the direction of the transferred quasimomentum ${\bf
q}={\bf Q}-({\bf K^{\prime}}-{\bf K})$ and the drag current ${\bf
j}_{\mathrm{DC}}$ as function of the angle between ${\bf q}$ and ${\bf
K^{\prime}}-{\bf K}$. In the absence of the Hall effect, the drag current
would be parallel to ${\bf q}$; c) Berry curvature as function of chemical
potential.

Essentially, the drag effect considered here represents the conversion of the
kinetic energy of the moving polar crystalline layer into the energy of a DC
electric current in graphene. This energy transfer from the moving crystalline
layer to the electrons in graphene results in the back-action force on the
layer, which can be formulated in terms of a renormalized shear viscosity
coefficient. If one adopts the Poiseuille equation for a crystalline layer
moving with a velocity $v$ between the two graphene layers, the following
relation is obtained $\Delta phL_{\perp}=8\eta vL_{\perp}L/h$, where $\Delta
p$ denotes the pressure difference at the ends of the interlayer channel of
the lengh $L$, $h$ denotes the distance between the graphene layers forming
the channel (see Fig. 1a), $v$ is the velocity of the crystalline layer, and
$\nu_{0}$ denotes the shear viscosity in the absence of the drag current. To
determine the change in shear viscosity due to the drag, consider the energy
dissipation during the motion of the crystalline layer. In the absence of the
drag current, the energy dissipation rate is given by the relation
$\dot{Q}_{0}=\Delta phL_{\perp}v=8\eta_{0}v^{2}L_{\perp}L/h.$ The drag current
results in the additional dissipation, given by
$\dot{Q}_{\mathrm{drag}}=2j_{d}^{2}L_{\perp}^{2}\rho_{g}L=f_{d}^{2}v^{2}L_{\perp}^{2}\rho_{g}L,$
where $j_{d}=vf_{d}$ is the absolute value of the drag current density in
graphene, $\rho_{g}$ denotes the resititivity of the graphene layer, and the
factor 2 accounts for the two layers forming the nano-channel. In the last
equation, we took into accont the proportionality of the drag current to the
velocity of the moving layer, as given by Eq. (7). The coefficient $f_{d}$ can
be read from Eq. (7),
$f_{d}=\frac{\pi^{4}}{e\hbar}\frac{U_{0}^{2}}{a_{0}\rho_{g}v_{F}|q|}\frac{|\mathcal{\bf
F}^{\mathrm{DC}}(\alpha,\phi)|}{|\sin\phi|}$. Therefore, the total energy
dissipation in presence of the drag reads

$\dot{Q}_{0}+\dot{Q}_{\mathrm{drag}}=\eta_{0}v^{2}\frac{8L_{\perp}L}{h}+f_{d}^{2}v^{2}L_{\perp}^{2}\rho_{g}L\equiv\eta
v^{2}\frac{8L_{\perp}L}{h},$ (9)

where the last equation defines the renormalized shear friction coefficient
$\eta$. Eq. (9) allows one to determine the shear friction coefficient due to
the drag effect, $\eta=\eta_{0}+\eta_{d}$, where
$\eta_{d}=f_{d}^{2}\rho_{g}L_{\perp}h/4$. The change in dynamic viscosity
provides an alternative way to detect the drag phenomenon. Specifically,
crossing the drag threshold by altering the chemical potential of graphene
would lead to a sharp decrease in the velocity of the crystalline layer, as
described by the relation
$\frac{v}{v_{0}}=\frac{\eta_{0}}{\eta_{0}+\eta_{d}}.$

In summary, we predict a Coulomb drag effect that should be observable in
experiments involving sliding polar crystalline layers through nano-
constrictions between graphene planes. This drag mechanism hinges on the
comparable lattice constants of graphene and the moving crystalline layer,
which facilitates intervalley and umklapp scattering induced by the potential
of the moving layer. The predicted effect can be experimentally detected
through the induction of drag current or drag voltage in graphene planes, as
well as by observing an increase in viscosity due to the motion of the
crystalline layer, resulting in a reduction of its velocity. Moreover, the
drag manisfests a nonlinear intervalley Hall effect – specifically, a
component of the drag current perpendicular to the transferred quasi-momentum.
This effect can be used for a momentum-resolved determination of the Berry
curvature. Achieving the drag effect requires significant doping of graphene
to create sufficiently large Fermi surfaces around each $K$-point. Our
estimations indicate that modern experimental doping levels are adequate for
observing this phenomenon.

The author acknowledges fruitful discussions with M. Trushin, A. Kamenev, A.
Levchenko, V. Cheianov, and T. Wehling.

## References

  * Algara-Siller  _et al._ [2015] G. Algara-Siller, O. Lehtinen, F. C. Wang, R. R. Nair, U. Kaiser, H. A. Wu, A. K. Geim, and I. V. Grigorieva, Square ice in graphene nanocapillaries, Nature 519, 443 (2015).
  * Negi  _et al._ [2022] S. Negi, A. Carvalho, M. Trushin, and A. H. C. Neto, Edge-driven phase transitions in 2d ice, The Journal of Physical Chemistry C 126, 16006 (2022).
  * Gao  _et al._ [2018] Z. Gao, N. Giovambattista, and O. Sahin, Phase diagram of water confined by graphene, Scientific Reports 8, 6228 (2018).
  * Li and Schmidt [2019] S. Li and B. Schmidt, Two-dimensional water in graphene nanocapillaries simulated with different force fields: Rhombic versus square structures, proton ordering, and phase transitions (2019), arXiv:1901.04236 [physics.chem-ph] .
  * Gopinadhan  _et al._ [2019] K. Gopinadhan, S. Hu, A. Esfandiar, M. Lozada-Hidalgo, F. C. Wang, Q. Yang, A. V. Tyurnina, A. Keerthi, B. Radha, and A. K. Geim, Complete steric exclusion of ions and proton transport through confined monolayer water, Science 363, 145 (2019).
  * Nair  _et al._ [2012] R. R. Nair, H. A. Wu, P. N. Jayaram, I. V. Grigorieva, and A. K. Geim, Unimpeded permeation of water through helium-leak–tight graphene-based membranes, Science 335, 442 (2012).
  * Katsnelson [2012] M. I. Katsnelson,  _Graphene: Carbon in Two Dimensions_ (Cambridge University Press, 2012).
  * Wallbank  _et al._ [2019] J. R. Wallbank, R. Krishna Kumar, M. Holwill, Z. Wang, G. H. Auton, J. Birkbeck, A. Mishchenko, L. A. Ponomarenko, K. Watanabe, T. Taniguchi, K. S. Novoselov, I. L. Aleiner, A. K. Geim, and V. I. Fal’ko, Excess resistivity in graphene superlattices caused by umklapp electron–electron scattering, Nature Physics 15, 32 (2019).
  * Parmenter [1953] R. H. Parmenter, The acousto-electric effect, Phys. Rev. 89, 990 (1953).
  * Eckstein [1964] S. G. Eckstein, Acoustoelectric Effect, Journal of Applied Physics 35, 2702 (1964).
  * Many and Balberg [1969] A. Many and I. Balberg, The acoustoelectric effect, in  _Electronic Structures in Solids: Lectures presented at the Second Chania Conference, held in Chania, Crete, June 30–July 14, 1968_, edited by E. D. Haidemenakis (Springer US, Boston, MA, 1969) pp. 385–416.
  * Fal’ko  _et al._ [1993] V. I. Fal’ko, S. V. Meshkov, and S. V. Iordanskii, Acoustoelectric drag effect in the two-dimensional electron gas at strong magnetic field, Phys. Rev. B 47, 9910 (1993).
  * Kalameitsev  _et al._ [2019] A. V. Kalameitsev, V. M. Kovalev, and I. G. Savenko, Valley acoustoelectric effect, Phys. Rev. Lett. 122, 256801 (2019).
  * Lapa  _et al._ [2020] P. N. Lapa, G. Kassabian, F. Torres, P. Salev, M.-H. Lee, J. del Valle, and I. K. Schuller, Acoustoelectric drag current in vanadium oxide films, Journal of Applied Physics 128, 155104 (2020).
  * [15] A. L. Chudnovskiy, Supplemental material, Supplemental material .
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Rev. Mod. Phys. 81, 109 (2009).
  * Rosenzweig  _et al._ [2019] P. Rosenzweig, H. Karakachian, S. Link, K. Küster, and U. Starke, Tuning the doping level of graphene in the vicinity of the van hove singularity via ytterbium intercalation, Phys. Rev. B 100, 035445 (2019).
  * Rosenzweig  _et al._ [2020] P. Rosenzweig, H. Karakachian, D. Marchenko, K. Küster, and U. Starke, Overdoping graphene beyond the van hove singularity, Phys. Rev. Lett. 125, 176403 (2020).
  * Deyo  _et al._ [2009] E. Deyo, L. E. Golub, E. L. Ivchenko, and B. Spivak, Semiclassical theory of the photogalvanic effect in non-centrosymmetric systems (2009), arXiv:0904.1917 .
  * Sodemann and Fu [2015] I. Sodemann and L. Fu, Quantum nonlinear hall effect induced by berry curvature dipole in time-reversal invariant materials, Phys. Rev. Lett. 115, 216806 (2015).

##  I Coulomb drag by motion of a monolayer polar crystal through graphene
nano-constriction: Supplemental Material

In this Supplemental Material we provide detailed derivations of the formulas
presented in the main text of the paper.

##  II Important length scales

Let us list the important length scales in graphene that will be used for
further estimations. Lattice vectors in graphene:

${\bf a}_{1}=\frac{a}{2}(3,\sqrt{3}),\,\,\,{\bf
a}_{2}=\frac{a}{2}(3,-\sqrt{3})$ (10)

Vectors connecting the nearest neighbors atoms in graphene

${\bm{\delta}}_{1}=\frac{a}{2}(1,\sqrt{3}),\,\,\,{\bm{\delta}}_{2}=\frac{a}{2}(1,-\sqrt{3}),{\bm{\delta}}_{3}=-a(1,0).$
(11)

$a\approx 1.42$ Å. Corresponding reciprocal-lattice vectors

${\bf b}_{1}=\frac{2\pi}{3a}(1,\sqrt{3}),\,\,\,{\bf
b}_{2}=\frac{2\pi}{3a}(1,-\sqrt{3})$ (12)

The $K$-points

${\bf K^{\prime}}=\frac{2\pi}{3a}(1,1/\sqrt{3}),\,\,\,{\bf
K}=\frac{2\pi}{3a}(1,-1/\sqrt{3}).$ (13)

There are also equivalent $K$ and $K^{\prime}$ points given by the vectors

${\bf K^{\prime}}<=>{\bf K^{\prime}}-{\bf b}_{1}<=>{\bf K^{\prime}}-{\bf
b}_{1}-{\bf b}_{2},\,\,\,{\bf K}<=>{\bf K}-{\bf b}_{2}<=>{\bf K}-{\bf
b}_{1}-{\bf b}_{2}.$ (14)

The lattice constant of the square ice crystal $a_{0}\sim 2.8$ Å. Furthermore
we put $\hbar=e=1$ in all intermediate formulas.

##  III Operator of electric current density in the tight-binding model of
graphene

Since spin plays no role in furhter considerations, we suppress spin-indexes
here and in what follows. The tight binding Hamiltonian of graphene (here the
sublattices $A,B$ correspond to pseudospin) is given by

$H=it\sum_{i}\sum_{j=1}^{3}(\psi^{+}_{A}({\bf r}_{i})\psi_{B}({\bf
r}_{i}+{\bf\delta}_{j})-\psi_{B}^{+}({\bf r}_{i}+{\bf\delta}_{j})\psi_{A}({\bf
r}_{i}))=\sum_{\bf k}(\psi_{A}^{+}({\bf k}),\psi_{B}^{+}({\bf
k}))\left(\begin{array}[]{cc}0&it\sum_{\langle j\rangle}e^{i{\bf
k}{\bm{\delta}}_{j}}\\\ -it\sum_{\langle j\rangle}e^{-i{\bf
k}{\bm{\delta}}_{j}}&0\end{array}\right)\left(\begin{array}[]{c}\psi_{A}({\bf
k})\\\ \psi_{B}({\bf k})\end{array}\right).$ (15)

Here $\psi^{+}_{A,B}({\bf r}_{i}),\psi_{A,B}({\bf r}_{i})$ denote fermion
creation and annihilation operators at the site $i$ of $A$ and $B$ sublattices
respectively, and the sum over $j$ runs over the nearest neighbor cites.
External vector potential induces phase shifts on links $t_{ij}\rightarrow
t_{ij}e^{i\phi_{ij}}$, where $\phi_{ij}={\bf A}_{i}\cdot{\bm{\delta}}_{j}$
(the scalar product of the vector potential at the site $i$ and the vector to
the next-nearest neighbor $j$). In the (real-time) action, the coupling to the
external vector potential is described by the following term

$S_{A}=-\int_{t}\sum_{i}\sum_{j=1}^{3}(\bar{\psi}_{A}({\bf
r}_{i}),\bar{\psi}_{B}({\bf
r}_{i}+{\bm{\delta}}_{j}))\left(\begin{array}[]{cc}0&ite^{i\phi_{ij}}\\\
-ite^{-i\phi_{ij}}&0\end{array}\right)\left(\begin{array}[]{c}\psi_{A}({\bf
r}_{i})\\\ \psi_{B}({\bf r}_{i}+{\bm{\delta}}_{j})\end{array}\right).$ (16)

Transforming to the ${\bf k}$-space according to

$\psi_{A}({\bf r})=\int_{\bf k}\psi_{A}({\bf k})e^{i{\bf k}{\bf
r}},\,\,\,\psi_{B}({\bf r}+{\bm{\delta}}_{3})=\int_{\bf k}\psi_{B}({\bf
k})e^{i{\bf k}{\bf r}},$ (17)

we obtain

$S_{A}=-\int_{t}\sum_{i}\sum_{j=1}^{3}\int_{\bf
k,k^{\prime}}\left(\bar{\psi}_{A}({\bf k^{\prime}}),\bar{\psi}_{B}({\bf
k^{\prime}})e^{-i{\bf
k^{\prime}}({\bm{\delta}_{j}}-{\bm{\delta}_{3}})}\right)e^{-i{\bf
k^{\prime}}{\bf r}_{i}}\left(\begin{array}[]{cc}0&ite^{i\phi_{ij}}\\\
-ite^{-i\phi_{ij}}&0\end{array}\right)e^{i{\bf k}{\bf
r}_{i}}\left(\begin{array}[]{c}\psi_{A}({\bf k})\\\ \psi_{B}({\bf k})e^{i{\bf
k}(\bm{\delta}_{j}-\bm{\delta}_{3})}\end{array}\right).$ (18)

The expression for the operator of the current density is obtained by the
variation of the action with respect to the vector potential.

${\bf j}({\bf r}_{i})=\frac{\delta S}{\delta{\bf A}({\bf
r}_{i})}=\sum_{\langle j\rangle}{\bm{\delta}}_{j}\frac{\delta
S}{\delta\phi_{ij}}\bigg{|}_{\phi_{ij}=0}=\int_{{\bf k},{\bf
k^{\prime}}}\sum_{j=1}^{3}e^{i({\bf k}-{\bf k^{\prime}}){\bf
r}_{i}}\left(\bar{\psi}_{A}({\bf k^{\prime}}),\bar{\psi}_{B}({\bf
k^{\prime}})\right)\left(\begin{array}[]{cc}0&-{\bm{\delta}}_{j}te^{i{\bf
k}({\bm{\delta}}_{j}-{\bm{\delta}}_{3})}\\\ -{\bm{\delta}}_{j}te^{-i{\bf
k^{\prime}}({\bm{\delta}}_{j}-{\bm{\delta}}_{3})}&0\end{array}\right)\left(\begin{array}[]{c}\psi_{A}({\bf
k})\\\ \psi_{B}({\bf k})\end{array}\right).$ (19)

###  III.1 Expansion of the current density operator for the states close to
$K$ and $K^{\prime}$ points.

Assume ${\bf k}={\bf K}+{\bf p}$, ${\bf k^{\prime}}={\bf K^{\prime}}+{\bf
p^{\prime}}$, ($|{\bf p}|,|{\bf p^{\prime}}|\ll K,K^{\prime}$). Then the
current density operator can be represented as a sum of intra- and inter-
valley components

${\bf j}({\bf r}_{i})={\bf j}_{KK}({\bf r}_{i})+{\bf
j}_{K^{\prime}K^{\prime}}({\bf r}_{i})+{\bf j}_{KK^{\prime}}({\bf r}_{i})+{\bf
j}_{K^{\prime}K}({\bf r}_{i}),$ (20)

where for instance

${\bf j}_{K^{\prime}K}({\bf r})=-t\sum_{j=1}^{3}\int_{{\bf p},{\bf
p^{\prime}}}e^{i({\bf K}-{\bf K^{\prime}}+{\bf p}-{\bf p^{\prime}}){\bf
r}}\left(\bar{\psi}_{A}({{\bf K^{\prime}}+\bf
p^{\prime}}),\bar{\psi}_{B}({{\bf K^{\prime}}+\bf
p^{\prime}})\right)\left(\begin{array}[]{cc}0&{\bm{\delta}}_{j}e^{i({\bf
K}+{\bf p})({\bm{\delta}}_{j}-{\bm{\delta}}_{3})}\\\
{\bm{\delta}}_{j}e^{-i({\bf K^{\prime}}+{\bf
p^{\prime}})({\bm{\delta}}_{j}-{\bm{\delta}}_{3})}&0\end{array}\right)\left(\begin{array}[]{c}\psi_{A}({{\bf
K}+\bf p})\\\ \psi_{B}({{\bf K}+\bf p})\end{array}\right).$ (21)

Let us introduce the pseudospin spinors for the states close to the $K$-points
$\Psi_{K}({\bf p})=(\psi_{A}({\bf K+p}),\psi_{B}({\bf K+p}))$,
$\Psi_{K^{\prime}}({\bf p^{\prime}})=(\psi_{A}({\bf
K^{\prime}+p^{\prime}}),\psi_{B}({\bf K^{\prime}+p^{\prime}}))$. Performing
the Fourier transform and denoting ${\bf q}={\bf p^{\prime}}-{\bf p}$, we can
write the Fourier component at the vector ${\bf K}-{\bf K^{\prime}}-{\bf q}$
as

$\displaystyle{\bf j}_{K^{\prime}K}({\bf q})=\int_{\bf
p}\bar{\Psi}_{K^{\prime}}({\bf p+q})\left(\begin{array}[]{cc}0&-t\sum_{\langle
j\rangle}{\bm{\delta}}_{j}e^{i({\bf K}+{\bf
p})({\bm{\delta}}_{j}-{\bm{\delta}}_{3})}\\\ -t\sum_{\langle
j\rangle}{\bm{\delta}}_{j}e^{-i({\bf K^{\prime}}+{\bf p}+{\bf
q})({\bm{\delta}}_{j}-{\bm{\delta}}_{3})}&0\end{array}\right)\Psi_{K}({\bf
p})$ (24)

Let us introduce the operators of partial current density in each valley,
which are defined as follows

${\bf j}_{K}({\bf p})=-t\sum_{j=1}^{3}{\bm{\delta}}_{j}e^{i({\bf K}+{\bf
p})({\bm{\delta}}_{j}-{\bm{\delta}}_{3})},\,\,\,\,{\bf j}_{K^{\prime}}({\bf
p})=-t\sum_{\langle j\rangle}{\bm{\delta}}_{j}e^{i({\bf K^{\prime}}+{\bf
p})({\bm{\delta}}_{j}-{\bm{\delta}}_{3})}.$ (25)

Then the inter- and intra-valley components of the current density operator
can be cast to the form

${\bf j}_{ab}({\bf q})=\int_{\bf p}(\bar{\Psi}_{a}({\bf p+q}){\bf J}_{ab}({\bf
p},{\bf p+q})\Psi_{b}({\bf p}),$ (26)

where

${\bf J}_{ab}({\bf p+q},{\bf p})=\left(\begin{array}[]{cc}0&{\bf j}_{b}({\bf
p})\\\ {\bf j}^{*}_{a}({\bf p}+{\bf q})&0\end{array}\right).$ (27)

Here the subscripts $a,b$ relate to the valley, taking the values
$K,K^{\prime}$.

##  IV Green functions near $K$ points in the Keldysh formalism

The Hamiltonian of graphene for the states close to the $K,K^{\prime}$-points
(${\bf K}=\frac{2\pi}{3a}(1,-1/\sqrt{3})$, ${\bf
K^{\prime}}=\frac{2\pi}{3a}(1,1/\sqrt{3})$) reads

$\displaystyle H_{K}=v_{F}\int_{{\bf p}}\Psi_{K}^{+}({\bf
p})(\sigma_{x}p_{x}+\sigma_{y}p_{y})\Psi_{K}({\bf p}),$ (28) $\displaystyle
H_{K^{\prime}}=v_{F}\int_{{\bf p}}\Psi_{K^{\prime}}^{+}({\bf p})({\bf
p})(\sigma_{x}p_{x}-\sigma_{y}p_{y})\Psi_{K^{\prime}}({\bf p}),$ (29)

where $v_{F}=\frac{3}{2}ta$ is the Fermi velocity. The Green functions in the
Keldysh formalism have a matrix structure given by

$\hat{G}(\epsilon,{\bf p})=\left(\begin{array}[]{cc}G^{R}(\epsilon,{\bf
p})&G^{K}(\epsilon,{\bf p})\\\ 0&G^{A}(\epsilon,{\bf p})\end{array}\right).$
(30)

Here the superscripts $R/A/K$ relate to the retarded, advanced, and Keldysh
components of the Green function respectively. According to the Hamiltonian
Eq. (28), close to the $K$-point, the explicit form of the Green functions in
equilibrium at temperature $T$ is given by

$\displaystyle G^{R/A}_{K}(\epsilon,{\bf p})=\frac{\epsilon{\bf
1}_{2}+v_{F}{\bf p}\cdot\bm{\sigma}}{\epsilon^{2}-v_{F}^{2}p^{2}\pm
io}=\frac{1}{\epsilon^{2}-v_{F}^{2}p^{2}\pm
io}\left(\begin{array}[]{cc}\epsilon&v_{F}(p_{x}-ip_{y})\\\
v_{F}(p_{x}+ip_{y})&\epsilon\end{array}\right)$ (33) $\displaystyle
G^{K}_{K}(\epsilon,{\bf
p})=\tanh\left(\frac{\epsilon}{2T}\right)\left(G^{R}(\epsilon,{\bf
p})-G^{A}(\epsilon,{\bf p})\right).$ (34)

In what follows we introduce complexified momentum $p=p_{x}+ip_{y}$, denoting
$\bar{p}=p_{x}-ip_{y}$, and $|p|=\sqrt{p_{x}^{2}+p_{y}^{2}}$. Then the
expressions for the Green functions can be represented in the form that
exhibits their pole structure explicitly

$G_{K}^{R/A}(\epsilon,{\bf p})=\frac{1}{2}\left[\frac{1}{\epsilon-v_{F}|p|\pm
io}+\frac{1}{\epsilon+v_{F}|p|\pm
io}\right]\left(\begin{array}[]{cc}1&v_{F}\bar{p}/\epsilon\\\
v_{F}p/\epsilon&1\end{array}\right),$ (35) $G^{K}_{K}(\epsilon,{\bf
p})=-i\pi\tanh\left(\frac{\epsilon}{2T}\right)\left(\begin{array}[]{cc}1&v_{F}p/\epsilon\\\
v_{F}\bar{p}/\epsilon&1\end{array}\right)[\delta(\epsilon-
v_{F}|p|)+\delta(\epsilon+v_{F}|p|)].$ (36)

The Green functions for momenta close to the $K^{\prime}$ point are determined
by the Hamiltonian Eq. (29). They read

$G^{R/A}_{K^{\prime}}(\epsilon,{\bf p})=\frac{\epsilon{\bf 1}_{2}+v_{F}{\bf
p}\cdot\bm{\sigma}^{*}}{\epsilon^{2}-v_{F}^{2}p^{2}\pm
io}=\frac{1}{\epsilon^{2}-v_{F}^{2}p^{2}\pm
io}\left(\begin{array}[]{cc}\epsilon&v_{F}(p_{x}+ip_{y})\\\
v_{F}(p_{x}-ip_{y})&\epsilon\end{array}\right),$ (37)

or, in the form explicitely exibiting the pole structure,

$G_{K^{\prime}}^{R/A}(\epsilon,{\bf p})=\frac{1}{2}\left[\frac{1}{\epsilon-
v_{F}|p|\pm io}+\frac{1}{\epsilon+v_{F}|p|\pm
io}\right]\left(\begin{array}[]{cc}1&v_{F}p/\epsilon\\\
v_{F}\bar{p}/\epsilon&1\end{array}\right),$ (38)
$G^{K}_{K^{\prime}}(\epsilon,{\bf
p})=-i\pi\tanh\left(\frac{\epsilon}{2T}\right)\left(\begin{array}[]{cc}1&v_{F}p/\epsilon\\\
v_{F}\bar{p}/\epsilon&1\end{array}\right)[\delta(\epsilon-
v_{F}|p|)+\delta(\epsilon+v_{F}|p|)].$ (39)

Since the Green function in different valleys differ only by the matrix
structure in the pseudospin space, it is convenient to introduce the following
representation

$\hat{G}_{\nu}^{R/A}(\epsilon,{\bf p})=G^{R/A}(\epsilon,{\bf
p})\hat{g}_{\nu}(\epsilon,{\bf p}),\,\,\,G_{\nu}^{K}(\epsilon,{\bf
p})=G^{K}(\epsilon,{\bf p})\hat{g}_{\nu}(\epsilon,{\bf p}),$ (40)

where

$\displaystyle G^{R/A}(\epsilon,{\bf p})=\frac{1}{2}\left[\frac{1}{\epsilon-
v_{F}|p|\pm io}+\frac{1}{\epsilon+v_{F}|p|\pm io}\right],$ (41) $\displaystyle
G^{K}(\epsilon,{\bf
p})=-i\pi\tanh\left(\frac{\epsilon}{2T}\right)[\delta(\epsilon-
v_{F}|p|)+\delta(\epsilon+v_{F}|p|)],$ (42)

denote the parts of the Green functions independent of their pseudospin
structure, and

$\displaystyle\hat{g}_{K^{\prime}}(\epsilon,{\bf
p})=\left(\begin{array}[]{cc}1&v_{F}p/\epsilon\\\
v_{F}\bar{p}/\epsilon&1\end{array}\right)=\sigma_{0}+\frac{v_{F}}{\epsilon}(p_{x}\sigma_{x}-p_{y}\sigma_{y}),$
(45) $\displaystyle\hat{g}_{K}(\epsilon,{\bf
p})=\left(\begin{array}[]{cc}1&v_{F}\bar{p}/\epsilon\\\
v_{F}p/\epsilon&1\end{array}\right)=\sigma_{0}+\frac{v_{F}}{\epsilon}(p_{x}\sigma_{x}+p_{y}\sigma_{y}),$
(48)

capture the structure of Green functions in the pseudospin space.

##  V Model of the drag potential by moving polar crystalline layer

In this section we provide detailed derivation of the interaction constants
$u_{i}({\bf q,p})$ used in the main text of the paper. As it is explained in
the main text, the dominant contirbution to the drag is provided to the
Fourier component of the crystalline potential that has a maximal projection
upon the direction of velocity of the crystalline layer. In what follows we
consider just that component, denoting its wave vector as ${\bf Q}$.

###  V.1 Modulation of on-site energies

We introduce the basis wave functions $\psi_{A,B}({\bf r})$ localized at the
positions of the atoms of the two ($A,B$) graphene sublattices. Thereby the
atoms of the $A$-sublattice are situated at ${\bf r}_{nm}=n{\bf a}_{1}+m{\bf
a}_{2}$, where $n,m$ are integers and ${\bf a}_{1,2}$ denote the lattice
vectors in graphene, whereas the positions of atoms of the $B$-sublattice are
given by ${\bf r}_{nm}+{\bm{\delta}}_{3}$. The action of the periodic
potential on the on-site energies of the atoms is given by

$\displaystyle H_{U}=U_{0}\sum_{n,m}\left\\{\cos({\bf Q}{\bf r}_{nm}-\omega
t)\bar{\psi}_{A}({\bf r}_{nm})\psi_{A}({\bf r}_{nm})+\cos({\bf Q}({\bf
r}_{nm}+{\bm{\delta}}_{3})-\omega t)\bar{\psi}_{B}({\bf
r}_{nm}+{\bm{\delta}}_{3})\psi_{B}({\bf r}_{nm}+{\bm{\delta}}_{3})\right\\}$

Perfoming Fourier-transform according to Eq. (17), we cast the external
potential to the following expression

$\displaystyle H_{U}=\frac{U_{0}}{2}e^{-i\omega t}e^{\frac{i}{2}{\bf
Q}{\bm{\delta}}_{3}}\int_{\bf k}\left\\{\cos\left(\frac{{\bf
Q}}{2}{\bm{\delta}}_{3}\right)\left[\bar{\psi}_{A}({\bf k}+{\bf
Q})\psi_{A}({\bf k})+\bar{\psi}_{B}({\bf k}+{\bf Q})\psi_{B}({\bf
k})\right]\right.$ $\displaystyle\left.-i\sin\left(\frac{{\bf
Q}}{2}{\bm{\delta}}_{3}\right)\left[\bar{\psi}_{A}({\bf k}+{\bf
Q})\psi_{A}({\bf k})-\bar{\psi}_{B}({\bf k}+{\bf Q})\psi_{B}({\bf
k})\right]\right\\}+$ $\displaystyle e^{i\omega t}e^{-\frac{i}{2}{\bf
Q}{\bm{\delta}}_{3}}\int_{\bf k}\left\\{\cos\left(\frac{{\bf
Q}}{2}{\bm{\delta}}_{3}\right)\left[\bar{\psi}_{A}({\bf k}-{\bf
Q})\psi_{A}({\bf k})+\bar{\psi}_{B}({\bf k}-{\bf Q})\psi_{B}({\bf
k})\right]\right.$ $\displaystyle\left.+i\sin\left(\frac{{\bf
Q}}{2}{\bm{\delta}}_{3}\right)\left[\bar{\psi}_{A}({\bf k}-{\bf
Q})\psi_{A}({\bf k})-\bar{\psi}_{B}({\bf k}-{\bf Q})\psi_{B}({\bf
k})\right]\right\\}.$ (49)

Fruthermore, using the pseudospin spinors $\Psi({\bf k})=(\psi_{A}({\bf
k}),\psi_{B}({\bf k}))^{T}$, the Hamiltonian $H_{U}$ can be written as

$\displaystyle H_{U}=\frac{U_{0}}{2}e^{-i\omega t}e^{\frac{i}{2}{\bf
Q}{\bm{\delta}}_{3}}\int_{\bf k}\left\\{\cos\left(\frac{{\bf
Q}}{2}{\bm{\delta}}_{3}\right)\bar{\Psi}({\bf k}+{\bf Q})\sigma_{0}\Psi({\bf
k})-i\sin\left(\frac{{\bf Q}}{2}{\bm{\delta}}_{3}\right)\bar{\Psi}({\bf
k}+{\bf Q})\sigma_{z}\Psi({\bf k})\right\\}+$ $\displaystyle e^{i\omega
t}e^{-\frac{i}{2}{\bf Q}{\bm{\delta}}_{3}}\int_{\bf
k}\left\\{\cos\left(\frac{{\bf Q}}{2}{\bm{\delta}}_{3}\right)\bar{\Psi}({\bf
k}-{\bf Q})\sigma_{0}\Psi({\bf k})+i\sin\left(\frac{{\bf
Q}}{2}{\bm{\delta}}_{3}\right)\bar{\Psi}({\bf k}-{\bf Q})\sigma_{z}\Psi({\bf
k})\right\\},$ (50)

where $\sigma_{0}={\bf 1}_{2}$ and $\sigma_{z}$ denotes the Pauli matrix.
Assuming ${\bf Q}$ to be close to ${\bf K^{\prime}}-{\bf K}$, we introduce the
parametrization ${\bf Q}={\bf K^{\prime}}-{\bf K}+{\bf q}$. The external
potential induces transitions between the electron states in different
valleys. Then, assuming $|q|\ll|Q|$, we can represent Eq. (50) as

$\displaystyle H_{U}=\frac{U_{0}}{2}e^{-i\omega t}e^{\frac{i}{2}({\bf
K^{\prime}}-{\bf K}+{\bf q}){\bm{\delta}}_{3}}\int_{\bf
p}\left\\{\cos\left(\frac{{\bf K^{\prime}}-{\bf K}+{\bf
q}}{2}{\bm{\delta}}_{3}\right)\bar{\Psi}_{K^{\prime}}({\bf p}+{\bf
q})\sigma_{0}\Psi_{K}({\bf p})-i\sin\left(\frac{{\bf K^{\prime}}-{\bf K}+{\bf
q}}{2}{\bm{\delta}}_{3}\right)\bar{\Psi}_{K^{\prime}}({\bf p}+{\bf
q})\sigma_{z}\Psi_{K}({\bf p})\right\\}$ $\displaystyle+e^{i\omega
t}e^{-\frac{i}{2}({\bf K^{\prime}}-{\bf K}+{\bf q}){\bm{\delta}}_{3}}\int_{\bf
p}\left\\{\cos\left(\frac{{\bf K^{\prime}}-{\bf K}+{\bf
q}}{2}{\bm{\delta}}_{3}\right)\bar{\Psi}_{K}({\bf
p})\sigma_{0}\Psi_{K^{\prime}}({\bf p}+{\bf q})+i\sin\left(\frac{{\bf
K^{\prime}}-{\bf K}+{\bf q}}{2}{\bm{\delta}}_{3}\right)\bar{\Psi}_{K}({\bf
p})\sigma_{z}\Psi_{K^{\prime}}({\bf p}+{\bf q})\right\\},$ (51)

Further simplification is reached by noticing that $({\bf K^{\prime}}-{\bf
K})\cdot{\bm{\delta}}_{3}=0$, which results in

$\displaystyle H_{U}=\frac{U_{0}}{2}e^{-i\omega t}e^{\frac{i}{2}{\bf
q}{\bm{\delta}}_{3}}\int_{\bf p}\left\\{\cos\left(\frac{{\bf
q}}{2}{\bm{\delta}}_{3}\right)\bar{\Psi}_{K^{\prime}}({\bf p}+{\bf
q})\sigma_{0}\Psi_{K}({\bf p})-i\sin\left(\frac{{\bf
q}}{2}{\bm{\delta}}_{3}\right)\bar{\Psi}^{K^{\prime}}({\bf p}+{\bf
q})\sigma_{z}\Psi_{K}({\bf p})\right\\}+$ $\displaystyle e^{i\omega
t}e^{-\frac{i}{2}{\bf q}{\bm{\delta}}_{3}}\int_{\bf
p}\left\\{\cos\left(\frac{{\bf
q}}{2}{\bm{\delta}}_{3}\right)\bar{\Psi}_{K}({\bf
p})\sigma_{0}\Psi_{K^{\prime}}({\bf p}+{\bf q})+i\sin\left(\frac{{\bf
q}}{2}{\bm{\delta}}_{3}\right)\bar{\Psi}_{K}({\bf
p})\sigma_{z}\Psi_{K^{\prime}}({\bf p}+{\bf q})\right\\}.$ (52)

###  V.2 Change of the hopping amplitude by the potential of the moving
lattice

We assume that the change of the hopping amplitude between the nearest
neighbor atoms is proportional to the local external potential in the middle
of the link connecting the two atoms in the graphene lattice. This assumption
leads to the following Hamiltonian for the change of the nearest-neighbor
hopping

$H_{h}=U_{1}\sum_{{\bm{\delta}}_{1},{\bm{\delta}}_{2},{\bm{\delta}}_{3}}\sum_{n,m}\left\\{\cos\left[{\bf
Q}\left({\bf r}_{nm}+\frac{{\bm{\delta}}_{\nu}}{2}\right)-\omega
t\right]\left[\bar{\psi}_{A}({\bf r}_{nm})\psi_{B}({\bf
r}_{nm}+{\bm{\delta}}_{\nu})+\bar{\psi}_{B}({\bf
r}_{nm}+{\bm{\delta}}_{\nu})\psi_{A}({\bf r}_{nm})\right]\right\\}$ (53)

Here the coupling constant $U_{1}$ is in general different from the coupling
$U_{0}$ that describes the change of the local on-site potential. Performing
the Fourier transformation in Eq. (53) into the ${\bf k}$ space and leaving
only the states close to $K$ and $K^{\prime}$ point, we cast the Fourier
transformed Hamiltonian $H_{h}$ to the form

$\displaystyle H_{h}=\frac{U_{1}}{2}\sum_{i=1}^{3}\int_{\bf
p}\left\\{e^{-i\omega t}e^{\frac{i}{2}{\bm{\delta}}_{3}({\bf K^{\prime}}-{\bf
K}+{\bf q})}\bar{\Psi}_{K^{\prime}}({\bf q}+{\bf
p})\left(\sigma_{x}\cos\left[\left(\frac{{\bf K}+{\bf K^{\prime}}+{\bf
q}}{2}+{\bf
p}\right)({\bm{\delta}}_{i}-{\bm{\delta}}_{3})\right]-\right.\right.$
$\displaystyle\left.\sigma_{y}\sin\left[\left(\frac{{\bf K}+{\bf
K^{\prime}}+{\bf q}}{2}+{\bf
p}\right)({\bm{\delta}}_{i}-{\bm{\delta}}_{3})\right]\right)\Psi_{K}({\bf
p})+$ $\displaystyle e^{i\omega t}e^{-\frac{i}{2}{\bm{\delta}}_{3}({\bf
K^{\prime}}-{\bf K}+{\bf q})}\bar{\Psi}_{K}({\bf
p})\left(\sigma_{x}\cos\left[\left(\frac{{\bf K}+{\bf K^{\prime}}+{\bf
q}}{2}+{\bf p}\right)({\bm{\delta}}_{i}-{\bm{\delta}}_{3})\right]-\right.$
$\displaystyle\left.\left.\sigma_{y}\sin\left[\left(\frac{{\bf K}+{\bf
K^{\prime}}+{\bf q}}{2}+{\bf
p}\right)({\bm{\delta}}_{i}-{\bm{\delta}}_{3})\right]\right)\Psi_{K^{\prime}}({\bf
p}+{\bf q})\right\\}.$ (54)

###  V.3 Total drag Hamiltonian

Upon projecting on the states close to the $K$ points, the drag potential
acquires the form

$H_{d}=\frac{U_{0}}{2}\int_{\bf p}\left\\{\Psi^{+}_{K^{\prime}}({\bf p}+{\bf
q})\left({\bf u}({\bf q},{\bf p})\cdot{\bm{\sigma}}\right)\Psi_{K}({\bf
p})e^{-i\omega t}+\Psi^{+}_{K}({\bf p})\left({\bf u}^{*}({\bf q},{\bf
p})\cdot{\bm{\sigma}}\right)\Psi_{K^{\prime}}({\bf p}+{\bf q})e^{i\omega
t}\right\\},$ (55)

where $\left({\bf u}({\bf
q})\cdot{\bm{\sigma}}\right)=\sum_{i=0}^{3}u_{i}({\bf q})\sigma_{i}$. The
components of the vector of coupling constants ${\bf u}$ are obtained by
comparison with Eqs. (52), (54) as follows

$\displaystyle u_{0}=\frac{1}{2}e^{\frac{i}{2}{\bf
q}{\bm{\delta}}_{3}}\cos\left(\frac{{\bf q}}{2}{\bm{\delta}}_{3}\right),$ (56)
$\displaystyle u_{3}=-i\frac{1}{2}e^{\frac{i}{2}{\bf
q}{\bm{\delta}}_{3}}\sin\left(\frac{{\bf q}}{2}{\bm{\delta}}_{3}\right),$ (57)
$\displaystyle u_{1}=\frac{U_{1}}{2U_{0}}e^{\frac{i}{2}{\bf
q}{\bm{\delta}}_{3}}\sum_{i=1}^{3}\cos\left[\left(\frac{{\bf K}+{\bf
K^{\prime}}+{\bf q}}{2}+{\bf
p}\right)({\bm{\delta}}_{i}-{\bm{\delta}}_{3})\right]=\frac{U_{1}}{2U_{0}}e^{\frac{i}{2}{\bf
q}{\bm{\delta}}_{3}}\left(1-\sum_{i=1}^{2}\cos\left[\left(\frac{{\bf
q}}{2}+{\bf p}\right)({\bm{\delta}}_{i}-{\bm{\delta}}_{3})\right]\right),$
(58) $\displaystyle u_{2}=-\frac{U_{1}}{2U_{0}}e^{\frac{i}{2}{\bf
q}{\bm{\delta}}_{3}}\sum_{i=1}^{3}\sin\left[\left(\frac{{\bf K}+{\bf
K^{\prime}}+{\bf q}}{2}+{\bf
p}\right)({\bm{\delta}}_{i}-{\bm{\delta}}_{3})\right]=\frac{U_{1}}{2U_{0}}e^{\frac{i}{2}{\bf
q}{\bm{\delta}}_{3}}\sum_{i=1}^{2}\sin\left[\left(\frac{{\bf q}}{2}+{\bf
p}\right)({\bm{\delta}}_{i}-{\bm{\delta}}_{3})\right],$ (59)

where in the second equalities we took into account $\frac{{\bf K}+{\bf
K^{\prime}}}{2}\cdot({\bm{\delta}}_{1}-{\bm{\delta}}_{3})=\frac{{\bf K}+{\bf
K^{\prime}}}{2}\cdot({\bm{\delta}}_{2}-{\bm{\delta}}_{3})=\pi$. Here we use
the coupling $U_{0}$ as the overall scale for the coupling strengths.

##  VI Projection on the conduction band

Due to the necessity of finite doping, the main contribution to the drag
current at low temperature is produced by transitions between the quantum
states in the conduction band. This observation allows simplified calculation
of the current by projection on the conduction band. To realize the
projection, we first employ the unitary rotation that diagonalizes the
Hamiltonian in the pseudospin sector. Let $\phi_{p}$ be the polar angle of the
wave vector ${\bf p}=(p_{x},p_{y})$. Then the Hamiltonian of the state with
the wave vector ${\bf p}$ is diagonalized by the unitary transformation

$H({\bf p})=\left(\begin{array}[]{cc}v_{F}|p|&0\\\
0&-v_{F}|p|\end{array}\right)=U_{\nu}^{+}(\phi_{p})H_{\nu}({\bf
p})U_{\nu}(\phi_{p}),$ (60)

where

$U_{K}({\bf p})=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}e^{-i\phi_{\bf
p}/2}&e^{-i\phi_{\bf p}/2}\\\ e^{i\phi_{\bf p}/2}&-e^{i\phi_{\bf
p}/2}\end{array}\right),\,\,\,\,U_{K^{\prime}}({\bf
p})=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}e^{i\phi_{\bf
p}/2}&e^{i\phi_{\bf p}/2}\\\ e^{-i\phi_{\bf p}/2}&-e^{-i\phi_{\bf
p}/2}\end{array}\right)$ (61)

Further projection on the conduction bend is achieved by the application of
the operator $\frac{1}{2}(1+\sigma_{z})$, so that the total projection
operator reads

$\hat{P}^{c}_{\nu}({\bf p})=U_{\nu}({\bf p})\left(\begin{array}[]{cc}1&0\\\
0&0\end{array}\right)$ (62)

Explicitly, close to the points $K$ and $K^{\prime}$, the eigenstates of the
conduction band are expessed through the eigenstates of $A$ and $B$
sublattices as follows

$|p\rangle_{c}^{K}=\frac{1}{\sqrt{2}}\left(e^{-i\phi_{\bf
p}/2}|A\rangle+e^{i\phi_{\bf
p}/2}|B\rangle\right),\,\,\,\,|p\rangle_{c}^{K^{\prime}}=\frac{1}{\sqrt{2}}\left(e^{i\phi_{\bf
p}/2}|A\rangle+e^{-i\phi_{\bf p}/2}|B\rangle\right)$ (63)

Let us now consider the kinematic restriction on the wave vectors of the
initial and final states due to the energy and quasi-momentum conservation in
detail. Since the energy transfer by the scattering processes considered here
is much smaller than the Fermi energy, we can approximately set the absolute
value of the wave vectors for the initial and final states equal to the Fermi
wave vector $p_{F}$. It turns out that this condition leaves only two possible
choices for the wave vector of the initial state ${\bf p}$ (the wave vector of
the final state is then fixed automatically to ${\bf p}+{\bf q}$) for the
fixed transferred wave vector ${\bf q}$ (see Fig. 4). Denote the polar angle
of the wave vector $\bf q$ as $\alpha$ and the angle between vectors ${\bf q}$
and ${\bf p}$ as $\phi$. Then the two possible choices are

Figure 4: Possible directions of initial and finite wave vectors for a given
transfer wave vector ${\bf q}$

$\displaystyle{\bf i)}\,\,\,{\phi_{\bf p}}=\alpha+\phi,\,\,\,{\phi_{\bf
p^{\prime}}}=\pi+\alpha-\phi,$ (64) $\displaystyle{\bf ii)}\,\,\,{\phi_{\bf
p}}=\alpha-\phi,\,\,\,{\phi_{\bf
p^{\prime}}}=\alpha+\phi-\pi<=>\pi+\alpha+\phi.$ (65)

It follows that all results for the case (ii) are obtained from (i) by change
$\phi\rightarrow-\phi$. The change of the sign of the trasnfer quasi-momentum
${\bf q}\rightarrow-{\bf q}$ corresponds to the change
$\alpha\rightarrow\alpha+\pi$.

Let us write down explicit expression for the components of current density
operators Eqs. (25) projected on the conduction band. We consider the geometry
as given in the case (i) above.

$\displaystyle{\bf j}_{K}({\bf p})=-ta\left\\{\left(\begin{array}[]{c}-1\\\
0\end{array}\right)+\left(\begin{array}[]{c}1/2\\\
\sqrt{3}/2\end{array}\right)e^{i\frac{2\pi}{3}}e^{i\sqrt{3}ap_{F}\cos(\alpha+\phi-\pi/6)}+\left(\begin{array}[]{c}1/2\\\
-\sqrt{3}/2\end{array}\right)e^{-i\frac{2\pi}{3}}e^{i\sqrt{3}ap_{F}\cos(\alpha+\phi+\pi/6)}\right\\}$
(72) $\displaystyle{\bf j}_{K^{\prime}}({\bf p}+{\bf q})={\bf
j}_{K^{\prime}}({\bf p^{\prime}})=-ta\left\\{\left(\begin{array}[]{c}-1\\\
0\end{array}\right)+\left(\begin{array}[]{c}1/2\\\
\sqrt{3}/2\end{array}\right)e^{-i\frac{2\pi}{3}}e^{-i\sqrt{3}ap_{F}\cos(\alpha-\phi-\pi/6)}+\left(\begin{array}[]{c}1/2\\\
-\sqrt{3}/2\end{array}\right)e^{i\frac{2\pi}{3}}e^{-i\sqrt{3}ap_{F}\cos(\alpha-\phi+\pi/6)}\right\\}$
(79) (80)

The projection of the current density operator on the conduction band is
calculated using Eqs. (72), (80) as follows

$\displaystyle\langle{\bf p}+{\bf q},K^{\prime}|{\bf J}_{K^{\prime}K}({\bf
p}+{\bf q},{\bf p})|{\bf
p},K\rangle_{c}=\frac{1}{2}\left(e^{-i\phi^{\prime}/2},e^{i\phi^{\prime}/2}\right)\left(\begin{array}[]{cc}0&{\bf
j}_{K}({\bf p})\\\ {\bf j}_{K^{\prime}}^{*}({\bf p}+{\bf
q})&0\end{array}\right)\left(\begin{array}[]{c}e^{-i\phi/2}\\\
e^{i\phi/2}\end{array}\right)=$ (85) $\displaystyle\frac{1}{2}\left\\{{\bf
j}_{K}({\bf p})e^{i(\phi-\phi^{\prime})/2}+{\bf j}_{K^{\prime}}^{*}({\bf
p}+{\bf q})e^{-i(\phi-\phi^{\prime})/2}\right\\},$ (86)
$\displaystyle\langle{\bf p},K|{\bf J}_{KK^{\prime}}({\bf p},{\bf p}+{\bf
q})|{\bf p}+{\bf q},K^{\prime}\rangle_{c}=\langle{\bf p}+{\bf
q},K^{\prime}|{\bf J}_{K^{\prime}K}({\bf p}+{\bf q},{\bf p})|{\bf
p},K\rangle_{c}^{*}=$ $\displaystyle\frac{1}{2}\left\\{{\bf j}^{*}_{K}({\bf
p})e^{-i(\phi-\phi^{\prime})/2}+{\bf j}_{K^{\prime}}({\bf p}+{\bf
q})e^{i(\phi-\phi^{\prime})/2}\right\\},$ (87)

where $\phi=\phi_{\bf p}$ and $\phi^{\prime}=\phi_{\bf p^{\prime}}$ as defined
by Eqs. (64) , (65).

For further calculations we also need the expressions for the projected
interaction vertices

$\displaystyle\langle{\bf p+q},K^{\prime}|u_{0}\sigma_{0}|{\bf
p},K\rangle_{c}=-\frac{1}{2}e^{-\frac{i}{2}|q|a\cos\alpha}\cos\left(\frac{|q|a}{2}\cos\alpha\right)\sin\alpha,$
(88) $\displaystyle\langle{\bf p+q},K^{\prime}|u_{3}\sigma_{z}|{\bf
p},K\rangle_{c}=\frac{1}{2}e^{-\frac{i}{2}|q|a\cos\alpha}\sin\left(\frac{|q|a}{2}\cos\alpha\right)\cos\alpha,$
(89) $\displaystyle\langle{\bf p+q},\\!K^{\prime}|u_{1}\sigma_{x}|{\bf
p},\\!K\rangle_{c}\\!=\\!\frac{U_{1}}{2U_{0}}e^{-\frac{i}{2}|q|a\cos\alpha}\sin\phi\left\\{1\\!-\\!\cos\left[p_{F}a\sin\phi\cos\left(\frac{\pi}{3}\\!-\\!\alpha\right)\right]\\!-\\!\cos\left[p_{F}a\sin\phi\cos\left(\frac{2\pi}{3}\\!-\\!\alpha\right)\right]\right\\},$
(90) $\displaystyle\langle{\bf p+q},K^{\prime}|u_{2}\sigma_{y}|{\bf
p},K\rangle_{c}=-\frac{U_{1}}{2U_{0}}e^{-\frac{i}{2}|q|a\cos\alpha}\cos\phi\left\\{\sin\left[p_{F}a\sin\phi\cos\left(\frac{\pi}{3}\\!-\\!\alpha\right)\right]\\!+\\!\sin\left[p_{F}a\sin\phi\cos\left(\frac{2\pi}{3}\\!-\\!\alpha\right)\right]\right\\}.$

##  VII Perturbative calculation of the AC drag current

Figure 5: Diagram for the lowest order contribution to the AC current.

The diagram in Fig. 5 corresponds to the lowest order contribution to the
current density due to the intervalley scattering $(KK^{\prime})$. The
mathematical expression reads

$\langle J_{K^{\prime}K}({\bf
q},\omega)\rangle=-iU_{0}\mathrm{Tr}\left\\{\gamma^{q}{\bf
J}_{K^{\prime}K}({\bf p}+{\bf q},{\bf p}){\hat{G}}_{K}({\bf
p},\epsilon)\left({\bf u}^{*}({\bf
q,p})\cdot{\bm{\sigma}}\right)\gamma^{\mathrm{cl}}{\hat{G}}_{K^{\prime}}({\bf
p}+{\bf q},\epsilon+\omega)\right\\}.$ (92)

Here $\gamma^{q}=\sigma_{1}$ and $\gamma^{\mathrm{cl}}=\sigma_{0}$ denote the
Pauli matrices acting in the Keldysh (RA) space.

Using the notations introduced in Eq. (40), we separate the parts reflecting
the pole structure of Green functions (the spectrum) and the part
corresponding to the matrix elements of the transtion in the pseudo-spin and
valley spaces

$\displaystyle\langle J_{K^{\prime}K}({\bf
q},\omega)\rangle=-iU_{0}\int_{\epsilon,{\bf
p}}I_{K^{\prime}K}\left\\{G^{K}({\bf p},\epsilon)G^{A}({\bf p}+{\bf
q},\epsilon+\omega)+G^{R}({\bf p},\epsilon)G^{K}({\bf p}+{\bf
q},\epsilon+\omega)\right\\}=$ $\displaystyle-iU_{0}\int_{\epsilon,{\bf
p}}I_{K^{\prime}K}\left\\{\tanh\left(\frac{\epsilon-\mu}{2T}\right)\left[G^{R}({\bf
p},\epsilon)-G^{A}({\bf p},\epsilon)\right]G^{A}({\bf p}+{\bf
q},\epsilon+\omega)+\right.$ $\displaystyle\left.G^{R}({\bf
p},\epsilon)\left[G^{R}({\bf p}+{\bf q},\epsilon+\omega)-G^{A}({\bf p}+{\bf
q},\epsilon+\omega)\right]\tanh\left(\frac{\epsilon+\omega-\mu}{2T}\right)\right\\}=$
$\displaystyle-iU_{0}\int_{\epsilon,{\bf
p}}I_{K^{\prime}K}\left\\{\left[\tanh\left(\frac{\epsilon-\mu}{2T}\right)-\tanh\left(\frac{\epsilon+\omega-\mu}{2T}\right)\right]G^{R}({\bf
p},\epsilon)G^{A}({\bf p}+{\bf q},\epsilon+\omega)\right\\},$ (93)

where in the last line we omitted the products of two retarded and two
advanced Green functions, and introduced the notation

$I_{K^{\prime}K}(\epsilon,{\bf p};\omega,{\bf q})=\mathrm{tr}\left\\{{\bf
J}_{K^{\prime}K}({\bf p},{\bf p}+{\bf q})g_{K}({\bf p},\epsilon)\left({\bf
u}^{*}({\bf q,p})\cdot{\bm{\sigma}}\right)g_{K^{\prime}}({\bf p}+{\bf
q},\epsilon+\omega)\right\\}$ (94)

Analogously, the contribution to the current by the scattering from the valley
$K^{\prime}$ to the valley $K$ is given by

$\displaystyle\langle J_{KK^{\prime}}({\bf
q},\omega)\rangle=-iU_{0}\int_{\epsilon,{\bf
p}}I_{KK^{\prime}}\left\\{\left[\tanh\left(\frac{\epsilon+\omega-\mu}{2T}\right)-\tanh\left(\frac{\epsilon-\mu}{2T}\right)\right]G^{R}({\bf
p+q},\epsilon+\omega)G^{A}({\bf p},\epsilon)\right\\},$ (95)

where

$I_{KK^{\prime}}(\epsilon,{\bf p};\omega,{\bf q})=\mathrm{tr}\left\\{{\bf
J}_{KK^{\prime}}({\bf p},{\bf p}+{\bf q})g_{K^{\prime}}({\bf
p+q},\epsilon+\omega)\left({\bf u}({\bf
q,p})\cdot{\bm{\sigma}}\right)g_{K}({\bf p},\epsilon)\right\\}.$ (96)

###  VII.1 Momentum and frequency integration

Closing the integration contour in the lower $\epsilon$ half-plane, we reduce
the integral to the sum over the residue of the poles of the retarded Green
function, which results in

$\displaystyle\langle J_{K^{\prime}K}({\bf q},\omega)\rangle\frac{-\pi
U_{0}}{2}\int_{\epsilon,{\bf
p}}I_{K^{\prime}K}\left[\tanh\left(\frac{\epsilon+\omega-\mu}{2T}\right)-\tanh\left(\frac{\epsilon-\mu}{2T}\right)\right]\times$
$\displaystyle\left\\{\delta(\epsilon-
v_{F}|p|)\left[\frac{1}{\omega+v_{F}|p|-v_{F}|p+q|-i/\tau}+\frac{1}{\omega+v_{F}|p|+v_{F}|p+q|-i/\tau}\right]+\right.$
$\displaystyle\left.\delta(\epsilon+v_{F}|p|)\left[\frac{1}{\omega-
v_{F}|p|-v_{F}|p+q|-i/\tau}+\frac{1}{\omega-
v_{F}|p|+v_{F}|p+q|-i/\tau}\right]\right\\}$ (97)

$\displaystyle\langle J_{KK^{\prime}}({\bf q},\omega)\rangle=\frac{-\pi
U_{0}}{2}\int_{\epsilon,{\bf
p}}I_{KK^{\prime}}\left[\tanh\left(\frac{\epsilon+\omega-\mu}{2T}\right)-\tanh\left(\frac{\epsilon-\mu}{2T}\right)\right]\times$
$\displaystyle\left\\{\delta(\epsilon+\omega-v_{F}|{\bf
p+q}|)\left[\frac{1}{-\omega+v_{F}|{\bf
p+q}|-v_{F}|p|-i/\tau}+\frac{1}{-\omega+v_{F}|{\bf
p+q}|+v_{F}|p|-i/\tau}\right]+\right.$
$\displaystyle\left.\delta(\epsilon+\omega+v_{F}|{\bf
p+q}|)\left[\frac{1}{-\omega-v_{F}|{\bf
p+q}|-v_{F}|p|-i/\tau}+\frac{1}{-\omega-v_{F}|{\bf
p+q}|+v_{F}|p|-i/\tau}\right]\right\\}$ (98)

Here we introduced small imaginary part at the poles
$io\rightarrow\frac{i}{2\tau}$, which accounts phenomenologically for
intrinsic scattering processes in graphene.

Furthermore, leaving only the pole contributions to the integrals over
momentum ${\bf p}$, we obtain

$\displaystyle\langle J_{K^{\prime}K}({\bf
q},\omega)\rangle=\frac{i\pi^{2}U_{0}}{2}\int_{{\bf
p}}\left\\{I_{K^{\prime}K}(\epsilon=v_{F}|p|-\omega)\left[\tanh\left(\frac{\omega+v_{F}|p|-\mu}{2T}\right)-\tanh\left(\frac{v_{F}|p|-\mu}{2T}\right)\right]\times\right.$
$\displaystyle\left[\delta(\omega+v_{F}|p|-v_{F}|p+q|)+\delta(\omega+v_{F}|p|+v_{F}|p+q|)\right]+$
$\displaystyle
I_{K^{\prime}K}(\epsilon=-v_{F}|p|)\left[\tanh\left(\frac{\omega-
v_{F}|p|-\mu}{2T}\right)-\tanh\left(\frac{-v_{F}|p|-\mu}{2T}\right)\right]\times$
$\displaystyle\left.\left[\delta(\omega-v_{F}|p|-v_{F}|p+q|)+\delta(\omega-
v_{F}|p|+v_{F}|p+q|)\right]\right\\}\approx$
$\displaystyle\frac{i\pi^{2}U_{0}}{2}\int_{\bf
p}I_{K^{\prime}K}(\epsilon=v_{F}|p|)\left[\tanh\left(\frac{\omega+v_{F}|p|-\mu}{2T}\right)-\tanh\left(\frac{v_{F}|p|-\mu}{2T}\right)\right]\delta(\omega+v_{F}|p|-v_{F}|p+q|),$
(99)

$\displaystyle\langle J_{KK^{\prime}}({\bf
q},\omega)\rangle=\frac{i\pi^{2}U_{0}}{2}\int_{{\bf
p}}\left\\{I_{KK^{\prime}}(\epsilon=v_{F}|{\bf
p+q}|-\omega)\left[\tanh\left(\frac{v_{F}|{\bf
p+q}|-\mu}{2T}\right)-\tanh\left(\frac{v_{F}|{\bf
p+q}|-\omega-\mu}{2T}\right)\right]\times\right.$
$\displaystyle\left[\delta(v_{F}|{\bf p+q}|-\omega-
v_{F}|p|)+\delta(-\omega+v_{F}|p|+v_{F}|{\bf p+q}|)\right]+$ $\displaystyle
I_{KK^{\prime}}(\epsilon=-\omega-v_{F}|{\bf
p+q}|)\left[\tanh\left(\frac{-v_{F}|{\bf
p+q}|-\mu}{2T}\right)-\tanh\left(\frac{-v_{F}|{\bf
p+q}|-\omega-\mu}{2T}\right)\right]\times$
$\displaystyle\left.\left[\delta(-\omega-v_{F}|{\bf
p+q}|-v_{F}|p|)+\delta(-\omega-v_{F}|{\bf
p+q}|+v_{F}|p|)\right]\right\\}\approx$
$\displaystyle\frac{i\pi^{2}U_{0}}{2}\int_{\bf
p}I_{KK^{\prime}}(\epsilon=v_{F}|p|)\left[\tanh\left(\frac{v_{F}|p|+\omega-\mu}{2T}\right)-\tanh\left(\frac{v_{F}|p|-\mu}{2T}\right)\right]\delta(v_{F}|{\bf
p+q}|-\omega-v_{F}|p|).$ (100)

In the last lines of Eqs. (99), (100) we left the leading contributions at
positive chemical potential and low temperatures. The combination of
$\tanh$-functions in that term together with the energy and quasi-momentum
conservation imposed by the $\delta$-function projects the integrals on the
states of conduction band.

It is convenient to perfom the two-dimensional integration over momenta using
polar coordinates, $p=(|{\bf p}|,\phi)$ where $\phi$ is defined as the angle
between vectors ${\bf p}$ and ${\bf q}$. Furthermore, we satisfy the
$\delta$-functions by perfoming the angular integration. The angular integrals
are calculated according to the formula

$\int d\phi\delta(f(\phi))F(\phi)=\frac{F(\phi_{0})}{|f^{\prime}(\phi_{0})|},$
(101)

where $f^{\prime}(\phi_{0})=\frac{df(\phi)}{d\phi}|_{\phi=\phi_{0}}$, and
$\phi_{0}$ satisfies the relation $f(\phi_{0})=0$. Here we show explicitely
the computation of the integral with $\delta(\omega+v_{F}|p|-v_{F}|p+q|)$ in
Eq. (99). The angular dependence in the $\delta$-function is contained in the
term

$|p+q|=\sqrt{p^{2}+q^{2}+2pq\cos\phi},$ (102)

therefore

$f(\phi)=\omega+v_{F}p-v_{F}\sqrt{p^{2}+q^{2}+2pq\cos\phi},\,\,\,\,f^{\prime}(\phi)=\frac{v_{F}pq\sin\phi}{\sqrt{p^{2}+q^{2}+2pq\cos\phi}}=\frac{v_{F}pq\sin\phi}{|p+q|}=\frac{v_{F}^{2}pq\sin\phi}{\omega+v_{F}p},$
(103)

where in the last equation we used the condition imposed by the
$\delta$-function.

Furhermore, in the physically relevant regime $v_{F}p\gg\omega$, we can
simplify

$f^{\prime}(\phi)\approx v_{F}q\sin\phi.$ (104)

At low temperature, $T\ll\omega$, the $\tanh$ can be replaced by a step
function. Then, the integration over the absolute value of momentum $|p|$ can
be performed by replacing $|p|\approx p_{F}$, whereas the difference of the
step-functions determines the integration range $\omega/v_{F}$. Thus, the low-
temperature approximation to the integration can be summarized by the
expression

$\int_{\bf
p}\left[\tanh\left(\frac{\omega+v_{F}|p|-\mu}{2T}\right)-\tanh\left(\frac{v_{F}|p|-\mu}{2T}\right)\right]\delta(v_{F}|{\bf
p+q}|-\omega-v_{F}|p|)...\approx\frac{\omega p_{F}}{v_{F}^{2}|q|}\int
d|p|\int\frac{d\phi}{\sin\phi}\delta(|p|-p_{F})\delta(\phi-\phi_{0})...,$
(105)

where $\phi_{0}$ is the angle between $\bf p$ and $\bf q$ that solves the
equation $|{\bf p+q|}=|p|=p_{F}$, which is explicitely given by

$\phi_{0}=\arccos\left(-\frac{q}{2p_{F}}\right).$ (106)

Finally, the AC drag current density, given by the sum of Eqs. (99) and (100)
becomes

$\displaystyle{\bf j}_{\mathrm{AC}}=\frac{i\pi^{2}U_{0}}{2}\int_{{\bf
p}}\left(I_{K^{\prime}K}(|p|,{\bf q},\phi)+I_{KK^{\prime}}(|p|,{\bf
q},\phi)\right)\left[\tanh\left(\frac{\omega+v_{F}|p|-\mu}{2T}\right)-\tanh\left(\frac{v_{F}|p|-\mu}{2T}\right)\right]\delta(\omega+v_{F}|p|-v_{F}|p+q|)=$
$\displaystyle\frac{i\pi^{2}U_{0}}{2v_{F}|q|\sin\phi_{0}}\int
pdp\left(I_{K^{\prime}K}(|p|,{\bf q},\phi_{0})+I_{KK^{\prime}}(|p|,{\bf
q},\phi_{0})\right)\left[\tanh\left(\frac{\omega+v_{F}|p|-\mu}{2T}\right)-\tanh\left(\frac{v_{F}|p|-\mu}{2T}\right)\right]\approx$
$\displaystyle\frac{i\pi^{2}U_{0}}{2v_{F}|q|\sin\phi_{0}}\frac{p_{F}\omega}{v_{F}}\left(I_{K^{\prime}K}(p_{F},{\bf
q},\phi_{0})+I_{KK^{\prime}}(p_{F},{\bf
q},\phi_{0})\right)=\frac{i\pi^{2}U_{0}\omega
p_{F}}{2v_{F}|q||\sin\phi_{0}|}\left(\mathcal{I}_{K^{\prime}K}(p_{F},{\bf
q},\phi_{0})+\mathcal{I}_{K^{\prime}K}(p_{F},{\bf q},\phi_{0})\right),$ (107)

where in the last line we introduced the dimensionless parital current
densities $\mathcal{I}_{ab}$ defined by the relation

$I_{ab}(p_{F},{\bf q},\phi_{0})=v_{F}\mathcal{I}_{ab}(p_{F},{\bf
q},\phi_{0}).$ (108)

The drag current given by Eq. (107) formally diverges at the threshold
$p_{F}=\frac{1}{2}(q-\omega/v_{F})$, which corresponds to the angle
$\phi_{0}=\pi$. That divergence is removed if one takes into account the
scattering processes that change momenum of quasipartcles. Therefore, close to
the threshold, one should take into account the finite life-time $\tau$ in
Eqs. (97), (98) explicitely, which amounts to the replacement
$\sin\phi_{0}\rightarrow\frac{1}{\sqrt{v_{F}q\tau}}$ in Eq. (107) at the
cutoff $\phi_{0}=\pi$. The exact dependence of the prefactor in Eq. (107) on
the angle $\phi$ exhibits a crossover from the constant $\sqrt{v_{F}q\tau}$
for the angles satisfying $\sin\phi\ll\frac{1}{\sqrt{v_{F}q\tau}}$ to
$1/\sin\phi$ for $\sin\phi\gg\frac{1}{\sqrt{v_{F}q\tau}}$.

###  VII.2 Drag current density

In terms of the matrix elements projected onto the conduction band, the drag
current density given by Eq. (107) can be cast in the form

$\displaystyle{\bf j}_{\mathrm{AC}}\approx\frac{i\pi^{2}U_{0}\omega
p_{F}}{2v_{F}|q||\sin\phi|}\left\\{\langle{\bf p}+{\bf q},K^{\prime}|{\bf
J}_{K^{\prime}K}({\bf p}+{\bf q},{\bf p})|K,{\bf p}\rangle\langle{\bf
p},K|{\bf u}^{*}{\bm{\sigma}}|{\bf p}+{\bf q},K^{\prime}\rangle+\right.$
$\displaystyle\left.\langle{\bf p},K|{\bf J}_{KK^{\prime}}({\bf p},{\bf
p}+{\bf q})|K^{\prime},{\bf p}+{\bf q}\rangle\langle{\bf p}+{\bf
q},K^{\prime}|{\bf u}{\bm{\sigma}}|{\bf p},K\rangle\right\\}$
$\displaystyle=\frac{i\pi^{2}U_{0}\omega
p_{F}}{2v_{F}|q||\sin\phi|}\Re\left[\left({\bf j}_{K}({\bf
p})e^{i(\phi-\phi^{\prime})/2}+{\bf j}_{K^{\prime}}^{*}({\bf p}+{\bf
q})e^{-i(\phi-\phi^{\prime})/2}\right)\langle{\bf p},K|{\bf
u}^{*}{\bm{\sigma}}|{\bf p}+{\bf q},K^{\prime}\rangle\right]$
$\displaystyle=\frac{i\pi^{2}U_{0}\omega p_{F}}{2v_{F}|q||\sin\phi|}{\bf
F}^{\mathrm{AC}}_{K^{\prime}K}(\alpha,\phi)$ (109)

where the dimensionless vector function ${\bf
F}^{\mathrm{AC}}_{K^{\prime}K}(\alpha,\phi)$ is determined by the matrix
elements of the current density operator and interaction vertices. It depends
only on the directions of the transferred quasi-momentum and the momenta of
the initial and final scattered states, parametrized by the angles
$\alpha,\phi$ as shown in Fig. 4).

##  VIII Perturbative calculation of the DC drag current

The lowest order contributions to the current density of the DC drag are shown
in Fig. 6. Here only the inra-valley current vertices contribute.

Figure 6: Diagram for the lowest order contribution to the DC current.

###  VIII.1 Calculation of the $J_{KK}$ term

The first diagram (with the vertex $j_{KK}$) in Fig. 6 corresponds to the
expression

$\displaystyle\left\langle J_{KK}({\bf
q},\omega)\right\rangle=-iU_{0}^{2}\mathrm{Tr}\left\\{\gamma^{\mathrm{q}}{\bf
J}_{KK}({\bf p},{\bf p})\hat{G}_{K}({\bf p},\epsilon)\gamma^{\mathrm{cl}}({\bf
u}^{*}({\bf q},{\bf p})\cdot{\bm{\sigma}})\hat{G}_{K^{\prime}}({\bf p}+{\bf
q},\epsilon+\omega)\gamma^{\mathrm{cl}}({\bf u}({\bf q},{\bf
p})\cdot{\bm{\sigma}})\hat{G}_{K}({\bf p},\epsilon)\right\\}=$ $\displaystyle-
iU_{0}^{2}\int_{\epsilon,{\bf p}}\mathrm{Tr}\left\\{{\bf J}_{KK}({\bf p},{\bf
p})\left[G_{K}^{K}(\epsilon,{\bf p})({\bf u}^{*}({\bf q},{\bf
p})\cdot{\bm{\sigma}})G_{K^{\prime}}^{A}(\epsilon+\omega,{\bf p}+{\bf q})({\bf
u}({\bf q},{\bf p})\cdot{\bm{\sigma}})G_{K}^{A}(\epsilon,{\bf
p})+\right.\right.$ $\displaystyle G_{K}^{R}(\epsilon,{\bf p})({\bf
u}^{*}({\bf q},{\bf
p})\cdot{\bm{\sigma}})G_{K^{\prime}}^{K}(\epsilon+\omega,{\bf p}+{\bf q})({\bf
u}({\bf q},{\bf p})\cdot{\bm{\sigma}})G_{K}^{A}(\epsilon,{\bf p})+$
$\displaystyle\left.\left.G_{K}^{R}(\epsilon,{\bf p})({\bf u}^{*}({\bf q},{\bf
p})\cdot{\bm{\sigma}})G_{K^{\prime}}^{R}(\epsilon+\omega,{\bf p}+{\bf q})({\bf
u}({\bf q},{\bf p})\cdot{\bm{\sigma}})G_{K}^{K}(\epsilon,{\bf
p})\right]\right\\}.$ (110)

Using the notations Eq. (40), one can represent the expression for the current
$\left\langle J_{KK}\right\rangle$ in the form

$\displaystyle\left\langle J_{KK}\right\rangle=-i\int_{\epsilon,{\bf
p}}\left\\{G^{R}(\epsilon,{\bf p})G^{A}(\epsilon,{\bf
p})\left[\tanh\left(\frac{\epsilon+\omega-\mu}{2T}\right)-\tanh\left(\frac{\epsilon-\mu}{2T}\right)\right]\left[G^{R}(\epsilon+\omega,{\bf
p}+{\bf q})-G^{A}(\epsilon+\omega,{\bf p}+{\bf q})\right]+\right.$
$\displaystyle\left.\tanh\left(\frac{\epsilon-\mu}{2T}\right)\left[\left(G^{R}(\epsilon,{\bf
p})\right)^{2}G^{R}(\epsilon+\omega,{\bf p}+{\bf q})-\left(G^{A}(\epsilon,{\bf
p})\right)^{2}G^{A}(\epsilon+\omega,{\bf p}+{\bf q})\right]\right\\}\times$
$\displaystyle\mathrm{tr}\left\\{{\bf J}_{KK}({\bf p},{\bf
p})\hat{g}_{K}(\epsilon,{\bf p})({\bf u}^{*}({\bf q},{\bf
p})\cdot{\bm{\sigma}})\hat{g}_{K^{\prime}}(\epsilon+\omega,{\bf p}+{\bf
q})({\bf u}({\bf q},{\bf p})\cdot{\bm{\sigma}})\hat{g}_{K}(\epsilon,{\bf
p})\right\\}.$ (111)

In the following calculation we neglect the term in the second line in Eq.
(111) that contains traces ot three retarded or three advanced Green
functions.

The trace in the pseudospin space leads to the sum of the terms proportional
to combinations of different components of the driving potential,

$\mathrm{tr}\left\\{{\bf J}_{KK}({\bf p},{\bf p})\hat{g}_{K}(\epsilon,{\bf
p})({\bf u}^{*}({\bf q},{\bf
p})\cdot{\bm{\sigma}})\hat{g}_{K^{\prime}}(\epsilon+\omega,{\bf p}+{\bf
q})({\bf u}({\bf q},{\bf p})\cdot{\bm{\sigma}})\hat{g}_{K}(\epsilon,{\bf
p})\right\\}=\sum_{\mu,\nu=0}^{3}u_{\mu}u^{*}_{\nu}{\bf
J}_{K}^{\mu\nu}(\epsilon,\omega;{\bf p},{\bf q}),$ (112)

where

${\bf J}_{K}^{\mu\nu}(\epsilon,\omega;{\bf p},{\bf q})=\mathrm{tr}\left\\{{\bf
J}_{KK}({\bf p},{\bf p})\hat{g}_{K}(\epsilon,{\bf
p})\sigma_{\mu}\hat{g}_{K^{\prime}}(\epsilon+\omega,{\bf p}+{\bf
q})\sigma_{\nu}\hat{g}_{K}(\epsilon,{\bf p})\right\\}$ (113)

Substituting explicit expressions for Green functions Eqs. (41), (42), we
obtain

$\displaystyle\left\langle J_{KK}\right\rangle=-\frac{\pi
U_{0}^{2}}{4}\int_{\epsilon,{\bf p}}\sum_{\mu,\nu=0}^{3}u_{\mu}u^{*}_{\nu}{\bf
J}_{K}^{\mu\nu}(\epsilon,\omega;{\bf p},{\bf
q})\left[\tanh\left(\frac{\epsilon+\omega-\mu}{2T}\right)-\tanh\left(\frac{\epsilon-\mu}{2T}\right)\right]\times$
$\displaystyle\left[\delta(\epsilon+\omega-
v_{F}|p+q|)+\delta(\epsilon+\omega+v_{F}|p+q|)\right]\left[\frac{1}{\epsilon-
v_{F}|p|+\frac{i}{2\tau}}+\frac{1}{\epsilon+v_{F}|p|+\frac{i}{2\tau}}\right]\left[\frac{1}{\epsilon-
v_{F}|p|-\frac{i}{2\tau}}+\frac{1}{\epsilon+v_{F}|p|-\frac{i}{2\tau}}\right]\approx$
$\displaystyle-\frac{\pi U_{0}^{2}}{4}\int_{\epsilon,{\bf
p}}\sum_{\mu,\nu=0}^{3}u_{\mu}u^{*}_{\nu}{\bf
J}_{K}^{\mu\nu}(\epsilon,\omega;{\bf p},{\bf
q})\left[\tanh\left(\frac{\epsilon+\omega-\mu}{2T}\right)-\tanh\left(\frac{\epsilon-\mu}{2T}\right)\right]\left[\delta(\epsilon+\omega-
v_{F}|p+q|)+\delta(\epsilon+\omega+v_{F}|p+q|)\right]\times$
$\displaystyle\left[\frac{1}{(\epsilon-
v_{F}|p|)^{2}+\frac{1}{4\tau^{2}}}+\frac{1}{(\epsilon+v_{F}|p|)^{2}+\frac{1}{4\tau^{2}}}\right]\approx$
$\displaystyle-\frac{\pi^{2}\tau U_{0}^{2}}{2}\int_{\epsilon,{\bf
p}}\sum_{\mu,\nu=0}^{3}u_{\mu}u^{*}_{\nu}{\bf
J}_{K}^{\mu\nu}(\epsilon,\omega;{\bf p},{\bf
q})\left[\tanh\left(\frac{\epsilon+\omega-\mu}{2T}\right)-\tanh\left(\frac{\epsilon-\mu}{2T}\right)\right]\times$
$\displaystyle\left[\delta(\epsilon+\omega-
v_{F}|p+q|)+\delta(\epsilon+\omega+v_{F}|p+q|)\right]\left[\delta(\epsilon-
v_{F}|p|)+\delta(\epsilon+v_{F}|p|)\right]\approx$
$\displaystyle-\frac{\pi^{2}\tau U_{0}^{2}}{2}\int_{\bf
p}\sum_{\mu,\nu=0}^{3}u_{\mu}u^{*}_{\nu}{\bf
J}_{K}^{\mu\nu}(\epsilon=v_{F}|p|,\omega;{\bf p},{\bf
q})\left[\tanh\left(\frac{v_{F}|p|+\omega-\mu}{2T}\right)-\tanh\left(\frac{v_{F}|p|-\mu}{2T}\right)\right]\delta(v_{F}|p|+\omega-
v_{F}|p+q|)$ $\displaystyle-\frac{\pi^{2}\tau U_{0}^{2}}{2}\int_{\bf
p}\sum_{\mu,\nu=0}^{3}u_{\mu}u^{*}_{\nu}{\bf
J}_{K}^{\mu\nu}(\epsilon=-v_{F}|p|,\omega;{\bf p},{\bf
q})\left[\tanh\left(\frac{\omega-
v_{F}|p|-\mu}{2T}\right)-\tanh\left(\frac{-v_{F}|p|-\mu}{2T}\right)\right]\delta(\omega-
v_{F}|p|+v_{F}|p+q|)=$ $\displaystyle-\frac{\pi^{2}\tau U_{0}^{2}}{2}\int_{\bf
p}\sum_{\mu,\nu=0}^{3}u_{\mu}u^{*}_{\nu}{\bf
J}_{K}^{\mu\nu}(\epsilon=v_{F}|p|,\omega;{\bf p},{\bf
q})\left[\tanh\left(\frac{v_{F}|p|+\omega-\mu}{2T}\right)-\tanh\left(\frac{v_{F}|p|-\mu}{2T}\right)\right]\delta(v_{F}|p|+\omega-
v_{F}|p+q|).$ (114)

Calculations in Eq. (114) are performed under the assumptions of a large mean-
free time, $1/\tau\ll v_{F}p_{F}$, which allowed replacing the Lorentz-
functions by $\delta$-functions. Furthermore, some $\delta$-functions are
eliminated by the condition $\omega\ll v_{F}|p|$. In the last line we assmed
positive frequency $\omega>0$.

Again, as in the case of AC current, the restrictions imposed by the energy
and momentum conservation together with Fermi-distributions amount to
projection onto the states in the conduction band at low temperatures. It
shows up mathematically in that the main contribution in Eq. (114) for $\mu>0$
is provided by the term with the combination
$\left[\tanh\left(\frac{v_{F}|p|+\omega-\mu}{2T}\right)-\tanh\left(\frac{v_{F}|p|-\mu}{2T}\right)\right]\delta(v_{F}|p+q|-v_{F}|p|-\omega)$.

###  VIII.2 Symmetries by change $K\leftrightarrow K^{\prime}$

The following symmetries between the $K$ and $K^{\prime}$ points emerge in
frame of the rotating wave approximation for the drag potential

$K\rightarrow
K^{\prime}\implies\omega\rightarrow-\omega,\,\,q\rightarrow-\bar{q},\,\,p\rightarrow\bar{p},$
(115)

where we used complex valued vectors $p=p_{x}+ip_{y}$, $q=q_{x}+iq_{y}$, and
$\bar{p},\bar{q}$ denote the complex conjugated vectors. For the angles
$\alpha$ and $\phi$, the transformation Eq. (115) implies

$\alpha\rightarrow\pi-\alpha,\,\,\,\phi\rightarrow\pi-\phi.$ (116)

###  VIII.3 Calculation of the $J_{K^{\prime}K^{\prime}}$ term

The contribution of the diagram with $J_{K^{\prime}K^{\prime}}$ vertex is
obtained from the one for $J_{KK}$ employing the symmetries between the
expressions for $K$ and $K^{\prime}$ valleys $K\rightarrow
K^{\prime}\implies\omega\rightarrow-\omega,\,\,q\rightarrow-\bar{q},\,\,p\rightarrow\bar{p}$.
Below we repeat the main formulas applying the symmetries to the equations for
$J_{KK}$.

$\displaystyle\left\langle J_{K^{\prime}K^{\prime}}({\bf
q},\omega)\right\rangle=-iU_{0}^{2}\mathrm{Tr}\left\\{\gamma^{\mathrm{q}}{\bf
J}_{K^{\prime}K^{\prime}}({\bf p},{\bf p})\hat{G}_{K^{\prime}}({\bf
p},\epsilon)\gamma^{\mathrm{cl}}({\bf u}({\bf q},{\bf
p})\cdot{\bm{\sigma}})\hat{G}_{K}({\bf p}-{\bf
q},\epsilon-\omega)\gamma^{\mathrm{cl}}({\bf u}^{*}({\bf q},{\bf
p})\cdot{\bm{\sigma}})\hat{G}_{K^{\prime}}({\bf p},\epsilon)\right\\}=$
$\displaystyle-iU_{0}^{2}\int_{\epsilon,{\bf p}}\mathrm{Tr}\left\\{{\bf
J}_{K^{\prime}K^{\prime}}({\bf p},{\bf
p})\left[G_{K^{\prime}}^{K}(\epsilon,{\bf p})({\bf u}({\bf q},{\bf
p})\cdot{\bm{\sigma}})G_{K}^{A}(\epsilon-\omega,{\bf p}-{\bf q})({\bf
u}^{*}({\bf q},{\bf p})\cdot{\bm{\sigma}})G_{K^{\prime}}^{A}(\epsilon,{\bf
p})+\right.\right.$ $\displaystyle G_{K^{\prime}}^{R}(\epsilon,{\bf p})({\bf
u}({\bf q},{\bf p})\cdot{\bm{\sigma}})G_{K}^{K}(\epsilon-\omega,{\bf p}-{\bf
q})({\bf u}^{*}({\bf q},{\bf
p})\cdot{\bm{\sigma}})G_{K^{\prime}}^{A}(\epsilon,{\bf p})+$
$\displaystyle\left.\left.G_{K^{\prime}}^{R}(\epsilon,{\bf p})({\bf u}({\bf
q},{\bf p})\cdot{\bm{\sigma}})G_{K}^{R}(\epsilon-\omega,{\bf p}-{\bf q})({\bf
u}^{*}({\bf q},{\bf p})\cdot{\bm{\sigma}})G_{K^{\prime}}^{K}(\epsilon,{\bf
p})\right]\right\\}.$ (117)

Reproducing the calculational steps for the calculation in $K$-valley, we
subsequently obtain

$\displaystyle\left\langle
J_{K^{\prime}K^{\prime}}\right\rangle=-i\int_{\epsilon,{\bf
p}}\left\\{G^{R}(\epsilon,{\bf p})G^{A}(\epsilon,{\bf
p})\left[\tanh\left(\frac{\epsilon-\omega-\mu}{2T}\right)-\tanh\left(\frac{\epsilon-\mu}{2T}\right)\right]\left[G^{R}(\epsilon-\omega,{\bf
p}-{\bf q})-G^{A}(\epsilon-\omega,{\bf p}-{\bf q})\right]+\right.$
$\displaystyle\left.\tanh\left(\frac{\epsilon-\mu}{2T}\right)\left[\left(G^{R}(\epsilon,{\bf
p})\right)^{2}G^{R}(\epsilon-\omega,{\bf p}-{\bf q})-\left(G^{A}(\epsilon,{\bf
p})\right)^{2}G^{A}(\epsilon-\omega,{\bf p}-{\bf q})\right]\right\\}\times$
$\displaystyle\mathrm{tr}\left\\{{\bf J}_{K^{\prime}K^{\prime}}({\bf p},{\bf
p})\hat{g}_{K^{\prime}}(\epsilon,{\bf p})({\bf u}({\bf q},{\bf
p})\cdot{\bm{\sigma}})\hat{g}_{K}(\epsilon-\omega,{\bf p}-{\bf q})({\bf
u}^{*}({\bf q},{\bf p})\cdot{\bm{\sigma}})\hat{g}_{K^{\prime}}(\epsilon,{\bf
p})\right\\}.$ (118)

$\mathrm{tr}\left\\{{\bf J}_{K^{\prime}K^{\prime}}({\bf p},{\bf
p})\hat{g}_{K^{\prime}}(\epsilon,{\bf p})({\bf u}({\bf q},{\bf
p})\cdot{\bm{\sigma}})\hat{g}_{K}(\epsilon-\omega,{\bf p}-{\bf q})({\bf
u}^{*}({\bf q},{\bf p})\cdot{\bm{\sigma}})\hat{g}_{K^{\prime}}(\epsilon,{\bf
p})\right\\}=\sum_{\mu,\nu=0}^{3}u^{*}_{\mu}u_{\nu}{\bf
J}_{K^{\prime}}^{\mu\nu}(\epsilon,\omega;{\bf p},{\bf q}),$ (119)

where

${\bf J}_{K^{\prime}}^{\mu\nu}(\epsilon,\omega;{\bf p},{\bf
q})=\mathrm{tr}\left\\{{\bf J}_{K^{\prime}K^{\prime}}({\bf p},{\bf
p})\hat{g}_{K^{\prime}}(\epsilon,{\bf
p})\sigma_{\mu}\hat{g}_{K}(\epsilon-\omega,{\bf p}-{\bf
q})\sigma_{\nu}\hat{g}_{K^{\prime}}(\epsilon,{\bf p})\right\\}$ (120)

Shifting the integration variables $\epsilon\rightarrow\epsilon+\omega$, ${\bf
p}\rightarrow{\bf p}+{\bf q}$ and neglecting the products of three retarded
and three advanced Green functions, we obtain

$\displaystyle\left\langle
J_{K^{\prime}K^{\prime}}\right\rangle=-i\int_{\epsilon,{\bf
p}}\left\\{G^{R}(\epsilon+\omega,{\bf p+q})G^{A}(\epsilon+\omega,{\bf
p+q})\left[\tanh\left(\frac{\epsilon-\mu}{2T}\right)-\tanh\left(\frac{\epsilon+\omega-\mu}{2T}\right)\right]\left[G^{R}(\epsilon,{\bf
p})-G^{A}(\epsilon,{\bf p})\right]\right\\}\times$
$\displaystyle\mathrm{tr}\left\\{{\bf J}_{K^{\prime}K^{\prime}}({\bf p+q},{\bf
p+q})\hat{g}_{K^{\prime}}(\epsilon+\omega,{\bf p+q})({\bf u}({\bf q},{\bf
p+q})\cdot{\bm{\sigma}})\hat{g}_{K}(\epsilon,{\bf p})({\bf u}^{*}({\bf q},{\bf
p+q})\cdot{\bm{\sigma}})\hat{g}_{K^{\prime}}(\epsilon,{\bf p+q})\right\\}=$
$\displaystyle-\frac{\pi U_{0}^{2}}{4}\int_{\epsilon,{\bf
p}}\sum_{\mu,\nu=0}^{3}u_{\mu}u^{*}_{\nu}{\bf
J}_{K^{\prime}}^{\mu\nu}(\epsilon+\omega,\omega;{\bf p+q},{\bf
q})\left[\tanh\left(\frac{\epsilon-\mu}{2T}\right)-\tanh\left(\frac{\epsilon+\omega-\mu}{2T}\right)\right]\left[\delta(\epsilon-
v_{F}|p|)+\delta(\epsilon+v_{F}|p|)\right]\times$
$\displaystyle\left[\frac{1}{\epsilon+\omega-v_{F}|{\bf
p+q}|+\frac{i}{2\tau}}+\frac{1}{\epsilon+\omega+v_{F}|{\bf
p+q}|+\frac{i}{2\tau}}\right]\left[\frac{1}{\epsilon+\omega-v_{F}|{\bf
p+q}|-\frac{i}{2\tau}}+\frac{1}{\epsilon+\omega+v_{F}|{\bf
p+q}|-\frac{i}{2\tau}}\right]\approx$ $\displaystyle-\frac{\pi
U_{0}^{2}}{4}\int_{\epsilon,{\bf p}}\sum_{\mu,\nu=0}^{3}u_{\mu}u^{*}_{\nu}{\bf
J}_{K^{\prime}}^{\mu\nu}(\epsilon+\omega,\omega;{\bf p+q},{\bf
q})\left[\tanh\left(\frac{\epsilon-\mu}{2T}\right)-\tanh\left(\frac{\epsilon+\omega-\mu}{2T}\right)\right]\left[\delta(\epsilon-
v_{F}|p|)+\delta(\epsilon+v_{F}|p|)\right]\times$
$\displaystyle\left[\frac{1}{(\epsilon+\omega-v_{F}|{\bf
p+q}|)^{2}+\frac{1}{4\tau^{2}}}+\frac{1}{(\epsilon+\omega+v_{F}|{\bf
p+q}|)^{2}+\frac{1}{4\tau^{2}}}\right]\approx$
$\displaystyle-\frac{\pi^{2}\tau U_{0}^{2}}{2}\int_{\epsilon,{\bf
p}}\sum_{\mu,\nu=0}^{3}u_{\mu}u^{*}_{\nu}{\bf
J}_{K^{\prime}}^{\mu\nu}(\epsilon+\omega,\omega;{\bf p+q},{\bf
q})\left[\tanh\left(\frac{\epsilon-\mu}{2T}\right)-\tanh\left(\frac{\epsilon+\omega-\mu}{2T}\right)\right]\left[\delta(\epsilon-
v_{F}|p|)+\delta(\epsilon+v_{F}|p|)\right]\times$
$\displaystyle\left[\delta(\epsilon+\omega-v_{F}|{\bf
p+q}|)+\delta(\epsilon+\omega+v_{F}|{\bf p+q}|)\right]\approx$
$\displaystyle-\frac{\pi^{2}\tau U_{0}^{2}}{2}\int_{{\bf
p}}\sum_{\mu,\nu=0}^{3}u_{\mu}u^{*}_{\nu}{\bf
J}_{K^{\prime}}^{\mu\nu}(v_{F}|p|+\omega,\omega;{\bf p+q},{\bf
q})\left[\tanh\left(\frac{v_{F}|p|-\mu}{2T}\right)-\tanh\left(\frac{v_{F}|p|+\omega-\mu}{2T}\right)\right]\delta(v_{F}|p|+\omega-
v_{F}|{\bf p+q}|)-$ $\displaystyle\frac{\pi^{2}\tau U_{0}^{2}}{2}\int_{{\bf
p}}\sum_{\mu,\nu=0}^{3}u_{\mu}u^{*}_{\nu}{\bf
J}_{K^{\prime}}^{\mu\nu}(-v_{F}|p|+\omega,\omega;{\bf p+q},{\bf
q})\left[\tanh\left(\frac{-v_{F}|p|-\mu}{2T}\right)-\tanh\left(\frac{-v_{F}|p|+\omega-\mu}{2T}\right)\right]\delta(\omega+v_{F}|{\bf
p+q}|-v_{F}|p|)=$ $\displaystyle-\frac{\pi^{2}\tau U_{0}^{2}}{2}\int_{{\bf
p}}\sum_{\mu,\nu=0}^{3}u_{\mu}u^{*}_{\nu}{\bf
J}_{K^{\prime}}^{\mu\nu}(v_{F}|p|+\omega,\omega;{\bf p+q},{\bf
q})\left[\tanh\left(\frac{v_{F}|p|-\mu}{2T}\right)-\tanh\left(\frac{v_{F}|p|+\omega-\mu}{2T}\right)\right]\delta(v_{F}|p|+\omega-
v_{F}|{\bf p+q}|)$ (121)

###  VIII.4 Total DC drag current

The total DC drag current density, given by the sum of Eqs. (114) and (121),
can be represented in the form

$\displaystyle{\bf j}^{\mathrm{DC}}=-\frac{\pi^{2}\tau U_{0}^{2}}{2}\int_{\bf
p}\sum_{\mu,\nu=0}^{3}\left[u_{\mu}({\bf p},{\bf q})u^{*}_{\nu}({\bf p},{\bf
q}){\bf J}_{K}^{\mu\nu}(\epsilon=v_{F}|p|,\omega;{\bf p},{\bf q})-u_{\mu}({\bf
p+q},{\bf q})u^{*}_{\nu}({\bf p+q},{\bf q}){\bf
J}_{K^{\prime}}^{\mu\nu}(\epsilon=v_{F}|p|+\omega,\omega;{\bf p+q},{\bf
q})\right]\times$
$\displaystyle\left[\tanh\left(\frac{v_{F}|p|+\omega-\mu}{2T}\right)-\tanh\left(\frac{v_{F}|p|-\mu}{2T}\right)\right]\delta(v_{F}|p|+\omega-
v_{F}|p+q|).$ (122)

###  VIII.5 Projection on the conduction band

The matrix elements of the current density operators between the states in
conduction band are given by

$\displaystyle\langle{\bf p},K|{\bf J}_{KK}({\bf p},{\bf p})|{\bf
p},K\rangle_{c}=\frac{1}{2}\left(e^{i\phi/2},e^{-i\phi/2}\right)\left(\begin{array}[]{cc}0&{\bf
j}_{K}({\bf p})\\\ {\bf j}_{K}^{*}({\bf
p})&0\end{array}\right)\left(\begin{array}[]{c}e^{-i\phi/2}\\\
e^{i\phi/2}\end{array}\right)=$ (127) $\displaystyle\frac{1}{2}\left\\{{\bf
j}_{K}({\bf p})e^{i\phi}+{\bf j}_{K}^{*}({\bf p})e^{-i\phi}\right\\},$ (128)

where $\phi=\phi_{\bf p}$. Analogously, for the valley $K^{\prime}$ we obtain

$\displaystyle\langle{\bf p^{\prime}},K^{\prime}|{\bf
J}_{K^{\prime}K^{\prime}}({\bf p^{\prime}},{\bf p^{\prime}})|{\bf
p^{\prime}},K^{\prime}\rangle_{c}=\frac{1}{2}\left(e^{-i\phi^{\prime}/2},e^{i\phi^{\prime}/2}\right)\left(\begin{array}[]{cc}0&{\bf
j}_{K^{\prime}}({\bf p^{\prime}})\\\ {\bf j}_{K^{\prime}}^{*}({\bf
p^{\prime}})&0\end{array}\right)\left(\begin{array}[]{c}e^{i\phi^{\prime}/2}\\\
e^{-i\phi^{\prime}/2}\end{array}\right)=$ (133)
$\displaystyle\frac{1}{2}\left\\{{\bf j}_{K^{\prime}}({\bf
p^{\prime}})e^{-i\phi^{\prime}}+{\bf j}_{K^{\prime}}^{*}({\bf
p^{\prime}})e^{i\phi^{\prime}}\right\\},$ (134)

where $\phi^{\prime}=\phi_{{\bf p}+{\bf q}}$ as given by Eqs. (64), (65).

Upon projection on the conduction band, the leading contributions to the DC
drag current density are given by

$\displaystyle\left\langle
J_{KK}\right\rangle\approx\frac{\pi^{2}\tau}{2}U_{0}^{2}\int_{{\bf
p}}\sum_{\mu,\nu=0}^{3}u_{\mu}u^{*}_{\nu}{\bf
J}_{K}^{\mu\nu}(\epsilon=v_{F}|p|,\omega;{\bf p},{\bf
q})\left[\tanh\left(\frac{v_{F}|p|+\omega-\mu}{2T}\right)-\tanh\left(\frac{v_{F}|p|-\mu}{2T}\right)\right]\delta(v_{F}|p+q|-v_{F}|p|-\omega)\approx$
$\displaystyle\frac{\pi^{2}U_{0}^{2}\tau\omega
p_{F}}{v_{F}|q||\sin\phi|}\langle{\bf p},K|{\bf J}_{KK}({\bf p},{\bf p})|{\bf
p},K\rangle_{c}\langle{\bf p},K|({\bf u}^{*}\cdot{\bm{\sigma}})|{\bf p}+{\bf
q},K^{\prime}\rangle_{c}\langle{\bf p}+{\bf q},K^{\prime}|({\bf
u}\cdot{\bm{\sigma}})|{\bf p},K\rangle_{c},$ (135) $\displaystyle\left\langle
J_{K^{\prime}K^{\prime}}\right\rangle\approx-\frac{\pi^{2}\tau
U_{0}^{2}}{2}\int_{{\bf p}}\sum_{\mu,\nu=0}^{3}u_{\mu}u^{*}_{\nu}{\bf
J}_{K^{\prime}}^{\mu\nu}(v_{F}|p|+\omega,\omega;{\bf p+q},{\bf
q})\left[\tanh\left(\frac{v_{F}|p|-\mu}{2T}\right)-\tanh\left(\frac{v_{F}|p|+\omega-\mu}{2T}\right)\right]\delta(v_{F}|p|+\omega-
v_{F}|{\bf p+q}|)\approx$ $\displaystyle\frac{\pi^{2}U_{0}^{2}\tau\omega
p_{F}}{v_{F}|q||\sin\phi|}\langle{\bf p+q},K^{\prime}|{\bf
J}_{K^{\prime}K^{\prime}}({\bf p+q},{\bf p+q})|{\bf
p+q},K^{\prime}\rangle_{c}\langle{\bf p},K|({\bf
u}^{*}\cdot{\bm{\sigma}})|{\bf p}+{\bf q},K^{\prime}\rangle_{c}\langle{\bf
p}+{\bf q},K^{\prime}|({\bf u}\cdot{\bm{\sigma}})|{\bf p},K\rangle_{c}$ (136)

Close to the threshold doping, at the angles $\phi$ such that
$|\sin\phi|\ll\frac{1}{\sqrt{v_{F}|q|\tau}}$, the divergent factor
$1/|\sin\phi|$ has to be replaced by the cutoff $\sqrt{v_{F}|q|\tau}$ due to
the intrinsic scattering processes in graphene.

###  VIII.6 Drag current at zero temperature

Projecting Eq. (122) on the states of conduction band, be obtain

$\displaystyle{\bf j}^{\mathrm{DC}}=\frac{\pi^{2}U_{0}^{2}\tau\omega
p_{F}}{v_{F}^{2}|q||\sin\phi|}\left[\langle{\bf p},K|{\bf J}_{KK}({\bf p},{\bf
p})|{\bf p},K\rangle_{c}-\langle{\bf p+q},K^{\prime}|{\bf
J}_{K^{\prime}K^{\prime}}({\bf p+q},{\bf p+q})|{\bf
p+q},K^{\prime}\rangle_{c}\right]\times$ $\displaystyle\langle{\bf p},K|({\bf
u}^{*}\cdot{\bm{\sigma}})|{\bf p}+{\bf q},K^{\prime}\rangle_{c}\langle{\bf
p}+{\bf q},K^{\prime}|({\bf u}\cdot{\bm{\sigma}})|{\bf p},K\rangle_{c}$ (137)

Finally, keeping the contribution $\left\langle J_{KK}\right\rangle$ only,
restoring the factors of $\hbar$, and expressing the life-time $\tau$ through
the Drude conductivity of graphene
$\rho_{g}=\frac{\pi\hbar^{2}}{e^{2}\mu\tau}$, where $\mu$ denotes the chemical
potential, we obtain the dc drag current at zero temperature in the form

${\bf j}^{\mathrm{DC}}=\frac{\pi^{3}U_{0}^{2}({\bf v}_{0}\cdot{\bf
Q})}{e\hbar\rho_{g}v_{F}^{2}|q|}\frac{{\bf
F}^{\mathrm{DC}}(\alpha,\phi)}{|\sin\phi|}.$ (138)

##  IX Berry curvature

The emergence of the Hall component of the drag current is analogous to the
nonlinear Hall effect explored in Refs. [19, 20]. For a spatially homogeneous
periodic electric field, the relationship between the DC drag current and the
Berry curvature is given by the formula (the elementary charge is set to 1)

${\bf
j}_{a}^{\mathrm{DC}}=\epsilon_{adc}\mathcal{E}_{b}\mathcal{E}^{*}_{c}\frac{\tau}{1+i\omega\tau}\int_{\bf
p}\left(\frac{\partial f_{0}({\bf p})}{\partial p_{b}}\right)\Omega_{d}.$
(139)

For the inter-valley nonlinear Hall effect, Eq. (139) is modified due to
additional constraints imposed by momentum and energy conservation
requirements in intervalley scattering

${\bf
j}_{a}^{\mathrm{DC}}=\epsilon_{adc}\mathcal{E}_{b}\mathcal{E}^{*}_{c}\tau\int_{\bf
p}\left(\frac{\partial f_{0}({\bf p})}{\partial p_{b}}\right)\delta(|{\bf
p+q|}-|{\bf
p}|-\omega)\Omega_{d}\approx\epsilon_{adc}\mathcal{E}_{b}\mathcal{E}^{*}_{c}\frac{\tau\omega}{v_{F}}\int_{\bf
p}\frac{p_{b}}{p_{F}}\delta(|{\bf p}|-p_{F})\delta(|{\bf
p+q|}-p_{F})\Omega_{d}$ (140)

where in the last equation we consider the limit $\omega\tau\ll 1$. Here
$f_{0}({\bf p})$ denotes the Fermi distribution, $\mathcal{E}$ denotes the
comlpex-valued amplitude of the electric field, which is related to the
periodic electric field by ${\bf E}=\mathcal{E}e^{i\omega t}$, and the indexes
$a,b,c,d$ take possible vallues $x,y,z$. $\Omega_{d}$ denotes the
$d$-component of the Berry curvature. In a two-dimensional system, the only
component of the Berry curvature entering Eq. (139) is $\Omega_{z}$. In the
case of Coulomb drag by a moving polar crystal, the external elelectric field
is determined by the gradient of the potential created by the moving lattice.
Transforming the drag Hamiltonian Eq. (55) back to coordinate space, we obtain

$H_{d}=\frac{U_{0}}{2}\int d{\bf R}d{\bf r}\Psi^{+}_{K^{\prime}}({\bf
R})\left({\bf u}({\bf R},{\bf r})\cdot\bm{\sigma}\right)\Psi_{K}({\bf R}-{\bf
r})e^{-i\omega t}+\mathrm{h.c.}$ (141)

Here the potential ${\bf u}({\bf R},{\bf r})$ is given by

${\bf u}({\bf R},{\bf r})=e^{i{\bf qR}}\int_{\bf p}{\bf u}({\bf p},{\bf
q})e^{i{\bf pr}}.$ (142)

Explicit expressions for ${\bf u}({\bf p},{\bf q})$ as given by Eqs. (56) –
(59) show that the components $u_{0}$ and $u_{3}$ are local, proportional to
$\delta({\bf r})$, while the components $u_{1}$ and $u_{2}$ are short-ranged
in ${\bf r}$. Therefore, for the determination of the electric field we
consider the potential as completely local and use the definition

${\mathcal{\bm{E}}}({\bf R})=-U_{0}\left(\nabla_{\bf R}{\bf u}({\bf R},{\bf
r})\cdot\bm{\sigma}\right)=-i{\bf q}e^{i{\bf qR}}\int_{\bf p}\left({\bf
u}({\bf p},{\bf q})\cdot\bm{\sigma}\right)e^{i{\bf pr}}$ (143)

Perfoming the Fourier transform with respect to ${\bf R},{\bf r}$, we obtain
the electric field as

${\mathcal{\bm{E}}}({\bf p},{\bf q})=-iU_{0}{\bf q}\left({\bf u}({\bf p},{\bf
q})\cdot\bm{\sigma}\right).$ (144)

For the electric field pointing in the direction of ${\bf q}$, Eq. (140)
relates the Berry curvature to the Hall component to the current, which is
perpendicular to ${\bf q}$

$j_{\mathrm{H}}^{\mathrm{DC}}=|\mathcal{E}_{q}|^{2}\frac{\tau\omega
p_{F}}{v_{F}q}\frac{\cos\phi}{|\sin\phi|}\Omega_{z},$ (145)

where $\phi$ denotes the angle between $\bf p$ and $\bf q$.

Let us now relate the Berry curvature to the component of the drag current
perpendicular to the applied field (the Hall component). Note, that according
to Eq. (137), the total drag current can be represented as a difference of the
partial currents associated to the valleys $K$ and $K^{\prime}$ as follows

${\bf j}^{\mathrm{DC}}={\bf j}^{\mathrm{DC}}_{K}-{\bf
j}^{\mathrm{DC}}_{K^{\prime}},$ (146)

where the partial currents are determined using Eqs. (128), (134)

$\displaystyle{\bf j}^{\mathrm{DC}}_{K}=\frac{\pi^{2}\tau\omega
p_{F}}{2v_{F}^{2}|q|^{3}|\sin\phi|}|{\mathcal{E}}_{q}|^{2}\left\\{{\bf
j}_{K}({\bf p})e^{i\phi_{\bf p}}+{\bf j}_{K}^{*}({\bf p})e^{-i\phi_{\bf
p}}\right\\},$ (147) $\displaystyle{\bf
j}^{\mathrm{DC}}_{K^{\prime}}=\frac{\pi^{2}\tau\omega
p_{F}}{2v_{F}^{2}|q|^{3}|\sin\phi|}|{\mathcal{E}}_{q}|^{2}\left\\{{\bf
j}_{K^{\prime}}({\bf p}+{\bf q})e^{i\phi_{\bf p+q}}+{\bf
j}_{K^{\prime}}^{*}({\bf p}+{\bf q})e^{-i\phi_{\bf p+q}}\right\\}.$ (148)

Comparing Eqs. (147), (148) with Eq. (145), we relate the Hall component of
the matrix element of the partial current to the Berry phase as follows

$\Omega_{z}=\frac{\pi^{2}\Re\left[j^{\mathrm{DC}}_{H}e^{i\phi_{\bf
p}}\right]}{v_{F}q^{2}\cos\phi}$ (149)

Note, that the angles $\phi$ and $\phi_{\bf p}$ are different, namely, $\phi$
denotes the angle between the vectors ${\bf p}$ and $\bf q$, whereas
$\phi_{\bf p}$ denotes the direction of ${\bf p}$. The relations between the
angles is are given by Eqs. (64), (65). Using Eqs. (64), (65) (see also Fig.
4), one can represent the Hall component of the current as

$j^{\mathrm{DC}}_{H}=|\mathcal{E}|^{2}\frac{\tau\omega
p_{F}}{v_{F}q}\cot\phi\left[\Omega^{K}_{z}({\bf p})+\Omega^{K}_{z}({\bf-
p-q})+\Omega^{K^{\prime}}_{z}({\bf p+q})+\Omega^{K^{\prime}}_{z}(-{\bf
p})\right].$ (150)

Therefore, the Hall component of the drag current measures the sum of the
Berry phases in $K$ and $K^{\prime}$ valleys , which we denote as

$\bar{\Omega}({\bf p},{\bf q})=\Omega^{K}_{z}({\bf p})+\Omega^{K}_{z}({\bf-
p-q})+\Omega^{K^{\prime}}_{z}({\bf p+q})+\Omega^{K^{\prime}}_{z}(-{\bf p})$
(151)

The sum of the Berry phases Eq. (151) can be calculated from the relation of
the longitudinal (with respect to ${\bf q}$ ) $j_{q}$ and Hall, $j_{H}$
components of the drag current

$\bar{\Omega}({\bf p},{\bf
q})=\frac{\pi^{2}}{v_{F}q^{2}\cos\phi}\frac{j^{\mathrm{DC}}_{H}}{j^{\mathrm{DC}}_{q}}\Re\left\\{j_{q,K}({\bf
p})e^{i\phi_{\bf p}}+j_{q,K^{\prime}}({\bf p+q})e^{i\phi_{\bf p+q}}\right\\}$
(152)

Here the components of the drag current $j^{\mathrm{DC}}_{H}$,
$j^{\mathrm{DC}}_{q}$ are given by experimental measurements, whereas the
matrix elements $j_{q,K}({\bf p})$, $j_{q,K^{\prime}}({\bf p+q})$ are
calculated by projecting Eqs. (72), (80) on the direction of ${\bf q}$.

##  X Estimations of the threshold doping level

To estimate the threshold doping level, we use the following numerical values:
the absolute value of the inverse lattice vector of the square ice
$Q=\frac{2\pi}{a_{0}}=\frac{2\pi}{2.8}10^{10}$m${}^{-1}\approx 2.24\cdot
10^{10}$m-1; the distance between the $K$ and $K^{\prime}$ points $|{\bf
K}-{\bf K^{\prime}}|=\frac{4\pi}{3\sqrt{3}a}=\frac{4\pi}{3\sqrt{3}1.42}\cdot
10^{10}$m${}^{-1}\approx 1.7\cdot 10^{10}$m-1. The threshold chemical
potential corresponds to the Fermi wave vector $p_{F}\approx q/2$. For the
parallel orientation of ${\bf q}$ and ${\bf K}-{\bf K^{\prime}}$, we obtain
$q=Q-|{\bf K}-{\bf K^{\prime}}|\approx 0.54\cdot 10^{10}$m-1, which
corresponds to the Fermi energy $E_{F}=\hbar v_{F}q/2\approx 1.7$eV. This
value is below the energy of the van Hove singularity at the M-point
$E_{vH}\approx 2$eV.