Stokes flow analogous to viscous electron current in graphene

TL;DR

This paper demonstrates the analogy between viscous electron flow in graphene and viscous fluid flow in microfluidic devices, showing vortex formation and validating a 2D theoretical model for 3D Stokes flow.

Contribution

It provides experimental evidence and numerical validation of vortex behavior in viscous flows, extending the analogy to electronic systems in graphene.

Findings

Vortices observed in viscous fluid flow analogous to electronic systems.

A threshold for vortex formation is identified.

A 2D model accurately captures 3D Stokes flow features.

Abstract

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.

Equations

$$P(x,y)=\eta u_0\Re \left(\frac{8i{R}^2}{\pi }\frac{\bar{z}}{({\bar{z}}^2+R^2{)}^2}\right),$$

$$\Phi (\bar{z})=\frac{1}{\eta }P+i\omega .$$

$$\kappa =\sqrt[4]{m}\text{ sn}\left(\frac{2K(m)}{\pi }{\sin }^{-1}\left(\frac{z}{\sqrt{a^2-1}}\right),m\right),$$

$$F(\mu ,\nu )=\sum _{n=0}^N\left(A_n{e}^{-n\mu }+B_n{e}^{n\mu }\right)\cos n\nu .$$

Full text

Electron transport in conducting materials is often hindered by lattice disorder, impurities, and interactions with phonons and electrons 1. In ultraclean two-dimensional (2D) materials, the transport is primarily affected by electron-phonon ($\gamma_{\mathrm{p}}$) and electron-electron ($\gamma_{\mathrm{ee}}$) scattering processes; $\gamma$ is the scattering rate of each process. In the limit $\gamma_{\mathrm{p}}\gg\gamma_{\mathrm{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 2 ($\gamma_{\mathrm{ee}}\gg\gamma_{\mathrm{p}}$). In this regime, strongly interacting electrons move in a neatly coordinated manner which resembles the flow of viscous fluids 3, 4. Thus, hydrodynamic equations can be employed to describe the transport behaviour of electron flow provided mean free path ($\ell_{\mathrm{ee}}$) for momentum-conserving electron-electron collisions is the shortest length scale in the problem, i.e., $\ell_{\mathrm{ee}}\ll w,\ell_{\mathrm{p}}$; where $w$ is the system size and $\ell_{\mathrm{p}}$ is the mean free path for momentum-nonconserving electron-phonon collisions 3, 4, 5, 6, 7. 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 8, 9.

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 10, 11, 12, 13. Experiments on viscous electron flow through narrow constrictions in doped single- and bi-layer graphene reveal anomalous (negative) local resistance 10. 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 14, 10, 2, 15. Concurrently, Levitov and Falkovich (L$\&$F) developed a theoretical model based on Stokes flow of strongly interacting electrons 2, 16. 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 ($\mathrm{Re}\ll 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 ($\mathrm{Re}$) defines the ratio of inertial and viscous stresses in fluids 4. 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 2.

A steady, viscous, and incompressible Newtonian fluid flow in a 3D slab of width $w$ and height $h$ at $\mathrm{Re}\ll 1$ is described as $\eta(\partial_{x}^{2}+\partial_{y}^{2}+\partial_{z}^{2})u_{i}=\nabla P$, where $u_{i}(x,y,z)=u_{i}(x,y)u_{z}(z)$, $\eta$ is the dynamic viscosity of the fluid and $P$ is the pressure field 3. Since a vertical parabolic velocity profile $u_{z}=z(h-z)/h^{2}$ is globally uniform, and $\nabla P(z)=\mathrm{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. 2 the following 2D linearized stationary Navier-Stokes equation after integration of the vertical velocity profile: $[\eta(\partial_{x}^{2}+\partial_{y}^{2})-12\eta/h^{2}]u_{i}(x,y)=6\nabla P$. The threshold value to realize a vortex in a rectangular slab of thickness $h$ is $\epsilon\equiv 12w^{2}/h^{2}\leq 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. $\epsilon<\epsilon_{c}$, is estimated numerically for the rectangular slab in ref. 2 (see also Supp. Info. therein).

The generation of vortices in a fluid flow is typically associated with high-$\mathrm{Re}$ flow and inertial effects. However, such inferences are likely based on the incorrect notion that low-$\mathrm{Re}$ flow is irrotational, which is only applicable to an ideal fluid without viscosity, where the Kelvin circulation theorem is valid 3, 17. 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 $\mathrm{Re}\ll 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 2.

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-$\mathrm{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 $\epsilon<\epsilon_{c}$ is satisfied, in agreement with the predictions 2 of L$\&$F. Further, we expand our observations to elliptical cavities of different eccentricities ($e$) and verify them with the analytical predictions 2 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. 1b (see also corresponding Supplementary Movie 1 18) illustrates the vortex flow at $\mathrm{Re}=0.07$ in the rectangular cavity region ($e\to 1$) of the device shown in Fig. 1a 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. 2 and discussed in the next section. For the rectangular cavity, the streamlines and velocity field, obtained from micro-particle image velocimetry ($\mu$PIV), of the vortex flow (Fig. 3a) 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 16. Furthermore, we obtain the location of the vortex center at $x_{0}\approx 0.96w$ (see Figs. 1b and 3a), which is in fair agreement with the analytical prediction 16 of $x_{0}\approx w$.

Elliptic cavities. Next, we investigated elliptical cavities with different eccentricities, $e=\sqrt{1-b^{2}/a^{2}}$, where $a$ and $b$ are the major and minor axes of an ellipse, respectively. Figures 2a,b show long-exposure images depicting the vortex flow for $e=0.95$ (at $\mathrm{Re}=0.08$) and $e=0.9$ (at $\mathrm{Re}=0.07$), respectively; see also Supplementary Movies 2$-$4 18. The corresponding streamlines and velocity field are shown in Figs. 3b,c. For $e=0.95$, the vortex center appears at $x_{0}\approx 0.85w$ (see Figs. 2a and 3b), and for $e=0.9$ the vortex is pushed towards the cavity edge and its center is located at $x_{0}\approx 0.95w$ (see Figs. 2b and 3c). However, for $e=0.75$ we do not observe a vortex (Fig. 2c and Supplementary Movie 4 18), 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 $\psi(x,y)$ via two-dimensional velocity $\mathbf{u}(x,y)=\mathbf{z}\times\nabla\psi$, which reduces the Stokes equation to the bi-harmonic equation 19, 20 $\left(\partial_{x}^{2}+\partial_{y}^{2}\right)^{2}\psi(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 $\mathbf{u}$, 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) 18). However, the real part of the pressure field

[TABLE]

where $\bar{z}=x-iy$, exhibits non-trivial patterns (Supplementary Fig. 1(c) 18).

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,\pm b)$. Since vortices did not appear for a disk ($e=0$), and did appear for an infinite rectangular strip ($e\to 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 $\omega=\nabla\times\mathbf{u}$, which is used to obtain the complex harmonic function

[TABLE]

Since $\frac{1}{\eta}P$ and $\omega$ form a Cauchy-Riemann pair, we can use conformal mapping to map $\Phi$ from a disk to an ellipse. The mapping formula is given by 21, 22:

[TABLE]

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 $\kappa$-complex plane. Then, from $\Phi$ in the elliptical geometry, the real part of the vorticity is extracted as $\omega=\Re(\Phi)$ and then an inverse Laplacian is applied on $\omega$ that yields $\psi$ up to a harmonic function $\psi=\psi_{0}+F$, where $\psi_{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 ($\mu,\nu$) that preserves the symmetries in the problem is given by:

[TABLE]

The coefficients $A_{n}$ and $B_{n}$ are determined by fitting $\psi$ 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 18. Results from numerics show that vortices exist for any ellipse with non-zero eccentricity, for instance, see Fig. 4b for $e=0.97$ and corresponding experimental results in Fig. 4a. In addition, we show velocity fields computed numerically and obtained experimentally for $e=0.9$ and 0.95 in Supplementary Figs. 2 and 3 18, respectively.

Next, we compute the distance of separatrix from the edge of the ellipse normalized with the cavity width, i.e. $\xi/w$, for different elliptical cavities. Figure 5 shows the variation of $\xi/w$ as a function of $e$, ranging from a disk ($e=0$) to an infinite rectangular strip ($e\to 1$), obtained numerically from the fitted harmonic function of different orders $N$ (open symbols). Based on Fig. 5, 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. 5 (solid symbols) exhibit the same trend as from the numerical simulations. There is though a significant quantitative discrepancy between the experimentally measured $\xi/w$ and the computed value.

Discussion

The quantitative discrepancy in the $\xi/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 $\sim 2$ mm). The latter flow can be approximated by a quasi-2D viscous flow with the friction term $\epsilon=12w^{2}/h^{2}$, corresponding to a mixed viscous-ohmic case.

As discussed by L&F 2, in the mixed viscous-ohmic regime, vortices will form when the resistance to viscosity ratio $\epsilon$ is smaller than some critical value $\epsilon_{c}$, which is geometry dependent. In our experimental system, $\epsilon<\epsilon_{c}$ for the rectangular cavity with $\epsilon=12$ (see Methods section) and for elliptic cavities, the value of $\epsilon_{c}$ is not known and probably further depends on the eccentricity $e$; $\epsilon<\epsilon_{c}$ holds true only for $e\to 1$. Thus, the experimental results for the rectangular strip $\epsilon/\epsilon_{c}=0.1\ll 1$ match quite well with the analytical predictions (Fig. 3a 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$, $\epsilon_{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. 4. For ellipse with $e\sim 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 $\epsilon$.

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 ($\nabla\cdot\mathbf{u}$) based on the experimental data. As can be seen in Fig. 6, 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 18.

In conclusion, we have numerically examined the two-dimensional flow of a Newtonian fluid by using an experimental, three-dimensional flow at $\mathrm{Re}\ll 1$; the latter can be reduced to a quasi-2D flow with the wall friction term, as suggested in ref. 2. We have shown that vortices do appear in a rectangular cavity when the criterion $\epsilon<\epsilon_{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 $\xi/w$. The observed quantitative discrepancies in $\xi/w$ as a function of the eccentricity may be explained by a reduction of the criterion value for the onset of the vortex observation $\epsilon_{c}$ with decreasing eccentricity, $e$, from unity for a rectangular cavity down to zero for a circular disk. The criterion $\epsilon<\epsilon_{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-$\mathrm{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, 23, 24, 25, 26 etc.

Methods

Experiments. The experimental system consists of a straight channel ($40\times 0.5\times 2$ mm3) endowed with a long rectangular cavity ($e\to 1$) at the centre of the channel, as illustrated in Fig. 1a. 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 $\epsilon=12$. An aqueous glycerol solution ($60\%$ by weight) of viscosity $\eta=11$ mPa$\cdot$s at $20\degree\mathrm{C}$ 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. 1a). 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 ($\bar{u}$) is estimated as $\bar{u}=\bar{Q}/(\rho wh)$, where time-averaged fluid discharge rate $\bar{Q}=\overline{\Delta W/\Delta t}$, fluid density $\rho=1156$ Kg m*-3*. Thus, the Reynolds number is defined as $\mathrm{Re}=w\rho\bar{u}/\eta$.

Imaging system. The fluid is seeded with $2\leavevmode\nobreak\ \mu 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 104 images with a spatial resolution of $1936\times 1456$ pixel and at a rate of 9 fps. We employ micro-particle image velocimetry ($\mu$PIV) to obtain the spatially-resolved velocity field $\mathbf{u}=(u_{y},u_{x})$ in the cavity region 27. An interrogation window of $32\times 32$ pixel2 ($55\times 55\leavevmode\nobreak\ \mu m^{2}$) with $50\%$ overlap is chosen to procure $\mathbf{u}$. Experiments are repeated on cavities of different $e$ values.