Stokes flow analogous to viscous electron current in graphene
Jonathan Mayzel Victor Steinberg Atul Varshney




^1Department of Physics of Complex Systems, Weizmann Institute of Science, Rehovot 76100,  Israel
^2The Racah Institute of Physics, Hebrew University of Jerusalem, Jerusalem 91904, Israel
^3Institute of Science and Technology Austria, Am Campus 1, 3400 Klosterneuburg, Austria





Electron transport in two-dimensional conducting materials such as graphene,  with dominant electron-electron interaction, exhibits unusual vortex flow that leads to a nonlocal current-field relation (negative resistance), distinct from the classical Ohm's law. The transport behavior of these materials is best described by low Reynolds number hydrodynamics, where the constitutive pressure-speed relation is Stoke's law. Here we report  evidence of such vortices observed in a viscous flow of Newtonian fluid in a microfluidic device consisting of a rectangular cavity-analogous to the electronic system.   We extend our experimental observations  to elliptic cavities of different eccentricities,  and validate them by  numerically solving  bi-harmonic equation obtained for the viscous flow with no-slip boundary conditions. We verify the existence of a predicted threshold at which vortices appear. Strikingly, we find that  a two-dimensional  theoretical model  captures the essential  features of  three-dimensional Stokes flow in experiments.


Stokes flow analogous to viscous electron  current in graphene
    Jonathan Mayzel,^1  Victor Steinberg^1,2 and Atul Varshney^1,3
    December 30, 2023
==================================================================




Electron transport in conducting materials is often hindered by  lattice disorder, impurities, and interactions with phonons and electrons <cit.>. In ultraclean two-dimensional (2D)  materials, the transport is primarily affected by electron-phonon (γ_p) and electron-electron (γ_ee) scattering processes; γ is the scattering rate of each process. In the limit γ_p≫γ_ee, a fast momentum relaxation of electrons results in a linear relationship between local current and applied electric field  called Ohm's law.  This relationship breaks down when the momentum exchange rate of electrons with each other is much faster than with the lattice <cit.> (γ_ee≫γ_p). In this regime, strongly interacting electrons move in a neatly coordinated manner which resembles  the flow of viscous fluids  <cit.>. Thus, hydrodynamic equations can be employed to describe the transport behaviour of electron   flow  provided mean free path (ℓ_ee) for momentum-conserving electron-electron collisions is the shortest length scale in the problem, i.e., ℓ_ee≪ w, ℓ_p; where w is the system size and ℓ_p is the mean free path for momentum-nonconserving electron-phonon collisions  <cit.>.  The first experimental signature of hydrodynamic electron flow  was obtained in  the measurements of differential resistance of electrostatically defined wires in a 2D electron gas in (Al,Ga)As heterostructures  <cit.>.



Recent development in the synthesis  of ultrapure  crystals has  facilitated   the investigation of the   viscous flow regime of electrons at elevated temperatures, so called viscous electronics  <cit.>. Experiments on viscous electron flow through narrow constrictions in doped single- and bi-layer graphene   reveal anomalous (negative) local resistance <cit.>. It was understood that the negative resistance may arise due to viscous shear flow which generates a vortex (whirlpool) and a backflow producing a reverse electric field that acts against  the applied field driving the source-drain current <cit.>.   Concurrently,   Levitov and Falkovich (L&F) developed a theoretical model based on Stokes flow of strongly interacting electrons <cit.>. They explored  three transport regimes of electron fluids, namely: Ohmic, mixed ohmic-viscous and  viscous, in an infinitely long 2D rectangular strip, with  point source and drain contacts located at the center, on  opposite sides  of the strip. Linearized 2D Navier-Stokes equation in the limit of low Reynolds number (Re≪1) is used to describe the transport regimes of incompressible electron fluids that yields bi-harmonic stream function, in contrast to the harmonic stream function for the Ohmic case. Here, the dimensionless Reynolds number (Re) defines the ratio of inertial and viscous stresses in fluids <cit.>. Furthermore,  they extended their  model  for   viscous  electron flow  in a three-dimensional (3D) conducting slab of small, finite thickness, and showed that the extra dimension in the Stokes equation translates into an effective resistance term  <cit.>.


