††thanks: Equal contribution††thanks: Equal contribution

# Electronic Poiseuille Flow in Hexagonal Boron Nitride Encapsulated Graphene
FETs

Wenhao Huang  Empa, Swiss Federal Laboratories for Materials Science and
Technology, Transport at Nanoscale Interfaces Laboratory, Überlandstrasse 129,
CH-8600 Dübendorf, Switzerland.  Department of Physics, University of Basel,
Klingelbergstrasse 82, CH-4056 Basel, Switzerland.  Tathagata Paul
<EMAIL_ADDRESS>Empa, Swiss Federal Laboratories for Materials Science
and Technology, Transport at Nanoscale Interfaces Laboratory, Überlandstrasse
129, CH-8600 Dübendorf, Switzerland.  Kenji Watanabe  Research Center for
Functional Materials, National Institute for Materials Science, 1-1 Namiki,
Tsukuba, Ibaraki 305-0044, Japan.  Takashi Taniguchi  International Center for
Materials Nanoarchitectonics, National Institute for Materials Science, 1-1
Namiki, Tsukuba, Ibaraki 305-0044, Japan.  Mickael L. Perrin
<EMAIL_ADDRESS>Empa, Swiss Federal Laboratories for Materials Science
and Technology, Transport at Nanoscale Interfaces Laboratory, Überlandstrasse
129, CH-8600 Dübendorf, Switzerland.  Department of Information Technology and
Electrical Engineering, ETH Zürich, 8092 Zürich, Switzerland.  Michel Calame
<EMAIL_ADDRESS>Empa, Swiss Federal Laboratories for Materials Science
and Technology, Transport at Nanoscale Interfaces Laboratory, Überlandstrasse
129, CH-8600 Dübendorf, Switzerland.  Department of Physics, University of
Basel, Klingelbergstrasse 82, CH-4056 Basel, Switzerland.  Swiss Nanoscience
Institute, University of Basel, Klingelbergstrasse 82, CH-4056 Basel,
Switzerland.

# Electronic Poiseuille Flow in Hexagonal Boron Nitride Encapsulated Graphene
FETs

Wenhao Huang  Empa, Swiss Federal Laboratories for Materials Science and
Technology, Transport at Nanoscale Interfaces Laboratory, Überlandstrasse 129,
CH-8600 Dübendorf, Switzerland.  Department of Physics, University of Basel,
Klingelbergstrasse 82, CH-4056 Basel, Switzerland.  Tathagata Paul
<EMAIL_ADDRESS>Empa, Swiss Federal Laboratories for Materials Science
and Technology, Transport at Nanoscale Interfaces Laboratory, Überlandstrasse
129, CH-8600 Dübendorf, Switzerland.  Kenji Watanabe  Research Center for
Functional Materials, National Institute for Materials Science, 1-1 Namiki,
Tsukuba, Ibaraki 305-0044, Japan.  Takashi Taniguchi  International Center for
Materials Nanoarchitectonics, National Institute for Materials Science, 1-1
Namiki, Tsukuba, Ibaraki 305-0044, Japan.  Mickael L. Perrin
<EMAIL_ADDRESS>Empa, Swiss Federal Laboratories for Materials Science
and Technology, Transport at Nanoscale Interfaces Laboratory, Überlandstrasse
129, CH-8600 Dübendorf, Switzerland.  Department of Information Technology and
Electrical Engineering, ETH Zürich, 8092 Zürich, Switzerland.  Michel Calame
<EMAIL_ADDRESS>Empa, Swiss Federal Laboratories for Materials Science
and Technology, Transport at Nanoscale Interfaces Laboratory, Überlandstrasse
129, CH-8600 Dübendorf, Switzerland.  Department of Physics, University of
Basel, Klingelbergstrasse 82, CH-4056 Basel, Switzerland.  Swiss Nanoscience
Institute, University of Basel, Klingelbergstrasse 82, CH-4056 Basel,
Switzerland.

###### Abstract

Electron-electron interactions in graphene are sufficiently strong to induce a
correlated and momentum-conserving flow such that charge carriers behave
similarly to the Hagen-Poiseuille flow of a classical fluid. In the current
work, we investigate the electronic signatures of such a viscous charge flow
in high-mobility graphene FETs. In two complementary measurement schemes, we
monitor differential resistance of graphene for different channel widths and
for different effective electron temperatures. By combining both approaches,
the presence of viscous effects is verified in a temperature range starting
from 178 K and extending up to room temperature. Our experimental findings are
supported by finite element calculations of the graphene channel, which also
provide design guidelines for device geometries that exhibit increased viscous
effects. The presence of viscous effects near room temperature opens up
avenues for functional hydrodynamic devices such as geometric rectifiers like
a Tesla valve and charge amplifiers based on electronic Venturi effect.

##  I Introduction

The dynamics governing the flow of charges in conductors is heavily dependent
on the ratios of the different physical length scales present in the system.
These length scales dictate how the charge carriers interact with each other,
with the system boundaries, and with defects or impurities. Microscopically,
the flow of charge carriers is resisted by scattering from defects and lattice
vibrations (phonons), characterized by the length scale, $l_{diff}$, which is
the average distance traveled by the charge carriers between two such
momentum-relaxing (inelastic) scattering events. The charge carriers also
interact with each other in a momentum-conserving (elastic) fashion, with a
length scale, $l_{ee}$. The third relevant length scale is the channel width,
$w$ (assuming channel length $l\gg w$), dictating how often charge carriers
interact with the boundary of the system.

In most conductors, diffusive scattering prevails ($l_{diff}<l_{ee},~{}w$) and
charge transport is Ohmic. Charge carriers traveling across the channel suffer
many momentum relaxing collisions leading to a constant drift velocity along
the direction of the applied electric field [1]. Alternatively, transport is
ballistic, when the channel dimensions are the smallest length scale in the
system. Under these circumstances, charge carriers travel collision-less from
the source to the drain terminal, dissipating energy only at the contacts [2].
However, a third, and relatively unexplored transport regime emerges, when
inter-particle scattering dominates ($l_{ee}<l_{diff},w$). These strong inter-
particle interactions lead to a correlated, viscous charge flow, similar to
the flow of fluids and governed by the theory of hydrodynamics [3, 4, 5, 6, 7,
8, 9, 10]. A direct consequence of this behavior is the non-uniformity of the
carrier velocity in the channel, which follows a Poiseuille-like flow profile,
similar to fluids (Fig. 1a). In the context of electrical transport, this
translates to a (non-uniform) position dependent current density, being
largest in the center of the channel and decreasing on approaching the channel
edges, following a parabolic profile.

Fulfilling the criteria of charge hydrodynamics requires high mobility devices
with reduced scattering from defects and phonons (large $l_{diff}$) and strong
inter-particle interactions (small $l_{ee}$, generally observed in materials
with a simple Fermi surface). Graphene fulfills both these requirements.
First, the electron-electron interaction strength is enhanced by the nearly
circular Fermi surface at low number densities (small $l_{ee}$) [11, 12].
Second, graphene has a weak electron-phonon coupling strength [13, 14, 15] and
advancements in fabrication techniques make it possible to limit scattering
from extrinsic disorders, in particular by encapsulation in hexagonal boron
nitride ($h$-BN). $h$-BN also serves as an ideal substrate for graphene, due
to its atomically flat surface and low lattice mismatch ($\sim 1.7\%$) [16]
(Fig. 1a). These lead to graphene heterostructures having outstanding
electronic properties and large $l_{diff}$.

Figure 1: Device structure and electrical performance. a. Schematic depiction
of the graphene heterostructure used in the current work. The Poiseuille flow
profile is depicted above the Top $h-$BN. The arrows and the curve represent
the particle velocity and the velocity profile expected for Poiseuille flow,
respectively. b. Transfer characteristics of a typical $GrFET$ device at
different values of $T$. c. Diffusive scattering length scale ($l_{diff}$) of
$GrFET$ device extracted from the transfer characteristics in subsection b.

In recent years, hydrodynamic transport of charge has attracted significant
theoretical attention [5, 17, 18, 19, 20], and several experimental signatures
have been reported [21, 22, 23, 6, 7, 8]. Of particular interest is the direct
observation of the parabolic flow profile. Recent investigations utilizing
direct imaging of Hall field [24] and local magnetic field [22, 25] have
hinted at the presence of viscous effects in charge transport in graphene and
layered Weyl semimetals. Additionally, unique transport properties like
formation of vortices, negative local resistance and superballisticity have
also been linked to a hydrodynamic charge flow [26, 6, 23]. Apart from the
Coulomb mediated electron-electron interactions [5, 17, 18, 19, 20], recent
reports also indicate the possibility of phonon mediated electronic
interactions generating a correlated charge flow at higher carrier densities
[26, 25]. Though electrical signatures of Poiseuille flow in graphene have
been reported [6, 23], a comparative study of the transport signatures from
multiple types of experiments and a determination of the temperature range for
the viscous effects is lacking.

In the current work, we perform electrical transport measurements in $h$-BN
encapsulated graphene field effect transistors ($GrFETs$) to investigate the
presence of electronic Poiseuille flow. For this purpose two separate
measurement schemes are used, 1. Current Bias measurements (Section III) and
2. Width dependent measurements (Section IV). For the current bias
measurements, we look into the scaling of channel resistance as a function of
DC bias current in a sample of constant width (schematic in Fig. 2a). The
width dependent measurements investigate the flow of charge carriers in
channels of different widths (schematic in Fig. 3a). We observe signatures of
correlated charge flow for temperatures starting from $\sim 178$ K and
persisting upto 300 K. Further, finite element calculations are performed to
verify the experimental observations (Section V). Our results demonstrate the
presence of strong electronic correlations near room temperature opening up
possibilities of utilizing the viscous electronic flow in graphene for
functional devices.

##  II Electrical performance of $GrFETs$

