Fluidity onset in graphene

Viscous electron fluids have emerged recently as a new paradigm of strongly-correlated electron transport in solids. Here we report on a direct observation of the transition to this long-sought-for state of matter in a high-mobility electron system in graphene. Unexpectedly, the electron flow is found to be interaction-dominated but non-hydrodynamic (quasiballistic) in a wide temperature range, showing signatures of viscous flows only at relatively high temperatures. The transition between the two regimes is characterized by a sharp maximum of negative resistance, probed in proximity to the current injector. The resistance decreases as the system goes deeper into the hydrodynamic regime. In a perfect darkness-before-daybreak manner, the interaction-dominated negative response is strongest at the transition to the quasiballistic regime. Our work provides the first demonstration of how the viscous fluid behavior emerges in an interacting electron system.

E lectron fluids, an exotic state of matter in which electron-electron (ee) interactions dominate transport, have been long anticipated theoretically [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15] but until recently they were far from experimental reality. This situation is currently changing owing to the discovery of new materials in which ee interactions are particularly strong or momentum relaxation due to disorder and phonons is weak. The inventory of experimental systems that can host viscous e-fluids, as we will call them for brevity, has been steadily growing in the last few years [16][17][18][19] , stimulating wide interest in their properties. E-fluids may exhibit new behaviors such as vortices 20,21 , whirlpools 16 , superballistic transport 22,23 , Poiseuille flow 10,11,13,14,18 , anomalous heat conduction 17 , and viscous magnetotransport 24,25 . The questions about the genesis of e-fluids, on the other hand, received relatively little attention. How does an electron system enter the fluid state? What happens when l ee becomes comparable or larger than the system dimensions? What is the relation between electric current and potential at the transition? All these questions are at present poorly understood: neither there exists a detailed theory treating both ballistic and viscous electron regimes on equal footing, nor any systematic experimental study of the transition has been performed. Searching for the fluidity onset is the subject of this work.
So far, the behavior of e-fluids was mostly discussed deep in the hydrodynamic regime, where the mean free path l ee was the shortest lengthscale of the system. However, the experimental conditions are usually such that l ee , tunable by varying temperature T, is either comparable or at most a few times smaller than the system dimensions, putting the experimentally investigated e-fluids close to the onset of fluidity. As we will show below, this regime hosts an interaction-dominated quasiballistic state, which exhibits a negative voltage response similar to that observed at not-too-high T in the ref. 16 . The negative response arises because ambient carriers, as a result of momentumconserving collisions with injected carriers, are blocked from reaching voltage probes. Furthermore, the negative response is enhanced by "memory effects", so that it may exceed the negative response in the viscous state 26 . Thus, the interaction-dominated quasiballistic state, while quite distinct from the viscous fluid state, can in some cases serve as a proxy for the latter.
Graphene offers a convenient venue for this study. First, due to their exceptional cleanness and weak electron-phonon (el-ph) coupling, state-of-the-art graphene devices support micrometerscale ballistic transport with respect to momentum-nonconserving collisions over a wide range of temperatures 27 , from liquid-helium to room T. Second, above the temperatures of liquid nitrogen, ee collisions become the dominant scattering mechanism, so that the behavior of the electron system resembles that of viscous fluids 16,23 . Third, l ee in graphene can be varied over a wide range 23 by changing the carrier density n and T. This enables a smooth transition (or, more precisely, a crossover) between single-particle ballistic and viscous transport regimes, allowing us to track how the electron system enters the collective fluid state.

