# Hydrodynamic electrons in Graphene: a viscous boundary-layer description.

Pedro Cosme<EMAIL_ADDRESS>Instituto de Plasmas e Fusão
Nuclear  Instituto Superior Técnico 1049-001 Lisboa, Portugal  João S. Santos
Instituto Superior Técnico 1049-001 Lisboa, Portugal  Hugo Terças  Instituto
de Plasmas e Fusão Nuclear  Instituto Superior Técnico 1049-001 Lisboa,
Portugal

(August 27, 2024)

###### Abstract

In this paper we dwell over the study of the boundary layer problem in a
hydrodynamical description of the electrons in gated graphene. It has been
verified experimentally that this fluid can display non-Poiseuille like flow
as reproduced in our numerical simulation. In fact, the velocity profile
displays a maximum value close to the boundary and then decreases as it
approaches the bulk of the graphene layer. This work aims to present a
satisfactory theoretical description of the boundary layer problem in
graphene. We found that by using the fluid equations and following a method
similar to that for deriving Blasius’ equation, a non-linear model can be
obtained whose solutions display the maximum values of velocity near the edges
of the graphene layer. We argue that such a non-monotonic model and behaviour
can shed some light on the subject of non-topological edge currents in
graphene.

##  I Introduction

The study of the collective state of carriers, be it electrons or holes, in
graphene is a difficult task to endeavour, particularly in a fully quantum
description with few works [1] dedicated to that goal. Instead, a growing
interest in semi-classical descriptions can be found. In particular, several
hydrodynamics descriptions of the mean field quantities have been brought
forward in the late years [2, 3, 4, 5, 6], with some experimental works also
taking the first steps [7, 8].

For such hydrodynamical models, a central point of discussion is the nature of
the flow, that is, what kind of velocity profile is foreseeable under a
specific set of conditions, as with regular hydrodynamics this is ought to be
primarily determined by the viscosity of the fluid. Nonetheless, to calculate
the viscosity of such system from first principles is a complex task, and even
if the results for undoped graphene are relatively consensual [9, 10] the same
does not happen for the doped scenario, seemingly the most relevant situation
for technological applications. Moreover, the boundary conditions for the
Fermi liquid at the edges of a graphene layer are still not completely
understood and the naive approach of a no-slip condition might not be the most
accurate one[11, 12].

Some authors have committed their analysis to a highly viscous regime, leading
to a Poiseuille flow [7, 13]. However, for such a parabolic flow to develop
there are two conditions required; namely, that the flow is quasi-stationary,
i.e. it has a small Womersley number [14], and that the transverse dimension
of the system is smaller than the boundary layer[15], i.e. in a narrow
channel. These condition might not be met in some conditions, either when
dealing with regimes with high frequency or wide channels, for instance in the
case of graphene nano-antennae operating at THz [2, 16, 17]. Thus, the
assessment of the boundary layer thickness for the electronic flow is of
central importance for study of the hydrodynamic models of conduction in
graphene.

In this work we present a description of the boundary layer for graphene
electrons, in the degenerate limit. We start by establishing the hydrodynamic
model for the conducting electrons in gated graphene [18] and, in the
following sections, proceed to derive a nonlinear model for the velocity
profile across the channel, comparing it with the simulations from which we
also extract scaling laws for the boundary layer thickness.

##  II Electron Fluid in Graphene

In the highly degenerate regime, the electrons in graphene are all described
by a hydrodynamic model resorting to the Drude mass [3, 5]

$m^{\star}=\frac{\hbar k_{F}}{v_{F}}=\frac{\hbar\sqrt{\pi n}}{v_{F}},$ (1)

where $\hbar k_{F}$ is the Fermi momentum and $n$ is the carrier number
density. It is important to note that, since the electronic fluid is
compressible, the effective mass is not a conserved quantity, contrary to
customary fluids.

From a kinetic description one can derive [18] the conservation laws for
number of particles, momentum and energy as