The quality of the fabricated $GrFETs$ was assessed by measuring the
conductance of the channel as a function of gate voltage, as called a transfer
characteristic at different temperatures (Fig. 1b). The devices were
fabricated using a standard dry transfer technique [27] and measurements
performed in a four-probe geometry to remove the contact resistance (See
Methods). $h$-BN encapsulation allows fabrication of edge-contacts, which are
known to outperform conventional surface contacts [28]. Fig. 1b demonstrates a
typical graphene transfer characteristic, with a maximum of resistance near
the Dirac point and a sharp decrease with increasing number density ($n$).
This corresponds to the transition from insulating behavior, caused by the
formation of charge puddles near the neutrality point, to metallic behavior
[29, 30]. The high quality of our FETs is demonstrated by a large carrier
mobility ($\mu\geq 2\times 10^{5}$) (Supplementary Fig. S1) and a large
$\l_{diff}\geq$ 1 $\mu$m (Fig. 1c), thereby relaxing the criteria for
hydrodynamic behavior ($l_{ee}<l_{diff},w$) for properly designed channels.
The value of $l_{diff}$ is computed considering a Drude model,

$l_{diff}=\frac{\sigma m^{*}v_{f}}{ne^{2}}$ (1)

where, $\sigma$ is the channel conductivity, $m^{*}$ the cyclotron mass [31]
and $v_{f}=10^{6}$ ms-1, the Fermi velocity of graphene. We observe a
reduction in the $l_{diff}$ with increasing temperature (especially at $T=$337
K) which can be attributed to increased momentum relaxing scattering events
[16]. Temperature dependence of $l_{d}iff$ is provided in Supplementary Fig.
S1.

Figure 2: Current bias measurements in $GrFETs$. a. Schematic depiction of a
typical $GrFET$ device with circuit connections for current bias measurements.
$I_{AC}$ and $I_{DC}$ represent the AC and DC bias current. $V_{g}$ is the
gate bias. b. Optical image of the device used for current bias measurements.
Red and yellow dot lines highlight the bottom and top $h$-BN flakes,
respectively. The scale bar is 5 $\mu$m. c. Color plots of d$V$/d$I$ as a
function of $I_{DC}$ and $n$ at $T=$ 9 K (left) and 150 K (right). d. Plots
depicting the $I_{DC}$ dependence of differential resistance at $T=$ 9 K for
three different $n$ values indicated by white dashed lines in subsection c..
The number density ranges for the different d$V$/d$I$ scaling behaviors are
indicated at the top of the color plot (left panel subsection c.) by solid
lines of the same color as the corresponding scatter plot in d. e. Evolution
of d$V$/d$I$ curves at $n=$ 1.6 $\times$ 1011cm-2 (white dashed line for $T=$
9 K and 150 K data in subsection c.) for increasing values $T$ from 9 K (top
panel) to 300 K (bottom panel). The Knudsen, Poiseuille and Diffusive
transport regimes are indicated by color bars in the top panel.

##  III Current Bias Measurements

With the high quality of our $GrFETs$ established, we turn to current bias
measurements which consist of measuring the device conductance for increasing
DC bias current. The role of the DC bias ($|I_{DC}|\leq 300$ $\mu$A) is to
increase the electronic temperature of the devices, thereby providing a knob
to increase the inter-particles interactions and reduce $l_{ee}$. At the same
time, the weak electron-phonon coupling in SLG ensures the lattice temperature
does not significantly increase, which would lead to an increase in diffusive
scattering processes, and hence a decrease in $l_{diff}$ [13, 14, 15]. The
decoupling of the two length scales therefore creates an ideal scenario for
hydrodynamic charge flow [9].

Fig. 2a presents the device geometry and circuit connections used for the
current bias measurements and Fig. 2b is an optical micrograph of the measured
device. The measured device has a channel length and width of 1.1 $\mu$m and
0.55 $\mu$m, respectively. The transfer characteristics, mobility and
$l_{diff}$ of the device are presented in Supplementary Fig. S2. Fig. 2c shows
color plots of d$V$/d$I$ as a function of $n$ and $I_{DC}$ recorded at 9 K and
150 K. At $T=9$ K, for changing number densities, three types of $I_{DC}$
curves are observed. The number densities in which these behaviors are
observed are indicated above the color plot. For each behavior, a typical
$I_{DC}$ curve is plotted in Fig. 2d, with the number densities highlighted by
the vertical dashed lines in Fig. 2c . For low $n$ (Scaling behavior i@), the
d$V$/d$I$ values follow a monotonic decreasing trend with increasing $I_{DC}$
until about $5\times 10^{-5}~{}$A ($J_{min}\approx 0.09$ mA/$\mu$m). For
larger values of $I_{DC}$, the curve is stable with only a minute increase in
resistance. This behavior results in a central bright spot in the color map.
With increasing $n$, we observe a non-monotonicity in the measured d$V$/d$I$,
as apparent from the central bright spot that splits into two ridges and a
valley appearing in the center (Scaling behavior ii@).

For Scaling behavior ii@, the $I_{DC}$ curves exhibit three distinct regions.
First, for low DC currents, the d$V$/d$I$ values increase with current.
d$V$/d$I$ then reaches a peak and in the second region, decreases for
increasing current (reaching a minimum of d$V$/d$I$), while in the third
region, d$V$/d$I$ steadily increases with DC current. The experimentally
observed non-monotonicity in the differential resistance hints towards the
presence of viscous effects ($n=1.6\times 10^{11}$ cm-2 in Fig. 2d), as we
will discuss later on. Finally, for the largest $n$ (Scaling behavior iii@),
the negative slope region vanishes and the d$V$/d$I$ is constant for small
$I_{DC}$ with a monotonic increase for higher current values.

The three different Scaling behaviors with changing $n$ can be rationalized as
follows. For Scaling behavior i@, electrical transport is dominated by the
formation of charge puddles, leading to an insulating behavior at low DC
currents, with a transition to metallic behavior for increasing DC current
[29, 30, 32, 33]. The initial drop in resistance is associated with electric
field assisted interband carrier generation [32]. We can calucate the charge
inhomogeneity ($\Delta n$) near the Dirac point from the minima in the
d$V$/d$I$. This gives $\Delta n=J_{min}/ev_{f}\approx 5\times 10^{10}$ cm-2,
which is similar to those reported in encapsulated graphene [33]. The metallic
behavior at higher $I_{DC}$ has been attributed to a Schwinger production of
electron hole plasma in graphene [33]. For Scaling behavior ii@, increasing
the DC current leads to three distinct regions in the $I_{DC}$ curve in which
the charge transport is quasi-ballistic/Knudsen, hydrodynamic/Poiseuille and
diffusive, respectively. For low DC currents, the charge carriers flow nearly
collisionless, leading to a quasi-ballistic transport, also known as the
electronic Knudsen regime. The transfer characteristics also indicate the
presence of quasi-ballistic transport (Supplementary Fig. S2). Upon increasing
$I_{DC}$, the electron temperature increases and consequently the inter-
particle collisions are increased. Even though such collisions are momentum
conserving and lead to no resistance, they make it more probable for the
carriers to reach and diffusively interact with the sample boundaries [9],
leading to an increase in the resistance (d(d$V$/d$I$)/d$I_{DC}>0$). A further
increase in $I_{DC}$ until $l_{ee}<w$ brings about a transition to the
Poiseuille regime. Here, the inter-particle collisions actively prevent a
large fraction of the charge carriers (ones near the center of the channel)
from seeing the boundaries, resulting in a slope reversal,
d(d$V$/d$I$)/d$I_{DC}<0$, and a minimum in the differential resistance. This
negative slope for intermediate $I_{DC}$ and minima of d$V$/d$I$ in clean
conductors are signatures of the electronic Gurzhi effect [34]. This has been
experimentally observed previously and is indicative of the presence of
viscous/Poiseuille flow [6, 9]. This regime lasts until diffusive scattering
processes within the sample gain prominence, restoring the positive slope of
d$V$/d$I$ as expected in standard metallic behavior or the diffusive regime
(d(d$V$/d$I$)/d$I_{DC}>0$). A similar behavior is observed in four different
samples (See Methods for details). Finally, for higher number densities
(Scaling behavior iii@), the inter-particle interactions are screened (larger
$l_{ee}$) [35, 36] while the diffusive scattering remains largely unaffected
(constant $l_{diff}$ in Supplementary Fig. S2 and Fig. 1c) [6]. A combination
of these effects leads to the vanishing of the Poiseuille flow (negative slope
of d$V$/d$I$). The d$V$/d$I$ values are nearly constant for low $I_{DC}$
(quasi-ballistic regime) and show an increasing trend once diffusive
scattering starts dominating the transport around DC currents of $3\times
10^{-5}~{}$A (0.05 mA/$\mu$m). Similar observations were also made for hole
transport in $GrFETs$ (negative $n$), and are presented in Supplementary Fig.
S3.

We also investigate the effect of the lattice temperature ($T$) on the Gurzhi
effect by recording the conductance maps at different temperatures. The color
plots at all temperatures are presented in Supplementary Fig. S4. First, we
compare the conductance map at 150 K (Fig. 2c (right panel)) to the one
obtained at 9 K (Fig. 2c (left panel)). Two major differences are observed.
First, the central bright spot has a reduced intensity and, second, the
Knudsen flow is absent for higher $T$ values. Color plots at intermediate
temperatures are presented in the Supplementary Fig. S4. The reduction in
intensity of the central bright spot is indicative of a reduced Dirac point
resistance at higher values of $T$ in $GrFETs$. This is an effect of the
enhanced particle-hole excitation, similar to our observations in Fig. 1b. The
disappearance of the Knudsen regime is better visualized when the $I_{DC}$
curves for an intermediate number density of $n=~{}$1.6$\times$1011 cm-2 are
investigated for all temperatures (Fig. 2e). At low $T$, all three regimes,
Knudsen, Poiseuille and diffusive are visible. However, with increasing $T$,
the Knudsen regime first reduces in size, until it disappears at $T=$ 150 K.
Thereafter, at $T=$ 300 K, the Poiseuille regime also vanishes and the
$GrFETs$ only demonstrate signatures of diffusive transport.