Results
Experimental data. We explore the onset of the hydrodynamic state by studying graphene devices in the so-called vicinity geometry 16 , illustrated in Fig. 1a: The current I is injected through a narrow contact into a wide graphene channel, and a local potential is probed at a small distance x from the injector. The main result of our study is that the vicinity resistance R v = V/I reaches an extreme negative value at the onset of fluidity. In particular, this behavior manifests itself most clearly through the temperature dependence of R v (Fig. 1b, c), with the quasiballistic and hydrodynamic regimes occurring at low and high T, respectively. We will show that the deep minimum at intermediate temperatures in the R v (T) dependences is the hallmark of the transition. Furthermore, we will demonstrate that this transition can be conveniently quantified by the electron Knudsen number taking values Kn ( 1 and Kn > 1 in the hydrodynamic and quasiballistic transport regimes, respectively, and approaching unity at the fluidity onset.
Importantly, the negative sign of R v , observed across the entire transition, signals that ee interactions dominate transport in both the quasiballistic and hydrodynamic regimes. The hydrodynamic regime, where theory predicts dR v /dT > 0 20,28 , occurs only at high enough temperatures and low enough carrier densities. This regime is preceded by an extended quasiballistic regime with dR v / dT < 0, discussed in detail below. The occurrence of two distinct interaction-dominated regimes in a 2D electron system is a surprising finding, which is of interest from a fundamental perspective and important for possible applications.
To explore the onset of the fluid state experimentally, we fabricated high-quality devices based on bilayer graphene (BLG) encapsulated between hexagonal boron nitride (for details, see Methods). The latter provides a clean environment for graphene's electron system ensuring micrometer-scale ballistic transport with respect to extrinsic momentum-non-conserving scatterering. The a Optical photograph of one of our devices on which the measurement geometry is indicated: current I is injected into the graphene channel through a 300 nm contact and the voltage drop is measured at a distance x from the injection point. Device width W is 2.3 μm. b, c Temperature dependence of the vicinity resistance measured experimentally and computed theoretically for different carrier densities n in bilayer graphene. The most negative value occurs at the fluidity onset, Kn $ 1, where Kn is the Knudsen number, (1) devices were shaped in a form of dual-gated multiterminal Hall bars (Fig. 1a), allowing us to study the distance-dependent potential anticipated at the transition upon varying the carrier densities n. The dual-gated design allowed us to maintain zero displacement between the graphene layers, so that one could tune the Fermi energy ε F in BLG without altering its band structure (opening the band gap). We have strategically chosen the BLG system because it ε F varies with n stronger than in monolayer graphene (MLG) (n vs. n 1/2 ). The standard dependence l ee = ℏv F ε F /(k B T) 2 translates into the scaling l ee ∼ n 3/2 , which is much faster than the n 1/2 dependence in MLG. This allowed us to explore a wider range of l ee than in MLG by varying the carrier density for a given T (see below), providing a convenient knob to tune the Kn value and probe the quasiballistic-to-hydrodynamic transition 29 .
Notably, the signal measured in the vicinity configuration contains a non-negligible offset due to momentum-nonconserving scattering (by phonons and/or disorder) which we further refer to as an Ohmic contribution. To distill the viscous contribution, we employed the approach introduced in the ref. 16 in which the Ohmic term, expressed as bρ, was subtracted from the measured vicinity signal, assuming the additive behavior of these contributions 28 . Here ρ = ρ(n, T) is the BLG sheet resistance measured in the conventional four-terminal geometry and b is the geometric factor that depends on sample dimensions and the distance between the injection point and the voltage probe 16,28 (for example, b = 0.1 for the measurementent configuration shown in Fig. 1a). As discussed below, the procedure of subtracting the Ohmic contribution, while somewhat ad hoc, can be justified for the geometry of our experiment. Below we refer to this adjusted vicinity resistance using the same notation R v unless stated otherwise. Figure 1b shows R v as a function of T measured in one of our BLG devices. Far away from the charge neutrality point (CNP) and at liquid helium T, R v is positive for all experimentally accessible n. When the temperature is increased, R v rapidly drops, reverses its sign, reaches a minimum and then starts to grow. Figure 2a details this observation by mapping R v on the (n, T)plane. The non-monotonic dependence R v vs. T is observed for all n, whereby the temperature at which R v dips, grows with increasing n (red dashed line).
To understand this nonmonotonic behavior, we first consider the limiting cases: the hydrodynamic regime l ee ( x, realized at large T, and the free-particle regime l ee ) W, realized at the lowest T (here W is the device width). In the hydrodynamic regime, negative R v arises as a result of viscous entrainment by the injected current of the fluid in adjacent regions 16,20,28 . In the free-particle regime, positive R v is expected from single-particle ballistic transport due to reflection of injected carriers from the opposite boundaries 30,31 . Therefore, the sign of R v must change from negative to positive upon lowering T, as indeed seen in the data shown in Fig. 1a. Furthermore, the hydrodynamic R v is proportional to viscosity 20,28 , giving the dependence R v $ l ee ðTÞ. The quantity l ee (T) increases as T decreases, leading to increasingly more negative R v . The non-monotonic temperature dependence R v (T), implied by these observations, is indeed seen in our measurements (Figs. 1b,2a).
Importantly, in between the free-particle regime l ee ) W and the hydrodynamic regime l ee ) x lies an interesting regime x < l ee < W that has hitherto been ignored in the literature. This intermediate regime, which for the lack of a better name will be called "quasiballistic", features an interaction-dominated response of a non-hydrodynamic nature, since the mean free path l ee is greater than the distance from the injector to the probe. Conspicuously, R v remains negative in this regime. However, since now R v $ 1=l ee ðTÞ, the sign of dR v /dT is reversed compared to the hydrodynamic regime. The negative sign of R v can be understood by considering injected carriers that travel over a large distance of the order of l ee > x and then scatter off ambient thermal carriers. After scattering, some of the injected carriers make it back into the probe, creating a positive contribution to R v . Simultaneously, some of the ambient carriers, through scattering off the injected carriers, are blocked from reaching the probe. This process creates a negative contribution to R v . Detailed analysis shows that the latter contribution dominates 26 , giving rise to negative R v . As T increases, R v grows progressively more negative until the point l ee = x, where the hydrodynamic behavior sets in and the sign of the T dependence is reversed. Interestingly, in the quasiballistic regime, the value |R v | decreases with n and grows with T, in qualitative agreement with the behavior of a MLG R v at not-too-high T found in the ref. 16 . This suggests a possible resolution of the conundrum posed by the findings of ref. 16 , in which a hydrodynamic-like negative R v was found to depend on n and T differently from what is expected in the hydrodynamic regime.
Theory and comparison with experiment. To capture all these different regimes in a single model, we employ the kinetic equation for quasiparticles in the graphene Fermi liquid. Transport in the geometry of Fig. 1 is described by solving the kinetic equation in an infinite strip of width W: −∞ < x < ∞, 0 < y < W, with diffuse boundary conditions at the strip edges y = 0, W. Current I is injected through a point-like source at x = y = 0 and is drained on the far left, x = −∞. We find the potential at (x, 0) by evaluating the particle flux entering the probe (for details, see Methods). At low temperatures the ee rate γ ee is small, and the ee collision term can be ignored 12 .
The model then describes ballistic particles bouncing between the strip edges, as illustrated in the upper inset of Fig. 3b. The net flux of particles into the probe then gives a positive value R v = V p (x)/I. At high T, on the other hand, the ee collision term dominates, and the distribution function approaches the local equilibrium. The resulting hydrodynamic behavior is then described by the Stokes equation that states the balance between the viscous friction and electric forces: eE/m = −ν▽ 2 v. (The latter follows directly from the Eq. (12) of Methods, multiplied by p, integrated over momenta and combined with an expression for the stress tensor obtained from 1/γ ee expansion.) In this case, we obtain R v $ η=ðnexÞ 2 , where η is the dynamic viscosity given by η ¼ 1 4 m Ã nv F l ee and m* is the carrier effective mass 20,22 . The single parameter γ ee allows us to explore both the ballistic and viscous regime through the dependence of R v on T and n. Carrier dynamics in the quasiballistic regime is shown schematically in the lower inset of Fig. 3b.
In Fig. 1b, c we compare the experimental data for R v vs. T with the results of our modeling, assuming the ee collision rate that depends on T and n as hγ ee $ T 2 e =ε F 12 . For bilayer graphene, the Fermi energy ε F is related to the carrier density as n = m * ε F /(πℏ 2 ), where m * = 0.033 m e . The two panels flaunt good qualitative agreement; namely, our theory captures the main experimental features: positive R v at small T that rapidly drops with increasing T and monotonically grows with n, so that the minima and sign changes in R v occur at higher T for larger n. Furthermore, our model reproduces some of the more subtle features of the data. For example, the nodes in R v vs. T shift to higher T and the minima to lower T, as the distance to the probe x increases, see Fig. 3. An overall agreement is also found for the full R v (n, T) maps shown in Fig. 2a, b that become near-identical after rescaling the T axis.

