# Plasmonic drag photocurrent in graphene at extreme nonlocality

Vladimir Silkin  Dmitry Svintsov  Laboratory of 2d Materials for
Optoelectronics, Moscow Institute of Physics and Technology, Dolgoprudny
141700, Russia

###### Abstract

It is commonly assumed that photocurrent in two-dimensional systems with
centrosymmetric lattice is generated at structural inhomogenities, such as p-n
junctions. Here, we study an alternative mechanism of photocurrent generation
associated with inhomogenity of the driving electromagnetic field, termed as
’plasmonic drag’. It is associated with direct momentum transfer from field to
conduction electrons, and can be characterized by a non-local non-linear
conductivity $\sigma^{(2)}({\bf q},\omega)$. By constructing a classical
kinetic model of non-linear conductivity with full account of non-locality, we
show that it is resonantly enhanced for wave phase velocity coinciding with
electron Fermi velocity. The enhancement is interpreted as phase locking
between electrons and the wave. We discuss a possible experiment where non-
uniform field is created by a propagating graphene plasmon, and find an upper
limit of the current responsivity vs plasmon velocity. This limit is set by a
competition between resonantly growing $\sigma^{(2)}({\bf q},\omega)$ and
diverging kinetic energy of electrons as the wave velocity approaches Fermi
velocity.

##  I Introduction

It is commonly believed that photocurrent generation in two-dimensional
systems without lattice inversion asymmetry occurs at the structural
inhomogenities, such as p-n junctions Song  _et al._ (2011); Muravev and
Kukushkin (2012); Tielrooij  _et al._ (2015) and contacts with metals Cai  _et
al._ (2014); Echtermeyer  _et al._ (2014). At the same time, the inversion
asymmetry of driving electromagnetic field itself can lead to emergence of
photocurrent. The physics beyond such photocurrent is the direct transfer of
electromagnetic field momentum to the electrons, termed in literature as
photon drag Gibson  _et al._ (1970), or dynamic Hall effect Karch  _et al._
(2010). As the optoelectronic structures become deeply sub-wavelength, their
electromagnetic response is dominated by near-fields. The mechanism of current
generation by field with finite momentum is thus termed as plasmonic drag
Vengurlekar and Ishihara (2005); Durach and Noginova (2016); Popov (2013),
implying that momentum is transferred via excitation of two-dimensional
plasmons.

Plasmonic drag photocurrent is ubiquitous to systems with uniform channel and
asymmetric electromagnetic environment. Experimental examples of such
environments are gate edges Song  _et al._ (2011), metal contacts, and grating
couplers without inversion symmetry Olbrich  _et al._ (2016); Popov  _et al._
(2015) (Fig. 1). Plasmonic drag can occur upon tilted illumination of fully
symmetric structure with identical contacts, or illumination of one of its
contacts which acts as photon-to-plasmon coupler Nikulin  _et al._ (2021).
Compared to photovoltaic and photo-thermoelectric effects, the plasmonic drag
does not need channel doping non-uniformity. We may therefore speculate that
it is the most omnipresent mechanism of photocurrent generation in two
dimensions 111The plasmonic drag effect is conceptually equivalent to the
distributed resistive self-mixing Sakowicz  _et al._ (2011) and Dyakonov-Shur
Sakowicz  _et al._ (2011) rectification. This can be seen by presenting the
photocurrent via linear-response electric field and its gradients, all three
mechanisms will be described by identical expressions..

Figure 1: Possible structures where plasmonic drag can be observed upon
diffraction of incident plane electromagnetic wave: (a) graphene under the
edge of a metal gate; (b) graphene contacted by metal (c) graphene under the
grating lacking inversion symmetry

As the size of the optoelectronic devices $L$ is shrinking, the characteristic
wave number $q$ of electromagnetic near fields scales up as $L^{-1}$. Down-
scaling of the gate-channel separation also leads to slow-down of plasmons
Chaplik (1972) and increase in characteristic wave vector at given frequency.
At very large wave vectors, the current-field relations become non-local and
reflect the dynamics of individual electrons. While linear non-local response
of two-dimensional systems is quite well studied Nehls  _et al._ (1996);
Bandurin  _et al._ (2021); Lundeberg  _et al._ (2017), the non-linear
processes of photocurrent generation at strong non-locality are almost
unexplored. Indeed, previous studies of plasmonic drag were limited either to
expansions linear in field wave vector Olbrich  _et al._ (2016), or to the
hydrodynamic description of electron systems Popov (2013). The latter is valid
only at small frequencies compared to that of carrier-carrier collisions, and
at long wavelengths much exceeding the mean free path Svintsov (2018).

