Peter Grünberg Institute 6 & 9, Forschungszentrum Jülich GmbH, 52425 Jülich,
Germany RWTH Aachen University] 2nd Institute of Physics B and JARA-FIT, RWTH
Aachen University, Otto-Blumenthal-Straße, 52074 Aachen, Germany AMO Gmbh]
Advanced Microelectronic Center Aachen (AMICA), AMO GmbH, Otto-Blumenthal-Str.
25, 52074 Aachen, Germany University of Wuppertal] University of Wuppertal,
42285 Wuppertal, Germany [AMO Gmbh] Advanced Microelectronic Center Aachen
(AMICA), AMO GmbH, Otto-Blumenthal-Str. 25, 52074 Aachen, Germany RWTH Aachen
University] 2nd Institute of Physics B and JARA-FIT, RWTH Aachen University,
Otto-Blumenthal-Straße, 52074 Aachen, Germany

# Local spots of viscous electron flow in graphene at room temperature and
moderate mobility

Sayanti Samaddar  2nd Institute of Physics B and JARA-FIT, RWTH Aachen
University, Otto-Blumenthal-Straße, 52074 Aachen, Germany  Jeff Strasdas  2nd
Institute of Physics B and JARA-FIT, RWTH Aachen University, Otto-Blumenthal-
Straße, 52074 Aachen, Germany  Kevin Janßen  2nd Institute of Physics B and
JARA-FIT, RWTH Aachen University, Otto-Blumenthal-Straße, 52074 Aachen,
Germany  Tjorven Johnsen  [  Zhenxing Wang  [  Daniel Neumaier  [  Marcus
Liebmann  [  Markus Morgenstern<EMAIL_ADDRESS>2nd Institute
of Physics B and JARA-FIT, RWTH Aachen University, Otto-Blumenthal-Straße,
52074 Aachen, Germany

###### Abstract

Dominating electron-electron scattering enables viscous electron flow
exhibiting hydrodynamic current density patterns such as Poiseuille profiles
or vortices. The viscous regime has recently been observed in graphene by non-
local transport experiments and mapping of Poiseuille profiles up to room
temperature. Here, we probe the current-induced surface potential maps of
graphene field effect transistors using scanning probe microscopy at room
temperature. We discover the appearance of $\mu$m large areas close to charge
neutrality, where the current induced electric field opposes the externally
applied field. By estimating the scattering lengths from the gate dependence
of local electric fields, we find that these areas exhibit a dominating
electron-electron scattering as expected for viscous flow. We map the
respective meandering electric fields and carefully rule out artifacts such as
by source-drain voltage induced local doping. Our results imply that viscous
electron flow is omnipresent in graphene devices, even at moderate mobility.

###### keywords:

graphene, electron viscosity, negative electric fields, field effect,
electrostatic force microscopy, Kelvin probe force microscopy

##  1 Introduction

Since electron-electron scattering is momentum conserving, if Umklapp
scattering is suppressed, the corresponding electric resistance is not related
to momentum relaxation, but to viscous properties of the electron liquid.1, 2,
3 Indeed, Navier-Stokes type equations have been employed at dominating
electron-electron scattering to calculate resistance and charge flow patterns
4, 5, 6, 1. Most notably, it has been shown that the resistance of a
constriction can drop below its ballistic Landauer-Büttiker-type value via
lateral drag as dubbed the Gurzhi effect,7, 8 that a Poisseuile flow implying
an inverted parabolic velocity profile appears across a ribbon 6, 9, 10, 11
and that vortices of current flow can develop 12 sideways from a current
injection point 6, 13, 14, 15, 16, 17 or within a disorder potential 18, 19.
Vortices at current injection points develop also without a dominating viscous
flow and are, hence, not a smoking gun.20, 15

Experimentally, first evidence of the Gurzhi effect has been found for GaAs
constrictions.21 More recently, indications of a dominating viscous electron
flow appeared in other 2D materials such as graphene 22, 23, 24, 25, 26, 27,
28, 29, 30 or PdCoO2 31 as well as in the 3D Dirac- and Weyl-type materials
PtSn432 and WP233, respectively. These encouraging results using electric or
heat transport experiments, partly at optical frequencies 25, 34, revealed the
viscous flow indirectly via its detailed parameter dependence. A real-space
visualization has been accomplished for graphene displaying the Poiseuille
charge flow profile 35 up to room temperature 36 and its transition to ohmic
37 or ballistic transport profiles 35. Moreover, a GaAs constriction in the
viscous regime has been probed by scanning gate microscopy without conclusive
interpretation of the observed patterns 38.

Here, we employ Kelvin probe force microscopy (KPFM)39 and electrostatic force
microscopy (EFM)40 under current flow as previously applied to static graphene
at ambient conditions.41, 42, 43, 44 To tune the charge carrier density, we
use a backgate voltage $V_{\rm gate}$ that is sufficiently screened to avoid a
direct influence onto the probing cantilever. This enables us to map the
current induced electric potentials at different $V_{\rm gate}$ and source-
drain voltages $V_{\rm SD}$. Surprisingly, we find areas, where the local
potential drop opposes the externally applied $V_{\rm SD}$. Such areas appear
increasingly, if the sample is tuned towards charge neutrality, i.e., the
Dirac point voltage $V_{\rm D}$. One explanation is a novel artifact in KPFM
and EFM, that we discovered during this study and dub source-drain voltage
induced local doping (SDILD). Its effect including opposing fields can be
quantitatively deduced from the measured map of electron concentration and,
hence, identified. For the remaining areas of inverted potential drops, we use
the gate dependence of the measured potentials to estimate a local electron-
disorder scattering length $l_{\rm dis}$ that turns out to be larger than the
local electron-electron scattering length $l_{\rm ee}$ for charge carrier
densities up to about $\pm 5\cdot 10^{16}$/m2 around $V_{\rm D}$. Thus,
consistently, areas of relatively short $l_{\rm ee}$ exhibit inverted lateral
electric fields as expected in the hydrodynamic regime.22, 45, 17, 30, 1 Since
the sample exhibits a moderate mobility $\mu=2500-4000$ cm2/Vs and the effects
are observed at $300$ K, our results imply that viscous electron flow is
ubiquitous in graphene devices.

##  2 Results and discussion

We use a graphene monolayer (Graphenea SE) deposited on a Si(100)/SiN(150 nm)
backgate with source and drain electrodes (Ni/Al (12/50 nm)) structured via
optical lithography (Supplementary Section 1).46 The large drain contact (Fig.
1a), set to ground during all measurements, protects the cantilever from
direct influences of $V_{\rm gate}$. A commercial atomic force microscope
(Bruker Dimensions Icon PT) enables EFM and KPFM with lateral resolution down
to 20 nm, while varying, both, $V_{\rm gate}$ and $V_{\rm SD}$. Both methods
are applied in lift mode. Firstly, one records a tapping mode image of the
surface with cantilever amplitude of 30–33 nm. Afterwards, the average tip-
sample distance is increased by $z_{\rm tip}\approx 30$ nm to determine maps
of the contact potential difference $V_{\rm CPD}(x,y)=W_{\rm tip}-W_{\rm
sample}(x,y)$. Here, $W_{\rm tip}$ and $W_{\rm sample}(x,y)$ are the work
function of the tip and of the sample area below the tip at position $(x,y)$,
respectively. Thus, both methods map the surface potential of graphene. At
$V_{\rm SD}=0$ V, it is directly related to the charge carrier density
$n_{0}(x,y)$, i.e., the density of excess electrons with respect to charge
neutrality. Deducing $n_{0}(x,y)$ from $V_{\rm CPD}(x,y)$ requires the contact
potential difference $V^{0}_{\rm CPD}$ at charge neutrality 47, 48. This is
determined straightforwardly by tuning $V_{\rm gate}$ to charge neutrality,
i.e., to the maximum of two-point resistance as function of $V_{\rm gate}$
(Fig. 1e) (Supplementary Section 5).

Figure 1: Measurement scheme, sample characterization and inverted electric
fields. (a) Setup to image lateral electric fields induced by charge transport
$E_{x}^{\rm meas}(x,y)$, $V_{\rm gate}$: gate voltage, $V_{\rm SD}$: source-
drain voltage, $V_{\rm tip}$: tip voltage, $f_{\rm drive}$: drive frequency of
cantilever deflection voltage $V_{\rm piezo}(t)$ for EFM and $V_{\rm tip}(t)$
for KPFM, serving as reference for lock-in detection. The large, grounded
drain pad (beige) protects the cantilever (violet) from influences of $V_{\rm
gate}$. Cantilever size and sample geometry are roughly to scale. (b) Red:
Measured phase of cantilever oscillation on graphene. Black: parabolic fit
with marked, deduced contact potential difference $V_{\rm CPD}=W_{\rm
tip}-W_{\rm sample}$ ($V_{\rm tip}$ at maximum). Blue dots: three phases
probed subsequently at the same location as the red data points. Such three
phases are recorded for each position $(x,y)$ to reconstruct the full parabola
and, hence, $V_{\rm CPD}(x,y)$ (Supplementary Section 3). (c) Blue, red:
$V_{\rm CPD}(x)$ along the same line, but at different $V_{\rm SD}$, $V_{\rm
gate}=85$ V (charge neutrality). Grey areas: topography along the same line
with graphene, source and drain electrode marked. (d) Blue: $\frac{V_{\rm
transport}}{V_{\rm SD}}(x):=\frac{[V_{\rm CPD}(x,V_{\rm SD})-V_{\rm
CPD}(x,V_{\rm SD}=0\hskip 1.42262pt{\rm V})]}{V_{\rm SD}}$ using the two
$V_{\rm CPD}(x,V_{\rm SD})$ from c. Red: deduced $E_{x}^{\rm meas}(x)$ via the
formula between c and d. Pink background highlights inverted $E_{x}^{\rm
meas}(x)$ areas. (e) Red: Spatially averaged $\langle E_{x}^{\rm
meas}(x,y)\rangle\cdot L/I_{\rm SD}$ representing the sample averaged
resistance as probed by KPFM, $L=18\,{\rm\mu}$m: graphene sample length along
$V_{\rm SD}$, $I_{\rm SD}$: source-drain current. Blue: simultaneously
measured two point resistance, $V_{\rm D}$: Dirac point. (f) Topography of
graphene area (tapping mode AFM). (g) $E_{x}^{\rm meas}(x,y)$ in the area of
the white rectangle in f probed by KPFM, $V_{\rm gate}=80$ V, $V_{\rm SD}=500$
mV, $z_{\rm tip}=\,$ 30 nm, amplitude of $V_{\rm tip}(t)$ for lock-in
detection: $V_{\rm tip,ac}=\,$ 1.1 V, $f_{\rm drive}\,=\,$ 61.34 kHz,
recording time of each $V_{\rm CPD}(x,y)$ image: 26 min.

In EFM, one drives the cantilever at frequency $f_{\rm drive}$ slightly below
the resonance frequency of the free cantilever ($\sim 62$ kHz) and detects the
phase lag $\phi_{\rm EFM}$ of the excitation with respect to the cantilever
oscillation. The resulting $\phi_{\rm EFM}$ as function of tip voltage $V_{\rm
tip}$ features an inverted parabola $\phi_{\rm EFM}(V_{\rm tip})$ (Fig. 1b).
This is due to the cantilever softening that is caused by the attractive
electrostatic interaction between the charged tip and the induced charge
within the sample. The maximum of the parabola is, hence, at vanishing
electrostatic force, i.e., at $V_{\rm tip}=V_{\rm CPD}$.40 To obtain $V_{\rm
CPD}(x,y)$ maps, it is sufficient to probe $\Delta\phi$ at three distinct
$V_{\rm tip}$ such that the fitted parabola (Fig. 1b, black line) reveals
$V_{\rm CPD}$ with an accuracy of about 2 mV (Supplementary Section 3.1). KPFM
is faster in recording, since adapting $V_{\rm tip}=V_{\rm CPD}(x,y)$ via
feedback control that minimizes electrostatic forces, but it is more prone to
voltage signals from the environment such as remaining electric fields from
$V_{\rm gate}$ penetrating to the cantilever.41, 40 This leads to a smooth
background in $V_{\rm CPD}(x,y)$ images that had to be subtracted
(Supplementary Section 2.2).

Figure 1c shows two measured profiles $V_{\rm CPD}(x)$ ($x$: lateral position)
along the same line of graphene at $V_{\rm SD}=0$ V and $V_{\rm SD}=0.25$ V
($V_{\rm gate}=V_{\rm D}=85$ V). The corresponding topography (grey areas)
reveals the positions of the Al electrodes, where $V_{\rm CPD}(x)$ exhibits a
work function difference to graphene. The $V_{\rm CPD}(x)$ lines show
fluctuations across the graphene, attributed to charge puddles 49. Moreover,
application of $V_{\rm SD}$ changes the slope of $V_{\rm CPD}(x)$ indicating
the current induced potential drop. To disentangle charge puddles and current
induced potential $V_{\rm transport}(x)$, we simply subtract the two curves.
Dividing $V_{\rm transport}(x)$ by the applied $V_{\rm SD}$ reveals that about
80 % of $V_{\rm SD}$ drop across graphene (Fig. 1d). Notably, $V_{\rm
transport}/V_{\rm SD}$ partially exhibits a negative slope (highlighted areas
in Fig. 1d), i.e., an inverted voltage drop with respect to $V_{\rm SD}$. The
resulting in-plane electric field $E_{x}^{\rm meas}(x)=-dV_{\rm
transport}(x)/dx$, hence, is also inverted with respect to the electric field
direction caused by the the applied $V_{\rm SD}$ (red line, Fig. 1d). We
consistently display $-E_{x}^{\rm meas}$ in graphs and maps such that the
inverted electric fields always appear negative. The electric field inversion
appears ubiquitously in maps of $E_{x}^{\rm meas}(x,y)$ if recorded close to
charge neutrality, i.e., at $V_{\rm gate}\simeq V_{\rm D}\pm 15\,{\rm V}$
(Fig. 1g). Areas of inverted $E_{x}^{\rm meas}(x,y)$ are partly correlated
with topographic features as at the diagonal fold starting at the lower left
in the topography map (Fig. 1f), where multiple inverted $E_{x}^{\rm
meas}(x,y)$ areas appear on its right (Fig. 1g, see also Fig. 5). We
crosschecked that the simultaneously recorded two-point resistance $R_{\rm
transport}$ agrees with the gate dependence of the spatially averaged
$E_{x}^{\rm meas}(x,y)$ after adequate scaling (Fig. 1e). The data nicely
match except of a small difference of 80-160 $\Omega$ that is nearly the same
as the metal-graphene contact resistance revealed by 4-point measurements of
identically prepared samples 46.

The inverted electric fields with respect to $V_{\rm SD}$ (Fig. 1g) imply a
complex charge redistribution by the current flow. Such charge redistribution
is known to appear in the hydrodynamic regime, e.g., via current vortices 1.
Hence, it is tempting to assume that the inverted $E_{x}^{\rm meas}(x,y)$ is
due to viscous electron flow.1, 2, 3

However, there is an artifact leading to an apparent inverted $E_{x}^{\rm
meas}(x,y)$ at $V_{\rm gate}\simeq V_{\rm D}$. It results from the local
doping change by the applied $V_{\rm SD}$29 (Fig. 2a) and has to be carefully
distinguished from a current induced inverted $E_{x}^{\rm meas}(x,y)$. To
understand the artifact, we firstly recall that $E_{x}^{\rm meas}(x,y)$ at
given $V_{\rm SD}$ is determined from recorded $V_{\rm CPD}(x,y)$ maps via

$E_{x}^{\rm meas}(x,y)=-\frac{d(V_{\rm CPD}(x,y,V_{\rm SD})-V_{\rm
CPD}(x,y,V_{\rm SD}=0))}{dx}$ (1)

.

Assuming that the applied $V_{\rm SD}$ drops rigidly linear across the
graphene, one straightforwardly obtains (Supplementary Section 6, eq. (S16))

$E_{x}^{\mathrm{meas}}(x,y)=-\frac{V_{\rm SD}}{L}-\frac{\hbar
v_{\mathrm{F}}\sqrt{\pi}}{2|e|}\left(\frac{1}{\sqrt{|n(x,y)|}}\frac{dn(x,y)}{dx}-\frac{1}{\sqrt{|n_{0}(x,y)|}}\frac{dn_{0}(x,y)}{dx}\right):=E_{x}^{\mathrm{SDILD}}(x,y).$
(2)

Here, $L=18.5$ $\mu$m is the length of $V_{\rm SD}$ drop, $v_{\rm F}\simeq
10^{6}$ m/s is graphene’s Fermi velocity and $n(x,y)$ ($n_{0}(x,y)$) is the
charge carrier density at applied $V_{\rm SD}$ (without $V_{\rm SD}$). We
include quantum capacitance to calculate $n(x,y)$, $n_{0}(x,y)$ (Supplementary
Section 5). The field is dubbed $E_{x}^{\mathrm{SDILD}}(x,y)$ with SDILD
abbreviating source-drain voltage induced local doping. Obviously,
$E_{x}^{\mathrm{SDILD}}(x,y)$ diverges at $n(x,y)=0$ and $n_{0}(x,y)=0$ with
sign depending on the spatial derivative of the corresponding charge carrier
density, and being opposite for $n(x,y)$ and $n_{0}(x,y)$. After normalizing
$E_{x}^{\mathrm{SDILD}}(x,y)$ to $V_{\rm SD}$ for easier comparison
(analogously $E_{x}^{\mathrm{meas}}(x,y)$):

$\widehat{E}_{x}^{\mathrm{SDILD}}(x,y)=\frac{E_{x}^{\mathrm{SDILD}}(x,y)}{V_{\rm
SD}},\hskip
19.91684pt\widehat{E}_{x}^{\mathrm{meas}}(x,y)=\frac{E_{x}^{\mathrm{meas}}(x,y)}{V_{\rm
SD}},$ (3)

the sign of the divergence also depends on the sign of $V_{\rm SD}$.

Figure 2: Inverted $\widehat{E}_{x}^{\rm meas}(x)$ by source-drain voltage
induced local doping (SDILD). (a) Sample sketch with exemplary charge carrier
density $n_{0}(x)$ (black line, $V_{\rm SD}=0$ V) and resulting $n(x)$ (red
line, $V_{\rm SD}>0$ V). Thin lines: $n_{0}(x)$, $n(x)$ without charge
puddles. (b) Smaller scale of $n_{0}(x)$, $n(x)$ for multiple $V_{\rm SD}$,
$V_{\rm gate}\simeq V_{\rm D}$. The $n_{0}(x)$ profile (black, $V_{\rm SD}=0$
V) corresponds to a Gaussian shaped potential that crosses $E_{\rm F}$, i.e.,
an electron puddle surrounded by holes. Vertical lines: positions of
$n_{0}(x)=0$ (black) and $n(x)=0$ (colored). (c) $\widehat{E}_{x}^{\rm
SDILD}(x)$ deduced from b according to eq. (2), (3) as indicated between b and
c ($\beta=\hbar v_{\mathrm{D}}\sqrt{\pi}/2|e|$). (d) Graphene topography
(tapping mode AFM). (e) Equilibrium doping $n_{0}(x,y)$ of the same area as in
d, deduced from $V_{\rm CPD}(x,y,V_{\rm SD}=0)$ as recorded by KPFM, $V_{\rm
gate}=V_{\rm D}=85\,$V. (f), (h) $\widehat{E}_{x}^{\rm SDILD}(x,y)$ deduced
from $n_{0}(x,y)$ of e via eq. (2), (3) after tilting $V_{\rm CPD}(x,y,V_{\rm
SD}=0)$ to get $n(x,y)$ according to eq. (S14) (Supplementary Section 5.1) at
$V_{\rm SD}\neq 0$ V as marked on the right of g, i. (g), (i) Measured
$\widehat{E}_{x}^{\rm meas}(x,y)$ at the same $V_{\rm SD}$ as f, h. (e), (g),
(i) $z_{\rm tip}=\,$ 30 nm, $V_{\rm ac}=\,$ 1.1 V, $f_{\rm drive}\,=\,$ 61.355
kHz, $V_{\rm CPD}(x,y)$ image sequence ($V_{\rm SD}$ (waiting time/recording
time)): -0.25 V (5’/22’) $\rightarrow$ 0 V (5’/22’) $\rightarrow$ 0.25 V
(5’/22’).

Since $V_{\rm SD}$ necessarily shifts the zero crossings of the charge carrier
density along $x$, parallel to $V_{\rm SD}$, (Fig. 2b) one usually gets two
divergences of opposite sign (Fig. 2c), one at $n_{0}(x)=0$ ($V_{\rm SD}=0$)
and one at the shifted $n(x)=0$ ($V_{\rm SD}\neq 0$). As a result, most of the
$n_{0}(x)=0$ points feature a double divergence of
$\widehat{E}_{x}^{\mathrm{SDILD}}(x)$ in opposite directions (Fig. 2c). Figure
2b and c show an example for an electron puddle $n_{0}(x)$ with diameter 1.4
$\mu$m embedded into a hole density. The puddle implies two zero crossing
$n_{0}(x)=0$ (Fig. 2b). Tilting the $n_{0}(x)$ profile by $V_{\rm SD}$ shifts
the zero crossings laterally (Fig. 2b). The resulting
$\widehat{E}_{x}^{\mathrm{SDILD}}(x)$ for $V_{\rm SD}>0$ V (Fig. 2c,
orange/red line) shows a negative dip for the left position $n_{0}(x)\simeq 0$
(left vertical black line) due to positive slope $dn_{0}(x)/dx$. The shifted
$n(x)=0$ (larger positive slope) exhibits a sharper, but positive bump shifted
to the right. Consistently, the signs of the (virtual) divergences are
inverted by changing either the sign of $V_{\rm SD}$ or the sign of
$dn_{0}(x)/dx$, e.g., at $x\simeq 1.4$ $\upmu$m in Fig. 2b, c. The divergences
are naturally softened by the finite raster size avoiding a precise
$n_{0}(x)=0$.

The SDILD effect is indeed observed in the experiment. Figure 2d shows the
topography and Fig. 2e the $n_{0}(x,y)$ map deduced from $V_{\rm
CPD}(x,y,V_{\rm SD}=0\,{\rm V})$ recorded at $V_{\rm gate}=V_{\rm D}$. An
electron puddle of size 5 $\mu$m is apparent with amplitude $2\times
10^{11}$/cm2. Figure 2f, h display the resulting
$\widehat{E}_{x}^{\mathrm{SDILD}}(x,y)$ for two opposite $V_{\rm SD}$ using
eq. (2), (3) and calculating $n(x,y)$ by tilting the potential and including
quantum capacitance (Supplementary Section 5.1, eq. (S14)). The
$\widehat{E}_{x}^{\mathrm{meas}}(x,y)$ maps recorded by measurement of $V_{\rm
CPD}(x,y,V_{\rm SD})$ and using eq. (1) are displayed in Fig. 2g, i. The
latter include both, current induced electric fields and SDILD effects.
Obviously, most details of the experiment are reproduced by the simulation
revealing that the apparent inverted $\widehat{E}_{x}^{\mathrm{meas}}(x,y)$
(red areas) is induced by SDILD and not by a spatially inhomogeneous current
flow. Thus, one must distinguish SDILD artifacts from real current induced
inverted $\widehat{E}_{x}^{\mathrm{meas}}(x,y)$.

Figure 3 shows an example, where this has been accomplished. Figure 3a
displays the charge carrier density $n_{0}(x,y)$ deduced from $V_{\rm
CPD}(x,y,V_{\rm SD}=0\,{\rm V})$. Figure 3c displays the resulting
$\widehat{E}_{x}^{\mathrm{SDILD}}(x,y)$ at $V_{\rm SD}=0.1$ V according to eq.
(2) and (3), and Fig. 3e the measured $\widehat{E}_{x}^{\mathrm{meas}}(x,y)$
at the same $V_{\rm SD}$. The $\widehat{E}_{x}^{\mathrm{meas}}(x,y)$ map
features an extended dipolar structure (total width: $\sim 4$ $\mu$m)
consisting of two lobes with opposite fields having the inverted electric
fields on the left (Fig. 3e). In contrast, the weak dipolar structure of
$\widehat{E}_{x}^{\mathrm{SDILD}}(x,y)$ is a factor of ten smaller in
amplitude (color bars), a factor of four smaller in $x$ extension and has the
wrong order with inverted field on the right (Fig. 3c). Since the doping
profile $n_{0}(x,y)$ temporarily fluctuates at ambient conditions, we firstly
minimized these fluctuations by adequate waiting times prior to the $V_{\rm
CPD}(x,y)$ recording and by optimizing the sequences to change $V_{\rm SD}$
and $V_{\rm gate}$ (Supplementary Section 4.2, Fig. S8). Moreover, we recorded
$V_{\rm CPD}(x,y,V_{\rm SD}=0\,{\rm V})$ prior and after $V_{\rm
CPD}(x,y,V_{\rm SD}\neq 0\,{\rm V})$ (eq. (1)). The resulting two image sets
consisting of $n_{0}(x,y)$, $\widehat{E}_{x}^{\mathrm{SDILD}}(x,y)$, and
$\widehat{E}_{x}^{\mathrm{meas}}(x,y)$ are compared in Fig. 3 revealing that
the small changes in $n_{0}(x,y)$ (Fig. 3a, b) barely change the dipolar
structure of $\widehat{E}_{x}^{\mathrm{meas}}(x,y)$ (Fig. 3e, f) that in both
cases strongly deviates from $\widehat{E}_{x}^{\mathrm{SDILD}}(x,y)$ (Fig. 3c,
d).

Figure 3: Current induced inverted electric fields not caused by
$\widehat{E}_{x}^{\rm SDILD}$. (a), (b) Graphene doping maps $n_{0}(x,y)$ as
deduced from $V_{\rm CPD}(x,y,V_{\rm SD}=0V)$ recorded by EFM. Both images
display the same sample area, but are recorded prior (a) and after (b)
recording $V_{\rm CPD}(x,y,V_{\rm SD}=0.1\,{\rm V})$ to determine
$\widehat{E}_{x}^{\rm meas}(x,y)$, $V_{\rm gate}=V_{\rm D}=\,$85 V, $z_{\rm
tip}=\,$ 20 nm, $f_{\rm drive}\,=\,$ 61.72 kHz, cantilever amplitude: ${A}=\,$
31 nm, recording time/image: 34 min. (c), (d) $\widehat{E}_{x}^{\rm
SDILD}(x,y)$ deduced from a, b, respectively, using eqs. (2), (3) and $V_{\rm
SD}=0.1$ V. (e), (f) $\widehat{E}_{x}^{\rm meas}(x,y)$ at $V_{\rm SD}=0.1$ V
via eq. (1) employing $V_{\rm CPD}(x,y,V_{\rm SD}=0$ V) from a, b,
respectively, as recorded $\sim 90$ min after (a) and $\sim 30$ min before (b)
the additionally required $V_{\rm CPD}(x,y,V_{\rm SD}=0.1\,{\rm V})$.

Consequently, we attribute the observed dipolar
$\widehat{E}_{x}^{\mathrm{meas}}(x,y)$ feature in Fig. 3e–f, including a large
area of inverted $\widehat{E}_{x}^{\mathrm{meas}}(x,y)$, to a spatially
inhomogeneous current flow. Inverted $\widehat{E}_{x}^{\mathrm{meas}}(x,y)$
areas can be rationalized by a passing current that reduces the charge carrier
density in a nearby area compared to the one that would be there by SDILD
only, i.e., charge puddles are sucked from the passing current without being
compensated by the forward electron flow. This has been observed in
simulations for viscous electron flow within a disorder potential18, 45, 30, 1
where a lateral viscous force pulls electrons out of a relatively quiet area
that is protected from electron flow, e.g., by an upstream obstacle.1, 14, 17

To corroborate the claim of decisive electron viscosity, we demonstrate that
the puddles of inverted $\widehat{E}_{x}^{\mathrm{meas}}(x,y)$ indeed exhibit
conditions favoring viscous current flow, namely $l_{\rm ee}<l_{\rm dis}$.1,
30 Figure 4a$-$f show $\widehat{E}_{x}^{\mathrm{meas}}(x,y)$ of graphene at
various $V_{\rm gate}$. Areas of inverted
$\widehat{E}_{x}^{\mathrm{meas}}(x,y)$ appear close to $V_{\rm gate}=V_{\rm
D}$, here at the electron side ($V_{\rm gate}>V_{\rm D}$). We note that the
measured $n_{0}(x,y)$ in this area does not exhibit any zeroes at $V_{\rm
gate}-V_{\rm D}\geq 3$ V, hence, SDILD is negligible. Firstly, we estimate
$l_{\rm ee}(x,y)$ as local property via $n(x,y)$ deduced from the measured
$n_{0}(x,y)$ (Supplementary Section 5.1, eq. (S14)) using 50, 51

$l_{\rm ee}(x,y)=\frac{4}{\pi}\left(\frac{\hbar v_{\rm F}}{k_{\rm
B}T}\right)^{2}\left<\sqrt{\pi|n(x,y)|}\frac{1}{\ln{\frac{2\hbar
v_{\mathrm{F}}\sqrt{\pi|n(x,y)|}}{k_{\mathrm{B}}T}}}\right>$ (4)

as reasonable estimate for $l_{\rm ee}$ within a factor of 1.5 except at
$|n|<2\cdot 10^{14}$/ m2 (quantum critical regime)52, 53 ($T=298$ K:
temperature, $k_{\rm B}$: Boltzmann constant, $\hbar$: reduced Planck’s
constant).

Secondly, we determine $l_{\rm dis}(x,y)$ from the measured
$E_{x}^{\mathrm{meas}}(x,y)$ far away from charge neutrality, where electron-
electron scattering is irrelevant, using 54

$l_{\rm
dis}(x,y)=\left(\frac{h}{2e^{2}}\right)\left<\frac{1}{\sqrt{\pi|n(x,y)|}\cdot\rho_{\mathrm{local}}(x,y)}\>\right>.$
(5)

Here, $\rho_{\mathrm{local}}(x,y)=E_{x}^{\mathrm{meas}}(x,y)\cdot W/I_{\rm
SD}$ ($e$: electron charge, $W=28$ $\mu$m: width of sample, $I_{\rm SD}$:
source-drain current) represents the local resistivity for areas where the
current flow is unidirectional, as at $|V_{\rm gate}-V_{\rm D}|>70$ V (Fig.
5a). For consistency, we subsequently average $l_{\rm ee}(x,y)$ and $l_{\rm
dis}(x,y)$ across areas larger than the scattering lengths (Supplementary
Section 7).

To extrapolate $l_{\rm dis}(x,y)$ towards $V_{\rm D}$, we have to deduce if
short-range or long-range scattering dominates.54 The inset in Fig. 4g shows
$\left<\rho_{\mathrm{local}}(x,y)\right>(V_{\rm gate}-V_{\rm D})$, averaged
across a larger area of $\sim 100$ $\mu$m2. The maximum is close to $V_{\rm
D}$ as expected. More importantly, a rather constant
$\left<\rho_{\mathrm{local}}(x,y)\right>$ appears at large hole doping
indicating a dominating short-range scattering.54 Hence, we use eq. (5) to
extrapolate $\left<l_{\rm dis}(x,y)\right>$ to $V_{\rm gate}-V_{\rm D}\simeq
0$ V using gate independent $\rho_{\mathrm{local}}(x,y)$.54 This formally
includes the electron-phonon contribution (Supplementary Section 7.2).54, 55.

Figure 4: Relation of inverted ${E}_{x}^{\mathrm{meas}}(x,y)$ areas to the
different scattering lengths. (a)-(f) $\widehat{E}_{x}^{\mathrm{meas}}(x,y)$
of a single area at different $V_{\rm gate}-V_{\rm D}$ (EFM), $V_{\rm D}=81$
V, $z_{\rm tip}=\,$ 30 nm, $f_{\rm drive}\,=\,$ 62.344 kHz, sequence of
$V_{\rm CPD}(x,y)$ images at each $V_{\rm gate}$ ($V_{\rm SD}$ (waiting
time/recording time)): 0 V (5’/20’) $\rightarrow$ 0.4 V (5’/20’). Black
rectangle in d highlights an area of inverted
$\widehat{E}_{x}^{\mathrm{meas}}(x,y)$. (g) Gate dependence of electron-
electron scattering length $l_{\rm ee}$ (blue, black, eq. (4), $T=298$ K) and
electron-disorder scattering length $l_{\rm dis}$ (pink, eq. (5)) averaged
across the black rectangle in d (large symbols) and across a larger area
($12.5$ $\mu{\rm m}$ $\times$ $8.33$ $\mu{\rm m}$, small symbols). The
divergence in $l_{\rm dis}$ at $V_{\rm gate}=V_{\rm D}$ (eq. (5)) and in
$l_{\rm ee}$ (eq. (4)) are avoided by using $n_{0}(x,y)=5\cdot 10^{14}$ m-2 at
$V_{\rm D}$ as naturally appearing by charge puddles (Figs. 2e, 3a–b). These
puddle fluctuations are larger than the quantum critical limit $n_{\rm
QC}=2\times 10^{14}$ m-2 53 or the thermal limit $n_{\rm Th}=4\times 10^{14}$
m-2 51, both principally bordering the low $n$ applicability of eq. (5).51
Inset: $\left<E_{x}^{\mathrm{meas}}(x,y)\right>$ spatially averaged across
$12.5$ $\mu{\rm m}$ $\times$ $8.33$ $\mu{\rm m}$ and scaled to represent the
local resistivity $\left<\rho_{\rm
local}(x,y)\right>=\left<{E}_{x}^{\mathrm{meas}}(x,y)\right>\cdot W/I_{\rm
SD}$ ($W$: sample width). Constant resistivity appears at large hole doping
and, thus, $\left<\rho_{\rm local}(x,y)\right>$ is attributed to short-range
disorder scattering.54

Figure 4g displays $\left<l_{\rm dis}(x,y)\right>$ for the black rectangle
area of Fig. 4d (large pink dots) and for a 100 $\mu$m2 area (small pink dots)
together with the more directly determined $\left<l_{\rm ee}(x,y)\right>$ (eq.
(4), blue, black dots). Obviously, the $V_{\rm gate}$ range with $l_{\rm
ee}<l_{\rm dis}$ is significantly larger for the area of inverted
$\widehat{E}_{x}^{\mathrm{meas}}(x,y)$ than for the extended area. An extended
$V_{\rm gate}$ area with $l_{\rm ee}<l_{\rm dis}$ is consistently observed in
most regions of inverted $\widehat{E}_{x}^{\mathrm{meas}}(x,y)$ (Supplementary
Section 7.3, Fig. S15f). Hence, the conditions for viscous electron flow are
generally realized at $V_{\rm gate}\simeq V_{\rm D}$ in our sample and appear
most pronounced in areas of inverted $\widehat{E}_{x}^{\mathrm{meas}}(x,y)$.
We regard this as strong evidence that hydrodynamic electron flow indeed
causes the experimentally found inverted
$\widehat{E}_{x}^{\mathrm{meas}}(x,y)$ areas (except if attributed to SDILD).
This central result is remarkable, since we operate at moderate mobility,
$\mu=2000-4000$ cm2/Vs, and ambient conditions similar to typical graphene
devices.56 It implies that hydrodynamic electron flow is relevant also in
corresponding graphene applications.

Moreover, our results establish EFM and KPFM as commercially available methods
to probe consequences of hydrodynamic electron flow with high spatial
resolution.45 This is particularly helpful for regimes, where relatively short
length scales of viscous patterns prohibit their detection by negative
vicinity resistance.17, 22 Moreover, EFM is significantly simpler than the
more elaborate scanning probe methods as, e.g., scanning of a carbon nanotube
quantum dot57, scanning NV centers35, 36, 37 or scanning tunneling
potentiometry58. Figure 5 showcases additional experimental observations
highlighting the opportunities for future studies.

Figure 5: Ubiquity of inverted $\widehat{E}_{\rm x}^{\rm meas}(x,y)$.
(a)$-$(f) $\widehat{E}_{x}^{\rm meas}(x,y)$ of another graphene area recorded
by KPFM, $V_{\rm gate}-V_{\rm D}$ marked, $V_{\rm SD}=\,0.5\,V$, $z_{\rm
tip}=\,$ 30 nm, $f_{\rm drive}\,=\,$ 61.34 kHz, $V_{\rm ac}=\,$ 1.1 V,
recording time per $V_{\rm CPD}(x,y)$ image: 21’, waiting time prior to
recording: 5’. (g) same as (d) in different representation: arrows: in-plane
electric field vectors ${\mathbf{E}}^{\rm meas}(x,y)=-\mathbf{\nabla}V_{\rm
transport}(x,y)$, arrow color: blue (red): $x$ component along (opposing)
$V_{\rm SD}$, background color: $n(x,y)$, a background subtraction has been
applied to $V_{\rm CPD}(x,y,V_{\rm SD})$ that slightly influences
${\mathbf{E}}^{\rm meas}(x,y)$ (Supplementary Section 2.2, 4). (h) Topography
of a different graphene area (tapping mode AFM). (i)$-$(j)
$\widehat{E}_{x}^{\rm SDILD}(x,y)$ maps of the area of h deduced from $V_{\rm
CPD}(x,y,V_{\rm SD}=0$ V) recorded by EFM at $V_{\rm gate}=V_{\rm D}$, $V_{\rm
SD}$ marked. (k)$-$(m) $\widehat{E}_{x}^{\rm meas}(x,y)$ maps recorded via EFM
of the same area as in h$-$j, $V_{\rm SD}$ marked, $V_{\rm gate}=V_{\rm D}$,
$z_{\rm tip}=\,$ 30 nm, $f_{\rm drive}\,=\,$ 62.34 kHz, sequence of $V_{\rm
CPD}(x,y)$ images ($V_{\rm SD}$ (waiting time/recording time)): 0 V (2
hours/21’) $\rightarrow$ 0.1 V (5’/21’) $\rightarrow$ 0.2 V (5’/21’)
$\rightarrow$ 0.1 V (5’/21’) $\rightarrow$ 0 V (5’/21’). The last image serves
as reference. (n)$-$(s) $\widehat{E}_{x}^{\rm meas}(x,y)$ recorded by KPFM,
$V_{\rm SD}$ marked, $V_{\rm gate}=V_{\rm D}$, $z_{\rm tip}=\,$ 30 nm, $V_{\rm
ac}=\,$ 1.1 V, $f_{\rm drive}\,=\,$ 62.34 kHz, sequence of $V_{\rm CPD}(x,y)$
images ($V_{\rm SD}$ (waiting time/recording time)): -1 V (2 hours/20’)
$\rightarrow$ -0.5 V (35’/20’) $\rightarrow$ -0.25 V (10’/20’) $\rightarrow$ 0
V (5’/20’) $\rightarrow$ 0.25 V (5’/20’) $\rightarrow$ 0.5 V (10’/20’)
$\rightarrow$ 1 V (35’/20’).

Figure 5a$-$f display $\widehat{E}_{x}^{\rm meas}(x,y)$ of an area recorded by
KPFM at different $V_{\rm gate}$. Again, inverted $\widehat{E}_{x}^{\rm
meas}(x,y)$ areas appear consistently at different $V_{\rm gate}$ around
$V_{\rm gate}\simeq V_{\rm D}$, but, in contrast to Fig. 4a$-$f, on the hole
side, i.e., at negative $V_{\rm gate}-V_{\rm D}$. Again, these areas cannot be
attributed to SDILD by using the probed $n_{0}(x,y)$. Comparison with the
topography of the area (Fig. 1f) reveals a strong influence of a topographic
fold ($\sim 1.4$ nm in height) on the current flow pattern. Figure 5g shows an
arrow plot of the in-plane electric field ${\mathbf{E}}^{\rm
meas}(x,y)=-\mathbf{\nabla}V_{\rm transport}(x,y)$ at charge neutrality
displayed on top of a color plot showcasing $n(x,y)$. This visualizes the
correlation between local charge carrier density and current induced
curvatures of ${\mathbf{E}}$ field patterns (Supplementary Section Fig. S6).
Moreover, Fig. 5h$-$m demonstrate the disentangling of SDILD and current
induced electric fields via comparing $\widehat{E}^{\rm meas}(x,y)$ and
$\widehat{E}^{\rm SDILD}(x,y)$ of the same area. Obviously, most of the
calculated SDILD patterns are reproduced by the measured $\widehat{E}^{\rm
meas}(x,y)$, such that additional, more extended features can be attributed
rather unambiguously to viscous electron flow. Finally, Fig. 5n$-$s
demonstrate how areas of inverted $\widehat{E}_{x}^{\rm meas}(x,y)$ develop
with applied $V_{\rm SD}$. We generally observe that they weaken with
increasing $|V_{\rm SD}|$ without a clear understanding yet.

In conclusion, we have discovered areas of inverted electric field with
respect to the applied source-drain voltage in graphene field-effect
transistors at room temperature and moderate mobility. Via carefully analyzing
artifacts of the measurement such as SDILD, we could attribute several of
these features to local viscous electron flow, in particular, by correlating
its appearance with strongly dominating electron-electron scattering compared
to electron-disorder scattering. This establishes viscous electron flow as
relevant for material parameters used for applications and provides a new
method to study these intriguing electron transport phenomena with high
spatial resolution.

##  3 Supplementary Information

##  S1 Sample Preparation and Transport Characterization

The field effect graphene devices were made from commercially obtained CVD
graphene grown on copper (Graphenia SE) that has been transferred to a SiN
(150 nm)/Si(100) substrate after wet chemical etching of copper with FeCl3,
using PMMA as a supporting layer for the transfer 59. The source and drain
electrodes of each device are connected to the graphene at the edges according
to a procedure described elsewhere 46. The contacts are defined by a polymer
resist AZ5214E via optical lithography followed by removal of graphene from
the exposed areas by oxygen plasma in a reactive ion etching chamber operated
at 100 W for 30 mins. The contact metals are deposited directly afterwards in
order to contact the open bonds of graphene. Firstly, 12 nm Ni is sputtered,
which facilitates edge contact with low contact resistance. Then, 50 nm Al is
deposited via e-beam evaporation followed by lift-off of the resist.
Subsequently, optical lithography is used to remove undesired graphene areas
by etching in oxygen plasma. Contact resistances obtained by this process are
$100-200\,{\rm\Omega}$ slightly depending on the graphene doping 46.

Figure S1a shows an optical image of a typical device with the AFM cantilever
on top. The device design ensures minimum exposure of the cantilever body to
the Si/SiN substrate and, hence, to the electric fields of the gate voltage
$V_{\rm gate}$. Only $2.7\,$ % of the cantilever are in line of sight of the
gate dielectric SiN as barely changing during scanning. Hence, gate voltages
as high as $\pm 100\,V$ can be applied, while requiring less than 3 V of
compensation voltage in KPFM mode as enabled by the Bruker instrument.
Furthermore, to reduce the influences of the source and drain voltages, we
restrict the measurements to graphene areas away from source and drain as
marked, e.g., by the central red rectangle in Fig.S1b.

Figure S1: Experimental Setup. (a) Optical image of the device in top view
with the cantilever above it. The graphene is the darker region below the
cantilever. The dashed rectangle marks the area of b. The $x$ and $y$
coordinates depict the scan directions. (b) Topography of the area marked in a
as recorded by tapping mode AFM. The central, red rectangle marks the region
analyzed in Fig. S7a–f. It is displayed at increased contrast with respect to
the color bar on top that is valid for the rest of the image. (c) Gate
dependent two terminal resistance for the device analyzed in Fig. 4, main
text. The curves are recorded after the time delays as marked, recording time
per curve: 8/40/40/1.5/13/2.8 min (from top (0 hrs) to bottom (3 days)), 0.5
s/measurement point. Mobilities $\mu$ and Dirac point voltages $V_{\rm D}$ are
deduced from fits by eq. (S.1) revealing, additionally, the residual doping
$n^{*}=(1.5\pm 0.3)\times 10^{16}\,{\rm m^{-2}}$ and the contact resistance
$R_{\rm series}=(200\pm 20)\,{\rm\Omega}$. The arrow indicates the direction
of $V_{\rm gate}$ sweep for all curves. (d) Two-point resistance of the same
device as c, but acquired during imaging the maps of Fig. 4a-f, main text .
Two $V_{\rm SD}$ (blue, red) are applied subsequently at each $V_{\rm gate}$,
20 min/measurement point. The arrow marks the direction of $V_{\rm gate}$
sweep. Parameters $\mu$, $V_{\rm D}$, $n^{*}$, and $R_{\rm series}$ result
from a fit by eq. (S.1).

Four devices from a single chip are probed by KPFM/EFM, each for about one
week and, at least, twice. Since the measurements were performed at ambient
conditions, the two terminal resistance $R_{\rm transport}$ of the devices
varied continuously with time concerning the voltage at charge neutrality
$V_{\rm D}$ and the mobility $\mu$. Moreover, a gate voltage hysteresis
appeared, where the sweeps in the reverse direction have been found to be more
reproducible. A series of such sweeps are shown for one device in Fig. S1c.
Consequently, the reverse sweep direction was exclusively used for all
measurements emplyoing sweep speeds not exceeding 0.1 V/s. Moreover, we took
care that $V_{\rm SD}$ was changed in the same direction prior to recording
images that are compared directly. In addition, a waiting time was established
after setting $V_{\rm gate}$ and $V_{\rm SD}$ for each EFM/KPFM image as given
in the respective figure captions. Finally, the two-point resistance was
continuously monitored during imaging revealing minor variations as depicted
in Fig. 1e, main text. Figure S1d shows $R_{\rm transport}(V_{\rm gate})$ as
measured during an image series. It nicely reproduces the fast $R_{\rm
transport}(V_{\rm gate})$ sweeps of the same device recorded without imaging
(Fig. S1c).

The device mobility $\mu$, the gate voltage at charge neutrality $V_{\rm D}$,
the residual doping $n^{*}$ at charge neutrality, and the contact resistance
$R_{\rm series}$ are deduced from fitting the two-terminal $R_{\rm
transport}(V_{\rm gate})$ with a relation developed for long range Coulomb
type disorder 60, 54 (constant mobility):

$R_{\rm transport}=R_{\rm series}+\frac{1}{\mu}\frac{L/W}{\sqrt{\left(C_{\rm
gate}(V_{\rm gate}-V_{\rm D})\right)^{2}+(n^{*}e)^{2}}}$ (S.1)

using the length $L=18.5$ $\mu$m and the width $W=28$ $\mu$m of the graphene
area. The gate capacitance per unit area reads $C_{\rm gate}=\kappa_{\rm
SiN}\epsilon_{0}/t$ with $\kappa_{\rm SiN}=7.5$ and $t=150\,$nm as dielectric
constant and thickness of the dielectric, respectively, and $\epsilon_{0}$ as
the vacuum dielectric constant. The four devices featured $\mu=2500-4000$
${\rm cm^{2}/Vs}$ and $n^{*}=(0.8-2.7)\times 10^{16}\,{\rm m^{-2}}$ with
contact resistance $R_{\rm series}=(160-300)\,{\rm\Omega}$ compatible with
previous measurements of identically produced devices 46.

.

##  S2 Amplitude Modulation KPFM

Figure S2: Kelvin Probe Force Microscopy. (a) Schematic of amplitude
modulation KPFM as implemented in the Dimension Icon AFM setup from Bruker 61.
The lower image shows a profile line $z(x)$ (black) recorded by the cantilever
(violet) in tapping mode AFM with amplitude $A^{\rm TM}_{f_{\rm drive}}$
excited via the voltage $V_{\rm drive}\sin{(2\pi f_{\rm drive}t)}$ applied to
the piezo drive of the cantilever. The middle image shows the profile line
from the lower image (full grey line) and the same line lifted by $z_{\rm
tip}$ (dashed grey line) as traced during recording $V_{\rm CPD}(x)$ (green
line). $V_{\rm CPD}$ is determined at each point by the feedback loop depicted
on top. It minimizes the cantilever amplitude $A_{f_{\rm drive}}$ at frequency
$f_{\rm drive}$, while the AC voltage $V_{\rm AC}$ is applied to the tip with
constant amplitude $V_{\rm exc}$ simultaneously with a DC voltage $V_{\rm
DC}$, that itself is optimized by the feedback. A remaining amplitude
$A_{2f_{\rm drive}}$ appears at frequency $2f_{\rm drive}$ (eq. (S.5)).
$\phi_{f_{\rm drive}}$ measures the relative phase between $V_{\rm AC}$
excitation and detected cantilever oscillation. (b) Topography of graphene
area for one of the measured devices recorded by tapping mode AFM. (c)
Corresponding $V_{\rm CPD}(x,y)$ map, acquired by KPFM at $V_{\rm gate}=V_{\rm
D}=\,$87 V, $V_{\rm SD}=0$ V. (d) Same map as in c after subtracting a second
order polynomial background for horizontal and vertical direction. 62 Image
scan parameters: $z_{\rm tip}=\,$ 30 nm, $A_{f_{\rm drive}}^{\rm TM}=\,$ 30
nm, $f_{\rm drive}\,=\,$ 61.355 kHz, $V_{\rm exc}=1.066$ V, time/image = 11
min.

One of the two scanning probe microscopy methods employed to image the current
induced electric fields is Kelvin Probe Force Microscopy with amplitude
modulation (AM-KPFM). Figure S2a illustrates its implementation in the
Dimension Icon AFM setup from Bruker as used here. It employs a two-path
process conducted for each scan line. In the first path, the topography of the
sample is acquired by tapping mode AFM. In the second path, called the lift
mode, the tip movement mimicks the sample topography at an offset height
$z_{\rm tip}$, while the contact potential difference of the tip with respect
to the sample, $V_{\rm CPD}$, is determined at each point. Hence, the tip
retraces the measured topography from the first path with an additional tip-
sample distance $z_{\rm tip}$. During tapping mode, the cantilever oscillation
is mechanically driven by a piezo-electric element while the cantilever scans
the sample surface employing a $z$-feedback that regulates the tip-sample
distance $z$ to maintain a constant amplitude of the cantilever oscillation
$A_{f_{\rm drive}}^{\rm TM}\approx\,$30 nm. The drive amplitude of the
piezoelectric element is $V_{\rm drive}=2-3$ V at a drive frequency $f_{\rm
drive}$ that is chosen at 50 Hz below the resonance frequency of the free
cantilever $f_{\rm res}$ and given in the respective figure captions.

In the second path, both the $z$-feedback and the piezo drive voltage are
switched off, while a potential feedback regulates the DC voltage applied to
the cantilever $V_{\rm DC}$. An AC modulation $V_{\rm AC}(t)=V_{\rm exc}{\rm
sin}(2\pi f_{\rm drive}t)$ of constant amplitude $V_{\rm exc}$ is superposed
to $V_{\rm DC}$ such that the total tip voltage reads $V_{\rm tip}(t)=V_{\rm
DC}+V_{\rm exc}{\rm sin}(2\pi f_{\rm drive}t)$. This causes the electrostatic
tip-sample interaction force $F_{\rm es}$ to be modulated as 63

$F_{\mathrm{es}}(t)=\frac{1}{2}\,\frac{dC_{\rm
ts}(z)}{dz}(V_{\mathrm{tip}}(t)-V_{\mathrm{CPD}})^{2}=F_{\mathrm{es}}^{0}+(F_{\mathrm{es}})_{f_{\rm
drive}}(t)+(F_{\mathrm{es}})_{2f_{\rm drive}}(t)$ (S.2)

where $C_{\rm ts}(z)$ is the distance dependent tip-sample capacitance and

$\displaystyle F_{\mathrm{es}}^{0}$ $\displaystyle=\frac{1}{2}\frac{dC_{\rm
ts}}{dz}\left[(V_{\mathrm{DC}}-V_{\mathrm{CPD}})^{2}+\frac{V_{\rm
exc}^{2}}{2}\right]$ (S.3) $\displaystyle(F_{\mathrm{es}})_{f_{\rm drive}}(t)$
$\displaystyle=\frac{dC_{\rm ts}}{dz}(V_{\mathrm{DC}}-V_{\mathrm{CPD}})V_{\rm
exc}\,{\rm sin}(2\pi f_{\rm drive}t)$ (S.4)
$\displaystyle(F_{\mathrm{es}})_{2f_{\rm drive}}(t)$
$\displaystyle=\frac{1}{4}\frac{dC_{\rm ts}}{dz}\,V_{\rm exc}^{2}\,{\rm
cos}(4\pi f_{\rm drive}t)$ (S.5)

It can be shown that the oscillating $F_{\rm es}(t)$ causes the cantilever to
vibrate with $z_{f_{\rm drive}}(t)\propto Q\cdot(F_{\rm es})_{f_{\rm
drive}}(t)/k$ at excitation frequency $f_{\rm drive}$, if $f_{\rm drive}$ is
close to $f_{\rm res}$.63 An additional oscillation $z_{2f_{\rm
drive}}(t)\propto(F_{\rm es})_{2f_{\rm drive}}(t)/k$ appears at frequency
$2f_{\rm drive}$. Here, $k\approx\,$3 N/m is the stiffness constant of the
cantilever and $Q\approx 250$ its quality factor (SCM-PIT-V2, Bruker 64). The
amplitude and phase lag of the cantilever oscillation with respect to $V_{\rm
AC}(t)$ is detected by a lock-in amplifier and passed to the potential
feedback that nullifies $z_{f_{\rm drive}}(t)$ by adjusting $V_{\rm DC}$.
Using $V_{\rm CPD}$ as the adjusted $V_{\rm DC}$ as recorded at each position
$(x,y)$, a map $V_{\rm CPD}(x,y)$ results as displayed in Fig. S2c.

Since $z_{f_{\rm drive}}(t)$ is rendered zero by the feedback, the cantilever
oscillates only with $z_{2f_{\rm drive}}(t)$. The amplitude can be estimated
by eq. S.5 using $V_{\rm exc}=1.066$ V, tip-sample distance $d=A_{f_{\rm
drive}}^{\rm TM}+z_{\rm tip}=60$ nm and assuming the simple model of a charged
sphere above an infinite metal surface for the tip-sample capacitance leading
to 65 $\frac{dC_{\rm ts}}{dz}=2\pi\epsilon_{0}\frac{R^{2}}{d(d+R)}$ with tip
radius $R=25$ nm. One obtains an amplitude $A_{2f_{\rm drive}}\approx 0.6$ pm
$\ll d$.

###  S2.1 Cantilever Contribution

The above description is simplified by assuming a homogeneous surface
potential and by neglecting the influence of the electrostatic interaction of
the cantilever body with the device. Both are relevant due to the long range
nature of electrostatic forces. More precisely, the recorded $V_{\rm
CPD}(x,y)$ is described at each point $(x,y)$ by a convolution of the surface
potential map with a point spread function (PSF) deduced from the tip geometry
at distance $z_{\rm tip}$ 66, 40. Since AM-KPFM nullifies the $f_{\rm drive}$
component of the tip sample force $F_{\rm es}$ (eq. (2)) and not the force
gradient $dF_{\rm es}/dz$ as EFM (section S3), the cantilever body as well as
tip areas more distant from the apex influence the measured $V_{\rm CPD}$ 40.
This impact is known to reduce the signal intensity, but barely the spatial
resolution 40. For the cantilever of our measurements (SCM-PIT-V2) and $z_{\rm
tip}=30\,$nm, the approximate PSF, as given for this cantilever type by Xu et
al.,40 reveals a spatial resolution of 60 nm and a reduction of signal
intensity by $20\,\%$ assuming a homogeneous surface potential.

This is compatible with our results (Fig. S3). The relative drop of transport
potential across the graphene as measured by AM-KPFM is only $\sim 80$ % of
the applied $V_{\rm SD}$, but it is 100 % using EFM (Fig. S3b) as expected
(section S3)).40 The relation for KPFM turns out to be independent on the
applied $V_{\rm SD}$ (Fig. S3a) and $V_{\rm gate}$ (Fig. S3b), but changes
with $z_{\rm tip}$ (Fig. S3b) as predicted via its PSF.40 Note that the
potential drop across graphene was determined 3 $\mu$m away from the lateral
metal-graphene interface avoiding influences of the potential inhomogeneities
there. Hence, our results confirm the considerations by PSF40 implying a
residual long range influence.

Figure S3: Transport Potential Drop in KPFM. (a) Transport potential drops
across graphene (colored lines), each averaged from 50 distinct profile lines
at the same $V_{\rm SD}$ and normalised by $V_{\rm SD}$ (right vertical
scale). Various $V_{\rm SD}$ differing by 0.5 V are color coded, $z_{\rm
tip}=30\,$nm, $V_{\rm gate}=85\,$V, $\Delta V:=V_{\rm transport}(x=25\,{\rm\mu
m})-V_{\rm transport}(x=0\,{\rm nm})$. A topography profile probed by tapping
mode AFM is added in grey (left vertical scale). (b) $\Delta V(V_{\rm SD})$
for different $V_{\rm gate}$, $z_{\rm tip}$ as marked (slopes: blue: 0.62,
red: 0.84, green: 0.82, black: 1.0). For the EFM data, $\Delta V$ is deduced
from shifts of the maxima of the EFM parabolas measured on the source
($x=\,41\,{\rm\mu m}$) while the drain is grounded (section S3).

###  S2.2 Background Subtraction

Figure S2c shows a KPFM map of a graphene area recorded at $V_{\rm SD}=$0 V
and $V_{\rm gate}=$87 V, i.e., at charge neutrality. The displayed area is, at
least, $3$ ${\rm\mu}$m away from the electrodes and from exposed SiN areas.
Thus, residual backgrounds cannot result from the tip apex rendering them long
range. Indeed, a dominating, long range contrast developing from the lower
right to the upper left appears. It is likely caused by remaining interactions
of the cantilever body with exposed SiN areas. Such long-range background is
only found in AM-KPFM maps but not in EFM maps corroborating that it is not
caused by the tip apex. The increase of $V_{\rm CPD}(x,y)$ towards the upper
left meets the expectation from the geometry of the exposed SiN areas (Fig.
S1a) and the positive $V_{\rm gate}$. Due to the inclination of $10\degree$ of
the cantilever towards the tip, the distance of the cantilever body part above
SiN decreases relative to the SiN by moving the cantilever into this direction
and, hence, the interaction force to the SiN increases.

To get rid of the long-range background, we apply a 2nd order polynomial
background subtraction using Gwyddion 62. The resulting map after subtraction
is displayed in Fig. S2d featuring an obvious correlation with the
simultaneously recorded topography (Fig. S2b), e.g., at the two nearly
vertically propagating wrinkles or at the two large bumps on the left that
obviously change the local doping. This justifies the subtraction method
featuring a spatial resolution of about 20 nm and a $V_{\rm CPD}$ resolution
of about 10 mV.

Of course, the background subtraction also removes the spatial average of
$V_{\rm CPD}$ and the linear potential drop induced by $V_{\rm SD}$. They are,
however, required for the doping maps and the electric field maps,
respectively. For the electric field maps, we do not apply the background
subtraction, but rely on the subtraction of two $V_{\rm CPD}(x,y)$ maps
recorded at the same $V_{\rm gate}$ (eq. (S.10)) and, hence, removing the
background, that is dominated by the gate voltage that penetrates to the
cantilever, automatically. For the doping maps (eq. (S.13)), we add to the
background subtracted $V_{\rm CPD}(x,y)$ an averaged value $\overline{V}_{\rm
CPD}$ as deduced from a straightforward capacitive charging model of the
graphene by $V_{\rm gate}$ (section S5). This procedure is corroborated by the
excellent agreement between the simulated electric field maps (Fig. 2f,h, main
text) as directly deduced from doping maps (section S5) with the
experimentally measured electric field maps (Fig. 2g,i, main text).

###  S2.3 Noise and Sensitivity

Two major sources of noise are known for KPFM.67 Firstly, thermal fluctuations
of the cantilever oscillation are present due to Brownian motion. Secondly,
sensor noise from detecting the optical beam deflection (OBD) contributes.
Using established formulas,67 we find the thermal noise to be $n_{\rm
thermal}=1.9\,{\rm pm}/\sqrt{{\rm Hz}}$ using our parameters
$Q\,=\,250,\,f_{\rm res}=62.34\,{\rm kHz}\,$ and $T=300\,$K. A typical value
for the OBD sensor noise is $n_{\rm OBD}=100\,{\rm fm}/\sqrt{{\rm Hz}}$,67
i.e., negligible. For the typical bandwidth of the feedback loop $B\approx
5\,$Hz, we get an amplitude noise at $f_{\rm drive}$ of $n_{\rm
thermal/OBD}\sqrt{B}\approx 4$ pm. This imposes a lower limit on the $V_{\rm
CPD}$ precision. At $z_{\rm tip}=30\,$nm and AC voltage amplitude $V_{\rm
AC}=1.066\,$V, we find this limit to be $\delta V_{\rm CPD}=10\,$mV.67 This is
roughly consistent with the noise in $V_{\rm CPD}(x,y)$ maps, e.g., in Fig.
S2d, i.e., we find RMS fluctuations within areas of $0.5\,\mu$m $\times
0.5\,\mu$m of 5–20 mV. It also roughly agrees with the data sheet from Bruker
promising a noise level of $\sim 10\,$mV and with other experiments in the
literature.47, 68, 63

##  S3 Electrostatic Force Microscopy

As second method, we employed Electrostatic Force Microscopy (EFM), well
established to map surface potentials with high resolution 69, 44. Like KPFM
(section S2), it uses a two step process acquiring the topography in a first
path. After lifting the tip for the second path, contrary to KPFM, the
cantilever is oscillated by an AC voltage applied to a piezo-electric element,
while a DC voltage $V_{\rm tip}$ is applied to the tip (Fig. S4a). This
$V_{\rm tip}$ induces an attractive electrostatic force between tip and sample
via image charges in the sample (eq. (S.2)) that causes a decrease in
resonance frequency of the cantilever $f_{0}$ with resulting shift of the
phase vs. frequency curve $\phi_{\rm EFM}(f)$ (Fig. S4e). Probing at the
excitation frequency $f_{\rm drive}<f_{0}$, the electrostatic force changes
the phase lag between cantilever oscillation and voltage oscillation at the
piezo (Fig. S4e) as well as the amplitude of the cantilever oscillation
$A_{\rm EFM}$. The measured $A_{\rm EFM}$ and $\phi_{\rm EFM}$ are both
related to the electrostatic tip-sample force gradient that must be nullified
to detect the local $V_{\rm CPD}$. Note that $\phi_{\rm EFM}$ depicts the
phase lag of the exciting oscillation with respect to the cantilever
oscillation after being nullified prior to each measurement at the resonance
frequency of the free cantilever oscillation.

Figure S4: Electrostatic Force Microscopy. (a) Sketch of the cantilever
(violet) during the second path (lift mode) including the detection scheme
(laser, photodiode). Full grey line: Topography $z(x)$ recorded during the
first path. Dashed grey line: Topography line after lifting the cantilever by
$z_{\rm tip}$ as tracked during the second path. Red line: phase of cantilever
oscillation $\phi_{\rm EFM}(x)$ as recorded during the second path. The
cantilever is continuously excited by $V_{\rm AC,piezo}$, applied to the piezo
element at frequency $f_{\rm drive}$ slightly below the resonance frequency of
the free cantilever $f_{\rm res}$. A DC tip voltage $V_{\rm tip}$ changes the
amplitude $A_{\rm EFM}$ and phase $\phi_{\rm EFM}$ of the cantilever
oscillation at $f_{\rm drive}$ as recorded by lock-in amplifier. (b–d) Three
exemplary $\phi_{\rm EFM}(x,y)$ maps of one graphene area (topography in Fig.
5h, main text) recorded at three different $V_{\rm tip}$ as marked. (e)
Simulation of $\phi_{\rm EFM}(f)$ for two distinct $V_{\rm tip}-V_{\rm CPD}$,
$f_{0}=75\,$kHz, $Q=240$ as present for the used SCM-PIT-V2 cantilever. The
downward arrow marks the decrease in $\phi_{\rm EFM}$ at $f_{\rm drive}$ due
to the electrostatic force induced by $V_{\rm tip}$. (f) Measured $A_{\rm
EFM}$ (blue, left axis) and $\phi_{\rm EFM}$ (red, right axis) as function of
$V_{\rm tip}$ for a centrally located point on graphene, $V_{\rm gate}=65\,$V,
$V_{\rm SD}=0\,$V. $V_{\rm CPD}$ marks the voltage at the $\phi_{\rm EFM}$
maximum. $\Delta\Phi$ marks the phase shift with respect to $V_{\rm CPD}$
induced by $V_{\rm tip}=0.8$ V. The vertical dashed lines with black dots
indicate the three $V_{\rm tip}$ as used for recording $\phi_{\rm EFM}(x,y)$
maps and subsequently deducing $V_{\rm CPD}(x,y)$ maps.

Assuming the cantilever as a damped, driven harmonic oscillator disturbed by a
force that is small compared to the restoring force of the cantilever and that
varies slowly on the scale of the oscillation amplitude, the following
expression applies for the induced resonance frequency shift $\Delta f_{0}$ 70
($k$: cantilever stiffness constant, $F_{\rm es}$: electrostatic force
according to eq. (S.2)).

$\Delta f_{0}=-\frac{f_{0}}{2k}\,\frac{\partial F_{\rm es}}{\partial z}.$
(S.6)

For $f_{\rm drive}\approx f_{0}$ and $Q\cdot\Delta f_{0}/f_{0}\ll 1$, the
corresponding phase shift reads:44

$\Delta\phi_{\rm EFM}(f_{\rm drive})=-\arcsin\left({\frac{2Q}{f_{0}}\,\Delta
f_{0}}\right)\approx-\frac{Q}{k}\,\frac{\partial F_{\rm es}}{\partial z}.$
(S.7)

Hence, the change of $\phi_{\rm EFM}(f_{\rm drive})$ is proportional to the
force gradient. This renders EFM more sensitive to local surface potentials
than AM KPFM, itself minimizing the force instead (section S2.1). Long range
forces acting on the cantilever body contribute much less to EFM via their
less steep gradient 40. Indeed, $\phi_{\rm EFM}(x,y)$ maps are found devoid
from any significant long range background (Fig. S4b–d, Fig. S5a). The spatial
resolution of EFM improves with decreasing $z_{\rm tip}$. For the selected
SCM-PIT-V2 cantilever and $z_{\rm tip}=20$ nm (Fig. 3, main text), the
resolution is approximately $1.1\cdot z_{\rm tip}+11$ nm $=33\,$nm.40

Substituting $F_{\rm es}$ from eq. (S.2) into eq. (S.7) implies a quadratic
dependence of $\Delta\phi_{\rm EFM}$ on $V_{\rm tip}$ reading

$\Delta\phi_{\rm EFM}(x,y)=-\frac{Q(x,y)}{2k}\,\frac{d^{2}C_{\rm
ts}(x,y)}{dz^{2}}\,\left(V_{\rm tip}-V_{\rm CPD}(x,y)\right)^{2},$ (S.8)

where the quality factor $Q$ and the tip sample capacitance $C_{\rm ts}$ are
position dependent as affected by local variations of dissipation and
screening, respectively, due to, e.g., adsorbates or doping fluctuations.
Nevertheless, the maximum of $\phi_{\rm EFM}(V_{\rm tip})$ is a direct measure
of $V_{\rm CPD}$. Indeed, the measured $\phi_{\rm EFM}(V_{\rm tip})$ at a
single point features a parabola as shown by the fit in Fig. 1b, main text.
Minor deviations at larger voltages are most likely caused by the influence of
quantum capacitance that is not captured by eq. (S.2). To determine $V_{\rm
CPD}$ for each position $(x,y)$, we probe $\phi_{\rm EFM}(x,y)$ at only three
$V_{\rm tip}=-1,0,1$ V (Fig. S4f) and use the three points for fitting a
parabolic function

$\phi_{\rm EFM}(x,y)=-\kappa(x,y)\left(V_{\rm tip}-V_{\rm
CPD}(x,y)\right)^{2}+\phi_{0}(x,y)$ (S.9)

via three parameters, namely the curvature $\kappa(x,y)$ and the maximum at
$(\phi_{0}$, $V_{\rm CPD})(x,y)$. Hence, we take the local nature of all three
parameters into account with $\phi_{0}(x,y)$ caused by local forces that are
not dependant on $V_{\rm tip}$. We choose the extrema of $V_{\rm tip}=\pm
1\,$V large enough to obtain high precision of the fit, but still small enough
to avoid deviations from the parabola due to quantum capacitance.

###  S3.1 Noise and Sensitivity

To illustrate the principle, Fig. S5a displays a $V_{\rm CPD}(x,y)$ map of an
electron doped graphene area as obtained by a parabolic fit using eight
$V_{\rm tip}$ for each position $(x,y)$ instead of three, hence, enabling
better precision. Features are discernable at several graphene folds, wrinkles
and point defects (compare with the topography in Fig. 5h, main text). These
topographic features obviously exhibit distinct surface potentials.
Importantly, no long range background had to be subtracted. Figure S5b shows
two selected sets of measured $\phi_{\rm EFM}(V_{\rm tip})$ at the eight
$V_{\rm tip}$ together with a parabolic fit and a 7th order polynomial fit.
The two fits reveal nearly identical maxima.

Figure S5c displays the variation of the curvature of the fitted parabolas
with $V_{\rm gate}$ (main) and with position (inset). While the $V_{\rm gate}$
dependence is rather irregular exhibiting fluctuations of about 5 % only, the
position dependence exhibits features at folds and wrinkles that likely
exhibit different dissipation strengths leading to different local $Q(x,y)$.
Variations in $\Phi_{0}$, related to the forces that do not depend on $V_{\rm
tip}$, are about 0.6∘ for, both, $V_{\rm gate}$ dependence and position
dependence, but without obvious correlations to topographic features (not
shown).

Figure S5: Parabolic fits of $\phi_{\rm EFM}(V_{\rm tip})$. (a) $V_{\rm
CPD}(x,y)$ map of graphene measured by EFM via parabolic fits of eight
recorded $\phi_{\rm EFM}(V_{\rm tip})$ values for each location $(x,y)$ (eq.
(S.9)), $V_{\rm gate}=100\,V,\,V_{\rm SD}=0\,V,\,z_{\rm tip}=30\,$nm. The
corresponding topography is shown in Fig. 5h, main text. (b) Variation of
$\phi_{\rm EFM}$ with $V_{\rm tip}$ for the two marked regions in (a). The
solid/dashed lines are fits to the measured data points with the specified
polynomial order. The black curves are shifted downwards for better
visibility. (c) Curvature $\kappa(V_{\rm gate})$ resulting from parabolic fits
as in (b) (eq. S.9) and displayed for a single location. Inset: $\kappa(x,y)$
of the same area as in (a), $V_{\rm gate}=100\,V,\,V_{\rm SD}=0\,V,\,z_{\rm
tip}=30\,$nm.

For $V_{\rm CPD}(x,y)$ mapping by EFM, we recorded only three $\phi_{\rm
EFM}(V_{\rm tip})$ at each location that are fitted by a parabola (eq. (S.9)).
The resulting noise is estimated in the following. Continuous sweeps of
$\phi_{\rm EFM}(V_{\rm tip})$ consisting of 512 points and recorded at
different $V_{\rm gate}$ revealed an average RMS deviation from the parabola
of $\delta\phi_{\rm EFM}=0.09\degree$ at a recording time per point of 1.5 ms.
The RMS deviation with respect to the parabola dropped continuously, if one
employs Gaussian averaging of $\phi_{\rm EFM}(V_{\rm tip})$ prior to fitting,
as expected. Moreover, fits of higher polynomial order only slightly decreased
$\delta\phi_{\rm EFM}$, e.g., by less than 5% up to a fit order of ten. This
corroborates that the parabola is a very good approximation of $\phi_{\rm
EFM}(V_{\rm tip})$. Since recording times of $3.4\,{\rm ms}$ are employed for
the three $\phi_{\rm EFM}(V_{\rm tip})$ used for $V_{\rm CPD}(x,y)$ maps, the
RMS error becomes $\delta\phi_{\rm
EFM}=0.09\degree\times\sqrt{1.5/3.4}=0.06\degree$. We used this
$\delta\phi_{\rm EFM}$ to deduce the resulting error in $V_{\rm CPD}$ by
numerical simulations revealing $\delta V_{\rm CPD}\approx 2$ mV at $V_{\rm
CPD}\approx 0.2$ V such as in Fig. 3, main text. This largely fits with the
uncorrelated noise that we observe in maps experimentally (section S4.1).

##  S4 Electric Fields due to Current Flow

To map the current induced electric fields, we measure the change of the
contact potential difference $V_{\rm CPD}(x,y)$ after applying a source-drain
voltage $V_{\rm SD}$. This is dubbed the transport voltage $V_{\rm
transport}(x,y,V_{\rm SD})$ defined as.

$V_{\rm transport}(x,y,V_{\rm SD})=V_{\rm CPD}(x,y,V_{\rm SD})-V_{\rm
CPD}(x,y,V_{\rm SD}=0\,V)$ (S.10)

The negative spatial gradient of $V_{\rm transport}(x,y,V_{\rm SD})$ is the
transport induced electric field reading

$\mathbf{E^{\rm meas}}(x,y)=-\nabla V_{\rm transport}(x,y).$ (S.11)

The derivation implicitly assumes that the static work-function fluctuations
of the two subtracted $V_{\rm CPD}(x,y)$ maps remain unchanged and therefore
cancel (section S6).

Figure S6a shows a $V_{\rm transport}(x,y)$ map derived from KPFM maps of
$V_{\rm CPD}(x,y)$ according to eq. (S.10) and after normalizing to $V_{\rm
SD}$. The image employs $V_{\rm SD}=0.5$ V at charge neutrality. The $x$ and
$y$ components of the resulting in-plane electric field $\mathbf{E^{\rm
meas}}(x,y)$ (eq. (S.11)) are displayed in Fig. S6b and c, respectively.
Figure S6g shows a vector representation of $\mathbf{E^{\rm meas}}(x,y)$
represented by arrows on the colored background of the deduced charge carrier
density $n(x,y)$ as present at $V_{\rm SD}=0.5$ V (eq. (S.14)).

Figure S6d-f and h show the same data, but deduced from the $V_{\rm CPD}(x,y)$
maps after 2nd order polynomial background subtraction (section S2.2). This
naturally removes the constant electric field along the applied $V_{\rm SD}$
and, thus, highlights deviations from this average electric field, but it
prohibits the direct identification of inverted electric fields. Correlations
between doping and electric field get more apparent, e.g., hole (electron)
doped regions in this particular area tend to reduce (increase) the electric
field in $x$ direction with respect to the average field. Moreover, it gets
more easy within the vector maps with background subtraction to identify
electric fields pointing in different $y$ directions (Fig. S6g). Finally, one
can identify source like areas (e.g. upper left area) or sink like areas (e.g.
upper fold area) of the current induced electric fields.

Figure S6: Current Induced Electric Fields after Background Subtraction. (a)
$V_{\rm transport}(x,y)/V_{\rm SD}$ map using KPFM for recording $V_{\rm
CPD}(x,y)$ at $V_{\rm SD}=0.5$ V and $V_{\rm SD}=0$ V as input for eq. (S.10),
$V_{\rm gate}-V_{\rm D}=-5\,$V, corresponding topography: Fig. 1f, main text.
(b, c) Maps of $x$ and $y$ component of the electric field, respectively, as
deduced by eq. (S.11) using a. Average electric fields $\left<E_{x,y}\right>$
and the RMS standard deviations $\sigma_{\rm RMS}$ are marked below the
images. The latter is determined within boxes of ($1\,\mu$m)2 and averaged
over all boxes afterwards. (d) Same as a, but after removal of a 2nd order
polynomial background from the two constituting $V_{\rm CPD}(x,y)$ maps
(section S2.2). (e, f) Maps of $x$ and $y$ component of the electric field
derived from d via eq. (S.11). (g) Vector plot of in-plane electric field
$\mathbf{E^{\rm meas}}(x,y)$ deduced via eq. (S.11) from a. The arrows are
colored blue (red), if pointing forward (backward) with respect to the applied
$V_{\rm SD}$. The background color shows the doping distribution $n(x,y)$ at
$V_{\rm SD}=0.5\,$V (eq. (S.14)). (h) Same as g, but deduced from d.

###  S4.1 Noise Filtering

Performing a nearest neighbor differentiation of $V_{\rm transport}(x,y)$ to
obtain $\mathbf{E^{\rm meas}}(x,y)$ maps (eq. S.11) revealed uncorrelated
noise for both components with strength $\delta E^{\rm meas}_{x}\approx\delta
E^{\rm meas}_{y}\approx 100\,{\rm mV}/\mu$m (Fig.S7a). This uncorrelated noise
obscured the observation of any feature. This is likely of electric origin at
the piezoelectric actuators reducing the accuracy to determine $V_{\rm CPD}$
in EFM (section S3).

To reduce the overwhelming noise, we employ two filtering processes. On one
hand, we apply a spatial Gaussian averaging with full width at half maximum
$\Gamma$ to $V_{\rm transport}(x,y)$ maps. On the other hand, we use locations
that are further apart from each other for the differentiation, i.e.,

$\displaystyle E_{x}^{\rm meas}(x,y)=\frac{V_{\rm transport}(x,y)-V_{\rm
transport}(x-\Delta x,y)}{\Delta x}$ $\displaystyle E_{y}^{\rm
meas}(x,y)=\frac{V_{\rm transport}(x,y)-V_{\rm transport}(x,y-\Delta
y)}{\Delta y}$ (S.12)

with $\Delta x$ and $\Delta y$ being the chosen distances in the respective
directions. The RMS fluctuations $\delta E_{x}^{\rm meas}$ and $\delta
E_{y}^{\rm meas}$ scale with $1/\Delta x$ and $1/\Delta y$, respectively,
according to error propagation. This naturally improves the signal to noise
ratio at the expense of spatial resolution.

This second approach is visualized for an experimentally measured $V_{\rm
transport}(x,y)$ map in Fig. S7a–c. Stable structures appear at $\Delta x\geq
0.91$ $\mu$m, i.e., using a distance of 35 measurement points for $\Delta x$.
The remaining uncorrelated noise exhibits a standard variation
$\sigma_{E_{x}}=5$ mV/$\mu$m as deduced by analyzing multiple boxes of
($0.3\,\mu$m)2, much smaller than the apparent features sizes. This can be
compared with the error $\delta V_{\rm CPD}=2$ mV resulting from the
$\Phi_{\rm EFM}(V_{\rm tip})$ noise (section S3.1). Error propagation implies
$\delta V_{\rm transport}=3$ mV (eq. (S.10)) and, respectively, $\delta
E_{x}^{\rm meas}=5$ mV/$\mu$m for $\Delta x=0.91\,\mu$m (eq. S.12) in very
good agreement with the measured noise of 5 mV/$\mu$m.

Figure S7g shows the same procedure of filtering via large $\Delta x$ for a
fictitious electric field $E_{x}(x)$ featuring a dipolar structure (black
line) that is mixed with uncorrelated noise in $V_{\rm transport}(x)$ of
strength $\delta V_{\rm transport}=16$ mV prior to using different $\Delta x$
to determine $E_{x}^{\rm meas}(x)$ according to eq. (S.12) (colored lines).
Obviously, the feature width and height are barely changed by the $\Delta x$
filtering, but the dipolar feature is slightly shifted to the right. In
contrast, the Gaussian averaging makes the features wider and weaker in
amplitude, while maintaining its center position (Fig. S7h). This is also
visible in the accordingly Gaussian filtered experimental images (Fig. S7d-f).

Since our main interest is the feature size and the feature strength, in
particular, during the analysis of Fig. 3, main text, we optimize $\Delta x$
(eq. S.12) towards the lowest possible signal/noise ratio, where features get
significantly stronger than the noise floor, but use only a mild Gaussian
averaging with $\Gamma\leq 100\,$nm. Table S1 summarizes the chosen $\Delta x$
and $\Gamma$ for all electric field maps presented in main text and
supplement.

Eventually, Fig. S7i and j show the rms electric field noise
$\sigma_{E_{x}}(\Delta x)$, determined within multiple boxes of size
$(0.3\,\mu$m)2 (average of all boxes within one image) for different $V_{\rm
gate}$ and $V_{\rm SD}$, respectively. The fitted slope in the double-
logarithmic plots (lines) consistently reveals $\sigma_{E_{x}}\propto\Delta
x^{-1}$, as expected from the discussion above. This evidences uncorrelated
$E_{x}^{\rm meas}(x,y)$ noise. In line, the correlation length of $E_{x}^{\rm
meas}(x,y)$ at smallest possible $\Delta x=26$ nm (eq. S.12) is $\xi\approx
30\,$nm only, i.e., the image resolution. Note that the noise depends barely
on $V_{\rm SD}$ and $V_{\rm gate}$.

Figure | Method | $\Gamma$ | $\Delta x$ | $\Delta y$
---|---|---|---|---
1d | KPFM | 92 nm | 1.17 $\mu$m | 0.12 $\mu$m
1g | KPFM | 81 nm | 0.43 $\mu$m | 0.12 $\mu$m
2f-i | KPFM | 93 nm | 0.39 $\mu$m | 0.08 $\mu$m
3c–d, | EFM | 86 nm | 0.47 $\mu$m | 0.03 $\mu$m
3e–f, S8m–n, | EFM | 62 nm | 0.91 $\mu$m | 0.03 $\mu$m
4 | KPFM | 31 nm | 0.39 $\mu$m | 0.05 $\mu$m
5a–f, S6b–c, S14 | KPFM | 61 nm | 0.52 $\mu$m | 0.12 $\mu$m
5g | KPFM | 61 nm | 0.52 $\mu$m | 0.58 $\mu$m
5i–m | EFM | 77 nm | 0.39 $\mu$m | 0.04 $\mu$m
5n–s | KPFM | 77 nm | 0.82 $\mu$m | 0.05 $\mu$m
S6e–f | KPFM | 81 nm | 0.52 $\mu$m | 0.93 $\mu$m
S6g–h | KPFM | 61 nm | 0.52 $\mu$m | 0.93 $\mu$m
S10c | KPFM | 77 nm | 0.33 $\mu$m | 0.05 $\mu$m
S10d | KPFM | 99 nm | 0.03 $\mu$m | 0.05 $\mu$m
Table S1: Table: Filter parameters for all electric field maps of main text
and supplement. Figure S7: Filtering procedure for electric field maps. (a–c)
Experimentally measured $\widehat{E}_{x}^{\rm meas}(x,y)$ (EFM) using eq.
(S.12) with different values of $\Delta x$ as indicated after mild Gaussian
averaging of $V_{\rm transport}(x,y)$ with $\Gamma=61.5$ nm, the upper left
area of b corresponds to Fig. 3e, main text. (d–f) $\widehat{E}_{x}^{\rm
meas}(x,y)$ using the same $V_{\rm transport}(x,y)$ map as in a–c, but using
Gaussian averaging with different $\Gamma$ before employing $\Delta x=26$ nm
for differentiation (eq. (S.12)). (g) Simulated dipolar feature of
$\widehat{E}_{x}(x)$ (black line) and the same feature after superposition
with white noise of strength $\delta V_{\rm transport}=16\,$mV in the
constituting $V_{\rm transport}(x)$ at $V_{\rm SD}=0.1\,$V and subsequent
application of eq. (S.12) with the indicated $\Delta x$ (colored lines). (h)
Same as g concerning simulated dipole (black line) and noise superposition,
but using subsequently Gaussian averaging with different $\Gamma$ as marked
(colored lines). (i,j) RMS noise of experimental $E_{x}^{\rm meas}(x,y)$ maps,
$\sigma_{E_{x}}(\Delta x)$, for different $V_{\rm gate}$ (i) and $V_{\rm SD}$
(j) (KPFM). The noise for each $\Delta x$ is deduced within boxes of
$0.3\,\mu{\rm m}\times 0.3\,\mu$m covering the whole image. Subsequently, the
noise of all boxes is averaged. The solid lines are linear fits revealing
$\sigma_{E_{x}}\propto\Delta x^{-1}$.

###  S4.2 Temporal Stability of Dopant Distribution

The transport potential $V_{\rm transport}(x,y)$ according to eq. (S.10)
quantifies the potential drop due to charge flow, if the work-function
fluctuations in the measured region remain the same while acquiring the two
$V_{\rm CPD}(x,y)$ maps, biased and unbiased. However, the ambient conditions
during recording partially lead to local charging and discharging depending on
the history of $V_{\rm gate}$ and $V_{\rm SD}$. An example is shown in Fig.
S8a–c, where changes in $V_{\rm CPD}(x,y)$ appear, albeit the maps are
recorded subsequently without changing $V_{\rm gate}$ or $V_{\rm SD}$ and a
full day after setting these voltages. However, changes after recording Fig.
S8c are much less pronounced (Fig. S8d–f), albeit $V_{\rm SD}$ was cycled in
between indicating that a first imaging itself leads to equilibration of the
doping pattern.

Hence, it is crucial to reduce the uncontrolled charging processes. Therefore,
we firstly monitored the two-terminal $R_{\rm transport}$ continuously. It
changed minimally after ramping $V_{\rm gate}$ to about $100\,$V and then ramp
it down slowly by $\sim 0.1\,V/s$ until the device resistance is maximized
signalling charge neutrality. Figure S8g shows a time trace of $R_{\rm
transport}$ after such stabilization revealing only small fluctuations of
about 5 % after a waiting time of roughly 5 hours. Using this procedure,
stability at the local scale has still to be ensured. This was more involved
and not always successful. There is no direct way to experimentally map the
doping distribution, while applying a finite $V_{\rm SD}$. Hence, we checked
the doping distribution by recording $V_{\rm CPD}(x,y,V_{\rm SD}=\,0\,$V) maps
before and after recording $V_{\rm CPD}(x,y,V_{\rm SD}\neq\,0\,$V) maps.
Figure S8h–n display characteristic features of a successful image sequence,
namely a minimum change in $V_{\rm transport}$ while cycling $V_{\rm SD}$ and
very similar spatially averaged $V_{\rm CPD}$ values for images recorded at
the same $V_{\rm SD}$ but after a distinct $V_{\rm SD}$ history (Fig. S8h).
Obviously, the first two images exhibit a rather different $\left<V_{\rm
CPD}(x,y)\right>$ in line with the changes in Fig. S8a–c, implying again that
the first maps within a certain area change the lateral doping distribution
more strongly, likely via the applied tip voltages. However, stable subsequent
images can be often recorded afterwards. Consequently, the $V_{\rm
transport}(x,y)$ maps deduced from a stable sequence are nearly identical, if
distinct reference images recorded at $V_{\rm SD}=\,0\,$V are employed (eq.
(S.10)). This is visible by comparing Fig. S8i and k as well as Fig. S8j and
l, that used image 6 and image 10 as reference, respectively. The similarity
naturally also applies for the resulting electric field maps in Fig. S8m–n.
These kind of images, selected by adequate monitoring, are, hence, reliably
attributed to consequences of the applied $V_{\rm SD}$ (discussion of Fig. 3,
main text).

Figure S8: Temporal Stability of Doping Distribution. (a–f) Subsequently
recorded $V_{\rm CPD}(x,y)$ maps using EFM, $V_{\rm SD}=0\,$V, $V_{\rm
gate}=V_{\rm D}=85\,$V, $z_{\rm tip}=20\,$nm, $f_{\rm drive}=61.7\,$kHz,
$V_{\rm tip}=2\,$V during tapping mode traces, recording time per image: 34
min, waiting time prior to each image: $>5$ min. $V_{\rm gate}$ and $V_{\rm
SD}$ was changed one full day prior to the acquisition of a. Before recording
d (f), $V_{\rm SD}$ was cycled to $+0.2$ V ($-0.2$ V) staying there for 10
min. (g) Two terminal resistance $R_{\rm transport}$ recorded between the
setting of $V_{\rm gate}$ and the data acquisition of Fig. 3, main text,
$V_{\rm gate}=V_{\rm D}=85\,$V. (h) Spatially averaged $V_{\rm CPD}$ (larger
dots, left axis) and simultaneously acquired temporarily averaged device
resistance (smaller dots connected by red line, right axis) for a series of
$V_{\rm CPD}(x,y)$ maps measured in order as numbered, just after finishing
the time trace in g. Same colors belong to the same $V_{\rm CPD}(x,y)$ map.
The current induced electric field map deduced from image 9 and image 6 (10)
is also shown in Fig. 3e (f), main text. (i–l) $V_{\rm transport}(x,y)/V_{\rm
SD}$ maps as deduced from the subtraction of the $V_{\rm CPD}(x,y)$ maps
labelled in h as marked below each image (eq. (S.10)). Note that different
reference maps at $V_{\rm SD}=0$ V are used for i, j and k, l. (m, n) Deduced
$\widehat{E}_{x}^{\rm meas}(x,y)$ map (eq. (S.11)) from j and l, respectively,
using the $V_{\rm CPD}(x,y)$ maps as marked below the image. The rectangular
area marked in l–n is displayed in Fig. 3e, f main text

##  S5 Doping Distribution from $V_{\rm CPD}(x,y)$ Maps

Spatial fluctuations of the contact potential difference between tip and
graphene are related to doping fluctuations of the graphene 48. The
corresponding doping density $n_{0}(x,y)$ can be deduced via $|e|\left(V_{\rm
CPD}(x,y)-V_{\rm CPD}^{0}\right)=E_{\rm F}-E_{\rm D}(x,y)$, where $E_{\rm
F}:=0\,$eV is the Fermi energy of graphene, $E_{\rm D}(x,y)$ is the local
Dirac point energy, and $V_{\rm CPD}^{0}$ is the contact potential difference
between the tip and charge neutral graphene. The resulting doping distribution
$n_{0}(x,y)$ at $V_{\rm SD}=0$ V reads

$n_{0}(x,y)=\frac{e^{2}}{\pi}\,{\rm sign}{\left[V_{\mathrm{CPD}}(x,y,V_{\rm
SD}=0\,{\rm
V})-V_{\mathrm{CPD}}^{\mathrm{0}}\right]}\,\left(\frac{V_{\mathrm{CPD}}(x,y,V_{\rm
SD}=0\,{\rm V})-V_{\mathrm{CPD}}^{\mathrm{0}}}{\hbar
v_{\mathrm{F}}}\right)^{2},$ (S.13)

whith the Fermi velocity of graphene $v_{\rm F}=1\times 10^{6}\,$m/s. Thus,
one needs $V_{\rm CPD}^{0}$. It can be deduced from $V_{\rm CPD}(x,y,V_{\rm
SD}=0\,$V) maps recorded by EFM at charge neutrality, i.e., at $V_{\rm
gate}=V_{\rm D}$, but not by KPFM.

For $V_{\rm CPD}(x,y)$ maps acquired by EFM, we then construct doping maps by
applying eq. (S.13) straightforwardly. However, for KPFM, we have to consider
that $V_{\rm CPD}(x,y)$ maps are obtained after subtracting a second order
polynomial background (section S2.2) that largely removes the spatial average
of $V_{\rm CPD}(x,y)$ (Fig. S9g). In order to restore the average, we
calculate the average contact potential difference for a particular $V_{\rm
gate}$ using $V_{\rm CPD}=V_{\rm CPD}^{0}+{\rm sign}(n_{\rm gate})\,\hbar
v_{\rm F}\sqrt{\pi|n_{\rm gate}|}$ with $n_{\rm gate}=C_{\rm gate}\left(V_{\rm
gate}-V_{\rm D}\right)/e$ as doping density induced by the gate via a
capacitive model (Fig. S9h). Subsequently, the histograms of the $V_{\rm
CPD}(x,y)$ maps (Fig. S9a–c), i.e., each value of the map, are shifted such
that the histogram maxima (Fig. S9g) are aligned with the calculated values of
Fig. S9h. These adjusted $V_{\rm CPD}(x,y)$ maps are then used to apply eq.
(S.13) resulting straightforwardly in the doping maps $n_{0}(x,y)$. Figure
S9d-f displays resulting doping maps of the same area at different $V_{\rm
gate}$. The whole map changes from hole doping to electron doping due to the
added $V_{\rm CPD}^{0}$. Moreover, folds and wrinkles (topography in Fig. 5h,
main text) charge less via gating, which is likely due to the larger distance
from the gate.

Figure S9: Doping Maps.(a–c) Histograms of $V_{\rm CPD}(x,y,V_{\rm SD}=\,$0 V)
as recorded by KPFM after background subtraction (section S2.2), $V_{\rm
gate}-V_{\rm D}$ as marked, $V_{\rm D}=81\,$V, red line in b: deduced $V_{\rm
CPD}^{0}=-1.3$ mV. (d–f) Doping maps $n_{0}(x,y)$ corresponding to the
adjacent histograms and deduced from $V_{\rm CPD}(x,y,V_{\rm SD}=\,$0 V) via
eq. (S.13) using the capacitively calculated background shift of (h),
corresponding topography: Fig. 5h, main text, black rectangle in e: region of
Fig. 4a–f, main text. (g) Maximum of $V_{\rm CPD}$ histograms for various
$V_{\rm gate}$ deduced from KPFM images after background subtraction. (h)
$\left<V_{\rm CPD}\right>$ as derived from a capacitive gating model ($V_{\rm
D}=81\,$V, text) and as added to the background subtracted $V_{\rm CPD}(x,y)$
maps prior to applying eq. (S.13) for calculating $n_{0}(x,y)$.

###  S5.1 Source Drain Induced Local Doping (SDILD)

As discussed in the main text, the doping distribution $n(x,y)$ is subject to
a linear gradient at finite $V_{\rm SD}$. The gate voltage reference at the
graphene is ground for the drain electrode and shifted by $V_{\rm SD}$ with
respect to ground for the source electrode (Fig. 2a, main text). Hence, we
define a local gate voltage $V_{\rm gate}^{\rm local}(x,V_{\rm SD})=V_{\rm
gate}-\frac{V_{\rm SD}}{L}(x-x_{0})$, where $x=x_{0}$ is the edge of the drain
electrode and $x=0$ the left side of the image (closest to the drain). We will
use $V_{\rm gate}^{\rm local}(x,V_{\rm SD})$ to calculate the resulting
$n(x,y)$ iteratively. Firstly, we regard the doping $n_{0}(x,y)$ at $V_{\rm
SD}=0$ V. Its fluctuations change the local charge neutrality point at finite
$V_{\rm gate}$ reading $V_{\rm D}^{\rm local}(x,y)=V_{\rm
gate}-e\frac{n_{0}(x,y)}{C_{\rm eff}(x,y,V_{\rm SD}=0\,{\rm V})}$. Here,
$C_{\rm eff}(x,y,V_{\rm SD}=0\,{\rm V})=\left(\frac{1}{C_{\rm
gate}}+\frac{1}{C_{\rm Q}(x,y,V_{\rm SD}=0\,{\rm V})}\right)^{-1}$ is the
effective capacitance consisting of geometric capacitance $C_{\rm gate}$ and
quantum capacitance $C_{\rm Q}(x,y,V_{\rm SD}=0\,{\rm
V})=\frac{2e^{2}}{\pi}\frac{\sqrt{\pi|n_{\rm 0}(x,y)|}}{\hbar v_{\rm F}}$.

In a second step, we set up the iterative loop for $n(x,y)$ at finite $V_{\rm
SD}$ reading

$\displaystyle n(x,y,V_{\rm SD})$ $\displaystyle=\frac{C_{\rm eff}(x,y,V_{\rm
SD})}{e}\left(V_{\rm gate}^{\rm local}(x,V_{\rm SD})-V_{\rm D}^{\rm
local}(x,y)\right)$ $\displaystyle=\frac{C_{\rm eff}(x,y,V_{\rm SD})}{C_{\rm
eff}(x,y,V_{\rm SD}=0\,{\rm V})}n_{0}(x,y)-\left(\frac{V_{\rm
SD}}{eL}\right)C_{\rm eff}(x,y,V_{\mathrm{SD}})(x-x_{0}),$ (S.14)

where $C_{\rm eff}(x,y,V_{\rm SD})$ is the effective capacitance for the
changed doping $n(x,y,V_{\rm SD})$ via $V_{\rm SD}$. Thus, both sides of the
equation depend on $n(x,y,V_{\rm SD})$ suggesting a self-consistent loop that
we applied for solving. Far away from charge neutrality, the second term is
negligible, since $|V_{\rm SD}|\ll V_{\rm gate}$ and $C_{\rm eff}(x,y,V_{\rm
SD})\approx C_{\rm gate}$ implying $n(x,y)\approx n_{0}(x,y)$, i.e., SDILD is
negligible. However, when $V_{\rm gate}\rightarrow V_{\rm D}$, the first term
is almost zero raising the importance of the second term. Consequently,
$V_{\rm SD}$ induced gating (SDILD) has to be considered carefully, in
particular, close to charge neutrality.

##  S6 Electric Field from SDILD

In section S4, we assigned the difference between biased and unbiased $V_{\rm
CPD}(x,y)$ maps to the current induced voltage drop $V_{\rm transport}(x,y)$
(eq. (S.10)) using the assumption that the doping distribution in the sample
remains unchanged. However, section S5.1 reveals that $V_{\rm SD}$ changes the
doping via SDILD (eq. (S.14)), most strongly close to charge neutrality. Even
for a linear potential drop along the sample by $V_{\rm SD}$, we get non-
linearities in $V_{\rm transport}(x,y)$ (eq. (S.10)) via the non-linear local
doping due to quantum capacitance (eq. (S.14)). We dub the resulting electric
field via SDILD $\mathbf{E}^{\rm SDILD}(x,y)$ (eq. S.11). Here, we present a
formalism to calculate it for known $V_{\rm SD}$ and $n_{0}(x,y)$ (eq. S.13).
A comparison of $\mathbf{E}^{\rm SDILD}(x,y)$ with $\mathbf{E}^{\rm
meas}(x,y)$ then allows to distinguish between non trivial transport features
and those due to SDILD.

The Fermi energy of graphene reads $E_{\rm F}-E_{\rm D}={\hbar v_{\rm F}\;\rm
sign}(n)\sqrt{\pi|n|}$ with $n(x,y,V_{\rm SD})$ being the doping distribution
including SDILD (eq. (S.14)). If we assume that the applied $V_{\rm SD}$ drops
linearly along the graphene, implying a potential $\frac{V_{\rm
SD}}{L}(x-x_{0})$ ($x_{0}$: edge of the drain electrode), we obtain

$\displaystyle V_{\rm CPD}(x,y,V_{\rm SD})$ $\displaystyle=V_{\rm
CPD}^{0}+\frac{\left(E_{\mathrm{F}}(x,V_{\rm SD})-E_{\rm D}(x,y,V_{\rm
SD})\right)}{|e|}+\frac{V_{\rm SD}}{L}(x-x_{0})$ $\displaystyle=V_{\rm
CPD}^{0}+\frac{\hbar v_{\rm F}}{|e|}\,{\rm sign}(n(x,y,V_{\rm
SD}))\sqrt{\pi|n(x,y,V_{\rm SD})|}+\frac{V_{\rm SD}}{L}(x-x_{0}).$ (S.15)

The contact potential difference in the unbiased case is (section S5)

$V_{\rm CPD}(x,y,0)=V_{\rm CPD}^{0}+{\frac{\hbar v_{\rm F}}{|e|}\rm
sign}(n_{0}(x,y))\sqrt{\pi|n_{0}(x,y)|}.$

The resulting $x$ component of the electric field (eqs. (S.10), (S.11))
becomes

$\displaystyle E_{x}^{\rm SDILD}(x,y,V_{\rm SD})$
$\displaystyle=-\frac{d}{dx}[V_{\rm CPD}(x,y,V_{\rm SD})-V_{\rm CPD}(x,y,0)]$
$\displaystyle=-\frac{V_{\rm SD}}{L}-\frac{\hbar v_{\rm
F}}{|e|}\frac{d}{dx}[{\rm sign}(n(x,y,V_{\rm SD}))\sqrt{\pi|n(x,y,V_{\rm
SD})|}-{\rm sign}(n_{0}(x,y))\sqrt{\pi|n_{0}(x,y)|}]$
$\displaystyle=E_{x}^{0}-\frac{\hbar v_{\rm
F}\sqrt{\pi}}{2|e|}\left(\frac{1}{\sqrt{|n(x,y,V_{\rm
SD})|}}\frac{dn(x,y,V_{\rm
SD})}{dx}-\frac{1}{\sqrt{|n_{0}(x,y)|}}\frac{dn_{0}(x,y)}{dx}\right),$ (S.16)

using $\frac{d}{dx}\left({\rm
sign}(n)\;\sqrt{\pi|n|}\right)=\frac{\sqrt{\pi}}{2}\frac{1}{\sqrt{|n|}}\frac{dn}{dx}$.
The term $E_{x}^{0}=-V_{\rm SD}/L$ is the trivial electric field by the linear
voltage drop due to $V_{\rm SD}$. Normalising the electric field to $V_{\rm
SD}$, we get

$\displaystyle\widehat{E}_{x}^{\rm SDILD}(x,y,V_{\rm SD})$
$\displaystyle=-\frac{1}{L}-\frac{\beta}{V_{\rm
SD}}\left(\frac{1}{\sqrt{|n(x,y,V_{\rm SD})|}}\frac{dn(x,y,V_{\rm
SD})}{dx}-\frac{1}{\sqrt{|n_{0}(x,y)|}}\frac{dn_{0}(x,y)}{dx}\right)$
$\displaystyle:=\widehat{E}_{x}^{0}-\frac{\beta}{V_{\rm SD}}\chi(x,y,V_{\rm
SD})$ (S.17)

with $\beta=\hbar v_{\mathrm{F}}\sqrt{\pi}/2|e|=6.425\cdot 10^{-10}\,V$m and
$\widehat{E}_{x}^{0}=-1/L\approx-0.055\,\mu{\rm m^{-1}}$.

Obviously, $\widehat{E}_{x}^{\rm SDILD}$ diverges, if $n_{0}(x,y)$ or
$n(x,y,V_{\rm SD})$ crosses zero as naturally appearing close to charge
neutrality. This implies pronounced local maxima and minima in the measured
electric field maps that are not caused by the current flow.

As an example, we consider a one-dimensional parabolic profile $n_{0}(x)$
along $x$ with negative curvature crossing zero twice and being located in the
center between source and drain electrode (black line, Fig. S10a). Applying
$V_{\rm SD}$ results in a vertical and horizontal shift (grey line, Fig.
S10a).The vertical shift is more pronounced than the horizontal one, since the
total length within Fig. S10a is much smaller than the distance to the drain
electrode. Since the drain is grounded, half of the $V_{\rm SD}$ induced
potential drops towards the displayed center region. Naturally, the two zero
crossings of the $n(x)$ parabola are shifted inwards with respect to the zero
crossings of $n_{0}(x)$ such that eq. (S.17) implies four distinct divergences
of $\widehat{E}_{x}^{\rm SDILD}(x)$ along $x$. The direction of divergence,
peak or dip, changes between the zeroes of $n(x)$ and the zeroes of $n_{0}(x)$
at the same slope of the parabola (eq. (S.17)). It also changes with direction
of slope for the same density, either $n(x)$ or $n_{0}(x)$. Consequently, a
quadrupole of $\widehat{E}_{x}^{\rm SDILD}(x)$ divergences appears for the
inverted $n_{0}(x)$ parabola with order dip-peak-dip-peak from left to right
(red line, Fig. S10a). Quantum capacitance barely changes this scenario (full
and dashed grey line, Fig. S10a).

Figure S10b shows a 2D plot of $\widehat{E}_{x}^{\rm SDILD}(x,y)$ for a
rotationally symmetric parabolic doping profile $n_{0}(x,y)$ with the same
apex and curvature as in Fig. S10a. It features four lobes, one for each of
the four divergences that we have already discussed in Fig. S10a. The inner
lobes appear more extended along $x$ than the outer lobes as consequence of
the weaker slopes $dn(x,y)/dx$ at the zeroes of $n(x,y)$ compared to the
stronger slopes $dn_{0}(x,y)/dx$ at the zeroes of $n_{0}(x,y)$.

Figure S10: Electric fields from SDILD. (a) Doping profiles (left axis):
parabolic $n_{0}(x)$ at $V_{\rm SD}=0\,$V (black line) using the equation
above with $\Delta=2.5$ $\mu$m, $\gamma=-8.3\cdot 10^{25}\,{\rm m^{-4}}$ and
resulting $n(x,V_{\rm SD}=0.1\,$V) (grey lines) according to eq. (S.14) either
neglecting quantum capacitance (dashed line) or including it (full line). In
addition, the deduced $\widehat{E}_{x}^{\rm SDILD}(x)$ (red line, right axis)
is shown using eq. (S.17) without quantum capacitance (dashed) or with it
(full). Dashed red line: trivial $\widehat{E}_{x}^{0}=-1/L=-0.055\,\mu{\rm
m^{-1}}$ resulting from the device geometry. (b) 2D color map of
$\widehat{E}_{x}^{\rm SDILD}(x,y,V_{\rm SD}=0.1\,$V) for a rotational
symmetric parabolic doping profile (formula and parameters as in a) located in
the center between source and drain electrode. The trivial field
$\widehat{E}_{x}^{0}$ has been subtracted. (c) $\widehat{E}_{x}^{\rm
SDILD}(x,y,V_{\rm SD})$ map calculated from a measured $V_{\rm CPD}(x,y,V_{\rm
SD}=0\,{\rm V})$ map (AM KPFM) at charge neutrality ($V_{\rm gate}=72\,$V) for
$V_{\rm SD}=-1\,$V. (d) Measured $\widehat{E}_{x}^{\rm meas}(x,y)$ of the same
area as c (same as Fig. 5p, main text). The red arrows in c,d mark features
present in both images and, hence, attributed to SDILD. The black arrow in d
points at an encircled feature not present in c, and, hence, assigned to the
current flow.

Figure S10c displays an $\widehat{E}_{x}^{\rm SDILD}(x,y)$ map deduced from a
measured $V_{\rm CPD}(x,y,V_{\rm SD}=0\,{\rm V})$ map by KPFM using eqs.
(S.13), (S.14) and (S.16). For comparison, the measured $\widehat{E}_{x}^{\rm
meas}(x,y)$ map of the same area is displayed in Fig. S10d. The features that
are similar in both images (red arrows) are attributed to artifacts from
SDILD, while additional features in the $\widehat{E}_{x}^{\rm meas}(x,y)$ map,
such as the encircled yellowish area of inverted electric field, can be
assigned to current induced features.

###  S6.1 Simulation of the Measured Electric Fields in Fig. 3 by SDILD

Figure 3e, f, main text, showcase a dipolar $\widehat{E}_{x}^{\rm meas}(x,y)$
structure as determined by EFM that exhibits a spatial extent and an intensity
of its two lobes much larger than the apparent electric field generated by
SDILD (Fig. 3c, d, main text). This structure is shown again in Fig. S12e.
Since we can not avoid temporal fluctuations of the doping profile
$n_{0}(x,y)$ completely (section S4.2), we also attempted to reproduce the
recorded dipolar $E_{x}^{\rm meas}(x,y)$ structure by an arbitrary, artificial
doping profile $n_{0}(x,y)$. Eventually, we found such a doping profile (black
line, Fig. S12d) that, however, appears to be impossible in reality. We
describe the reasoning in the following after recalling that the doping
profiles recorded prior and after the $V_{\rm CPD}(x,y,V_{\rm SD}=0.1\,{\rm
V})$ map necessary for $E_{x}^{\rm meas}(x,y)$ are, both, not capable to
reproduce $E_{x}^{\rm meas}(x,y)$ via SDILD (Fig. 3, main text).

We firstly consider a one dimensional doping profile $n_{0}(x)$ and later
extend the analysis to two dimensions. To observe a divergence of the electric
field via SDILD, either $n_{0}(x)$ or $n(x)$ must cross zero (eq. (S.17)). As
shown in Fig. S10a–b, $\widehat{E}_{x}^{\rm SDILD}$ rapidly decreases away
from such crossings. Hence, we firstly focus on the area close to such
crossings. If, both, $n_{0}(x)$ and the shifted $n(x)$ cross zero with the
same direction of slope, one gets a dipolar divergence of
$\widehat{E}_{x}^{\rm SDILD}(x)$ (Fig. S10a). Around the crossing point, we
apply a Taylor expansion of $n_{0}(x)$ and $n(x)$, here given for $n_{0}(x)$
around a crossing at $x=0$, to deduce leading terms reading

$n_{0}(x)=\sum_{m}\alpha_{m}\left(\frac{x}{L}\right)^{m},\,\,\,m=1,2,3.....$
(S.18)

Here, $m$ is an integer exponent, $\alpha_{m}$ the corresponding prefactor and
$L$ a constant length. Considering eq. (S.17), we find for the individual
terms of the expansion

$\frac{dn_{0}(x)/dx}{\sqrt{n_{0}(x)}}=\frac{m\sqrt{\alpha_{m}}}{L}\cdot\left(\frac{x}{L}\right)^{m/2-1}.$
(S.19)

Consequently, the linear term of the expansion ($m=1$) produces the
divergence, while the other terms contribute (in first order) by zero ($m>2$)
or $m\sqrt{\alpha_{m}}/L$ ($m=2$) to $\widehat{E}_{x}^{\rm SDILD}(x)$ at the
crossing point. Thus, the width of the prominent SDILD feature around
$n_{0}(x)=0$ is largely dominated by the linear term of the Taylor expansion
(see discussion below for the influence of higher order terms). A resulting
dipolar $\widehat{E}_{x}^{\rm SDILD}(x)$ for a linear $n_{0}(x)$ term, as
calculated numerically via eq. (S.17), is shown in Fig. S11a. It exhibits a
very small width $w_{\rm lobe}$ of its two lobes at the required strength
$|\widehat{E}_{x}^{\rm SDILD}(x)|>0.25/{\mathrm{\mu}}$m that has been found in
the experiment (Fig. S12e). This width $w_{\rm lobe}$ in $x$ direction, where
$|\widehat{E}_{x}^{\rm SDILD}(x)|>0.25/{\mathrm{\mu}}$m, is plotted as a
function of the linear prefactor $\alpha_{1}$ in Fig. S11f (pink lines)
revealing a maximum of $w_{\rm lobe}$ for both lobes (dashed, full line).
Thus, albeit the extremal $|\widehat{E}_{x}^{\rm SDILD}(x)|$ increases
monotonously with $\alpha_{1}$ (Fig. S11e), $w_{\rm lobe}$ is non-monotonous
with a maximum at $w_{\rm lobe}\approx 0.1\,\mathrm{\mu m}$, exactly for the
$n_{0}(x)$ slope as used in Fig. S11a. Thus, a linear zero crossing of
$n_{0}(x)$ can not reproduce the width of the experimentally observed
$\widehat{E}_{x}^{\rm meas}(x)$ dipolar feature with $w_{\rm lobe}\simeq
1.7\,{\mathrm{\mu}}$m. Numerical tests of various shapes around the
$n_{0}(x)=0$ crossing corroborated this result for more general crossings
within reasonable limits.

Figure S11: $\widehat{E}_{x}^{\rm SDILD}(x,y)$ for different exponents of
$n_{0}(x,y)$. (a–d) Simulated doping profiles $n_{0}(x)$ (black lines)
according to the formula on top with $\alpha_{m}=1\cdot 10^{15}$/m2 and $n(x)$
at $V_{\rm SD}=0.1$ V (grey lines) according to eq. (S.14) for different
exponents $m$ as marked (left axis). The resulting electric field from SDILD
$\widehat{E}_{x}^{\rm SDILD}(x)$ according to eq. (S.17) is drawn in red
(right axis). Full (dashed) lines include (do not include) quantum
capacitance. (e) Extremal electric field within the two lobes of the simulated
dipolar structures (as shown in a, c, d) for different $m$. Solid and dashed
lines refer to the lobe with enhanced and inverted electric field,
respectively. (f) Width along $x$ of the simulated two lobes with
$|\widehat{E}_{x}^{\rm SDILD}(x,y)|\geq 0.25\,{\mathrm{\mu}}$m. Large pink dot
marks the width deduced from a. (g) $\widehat{E}_{x}^{\rm SDILD}(x)$ for
different $m$ using $n_{0}(x)$ as given above a, b, $\alpha_{m}=1\cdot
10^{15}$/m2. (h–l) $\widehat{E}_{x}^{\rm SDILD}(x,y)-\widehat{E}_{x}^{0}$ (eq.
(S.17)) for different exponents $m$ as marked using a rotational symmetric
$n_{0}(|\vec{x}|)$ with radial dependence as given above a, b,
$\alpha_{m}=5\cdot 10^{18}$/m2 (b), h–j have the same size, k–l have the same
size. All images use $V_{\rm SD}=0.1$ V.

As alternative, we consider $n_{0}(x)$ functions that are touching
$n_{0}(x)=0$ with zero slope $dn_{0}(x)/dx=0$ (argument applies analogously
for $n(x)$). Again, we apply a Taylor expansion of $n_{0}(x)$ around
$n_{0}(x)=0$. This also adresses the higher order terms of the zero crossings
on the same footing. We realize that the odd powers of the Taylor expansion
with $m>1$ cannot produce a dipolar structure at all (Fig. S11c), while for
$m=1$ the same argument as above applies. As illustrated for $m=5$ in Fig.
S11c, the functions with odd exponent $m>1$ produce a strong divergence of
$\widehat{E}_{x}^{\rm SDILD}(x)$ at the shifted $n(x)=0$ due to its large
slope $|dn/dx|$. In addition, two weaker features appear that are symmetric
around the touching point $n_{0}(x)=0$. These two weaker features originate
from $n_{0}(x)$ and naturally exhibit the same strength on both sides of
$n_{0}(x)=0$ and the same distance from $x=0$ (eq. (S.17)). Thus, each odd
power produces either a monopolar or a tripolar structure, but not a dipolar
one as crosschecked numerically. In contrast, the even powers of $n_{0}(x)$
result in an antisymmetric $\widehat{E}_{x}^{\rm SDILD}(x)$ feature around the
touching point provided that $n(x)$ is moved completely away from $n(x)=0$ via
$V_{\rm SD}$ (Fig. S11b,d). Hence, a dipolar structure naturally appears for
an even power function of $n_{0}(x)$ ($n(x)$) with negative (positive)
curvature in case of the applied positive $V_{\rm SD}$. However, numerically,
it turns out that also the lobes of these dipolar $\widehat{E}_{x}^{\rm
SDILD}(x)$ structures exhibit a maximum width of 0.3 $\mu$m at the required
$\widehat{E}_{x}^{\rm SDILD}>0.25/\mu$m (Fig. S11f). This is again much too
small to reproduce the experiment (Fig. S12e).

However, the numerical analysis of the different exponents reveals that the
extrema of the dipolar (even $m$) and tripolar (odd $m>1$)
$\widehat{E}_{x}^{\rm SDILD}(x)$ structures move continuously outwards with
increasing $m$. This is displayed in Fig. S11g showing $\widehat{E}_{x}^{\rm
SDILD}(x)$ profiles for each exponent $m$, in this case using
$\alpha_{m}=1\cdot 10^{15}$/m2 independent of $m$. The systematic shift to
larger distances with $m$ is due to the increasingly flat $n_{0}(x)$ part in
the center (Fig. S11a–d). While eq. (S.19) indicates a continuously increasing
$\widehat{E}_{x}^{\rm SDILD}(x)$ with increasing $|x|$ and increasing $m$,
this increasing strength due to $n_{0}(x)$ gets increasingly compensated by
the additional $n(x)$ term in eq. (S.17). Basically, the offset between
$n_{0}(x)$ and $n(x)$ within the square root denominators looses its
importance. Very roughly, the extrema appear at the $x$ positions, where
$n_{0}(x)$ has a similar value than $n(x=0)$. This naturally explains that the
extrema move outwards with increasing $m$. Nevertheless, its width $w_{\rm
lobe}$ always remains far below the experimental $w_{\rm lobe}\simeq
1.7\mathrm{\mu}$m (Fig. S11f). In addition, Fig. S11g nicely shows the
alternating tripolar and dipolar structures for increasing $m$.

The extension to 2D isotropic profiles does not change the above arguments as
shown exemplarily in Fig. S11h–l.

Figure S12: Patching doping profiles to reproduce $\widehat{E}_{x}^{\rm
meas}(x,y)$ by SDILD. (a) $\widehat{E}_{x}^{\rm SDILD}(x)$ at $V_{\rm SD}=0.1$
V for different exponents $m$ of the employed $n_{0}(x)$ function given on top
using eqs. (S.14), (S.17). The vertical lines mark $x$ sections that are
patched with a singular exponent $m$ in eq. (S.20) with $m$ increasing
symmetrically from the center to the rims of the $x$ axis. In particular, the
section at $x=0$ exhibits $m=1$, the two adjacent ones $m=2$, the two next
outer ones $m=4$, etc.. (b) Doping distribution $n_{0}(x,y)$ using eq. (S.20)
and the patched $m$ sequence as marked/described in a, but rotational
symmetric. (c–d) $\widehat{E}_{x}^{\rm SDILD}(x,y)$ using $n_{0}(x,y)$ from b,
but with different constant $\beta_{m}$ as marked, according to eqs. (S.14),
(S.17). Black and red curves are cross sections through the doping profiles
$n_{0}(x,y=0)$ and $n(x,y=0,V_{\rm SD})$ along $x$ at the center of the image
(right axis). (e) $\widehat{E}_{x}^{\rm meas}(x,y)$ at $V_{\rm SD}=0.1$ V via
eq. (S.12), (S.10) (same as Fig.3e, main text), same colorbar for c–e.

Since the single components of the Taylor series are not able to reproduce
inverted field lobes of intensity above $0.25\,\mu{\rm m}^{-1}$ and width
$w_{\rm lobe}>1.7\,\mu$m, we have to patch $n_{0}(x)$ piecewise in order to
get a large enough width of the dipolar structure. Since a patched sequence of
crossings would produce a sequence of dipolar structures via their dominating
linear term instead of a single extended dipolar structure, the option of
patching crossing points does not exist.

The fact that the extrema of the dipolar structure for each even $m$ shift
outwards with increasing $m$ (Fig. S12a) suggests a natural way of patching by
choosing $x$ ranges of extremal $\widehat{E}_{x}^{\rm SDILD}$ for each even
$m$ until the extrema position of the largest $m$ matches the $x$ extension of
the experimental dipolar $\widehat{E}_{x}^{\rm meas}(x,y)$ structure. More
formally, we use

$n_{0}(x)=-\beta_{m}\cdot\left(\frac{x}{L}\right)^{m(x)}$ (S.20)

with $m(x)$ being a step function that is a piecewise constant integer along
$x$, featuring subsequently $m=1,2,4,6,8,...$ and using $L=2.5\,\mu$m (Fig.
S12a,b). Increasing the parameter $\beta_{m}$ tunes the intensity of the
dipolar structure (Fig. S11e), but increases the number of required patches
simultaneously by decreasing $w_{\rm lobe}$ (Fig. S11f). Using this
construction, we found that $\sim 50$ patches are required to reproduce the
width and strength of the experimentally observed dipolar structure (Fig.
S12d, e). However, this necessarily requires that each area of different $m$
is adjusted to $n_{0}(x)=0$/m2 in the center (at $x=0$) implying jumps of
$n_{0}(x)$ as displayed by the black line in Fig. S12d. It implies $\sim 50$
jumps in $n_{0}(x)$ on a width of $~{}3\,\mu$m (lines in Fig. S12d). Each jump
returns $n_{0}(x)$ basically back to zero. Such doping profiles are
experimentally very unlikely, in particular, to be present in a certain region
but not in its surrounding.

Importantly, albeit different patching profiles might reproduce the
experimental dipolar structure, the patching method with its multiple jumps
back to $n_{0}(x)\approx 0$ cannot by avoided. Any type of a smooth return to
$n_{0}(x)\approx 0$ would lead to the opposite electric field such that
instead of an extended lobe, one would get multiple dipolar structures within
the lobe area. The resulting requirement of the jumps strongly excludes a
physical possibility that the observed dipolar structure is caused by SDILD,
even in the unlikely case that the doping profile during recording $V_{\rm
CPD}(x,y,V_{\rm SD})$ is substantially different from the one prior and after
the recording.

##  S7 Estimates of Local Scattering Lengths

###  S7.1 Electron Electron Scattering Length $l_{\rm ee}$

The electron-electron scattering length $l_{\rm ee}$ is associated with the
quasiparticle decay rate $\tau_{\rm ee}^{-1}$ due to inelastic electron
electron scattering via $l_{\rm ee}=v_{\rm F}\tau_{\rm ee}$. This is directly
related to the imaginary part of the retarded self energy $\Sigma_{s}^{\rm
R}(\mathbf{k},\xi_{\rm\mathbf{k},s})$ of monolayer graphene as 71

$\frac{1}{\tau_{\rm ee}}=2\,{\rm Im}\left[\Sigma_{s}^{\rm
R}(\mathbf{k},\xi_{\rm\mathbf{k},s}),\right]$ (S.21)

where $s=+/-$ is the conduction/valence band index and
$\xi_{\rm\mathbf{k},s}=s\hbar v_{\rm F}|\mathbf{k}|-\mu$ is the single
particle band energy relative to the non-interacting chemical potential $\mu$
at finite temperature. Since graphene on SiN is a weakly correlated material
(interaction parameter $r_{\rm s}=0.51$)54, the $GW$ approximation is adequate
50 leading to 72

$\displaystyle{\rm Im}[\Sigma_{s}^{\rm
R}(k,\omega)]=\sum_{\mathbf{q},s^{\prime}=\pm}V_{q}\left[n_{\rm
B}\left(\xi_{\mathbf{k+q},\,s^{\prime}}-\hbar\omega\right)+n_{\rm
F}\left(\xi_{\mathbf{k+q},\,s^{\prime}}-\hbar\omega\right)\right]$
$\displaystyle\times\left(1+ss^{\prime}\cos{\theta}\right){\rm
Im}\left[\frac{1}{\epsilon\left(q,\xi_{\mathbf{k+q},\,s^{\prime}}/\hbar-\omega\right)}\right],$
(S.22)

where $V_{q}=\frac{2\pi e^{2}}{\kappa q}$ is the Coulomb interaction for
momentum transfer $q$, $\kappa=4.25$ is the background dielectric constant
combining SiN and vacuum, $n_{\rm B}$ and $n_{\rm F}$ are the Bose and Fermi
distribution functions, respectively, and $\theta$ is the angle between
electron wave vectors $\mathbf{k}$ and $\mathbf{k+q}$. The finite temperature
dynamic dielectric function $\epsilon(q,\omega)$ reads within Random Phase
Approximation (RPA) $\epsilon(q,\omega)=1+V_{q}\Pi(q,\omega)$, where
$\Pi(q,\omega)$ is the irreducible polarizability.

In the limit of small $\xi_{\rm k}$ and low temperature, $\xi_{\mathbf{k}}\ll
k_{\rm B}T\ll E_{\rm F}$, where $E_{\rm F}=\hbar v_{\rm F}\cdot{\rm
sign}(n)\sqrt{\pi|n|}$ is the Fermi energy for charge carrier density $n$, an
asymptotic form of eq. (S.22) has been derived by Li et al. 55 (eq. (10),
there), that by substituting into eq. (S.21) results in an electron-electron
scattering length $l_{\rm ee}$ 51, 52, 26:

$l_{\rm ee}=\frac{4}{\pi}\left(\frac{\hbar v_{\rm F}}{k_{\rm
B}T}\right)^{2}\sqrt{\pi|n|}\frac{1}{\ln{\left(\frac{2E_{\rm F}}{k_{\rm
B}T}\right)}}$ (S.23)

as also given as eq. (4), main text. For the gate voltages of Fig. 4g, main
text, where $l_{\rm ee}(V_{\rm gate})$ is plotted for $T=298\,$K, we have
$k_{\rm B}T/E_{\rm F}\in[0.08-0.24]$, except at charge neutrality, such that
the required limit of eq. (S.23) applies.55.

The formalism is, moreover, not applicable in the quantum critical regime that
appears at room temperature for carrier concentrations $n<1.5\times
10^{14}\,{\rm m^{-2}}$ 53, far lower than $n$ for all data points in Fig. 4g,
main text, and Fig. S14e, except of charge neutrality. But even at charge
neutrality, the lateral charge fluctuations $\delta n_{\rm dis}\approx 1\times
10^{15}\,{\rm m^{-2}}$ (puddles) are one order of magnitude larger than the
threshold for quantum criticality. Finally, the threshold carrier
concentration $n_{\rm thermal}$ above which eq. (S.23) is valid ( $E_{\rm
F}(n_{\rm thermal})=k_{\rm B}T$) is $n_{\rm thermal}\approx 5\times
10^{14}\,{\rm m^{-2}}$, i.e., larger than the quantum critical threshold, but
lower than $\delta n_{\rm dis}$.

###  S7.2 Electron Disorder Scattering Length $l_{\rm dis}$

Figure S13: Disorder Scattering Length. (a) Map of local resistivity
$\rho_{\rm loc}(x,y)=\left<E_{x}^{\rm meas}(x,y,V_{\rm gate},V_{\rm
SD}=0.4\,{\rm V})\right>_{V_{\rm gate}}/(I_{\rm SD}/W)$ with $E_{x}^{\rm
meas}(x,y)$ averaged for $V_{\rm gate}=-141--91\,$V (rectangle in b), $I_{\rm
SD}$: source-drain current during imaging, $W$: sample width, respective
topography: Fig. 5h, main text. Black rectangles mark the spatial averaging
areas for the data in b. (b) Spatially averaged $\overline{\rho}_{\rm
loc}(V_{\rm gate})=\left<E_{x}^{\rm meas}(x,y,V_{\rm SD}=0.4\,{\rm V},V_{\rm
gate})\right>_{x,y}/(I_{\rm SD}/W)$. Full colored lines: areas marked in a,
black line: whole region of a, red dashed line: two-terminal device
resistivity, $L$: sample length. (c) Simulated $V_{\rm gate}$ dependence (eq.
S.24) of the electron-phonon scattering length $l_{\rm e-ph}$ due to gate
independent resistivity $\rho_{\rm A}=30\,{\Omega}$ from longitudinal
accoustic phonon scattering73, 55 (blue) and the complete electron disorder
scattering length $l_{\rm dis}$ deduced from the measured $\rho_{\rm
local}=500\,{\Omega}$ at large $|V_{\rm gate}|$ (dashed black) as well as the
resulting short range disorder scattering length $l_{\rm SR}$ (magenta)
employing the Matthiesen rule shown as inset. The upper electron density axis
($n$) in b and c is deduced from a capacitive model.

All scattering processes of electrons with external perturbations as point
defects, Coulomb type impurities, sample edges, or phonons change the momentum
of the electron system. We dub the respective scattering length $l_{\rm dis}$
(disorder scattering length). A straightforward approach to determine $l_{\rm
dis}$ uses the resistivity of graphene at high carrier densities where
electron-electron scattering is negligible (eq. (S.23)), such that the
effective mean free path reading $l_{\rm MFP}=\left(l_{\rm ee}^{-1}+l_{\rm
dis}^{-1}\right)^{-1}$ is dominated by $l_{\rm dis}$.

Although we are not able to directly map the resistivity with KPFM or EFM, the
recorded electric field maps at high doping give a reasonable estimate
assuming homogeneous current density $I_{\rm SD}/W$ ($I_{\rm SD}$: source-
drain current, $W$: sample width) in first order, since non-trivial viscous
properties disappear at large $|V_{\rm gate}|$ (Fig. 5a, main text). This
allows to approximately map local variations in resistivity $\rho_{\rm
loc}(x,y)=\frac{\left<E_{x}^{\rm meas}(x,y,V_{\rm gate})\right>_{V_{\rm
gate}}}{I_{\rm SD}/W}$ by averaging $E_{x}^{\rm meas}(x,y)$ across $V_{\rm
gate}=-141--91$ V for each position (Fig. S13a). Averaging instead across a
spatial area results in $\overline{\rho}_{\rm local}(V_{\rm
gate})=\frac{\left<E_{x}^{\rm meas}(x,y,V_{\rm gate})\right>_{x,y}}{I_{\rm
SD}/W}$ as shown in Fig. S13b for a few areas and compared to the measured two
probe resistivity of the device. While differences appear around charge
neutrality, all curves display a largely gate independent similar resistivity
for $n\leq-2\times 10^{17}\,{\rm m^{-2}}$.

Independence of resistivity from $V_{\rm gate}$, respectively from carrier
concentration $n$, is indicative of scattering dominated by short range
disorder, while $\rho\propto n^{-1}$ indicates dominant long range disorder
scattering 54. Therefore, we deduce that the former dominates at large hole
density. Electron-phonon scattering at longitudinal accoustic phonons results
as well in $n$-independent resistivity 55, but can be estimated as $\rho_{\rm
A}\approx 30\,{\rm\Omega}$ 73, 55, i.e., much lower than the measured
resistivity of $\rho=500\,\Omega$. Finally, scattering at remote interfacial
phonons (polar optical phonons of the substrate) can play a role for graphene
at room temperature 54, 73 exhibiting a resistivity contribution $\rho_{\rm
remote}\propto\frac{\hbar\omega_{\rm s}}{e^{\hbar\omega_{\rm s}/k_{\rm
B}T}-1}$, 74 where $\omega_{\rm s}$ is the frequency of the remote surface
polar optical phonon being $110\,{\rm meV}$ for ${\rm Si_{3}N_{4}}$ 75
($59\,{\rm meV}$ for ${\rm SiO_{2}}$) 73, 74. Since the remote phonon
contribution on ${\rm SiO_{2}}$ has been measured as $\rho_{\rm remote}\approx
26\,{\rm\Omega}$,73 we get for ${\rm Si_{3}N_{4}}$: $\rho_{\rm remote}\approx
26\,{\rm\Omega}\left(\frac{110\,{\rm meV}}{59\,{\rm
meV}}\right)\frac{e^{59\,{\rm meV}/k_{\rm B}T}-1}{e^{110\,{\rm meV}/k_{\rm
B}T}-1}=5\,{\Omega}\ll\rho_{\rm A}\ll\rho_{\rm loc}$. Thus, the remote phonon
contribution can be neglected as well. This implies a dominating scattering at
short range defects for large hole doping with scattering length $l_{\rm SR}$.

To calculate the gate dependence of the resulting $l_{\rm dis}(E_{\rm F})$, we
use the semiclassical Boltzmann transport equation for the scattering time
$\tau(E_{\rm F})=l_{\rm dis}(E_{\rm F})/v_{\rm F}$ reading54

$\frac{1}{\rho}=\frac{e^{2}v_{\rm F}^{2}}{2}D(E_{\rm F})\tau(E_{\rm
F})=\frac{e^{2}v_{\rm F}}{2}D(E_{\rm F})l_{\rm dis}(E_{\rm F})$ (S.24)

with density of states of graphene $D(E_{\rm F})=\frac{4\sqrt{\pi|n(E_{\rm
F})|}}{hv_{\rm F}}$. Substituting $D(E_{\rm F})$ into eq. (S.24) and solving
for $l_{\rm dis}$, we obtain eq. (5), main text, when using $\rho=\rho_{\rm
local}(x,y)$ in order to maintain the spatial character of $l_{\rm dis}(x,y)$
and $n(x,y)$. Figure S13c displays the resulting $l_{\rm dis}(n)$,
respectively $l_{\rm dis}(V_{\rm gate})$ using a capacitive model together
with its contributions $l_{\rm SR}$ and the electron-phonon contribution
$l_{\rm e-ph}\propto 1/\rho_{\rm A}$ as discussed above. Importantly, $l_{\rm
dis}$ increases towards charge neutrality ($V_{\rm D}$), oppositely to the
behaviour of $l_{\rm ee}$ (eq. (S.23)), such that dominating electron-electron
scattering naturally results at low doping.

###  S7.3 Comparing Local Scattering Lengths

Figure S14: Comparison of Different Scattering Lengths. (a–d)
$\widehat{E}_{x}^{\rm meas}(x,y,V_{\rm SD}=0.5\,$V) at the mentioned gate
voltages (partly shown in Fig. 5, main text), $V_{\rm D}=85\,$V. The
respective topography is shown in Fig.1f, main text. The white rectangle
highlights an area including inverted electric field regions. (e) Dependence
on $V_{\rm gate}$ (bottom axis) and carrier concentration (top axis) of
electron-electron scattering length $l_{\rm ee}$ (dark lines) and electron-
disorder scattering length $l_{\rm dis}$ (pink lines) for the marked area in b
(area 1, thick symbols) and the whole area of a–d (thin symbols). Dashed lines
labelled $V_{\rm gate}^{l_{\rm ee}<l_{\rm dis}}$ mark the range where $l_{\rm
ee}<l_{\rm dis}$ for the thick symbols. Inset: $V_{\rm gate}$ dependence of
$\overline{\rho}_{\rm loc}$ for area 1. (f) Map of $V_{\rm gate}^{l_{\rm
ee}<l_{\rm dis}}(x,y)$ for the same area as displayed in a–d using a spatial
average of $0.2\,\mu{\rm m}\times 0.2\,\mu{\rm m}$ for each point $(x,y)$. The
rectangle marks the same area as in b.

The ability to map, within the same area, the doping distribution $n_{0}(x,y)$
and the transport induced electric fields $E_{x}^{\rm meas}(x,y)$, leading to
$\rho_{\rm local}(x,y)$ at high hole doping, allows us to derive corresponding
maps of electron-electron scattering length $l_{\rm ee}(x,y)$ and electron-
disorder scattering length $l_{\rm dis}(x,y)$ via eq. (4) and eq. (5), main
text, respectively. Hence, we are able to locally compare these length scales
as shown in Fig. 4g, main text, revealing that inverted electric fields appear
preferentially in areas where the $V_{\rm gate}$ range with $l_{\rm ee}<l_{\rm
dis}$ is largest. In the particular case of Fig. 4, the inverted electric
fields appeared on the electron side close to charge neutrality. Here, we add
an example, where the inverted fields are observed on the hole side close to
charge neutrality (Fig. S14a–d), in particular, at the right of the fold
penetrating the whole image (topography: Fig. 1f, main text). This area also
shows lowest $E_{x}^{\rm meas}(x,y)$ far away from charge neutrality (Fig.
S14a) implying low $\rho_{\rm loc}(x,y)$ and, thus, large $l_{\rm dis}(x,y)$
(eq. (S.24)), hence, favoring electron viscosity around charge neutrality. In
line, again the areas with largest $V_{\rm gate}$ range, where $l_{\rm ee}$
dominates with respect to $l_{\rm dis}$ (Fig. S14e), showcase inverted
electric fields as more directly corroborated by mapping the $V_{\rm gate}$
range with $l_{\rm ee}<l_{\rm dis}$ directly, called $V_{\rm gate}^{l_{\rm
ee}<l_{\rm dis}}(x,y)$ (Fig. S14f) and comparing it with the $E_{x}^{\rm
meas}(x,y)$ maps of the identical region (Fig. S14a–d).

The inset in Fig. S14e displays $\overline{\rho}_{\rm loc}(V_{\rm gate})$ for
the highlighted area 1 exhibiting a wide $V_{\rm gate}$ range with constant
$\overline{\rho}_{\rm loc}$. This confirms our previous conclusion of
dominating short range disorder scattering for $l_{\rm dis}$, particularly for
the area with inverted electric fields.

The authors thank for helpful discussions with K. Sotthewes, C. B. Winkelmann,
O. Kolosov, V. Falko, and A. Tzalenchuk. This project has received funding
from the European Union’s Horizon 2020 research and innovation programme under
grant agreement number 881603 (Graphene Core3) as well as the Humboldt
foundation via a grant of S.S.

## References

  * Polini and Geim 2020 Polini, M.; Geim, A. K. Viscous electron fluids. _Phys. Today_ 2020, _73_ , 28–34
  * Mayzel et al. 2019 Mayzel, J.; Steinberg, V.; Varshney, A. Stokes flow analogous to viscous electron current in graphene. _Nat. Commun._ 2019, _10_ , 937
  * Schäfer and Teaney 2009 Schäfer, T.; Teaney, D. Nearly perfect fluidity: from cold atomic gases to hot quark gluon plasmas. _Rep. Prog. Phys._ 2009, _72_ , 126001
  * Narozhny et al. 2017 Narozhny, B. N.; Gornyi, I. V.; Mirlin, A. D.; Schmalian, J. Hydrodynamic Approach to Electronic Transport in Graphene. _Ann. Phys._ 2017, _529_ , 1700043
  * Narozhny 2019 Narozhny, B. N. Electronic hydrodynamics in graphene. _Ann. Phys._ 2019, _411_ , 167979
  * Torre et al. 2015 Torre, I.; Tomadin, A.; Geim, A. K.; Polini, M. Nonlocal transport and the hydrodynamic shear viscosity in graphene. _Phys. Rev. B_ 2015, _92_ , 165433
  * Gurzhi 1968 Gurzhi, R. N. Hydrodynamic effects in solids at low temperature. _Sov. Phys. Usp._ 1968, _11_ , 255–270
  * Govorov and Heremans 2004 Govorov, A. O.; Heremans, J. J. Hydrodynamic Effects in Interacting Fermi Electron Jets. _Physical Review Letters_ 2004, _92_ , 026803
  * Guo et al. 2017 Guo, H.; Ilseven, E.; Falkovich, G.; Levitov, L. S. Higher-than-ballistic conduction of viscous electron flows. _Proc. Nat. Acad. Sci._ 2017, _114_ , 3068–3073
  * Moessner et al. 2019 Moessner, R.; Morales-Durán, N.; Surówka, P.; Witkowski, P. Boundary-condition and geometry engineering in electronic hydrodynamics. _Phys. Rev. B_ 2019, _100_ , 155115
  * Holder et al. 2019 Holder, T.; Queiroz, R.; Scaffidi, T.; Silberstein, N.; Rozen, A.; Sulpizio, J. A.; Ella, L.; Ilani, S.; Stern, A. Ballistic and hydrodynamic magnetotransport in narrow channels. _Phys. Rev. B_ 2019, _100_ , 245305
  * Mohseni et al. 2005 Mohseni, K.; Shakouri, A.; Ram, R. J.; Abraham, M. C. Electron vortices in semiconductors devices. _Phys. Fluids_ 2005, _17_ , 100602
  * Levitov and Falkovich 2016 Levitov, L.; Falkovich, G. Electron viscosity, current vortices and negative nonlocal resistance in graphene. _Nat. Phys._ 2016, _12_ , 672–676
  * Danz and Narozhny 2020 Danz, S.; Narozhny, B. N. Vorticity of viscous electronic flow in graphene. _2D Materials_ 2020, _7_ , 035001
  * Chandra et al. 2019 Chandra, M.; Kataria, G.; Sahdev, D.; Sundararaman, R. Hydrodynamic and ballistic AC transport in two-dimensional Fermi liquids. _Phys. Rev. B_ 2019, _99_ , 165409
  * Lent 1990 Lent, C. S. Ballistic current vortex excitations in electron waveguide structures. _Appl. Phys. Lett._ 1990, _57_ , 1678–1680
  * Pellegrino et al. 2016 Pellegrino, F. M. D.; Torre, I.; Geim, A. K.; Polini, M. Electron hydrodynamics dilemma: Whirlpools or no whirlpools. _Phys. Rev. B_ 2016, _94_ , 155414
  * Mendoza et al. 2011 Mendoza, M.; Herrmann, H. J.; Succi, S. Preturbulent Regimes in Graphene Flow. _Phys. Rev. Lett._ 2011, _106_ , 156601
  * Li et al. 2020 Li, S.; Levchenko, A.; Andreev, A. V. Hydrodynamic electron transport near charge neutrality. _Phys. Rev. B_ 2020, _102_ , 075305
  * Hui et al. 2020 Hui, A.; Lederer, S.; Oganesyan, V.; Kim, E.-A. Quantum aspects of hydrodynamic transport from weak electron-impurity scattering. _Phys. Rev. B_ 2020, _101_ , 121107
  * de Jong and Molenkamp 1995 de Jong, M. J. M.; Molenkamp, L. W. Hydrodynamic electron flow in high-mobility wires. _Phys. Rev. B_ 1995, _51_ , 13389–13402
  * Bandurin et al. 2016 Bandurin, D. A.; Torre, I.; Kumar, R. K.; Shalom, M. B.; Tomadin, A.; Principi, A.; Auton, G. H.; Khestanova, E.; Novoselov, K. S.; Grigorieva, I. V.; Ponomarenko, L. A.; Geim, A. K.; Polini, M. Negative local resistance caused by viscous electron backflow in graphene. _Science_ 2016, _351_ , 1055–1058
  * Crossno et al. 2016 Crossno, J.; Shi, J. K.; Wang, K.; Liu, X.; Harzheim, A.; Lucas, A.; Sachdev, S.; Kim, P.; Taniguchi, T.; Watanabe, K.; Ohki, T. A.; Fong, K. C. Observation of the Dirac fluid and the breakdown of the Wiedemann-Franz law in graphene. _Science_ 2016, _351_ , 1058–1061
  * Ghahari et al. 2016 Ghahari, F.; Xie, H.-Y.; Taniguchi, T.; Watanabe, K.; Foster, M. S.; Kim, P. Enhanced Thermoelectric Power in Graphene: Violation of the Mott Relation by Inelastic Scattering. _Phys. Rev. Lett._ 2016, _116_ , 136802
  * Gallagher et al. 2019 Gallagher, P.; Yang, C.-S.; Lyu, T.; Tian, F.; Kou, R.; Zhang, H.; Watanabe, K.; Taniguchi, T.; Wang, F. Quantum-critical conductivity of the Dirac fluid in graphene. _Science_ 2019, _364_ , 158–162
  * Kumar et al. 2017 Kumar, R. K. et al. Superballistic flow of viscous electron fluid through graphene constrictions. _Nat. Phys._ 2017, _13_ , 1182–1185
  * Berdyugin et al. 2019 Berdyugin, A. I.; Xu, S. G.; Pellegrino, F. M. D.; Kumar, R. K.; Principi, A.; Torre, I.; Shalom, M. B.; Taniguchi, T.; Watanabe, K.; Grigorieva, I. V.; Polini, M.; Geim, A. K.; Bandurin, D. A. Measuring Hall viscosity of graphene’s electron fluid. _Science_ 2019, _364_ , 162–165
  * Bandurin et al. 2018 Bandurin, D. A.; Shytov, A. V.; Levitov, L. S.; Kumar, R. K.; Berdyugin, A. I.; Shalom, M. B.; Grigorieva, I. V.; Geim, A. K.; Falkovich, G. Fluidity onset in graphene. _Nat. Commun._ 2018, _9_ , 4533
  * Geurs et al. 2020 Geurs, J.; Kim, Y.; Watanabe, K.; Taniguchi, T.; Moon, P.; Smet, J. H. Rectification by hydrodynamic flow in an encapsulated graphene Tesla valve. _arXiv:_ 2020, _2008.04862_
  * Lucas and Fong 2018 Lucas, A.; Fong, K. C. Hydrodynamics of electrons in graphene. _J. Phys.: Cond. Matt._ 2018, _30_ , 053001
  * Moll et al. 2016 Moll, P. J. W.; Kushwaha, P.; Nandi, N.; Schmidt, B.; Mackenzie, A. P. Evidence for hydrodynamic electron flow in PdCoO2. _Science_ 2016, _351_ , 1061–1064
  * Fu et al. 2020 Fu, C.; Guin, S. N.; Scaffidi, T.; Sun, Y.; Saha, R.; Watzman, S. J.; Srivastava, A. K.; Li, G.; Schnelle, W.; Parkin, S. S. P.; Felser, C.; Gooth, J. Largely Suppressed Magneto-Thermal Conductivity and Enhanced Magneto-Thermoelectric Properties in PtSn4. _Research_ 2020, _2020_ , 4643507
  * Gooth et al. 2018 Gooth, J.; Menges, F.; Kumar, N.; Sü$\upbeta$, V.; Shekhar, C.; Sun, Y.; Drechsler, U.; Zierold, R.; Felser, C.; Gotsmann, B. Thermal and electrical signatures of a hydrodynamic electron fluid in tungsten diphosphide. _Nat. Commun._ 2018, _9_ , 4093
  * Block et al. 2020 Block, A.; Principi, A.; Hesp, N. C. H.; Cummings, A. W.; Liebel, M.; Watanabe, K.; Taniguchi, T.; Roche, S.; Koppens, F. H. L.; van Hulst, N. F.; Tielrooij, K.-J. Observation of giant and tunable thermal diffusivity of Dirac fluid at room temperature. _arXiv:_ 2020, _2008.04189_
  * Sulpizio et al. 2019 Sulpizio, J. A. et al. Visualizing Poiseuille flow of hydrodynamic electrons. _Nature_ 2019, _576_ , 75–79
  * Ku et al. 2020 Ku, M. J. H. et al. Imaging viscous flow of the Dirac fluid in graphene. _Nature_ 2020, _583_ , 537–541
  * Jenkins et al. 2020 Jenkins, A.; Baumann, S.; Zhou, H.; Meynell, S. A.; Yang, D.; Watanabe, K.; Taniguchi, T.; Lucas, A.; Young, A. F.; Jayich, A. C. B. Imaging the breakdown of ohmic transport in graphene. _arXiv:_ 2020, _2002.05065_
  * Braem et al. 2018 Braem, B. A.; Pellegrino, F. M. D.; Principi, A.; Röösli, M.; Gold, C.; Hennel, S.; Koski, J. V.; Berl, M.; Dietsche, W.; Wegscheider, W.; Polini, M.; Ihn, T.; Ensslin, K. Scanning gate microscopy in a viscous electron fluid. _Phys. Rev. B_ 2018, _98_ , 241304
  * Melitz et al. 2011 Melitz, W.; Shen, J.; Kummel, A. C.; Lee, S. Kelvin probe force microscopy and its application. _Surf. Sci. Rep._ 2011, _66_ , 1–27
  * Xu et al. 2018 Xu, J.; Chen, D.; Li, W.; Xu, J. Surface potential extraction from electrostatic and Kelvin-probe force microscopy images. _J. Appl. Phys._ 2018, _123_ , 184301
  * Panchal et al. 2013 Panchal, V.; Pearce, R.; Yakimova, R.; Tzalenchuk, A.; Kazakova, O. Standardization of surface potential measurements of graphene domains. _Sci. Rep._ 2013, _3_ , 2597
  * Burnett et al. 2011 Burnett, T.; Yakimova, R.; Kazakova, O. Mapping of Local Electrical Properties in Epitaxial Graphene Using Electrostatic Force Microscopy. _Nano Lett._ 2011, _11_ , 2324–2328
  * Datta et al. 2009 Datta, S. S.; Strachan, D. R.; Mele, E. J.; Johnson, A. T. C. Surface Potentials and Layer Charge Distributions in Few-Layer Graphene Films. _Nano Lett._ 2009, _9_ , 7–11
  * Altvater et al. 2019 Altvater, M. A.; Wu, S.; Zhang, Z.; Zhu, T.; Li, G.; Watanabe, K.; Taniguchi, T.; Andrei, E. Y. Electrostatic imaging of encapsulated graphene. _2D Mat._ 2019, _6_ , 045034
  * Falkovich and Levitov 2017 Falkovich, G.; Levitov, L. Linking Spatial Distributions of Potential and Current in Viscous Electronics. _Phys. Rev. Lett._ 2017, _119_ , 066601
  * Shaygan et al. 2017 Shaygan, M.; Otto, M.; Sagade, A. A.; Chavarin, C. A.; Bacher, G.; Mertin, W.; Neumaier, D. Low Resistive Edge Contacts to CVD-Grown Graphene Using a CMOS Compatible Metal. _Ann. Phys._ 2017, _529_ , 1–7
  * Yu et al. 2009 Yu, Y.-J.; Zhao, Y.; Ryu, S.; Brus, L. E.; Kim, K. S.; Kim, P. Tuning the graphene work function by electric field effect. _Nano Lett._ 2009, _9_ , 3430–3434
  * Samaddar et al. 2016 Samaddar, S.; Coraux, J.; Martin, S. C.; Grévin, B.; Courtois, H.; Winkelmann, C. B. Equal variations of the Fermi level and work function in graphene at the nanoscale. _Nanoscale_ 2016, _8_ , 15162–15166
  * Martin et al. 2007 Martin, J.; Akerman, N.; Ulbricht, G.; Lohmann, T.; Smet, J. H.; von Klitzing, K.; Yacoby, A. Observation of electron-hole puddles in graphene using a scanning single-electron transistor. _Nat. Phys._ 2007, _4_ , 144–148
  * Giuliani and Vignale 2005 Giuliani, G.; Vignale, G. _Quantum Theory of the Electron Liquid_ ; Cambridge University Press, 2005; pp 432–439
  * Polini and Vignale 2016 Polini, M.; Vignale, G. _No-nonsense Physicist: An Overview of Gabriele Giuliani’s Work and Life_ ; Edizioni della Normale, 2016; Vol. 2; pp 107–124
  * Kim et al. 2020 Kim, M. et al. Control of electron-electron interaction in graphene by proximity screening. _Nat. Commun._ 2020, _11_ , 2339
  * Sheehy and Schmalian 2007 Sheehy, D. E.; Schmalian, J. Quantum critical scaling in graphene. _Phys. Rev. Lett._ 2007, _99_ , 226803
  * Sarma et al. 2011 Sarma, S. D.; Adam, S.; Hwang, E. H.; Rossi, E. Electronic transport in two-dimensional graphene. _Rev. Mod. Phys._ 2011, _83_ , 407–470
  * Li and Sarma 2013 Li, Q.; Sarma, S. D. Finite temperature inelastic mean free path and quasiparticle lifetime in graphene. _Phys. Rev. B_ 2013, _87_ , 085406
  * Neumaier et al. 2019 Neumaier, D.; Pindl, S.; Lemme, M. C. Integrating graphene into semiconductor fabrication lines. _Nat. Mater._ 2019, _18_ , 525–529
  * Ella et al. 2019 Ella, L.; Rozen, A.; Birkbeck, J.; Ben-Shalom, M.; Perello, D.; Zultak, J.; Taniguchi, T.; Watanabe, K.; Geim, A. K.; Ilani, S.; Sulpizio, J. A. Simultaneous voltage and current density imaging of flowing electrons in two dimensions. _Nat. Nanotechnol._ 2019, _14_ , 480–487
  * Sinterhauf et al. 2020 Sinterhauf, A.; Traeger, G. A.; Pakdehi, D. M.; Schädlich, P.; Willke, P.; Speck, F.; Seyller, T.; Tegenkamp, C.; Pierz, K.; Schumacher, H. W.; Wenderoth, M. Substrate induced nanoscale resistance variation in epitaxial graphene. _Nat. Commun._ 2020, _11_ , 555
  * Li et al. 2009 Li, X.; Zhu, Y.; Cai, W.; Borysiak, M.; Han, B.; Chen, D.; Piner, R. D.; Colombo, L.; Ruoff, R. S. Transfer of large-area graphene films for high-performance transparent conductive electrodes. _Nano Lett._ 2009, _9_ , 4359–4363
  * Adam et al. 2007 Adam, S.; Hwang, E. H.; Galitski, V. M.; Sarma, S. D. A self-consistent theory for graphene transport. _Proc. Nat. Acad. Sci._ 2007, _104_ , 18392––18397
  * 61 https://www.bruker.com/products/surface-and-dimensional-analysis/atomic-force-microscopes/dimension-icon/overview.html
  * Nečas and Klapetek 2012 Nečas, D.; Klapetek, P. Gwyddion: an open-source software for SPM data analysis. _Central European Journal of Physics_ 2012, _10_ , 181–188
  * Zerweck et al. 2005 Zerweck, U.; Loppacher, C.; Otto, T.; Grafström, S.; Eng, L. M. Accuracy and resolution limits of Kelvin probe force microscopy. _Phys. Rev. B_ 2005, _71_ , 125424
  * 64 https://www.brukerafmprobes.com/p-3950-scm-pit-v2.aspx
  * S.Hudlet et al. 1998 S.Hudlet,; Jean, M. S.; Guthmann, C.; Berger, J. Evaluation of the capacitive force between an atomic force microscopy tip and a metallic surface. _Eur. Phys. J. B_ 1998, _2_ , 5–10
  * Strassburg et al. 2005 Strassburg, E.; Boag, A.; Rosenwaks, Y. Reconstruction of electrostatic force microscopy images. _Rev. Sci. Instr._ 2005, _76_ , 083705
  * Li et al. 2012 Li, G.; Mao, B.; Lan, F.; Liu, L. Practical aspects of single-pass scan Kelvin probe force microscopy. _Rev. Sci. Instrum_ 2012, _83_ , 113701
  * Willke et al. 2016 Willke, P.; Möhle, C.; Sinterhauf, A.; Kotzott, T.; Yu, H. K.; Wodtke, A.; Wenderoth, M. Local transport measurements in graphene on SiO2 using Kelvin probe force microscopy. _Carbon_ 2016, _102_ , 470–476