This $T$-dependence can be understood by considering the variation of $l_{ee}$
and $l_{diff}$ with lattice temperature. At $T=9$ K, $l_{ee},~{}l_{diff}\geq
w$, which results in a quasi-ballistic or Knudsen flow persisting upto large
values of $I_{DC}$ (Supplementary Fig. S2). A transition from Knudsen to
Poiseuille flow occurs when $l_{ee}<w$. With increasing $T$, the Poiseuille
flow condition is satisfied at smaller values of $I_{DC}$, leading to a
diminishing Knudsen flow regime. This trend persists until the lattice
temperature is high enough for the sample to enter the Poiseuille flow regime
at $I_{DC}=$ 0\. This is indicated by the absence of the Knudsen regime and is
observed at $T=$ 150 K (Supplementary Fig. S2). The sample remains in the
Poiseuille regime until diffusive scattering processes gain prominence
($l_{diff}<l_{ee}$) and the differential resistance shows a monotonic positive
slope (300 K data in Fig. 2e). Our differential resistance measurements, thus,
hint towards the presence of a viscous Poiseuille flow in $GrFETs$ at
temperatures between 150 and 300 K.

Figure 3: Width dependent measurements in $GrFETs$ a. Left Panel, schematic
overview of $GrFET$ design and circuit connection used for width dependent
measurements. Regions of different widths ($w$) are fabricated on the same
heterostructure using EBL and RIE (See Methods). Top right, false colored
scanning electron microscopy (SEM) image, used for measuring the channel
dimensions. Scale bar is 1 $\mu$m. The SEM image is from a region depicted by
a red dotted line in the left panel. Bottom right, schematic of the device
cross-section depicting the correction introduced in the calculation of $w$
due to the etching angle. b. Scaling of $\sigma$ as a function of $w$ at
$T=200$ K. Solid line is a fit to the measured values (closed circles)
following Eq. 5. c. d. and e. represent $l_{diff}$, $D_{\nu}$ and $l_{v}$
values, respectively, as a function of $T$, extracted from the fits in
subsection b. Red line (both solid and dotted) and shaded region in e.
represent theoretically predicted scaling of $l_{\nu}$, following Eq. 7 and
the temperature range for Poiseuille efffects in $GrFETs$ from width dependent
measurements, respectively. The red colored bar above subsection e. represents
the temperature range for Poiseuille effects observed from the current bias
measurements (Fig. 2 and Section III). Error bars in subsection c., d. and e.
are calculated from the fitting error in b.

##  IV Width Dependent Measurements

In addition to differential resistance measurements, we also verify the
presence of Poisuille flow using geometric effects in the channel conductivity
($\sigma$). For diffusive transport, $\sigma$ is a material parameter and
therefore the current density (J) is uniform across the device cross-section.
However, in the Poiseuille regime, the parabolic profile of J (or carrier
velocity) (Fig. 1a) manifests itself as a strong width-dependence in $\sigma$.
For a 2D Fermi Liquid in the Poiseuille regime, the transport equations are
given by [20]

$\mathbf{\nabla}\cdot\mathbf{J(r)}=0$ (2)
$\frac{\sigma_{0}}{e}\mathbf{\nabla}\phi(r)+D_{\nu}^{2}\nabla^{2}\mathbf{J(r)}=\mathbf{J(r)}$
(3)

where, Eq. 2 and Eq. 3 are the electrical analogues of the continuity and
Navier-Stokes equations, respectively. $\sigma_{0}$ is the Drude conductivity,
$\phi(r)$ the electric potential and $e$ the electronic charge. $D_{\nu}$, the
Gurzhi length, is a charcteristic length scale determining the strength of
viscous effects. Its functional dependence is given by

$D_{\nu}=\sqrt{\frac{\nu\l_{diff}}{v_{f}}}$ (4)

where, $\nu$ is the kinematic viscosity. The second term of Eq. 3 is the
viscous term and is responsible for the geometry dependence of $\sigma$. An
analytical relation for $\sigma$ is obtained by solving Eq. 2 and Eq. 3 with a
diffusive/no-slip boundary condition. The resulting functional dependence is

$\sigma=\sigma_{0}\left[1-2\frac{D_{\nu}}{w}\mathrm{tanh}\left(\frac{w}{2D_{\nu}}\right)\right]$
(5)

From Eq. 5, we can define an asymptotic limit of strong Poiseuille flow using
the condition $D_{\nu}>>w$. In this case, we obtain,
$\sigma=\sigma_{0}w^{2}/12D_{\nu}$. This $w^{2}$ scaling of $\sigma$ has been
reported previously in the Fermi liquid of Weyl semimetals [8]. However,
comparing the values of $w$ used in the current work ($0.6<w<5~{}\mu$m) and
previously reported values of $D_{\nu}\approx 0.2-0.4~{}\mu$m in graphene [6],
we expect that the measured $GrFETs$ do not satisfy the conditions for strong
Poiseuille flow. Hence, unlike previous reports [8], we find it suitable to
fit the experimental $\sigma$ with the non-asymptotic relation given by Eq. 5.

Fig. 3a shows a schematic of the $GrFETs$ used for our width-dependent
measurements. To fabricate these devices, the $h$-BN encapsulated graphene
heterostructure was shaped into channels of different widths using reactive
ion etching (RIE) (See Methods). Channel widths ($w_{1}$) and lengths ($l$)
were confirmed by scanning electron micrographs (top right panel of Fig. 3a).
The actual channel width ($w$) was estimated from the observed width
($w_{1}$), plus taking into consideration the height of the top $h$-BN and the
45 degree angle generated by the RIE process (bottom right panel of Fig. 3a)
[28]. The width-dependent fitting allows us to extract several flow
parameters, providing further insight into the viscous charge transport in
$GrFETs$. A detailed discussion of the process is presented below.

The transfer characteristic and $l_{diff}$ of the width dependent device is
presented in Fig. 1b and c, respectively. Fig. 3b shows the fitting of
$\sigma$ for lattice temperature $T=200~{}$K and $n=7\times 10^{11}$ cm-2.
This particular choice of $T$ and $n$ is guided by previous reports of viscous
effects in graphene at similar values of these parameters [6]. The closed
circles are measured values of $\sigma$, while the solid line is a fit
following Eq. 5, with $\sigma_{0}$ and $D_{\nu}$ as fitting parameters. From
the fit, one can directly extract the values of $\sigma_{0}$ and $D_{\nu}$.
The extracted value of $\sigma_{0}$ enables us to compute the diffusive mean
free path $l_{diff}$ using Eq. 1. We find $l_{diff}\sim 1~{}\mu$m, which is a
common occurrence in high quality graphene devices (Fig. 3c) [6, 28]. The
value of $D_{\nu}\sim 0.2~{}\mu$m (Fig. 3d), obtained from the fit is also in
close agreement with previous reports [6]. Substituting $l_{diff}$ and
$D_{\nu}$ in Eq. 4 gives us an estimate of the viscosity ($\nu$) of the charge
carriers. We obtain a high kinematic viscosity, $\nu\approx 0.04~{}$m2s-1
($T=200~{}$K), making the flow more viscous than honey, a phenomenon which has
been reported previously (Supplementary Fig. S5) [6]. This high viscosity is
also one of the assumptions for using the linearized Navier-Stokes equation in
Eq. 3. Additionally, the flow viscosity, which arises due to the inter-
particle scattering, can be directly related to the viscous mean free path,
$l_{\nu}$, a parameter which is closely related to $l_{ee}$ [6, 19]. The exact
relationship at low excitation frequencies ($f<<10^{12}$ Hz) is given by

$l_{\nu}=\frac{4m^{*}\nu}{\hbar k_{f}}$ (6)

The extracted $l_{\nu}\approx 0.14~{}\mu$m at $T=200$ K, is the smallest
length scale in the system ($l_{\nu}<<w,~{}l_{diff}$), indicative of the
presence of Poiseuille flow in $GrFETs$ at $T=200$ K (Fig. 3e). To explore
further, we repeat the width-dependent measurements at different lattice
temperatures and explore the temperature dependence of $l_{diff}$, $D_{\nu}$,
$l_{\nu}$ and $\nu$ in Fig. 3c, d, e and Supplementary Fig. S5, respectively.
The flow parameters demonstrate a $T-$dependence only in an intermediate
temperature range of $178<T<300$ K and saturate (except for $l_{diff}$) for
$T<178$ K and $T>300$ K.

To understand the physical relevance of these values, we compare the extracted
$l_{\nu}$ with its predicted $T-$dependence for the Dirac fermionic liquid in
a doped graphene sheet [19]. The analytical relation for $l_{\nu}$ is given by

$l_{\nu}=\frac{45v_{f}(1+N_{f}\alpha_{ee})^{2}\hbar\epsilon_{f}}{64\pi\sqrt{2}N_{f}\alpha_{ee}^{2}(k_{B}T)^{2}}$
(7)

where, $N_{f}=4$ is the number of Fermionic flavours in graphene,
$\alpha_{ee}=7.3$ the fine structure constant and $k_{B}$ the Boltzmann
constant. $\alpha_{ee}$ is computed from the relation [19, 37]

$\alpha_{ee}=\frac{e^{2}}{\epsilon\hbar v_{f}}$ (8)

where $\epsilon$ is the electrical permittivity. For the current work, we have
used the permittivity of $h-$BN considering a dielectric constant of 4 [38].
Eq. 7 (solid line in Fig. 3e) has no fitting parameters and matches the
extracted values of $l_{\nu}$(closed circles in Fig. 3e) for the intermediate
temperature range, $178<T<$ 300 K (shaded region in Fig. 3e). The extracted
flow parameters are thus physically relevant only for these intermediate
temperatures. Additionally, we note that the assumptions of viscous flow made
in the formulation of the transport equations (Eq. 2 & 3) limit their
validity, and as a result, the validity of the geometry dependent analysis for
values of $T$ where Poiseuille flow is present. In combination, the physical
relevance of the extracted flow parameters and the validity of the viscous
flow equations confirm the presence of Poiseuille flow for $178<T<300$ K in
$GrFETs$. Finally, this observation is strengthened by the fact that similar
$T$ values are obtained for both the current bias measurements (Section III)
and width dependent measurements (Fig. 2e and Fig. 3e). We also look at the
number density dependence of the flow parameters in Suppementary Fig. S6. The
geometry dependent measurements indicate the presence of hydrodynamic effects
for number densities $>10^{11}$ cm-2. Further discussions are provided in
Suppementary Fig. S6.