$\frac{\partial n}{\partial
t}+\bm{\nabla}\\!\cdot\\!\frac{\mathbf{p}}{m^{\star}}=0,$ (2a)
$\frac{\partial\mathbf{p}}{\partial
t}+\bm{\nabla}\\!\cdot\\!\left(\frac{\mathbf{p}\otimes\mathbf{p}}{nm^{\star}}+\mathbb{P}\right)+en\bm{\nabla}\phi=0.$
(2b) $\frac{\partial\mathcal{E}}{\partial
t}+\bm{\nabla}\cdot\left(\frac{3}{2}\mathcal{E}\mathbf{v}\right)+e\frac{\mathbf{p}}{m^{\star}}\cdot\bm{\nabla}\phi=0,$
(2c)

respectively. Where $\mathbf{p}=nm^{\star}\mathbf{v}$ is the fluid two-
dimensional momentum density, $\mathcal{E}$ the energy density, $\mathbb{P}$
the pressure stress tensor, and $\phi$ the electrical potential. In the
present work however, we will not make use of the energy conservation equation
(2c) since the continuity (2a) and Cauchy (2b) equations will form a closed
set of equations by themselves as will become evident.

The two-dimensional pressure term in (2b) can be, as customary, parcelled out
between the hydrostatic and viscous terms and cast as

$\bm{\nabla}\cdot\mathbb{P}=\bm{\nabla}P-\eta_{s}\nabla^{2}\mathbf{v}-\zeta\bm{\nabla}\left(\bm{\nabla}\cdot\mathbf{v}\right),$
(3)

with the pressure in the 2D Fermi-Dirac system [19, 20, 21]

$P=\frac{\hbar v_{F}}{3\pi}\big{(}\pi n\big{)}^{\frac{3}{2}},$ (4)

and where the coefficients $\eta^{s}$ and $\zeta$ represent the shear and bulk
viscosity, respectively. Nonetheless, the bulk viscosity in graphene is
extremely low [9] and will, therefore, be neglected in this discussion.

Also, in the gradual channel approximation, the Fermi level varies slowly
along the graphene, and so the field can be considered uniform and akin to a
plane capacitor; with a capacitance per area $C\approx\varepsilon/d_{0}$ where
$d_{0}$ is the graphene/gate separation and $\varepsilon$ the medium
permittivity. Thus, the applied bias potential $U_{\text{gate}}$ is directly
linked to the carrier density as $U_{\text{gate}}=en/C$ and so, the force term
in (2c) is simply

$\bm{\nabla}\phi=\bm{\nabla}U_{\text{gate}}=\frac{ed_{0}}{\varepsilon}\bm{\nabla}n.$
(5)

Thus, since both the pressure and force terms are functions of the number
density we group them as an effective pressure $\mathcal{P}$ such that

$\bm{\nabla}\mathcal{P}=\bm{\nabla}P+\frac{e^{2}d_{0}n}{\varepsilon}\bm{\nabla}n$
(6)

##  III Boundary Layer Problem

Although several works have been devoted to the viscous flow of electrons on
graphene, akin to Poiseuille flow [7, 22, 8], for the flow to be in such
conditions, of full development and entirely under the viscous effects, the
longitudinal dimension of the flow must be significantly larger than the
cross-section, so that it might be larger than the entrance length[15]; that
is, conform to the narrow channel approximation. In wider channels, or less
viscous situations, we can expect that the flow is separated in two distinct
regions, the uniform bulk and the boundary layer. Boundary layer is defined as
the region of the flow where the velocity field is not uniform due to the
influence of the viscous stress exerted on the fluid from the interface; it is
bounded by a typical thickness $\delta$[15], as schematized on Fig. 1.

$x$$y$$\delta(x)$$U_{0}$$U(x)$$u(x,y)\quad v(x,y)$

Figure 1: Schematic representation of the boundary layer after the injection
of a flow $U_{0}$ near an interface. The flow is divided between two regions
at a distance $\delta(x)$. Outside the boundary layer the flow
$U(x)\rightarrow U_{\infty}$ is assumed to be unperturbed and parallel to the
interface whereas inside the boundary layer the flow is described by
$\mathbf{v}=u(x,y)\mathbf{\hat{x}}+v(x,y)\mathbf{\hat{y}}$