In this Letter, we theoretically explore the limits of plasmonic drag in
graphene at arbitrarily strong non-locality, i.e. at arbitrarily large field
wave vector $q$. Our choice of graphene is dictated by recent observations of
ultra-confined acoustic graphene plasmons Lundeberg  _et al._ (2017); Iranzo
_et al._ (2018); Epstein  _et al._ (2020); Bylinkin  _et al._ (2019) which
spectrum is governed by non-locality, and their potential for photocurrent
harvesting Bandurin  _et al._ (2018a). We derive the non-local non-linear
conductivity of graphene $\sigma^{(2)({\bf q},\omega)}$ being the
proportionality coefficient between photocurrent density and squared ac
electric field, ${\bf j}^{(2)}={\bf n}_{\bf q}\sigma^{(2)({\bf
q},\omega)}E^{2}$. It possesses a square-root singularity at phase velocity
$\omega/q$ approaching the Fermi velocity of 2d electrons $v_{0}$. We
interpret this singularity as phase locking between dragged electrons and
electromagnetic field. Second, we quantify the electromagnetic energy flux
$S_{\rm pl}$ carried by 2d plasmons at large $q$. The ratio of these
quantities is nothing but current responsivity measured in photodetetion
experiments. The energy flux has a counter-balancing singularity to that in
$\sigma^{(2)}$, it appears due to large kinetic energy of charge carriers in
electromagnetic wave. As a result, the photocurrent responsivity has a
universal maximum order of $0.25e/E_{F}$ achieved at $\omega\approx
1.4qv_{0}$.

The starting point for evaluation of photocurrent in non-uniform field ${\bf
E}({\bf r},t)$ is the classical kinetic equation for electron distribution
function $f({\bf r},{\bf p},t)$:

$\frac{\partial f}{\partial t}+{\bf v_{p}}\frac{\partial f}{\partial{\bf
r}}-e{\bf E}({\bf r},t)\frac{\partial f}{\partial{\bf
p}}=-\frac{f-f_{0}}{\tau_{p}},$ (1)

where ${\bf v}_{\bf p}$ is the electron velocity, $\tau_{p}$ is the momentum
relaxation time, and $f_{0}$ is the equilibrium (Fermi) distribution function.
We specify the field in the form of running wave ${\bf E}({\bf r},t)={\bf
E}_{0}\>e^{i({\bf q}{\bf r}-\omega t)}+{\rm h.c.}$, where ${\rm h.c.}$ stands
for complex conjugate, and limit ourselves to purely longitudinal fields ${\bf
E}_{0}\parallel{\bf q}$. We adopt successive approximations to distribution
function in powers of the field, $f=f_{0}+f_{1}+f_{2}$, where $f_{n}\propto
E_{0}^{n}$.

Figure 2: Color map of non-linear conductivity $\sigma^{(2)}(q,\omega)$ vs
frequency $\omega/2\pi$ and wave vector $q$. Vertical line (a) corresponds to
constant wave vector $q=q_{F}$; line (b) shows the dispersion of graphene
plasmons. Momentum relaxation time $\tau_{p}=10^{-13}$ s, carrier density
$n=10^{10}$ cm-2.

The quantity of interest for calculation of dc photocurrent is the to time-
and space-averaged second-order distribution function $\langle
f_{2}\rangle_{t,x}$. Performing lengthy yet common calculations, we find it in
the form

${\langle
f_{2}\rangle}_{t,x}=\frac{e^{2}E_{0}^{2}}{2}\frac{\partial}{\partial{\bf
p}}\Bigg{\\{}\frac{\partial f_{0}/\partial{\bf p}}{(\omega-{\bf
qv_{p}})^{2}+{\tau_{p}}^{-2}}\Bigg{\\}}.$ (2)

The rectified current is obtained by integration over momentum space ${\bf
j}^{(2)}=g\sum\limits_{\bf p}{{\langle f_{2}\rangle}_{t,x}}$, where $g=4$ is
the electron degeneracy factor in graphene. Apparently, the rectified current
is proportional to $E_{0}^{2}$ and directed along ${\bf n_{q}}={\bf q}/q$.
Introducing the nonlinear conductivity ${\bf
j}^{(2)}=\sigma^{(2)}(q,\omega)E_{0}^{2}{\bf n_{q}}$, we find (see Appendix A
for details):