Figure 4: Finite element modelling of graphene channel. a. Plots of $J_{norm}$
as a function of $w$ and $Position~{}(norm.)$ at $T=178~{}$K depicting a
parabolic current flow profile for small values of $w$ and a uniform
$J_{norm}$ for the largest simulated widths. b. Flow curvature (solid line) as
a function of $w$ computed from flow profiles in a. at $Position~{}(norm.)=0$.
Filled circles indicate the curvatures for measured values of $w$ in the
current work. The dashed line represents curvature for an ideal parabolic
flow. Shaded region indicates $w$ values for strong Poiseuille flow (from
inset of subsection d.). c.,d.,e., and f. Experimental (filled circles) and
simulated (solid line) values of $\sigma$ as a function of $w$ for
$T=40,178,250$ and $337$ K, respectively. Dashed line in subsection d. and e.
indicates the $w^{2}$ scaling of $\sigma$ expected for strong Poiseuille
effects. Inset of d. compares the $w^{2}$ scaling and simulated $\sigma$ for
$w\leq 0.5~{}\mu$m. Shaded region represents the predicted $w$ values for
strong Poiseuille flow. The finite element modelling is performed at
$n=7\times 10^{11}$ cm-2.

##  V Finite element methods

Finally, we turn to finite element calculations (implemented in Comsol
Multiphysics 5.6) to model the Poiseuille flow in $GrFETs$. This process
involves numerically solving Eq. 2 and 3 with an uniform electric field
($\nabla\phi$) and no-slip boundary conditions for different values of $w$.
The input parameters, $D_{\nu}$ and $\sigma_{0}$, are calculated from the
theoretical predictions of $l_{\nu}$ (Eq. 7) and the field effect mobility
($\mu$). Our channels are oriented along the x-axis while their width follows
the y-axis. An example of a simulated flow profile is depicted in Fig. 4a for
$T=178~{}$K. We observe strong Poiseuille effects for small channel widths,
indicated by the parabolic dependence of the normalized current density
($J_{norm}=J_{x}(y)/J_{x}(y=0)$) on the normalized channel position
($Position~{}(norm.)=y/(w/2)$), and a weak Poiseuille flow, indicated by a
constant $J_{norm}$, for larger values of $w$. The curvature of the flow
profile computed at the center of the channel ($Position~{}(norm.)=0$) also
exhibits a width dependence (solid line in Fig. 4b). It saturates to a value
of 2 for the smallest widths with a parabolic flow profile ($w=0.1~{}\mu$m in
Fig. 4a) and reduces to 0 for larger widths with a constant $J_{norm}$
($w=5~{}\mu$m in Fig. 4a). The filled circles in Fig. 4b correspond to the
measured values of $w$ in the current work. Further, we extract the channel
conductivity ($\sigma$) from the simulated flow profile. In Fig. 4c to f, the
simulated (solid red line) and experimental (filled circles) values of
$\sigma$ are compared for four different values of $T$. The simulated results
match experimental data for temperatures in the Poiseuille flow regime
($T=178$ and $250$ K), with mismatches at lower ($T=40$ K) and higher ($T=337$
K) temperatures, where our experiments indicate the absence of hydrodynamic
effects.

In what follows, we quantify our observations and predict the values of $w$
necessary for a strong Poiseuille flow. As discussed before, the strong
Poiseuille regime is characterized by a $w^{2}$ scaling of $\sigma$. This is
indicated by the dashed blue lines in Fig. 4d and e, computed with the same
parameter values as in the simulation (solid red line). Indeed, the
experimentally measured widths in our $GrFETs$ devices are too large for
observing this effect as the flow is in the weak Poiseuille regime. The strong
Poiseuille regime is determined by computing the relative difference between
the simulated (solid red line) and strong Poiseuille (dashed blue line) values
of $\sigma$. Values of $w$ for which the relative difference remains below 10%
are considered to be in the strong Poiseuille regime. The criteria is
satisfied for $w<0.15~{}\mu$m, indicated by the shaded region in the inset of
Fig. 4d. Another indicator of the strong Poiseuille regime is the saturation
of the flow curvature to a value of 2 (shaded region in Fig. 4b). In
conclusion, our finite element model confirms our $GrFETs$ are in the weak
Poiseuille regime and predicts strong Poiseuille flow for $w<0.15~{}\mu$m in
$GrFETs$. Comparing the measured values of $l_{\nu}$ (Fig. 3e) with the
computed length scale for strong hydrodynamics, we would expect this behavior
to be present at a slightly higher temperature than the range indicated in the
current work in order to satisfy the criteria $l_{ee}<w$.

##  VI Discussions and Conclusions

Strikingly, we observe hydrodynamic effects at different charge number
densities in the geometry dependence study ($n>1\times 10^{11}$ cm-2)
(Supplementary Fig. S6) and the current bias measurements ($1.6\times
10^{11}\lesssim n\lesssim 2.4\times 10^{11}$ cm-2), for which we have no
clear, unique explanation. Similar behavior has been seen in previous reports,
with current bias measurements showing viscous effects for lower number
density ranges compared to negative vicinity resistance [6]. This could
indicate either a number density selectivity in the transport probes of charge
hydrodynamics or a number density variability of the hydrodynamic regime in
graphene.

In conclusion, we demonstrate the electrical signatures of Poiseuille flow in
charge transport through graphene channels using the Gurzhi effect and
geometry dependent measurements. We observe a negative slope in the
differential resistance measurements and a strong geometry dependence of the
channel conductivity, both indicative of the presence of viscous effects. The
two different measurements types, in conjunction, enable us to define a
temperature range for the viscous effects to occur, which reaches as high as
300 K. Moreover, we provide an electrical transport-based characterization
framework for detecting Poiseuille flow in conductors that can easily be
extended to other conductors. Our experimental observations are corroborated
by finite-element calculations that shed light on the channel geometries
required for a strong Poiseuille flow. Our findings offer promising prospects
for functional hydrodynamic graphene devices in the near future, such as
geometric rectifiers like Tesla valves and charge amplifiers based on the
electronic Venturi effects [39, 40, 41].

##  VII Methods

###  VII.1 Device Fabrication

$GrFETs$ measured in the current work are fabricated using the following
process: the heterostructure assembly begins with the mechanical exfoliation
of individual layers from bulk crystals (the graphite is obtained commercially
from NGS Trading & Consulting GmbH and the $h$-BN is from our co-authors in
National Institute for Materials Science, Japan). Exfoliated flakes are
optically screened using an optical microscope (Zeiss Axio Imager M2m). The
layer number and quality of grphene flakes are determined using Raman
spectroscopy (See Supplementary Fig. S7). We use a Raman system with a
backscattering geometry confocal Raman microscope (WITec, Alpha 300 R). An
excitation laser with a wavelength of 532 nm from a diode laser is used for
all Raman measurements. The thickness of $h-$BN layers is characterized using
atomic force microscopy in the tapping mode (Bruker Icon3 AFM). Selected
flakes are aligned and stacked layer by layer using commercially available
micromanipulator (hq Graphene 2D Heterostructure Transfer System). The pick-up
process is performed using PDMS backed polycarbonate (PC) stamps. PC shows
good adherence to $h$-BN for $330~{}$K$\leq T\leq 350~{}$K. After successful
pick-up of the top $h$-BN layer, the middle graphene layer and bottom $h$-BN
layer are picked up using the van der Waals interaction between the different
two-dimensional materials. Finally, the heterostructure is deposited on a pre-
patterned Si++/ SiO2 (285 nm) wafer at $T=450~{}$K. The PC is dissolved in
dchlorometane (DCM). The transfer process followed is similar to that in Ref.
[27]. Following this, a thermal annealing step is performed at 650 K with
H2/Ar 35/200 sccm for 3 hours to improve the quality of the heterostructures
by reducing interfacial bubbles, contaminants etc. (See Supplementary Fig. S8
and S9). The small electrodes are defined using EBL and RIE (CHF${}_{3}=40$
sccm, O${}_{2}=4$ sccm, $P=60$ mTorr, Power = 60 W). Next, we deposit 5/25 nm
Cr/Au directly using ebeam evaporator for edge-contacted electrodes and the
pattern is lifted off using acetone for 45 min [28]. After that, the second
EBL and metals deposition (5/65 nm Cr/Au) are performed for contact pads.
Finally, we perform the third ebeam processing, followed by the second RIE for
defining the channel shape. The morphology of the device is recorded by the
same optical microscope mentioned above. Channel widths of devices are
acquired from SEM (Hitachi S-4800). For the width dependent device the channel
width ($w$) was varied while keeping the channel length constant.

For the current bias measurements, four devices on four different chips were
fabricated. All samples demonstrate a similar scaling of d$V$/d$I$ at base
temperature. We observe the current biasing measurements have a tendency to
generate irreversible hysteresis and electron doping (shifting of the Dirac
point to negative gate voltages) in devices. Due to this we could perform the
temperature characterization in one device and partial temperature
characterization in a second device.

For the width dependent measurements, the largest four values of $w$ are from
a single large graphene channel which was shaped into channels of different
widths, while the smallest $w$ value is obtained from a second device on a
different chip made following similar fabrication protocols.

###  VII.2 Measurement Details

The four-probe resistance of $GrFETs$ was measured using a lock-in technique,
with an AC source-drain current ( of $I_{AC}=100$ nA (231.45 Hz); $I_{DC}=0$)
while measuring the resulting voltage drop (Lockins used are SRS 830 or EG&G
7265). The SRS CS580 voltage controlled current source is used to supply
$I_{AC}$ to the sample. The gate voltage is supplied by an AdWin Gold ii@
connected to a voltage amplifier (Physics Basel SP 908). From the applied gate
voltage, the number density is computed using $n=C_{ox}(V_{g}-V_{d})/e$ [42],
where, $C_{ox}$ is the oxide capacitance (series combination of SiO2 (285 nm)
and bottom $h$-BN), $e$ is the electron charge, $V_{g}$ and $V_{d}$ are the
gate and the Dirac point voltage, respectively.