While for incompressible fluids with a no-slip boundary condition Blasius
theory [15, 23] provides us a good description of the characteristic thickness
of the boundary layer, the situation of Fermi liquid herein considered is more
demanding. Not only is the fluid compressible due to the non-constant mass,
but also the fluid at the edges may not be completely at rest and, on the
contrary, slip along the boundary[11]. Moreover, our simulations of the flow
have showed that the velocity profile is not monotonic, displaying a maximum
value near the edge before converging to the free flow limit, as can be seen
in Fig.2. This can be exploited to easily characterize the distinct regions of
the flow, identifying the boundary layer thickness with the position of
maximum velocity $\delta_{\rm max}$.

Figure 2: Velocity profile of viscous electronic flow
($\nu_{s}\equiv\eta_{s}/m^{\star}n=4/Lv_{0}$). Transition from Poiseuille to
non-convex flow width the varying aspect ratio.

###  III.1 Heuristic scaling law for $\delta_{\rm max}$

From a dimensional analysis approach one can argue that the boundary layer
thickness should follow a scaling law

$\delta_{\rm max}\sim\frac{x}{{\rm Re}_{x}^{\alpha}},$ (7)

with ${\rm Re}_{x}$ the Reynolds number of the flow at position $x$ and
$\alpha$ an exponent to be determined. Solving numerically the model (2c),
with no slip boundary condition and sufficiently wide aspect ratio, provides
the velocity profiles to which (7) can the fitted to retrieve $\alpha$. Our
findings show that it increases marginally with kinematic viscosity and Mach
number as patent on Fig. 3.

Figure 3: Fitted exponent parameter for $\delta_{\rm max}$ for varying
kinematic viscosity and Mach number.

Although useful for empirical usage this numerical approach does not grant us
full understanding of the velocity profile, thus, an analytical description of
such profile was pursuit.

###  III.2 Analytical description of velocity profile

In order to obtain a more in depth description of the flow we followed the
usual procedure to determine the laminar boundary layers equations [15, 23].
From (2c) one gets,

$\frac{u}{2}\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial
y}-\frac{u}{2}\frac{\partial v}{\partial
y}-\frac{\eta_{s}}{nm^{\star}}\frac{\partial^{2}u}{\partial
y^{2}}+\frac{1}{nm^{\star}}\frac{\partial\mathcal{P}}{\partial x}=0,$ (8a)
$\text{and}\quad\frac{1}{nm^{\star}}\frac{\partial\mathcal{P}}{\partial y}=0.$
(8b)

Moreover, the introduction of a potential field $\psi$, such that
$nu=\partial_{y}\psi$ and $nv=-\partial_{x}\psi$, automatically solves the
steady-state of (2a).

Equation (8b) imposes that the pressure and, consequently trough (4), the
number density are solely functions of $x$. Thus, the Dorodnitsyn-Howarth [15]
transformation for the $y$ coordinate yields

$\eta=\frac{n(x)}{n_{\infty}}\sqrt{\frac{U_{\infty}}{2\nu x}}y$ (9)

which, along with the scaling of the potential field as

$\frac{\psi}{\sqrt{2\nu xU_{\infty}}n_{\infty}}=f(\eta),$ (10)

converts, after some algebra, Eq. (8a) to the nonlinear third-order equation

$\frac{n^{2}}{n_{\infty}^{2}}f^{\prime\prime\prime}+ff^{\prime\prime}=\frac{x}{n}\frac{dn}{dx}\left({f^{\prime}}^{2}+2\frac{S^{2}}{U_{\infty}^{2}}\right)$
(11)

with
$S=\left(v_{F}^{2}/2+e^{2}dv_{F}\sqrt{n_{0}}/\varepsilon\hbar\sqrt{\pi}\right)^{1/2}$
the sound speed of the plasmonic fluid. This equation ought to be supplemented
with the appropriate boundary conditions $f^{\prime}(0)=0$ and
$f^{\prime}(\infty)=1$ representing the non-slip condition and asymptotic
value of the velocity, respectively. Contrarily to regular boundary layers,
Eq. (11) can have overshoot solutions, i.e. with $f^{\prime}>1$ in some
region, correctly reproducing the behavior observed in the simulations.