A steady, viscous, and incompressible Newtonian fluid flow in a 3D slab of width w and height h at Re≪1 is described as η(∂_x^2+∂_y^2+∂_z^2)u_i=∇ P, where u_i(x,y,z)=u_i(x,y)u_z(z), η is the dynamic viscosity of the fluid and P is the pressure field <cit.>. Since a vertical parabolic velocity profile  u_z=z(h-z)/h^2 is globally uniform, and ∇ P(z)=const due to the uniformity of the velocity in both  spanwise (x) and streamwise (y) directions, one gets similar to what is suggested in ref.  <cit.> the following 2D linearized stationary Navier-Stokes equation after integration of the vertical  velocity profile: [η(∂_x^2+∂_y^2)-12η/h^2]u_i(x,y)=6∇ P. The threshold value to realize a vortex in a rectangular slab of thickness h is ϵ≡12w^2/h^2≤ 120 and it arises when the viscous shear force exceeds the wall friction (Ohmic) due to the boundaries. Therefore, in a three-dimensional system the viscous effect will be more pronounced when the system is thicker, at the same width, since the  friction arising from the top and bottom walls will be less significant.  The criterion value to observe vortices, i.e. ϵ<ϵ_c,  is estimated numerically for the rectangular slab in ref. <cit.> (see also Supp. Info. therein).

The generation of vortices in a fluid flow is typically associated with high-Re flow and inertial effects. However, such inferences are likely based on the incorrect notion that  low-Re flow is irrotational, which is only applicable to an ideal fluid without viscosity, where the Kelvin circulation theorem  is valid <cit.>. Thus, it is evident that a non-potential or rotational flow bears vorticity, but  to produce vortices in Stokes flow requires substantial efforts due to a strong dissipation of vorticity at Re≪1. Strikingly,  in a wall-dominated microfluidic channel flow,  vortices could be generated  at   low aspect ratio, w/h,  of the channel  despite a significant wall friction and viscous dissipation  <cit.>.


Here, motivated by the observation of   a distinct vortex flow in a strongly interacting electron system discussed above, we perform experiments  on a viscous flow of ordinary Newtonian fluid,  at low-Re, in a microfluidic  device consisting of a rectangular cavity, analogous to the 2D electronic system. Indeed, we observe a pair of symmetric vortices in the cavity region, when the geometrical criterion for the vortex observation ϵ< ϵ_c is satisfied, in agreement with the predictions <cit.> of L&F. Further, we expand our observations to elliptical cavities of different eccentricities (e)  and verify them with the analytical predictions <cit.> by  numerically solving the bi-harmonic equation obtained for the viscous flow with no-slip boundary conditions.




Results

Rectangular cavity. A long-exposure particle streak image in Fig. <ref>b (see also corresponding Supplementary Movie 1 <cit.>) illustrates the vortex flow at Re=0.07 in the rectangular cavity region (e → 1) of the device shown in Fig. <ref>a and described in Methods section. A pair of symmetric vortices appears in the cavity regions; here we show only  one side of the cavity. The flow in elliptic cavities is shown in Fig. <ref> and discussed in the next section. For the rectangular cavity, the streamlines and velocity field, obtained from micro-particle image velocimetry (μPIV), of the vortex flow (Fig. <ref>a) resemble the structures of the current streamlines and potential map as analysed by L&F for viscous flow in an infinite long 2D rectangular strip <cit.>. Furthermore, we obtain the location of the vortex center at x_0≈0.96w (see Figs. <ref>b and  <ref>a), which is in fair agreement with the analytical prediction <cit.> of x_0≈ w.










Elliptic cavities. Next, we investigated elliptical cavities with different eccentricities, e=√(1-b^2/a^2),  where a and b are the major and minor axes of an ellipse, respectively. Figures <ref>a,b show long-exposure images depicting the vortex flow for e=0.95 (at Re=0.08) and e=0.9 (at Re=0.07), respectively; see also  Supplementary Movies 2-4  <cit.>. The corresponding streamlines and velocity field are shown in Figs. <ref>b,c.  For e=0.95, the vortex center appears at x_0≈0.85w (see Figs. <ref>a and <ref>b), and for e=0.9 the vortex is pushed towards the cavity edge and its center is located at x_0≈0.95w (see Figs. <ref>b and <ref>c). However, for e=0.75 we do not observe a vortex (Fig. <ref>c and Supplementary Movie 4 <cit.>), which is discussed further.