Discussion
In our analysis, for simplicity, we disregarded the Ohmic effects due to the el-ph scattering. This is a reasonable starting point since the el-ph scattering mean free path l el−ph is considerably larger than l ee at the temperatures of interest (for details, see Methods). However, the flow can be distorted by the Ohmic effects at the lengthscales set by ξ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi which lies between l ee and l el−ph 20,21 . Thus caution must be exercised even when the el-ph scattering is weak. The procedure of extracting the viscous contribution by subtracting the Ohmic contribution is expected to work well so long as the Ohmic effects  Fig. 1a. Temperature enters through the ee scattering rate hγ ee ' T 2 e =ε F 17 . The top axis corresponds to the electron Knudsen number Kn = l ee /x calculated for the blue curve. The purple rectangles and arrows mark, respectively, the sign change and minima in the R v (T) dependences. Upper inset: Schematics of electron transport in the ballistic regime. The potential at the source is transported by carriers throughout the sample and into the probe. Lower inset: Schematics of the voltage sign change in the quasiballistic regime. Collisions between injected carriers (red arrows) and ambient thermal carriers (green arrows) diminish the number of thermal carriers reaching the probe. This reduces the potential, which reverses its sign and becomes negative. Arrows link the insets with the corresponding portions of the R v (T) dependence do not distort the current flow at the lengthscales which are being probed, i.e. when ξ exceeds the distance to the probe x ≈ 1 μm. Estimates show that the inequality ξ ) x holds at not-too-high temperatures, i.e. in the quasiballistic regime. At the fluidity onset, identified above as the turning point in the R v (T) dependence, for the estimated typical values l ee ≲0:2 μm and l elÀph $ 3 μm, the lengthscale ξ can become comparable to x. However, an analysis based on the Stokes equation indicates that, for the geometry of our experiment, the Ohmic and viscous contributions remain approximately additive even for ξ < x (for details, see Methods). We therefore believe that the subtraction procedure provides a reasonable approximation in the entire range of temperatures and dopings.
We also note that Figs. 2, 3 exhibit some discrepancy between the values of T at which theoretical and experimental R v reach the minimum. This is not particularly surprising given the simplistic expression of γ ee $ T 2 used in the model. Since γ ee is the only relevant temperature-dependent parameter in the model, the quantitative agreement can be improved through revising the dependence γ ee vs. T. Indeed, there are various effects that can give rise to deviation from the standard Fermi-liquid T 2 dependence. One is the logarithmic enhancement of the quasiparticle decay rate due to collinear ee collisions [32][33][34] . However, it is probably an unlikely culprit, since collinear collisions do not lead to angular relaxation. At the same time, recent analysis 35 indicates that the effective γ ee that determines electron viscosity depends on the lifetimes of the odd-m angular harmonics, m = ±3, ±5,..., which relax considerably slower than the Fermi-liquid T 2 estimate would suggest. Accounting for this effect could, effectively, extend the quasiballistic behavior to higher temperatures, which would improve the agreement with the observed dependence R v (T). Detailed analysis of these rates and of their impact on R v is beyond the scope of this work.
The experimental and theoretical R v (T) exhibit two prominent features: R v first changes sign from positive to negative and then passes through a deep minimum. Should the sign change or the minimum be taken as the signature of the onset of fluidity? That question can be answered with the help of the data presented in Fig. 3, demonstrating that R v is a non-trivial function of both l ee / W and l ee /x. We note in that regard that the sign reversal of theoretically computed R v occurs at Kn ) 1, that is inside the quasiballistic regime, for all values of x (Fig. 3b). Indeed, R v in Fig. 3b changes sign at T = 20 K which for a given n translates into l ee ≈ 10 μm, a length scale significantly greater than the values x $ 1 À 2 μm for this device. On the other hand, the most negative R v in Fig. 3b is found at Kn = 1-3, which corresponds to x $ l ee <W. Since in the hydrodynamic regime R v is proportional to η and thus should drop with increasing T, we infer that it is the condition Kn $ 1 (where R v is most negative) that describes the fluidity onset. Furthermore, R v is expected to be negative in the quasiballistic regime 26 when Kn > 1, so it is indeed the drop of | R v | with temperature, rather than the sign reversal, that marks the onset of the viscous flow.
Experimental observation of this anomalous behavior at the onset of the fluid state enables a direct electrical measurement of the mean free path l ee and electron viscosity. Good qualitative agreement of the experimental data and our theoretical model suggests further opportunities to study the physics of e-fluids, in particular the electron transport in the presence of magnetic field and/or confining potential, obstacles, funnels and electron pumps. Our work clearly shows that the initial deviation from the ballistic behavior observed experimentally in different systems 13,14,16,18,23 may be due to an entry into the interactiondominated "quasiballistic" regime rather than the true onset of electron fluidity. It requires higher temperatures and the observation of the behavior consistent with viscosity gradually decreasing with increasing T to ascertain that the Navier-Stokes description can be applied.