$\sigma^{(2)}(q,\omega)=-\frac{ge^{3}}{2\pi\hbar^{2}v_{0}}\frac{f_{0}(0)}{q^{2}}\left(\frac{2\,\widetilde{v}_{ph}^{2}-1}{2\,\sqrt{\widetilde{v}_{ph}^{2}-1}}-\widetilde{v}_{ph},\right)$
(3)

where $\widetilde{v}_{ph}=(\omega+i\tau_{p}^{-1})/(qv_{0})$ is the
dimensionless phase velocity, and $f_{0}(0)$ is the equilibrium electron
distribution function evaluated at zero energy. The photocurrent due to drag
of holes is additive with opposite sign, thus account of both sorts of
carriers amounts to replacement $f_{0}(0)\rightarrow 2f_{0}(0)-1$.

The most remarkable property of nonlinear conductivity (3) is the presence of
square-root singularity as the wave phase velocity $\omega/q$ approaches the
electron Fermi velocity from either side, as shown in Fig. 2. A detailed
inspection shows that singular contribution to current comes from electrons
moving in phase with the wave, i.e. at angles
$\cos\theta\approx\omega/(qv_{0})$. We may say that these electrons are
trapped in the minima of harmonic potential induced by the wave, and move
synchronously with the wave velocity.

The effect of finite scattering rate $\tau_{p}^{-1}$ on nonlinear conductivity
is highly non-trivial. It depends critically on the ratio of phase velocity
$\omega/q$ and Fermi velocity $v_{0}$, which is illustrated in Fig. 3 (a). For
slow waves with $\omega/qv_{0}<1$, the non-linear conductivity $\sigma^{(2)}$
is approximately proportional to relaxation time. This result is interpreted
as follows: the amount of momentum transferred from field to electrons is
proportional to the rate of Landau damping. Being a collisionless process, it
does not depend on $\tau_{p}$. The current established at given momentum
transfer rate (i.e. at given force density) is inversely proportional to
electron scattering rate.

Figure 3:  Dependence of the nonlinear conductivity $\sigma^{(2)}$ on
normalized phase velocity at (a) constant wave vector $q=q_{F}=20$ $\mu$m-1,
corresponding to the line cut (a) in Fig. 2 (b) $\omega$ and $q$ bound by
plasmon dispersion at various carrier densities $n=10^{10}$, $10^{11}$ and
$10^{12}$ cm-2. Momentum relaxation time $\tau_{p}=10^{-13}$ s.

The situation for fast waves, $\omega/qv_{0}>1$, is different. Such waves can
not induce intraband Landau damping by the virtue of momentum conservation,
and finite scattering rate is required to soften the momentum constraint. At
the same time, scattering acts to dissipate the generated current. As a
result, $\sigma^{(2)}$ becomes independent of scattering time for relatively
fast phase velocities $\omega/qv_{0}\gg 1$. This rule breaks down as we
approach the singularity at $\omega=qv_{0}$ which is softened by scattering at
both sides.

The above description of rectified photocurrent was performed at given field
${\bf E}_{0}({\bf q},\omega)$ in the plane of graphene. Determination of this
field requires a solution of linear diffraction problem for a wave scattered
by graphene and its electromagnetic environment (gratings Fateev  _et al._
(2017), tips McLeod  _et al._ (2014), contacts Nikulin  _et al._ (2021),
etc.). According to the results of analytical and numerical studies, the field
in the plane of graphene is dominated by 2d plasmons. Their wave vector is
bound to frequency via dispersion relation Ryzhii (2006):

$\omega_{\rm pl}(q)=qv_{0}\frac{1+g\alpha q_{F}/q}{\sqrt{1+2g\alpha
q_{F}/q}},$ (4)

where $\alpha=e^{2}/(\kappa\hbar v_{0})$ is the coupling constant in graphene,
$\kappa$ is the background dielectric constant, and $q_{F}$ is the electron
Fermi wave vector.

A distinctive property of dispersion (4) lies in complete absence of Landau
damping. This is guaranteed by the phase velocity of graphene plasmon that
always exceeds Fermi velocity $v_{0}$. It is thus tempting to see whether the
the singularity in non-linear conductivity, Eq. (3), can be probed with
graphene plasmons.