For the current bias measurements, the DC bias is supplied by an AdWin Gold
ii@ and the AC bias by a lockin. The two bias voltages are added using a
summing amplifier SIM980. The current source CS580 converts the summed voltage
to the supplied current bias. An $I_{AC}=100$ nA (231.5 Hz) is used for the
measurements and the output voltage from the sample measured using the lockin.

The width dependent measurements follow similar measurement technique as the
transfer characteristic measurements.

The current bias measurements are carried out on a closed cycle manual probe
station (Lakeshore). For width dependence measurements, the samples are wire
bonded to a chip carrier and measured in a commercially available closed cycle
cryostat from Advanced Research Systems (ARS).

##  VIII Acknowledgement

W.H., T.P. and M.C. acknowledge funding from the Swiss National Science
Foundation under the Sinergia grant no. 189924 (Hydronics). K.W. and T.T.
acknowledge support from the JSPS KAKENHI (Grant Numbers 19H05790 and
20H00354). M.L.P. acknowledges funding from the Swiss National Science
Foundation under Spark grant no. 196795. and the Eccellenza Professorial
Fellowship no. PCEFP2_203663, as well as supported by the Swiss State
Secretariat for Education, Research and Innovation (SERI) under contract
number MB22.00076. The authors acknowledge support from the Multiphysics Hub @
Empa for the COMSOL Multiphyics calculations. We thank the Cleanroom
Operations Team of the Binnig and Rohrer Nanotechnology Center (BRNC) for
their help and support. We also thank our project partners Dr. Bernd Gotsmann,
Prof. Ilaria Zardo, Prof. Nicola Marzari, Dr. Ivan Shorubalko for discussions.
W.H. would like to thank Prof. Ilaria Zardo for her mentorship.

##  IX Author contributions statement

W.H. and T.P. contributed equally to this work. W.H., T.P., M.L.P. and M.C.
conceived and designed the experiments. W.H. and T.P. prepared the devices and
performed the electrical measurements. W.H. and T.P. performed the annealing,
Raman, SEM and AFM measurements. K.W. and T. T. synthesized the hBN flakes.
W.H., T.P., M.L.P., and M.C. did the data analysis. M.L.P. performed the
finite element calculations, with the help of T.P.. W.H., T.P., M.L.P., and
M.C. discussed the figures and wrote the manuscript. M.L.P. and M.C.
supervised the study. All authors discussed the results and implications and
commented on the manuscript.

##  X Competing interests

The authors declare that there are no competing interests.

## References

  * Ashcroft and Mermin [1976] Ashcroft, N. W.; Mermin, N. D. _Solid state physics_ ; Holt, Rinehart and Winston, New York London, 1976.
  * Datta [1997] Datta, S. _Electronic transport in mesoscopic systems_ ; Cambridge University Press, 1997.
  * Landau and Lifshitz [2013] Landau, L. D.; Lifshitz, E. M. _Fluid Mechanics_ ; Elsevier, 2013; Vol. 6.
  * Batchelor and Batchelor [2000] Batchelor, C. K.; Batchelor, G. _An introduction to fluid dynamics_ ; Cambridge university press, 2000.
  * Ho et al. [2018] Ho, D. Y.; Yudhistira, I.; Chakraborty, N.; Adam, S. Theoretical determination of hydrodynamic window in monolayer and bilayer graphene from scattering rates. _Phys. Rev. B_ 2018, _97_ , 121404.
  * Bandurin et al. [2016] Bandurin, D.; Torre, I.; Kumar, R. K.; Shalom, M. B.; Tomadin, A.; Principi, A.; Auton, G.; Khestanova, E.; Novoselov, K.; Grigorieva, I.; Ponomarenko, L. A.; Geim, A. K.; Polini, M. Negative local resistance caused by viscous electron backflow in graphene. _Science_ 2016, _351_ , 1055–1058.
  * Moll et al. [2016] Moll, P. J.; Kushwaha, P.; Nandi, N.; Schmidt, B.; Mackenzie, A. P. Evidence for hydrodynamic electron flow in PdCoO2. _Science_ 2016, _351_ , 1061–1064.
  * Gooth et al. [2018] Gooth, J.; Menges, F.; Kumar, N.; Sü$\beta$, 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_ , 1–8.
  * De Jong and Molenkamp [1995] De Jong, M.; Molenkamp, L. Hydrodynamic electron flow in high-mobility wires. _Phys. Rev. B_ 1995, _51_ , 13389.
  * Keser et al. [2021] Keser, A. C.; Wang, D. Q.; Klochan, O.; Ho, D. Y. H.; Tkachenko, O. A.; Tkachenko, V. A.; Culcer, D.; Adam, S.; Farrer, I.; Ritchie, D. A.; Sushkov, O. P.; Hamilton, A. R. Geometric Control of Universal Hydrodynamic Flow in a Two-Dimensional Electron Fluid. _Phys. Rev. X_ 2021, _11_ , 031030.
  * Oka et al. [2019] Oka, T.; Tajima, S.; Ebisuoka, R.; Hirahara, T.; Watanabe, K.; Taniguchi, T.; Yagi, R. Ballistic transport experiment detects Fermi surface anisotropy of graphene. _Phys. Rev. B_ 2019, _99_ , 035440.
  * Roldán et al. [2008] Roldán, R.; López-Sancho, M. P.; Guinea, F. Effect of electron-electron interaction on the Fermi surface topology of doped graphene. _Phys. Rev. B_ 2008, _77_ , 115410.
  * Johannsen et al. [2013] Johannsen, J. C.; Ulstrup, S.; Bianchi, M.; Hatch, R.; Guan, D.; Mazzola, F.; Hornekær, L.; Fromm, F.; Raidel, C.; Seyller, T., et al. Electron–phonon coupling in quasi-free-standing graphene. _J. Phys. Condens. Matter_ 2013, _25_ , 094001.
  * Si et al. [2013] Si, C.; Liu, Z.; Duan, W.; Liu, F. First-principles calculations on the effect of doping and biaxial tensile strain on electron-phonon coupling in graphene. _Phys. Rev. Lett._ 2013, _111_ , 196802.
  * Calandra and Mauri [2007] Calandra, M.; Mauri, F. Electron-phonon coupling and electron self-energy in electron-doped graphene: Calculation of angular-resolved photoemission spectra. _Phys. Rev. B_ 2007, _76_ , 205411.
  * Dean et al. [2010] Dean, C. R.; Young, A. F.; Meric, I.; Lee, C.; Wang, L.; Sorgenfrei, S.; Watanabe, K.; Taniguchi, T.; Kim, P.; Shepard, K. L.; Hone, J. Boron nitride substrates for high-quality graphene electronics. _Nat. Nanotechnol._ 2010, _5_ , 722.
  * Pellegrino et al. [2016] Pellegrino, F. M.; Torre, I.; Geim, A. K.; Polini, M. Electron hydrodynamics dilemma: Whirlpools or no whirlpools. _Phys. Rev. B_ 2016, _94_ , 155414.
  * Svintsov [2018] Svintsov, D. Hydrodynamic-to-ballistic crossover in Dirac materials. _Phys. Rev. B_ 2018, _97_ , 121405.
  * Principi et al. [2016] Principi, A.; Vignale, G.; Carrega, M.; Polini, M. Bulk and shear viscosities of the two-dimensional electron liquid in a doped graphene sheet. _Phys. Rev. B_ 2016, _93_ , 125410.
  * 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.
  * Bandurin et al. [2018] Bandurin, D. A.; Shytov, A. V.; Levitov, L. S.; Kumar, R. K.; Berdyugin, A. I.; Ben Shalom, M.; Grigorieva, I. V.; Geim, A. K.; Falkovich, G. Fluidity onset in graphene. _Nat. Commun._ 2018, _9_ , 1–8.
  * Ku et al. [2020] Ku, M. J.; Zhou, T. X.; Li, Q.; Shin, Y. J.; Shi, J. K.; Burch, C.; Anderson, L. E.; Pierce, A. T.; Xie, Y.; Hamo, A., et al. Imaging viscous flow of the Dirac fluid in graphene. _Nature_ 2020, _583_ , 537–541.
  * Kumar et al. [2017] Kumar, R. K.; Bandurin, D.; Pellegrino, F.; Cao, Y.; Principi, A.; Guo, H.; Auton, G.; Shalom, M. B.; Ponomarenko, L. A.; Falkovich, G., et al. Superballistic flow of viscous electron fluid through graphene constrictions. _Nat. Phys._ 2017, _13_ , 1182–1185.
  * Sulpizio et al. [2019] Sulpizio, J. A.; Ella, L.; Rozen, A.; Birkbeck, J.; Perello, D. J.; Dutta, D.; Ben-Shalom, M.; Taniguchi, T.; Watanabe, K.; Holder, T., et al. Visualizing Poiseuille flow of hydrodynamic electrons. _Nature_ 2019, _576_ , 75–79.
  * Vool et al. [2021] Vool, U.; Hamo, A.; Varnavides, G.; Wang, Y.; Zhou, T. X.; Kumar, N.; Dovzhenko, Y.; Qiu, Z.; Garcia, C. A.; Pierce, A. T., et al. Imaging phonon-mediated hydrodynamic flow in WTe2. _Nat. Phys._ 2021, _17_ , 1216–1220.
  * Aharon-Steinberg et al. [2022] Aharon-Steinberg, A.; Völkl, T.; Kaplan, A.; Pariari, A. K.; Roy, I.; Holder, T.; Wolf, Y.; Meltzer, A. Y.; Myasoedov, Y.; Huber, M. E.; Yan, B.; Falkovich, G.; Levitov, L. S.; Hücker, M.; Zeldov, E. Direct observation of vortices in an electron fluid. _Nature_ 2022, _607_ , 74–80.
  * Zomer et al. [2014] Zomer, P.; Guimarães, M.; Brant, J.; Tombros, N.; Van Wees, B. Fast pick up technique for high quality heterostructures of bilayer graphene and hexagonal boron nitride. _Appl. Phys. Lett._ 2014, _105_ , 013101\.
  * Wang et al. [2013] Wang, L.; Meric, I.; Huang, P.; Gao, Q.; Gao, Y.; Tran, H.; Taniguchi, T.; Watanabe, K.; Campos, L.; Muller, D.; Guo, J.; Kim, P.; Hone, J.; Shepard, K. L.; Dean, C. R. One-dimensional electrical contact to a two-dimensional material. _Science_ 2013, _342_ , 614–617.
  * Martin et al. [2008] Martin, J.; Akerman, N.; Ulbricht, G.; Lohmann, T.; Smet, J. v.; Von Klitzing, K.; Yacoby, A. Observation of electron–hole puddles in graphene using a scanning single-electron transistor. _Nat. Phys._ 2008, _4_ , 144–148.
  * Xue et al. [2011] Xue, J.; Sanchez-Yamagishi, J.; Bulmash, D.; Jacquod, P.; Deshpande, A.; Watanabe, K.; Taniguchi, T.; Jarillo-Herrero, P.; LeRoy, B. J. Scanning tunnelling microscopy and spectroscopy of ultra-flat graphene on hexagonal boron nitride. _Nat. Mater._ 2011, _10_ , 282–285.
  * Neto et al. [2009] Neto, A. C.; Guinea, F.; Peres, N. M.; Novoselov, K. S.; Geim, A. K. The electronic properties of graphene. _Rev. Mod. Phys._ 2009, _81_ , 109.
  * Vandecasteele et al. [2010] Vandecasteele, N.; Barreiro, A.; Lazzeri, M.; Bachtold, A.; Mauri, F. Current-voltage characteristics of graphene devices: Interplay between Zener-Klein tunneling and defects. _Phys. Rev. B_ 2010, _82_ , 045416.
  * Berdyugin et al. [2022] Berdyugin, A. I.; Xin, N.; Gao, H.; Slizovskiy, S.; Dong, Z.; Bhattacharjee, S.; Kumaravadivel, P.; Xu, S.; Ponomarenko, L.; Holwill, M., et al. Out-of-equilibrium criticalities in graphene superlattices. _Science_ 2022, _375_ , 430–433.
  * Gurzhi [1963] Gurzhi, R. Minimum of resistance in impurity-free conductors. _J Exp Theor Phys_ 1963, _17_ , 521.
  * Giuliani and Vignale [2005] Giuliani, G.; Vignale, G. _Quantum theory of the electron liquid_ ; Cambridge University Press, 2005.
  * Razeghi [2006] Razeghi, M. _Fundamentals of solid state engineering_ ; Springer, 2006\.
  * Kotov et al. [2012] Kotov, V. N.; Uchoa, B.; Pereira, V. M.; Guinea, F.; Neto, A. C. Electron-electron interactions in graphene: Current status and perspectives. _Rev. Mod. Phys._ 2012, _84_ , 1067.
  * Kim et al. [2012] Kim, K. K.; Hsu, A.; Jia, X.; Kim, S. M.; Shi, Y.; Dresselhaus, M.; Palacios, T.; Kong, J. Synthesis and characterization of hexagonal boron nitride film as a dielectric layer for graphene devices. _ACS nano_ 2012, _6_ , 8583–8590.
  * Taubert et al. [2011] Taubert, D.; Schinner, G. J.; Tomaras, C.; Tranitz, H.-P.; Wegscheider, W.; Ludwig, S. An electron jet pump: The Venturi effect of a Fermi liquid. _J. Appl. Phys._ 2011, _109_ , 102412.
  * Taubert et al. [2010] Taubert, D.; Schinner, G.; Tranitz, H.; Wegscheider, W.; Tomaras, C.; Kehrein, S.; Ludwig, S. Electron-avalanche amplifier based on the electronic Venturi effect. _Phys. Rev. B_ 2010, _82_ , 161416.
  * Szelong et al. [2017] Szelong, M.; Ludwig, A.; Wieck, A. D.; Kunze, U. Positive centre voltage in T-branch junctions on n-type GaAs/AlGaAs based on hydrodynamics. _Semicond. Sci. Technol._ 2017, _32_ , 105005.
  * Novoselov et al. [2004] Novoselov, K. S.; Geim, A. K.; Morozov, S. V.; Jiang, D.; Zhang, Y.; Dubonos, S. V.; Grigorieva, I. V.; Firsov, A. A. Electric Field Effect in Atomically Thin Carbon Films. _Science_ 2004, _306_ , 666–669.