Methods
Device fabrication. Our devices were made of bilayer graphene encapsulated between ≈50 nm-thick crystals of hexagonal boron nitride (hBN). The hBNgraphene-hBN heterostructures were assembled using the dry-peel technique described elsewhere 27,36 and deposited on top of an oxidized Si wafer (290 nm of SiO 2 ) which served as a back gate. After this, a PMMA mask was fabricated on top of the hBN-graphene-hBN stack by electron-beam lithography. This mask was used to define contact areas to graphene, which was done by dry etching with fast selective removal of hBN 37 . Metallic contacts (usually, 5 nm of chromium followed by 50 nm gold) were then deposited onto exposed graphene edges that were a few nm wide. As the next step, another round of electron-beam lithography was used to prepare a thin metallic mask (40 nm Al) which defined a multiterminal Hall bar. After this, reactive ion plasma etching translated the shape of the metallic mask into encapsulated graphene. The Al mask also served as a top gate, in which case Al was wet-etched near the contact leads to remove the electrical contact to graphene.
Distilling the hydrodynamic contribution in the presence of Ohmic effects.
Here we assess the accuracy of the approach used in the main text to separate the viscous and Ohmic contributions to the R v signal. In this approach, it was assumed that the contributions are approximately additive, and thus the viscous contribution can be distilled by subtracting the (suitably scaled) Ohmic resistivity measured in a four-probe setup.
The validity of the additivity assumption can be verified using an exact solution of the hydrodynamic equations for current injected in a halfplane. The hydrodynamic approach applies when the ee mean free path is smaller than the elph scattering mean free path, l ee ( l elÀph . At the scales larger than l ee the electron flow satisfies the Stokes equation with an Ohmic term added to describe momentum relaxation: Taking a curl and defining κ 2 = ρ(en) 2 /η = 1/ξ 2 , we obtain the equation on the stream function: where v = ∇ × (ψz). Following 21 , we consider the flow in a half-plane y > 0 generated by the a point source on the boundary at x = 0: ψ x (x, 0) = δ(x)I/ne. The stream function in this case has the form Ae Àjkjy þ ð1 À AÞe Àqy where we defined q ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi k 2 þ κ 2 p >0. The stream function can be used to evaluate the potential. Plugging Eq. (4) in Eq. (2), we see that only the first (harmonic) term in the stream function contributes the potential: ð5Þ The yet-undetermined quantity A(k) depends on the type of boundary condition. The no-stress boundary condition at y = 0, which reads ψ yy (x, 0) = 0, yields A(k) = q 2 /κ 2 = 1 + k 2 /κ 2 . Remarkably, the exact potential is a sum of the viscous and Ohmic contributions, with each contribution unaffected by the presence of the other contribution in this case: where L is the system size. The subtraction procedure employed in analyzing the measurements is exact at all distances for the no-stress boundary condition. For the no-slip boundary condition, on the other hand, the additivity is only an approximate property. In this case, ψ y (x, 0) = 0 gives The last term in this expression gives a contribution which depends both on viscosity and resistivity. As illustrated in Fig. 4, this contribution is non-negligible at distances x ' ξ, where R v changes sign. Its magnitude, however, is small (under 10-15% of the total potential). Therefore, disregarding this contribution should provide a reasonably good approximation. Yet, this conclusion is almost certainly geometry-sensitive, being valid for the point source at a halfplane edge but not necessarily for other geometries.
Estimates of the electron-phonon scattering mean free path. Electron-phonon scattering rate in graphene was discussed mostly for the single-layer case [38][39][40] . Here we modify this analysis for the bilayer case. The value of the mean free path l el−ph is used in the main text to determine the lengthscales at which the the el-ph scattering does not distort the carrier flow.
We use the standard deformation potential Hamiltonian H elÀph ¼ Z d 2 rψ y ðr; tÞD∇uðr; tÞψðr; tÞ; where u(r, t) is the lattice displacement vector, D is the deformation potential coupling constant, ω k = s|k| is the phonon frequency, and ρ is the surface mass density of graphene sheet. Plugging these quantities into the Golden Rule for the elph emission rate gives where θ is the angle parameterizing the Fermi surface, and the deformation potential matrix element equals V fi ¼ , with the overlap 〈ψ f | ψ i 〉 = cos(θ p′ − θ p ) accounting for the chirality of charge carriers. Here p and p′ are electron momenta, and k = p − p′. (Parenthetically, for monolayer graphene, the cos factor is to be replaced with cos((θ p′ − θ p )/2).) The density of final states equals ν = m * /(2πℏ 2 ), where m * = 0.033 m e is the carrier effective mass; since electron-phonon scattering preserves carrier spin and valley index, the relevant degeneracies are not included in ν.
Phonon absorption is described by a similar expression with N ph (k) + 1 replaced by N ph (k). Since temperatures of interest are considerably larger than the Bloch-Gruneisen temperature T BG = ℏsk F , we can approximate the Bose factors N ph (k) and N ph (k) + 1 as T/ℏω k . Plugging N ph (k) + N ph (k) + 1 ≈ 2T/ℏω k in the expression for dΓ and replacing k with p − p′, gives dΓ ¼ cos 2 ðθÞ dθ 2π Then the transport scattering rate equals The electron-phonon mean free path is given by l el−ph = v/Γ tr , where v = ℏk F /m * is the carrier velocity. For bilayer graphene, we assume surface mass density ρ = 2 × 7.6 × 10 −7 kg/m 2 , the speed of sound s = 2 × 10 4 m/s. In single-layer graphene, transport measurements are consistent with deformation potential D of the order of 20 eV, see, e.g., ref. 41 . For bilayer graphene, ab-initio calculations 42 yield D = 15 eV. Assuming D in the range of 15-20 eV, we arrive at l el−ph of the order of 3 μm for typical experimental conditions. Table 1 provides a summary of the results for the single-layer and bilayer graphene. These estimates are in agreement with the el-ph scattering rates extracted from the temperature dependence of the four-probe resistance reported in the ref. 16 .
Details of the theoretical model. To describe the ballistic and viscous regimes on equal footing and provide a link between them, we use the kinetic Boltzmann equation for quasiparticles at the Fermi surface. Expanded to linear order in the deviation δf from the equilibrium Fermi-Dirac distribution, Boltzmann equation reads v∇ x δf ðp; xÞ À I ee ðδf ðp; xÞÞ ¼ Jðp; xÞ: The collision operator I ee in (12) describes scattering between single-particle states via momentum-conserving ee collisions. Near the Fermi surface, the distribution can be parameterized by the standard ansatz δf ðpÞ ¼ À ∂f 0 ∂ε χðpÞ, where the energy dependence in χ can be ignored on the account of fast quasiparticle thermalization by collinear scattering at the 2D Fermi surface 32,33 . We analyze the angular dependence χ(θ), where the angle θ parameterizes the Fermi surface and p ¼ ðcosθ; sinθÞ is the unit vector along the carrier momentum. We assume that all non-conserved angular harmonics of χ(θ) relax with equal rates γ ee = v F /l ee ,whereas the three angular harmonics corresponding to the conserved net momentum and particle number do not relax. The operator I ee , linearized in χ, therefore takes the form 13,22 : The angular brackets denote angular averaging over θ′.   12) is to be furnished with the boundary conditions describing momentum relaxation at the strip edges. We assume that particles are scattered diffusely, following Lambert's law. Hence the edges y = 0 and y = W effectively become isotropic current sources. At y = 0, we write χðθ>0; xÞ ¼ Jðθ; xÞ þ 1 2 The choice of the coefficient 1/2 in the second term is dictated by current conservation. Indeed, for an isotropic distribution of outgoing particles, χ(θ > 0) = χ 0 , the outgoing particle flux, ν R π 0 v F sinθχðθÞdθ=2π, is given by νv F χ 0 /π. Here ν is the density of particle states, and v F is the Fermi velocity. In the absence of current injection, this quantity must be equal to the incoming flux which is given by the integral in the second term. Similarly, an isotropic current source attached to the boundary, I(x, 0), is described via J(θ, x) = πI(x)/(eνv F ) in the Eq. (14). For the opposite orientation of the boundary, y = W, positive and negative angle values in the Eq. (14) are to be interchanged.
In general, distribution of particles in the Knudsen regime is not represented by a local equilibrium Fermi function, and a local chemical potential cannot be introduced. This poses a difficulty in relating the signal on a probe contact to the distribution function. To resolve this, we adopt the model of a probe which is commonly used to describe leads in mesoscopic circuits, see e.g., 43 . A probe is a perfect absorber for nonequilibrium carriers, which are equilibrated inside the probe and subsequently re-emitted into the fluid with an isotropic angular distribution. If the open-circuit condition is maintained in the probing circuit, the potential on the probe is proportional to the influx F of charge carriers into the probe, Since outgoing charge carriers are in equilibrium with the probe potential V p , they are characterized by the distribution function χ = eV p , so that the outgoing flux is νeV p /π. Balancing these fluxes, one finds the probe potential In the hydrodynamic regime, the distribution function is given by an equilibrium expression, which can be related to the local electric potential, χ(θ, x) ≈ eϕ(x). In this limit, one finds V p (x) = ϕ(x). For a generic nonequilibrium distribution, however, the relation between the local potential ϕ(x) and the probe signal V p (x) is less straightforward. In particular, in the ballistic limit (l ee = ∞) the probe attached to the edge of the sample does not register particles grazing along the edge. This suppresses the space charge effect, however this suppression is not a universal phenomenon, and should be viewed as an approximation.
Numerical modeling. To model the experimental geometry of Fig. 1a, we analyze the flow induced in a strip of width W, 0 < y < W by a point source on its edge at the point (0, 0) and a drain at x = −∞, see Fig. 4. Such a flow can be represented as a superposition of a symmetric flow emitted by the source with a uniform flow directed to the drain electrode. Both flows can be analysed numerically via the approach described below.
First, we pass to the Fourier representation with respect to the coordinate x along the strip, and discretize the transverse coordinate: y n = nh, where h = W/N y is the step size, n = 0,…, N y − 1. We also discretize the momentum direction as θ i = π(i + 1/2)/N θ , i = 0,…2N θ − 1. Hence the distribution function becomes a function of the wavevector k and two discrete coordinates, We employ the following finite-difference representation of the kinetic equation For numerical stability, the scheme is made "upwind": the form on the first line should be applied for upward-going particles (0 < θ i < π), and the form on the second line describes particles propagating downwards. Due to the choice of the exponential factors, the exact solution of (12) in the collisionless limit (γ ee = 0), χðk; y; θÞ / expðÀiky cotanθÞ, satisfies the discretized equation.
Thus, discretization of the advection term in Boltzmann equation, (v∇)δf, links the values of the distribution function at nearby sites. The discretized form of the collision integral (13) mixes propagation angles within the same site. Therefore, the above finite-difference system, together with the boundary conditions (14) can be recast into the well-known three-diagonal form in which only blocks on three adjacent sites n and n ± 1 are coupled: A n;ij ðkÞχ nÀ1;j ðkÞ þ B n;ij χ n;j ðkÞ þ C n;ij ðkÞχ nþ1;j ¼ b n;i ðkÞ: ð19Þ Here A n,ij (k), B n,ij (k) and C n,ij (k) are matrix operators acting on the angular index i describing propagation of particles and scattering between different momentum directions. The right-hand side b n,i (k) describes external sources of particles. Such a system can be efficiently solved via the standard three-diagonal matrix algorithm 44 .
The point source was represented as a source term in the Eq. (14), with Fourier image I(k) = 1. The uniform Poiseuille-like flow can be obtained by analyzing the k = 0 limit of the three-diagonal system (19), in which the flow is dragged by an external bias field. The bias field is incorporated into the Eq. (12) via the term −eEcosθ. The value of the bias field E is then obtained by normalizing the solution to the total current of 1/2.
To make sure that the details of the boundary layer near the edges are simulated properly, we have chosen a rather fine grid, N y = 5000. The propagation angles were discretized with N θ = 50, which corresponds to 3.6°step in θ. The particle distribution χ n,i (k) was calculated for |k|W < 50, which gives a satisfactory approximation to the distances of interest, 0.1W < x < W. The probe signal was then calculated as the particle flux (16), giving the results shown in Figs. 1c, 2b, 3b.

Data availability
The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.