To test this, we bind $q$ and $\omega$ via plasmon dispersion and plot
$\sigma^{(2)}$ parametrically via $v_{\rm ph}$ (which now cannot fall below
unity). The result is shown in Fig. 3 (b). Once the phase velocity is large,
$v_{\rm ph}\gg 1$, the nonlinear conductivity rolls down. Indeed, large phase
velocities correspond to initial part of plasmon dispersion curve where the
characteristic wave vector (transferred momentum) is small. Once the plasmon
phase velocity approaches Fermi velocity, $\sigma^{(2)}$ becomes small again.
This fact may look inconsistent with the presence of singularity in
$\sigma^{(2)}$, but it is. Indeed, synchronism of wave velocity and Fermi
velocity is achieved only at very high frequencies, where inertial electrons
do not keep up with rapidly oscillating electric field.

The only way to probe the singular non-linear conductivity is thus to bring
the whole plasmon dispersion closer to the singular line $\omega=qv_{0}$. This
can be achieved via reduction of carrier density and/or increase in background
dielectric constant. The result is illustrated in Fig. 3 (b), where reducing
$n=k_{F}^{2}/\pi$ from $10^{12}$ cm-2 to $10^{10}$ cm-2 leads to an order of
magnitude enhancement of non-linear conductivity.

Figure 4: Dependence of the current responsivity $R_{J}$ on frequency
$\omega/2\pi$ and wave vector $q$. Cut lines (a) and (b) correspond to $q={\rm
const}$ and plasmon dispersion relation at carrier density $n=10^{10}$ cm-2.
Plasmonic drag responsivity evaluated at plasmon dispersion relation at
various carrier densities

Experimental measurement of non-linear conductivity at finite wave vector $q$
is very challenging as the amplitude of electric field $E_{0}$ in the 2d plane
is different from the incident field $E_{\rm inc}$. This difference stems from
strong self-consistent field effects which, in fact, are responsible for
launching of plasmons. A more common measurable quantity is the photocurrent
responsivity $R_{J}=j^{(2)}/S$, where $S$ is the incident power density.
Below, we provide an upper bound for plasmonic drag responsivity from energy
balance considerations.

Under perfect matching conditions, the power flow of incident electromagnetic
wave $S$ is fully transformed into the power flow carried by 2d plasmon
$S_{\rm pl}$. A nearly perfect conversion is attainable under proper design of
grating couplers. In simpler systems, such as metal edges, the conversion
coefficient is well below unity, $S_{\rm pl}/S\ll 1$. Approximating $S_{\rm
pl}\approx S$, we obtain a natural upper bound of plasmonic drag responsivity.

The power flow density ${\bf s}$ carried by plasmon in a spatially dispersive
2d system is the sum of electromagnetic and kinetic contributions Landau  _et
al._ (2013):

${\bf s}=\frac{c}{8\pi}\left[{\bf
E}\times\textbf{H}^{*}\right]+\frac{1}{4}\frac{\partial{\rm
Im}\sigma^{(1)}(q,\omega)}{\partial\textbf{q}}\left(\textbf{E},\textbf{E}^{*}\right)\delta(z),$
(5)

where $\sigma^{(1)}(q,\omega)$ is the linear conductivity of graphene given in
the classical limit by

$\sigma^{(1)}(q,\omega)=ig\frac{e^{2}}{\hbar}\frac{E_{F}}{2\pi\hbar}\frac{\omega}{q^{2}v_{0}^{2}}\Bigg{[}\frac{\omega}{\sqrt{\omega^{2}-q^{2}v_{0}^{2}}}-1\Bigg{]}$
(6)

Integrating the flow density (5) over the vertical coordinate $z$, we obtain
the full power flow $S_{\rm pl}=\int{s_{x}dz}$ in the form

$S_{\rm pl}=\left.\frac{E_{0}^{2}}{4}\frac{1}{q}\frac{\partial}{\partial
q}\Big{[}q\,\sigma^{(1)}(q,\omega)\Big{]}\right|_{{\omega=\omega_{\rm
pl}(q)}}.$ (7)

It is now apparent that the the power flow carried by a plasmon is diverging
as its velocity approaches Fermi velocity. Physically, it comes from very
large contribution of carrier kinetic energy to the net energy flow. Formally,
it comes from differentiating singular conductivity ${\rm
Im}\sigma^{(1)}\propto[\omega^{2}-q^{2}v_{0}^{2}]^{-1/2}$; differentiation
enhances the strength of singularity.