# PINNtomo: Seismic tomography using physics-informed neural networks

Umair bin Waheed
Department of Geosciences
King Fahd University of Petroleum and Minerals
Dhahran 31261, Saudi Arabia.
<EMAIL_ADDRESS>
&Tariq Alkhalifah
Physical Sciences and Engineering Division
King Abdullah University of Science and Technology
Thuwal 23955, Saudi Arabia. &Ehsan Haghighat
Department of Civil Engineering
Massachusetts Institute of Technology
MA 02139, USA. &Chao Song
Physical Sciences and Engineering Division
King Abdullah University of Science and Technology
Thuwal 23955, Saudi Arabia. &Jean Virieux
ISTERRE
Université Grenoble Alpes
Saint-Martin-d’Heres 38400, France

###### Abstract

Seismic traveltime tomography using transmission data is widely used to image
the Earth’s interior from global to local scales. In seismic imaging, it is
used to obtain velocity models for subsequent depth-migration or full-waveform
inversion. In addition, cross-hole tomography has been successfully applied
for a variety of applications, including mineral exploration, reservoir
monitoring, and CO2 injection and sequestration. Conventional tomography
techniques suffer from a number of limitations, including the use of a
smoothing regularizer that is agnostic to the physics of wave propagation.
Here, we propose a novel tomography method to address these challenges using
developments in the field of scientific machine learning. Using seismic
traveltimes observed at seismic stations covering part of the computational
model, we train neural networks to approximate the traveltime factor and the
velocity fields, subject to the physics-informed regularizer formed by the
factored eikonal equation. This allows us to better compensate for the ill-
posedness of the tomography problem compared to conventional methods and
results in a number of other attractive features, including computational
efficiency. We show the efficacy of the proposed method and its capabilities
through synthetic tests for surface seismic and cross-hole geometries.
Contrary to conventional techniques, we find the performance of the proposed
method to be agnostic to the choice of the initial velocity model.

_K_ eywords Tomography $\cdot$ Inversion $\cdot$ Traveltimes $\cdot$ Neural
networks $\cdot$ Machine learning

##  1 Introduction

Seismic tomography has been used over the years as a pre-eminent tool for
subsurface model building at various scales ranging from global and regional
scales in earthquake seismology [1] to local scales in exploration seismology
[2]. Building a velocity macro-model is a crucial step for the success of
depth migration [3] and full-waveform inversion [4] for high-resolution
imaging of the Earth’s crust. First arrival traveltime tomography based on
refraction data or diving waves has been successfully used to build such
initial models [5]. In particular, the method is attractive for land seismic
data processing with surface acquisition, since it is often difficult or even
impossible to identify reflections in the data. Moreover, cross-hole seismic
tomography, which uses direct arrival times, has been around for more than
three decades in oil/gas [6] and mineral exploration [7]. It has also been
successful for applications in void and tunnel detection [8], reservoir
characterization [9], fracture detection [10], hydrological parameter
estimation [11], geotechnical site investigations [12], and time-lapse studies
related to carbon capture and sequestration [13].

Seismic tomography is typically solved as an inverse problem that minimizes
the misfit between a set of observed arrival times on the receivers and those
synthetically generated using an estimate of the velocity model. Minimization
of this misfit function requires a nonlinear optimization procedure. However,
the conventional ray tomography approach linearizes the tomography operator,
which requires the computation of the Fréchet derivatives. Then the linearized
tomography operator is inverted iteratively. For modern seismic surveys, the
requirement of explicitly computing the Fréchet derivatives is challenging to
handle in terms of the computation cost and memory requirements. This gave way
to the adjoint-state method [14, 15] that formulates tomographic inversion as
a nonlinear optimization process by directly computing the gradient of the
misfit function.

Nevertheless, these conventional methods still suffer from a number of
limitations. Usually, they use some form of smoothing regularization to
compensate for the ill-posedness of the problem. This ends up limiting the
resolution of the inverted velocity model. In addition, these methods
typically need an initial model with some general background features of the
Earth represented, like a constant depth gradient. The choice of the initial
model may affect the final solution and is usually not obvious prior to
inversion. Furthermore, for models with irregular topography, considerable
grid and algorithmic adaptions are needed to account for the free-surface
topography [16].

Therefore, in this work, we propose a novel algorithm for the seismic
tomography problem based on developments in the field of scientific machine
learning. In particular, we use the emerging paradigm of physics-informed
neural networks (PINNs) that overcomes the limitation of deep learning
associated with sparse data by incorporating the governing partial
differential equation (PDE) into the neural network’s loss function. PINNs
have already demonstrated success in solving a number of forward and inverse
problems in other scientific disciplines [17, 18]. Recently, PINNs have also
shown remarkable success in overcoming limitations associated with
conventional techniques in modeling seismic traveltimes [19, 20] and
wavefields [21, 22].

Here, we develop a PINN-based tomography (PINNtomo) algorithm to invert for
the velocity model. Given traveltimes at seismic stations covering part of the
computational domain, we use neural networks to approximate the traveltime
factor and the velocity fields, subject to the physics-informed regularizer
based on the factored eikonal equation. Doing so allows us to better
compensate for the poorly determined aspects of the velocity model compared to
conventional physics-agnostic smoothing regularizers. Also, we find the
performance of the method to be independent of the initial velocity model.
Moreover, since the method is mesh-free, it is easily adaptable to models with
irregular topography without modifications. Additional advantages of the
method include ease of deployment across a variety of platforms (CPUs, GPUs)
and architectures (desktops, clusters) without any modification.