Let Eq. (11) be recasted as

$af^{\prime\prime\prime}+ff^{\prime}=\beta\left(f^{\prime 2}+2m^{2}\right),$
(12)

with the obvious identifications, to ease its analysis. Note that this
equation reduces to Blasius equation for $a=1$ and $\beta=0$. All three
parameters can, in first order, be considered constant as they only vary along
the length of the channel. Regarding its values one expects, $a\simeq 1$, as
there should not occur a sizable variation of number density across the
channel width. In the case of $\beta$, it encodes the density drop along the
channel; thus, should be negative, as the density should not increase, but it
is also forseeable that $|\beta|\ll 1$. Finally, $m$, inverse Mach number, can
be tuned between about 10 and 100, due to the properties of the graphene
electron fluid in study[18].

In order to analyse (12) and to obtain some insight on its solutions, the
differential transformation method (DTM) [24] can be performed, from which one
obtains a recursive formula for the coefficients of the formal power series
expansion of $f(\eta)=\sum_{k=0}^{\infty}F_{k}\eta^{k}$ yielding

$F_{k+3}=-\frac{k!}{a(k+3)!}\Bigg{[}\sum_{q=0}^{k}(k-q+2)(k-q+1)F_{q-1}F_{k-q}\\\
-2\beta m^{2}-\beta\sum_{q=0}^{k}(q+1)(k-q+1)F_{q+1}F_{k-q+1}\Bigg{]}.$ (13)

As the condition $f^{\prime}(\infty)=1$ is incompatible with an expansion
about the origin it was replaced by $f^{\prime\prime}(0)=\tau$, i.e. fixing
the shear at the edge. Consequently, the first coefficients can be retrieved
as

$\displaystyle F_{0}$ $\displaystyle=0$ $\displaystyle,\quad F_{2}$
$\displaystyle=\frac{\tau}{2}$ $\displaystyle,\quad F_{4}$
$\displaystyle=\frac{m^{2}\beta}{12a}$ $\displaystyle F_{1}$ $\displaystyle=0$
$\displaystyle,\quad F_{3}$ $\displaystyle=\frac{m^{2}\beta}{3a}$
$\displaystyle,\quad F_{5}$
$\displaystyle=-\frac{\tau^{2}}{120}+\frac{\beta(2m^{2}+\tau^{2})}{60a}.$

It becomes clear that $m^{2}\beta$ is a figure of merit in the general
behaviour of the solutions, setting them apart from the classical boundary
layer profiles.

Solving, numerically, the boundary layer equation (12) we are able to
recapture the characteristic non-monotonic behaviour observed in the fluid
simulations. This aspect is clear in Fig. 5 where the influence of parameters
variation can also be observed.

Figure 4: Argument of maximum of $f^{\prime}$, $\eta_{\rm
max}\propto\delta_{\rm max}$, as function of parameters $\beta$ and $m$.

Figure 5:  Solutions of boundary layer equation
$af^{\prime\prime\prime}+ff^{\prime}=\beta\left(f^{\prime 2}+2m^{2}\right)$.
Sweeping parameters $a$ (left panel), $\beta$ (centre panel) and $m$ (right
panel). Reference values of parameters when not being varied: $a=1$;
$\beta=-2\\!\times\\!10^{-5}$; $m=30$.

##  IV Discussion

Although not expected for regular fluids, the peaked profile of the velocity
field is a prevailing feature of our simulations of graphene intermediate
viscous regimes; signalling a transition from Ohmic to fully viscous flow.
Thus, while the non-monotonicity of $f^{\prime}(\eta)$ impedes the usual
technique of characterizing the boundary layer by the displacement or momentum
thicknesses, $\delta_{\rm max}$ can be consistently used as $\delta(x)$
indicator.