It is now becoming clear that the current responsivity $R_{J}=j^{(2)}/S_{\rm
pl}$ is bounded from above as the plasmon phase velocity approaches the Fermi
velocity. Indeed, a singular growth in $\sigma^{(2)}$ at $\omega/q\rightarrow
v_{0}$ is overwhelmed by a faster growth in $S_{\rm pl}$. This situation is
illustrated in Fig. 4. Panel (a) shows $R_{J}$ as a function of independent
$q$ and $\omega$, naturally, the maximum of $R_{J}$ lies above the singular
line $\omega=qv_{0}$. Panel (b) shows the responsivity evaluated at $\omega$
and $q$ bound by dispersion (4); again, this function has a pronounced
maximum.

Instructively, it is possible to present the current responsivity in a
universal form at $T=0$:

$R_{J}=\frac{e}{E_{F}}f\left(\widetilde{v}_{ph}\right),$ (8)

where $f\left(\widetilde{v}_{ph}\right)$ is the dimesnionless function of a
single dimensionless parameter, the scaled phase velocity
$\widetilde{v}_{ph}=\omega/qv_{0}$. The function $f(\widetilde{v}_{ph})$
reaches a maximum value of $0.243$ at $\omega=1.38qv_{0}$. Remarkably, the
result does not depend on carrier density and dielectric environment.

Figure 5: Plasmonic drag responsivity vs Fermi energy evaluated at different
frequencies of incoming radiation. The wave vector is obtained from the
plasmon dispersion relation. Temperature $T=300$ K.

The ultimate plasmonic drag responsivity (8) can be compared to that of a
perfect photovoltaic cell $R_{\rm pv}<e/\hbar\omega$. Naturally, the drag
responsivity is below the photovoltaic limit, as all our calculations were
performed in the classical domain $E_{F}\gg\hbar\omega$. Nevertheless, the
maximum responsivities of $2...20$ A/W in Fig. 4 are large compared to those
of typical graphene photodetectors operating in the terahertz Bandurin  _et
al._ (2018b); Vicarelli  _et al._ (2012) and infrared Badioli  _et al._ (2014)
frequency ranges. It should be also noted that plasmonic drag mechanism
provides fast photoresponse Muravev  _et al._ (2016), as the timescale for
decay of photocurrent is the momentum relaxation time in the bulk.

All previous calculations carried out at $T=0$ have indicated that plasmonic
drag responsivity benefits from low carrier density. At finite temperature,
the density is limited by thermal excitation of carriers. The plasmon velocity
is sensitive to carrier density, and also cannot become very close to Fermi
velocity $v_{0}$. Account of finite temperature amounts to a simple
replacement of ’effective Fermi energy’ in expressions for linear conductivity
and plasmon dispersion

$E_{F}\rightarrow
kT\ln\left(1+e^{E_{F}/kT}\right)\left(1+e^{-E_{F}/kT}\right).$ (9)

It is possible to evaluate $R_{J}$ at finite temperature with full account of
thermally excited carriers, which is done in Fig. 5 at $T=300$ K. The obtained
dependence of plasmonic drag responsivity on Fermi energy is an anti-symmetric
function of $E_{F}$, with a maximum located at $E_{F}\sim kT$. The functional
dependence of $R_{J}$ on frequency is quite peculiar. Namely, the responsivity
grows with increasing frequency at large Fermi energies. This growth is
associated with increased plasmon wave vector at higher frequencies, and
higher average momentum transferred to an electronic system.