Through tests on realistic surface seismic and cross-hole geometries, we
demonstrate the efficacy of the proposed algorithm in solving the tomography
problem. This is done by obtaining a velocity model that produces traveltimes
matching those observed at seismic stations while honoring the physics of wave
propagation by minimizing the residual of the eikonal equation at selected
grid points in the computational domain.

##  2 Theory

In an isotropic medium, the eikonal equation relates the gradient of the
traveltime surfaces to the velocity of the wavefront through the relation:

$\displaystyle|\nabla T(\mathbf{x})|^{2}$
$\displaystyle=\frac{1}{v^{2}(\mathbf{x})},\qquad\forall\,\mathbf{x}\,\in\,\Omega,$
(1) $\displaystyle T(\mathbf{x}_{s})$ $\displaystyle=0,$

where $\Omega$ is a domain in $\mathbb{R}^{d}$ with $d$ as the space
dimension, $T(\mathbf{x})$ is the traveltime or Euclidean distance to any
point $\mathbf{x}$ from the point-source $\mathbf{x}_{s}$, $v(\mathbf{x})$ is
the velocity defined on $\Omega$, and $\nabla$ denotes the gradient operator.

Since equation (1) contains singularity at the point-source location,
traveltime modeling studies [19, 20] have shown that a factored form of the
eikonal equation is easier to train using PINNs. Therefore, we factorize the
traveltime $T(\mathbf{x})$ into two multiplicative functions [23], i.e.,

$T(\mathbf{x})=T_{0}(\mathbf{x})\,\tau(\mathbf{x}),$ (2)

where $T_{0}(\mathbf{x})$ is the known function which is computed
analytically, leaving $\tau(\mathbf{x})$ as the unknown traveltime factor. By
substituting the above in equation (1), we get the factored eikonal equation:

$\displaystyle T_{0}^{2}\,|\nabla\tau|^{2}+\tau^{2}\,|\nabla T_{0}|^{2}$
$\displaystyle+2\,T_{0}\,\tau\,(\nabla
T_{0}.\nabla\tau)=\frac{1}{v^{2}(\mathbf{x})},$ (3)
$\displaystyle\tau(\mathbf{x}_{s})$ $\displaystyle=1.$

The known traveltime $T_{0}$ is computed analytically using the expression:

$T_{0}(\mathbf{x})=\frac{|\mathbf{x}-\mathbf{x}_{s}|}{v(\mathbf{x}_{s})},$ (4)

where $v(\mathbf{x}_{s})$ is the velocity at the source location. This ensures
that $T_{0}$ captures the point-source singularity leaving $\tau$ as a smooth
function in the source neighborhood.

We can re-write the factored eikonal equation in its residual form as:

$L(\mathbf{x}):T_{0}^{2}\,|\nabla\tau|^{2}+\tau^{2}\,|\nabla
T_{0}|^{2}+2\,T_{0}\,\tau\,(\nabla
T_{0}.\nabla\tau)-\frac{1}{v^{2}(\mathbf{x})}=0.$ (5)

To invert for the unknown velocity model, we consider two multilayer neural
networks – one to approximate the unknown traveltime factor for an arbitrary
source location $\mathbf{x}_{s}$, $\tau(\mathbf{x}_{s},\mathbf{x})$, and the
other for the velocity, $v(\mathbf{x})$, i.e.,

$\displaystyle\tau(\mathbf{x}_{s},\mathbf{x})\approx\hat{\tau}(\mathbf{x}_{s},\mathbf{x})$
$\displaystyle=\mathcal{N}_{\tau}(\mathbf{x}_{s},\mathbf{x};\bm{\theta}_{\tau}),$
(6) $\displaystyle v(\mathbf{x})\approx\hat{v}(\mathbf{x})$
$\displaystyle=\mathcal{N}_{v}(\mathbf{x};\bm{\theta}_{v}),$

where $\mathcal{N}_{\tau}$ and $\mathcal{N}_{v}$ are the neural networks with
trainable parameters $\bm{\theta}_{\tau}$ and $\bm{\theta}_{v}$, respectively.
Since traveltime data from multiple sources are needed to obtain a reliable
velocity model, the traveltime factor network also takes shot locations as
input, in addition to the spatial coordinates. On the contrary, the velocity
model only takes spatial coordinates as input since it does not vary with the
source locations.

Since both the traveltime factor and velocity are strictly positive
quantities, we pass the output of the networks through a sigmoid function,
$\sigma()$, and multiply these by scaling coefficients, i.e.,

$\displaystyle\hat{\tau}(\mathbf{x}_{s},\mathbf{x})$
$\displaystyle=\sigma\left(\mathcal{N}_{\tau}(\mathbf{x}_{s},\mathbf{x};\bm{\theta}_{\tau})\right)\,\tau_{\mathrm{peak}},$
(7) $\displaystyle\hat{v}(\mathbf{x})$
$\displaystyle=\sigma\left(\mathcal{N}_{v}(\mathbf{x};\bm{\theta}_{v})\right)\,v_{\mathrm{peak}},$

where $\tau_{\mathrm{peak}}$ and $v_{\mathrm{peak}}$ are the peak values that
can be obtained from the traveltime factor network and the velocity network,
respectively. These scaling factors should be chosen such that they are larger
than the expected maximum values to avoid clipping of the output.

Finally, we use a single loss function to train both networks simultaneously.
The loss function is given as:

$\displaystyle\mathfrak{J}(\bm{\theta}_{\tau}$
$\displaystyle,\bm{\theta}_{v})=\frac{1}{N_{s}\,N_{r}}\sum_{n=1}^{N_{s}}\sum_{i=1}^{N_{r}}\left(T_{0}(\mathbf{x}_{n,i})\,\tau(\mathbf{x}_{n,i})-\tilde{T}_{n,i}\right)^{2}$
(8) $\displaystyle+$
$\displaystyle\frac{1}{N_{s}}\sum_{n=1}^{N_{s}}\left(\hat{\tau}(\mathbf{x}_{n,s})-1\right)^{2}+\frac{1}{N_{s}\,N_{t}}\sum_{n=1}^{N_{s}}\sum_{i=1}^{N_{t}}\left(L(\mathbf{x}_{n,i})\right)^{2},$

where $N_{s}$ denotes the total number of sources, $N_{r}$ is the total number
of receivers, and $N_{t}$ is the number of training (collocation) points from
the computational domain. The first term minimizes the misfit between
traveltimes predicted by the neural network and the observed traveltimes,
$\tilde{T}_{n,i}$, over all sources and receivers. The second term enforces
the boundary condition for all source positions, while the third term ensures
that the outputs of the neural networks minimize the residual of the eikonal
equation over all sources and training points. Figure 1 summarizes the
proposed loss functions and the neural networks used.

The network parameters $\bm{\theta}_{\tau}$ and $\bm{\theta}_{v}$ are then
identified by solving the following minimization problem:

$\arg\min_{\bm{\theta}_{\tau},\bm{\theta}_{v}}\;\mathfrak{J}(\bm{\theta}_{\tau},\bm{\theta}_{v}).$
(9)

Figure 1:  Illustration of the PINNtomo algorithm: We use two neural networks
to approximate the traveltime factor $\hat{\tau}$ and the velocity $\hat{v}$.
The loss function used to train these networks minimizes the traveltime
mismatch at the receivers (first term), honors the boundary condition (second
term), and minimizes the residual of the eikonal equation (third term) at
selected training points within the computational domain. These quantities are
minimized for all source locations $\mathbf{x}_{s}=(x_{s},z_{s})$.

##  3 Numerical Tests

In this section, we test the proposed tomography workflow on cross-hole and
surface acquisition geometries. For both tests, we use neural networks
containing 10 hidden layers. For the traveltime factor network, each layer
contains 20 neurons, whereas, for the velocity network, each layer contains 10
neurons. We use a locally adaptive exponential linear unit (l-ELU) as the
activation function for the hidden layers. Locally adaptive activation
functions have been shown to achieve superior performance and convergence
speed over base methods [24]. We train the networks first using the Adam
optimizer with a batch size of 1024 for 500 epochs and then using the L-BFGS-B
optimizer until convergence. These hyper-parameters are chosen based on some
initial tests and kept fixed throughout the study to avoid the need for
tuning. The PINN framework is implemented using the SciANN package [25] – a
high level Tensorflow wrapper for scientific computations.

###  3.1 Example 1: Cross-hole tomography

First, we apply the tomography workflow on a $1\times 1$ km2 computational
domain with the velocity model shown in Figure 2. We consider 11 sources on
the left boundary of the model ($x=0$ km) uniformly spaced with an interval of
100 m, and 51 receivers on the right boundary of the model ($x=1$ km) with a
uniform spacing of 20 m. The observed traveltimes are computed through a
first-order factored eikonal solver using the fast sweeping method [23]. The
initial velocity model for this test is shown in Figure 2. It is worth noting
that, contrary to conventional tomography methods, here, the initial velocity
model is automatically selected based on the initialization of the parameters
of the velocity network. While the choice of the initial velocity model is
critical to the success of ray tomography, we find our formulation to be
agnostic to this choice. For illustration, the initial velocity model is
plotted using the same color scale as the true model.

Figure 3 shows the inverted velocity model, indicating that the long-
wavelength features of the true model have been well-recovered, as expected
from traveltime tomography. Figure 4 compares the velocity profiles at
$x\leavevmode\nobreak\ =\leavevmode\nobreak\ 0.4$ km and
$x\leavevmode\nobreak\ =\leavevmode\nobreak\ 0.8$ km between the initial,
true, and inverted velocities. We observe that the initial velocity model
evolved significantly to match the long-wavelength trend of the true velocity
model. Finally, in Figure 5 we show the final fit between the observed and the
predicted traveltimes at all receivers for the source at
$(x_{s},z_{s})=(0\leavevmode\nobreak\ \mathrm{km},0.5\leavevmode\nobreak\
\mathrm{km})$, indicating good agreement between the two.