Analytic solution and numerical simulations. Using  the Stokes equation, the incompressibility condition, and the fact that the system is two-dimensional, we introduce the streamline function ψ(x,y) via  two-dimensional velocity 𝐮(x,y)=𝐳×∇ψ, which reduces the Stokes equation to the bi-harmonic equation <cit.> (∂_x^2 + ∂_y^2)^2 ψ(x,y) = 0. We solve  the bi-harmonic equation analytically for disk geometry and numerically  for elliptical  geometries with no-slip boundary conditions u_y = u_x = 0, where u_y and u_x are transverse and longitudinal components of 𝐮, respectively. For a disk geometry (e=0) of radius R centered at the origin and the flow is injected (collected) with speed u_0 at y=R (y=-R), we do not observe vortices in the flow field (Supplementary Fig. 1(a),(b) <cit.>). However, the real part of the pressure field

    P(x,y) = η u_0 ( 8 i R^2/πz̅/(z̅^2 + R^2)^2),

where z̅ = x-iy,  exhibits non-trivial  patterns (Supplementary Fig. 1(c) <cit.>).


Further, we consider an ellipse with major axis a and minor axis b, where flow is injected and collected with a speed u_o at (0,± b).  Since vortices did not appear for a disk (e = 0), and did appear for an infinite rectangular strip (e → 1) we expect that at some critical value of eccentricity vortices would start to emerge.
The streamline function for the ellipse is obtained using a multi-step computational approach: First we define the vorticity of the flow field ω = ∇×𝐮, which is used to obtain the  complex harmonic function

    Φ(z̅) = 1/η P + i ω.

 Since 1/η P and ω form a Cauchy-Riemann pair, we can use conformal mapping to map Φ from a disk to an ellipse. The mapping formula is given by <cit.>:

    κ = √(m) sn(2K(m)/πsin^-1( z/√(a^2 - 1)), m),

where sn is the Jacobian elliptic sine, and K(m) is the complete elliptic integral of the first kind with modulus m. This formula maps the interior of an ellipse in the z-complex plane to the interior of a disk in the κ-complex plane. Then, from Φ in the elliptical geometry,  the real part of the vorticity is extracted as ω = (Φ) and then an inverse Laplacian is applied  on ω  that yields ψ up to a harmonic function ψ = ψ_0 + F, where ψ_0 is the result of the inverse Laplacian and F is a general harmonic function  determined by boundary conditions. The most general form of a harmonic function in elliptic coordinates (μ,ν) that preserves the symmetries in the problem is given by:

    F(μ, ν) = ∑_n=0^N(A_n e^-n μ + B_n e^n μ) cos nν.

The coefficients A_n and B_n are determined by fitting ψ to the no-slip boundary conditions, and N is the highest order for which the coefficients are computed. A detailed numerical scheme is delineated in Supplementary Notes <cit.>.  Results from numerics show that vortices exist for any ellipse with non-zero eccentricity, for instance, see Fig. <ref>b for e=0.97 and corresponding experimental results in Fig. <ref>a. In addition, we show  velocity fields  computed numerically and obtained experimentally for e=0.9 and 0.95 in Supplementary Figs. 2 and 3 <cit.>, respectively.


Next, we compute the distance of separatrix from the edge of the ellipse normalized with the cavity width, i.e. ξ/w, for different elliptical cavities.  Figure <ref> shows the variation of ξ/w as a function of e, ranging from a disk (e=0)  to an infinite rectangular strip (e→1), obtained numerically from the fitted harmonic function of different orders N (open symbols). Based on Fig. <ref>, it suggests that, for a disk, a vortex “lies" on the disk edge, and as  the disk is deformed into an elliptical geometry the vortex escapes the edge and develops inside the ellipse. The experimental results in Fig. <ref> (solid symbols) exhibit the same trend as from the numerical simulations. There is though a significant quantitative discrepancy between the experimentally measured ξ/w and the computed value.



Discussion