## References

  * Song  _et al._ (2011) J. C. W. Song, M. S. Rudner, C. M. Marcus, and L. S. Levitov, Nano Letters 11, 4688 (2011), arXiv:1105.1142 .
  * Muravev and Kukushkin (2012) V. M. Muravev and I. V. Kukushkin, Applied Physics Letters 100, 10 (2012).
  * Tielrooij  _et al._ (2015) K. J. Tielrooij, L. Piatkowski, M. Massicotte, A. Woessner, Q. Ma, Y. Lee, K. S. Myhro, C. N. Lau, P. Jarillo-Herrero, N. F. Van Hulst, and F. H. Koppens, Nature Nanotechnology 10, 437 (2015), arXiv:1504.06487 .
  * Cai  _et al._ (2014) X. Cai, A. B. Sushkov, R. J. Suess, M. M. Jadidi, G. S. Jenkins, L. O. Nyakiti, R. L. Myers-Ward, S. Li, J. Yan, D. K. Gaskill, T. E. Murphy, H. D. Drew, and M. S. Fuhrer, Nature Nanotechnology 9, 814 (2014).
  * Echtermeyer  _et al._ (2014) T. J. Echtermeyer, P. S. Nene, M. Trushin, R. V. Gorbachev, A. L. Eiden, S. Milana, Z. Sun, J. Schliemann, E. Lidorikis, K. S. Novoselov, and A. C. Ferrari, Nano Letters 14, 3733 (2014), arXiv:1402.1266 .
  * Gibson  _et al._ (1970) A. F. Gibson, M. F. Kimmitt, and A. C. Walker, Applied Physics Letters 17, 75 (1970), https://doi.org/10.1063/1.1653315 .
  * Karch  _et al._ (2010) J. Karch, P. Olbrich, M. Schmalzbauer, C. Zoth, C. Brinsteiner, M. Fehrenbacher, U. Wurstbauer, M. M. Glazov, S. A. Tarasenko, E. L. Ivchenko, D. Weiss, J. Eroms, R. Yakimova, S. Lara-Avila, S. Kubatkin, and S. D. Ganichev, Phys. Rev. Lett. 105, 227402 (2010).
  * Vengurlekar and Ishihara (2005) A. S. Vengurlekar and T. Ishihara, Applied Physics Letters 87, 091118 (2005), https://doi.org/10.1063/1.2037851 .
  * Durach and Noginova (2016) M. Durach and N. Noginova, Phys. Rev. B 93, 161406 (2016).
  * Popov (2013) V. V. Popov, Applied Physics Letters 102, 253504 (2013), https://doi.org/10.1063/1.4811706 .
  * Olbrich  _et al._ (2016) P. Olbrich, J. Kamann, M. König, J. Munzert, L. Tutsch, J. Eroms, D. Weiss, M. H. Liu, L. E. Golub, E. L. Ivchenko, V. V. Popov, D. V. Fateev, K. V. Mashinsky, F. Fromm, T. Seyller, and S. D. Ganichev, Physical Review B 93, 1 (2016), arXiv:1510.07946 .
  * Popov  _et al._ (2015) V. V. Popov, D. V. Fateev, E. L. Ivchenko, and S. D. Ganichev, Physical Review B 91, 235436 (2015), arXiv:1505.06847 .
  * Nikulin  _et al._ (2021) E. Nikulin, D. Mylnikov, D. Bandurin, and D. Svintsov, Phys. Rev. B 103, 085306 (2021).
  * Note (1) The plasmonic drag effect is conceptually equivalent to the distributed resistive self-mixing Sakowicz  _et al._ (2011) and Dyakonov-Shur Sakowicz  _et al._ (2011) rectification. This can be seen by presenting the photocurrent via linear-response electric field and its gradients, all three mechanisms will be described by identical expressions.
  * Chaplik (1972) A. Chaplik, Sov. Phys. JETP 35, 395 (1972).
  * Nehls  _et al._ (1996) J. Nehls, T. Schmidt, U. Merkt, D. Heitmann, A. G. Norman, and R. A. Stradling, Phys. Rev. B 54, 7651 (1996).
  * Bandurin  _et al._ (2021) D. Bandurin, E. Mönch, K. Kapralov, I. Phinney, K. Lindner, S. Liu, J. Edgar, I. Dmitriev, P. Jarillo-Herrero, D. Svintsov, _et al._ , arXiv preprint arXiv:2106.02117 (2021).
  * Lundeberg  _et al._ (2017) M. B. Lundeberg, Y. Gao, R. Asgari, C. Tan, B. V. Duppen, M. Autore, P. Alonso-gonzález, A. Woessner, K. Watanabe, T. Taniguchi, R. Hillenbrand, J. Hone, M. Polini, and F. H. L. Koppens, Science 357, eaan2735 (2017), arXiv:1704.05518 .
  * Svintsov (2018) D. Svintsov, Physical Review B 97, 121405 (2018), arXiv:arXiv:1710.05054v1 .
  * Iranzo  _et al._ (2018) D. A. Iranzo, S. Nanot, E. J. Dias, I. Epstein, C. Peng, D. K. Efetov, M. B. Lundeberg, R. Parret, J. Osmond, J. Y. Hong, J. Kong, D. R. Englund, N. M. Peres, and F. H. Koppens, Science 360, 291 (2018), arXiv:1804.01061 .
  * Epstein  _et al._ (2020) I. Epstein, D. Alcaraz, Z. Huang, V.-V. Pusapati, J.-P. Hugonin, A. Kumar, X. M. Deputy, T. Khodkov, T. G. Rappoport, J.-Y. Hong, N. M. R. Peres, J. Kong, D. R. Smith, and F. H. L. Koppens, Science 368, 1219 (2020), arXiv:2002.00366 .
  * Bylinkin  _et al._ (2019) A. Bylinkin, E. Titova, V. Mikheev, E. Zhukova, S. Zhukov, M. Belyanchikov, M. Kashchenko, A. Miakonkikh, and D. Svintsov, Physical Review Applied 11, 054017 (2019).
  * Bandurin  _et al._ (2018a) D. A. Bandurin, D. Svintsov, I. Gayduchenko, S. G. Xu, A. Principi, M. Moskotin, I. Tretyakov, D. Yagodkin, S. Zhukov, T. Taniguchi, K. Watanabe, I. V. Grigorieva, M. Polini, G. N. Goltsman, A. K. Geim, and G. Fedorov, Nature Communications 9, 5392 (2018a), arXiv:1807.04703 .
  * Fateev  _et al._ (2017) D. V. Fateev, K. V. Mashinsky, and V. V. Popov, Applied Physics Letters 110 (2017), 10.1063/1.4975829.
  * McLeod  _et al._ (2014) A. S. McLeod, P. Kelly, M. D. Goldflam, Z. Gainsforth, A. J. Westphal, G. Dominguez, M. H. Thiemens, M. M. Fogler, and D. N. Basov, Phys. Rev. B 90, 085136 (2014).
  * Ryzhii (2006) V. Ryzhii, Japanese Journal of Applied Physics, Part 2: Letters 45, 923 (2006).
  * Landau  _et al._ (2013) L. D. Landau, L. Pitaevskii, and E. Lifshitz,  _Electrodynamics of continuous media_ , Vol. 8 (elsevier, 2013).
  * Bandurin  _et al._ (2018b) D. A. Bandurin, I. Gayduchenko, Y. Cao, M. Moskotin, A. Principi, I. V. Grigorieva, G. Goltsman, G. Fedorov, and D. Svintsov, Applied Physics Letters 112, 141101 (2018b).
  * Vicarelli  _et al._ (2012) L. Vicarelli, M. S. Vitiello, D. Coquillat, A. Lombardo, A. C. Ferrari, W. Knap, M. Polini, V. Pellegrini, and A. Tredicucci, Nature Materials 11, 865 (2012), arXiv:1203.3232 .
  * Badioli  _et al._ (2014) M. Badioli, A. Woessner, K. J. Tielrooij, S. Nanot, G. Navickaite, T. Stauber, F. J. García De Abajo, and F. H. Koppens, Nano Letters 14, 6374 (2014).
  * Muravev  _et al._ (2016) V. M. Muravev, V. V. Solov’ev, A. A. Fortunatov, G. E. Tsydynzhapov, and I. V. Kukushkin, JETP Letters 103, 792 (2016).
  * Sakowicz  _et al._ (2011) M. Sakowicz, M. B. Lifshits, O. A. Klimenko, F. Schuster, D. Coquillat, F. Teppe, and W. Knap, Journal of Applied Physics 110, 054512 (2011).