However, the simulation approach to determine $\delta_{\rm max}$ shows that de
boundary layer lacks universality; since $\delta_{\rm max}$ power law has a
small, yet noticeable, dependence on $S/v_{0}$ and $\nu$ itself. Nonetheless,
we can estimate $0.4\lesssim\alpha\lesssim 0.6$, in line with the expected
square root signature. Furthermore, the functional dependence of $\delta_{\rm
max}$ with flow parameters becomes evident with our model (11). Not only the
Mach number appears explicitly, but also, the density drop $\beta$ must be
dependent on the Mach and Reynolds numbers, albeit in a still unknown way.

Regarding the generalized boundary layer model, we found good agreement
between the solutions of (12) and the simulated velocity profiles,
particularly in replicating the velocity maximum. Showing that, even a minute
pressure drop, around $\beta\approx-10^{-4}$, can be responsible for the non-
convex flow profile.

##  V Conclusions

In this work, we derive a generalization of boundary layer theory for viscous
flow in a hydrodynamical description of the massless electrons in graphene.
This model, characterized by the parameters $a=n/n_{\infty}$, $\beta=d\log
n/d\log x$ and $m=S/U_{\infty}$, accurately predicts the non-monotonic
velocity profile observed in our performed fluid simulations and the
experimental observations reported in the literature. Our results, pointing
towards the existence of a velocity peak near the edge, show that edge
currents may show up in the system. This may lead to a better understanding of
the topic of non-topological edge currents in graphene, which have been
recently found in the literature [25].

Moreover, our findings, both from the analytical model and the simulations,
provide a valuable way to estimate the boundary layer thickness in
experimental configurations as $\delta_{\rm max}\sim x/{\rm Re}_{x}^{\alpha}$
with $0.4\lesssim\alpha\lesssim 0.6$. Such scaling law can be put in use when
designing or analysing graphene circuitry.

Looking towards the future, the effects of a slip length type boundary, i.e.
where $f^{\prime}(0)\neq 0$, may introduce new velocity profile contours. We
are confident that the model presented in this work is able to cope with such
conditions and remain faithful to the experimental and numerical results;
still, further investigation in such scenario is required.

Finally, the inclusion of odd viscosity in situations with broken parity
symmetry (namely with an applied magnetic field) may also be a feature of
future interest.

###### Acknowledgements.

The authors acknowledge the funding provided by Fundação para a Ciência e a
Tecnologia (FCT-Portugal) through Grant No. PD/BD/150415/2019 and through the
exploratory project UTA-EXPL/NPN/0038/2019.