The quantitative discrepancy in the ξ/w  dependence on e  can likely be explained by the different flow regimes considered in the experiment and numerical simulations. The numerical simulations were conducted for the 2D viscous flow without any wall friction and with zero size inlet-outlet channels (in  analogy with the electrical transport corresponding to a pure viscous case with a zero contact size), whereas the experimental flow is 3D with finite-size inlet-outlet channels (0.5 mm channel width compared to the cavity width of ∼2 mm). The latter flow can be approximated by a quasi-2D viscous flow with the friction term ϵ=12w^2/h^2, corresponding to a mixed viscous-ohmic case.

 As discussed by L&F <cit.>, in the mixed viscous-ohmic regime, vortices will form when the  resistance to viscosity ratio ϵ is smaller than some critical value ϵ_c, which is geometry dependent. In our experimental system, ϵ<ϵ_c for the rectangular cavity with ϵ=12 (see Methods section) and for  elliptic cavities,  the value of ϵ_c is not known and  probably further depends on the eccentricity e; ϵ<ϵ_c holds true only for e→1. Thus, the experimental results for the rectangular strip ϵ/ϵ_c = 0.1 ≪ 1 match quite well with the analytical predictions (Fig. <ref>a and related text). For ellipses with varying eccentricity e, this critical threshold is probably lower and  reduces further with  decreasing e. For e > 0.9, ϵ_c is possibly still higher than 12, so we do observe vortices, even though their size and streamlines do not match with the numerical solution of the purely viscous case. A comparison between the numerical and experimental results of the velocity magnitude and flow field for ellipse with e=0.97 can be seen in Fig. <ref>. For ellipse with e ∼ 0.75, the critical threshold is probably lower than 12, therefore we do not observe vortices. This issue requires further investigation by solving the  equation for the mixed viscous-ohmic case for an ellipse with no-slip boundary conditions and by conducting experiments for a fixed value of the eccentricity e and various values of ϵ.




By using the condition for the vortex generation, which is derived in the quasi-2D approximation of a 3D channel flow, we assume that the 2D model is applicable to characterize the 3D  flow. The results of vortex observation in 3D flow indicate the validity of this assumption and its applicability.  To further validate  this assumption we compute the divergence of the 2D velocity field  close to the mid-height of the cavity (∇·𝐮) based on the experimental data. As can be seen in Fig. <ref>, the divergence field is  close to zero in the cavity regions which suggests a two-dimensional nature of the flow. The non-zero deviations in the divergence field are produced either by velocity fluctuations since they appear in plus/minus pairs or by the lower accuracy of PIV measurements in the region close to the inlet and outlet, where the velocities are high. A comparison between the divergence fields computed from  experimentally and numerically obtained velocity fields for e=0.97 is shown in Supplementary Fig. 4 <cit.>.


In conclusion, we have numerically examined  the two-dimensional  flow of a Newtonian fluid  by using an experimental, three-dimensional flow  at Re≪ 1; the latter  can be reduced to a quasi-2D flow with the wall friction term, as suggested in ref. <cit.>.   We have shown that vortices do appear in a rectangular cavity when the criterion ϵ<ϵ_c is satisfied, verifying the analytical prediction of L&F with a good accuracy. Furthermore, we expand the analytical predictions to disk and elliptic geometries, and validate the latter by experimental observations of vortices and the normalized vortex separatrix distance ξ/w. The observed quantitative discrepancies in  ξ/w as a function of the eccentricity may be explained by a reduction of the criterion value for the onset of the vortex observation ϵ_c with decreasing  eccentricity, e, from unity for a rectangular cavity down to zero for a circular disk.  The criterion ϵ<ϵ_c of obtaining a vortex breaks down at some value of e, which leads to a  disagreement with the theory. Thus, conducting experiments in low-Re Newtonian fluid flow allows us to inspect unusual vortex flow  properties which are experimentally unobservable in graphene or in other 2D electronic systems.
It is important to notice that in a Newtonian fluid flow the boundary conditions are  known,  simple, and verified, whereas in graphene they are less clear and may evolve due to edge currents, partial slip, <cit.> etc.



Methods