##  Appendix A Сalculation of nonlinear conductivity

We calculate $\sigma^{(2)}={\langle j_{2}\rangle}_{t,x}/E_{0}^{2}$ through the
integral of the current.

${\langle
j_{2}\rangle}_{t,x}=-\frac{g}{(2\pi\hbar)^{2}}e\iint\limits_{-\infty}^{+\infty}v_{x}{\langle
f_{2}\rangle}_{t,x}dp_{x}dp_{y}$ (10)

After Revealing the derivative in ${\langle f_{2}\rangle}_{t,x}$ we use polar
system of coordinates, where:

$\begin{gathered}dp_{x}dp_{y}=pdpd\phi\\\ \frac{\partial f_{0}}{\partial
p_{x}}=\frac{\partial f_{0}}{\partial p}\cos{\phi}\\\
\frac{\partial^{2}f_{0}}{\partial
p_{x}^{2}}=\cos^{2}{\phi}\frac{\partial^{2}f_{0}}{\partial
p^{2}}+\frac{\sin^{2}{\phi}}{p}\frac{\partial f_{0}}{\partial p}\\\
\frac{\partial v_{x}}{\partial
p_{x}}=\frac{\sin^{2}{\phi}}{p}v_{0}\end{gathered}$ (11)

Momentum integrals yield to:

$\begin{gathered}\int_{0}^{+\infty}\frac{\partial f_{0}}{\partial
p}\,dp=-f_{0}(0)\\\ \int_{0}^{+\infty}\frac{\partial^{2}f_{0}}{\partial
p^{2}}p\,dp=f_{0}(0)\end{gathered}$ (12)

so, the $\sigma^{(2)}$ yields:

$\boxed{\sigma^{(2)}=-\frac{g}{(2\pi\hbar)^{2}}\frac{e^{3}f_{0}(0)}{q^{2}v_{0}}I_{\phi}(a,b)}$
(13)

Where:

$I_{\phi}(a,b)=\int\limits_{0}^{2\pi}\Bigg{\\{}\frac{\cos{\phi}(\cos^{2}{\phi}-1/2)}{(\cos{\phi}-a)^{2}+b^{2}}+\frac{\cos^{2}{\phi}\sin^{2}{\phi}(\cos{\phi}-a)}{\big{[}(\cos{\phi}-a)^{2}+b^{2}\big{]}^{2}}\Bigg{\\}}d\phi$
(14)

with $a=\omega/(qv_{0})$ and $b=1/(qv_{0}\tau)$ and that can be analitically
calculated from the theory of resudues. We need replace a variable to complex
units:

$\begin{gathered}z=e^{i\phi}\\\ \cos{\phi}=\frac{1+z^{2}}{2z}\\\
d\phi=\frac{1}{iz}dz\\\ \end{gathered}$ (15)

Inside of circle with radius 1 we have 3 poles:

$\begin{gathered}z=a-ib-\sqrt{(a-ib)^{2}-1}\;\;(second\>order)\\\
z=a+ib-\sqrt{(a+ib)^{2}-1}\;\;(second\>order)\\\
z=0\;\;(first\>order)\end{gathered}$ (16)

Thus, we can find residues of sub-integral function and find exact form for
the integral (14).

$\boxed{I_{\phi}(a,b)=2\pi
i\Big{(}ia+\frac{(a-ib)\sqrt{(a-ib)^{2}-1}}{4b}-\frac{(a+ib)\sqrt{(a+ib)^{2}-1}}{4b}\Big{)}}$
(17)

It is remarkable that second and third terms are complex conjugates, so the
sum of them is imaginary. After multiplying the entire bracket by $2\pi i$, it
become real. The integral is real for all $a$ and $b$.

The expression (17) can be interpreted as real function by allocation of the
main or secondary branch of the square root of complex function inside. It
leads to:

$\begin{gathered}I_{\phi}(a,b)=2\pi\Big{[}-a+\frac{\sqrt{|z|}}{2}\Big{(}\frac{a}{b}\sin{\frac{\theta}{2}}+\cos{\frac{\theta}{2}}\Big{)}\Big{]}\\\
\theta=\arccos{\Big{(}\frac{a^{2}-b^{2}-1}{|z|}\Big{)}}\\\
z=(a+ib)^{2}-1\end{gathered}$ (18)

It is worth noting that the value range of subcortical function $z$ (15) lies
in the upper part of the complex plane. For this reason we use the function
"$\arccos$" and the expression (18) turns out to be correctly defined for all
values of $a$ and $b$.

The system (18) is resource-intensive for calculating in computer programs, so
we are are recommend to use the equation (17), besides it contains a multi-
valued complex function, that does not affect the computation.

In collisionless electron plasma ($b\rightarrow 0$ and therefore
$\tau\rightarrow\infty$) the expression (14). leads to:

$\boxed{I_{\phi}^{(approximate)}(a)=2\pi\Bigg{(}\frac{2a^{2}-1}{2\sqrt{a^{2}-1}}-a\Bigg{)}}$
(19)

So, the current (10) leads to the expession of $\sigma^{(2)}$ with angle
integral that is easy calculated from (14) for all $a$ and $b$:

The nonlinear conductivity (13) in the approximation of fast waves
($q\rightarrow 0$ and low damping ($\omega\tau\rightarrow\infty$) leads to
next expresson:

$\sigma^{(2)}=-\frac{g}{(2\pi\hbar)^{2}}\frac{e^{3}f_{0}(0)}{v_{0}^{2}}2\pi\frac{\omega}{q^{3}}$
(20)