## References

  * [1] J. L. Figueiredo, J. P. S. Bizarro, and H. Terças, “Wigner description of massless Dirac plasmas,” pp. 1–12, 12 2020.
  * [2] P. Cosme and H. Terças, “Terahertz Laser Combs in Graphene Field-Effect Transistors,” ACS Photonics, vol. 7, pp. 1375–1381, 6 2020.
  * [3] A. Lucas and K. C. Fong, “Hydrodynamics of electrons in graphene,” Journal of Physics: Condensed Matter, vol. 30, p. 053001, 2 2018.
  * [4] B. N. Narozhny, I. V. Gornyi, A. D. Mirlin, and J. Schmalian, “Hydrodynamic Approach to Electronic Transport in Graphene,” Annalen der Physik, vol. 529, p. 1700043, 11 2017.
  * [5] D. Svintsov, V. Vyurkov, V. Ryzhii, and T. Otsuji, “Hydrodynamic electron transport and nonlinear waves in graphene,” Physical Review B - Condensed Matter and Materials Physics, vol. 88, no. 24, pp. 28–32, 2013.
  * [6] L. Barletti, “Hydrodynamic equations for electrons in graphene obtained from the maximum entropy principle,” Journal of Mathematical Physics, vol. 55, no. 8, 2014.
  * [7] J. A. Sulpizio, L. Ella, A. Rozen, J. Birkbeck, D. J. Perello, D. Dutta, M. Ben-Shalom, T. Taniguchi, K. Watanabe, T. Holder, R. Queiroz, A. Principi, A. Stern, T. Scaffidi, A. K. Geim, and S. Ilani, “Visualizing Poiseuille flow of hydrodynamic electrons,” Nature, vol. 576, pp. 75–79, 12 2019\.
  * [8] M. J. H. Ku, T. X. Zhou, Q. Li, Y. J. Shin, J. K. Shi, C. Burch, L. E. Anderson, A. T. Pierce, Y. Xie, A. Hamo, U. Vool, H. Zhang, F. Casola, T. Taniguchi, K. Watanabe, M. M. Fogler, P. Kim, A. Yacoby, R. L. Walsworth, M. M. Fogler, P. Kim, A. Yacoby, and R. L. Walsworth, “Imaging viscous flow of the Dirac fluid in graphene,” Nature, vol. 583, pp. 537–541, 7 2020\.
  * [9] A. Principi, G. Vignale, M. Carrega, and M. Polini, “Bulk and shear viscosities of the two-dimensional electron liquid in a doped graphene sheet,” Physical Review B, vol. 93, p. 125410, 3 2016.
  * [10] B. N. Narozhny and M. Schütt, “Magnetohydrodynamics in graphene: Shear and Hall viscosities,” Physical Review B, vol. 100, p. 035125, 7 2019\.
  * [11] E. I. Kiselev and J. Schmalian, “Boundary conditions of viscous electron flow,” Physical Review B, vol. 99, p. 35430, 1 2019.
  * [12] S. Danz and B. N. Narozhny, “Vorticity of viscous electronic flow in graphene,” 2D Materials, vol. 7, p. 035001, 4 2020.
  * [13] M. Polini and A. K. Geim, “Viscous electron fluids,” Physics Today, vol. 73, pp. 28–34, 6 2020.
  * [14] R. Moessner, P. Surówka, and P. Witkowski, “Pulsating flow and boundary layers in viscous electronic hydrodynamics,” Physical Review B, vol. 97, p. 161112, 4 2018.
  * [15] H. Schlichting and K. Gersten, Boundary-Layer Theory.  Berlin, Heidelberg: Springer Berlin Heidelberg, 6 2017.
  * [16] F. Bianco, V. Miseikis, D. Perenzoni, C. Coletti, M. Perenzoni, and A. Tredicucci, “Antenna-coupled graphene field-effect transistors as a terahertz imaging array,” IEEE Transactions on Terahertz Science and Technology, vol. XX, no. X, p. 1, 2020.
  * [17] S. Dash and A. Patnaik, “Behavior of graphene based planar antenna at microwave and terahertz frequency,” Photonics and Nanostructures - Fundamentals and Applications, vol. 40, p. 100800, 7 2020.
  * [18] P. Cosme and H. Terças, “Hydrodynamical study of terahertz emission in magnetized graphene field-effect transistors,” Applied Physics Letters, vol. 118, p. 131109, 3 2021.
  * [19] L. D. Landau and E. M. Lifshitz, Courses on theoretical physics: Statistical Physics, Part 1, vol. 5.  Pergamon Press, 3rd ed., 1980.
  * [20] G. Giuliani and G. Vignale, Quantum Theory of the Electron Liquid.  Cambridge University Press, 2005.
  * [21] A. J. Chaves, N. M. Peres, G. Smirnov, and N. Asger Mortensen, “Hydrodynamic model approach to the formation of plasmonic wakes in graphene,” Physical Review B, vol. 96, no. 19, pp. 1–18, 2017.
  * [22] I. Torre, A. Tomadin, A. K. Geim, and M. Polini, “Nonlocal transport and the hydrodynamic shear viscosity in graphene,” Physical Review B - Condensed Matter and Materials Physics, vol. 92, pp. 1–12, 10 2015.
  * [23] C. F. Dewey and J. F. Gross, “Exact Similar Solutions of the Laminar Boundary-Layer Equations,” Advances in Heat Transfer, vol. 4, no. C, 1967\.
  * [24] C. Bervillier, “Status of the differential transformation method,” Applied Mathematics and Computation, vol. 218, no. 20, pp. 10158–10170, 2012\.
  * [25] A. Aharon-Steinberg, A. Marguerite, D. J. Perello, K. Bagani, T. Holder, Y. Myasoedov, L. S. Levitov, A. K. Geim, and E. Zeldov, “Long-range nontopological edge currents in charge-neutral graphene,” Nature, vol. 593, no. May, 2020.