Figure 2: The true (a) and initial (b) velocity models used for the cross-hole
tomography test. For this test, we use 11 equispaced sources on the left
boundary ($x=0$ km) of the model and 51 equispaced receivers on the right
boundary ($x=1$ km) of the model. Figure 3:  The inverted velocity model for
the cross-hole tomography test indicating reliable reconstruction of the long-
wavelength features of the true velocity model.

Figure 4: A comparison of the velocity profiles at $x\leavevmode\nobreak\
=\leavevmode\nobreak\ 0.4$ km (a) and $x\leavevmode\nobreak\
=\leavevmode\nobreak\ 0.8$ km (b) between the initial (dotted blue), true
(solid black), and inverted (dashed red) velocity models. Figure 5:  A
comparison between the observed traveltimes (solid black) and those predicted
using the neural network (dashed red) at the receivers located on the right
boundary of the model for the source at $(x_{s},z_{s})=(0\leavevmode\nobreak\
\mathrm{km},0.5\leavevmode\nobreak\ \mathrm{km})$.

###  3.2 Example 2: Surface tomography

Next, we test the performance of PINNtomo using a surface acquisition
geometry. We consider a $1\times 5$ km2 model with the velocity distribution
as shown in Figure 6. We consider 51 sources with a uniform spacing of 100 m
and 126 receivers with a uniform spacing of 40 m, both on the surface of the
model ($z=0$ km). Like earlier, the observed traveltimes are computed through
a first-order factored eikonal solver using the fast sweeping method. Again,
the initial velocity is a consequence of the random initialization of the
velocity network’s parameters. For illustration, the initial velocity model is
shown in Figure 6 using the same color scale as the true model.

Figure 7 shows the inverted velocity model, indicating that the long-
wavelength features, in particular for the velocity bump between $x=2-3$ km
have been well-recovered. We also observe a good fit between the observed and
predicted traveltimes at the receivers, as shown in Figure 8. It is worth
mentioning that for conventional tomography, a depth gradient for the initial
velocity is necessary, and the choice of it may bear some consequences on the
inverted velocities. However, as highlighted earlier, PINNtomo is agnostic to
the initial velocity distribution and does not require an initial model with a
depth gradient, even for surface tomography.

Figure 6: The true (a) and initial (b) velocity models used for the surface
tomography test. For this test, we use 51 equispaced sources and 126
equispaced receivers on the surface ($z=0$ km) of the model. Figure 7:  The
inverted velocity model for the surface tomography test indicating reliable
reconstruction of the long-wavelength features of the true model.  Figure 8:
A comparison between the observed traveltimes (solid black) and those
predicted using the neural network (dashed red) at the receivers located on
the surface for the source at $x\leavevmode\nobreak\ =3$ km.

##  4 Discussion and conclusions

In this work, we presented a novel approach to the traveltime tomography
problem. To this end, we use neural networks to approximate the traveltime
factor and the velocity fields, subject to the physics-informed regularization
based on the eikonal equation. Through tests on cross-hole and surface
acquisition geometries, we observe that the method is capable of reliably
recovering the long wavelength features of the velocity field. Conventional
techniques require an initial depth gradient for velocity that may affect the
inversion result. On the contrary, the performance of the proposed approach is
largely agnostic to the initial velocity model.

The proposed approach enjoys several advantages compared to conventional
tomography methods. Typically, conventional methods use some form of smoothing
regularization to compensate for the poor illumination of the model space due
to the acquisition geometry. These physics-agnostic regularizers end up
limiting the resolution of the inverted velocity model. On the contrary, using
the residual of the eikonal equation in the loss function, PINNtomo enforces a
physics-informed regularizer, which uses the actual wave propagation physics
to address the ill-posedness of the problem. Moreover, unlike conventional
methods, the proposed approach is mesh-free, and therefore, can be easily used
for models with irregular topography and can also accommodate source and
receiver points that are not on a regular grid. Furthermore, through the use
of transfer learning, the approach is well-suited to study temporal variations
in the near-surface using surface seismic tomography or at larger depths using
cross-hole tomography.

The method also outperforms conventional techniques in terms of both memory
and computational requirements. While the adjoint-state tomography method
overcame the memory limitation of the ray method by avoiding the need to
explicitly compute the Fréchet derivative matrix, it still depends on the size
of the velocity model, which could still be cumbersome for large 3D surveys.
Here, the memory required is dependent on the optimization batch size, which
is usually much lower than the entire model size. Moreover, since PINNtomo
uses Tensorflow at the backend, it allows easy deployment of computations
across a variety of platforms (CPUs, GPUs) and architectures (desktops,
clusters). Therefore, in this study, we used an NVIDIA Tesla P100 GPU that
required only $\sim 6$ minutes to invert the velocity model for the cross-hole
example and $\sim 20$ minutes for the surface tomography example.

## References

  * [1] Nicholas Rawlinson, S Pozgay, and S Fishwick.  Seismic tomography: a window into deep earth.  Physics of the Earth and Planetary Interiors, 178(3-4):101–135, 2010.
  * [2] Robert R Stewart.  Exploration seismic tomography: Fundamentals.  Society of Exploration Geophysicists, 1991.
  * [3] John Etgen, Samuel H Gray, and Yu Zhang.  An overview of depth imaging in exploration geophysics.  Geophysics, 74(6):WCA5–WCA17, 2009.
  * [4] Jean Virieux and Stéphane Operto.  An overview of full-waveform inversion in exploration geophysics.  Geophysics, 74(6):WCC1–WCC26, 2009.
  * [5] Cédric Taillandier, Mark Noble, Hervé Chauris, and Henri Calandra.  First-arrival traveltime tomography based on the adjoint-state method.  Geophysics, 74(6):WCB1–WCB10, 2009.
  * [6] ND Bregman, RC Bailey, and CH Chapman.  Crosshole seismic tomography.  Geophysics, 54(2):200–215, 1989.
  * [7] Joe Wong.  Crosshole seismic imaging for sulfide orebody delineation near Sudbury, Ontario, Canada.  Geophysics, 65(6):1900–1907, 2000.
  * [8] Richard D Rechtien, Roy J Greenfield, and Robert F Ballard Jr.  Tunnel signature prediction for a cross-borehole seismic survey.  Geophysics, 60(1):76–86, 1995.
  * [9] Tokuo Yamamoto.  Imaging the permeability structure within the near-surface sediments by acoustic crosswell tomography.  Journal of applied geophysics, 47(1):1–11, 2001.
  * [10] Hansruedi Maurer and Alan G Green.  Potential coordinate mislocations in crosshole tomography: Results from the Grimsel test site, Switzerland.  Geophysics, 62(6):1696–1709, 1997.
  * [11] David W Hyndman, Jerry M Harris, and Steven M Gorelick.  Inferring the relation between seismic slowness and hydraulic conductivity in heterogeneous aquifers.  Water Resources Research, 36(8):2121–2132, 2000.
  * [12] Yannick Choy Hing Ng, William Danovan, and Taeseo Ku.  The potential of seismic cross-hole tomography for geotechnical site investigation.  In E3S Web of Conferences, volume 92, page 18006. EDP Sciences, 2019\.
  * [13] JB Ajo-Franklin, J Peterson, J Doetsch, and TM Daley.  High-resolution characterization of a CO2 plume using crosswell seismic tomography: Cranfield, MS, USA.  International Journal of Greenhouse Gas Control, 18:497–509, 2013\.
  * [14] Alain Sei and William W Symes.  Gradient calculation of the traveltime cost function without ray tracing.  In SEG Technical Program Expanded Abstracts 1994, pages 1351–1354. Society of Exploration Geophysicists, 1994.
  * [15] Shingyu Leung, Jianliang Qian, et al.  An adjoint state method for three-dimensional transmission traveltime tomography using first-arrivals.  Communications in Mathematical Sciences, 4(1):249–266, 2006.
  * [16] Ting Ma and Zhongjie Zhang.  Topography-dependent eikonal traveltime tomography for upper crustal structure beneath an irregular surface.  Pure and Applied Geophysics, 172(6):1511–1529, 2015.
  * [17] Maziar Raissi, Alireza Yazdani, and George Em Karniadakis.  Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations.  Science, 367(6481):1026–1030, 2020.
  * [18] Francisco Sahli Costabal, Yibo Yang, Paris Perdikaris, Daniel E Hurtado, and Ellen Kuhl.  Physics-informed neural networks for cardiac activation mapping.  Frontiers in Physics, 8:42, 2020.
  * [19] Umair Waheed, Ehsan Haghighat, Tariq Alkhalifah, Chao Song, and Qi Hao.  Eikonal solution using physics-informed neural networks.  arXiv preprint arXiv:2007.08330, 2020.
  * [20] Jonathan D Smith, Kamyar Azizzadenesheli, and Zachary E Ross.  Eikonet: Solving the eikonal equation with deep neural networks.  IEEE Transactions on Geoscience and Remote Sensing, 2020.
  * [21] Ben Moseley, Andrew Markham, and Tarje Nissen-Meyer.  Solving the wave equation with physics-informed deep learning.  arXiv preprint arXiv:2006.11894, 2020.
  * [22] Chao Song, Tariq Alkhalifah, and Umair Bin Waheed.  Solving the frequency-domain acoustic VTI wave equation using physics-informed neural networks.  Geophysical Journal International, 2021.
  * [23] Sergey Fomel, Songting Luo, and Hongkai Zhao.  Fast sweeping method for the factored eikonal equation.  Journal of Computational Physics, 228(17):6440–6455, 2009.
  * [24] Ameya D Jagtap, Kenji Kawaguchi, and George Em Karniadakis.  Locally adaptive activation functions with slope recovery for deep and physics-informed neural networks.  Proceedings of the Royal Society A, 476(2239):20200334, 2020.
  * [25] Ehsan Haghighat and Ruben Juanes.  Sciann: A Keras/TensorFlow wrapper for scientific computations and physics-informed deep learning using artificial neural networks.  Computer Methods in Applied Mechanics and Engineering, 373:113552, 2021.