Experiments. The experimental system consists of a straight channel (40×0.5×2 mm^3) endowed with a long rectangular cavity (e→1) at the centre of the channel, as illustrated in Fig. <ref>a. The devices are prepared from transparent acrylic glass [poly(methyl methacrylate)]. The width and thickness of the cavity are w=h=2 mm, which gives ϵ=12. An aqueous glycerol solution (60% by weight) of  viscosity η=11 mPa·s  at 20C is used as a working fluid in the experiments. A smooth gravity-driven flow of the fluid is injected via the inlet into the straight channel, and its flow rate is varied by changing the fluid column height (see Fig. <ref>a).  The  fluid exiting from the outlet of the channel is weighed instantaneously W(t) as a function of time t by a  PC-interfaced weighing balance (BA-210S; Sartorius) with a sampling rate of 5 Hz and a resolution of 0.1 mg. The flow speed (u̅) is estimated as u̅=Q̅/(ρ wh), where time-averaged fluid discharge rate Q̅=Δ W/Δ t, fluid density  ρ=1156 Kg m^-3. Thus, the Reynolds number is defined as Re=wρu̅/η.

Imaging system. The fluid is seeded with 2 μ m sized fluorescent particles (G0200; Thermo Scientific) for  flow visualization.  The cavity region is imaged in the midplane directly via a microscope (Olympus IMT-2), illuminated uniformly with a light-emitting diode (Luxeon Rebel) at 447.5 nm wavelength, and a CCD camera (GX1920; Prosilica) attached to the microscope records 10^4 images with a spatial resolution of 1936×1456  pixel and at a rate of 9 fps. We employ micro-particle image velocimetry (μPIV) to obtain the spatially-resolved velocity field 𝐮=(u_y, u_x )  in the cavity region <cit.>. An interrogation window of 32× 32 pixel^2 (55× 55 μ m^2) with 50% overlap is chosen to procure 𝐮. Experiments are repeated on  cavities of different  e values.






10

ashcroft_solid_1976
Ashcroft, N. W. and Mermin, N. D.
Solid State Physics.
Cengage Learning, 1 edition,  (1976).

levitov_electron_2016
Levitov, L. and Falkovich, G.
Electron viscosity, current vortices and negative nonlocal resistance
  in graphene.
Nature Phys. 12, 672–676 (2016).

landau_fluid_1987
Landau, L. D. and Lifshitz, E. M.
Course of Theoretical Physics: Fluid Mechanics.
Butterworth-Heinemann, 2nd edition,  (1987).

batchelor
Batchelor, G. K.
An Introduction to Fluid Dynamics.
Cambridge University Press,  (1973).

huang_statistical_1987
Huang, K.
Statistical Mechanics.
Wiley, 2nd edition,  (1987).

andreev_hydrodynamic_2011
Andreev, A. V., Kivelson, S. A., and Spivak, B.
Hydrodynamic Description of Transport in Strongly Correlated
  Electron Systems.
Phys. Rev. Lett. 106, 256804 (2011).

lucas_hydrodynamics_2018
Lucas, A. and Fong, K. C.
Hydrodynamics of electrons in graphene.
J. Phys.: Condens. Matter 30, 053001 (2018).

molenkamp_1994
Molenkamp, L. W. and de Jong, M. J. M.
Observation of Knudsen and Gurzhi transport regimes in a
  two-dimensional wire.
Solid-State Electron. 37, 551–553 (1994).

de_jong_hydrodynamic_1995
de Jong, M. J. M. and Molenkamp, L. W.
Hydrodynamic electron flow in high-mobility wires.
Phys. Rev. B 51, 13389–13402 (1995).

bandurin_negative_2016
Bandurin, D. A., Torre, I., Kumar, R. K., Shalom, M. B., Tomadin, A., Principi,
  A., Auton, G. H., Khestanova, E., Novoselov, K. S., Grigorieva, I. V.,
  Ponomarenko, L. A., Geim, A. K., and Polini, M.
Negative local resistance caused by viscous electron backflow in
  graphene.
Science 351, 1055–1058 (2016).

moll_evidence_2016
Moll, P. J. W., Kushwaha, P., Nandi, N., Schmidt, B., and Mackenzie, A. P.
Evidence for hydrodynamic electron flow in PdCoO_2.
Science 351, 1061–1064 (2016).

kumar_superballistic_2017
Kumar, R. K., Bandurin, D. A., Pellegrino, F. M. D., Cao, Y., Principi, A.,
  Guo, H., Auton, G. H., Shalom, M. B., Ponomarenko, L. A., Falkovich, G.,
  Watanabe, K., Taniguchi, T., Grigorieva, I. V., Levitov, L. S., Polini, M.,
  and Geim, A. K.
Superballistic flow of viscous electron fluid through graphene
  constrictions.
Nature Phys. 13, 1182–1185 (2017).

braem_scanning_2018
	Braem, B. A., Pellegrino, F. M. D., Principi, A., Röösli, M., Gold, C., Hennel, S.,
	Koski, J. V., Berl, M., Dietsche, W., Wegscheider, W., Polini, M., Ihn, T.,
	and Ensslin, K.
	Scanning Gate Microscopy in a Viscous Electron Fluid.
	Phys. Rev. B 98, 241304(R) (2018).

torre_nonlocal_2015
Torre, I., Tomadin, A., Geim, A. K., and Polini, M.
Nonlocal transport and the hydrodynamic shear viscosity in graphene.
Phys. Rev. B 92, 165433 (2015).

pellegrino_electron_2016
Pellegrino, F. M. D., Torre, I., Geim, A. K., and Polini, M.
Electron hydrodynamics dilemma: Whirlpools or no whirlpools.
Phys. Rev. B 94, 155414 (2016).

falkovich_linking_2017
Falkovich, G. and Levitov, L.
Linking Spatial Distributions of Potential and Current in
  Viscous Electronics.
Phys. Rev. Lett. 119, 066601 (2017).

white
White, F. M.
Viscous Fluid Flow.
McGraw-Hill, 2 edition,  (1991).

sm see Supplemental Material at <https://nature.com/articles/s41467-019-08916-5#Sec10>.

tikhonov_equations_2011
Tikhonov, A. N. and Samarskii, A. A.
Equations of Mathematical Physics.
Dover Publications,  (2011).

selvadurai_partial_2000
Selvadurai, A. P. S.
Partial Differential Equations in Mechanics 2: The
  Biharmonic Equation, Poisson's Equation.
Springer-Verlag,  (2000).

nehari_conformal_2011
Nehari, Z.
Conformal Mapping.
Dover Publications,  (2011).

karageorghis_conformal_2008
Karageorghis, A. and Smyrlis, Y.-S.
Conformal mapping for the efficient MFS solution of Dirichlet
  boundary value problems.
Computing 83, 1–24 (2008).

allen_spatially_2016
Allen, M. T., Shtanko, O., Fulga, I. C., Akhmerov, A. R., Watanabe, K.,
  Taniguchi, T., Jarillo-Herrero, P., Levitov, L. S., and Yacoby, A.
Spatially resolved edge currents and guided-wave electronic states in
  graphene.
Nature Phys. 12, 128–133 (2016).

wagner_boundary_2015
Wagner, G.
Boundary Conditions for Electron Flow in Graphene in the
  Hydrodynamic Regime.
arXiv [cond-mat]: 1509.07113 (2015).

maier_fluctuations_2005
Maier, S., Sang, Y., Filleter, T., Grant, M., Bennewitz, R., Gnecco, E., and
  Meyer, E.
Fluctuations and jump dynamics in atomic friction experiments.
Phys. Rev. B 72, 245418 (2005).

filleter_friction_2009
Filleter, T., McChesney, J. L., Bostwick, A., Rotenberg, E., Emtsev, K. V.,
  Seyller, T., Horn, K., and Bennewitz, R.
Friction and Dissipation in Epitaxial Graphene Films.
Phys. Rev. Lett. 102, 086102 (2009).

thielicke_pivlab_2014
Thielicke, W. and Stamhuis, E.
PIVlab-Towards User-friendly, Affordable and Accurate
  Digital Particle Image Velocimetry in MATLAB.
J. Open Research Software 2 (2014).




   	
Acknowledgments

J.M. thanks Gregory Falkovich for introducing the problem and his guidance during the studies.  We thank Guy Han for technical support. A.V. acknowledges support from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 754411. This work was partially supported by the Israel Science Foundation (ISF; grant #882/15) and the Binational USA-Israel Foundation (BSF; grant #2016145).



