pt

thanks: Contributed equally to this paper
jorgestr@ucm.es
thanks: Contributed equally to this paper
jorgestr@ucm.es

Superballistic conduction in hydrodynamic antidot graphene superlattices

Jorge Estrada-Álvarez GISC, Departamento de Física de Materiales, Universidad Complutense, 28040 Madrid, Spain.    Juan Salvador-Sánchez Nanotechnology Group, USAL–Nanolab, University of Salamanca
Plaza de la Merced, Edificio Trilingüe, 37008, Salamanca, Spain.
   Ana Pérez-Rodríguez Nanotechnology Group, USAL–Nanolab, University of Salamanca
Plaza de la Merced, Edificio Trilingüe, 37008, Salamanca, Spain.
   Carlos Sánchez-Sánchez Nanotechnology Group, USAL–Nanolab, University of Salamanca
Plaza de la Merced, Edificio Trilingüe, 37008, Salamanca, Spain.
   Vito Clericò Nanotechnology Group, USAL–Nanolab, University of Salamanca
Plaza de la Merced, Edificio Trilingüe, 37008, Salamanca, Spain.
   Daniel Vaquero Nanotechnology Group, USAL–Nanolab, University of Salamanca
Plaza de la Merced, Edificio Trilingüe, 37008, Salamanca, Spain.
   Kenji Watanabe Research Center for Functional Materials, National Institute for Materials Science
1-1 Namiki, Tsukuba 305-0044, Japan.
   Takashi Taniguchi International Center for Materials Nanoarchitectonics, National Institute for Materials Science
1-1 Namiki, Tsukuba 305-0044, Japan.
   Enrique Diez Nanotechnology Group, USAL–Nanolab, University of Salamanca
Plaza de la Merced, Edificio Trilingüe, 37008, Salamanca, Spain.
   Francisco Domínguez-Adame GISC, Departamento de Física de Materiales, Universidad Complutense, 28040 Madrid, Spain.    Mario Amado Nanotechnology Group, USAL–Nanolab, University of Salamanca
Plaza de la Merced, Edificio Trilingüe, 37008, Salamanca, Spain.
   Elena Díaz GISC, Departamento de Física de Materiales, Universidad Complutense, 28040 Madrid, Spain.
(July 5, 2024)
Abstract

Viscous electron flow exhibits exotic signatures such as superballistic conduction. Bending the geometry of the device is a must to observe hydrodynamic effects. To this end, we build three antidot graphene superlattices with different hole diameters. We measure their electrical properties at various temperatures and under the effect of a perpendicular magnetic field. We find an enhanced superballistic effect, suggesting the effectiveness of the geometry at bending the electron flow. In addition, superballistic conduction behaves non-monotonically with the magnetic field, which is related with the ballistic-hydrodynamic transition. We also analyze the device resistance as a function of the size of the antidot superlattice to find characteristic scaling laws describing the different transport regimes. We prove that the antidot superlattice is a convenient geometry for realizing hydrodynamic flow, and the experiment provides valuable explanations for the technologically relevant effects of superballistic conduction and scaling laws.

I Introduction

Collisions against impurities and phonons dominate electron-electron collisions in ordinary metals in most cases. In 1963, however, Gurzhi considered the opposite situation and claimed that a decrease in the electrical resistance with increasing temperature might appear in ultra-clean metals at moderate temperatures  [1]. This author attributed the phenomenon to the realization of a hydrodynamic regime of charge transport, where highly-correlated electrons behave collectively in a similar way to molecules in conventional viscous fluids [2]. This effect, also known as superballistic conduction, constitutes one of the archetypal hydrodynamic signatures [3, 4, 5, 6, 7, 8, 9, 10]. Moreover, since the collective motion of electrons leads to a resistance below the ballistic limit, superballistic conduction is a convenient property for low-power consumption devices [11]. In this regard, a reduction of the resistance up to 16%percent1616\%16 % has been reported in point contacts [9, 12, 13, 14, 15] and, up to 4%percent44\%4 %, in crenellated channels [16].

Another well-recognized exotic hydrodynamic signature is the experimental evidence of a Poiseuille flow [8]. In traditional fluids, the Poiseuille law is one of the fundamental principles that affect any fluid network. Its most immediate consequence is a resistance that scales as R1/ds4proportional-to𝑅1superscriptsubscript𝑑𝑠4R\propto 1/d_{s}^{4}italic_R ∝ 1 / italic_d start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT where dssubscript𝑑𝑠d_{s}italic_d start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT is the diameter of a single pipe carrying the fluid [17, 18, 19, 20]. If, instead of a single pipe, the space is filled up with several pipes of diameter d𝑑ditalic_d the Poiseuille law can be demonstrated to read R1/d2proportional-to𝑅1superscript𝑑2R\propto 1/d^{2}italic_R ∝ 1 / italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT [18]. Notice that this situation is the equivalent of an antidot superlattice in a two dimensional (2D) system. However, to the best of our knowledge, the equivalent of such scaling law for electrons is not understood in detail.

In the last decade, the development of 2D materials has multiplied the experimental realization of electron hydrodynamics. In this transport regime, electrons travel long distances, larger than the size of typical devices, before scattering against impurities or phonons. This is the case of ultra-pure PdCoO2 [5, 21], Weyl semimetals [22], (Al,Ga)As heterostructures [16] or graphene [8, 9, 6, 7], whose electronic properties dramatically improve when it is encapsulated between two layers of hexagonal boron nitride (hBN) [23]. However, direct visualization of the hydrodynamic flow often requires a complex microscopy setup [8, 10, 24, 22], making it difficult to establish an ubiquitous criterion to establish the occurrence of hydrodynamic transport [8, 7]. Due to the potential applications in electronic design and the need to easily explore new materials [4], it is desirable to look for novel platforms for viscous electron flow.

In this work, we demonstrate that electron trajectories become naturally bent after geometrically engineering the device. This boosts the hydrodynamic signatures, which highly depend on the device size and the non-uniformity of the electron flow. In particular, in this work, we build different antidot graphene superlattices to study superballistic conduction as an indicator of hydrodynamic flow. The experimental setup also allows us to study scaling laws that resemble the traditional fluid Poiseuille law. We also demonstrate a non-monotonic electrical response as a function of the magnetic field, which requires the development of an alternative theoretical understanding. Our experimental evidences are perfectly supported by detailed theoretical simulations that offer a better insight into the fundamentals of superballistic conduction.

II Geometrically-engineered devices for hydrodynamics

Despite the great quality of graphene heterostructures, initial experiments did not find any sign of the expected hydrodynamics regime. Bandurin et al. took an important step forward [6] by noting that the kinematic viscosity enters the Navier-Stokes equation as a coefficient in front of the second spatial derivative of velocity. Therefore, in order to maximize the viscous flow, a current flow as inhomogeneous as possible is needed. Together with the geometry, the edge scattering also needs to be controlled. Recently, the control of the hydrodynamic flow in a 2D electron gas of GaAs quantum wells with smooth boundaries was demonstrated in crenellated channels [16], with reductions of the resistance up to 4%percent44\%4 %.

In the present work, we design an optimized antidot geometry (see Appendix A) to bend the electron flow and uncover the collective hydrodynamic response of the electrons on a fully encapsulated graphene heterostructure. We perform a cryo-etching process [25] to favor smooth edges that ensure almost specular reflection at the edges and use a graphite bottom gate which enhances the charge mobility. In addition the latter screens undesired spurious effects arising from charged defects present in the underlying doped silicon substrate [26]. This engineered geometry helps us to disentangle viscous and ballistic effects.

Refer to caption
Figure 1: Graphene antidot superlattice. (a) Optical image and schematics of the electronic device. (b) Electron micrograph of a hBN flake showing the same antidot geometry lithographed onto the final device. Three regions can be found where the antidot diameter reaches d=100, 200𝑑100200d=100,\,200italic_d = 100 , 200, and 300nm300nm300\,\rm nm300 roman_nm. The center-to-center antidot distance is 2d2𝑑2d2 italic_d and they were arranged into a square lattice. (c) Close electron micrograph for d=100𝑑100d=100\,italic_d = 100nm antidots displaying smooth edges due to the cryo-etching process. (d) Simulation of the Boltzmann transport equation, where colors account for the potential and streamlines are average electron bent trajectories.

The heterostructure was fabricated by means of the standard mechanical exfoliation on pristine crystals of hBN and graphite, subsequently deposited onto silicon oxide wafer with top p-doped silicon oxide. Then, the stack was finally shaped using e-beam lithography into a typical 10-terminal Hall bar with a longitudinal contact-to-contact distance L=4μ𝐿4𝜇L=4\,\muitalic_L = 4 italic_μm, shown in Fig. 1(a). Once the final Hall bar was produced a last step of electron beam lithography + cryo etching [27] process was performed in order to define the antidot patterns as shown in Fig. 1(b) and (c). Three separated regions are shown, consisting of three antidot meshes of different diameters, namely d=100, 200𝑑100200d=100,\,200italic_d = 100 , 200, and 300nm300nm300\,\rm nm300 roman_nm. In all three regions, the antidots appear in a square lattice with a center-to-center distance of 2d2𝑑2d2 italic_d, defining our final periodic geometry. Fig. 1(d) displays the bent electron trajectories in this geometry. Different dose array tries were performed in order to define pristine structures with an enhanced definition of the antidots (see Appendix B for more details).

III Transport regimes

Viscous electron flow in 2D materials is nothing but a collective motion of conduction electrons that is expected to reduce the device resistance according to Gurzhi’s prediction [1, 2]. In particular, superballistic conduction involves the transition from a ballistic non-collective regime to a collective hydrodynamic one. Consequently, we must use a kinetic model that covers both regimes: the semiclassical Boltzmann transport equation  [28, 29] (see Appendix C for a detailed description). Certainly, it is more complex than a hydrodynamic model, but it is the only feasible way to rigorously explore all possible regimes of electronic transport, including the intermediate ones. In our case of interest the Boltzmann equation describing the distribution g(𝒓,θ)𝑔𝒓𝜃g(\bm{r},\theta)italic_g ( bold_italic_r , italic_θ ), which counts the electrons at position 𝒓𝒓\bm{r}bold_italic_r moving in the direction of θ𝜃\thetaitalic_θ, obeys

𝒖^𝒌𝒓(geV(𝒓)mvF)+θglBsubscript^𝒖𝒌subscript𝒓𝑔𝑒𝑉𝒓𝑚subscriptv𝐹subscript𝜃𝑔subscript𝑙𝐵\displaystyle\hat{\bm{u}}_{\bm{k}}\cdot\nabla_{{\bm{r}}}\left(g-\frac{eV(\bm{r% })}{m{\rm v}_{F}}\right)+\frac{\partial_{\theta}g}{l_{B}}over^ start_ARG bold_italic_u end_ARG start_POSTSUBSCRIPT bold_italic_k end_POSTSUBSCRIPT ⋅ ∇ start_POSTSUBSCRIPT bold_italic_r end_POSTSUBSCRIPT ( italic_g - divide start_ARG italic_e italic_V ( bold_italic_r ) end_ARG start_ARG italic_m roman_v start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT end_ARG ) + divide start_ARG ∂ start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT italic_g end_ARG start_ARG italic_l start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT end_ARG
+gle+gevengeeevenleeeven+goddgeeoddleeodd𝑔subscript𝑙𝑒superscript𝑔evensubscriptsuperscript𝑔even𝑒𝑒superscriptsubscript𝑙𝑒𝑒evensuperscript𝑔oddsubscriptsuperscript𝑔odd𝑒𝑒superscriptsubscript𝑙𝑒𝑒odd\displaystyle+\frac{g}{l_{e}}+\frac{g^{\mathrm{even}}-g^{\mathrm{even}}_{ee}}{% l_{ee}^{\mathrm{even}}}+\frac{g^{\mathrm{odd}}-g^{\mathrm{odd}}_{ee}}{l_{ee}^{% \mathrm{odd}}}+ divide start_ARG italic_g end_ARG start_ARG italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT end_ARG + divide start_ARG italic_g start_POSTSUPERSCRIPT roman_even end_POSTSUPERSCRIPT - italic_g start_POSTSUPERSCRIPT roman_even end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT end_ARG start_ARG italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_even end_POSTSUPERSCRIPT end_ARG + divide start_ARG italic_g start_POSTSUPERSCRIPT roman_odd end_POSTSUPERSCRIPT - italic_g start_POSTSUPERSCRIPT roman_odd end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT end_ARG start_ARG italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_odd end_POSTSUPERSCRIPT end_ARG =0.absent0\displaystyle=0\ .= 0 . (1)

where 𝒖^𝒌=(cosθ,sinθ)subscript^𝒖𝒌𝜃𝜃\hat{\bm{u}}_{\bm{k}}=(\cos\theta,\sin\theta)over^ start_ARG bold_italic_u end_ARG start_POSTSUBSCRIPT bold_italic_k end_POSTSUBSCRIPT = ( roman_cos italic_θ , roman_sin italic_θ ), vF106ms1subscriptv𝐹superscript106superscriptms1{\rm v}_{F}\approx 10^{6}\,\rm ms^{-1}roman_v start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT ≈ 10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT roman_ms start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT is the Fermi velocity, and m=kF/vF𝑚Planck-constant-over-2-pisubscript𝑘𝐹subscriptv𝐹m=\hbar k_{F}/{\rm v}_{F}italic_m = roman_ℏ italic_k start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT / roman_v start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT the cyclotron mass, being kF=π|n|subscript𝑘𝐹𝜋𝑛k_{F}=\sqrt{\pi|n|}italic_k start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT = square-root start_ARG italic_π | italic_n | end_ARG the Fermi wavenumber in graphene for a density of carriers n𝑛nitalic_n. The electrons propagate under the effect of a potential V(𝒓)𝑉𝒓V(\bm{r})italic_V ( bold_italic_r ) and a perpendicular magnetic field B𝐵Bitalic_B with a cyclotron radius lB=mvF/eBsubscript𝑙𝐵𝑚subscriptv𝐹𝑒𝐵l_{B}=m{\rm v}_{F}/eBitalic_l start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT = italic_m roman_v start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT / italic_e italic_B. In order to reproduce the experimental scenario Eq. 1 is solved numerically, using a finite element method (Appendix D) to compute the drift velocity as 𝒖=(1/π)02π𝒖^𝒌g(𝒓,θ)dθ𝒖1𝜋superscriptsubscript02𝜋subscript^𝒖𝒌𝑔𝒓𝜃differential-d𝜃\bm{u}=(1/\pi)\int_{0}^{2\pi}\hat{\bm{u}}_{\bm{k}}g(\bm{r},\theta)\,{\rm d}\thetabold_italic_u = ( 1 / italic_π ) ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 italic_π end_POSTSUPERSCRIPT over^ start_ARG bold_italic_u end_ARG start_POSTSUBSCRIPT bold_italic_k end_POSTSUBSCRIPT italic_g ( bold_italic_r , italic_θ ) roman_d italic_θ. Finally, the associated current density is integrated to find the current and the longitudinal resistance R𝑅Ritalic_R, which is the electrical property of interest in the experiments.

Equation (1) accounts for different mechanisms of electron collisions that largely affect the electrical response of any material. Collisions with crystal imperfections or lattice vibrations alter the total momentum of the electrons, characterized by the mean free path lesubscript𝑙𝑒l_{e}italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT. We estimate le=700nmsubscript𝑙𝑒700nml_{e}=700\,\rm nmitalic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT = 700 roman_nm for a gate voltage, referred to the Dirac peak, of Vg=3.9Vsubscript𝑉𝑔3.9VV_{g}=3.9\,\rm Vitalic_V start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT = 3.9 roman_V, which is a reasonable value after the nanostructuring process of the device [30]. Moreover, we also take into account scattering due to phonons [31, 32] at increasing temperatures (Appendix E). Contrary to impurities and phonons, electron-electron collisions conserve the total momentum, with a mean free path leesubscript𝑙𝑒𝑒l_{ee}italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT that can be computed for graphene [9, 33]. In this work, we go beyond the conventional Callaway’s ansatz [34] and define two different relaxation rates for the even (leeevensuperscriptsubscript𝑙𝑒𝑒evenl_{ee}^{\mathrm{even}}italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_even end_POSTSUPERSCRIPT) and odd-parity (leeoddsuperscriptsubscript𝑙𝑒𝑒oddl_{ee}^{\mathrm{odd}}italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_odd end_POSTSUPERSCRIPT) modes in the expansion of g𝑔gitalic_g. Indeed recent theoretical investigation [35, 36] proves that leeeven=leesuperscriptsubscript𝑙𝑒𝑒evensubscript𝑙𝑒𝑒l_{ee}^{\mathrm{even}}=l_{ee}italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_even end_POSTSUPERSCRIPT = italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT while leeoddleemuch-greater-thansuperscriptsubscript𝑙𝑒𝑒oddsubscript𝑙𝑒𝑒l_{ee}^{\mathrm{odd}}\gg l_{ee}italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_odd end_POSTSUPERSCRIPT ≫ italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT, resulting in the so-called tomographic description, which, for the sake of higher accuracy, we follow in this article. Last, there are also collisions against the edges, which are described with the appropriate boundary condition [37]. In our case, the cryo-etching technique produces a smooth boundary, so we assume a perfect slip boundary condition, where electrons have specular reflection at the edges.

The most commonly accepted route towards viscous electron flow is favoring electron-electron collisions to achieve the condition lee<Wsubscript𝑙𝑒𝑒𝑊l_{ee}<Witalic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT < italic_W or, more strictly by adding the requirement lee<lesubscript𝑙𝑒𝑒subscript𝑙𝑒l_{ee}<l_{e}italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT < italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT too, where W𝑊Witalic_W is the typical size of the device. However, recent works demonstrated these traditional requirements are too restrictive and there exist alternative pathways to achieve the hydrodynamic regime [28].

Although the Boltzmann equation is essential to describe the ballistic-hydrodynamic transition, a hydrodynamic Navier-Stokes equation is the common reference to viscous electron flow [6, 3]. The hydrodynamic equation is expected to provide accurate predictions in the hydrodynamic regime. Conversely, it gives wrong predictions in the ballistic one, so its error can be used to classify the nature of transport  [28] and in particular to distinguish between the ballistic and hydrodynamic transport. Regarding the current study, we want to highlight that the superballistic effect will be theoretically and experimentally proven when a ballistic-hydrodynamic transition occurs.

IV Results and discussion

Refer to caption
Figure 2: Enhanced superballistic conduction. Experimental resistance in the three regions of the antidot lattice with d=100, 200𝑑100200d=100,\,200italic_d = 100 , 200, and 300nm300nm300\,\rm nm300 roman_nm as a function of temperature for (a) electrons (Vg=3.9Vsubscript𝑉𝑔3.9VV_{g}=3.9\,\rm Vitalic_V start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT = 3.9 roman_V) and (b) holes (Vg=3.9Vsubscript𝑉𝑔3.9VV_{g}=-3.9\,\rm Vitalic_V start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT = - 3.9 roman_V), Vgsubscript𝑉𝑔V_{g}italic_V start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT defined respect to the Dirac point, and (c) theoretical resistance calculated by Boltzmann equation simulations. (d) Experimental magnitude of dR/dTd𝑅d𝑇\text{d}R/\text{d}Td italic_R / d italic_T with increasing temperature and various Vgsubscript𝑉𝑔V_{g}italic_V start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT for the d=200nm𝑑200nmd=200\,\rm nmitalic_d = 200 roman_nm region. (e) Gurzhi ratio 𝒢𝒢\mathcal{G}caligraphic_G as a function of Vgsubscript𝑉𝑔V_{g}italic_V start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT for the three antidot regions.

IV.1 Enhanced superballistic effect

First, we investigate the superballistic effect of the antidot superlattices. Experimental results in Fig. 2(a) and (b) clearly exhibit a decrease in the electrical resistance below its ballistic limit at low temperature in the three antidot regions. Thus, our results demonstrate the existence of the superballistic effect and constitute a signature of hydrodynamic transport. This supports the effectiveness of the antidot superlattice at bending the electron flow. Figure 2(c) shows the predictions of the Boltzmann equation for the same antidot geometries, showing a very good agreement with the experiment. In our simulations, we set the same value le=700subscript𝑙𝑒700l_{e}=700\,italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT = 700nm for the three regions for an easier comparison. Moreover, we have discussed (Appendix F) the influence of possible charge inhomogeneity effects to explain the better agreement with the experiment for hole conduction [see results for d=100nm𝑑100nmd=100\,\rm nmitalic_d = 100 roman_nm in Figs. 2(b) and (c)]. Also, the latter shows an enhanced superballistic conduction in the region of d=100nm𝑑100nmd=100\,\rm nmitalic_d = 100 roman_nm that may survive near room temperature. This result, different from the one observed in the regions of d=200nm𝑑200nmd=200\,\rm nmitalic_d = 200 roman_nm and 300nm300nm300\,\rm nm300 roman_nm, had not been reported in point contacts [9]. Nevertheless, it is accurately predicted by the Boltzmann equation. Indeed, both the bending of the electron flow and the smaller d/le𝑑subscript𝑙𝑒d/l_{e}italic_d / italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT restrain the detrimental impact of phonon scattering. Similar to the formation of current whirlpools recently observed in graphene [24], our work is consistent with hydrodynamic effects at room temperature as well. This brings in a paradigm where the technological advantages of reduced electrical resistance can be further exploited. For completeness, Fig. 2(d) monitors the superballistic effect by means of the dR/dT𝑑𝑅𝑑𝑇dR/dTitalic_d italic_R / italic_d italic_T measured in the region of d=200nm𝑑200nmd=200\,\rm nmitalic_d = 200 roman_nm when the corrected gate voltage Vgsubscript𝑉𝑔V_{g}italic_V start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT and T𝑇Titalic_T are varied. Notice that hereafter we will refer to Vgsubscript𝑉𝑔V_{g}italic_V start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT as the gate potential defined respect to the Dirac point. Last, Fig. 2(e) shows the magnitude of the superballistic effect estimated as the Gurzhi ratio 𝒢=1R(50K)/R(1.5K)𝒢1𝑅50K𝑅1.5K\mathcal{G}=1-R(50{\,\rm K})/R(1.5{\,\rm K})caligraphic_G = 1 - italic_R ( 50 roman_K ) / italic_R ( 1.5 roman_K ). This quantity allows us to distinguish the true superballistic effect from other physical scenarios occurring at higher temperatures near the charge neutrality point (within the dashed lines in Fig. 2(d) and (e)), as discussed in Appendix G. Reached ratios of 𝒢20%greater-than-or-equivalent-to𝒢percent20\mathcal{G}\gtrsim 20\%caligraphic_G ≳ 20 % improve those previously reported in point contacts [9]. Once again, this enlarged 𝒢𝒢\mathcal{G}caligraphic_G results as a consequence of the enhanced bending of the electron trajectories by the antidot superlattice. Numerical calculations offer a better insight into the fundamentals of superballistic conduction. Contrary to what was expected, the condition leedlemuch-less-thansubscript𝑙𝑒𝑒𝑑much-less-thansubscript𝑙𝑒l_{ee}\ll d\ll l_{e}italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT ≪ italic_d ≪ italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT is not essential for the Gurzhi effect to occur [1, 2]. Notice that leesubscript𝑙𝑒𝑒l_{ee}italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT is indeed very large in the low temperature limit. For example, in our experiment, we demonstrate superballistic conduction even at T=50K𝑇50KT=50\,\rm Kitalic_T = 50 roman_K when lee1400nmsubscript𝑙𝑒𝑒1400nml_{ee}\approx 1400\,\rm nmitalic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT ≈ 1400 roman_nm while le700nmsubscript𝑙𝑒700nml_{e}\approx 700\,\rm nmitalic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT ≈ 700 roman_nm. Furthermore, such a result is supported by our simulations that include the particular geometrical details beyond the approximated anti-Matthiessen rule [9, 15]. The experimental evidence violates the criteria leedmuch-less-thansubscript𝑙𝑒𝑒𝑑l_{ee}\ll ditalic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT ≪ italic_d and leelemuch-less-thansubscript𝑙𝑒𝑒subscript𝑙𝑒l_{ee}\ll l_{e}italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT ≪ italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT, so another reasoning must be developed. We conclude that not only the values of lesubscript𝑙𝑒l_{e}italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT and leesubscript𝑙𝑒𝑒l_{ee}italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT determine the collectivization of the electron flow, but mostly their decrease rate with temperature since dR/dT=(R/le)(le/T)+(R/lee)(lee/T)d𝑅d𝑇𝑅subscript𝑙𝑒subscript𝑙𝑒𝑇𝑅subscript𝑙𝑒𝑒subscript𝑙𝑒𝑒𝑇\text{d}R/\text{d}T=(\partial R/\partial l_{e})(\partial l_{e}/\partial T)+(% \partial R/\partial l_{ee})(\partial l_{ee}/\partial T)d italic_R / d italic_T = ( ∂ italic_R / ∂ italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT ) ( ∂ italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT / ∂ italic_T ) + ( ∂ italic_R / ∂ italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT ) ( ∂ italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT / ∂ italic_T ). Thus, the existence and the magnitude of the superballistic effect mainly depends of such rates.

In our theoretical study, we also paid attention to the relevance of the boundary conditions to reproduce the experimental results [37] (see Appendix H). Accordingly, we determine that an almost perfect slip condition typical of a specular boundary is compatible with our measurements. This is consistent with the high crystalline quality attained with the cryo-etching technique used in the sample preparation [25] since smooth edges favor the momentum conservation at the boundaries. Similarly to previous works [16], we can conclude that the demonstrated superballistic effect is almost universal since it relies on the geometry of the device and not on the particular considered edge scattering mechanism.

We also check the superballistic conduction in another device where the experimental measurements agree with these results (see Appendix I). More importantly, in the additional device, we demonstrate that there is no decrease in the resistance for a region where no antidots were fabricated. This further confirms that geometrical-engineering is responsible for the observed superballistic conduction.

Refer to caption
Figure 3: Intermittent superballistic conduction. (a) Longitudinal resistance for the three regions at Vg=3.9Vsubscript𝑉𝑔3.9VV_{g}=3.9\,\rm Vitalic_V start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT = 3.9 roman_V and T=1.5K𝑇1.5KT=1.5\,\rm Kitalic_T = 1.5 roman_K shows weak localization (WL), commensurability effect (C), the hydrodynamic region (H), and oscillations associated with the quantum Hall effect (QHE) . (b) Superballistic effect for several magnetic fields. (c) Experimental measurements and (d) Boltzmann equation simulations of the resistance as a function of the magnetic field, normalized to the commmensurability fields BCsubscript𝐵𝐶B_{C}italic_B start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT of the peak. (e) Gurzhi ratio 𝒢𝒢\mathcal{G}caligraphic_G, symbols are experimental observations and the lines are simulated values. (f) Qualitative explanation of the intermittent effect, where W𝑊Witalic_W is the typical width of a uniform device.

IV.2 Intermittent superballistic effect

Here we explore the superballistic effect in more detail by considering the electrical response of the antidot lattice in the presence of a magnetic field. First, let us analyze the low-temperature resistance as a function of the applied magnetic field, as shown in Fig. 3(a). Evidence of two quantum effects is shown in our measurements: i) the weak localization peak (WL) at B20mTless-than-or-similar-to𝐵20mTB\lesssim 20\,\rm mTitalic_B ≲ 20 roman_mT, due to quantum interference [38], and enhanced by intervalley scattering against the superlattice [39] and ii) the quantum Hall effect (QHE), whose peaks flatten for smaller values of d𝑑ditalic_d, suggesting a prominent role of the antidot geometry. More relevant for the matter of interest are the peaks at the particular field BC1.05πn/edsubscript𝐵𝐶1.05Planck-constant-over-2-pi𝜋𝑛𝑒𝑑B_{C}\approx 1.05\,\hbar\sqrt{\pi n}/editalic_B start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT ≈ 1.05 roman_ℏ square-root start_ARG italic_π italic_n end_ARG / italic_e italic_d related to the commensurability effect (C), occurring when the cyclotron radius is commensurate to the antidot lattice size [40, 30, 41]. Consequently, these peaks are shifted to higher magnetic fields for increasing d𝑑ditalic_d as indicated by the dotted line in Fig. 3(a). For d=200𝑑200d=200italic_d = 200 and 300nm300nm300\,\rm nm300 roman_nm, the hydrodynamic negative magnetoresistance (H) that follows this peak, due to the collectivization of electron flow, is also visible. The region of commensurability is zoomed out in Fig. 3(c) and numerically represented in Fig. 3(d) where the applied magnetic field is expressed in units of BCsubscript𝐵𝐶B_{C}italic_B start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT.

Regarding the superballistic conduction appearing when the temperature is increased, Fig. 3(b) shows that the effect arises for some magnetic fields and disappears for others, in an intermittent pattern. This is clearly revealed in Fig. 3(e) where the Gurzhi ratio 𝒢𝒢\mathcal{G}caligraphic_G is represented as a function of B/BC𝐵subscript𝐵𝐶B/B_{C}italic_B / italic_B start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT. Notice the agreement between the experimental measurements and the Boltzmann equation simulations in Figs. 3(c)-(e). Since the considered sizes d𝑑ditalic_d are much larger than the Fermi wavelength in these devices, the impact of quantum effects is expected to be reduced [40]. Also, the role of the tomographic approach in describing electron transport and its comparison against the hydrodynamic approach is explored in Appendix J, showing that experimental data supports the tomographic approach. The experimental results demonstrate the existence of superballistic transport not only at zero magnetic field as discussed in Sec. IV.1, but also when the magnetic field reaches the commensurability condition. Both cases are consistent with physical scenarios that support a ballistic-hydrodynamic transition as schematically summarized in Fig. 3(f). There we show [28] how the nature of the transport regime is established by way of the deviation of the hydrodynamic model results from the Boltzmann equation. In particular, Fig. 3(f) shows the transition between the ballistic (B) and the hydrodynamic (H) regime when the temperature or the magnitude of a perpendicular field is increased in a uniform graphene channel. Regarding the current study, we want to highlight that the superballistic effect is theoretically and experimentally proven when a ballistic-hydrodynamic transition occurs (notice blue solid lines in Fig. 3(f)). In summary, the superballistic conduction can be used to classify the ballistic and hydrodynamic transport regimes. The conventional formalism based on the Anti-Matthiessen rule [15, 9] does not reproduce the intermittent effect with magnetic field, not even qualitatively, and it had to be replaced by the Boltzmann equation description. Moreover, the classification of the hydrodynamic transition with a magnetic field, which only shows one hydrodynamic region [42], had to be replaced by another one with two hydrodynamic regions, in better agreement with Ref. [28].

IV.3 Quasi-Poiseuille law

Refer to caption
Figure 4: (a) Ratio between the measured resistances R1subscript𝑅1R_{1}italic_R start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT and R2subscript𝑅2R_{2}italic_R start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT in the d=100nm𝑑100nmd=100\,\rm nmitalic_d = 100 roman_nm and 200nm200nm200\,\rm nm200 roman_nm regions, respectively. The gate voltage is referred to the charge neutrality point. Notice that R1/R2>2α>1subscript𝑅1subscript𝑅22𝛼1R_{1}/R_{2}>2\Rightarrow\alpha>1italic_R start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / italic_R start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT > 2 ⇒ italic_α > 1. (b) Log-log plot for the resistance R𝑅Ritalic_R versus the antidot size d𝑑ditalic_d at Vg=3.9Vsubscript𝑉𝑔3.9VV_{g}=3.9\,\rm Vitalic_V start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT = 3.9 roman_V. Symbols are experimental results and lines are simulations of the Boltzmann equation. Notice the slope of the curves accounts for |α|𝛼|\alpha|| italic_α |, see the guidelines for α=1𝛼1\alpha=1italic_α = 1 and 2222. Different curves have been artificially shifted for a proper visualization. Panels (c)-(e) represent the potential color map inside the sample with the streamlines for the electron fluid and the velocity profile across a transversal section in different transport conditions. (c) Ballistic transport regime le=lee=5dsubscript𝑙𝑒subscript𝑙𝑒𝑒5𝑑l_{e}=l_{ee}=5ditalic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT = italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT = 5 italic_d where electrons move in a channel dodging holes. (d) Viscous transport le=5d,lee=0.2dformulae-sequencesubscript𝑙𝑒5𝑑subscript𝑙𝑒𝑒0.2𝑑l_{e}=5d,l_{ee}=0.2ditalic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT = 5 italic_d , italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT = 0.2 italic_d where the profiles are smoother. (e) Diffusive transport le=0.2d,lee=5dformulae-sequencesubscript𝑙𝑒0.2𝑑subscript𝑙𝑒𝑒5𝑑l_{e}=0.2d,l_{ee}=5ditalic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT = 0.2 italic_d , italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT = 5 italic_d where the antidot geometry loses its relevance and the electrons move through the whole sample.

The electronic equivalent of the Poiseuille law [20] provides relevant information about the regimes of transport. Let us assume that the resistance of the device scales as R1/dαproportional-to𝑅1superscript𝑑𝛼R\propto 1/d^{\alpha}italic_R ∝ 1 / italic_d start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT such that we could define αlogR/logd𝛼𝑅𝑑\alpha\equiv-{\partial\log R}/{\partial\log d}italic_α ≡ - ∂ roman_log italic_R / ∂ roman_log italic_d. For instance, α=0𝛼0\alpha=0italic_α = 0 in a purely diffusive regime (ledmuch-less-thansubscript𝑙𝑒𝑑l_{e}\ll ditalic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT ≪ italic_d), and α=1𝛼1\alpha=1italic_α = 1 is the value predicted by the Landauer formula [43, 9]. The physics is more complex when a hydrodynamic description applies in electron systems. Indeed, in Appendix K we find

α21+vFd2βνle,similar-to𝛼21subscriptv𝐹superscript𝑑2𝛽𝜈subscript𝑙𝑒\alpha\sim\frac{2}{1+\frac{\mathrm{v}_{F}d^{2}}{\beta\nu l_{e}}}\ ,italic_α ∼ divide start_ARG 2 end_ARG start_ARG 1 + divide start_ARG roman_v start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG italic_β italic_ν italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT end_ARG end_ARG , (2)

where β4.4𝛽4.4\beta\approx 4.4italic_β ≈ 4.4 is the bending ratio for the antidot geometry and ν𝜈\nuitalic_ν is the viscosity. The expression assumes specular edges, which is a reasonable hypothesis with the cryo-etching technique [25]. Notice that the hydrodynamic description supports all values 0<α<20𝛼20<\alpha<20 < italic_α < 2. This is different from the Poiseuille law in conventional fluids, leading to α=2𝛼2\alpha=2italic_α = 2 [18]. Indeed, scattering against impurities and phonons strongly affects the values of α𝛼\alphaitalic_α, and any finite value of lesubscript𝑙𝑒l_{e}italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT results in α<2𝛼2\alpha<2italic_α < 2. Therefore, it is not α=2𝛼2\alpha=2italic_α = 2, but rather 1<α<21𝛼21<\alpha<21 < italic_α < 2 that is assumed as a landmark in electron systems [20]. Nonetheless, the maximum values of α𝛼\alphaitalic_α occur for shorter d𝑑ditalic_d, longer lesubscript𝑙𝑒l_{e}italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT, and longer ν𝜈\nuitalic_ν, namely, at low temperatures. Because of the inevitable collisions against impurities and phonons, reducing the viscosity ν𝜈\nuitalic_ν to enter the hydrodynamic regime does not increase α𝛼\alphaitalic_α, but, on the contrary, diminishes it. Hence, we expect the largest values of α𝛼\alphaitalic_α near the ballistic regime. Indeed, this explains the observations in Fig. 4(a), where we plot the ratio between the resistances R1/R22αsimilar-tosubscript𝑅1subscript𝑅2superscript2𝛼R_{1}/R_{2}\sim 2^{\alpha}italic_R start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / italic_R start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ∼ 2 start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT in the d=100𝑑100d=100italic_d = 100 and 200nm200nm200\,\rm nm200 roman_nm regions. The maximum occurs at low temperatures and intermediate gate voltages, where the mean free paths of graphene are larger [6].

We study the scaling laws with a Log-Log plot in Fig. 4(b) and compare it with simulations of the Boltzmann equation. Yet simulations improve the accuracy, they show the same qualitative behavior as Eq. (2). Namely, the ballistic regime, whose flow profiles for a paradigmatic set of parameters is shown in Fig. 4(c), gives the largest values of α𝛼\alphaitalic_α. This regime is approached for low temperatures and shorter d𝑑ditalic_d. Increasing the temperature favors a transition to a more hydrodynamic regime, for which an archetypal example is shown in Fig. 4(d). Eventually, a diffusive regime, akin to the one in Fig. 4(e), is attained, and α0𝛼0\alpha\to 0italic_α → 0. Last, a magnetic field decreases the value of α𝛼\alphaitalic_α, as it also induces a transition from ballistic to hydrodynamic [28]. In summary, the adapted quasi-Poiseuille law helps us discuss transport regimes. Given the enormous dependence of the electrical properties with d𝑑ditalic_d, the scaling laws should be taken into account when designing any device.

V Conclusions

We found that the antidot superlattice is a convenient geometry for the realization of electron hydrodynamics. Certainly, the achievement is corroborated by the observation of the superballistic conduction with remarkable Gurzhi ratios (𝒢>20%𝒢percent20\mathcal{G}>20\%caligraphic_G > 20 %). The antidot geometry bends the electron flow so much that it overshadows the role of edge scattering. So, together with the control of the edge smoothness provided by the cryo-etching technique [25], we obtain an almost universal electron flow [16].

The observed hydrodynamic effects contribute to a better understanding of the hydrodynamic signatures. The simulations of the Boltzmann transport equation account for the ballistic-hydrodynamic transition, and enable us to study the tomographic dynamics of electrons [35]. Thereafter, we improve the conventional hydrodynamic description, and, in particular, the anti-Matthiessen rule that cannot account for our experimental results [15]. Particularly our formalism based on the Boltzmann equation is crucial to reproduce the intermittent superballistic effect experimentally found as a function of the applied magnetic field. Our description also sheds light on the classification of the hydrodynamic and ballistic transport regimes, by showing that two collective and hydrodynamic-like regions are achieved in presence of a magnetic field. The existence of collective regions right before and after the commensurability condition is in perfect agreement with the classification of transport regimes in Ref. [28]. Insight into the transport regimes is also attained with the scaling laws.

Both the advantages of the superlattices and the feasibility of their fabrication with lithography open an avenue to further optimize the geometries, look for novel hydrodynamic signatures, and explore materials to reduce the resistance in technological devices.

Acknowledgements.
We wish to acknowledge R. Brito and A. Hamilton for discussions. This work was supported by the “(MAD2D-CM)-UCM” project funded by Comunidad de Madrid, by the Recovery, Transformation and Resilience Plan, and by NextGenerationEU from the European Union, Agencia Estatal de Investigación of Spain (Grants PID2019-106820RB-C21/22 and PID2022-136285NB-C31/C32) and FEDER/Junta de Castilla y León Research Grant number SA106P23. J. E. acknowledges support from the Spanish Ministerio de Ciencia, Innovación y Universidades (Grant FPU22/01039). J. S.-S. acknowledges financial support from the Consejería de Educación, Junta de Castilla y León, and ERDF/FEDER. A. P.-R. acknowledges the financial support received from the Marie Skłodowska Curie-COFUND program under the Horizon 2020 research and innovation initiative of the European Union, within the framework of the USAL4Excellence program (grant agreement 101034371).

Appendix A Geometry optimization

The geometry of a device determines its hydrodynamic properties. Indeed, bending the electron flow enhances the hydrodynamic signatures. Figure S1 shows four antidot superlattices with different configurations and hole shapes, where the device current flows from left to right. These simulations were performed with a hydrodynamic model based on the Navier-Stokes equation for the characteristic lengths le=6dsubscript𝑙𝑒6𝑑l_{e}=6ditalic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT = 6 italic_d and lee=dsubscript𝑙𝑒𝑒𝑑l_{ee}=ditalic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT = italic_d, and perfect slip boundary conditions. In order to quantify the bending of the electron flow in every superlattice, we considered the scaling laws analyzed in Sec. IV, and evaluated the ratio R1/R2subscript𝑅1subscript𝑅2R_{1}/R_{2}italic_R start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / italic_R start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT, where R1(R2)subscript𝑅1subscript𝑅2R_{1}\leavevmode\nobreak\ (R_{2})italic_R start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_R start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) is the resistance of a device with antidot size d(2d)𝑑2𝑑d\leavevmode\nobreak\ (2d)italic_d ( 2 italic_d ), see Appendix K for a detailed description. Note that R1/R2>2subscript𝑅1subscript𝑅22R_{1}/R_{2}>2italic_R start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / italic_R start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT > 2 is a hydrodynamic signature, so we look for the largest R1/R2subscript𝑅1subscript𝑅2R_{1}/R_{2}italic_R start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / italic_R start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT as a quantitative criterion to support a maximized bending of the electron flow. Square antidots have sharp corners, so the electronic fluid follows almost straight trajectories without bending too much, as shown in Fig. S1(a). As a consequence R1/R2subscript𝑅1subscript𝑅2R_{1}/R_{2}italic_R start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / italic_R start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT is small. Moreover, given the difficulties in building and simulating samples with sharp corners, we avoid using square lattices. Superlattices with circles as antidots show a higher R1/R2subscript𝑅1subscript𝑅2R_{1}/R_{2}italic_R start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / italic_R start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT, especially when they are aligned at 45superscript4545^{\circ}45 start_POSTSUPERSCRIPT ∘ end_POSTSUPERSCRIPT respect to the current flow, see Fig. S1(d). In fact, misalignment with the latter avoids straight trajectories, which are less prone to bending around.

Refer to caption
Figure S1: Color maps of the electrical potential inside the sample with the streamlines for the electron fluid to visualize the bending of electron flow in antidot superlattices. Square antidots of side d𝑑ditalic_d with a separation of d𝑑ditalic_d yields (a) R1/R2=2.46subscript𝑅1subscript𝑅22.46R_{1}/R_{2}=2.46italic_R start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / italic_R start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = 2.46 in an aligned lattice and (b) R1/R2=2.52subscript𝑅1subscript𝑅22.52R_{1}/R_{2}=2.52italic_R start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / italic_R start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = 2.52 in a 45superscript4545^{\circ}45 start_POSTSUPERSCRIPT ∘ end_POSTSUPERSCRIPT rotated lattice with respect to the current flow. Circular antidots of diameter d𝑑ditalic_d separated by d𝑑ditalic_d yields (c) R1/R2=2.58subscript𝑅1subscript𝑅22.58R_{1}/R_{2}=2.58italic_R start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / italic_R start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = 2.58 in an aligned lattice and (d) R1/R2=2.63subscript𝑅1subscript𝑅22.63R_{1}/R_{2}=2.63italic_R start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / italic_R start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = 2.63 in a 45superscript4545^{\circ}45 start_POSTSUPERSCRIPT ∘ end_POSTSUPERSCRIPT rotated lattice with respect to the current flow.

Appendix B Sample preparation

The final heterostructure, consisting of a fully encapsulated single layer graphene with a graphite bottom gate was fabricated by means of the standard mechanical exfoliation on pristine crystals of hexagonal boron nitride (hBN) and graphite, subsequently deposited onto a silicon oxide wafer with top p-doped silicon oxide with nominal thickness of 290 nm. Top hBN and bottom hBN flakes had thicknesses of 888\,8nm and 606060\,60nm respectively and the graphite back gate consisted on a 151515\,15nm-thick layer. The thicknesses of these three layers of the stack were accurately measured with a profilometer model Bruker DekTakXT. For the characterization of the graphene monolayer, micro-Raman spectroscopy for the single-layer material was performed. For the stacking process of the heterostructure, a polycarbonate film was fabricated and deposited on polydimethylsiloxane. Top hBN was picked up at 50-60 and deposited on the graphene monolayer at 190. Employing the same technique, the hBN bottom flake was deposited onto the graphite back gate. Subsequently, the latter stack was annealed in vacuum at 350 to eliminate potential residues. The hBN top and graphene were finally picked up and deposited on the hBN bottom and graphite. A thorough micro-Raman map of the final heterostructure was performed to select the cleanest and defect-free region where the electron-beam lithography was performed.

Once the final heterostructure was fabricated (Fig. S2(a)) a premask process was carried out by EBL-SEM to remove excess flakes around it, thus avoiding possible electric shorts between pads and in order to facilitate further processing. A spin coating process was performed using homemade PMMA resist 5%percent\%% (by weight) in chlorobenzene. The premask was attacked with a cryo-etching system (Fig. S2(b)). Then, the stack was finally shaped by means of e-beam lithography into the final 10-terminal Hall bar (Fig. S2(c)). A subsequent last step of electron beam lithography + cryo etching process was carried out to define the antidot patterns within the Hall bar (Fig. S2(d)). Lastly, 10/55 nm Cr/Au contacts were deposited by e-beam evaporation (Fig. S2(e)) and the device was bonded on a LCC20 chip carrier for electrical characterization, finalizing the sample fabrication.

For the crucial lithographic step to define the periodic antidot lattice in the structure, a previous dose array process was carried out into a sacrificial hBN flake replicating the final design as shown in Fig. S3(a). There, doses ranging from 400 to 575 μ𝜇\muitalic_μC/cm2 were used showing that 550 μ𝜇\muitalic_μC/cm2 generates the optimal resolution for all three different antidot areas as shown in Figs. S3(b)-(d).

Refer to caption
Figure S2: Different steps in the device fabrication process and the final heterostructure (a) the premask after etching in (b) with the Hall bar finally defined in the heterostructure in panel (c). Panel (d) already displays the three different antidot regions with diameters of d=300, 200𝑑300200d=300,\,200italic_d = 300 , 200, and 100nm100nm100\,\rm nm100 roman_nm from the bottom region to the uppermost part of the panel. Finally, in panel (e) the final device after the evaporation of the Ohmic side contacts is shown.
Refer to caption
Figure S3: SEM electron micrograph for the dose array optimization process. (a) Antidot patterns obtained from doses ranging from 400 μ𝜇\muitalic_μC/cm2 (leftmost structure) to 575 μ𝜇\muitalic_μC/cm2. (b)-(d) From left to right: Antidot final structure onto a 8nm hBN flake displaying a great resolution and smooth edges for a dose of 550 μ𝜇\muitalic_μC/cm2 within d=300, 200𝑑300200d=300,\,200italic_d = 300 , 200, and 100nm100nm100\,\rm nm100 roman_nm regions.

Appendix C Theoretical model

Let us describe electrons as semiclassical particles [43, 29, 15], moving in a 2D device with well-defined position 𝒓=(x,y)𝒓𝑥𝑦{\bm{r}}=(x,y)bold_italic_r = ( italic_x , italic_y ) and wave vector 𝒌=k𝒖^𝒌(θ)𝒌𝑘subscript^𝒖𝒌𝜃{\bm{k}}=k\hat{\bm{u}}_{\bm{k}}(\theta)bold_italic_k = italic_k over^ start_ARG bold_italic_u end_ARG start_POSTSUBSCRIPT bold_italic_k end_POSTSUBSCRIPT ( italic_θ ), where 𝒖^𝒌(θ)(cosθ,sinθ)subscript^𝒖𝒌𝜃𝜃𝜃\hat{\bm{u}}_{\bm{k}}(\theta)\equiv\left(\cos\theta,\sin\theta\right)over^ start_ARG bold_italic_u end_ARG start_POSTSUBSCRIPT bold_italic_k end_POSTSUBSCRIPT ( italic_θ ) ≡ ( roman_cos italic_θ , roman_sin italic_θ ). Let 𝒗𝒗\bm{v}bold_italic_v denote the velocity of the electrons. The semiclassical description ignores quantum effects but accounts for the role of the geometry in the hydrodynamic effects. Then the Boltzmann transport equation reads

tf^+𝒗𝒓f^e(V^+𝒗×𝑩)𝒌f^=Γ[f^],subscript𝑡^𝑓𝒗subscript𝒓^𝑓𝑒Planck-constant-over-2-pi^𝑉𝒗𝑩subscript𝒌^𝑓Γdelimited-[]^𝑓\partial_{t}\hat{f}+{\bm{v}}\cdot\nabla_{{\bm{r}}}\hat{f}-\frac{e}{\hbar}\,% \big{(}-\nabla\hat{V}+{\bm{v}}\times{\bm{B}}\big{)}\cdot\nabla_{{\bm{k}}}\hat{% f}=\Gamma\left[\hat{f}\right]\ ,∂ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT over^ start_ARG italic_f end_ARG + bold_italic_v ⋅ ∇ start_POSTSUBSCRIPT bold_italic_r end_POSTSUBSCRIPT over^ start_ARG italic_f end_ARG - divide start_ARG italic_e end_ARG start_ARG roman_ℏ end_ARG ( - ∇ over^ start_ARG italic_V end_ARG + bold_italic_v × bold_italic_B ) ⋅ ∇ start_POSTSUBSCRIPT bold_italic_k end_POSTSUBSCRIPT over^ start_ARG italic_f end_ARG = roman_Γ [ over^ start_ARG italic_f end_ARG ] , (S1)

where 𝑩𝑩\bm{B}bold_italic_B is a perpendicular magnetic field and Γ[f]Γdelimited-[]𝑓\Gamma[f]roman_Γ [ italic_f ] is a collision operator

Γ[f]=ffele+fevenfeeevenleeeven+foddfeeoddleeodd.Γdelimited-[]𝑓𝑓subscript𝑓𝑒subscript𝑙𝑒superscript𝑓evensubscriptsuperscript𝑓even𝑒𝑒superscriptsubscript𝑙𝑒𝑒evensuperscript𝑓oddsubscriptsuperscript𝑓odd𝑒𝑒superscriptsubscript𝑙𝑒𝑒odd\Gamma[f]=\frac{f-f_{e}}{l_{e}}+\frac{f^{\mathrm{even}}-f^{\mathrm{even}}_{ee}% }{l_{ee}^{\mathrm{even}}}+\frac{f^{\mathrm{odd}}-f^{\mathrm{odd}}_{ee}}{l_{ee}% ^{\mathrm{odd}}}\ .roman_Γ [ italic_f ] = divide start_ARG italic_f - italic_f start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT end_ARG start_ARG italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT end_ARG + divide start_ARG italic_f start_POSTSUPERSCRIPT roman_even end_POSTSUPERSCRIPT - italic_f start_POSTSUPERSCRIPT roman_even end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT end_ARG start_ARG italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_even end_POSTSUPERSCRIPT end_ARG + divide start_ARG italic_f start_POSTSUPERSCRIPT roman_odd end_POSTSUPERSCRIPT - italic_f start_POSTSUPERSCRIPT roman_odd end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT end_ARG start_ARG italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_odd end_POSTSUPERSCRIPT end_ARG . (S2)

This includes scattering against impurities and phonons, towards the equilibrium distribution fesubscript𝑓𝑒f_{e}italic_f start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT with a mean free path lesubscript𝑙𝑒l_{e}italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT and electron-electron scattering towards an equilibrium distribution feesubscript𝑓𝑒𝑒f_{ee}italic_f start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT that moves with the fluids drift velocity. We split this last term with two relaxation times for the even and odd parts of the collision operator [35]. Thus, if f𝑓fitalic_f is expanded in angular harmonics f=nfnenθ𝑓subscript𝑛subscript𝑓𝑛superscript𝑒𝑛𝜃f=\sum_{n}f_{n}e^{n\theta}italic_f = ∑ start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT italic_f start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT italic_e start_POSTSUPERSCRIPT italic_n italic_θ end_POSTSUPERSCRIPT, fevensuperscript𝑓evenf^{\rm even}italic_f start_POSTSUPERSCRIPT roman_even end_POSTSUPERSCRIPT is the sum of the terms with even n𝑛nitalic_n and foddsuperscript𝑓oddf^{\rm odd}italic_f start_POSTSUPERSCRIPT roman_odd end_POSTSUPERSCRIPT includes the terms with odd n𝑛nitalic_n. They have different mean free paths, namely leeevensuperscriptsubscript𝑙𝑒𝑒evenl_{ee}^{\rm even}italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_even end_POSTSUPERSCRIPT and leeoddsuperscriptsubscript𝑙𝑒𝑒oddl_{ee}^{\rm odd}italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_odd end_POSTSUPERSCRIPT. Now, let us consider an isotropic conduction band with kFsubscript𝑘𝐹k_{F}italic_k start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT the Fermi wavenumber, vFsubscriptv𝐹{\rm v}_{F}roman_v start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT the Fermi velocity and m=kF/vF𝑚Planck-constant-over-2-pisubscript𝑘𝐹subscriptv𝐹m=\hbar k_{F}/{\rm v}_{F}italic_m = roman_ℏ italic_k start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT / roman_v start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT the cyclotron mass. We also assume a constant density of carriers n𝑛nitalic_n such that kF=πnsubscript𝑘𝐹𝜋𝑛k_{F}=\sqrt{\pi n}italic_k start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT = square-root start_ARG italic_π italic_n end_ARG, considering valley and spin degeneracy. Importantly, we assume that the relevant phenomena happen near the Fermi surface, so transport can be described just in terms of the θ𝜃\thetaitalic_θ direction. Thus, we define

g(𝒓,θ)𝑔𝒓𝜃\displaystyle g({\bm{r}},\theta)italic_g ( bold_italic_r , italic_θ ) =m0[f(𝒓,𝒌)fe(𝒌)]dk,absentPlanck-constant-over-2-pi𝑚superscriptsubscript0delimited-[]𝑓𝒓𝒌superscript𝑓𝑒𝒌differential-d𝑘\displaystyle=\frac{\hbar}{m}\int_{0}^{\infty}\Big{[}f({\bm{r}},{\bm{k}})-f^{e% }({\bm{k}})\Big{]}\,{\rm d}k\ ,= divide start_ARG roman_ℏ end_ARG start_ARG italic_m end_ARG ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT [ italic_f ( bold_italic_r , bold_italic_k ) - italic_f start_POSTSUPERSCRIPT italic_e end_POSTSUPERSCRIPT ( bold_italic_k ) ] roman_d italic_k , (S3a)
gee(𝒓,θ)subscript𝑔𝑒𝑒𝒓𝜃\displaystyle g_{ee}({\bm{r}},\theta)italic_g start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT ( bold_italic_r , italic_θ ) =m0[fee(𝒓,𝒌)fe(𝒌)]dk.absentPlanck-constant-over-2-pi𝑚superscriptsubscript0delimited-[]subscript𝑓𝑒𝑒𝒓𝒌superscript𝑓𝑒𝒌differential-d𝑘\displaystyle=\frac{\hbar}{m}\int_{0}^{\infty}\Big{[}f_{ee}({\bm{r}},{\bm{k}})% -f^{e}({\bm{k}})\Big{]}\,{\rm d}k\ .= divide start_ARG roman_ℏ end_ARG start_ARG italic_m end_ARG ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT [ italic_f start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT ( bold_italic_r , bold_italic_k ) - italic_f start_POSTSUPERSCRIPT italic_e end_POSTSUPERSCRIPT ( bold_italic_k ) ] roman_d italic_k . (S3b)

It is not difficult to show that gee(𝒓,θ)ux(𝒓)cosθ+uy(𝒓)sinθsimilar-to-or-equalssubscript𝑔𝑒𝑒𝒓𝜃subscript𝑢𝑥𝒓𝜃subscript𝑢𝑦𝒓𝜃g_{ee}({\bm{r}},\theta)\simeq u_{x}({\bm{r}})\cos\theta+u_{y}({\bm{r}})\sin\thetaitalic_g start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT ( bold_italic_r , italic_θ ) ≃ italic_u start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT ( bold_italic_r ) roman_cos italic_θ + italic_u start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT ( bold_italic_r ) roman_sin italic_θ, where the drift velocity is obtained as 𝒖(𝒓)=(1/π)02πg(𝒓,θ)𝒖^(θ)dθ𝒖𝒓1𝜋superscriptsubscript02𝜋𝑔𝒓𝜃^𝒖𝜃differential-d𝜃{\bm{u}}({\bm{r}})=(1/\pi)\int_{0}^{2\pi}g({\bm{r}},\theta)\hat{\bm{u}}(\theta% )\,{\rm d}\thetabold_italic_u ( bold_italic_r ) = ( 1 / italic_π ) ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 italic_π end_POSTSUPERSCRIPT italic_g ( bold_italic_r , italic_θ ) over^ start_ARG bold_italic_u end_ARG ( italic_θ ) roman_d italic_θ. We restrict ourselves to steady-state conditions, such that the non-equilibrium distribution function f(𝒓,𝒌)𝑓𝒓𝒌f({\bm{r}},{\bm{k}})italic_f ( bold_italic_r , bold_italic_k ) becomes independent of time. Hence, as described in Ref. [28], the Boltzmann equation for the non-equilibrium distribution function in a potential V(𝒓)𝑉𝒓V(\bm{r})italic_V ( bold_italic_r ) and a perpendicular magnetic field 𝑩𝑩{\bm{B}}bold_italic_B reduces to Eq.(1).

The geometry and edge scattering plays a crucial role in viscous electron flows, so the Boltzmann equation must be supplemented with the appropriate boundary condition. Two common choices [37, 28] are the diffusive (DF) edge that assumes g(θ)=0𝑔𝜃0g(\theta)=0italic_g ( italic_θ ) = 0 for all reflected electrons and the partially specular (PS) edge that reads

g(θ)𝑔𝜃\displaystyle g(\theta)italic_g ( italic_θ ) =g(θ)+𝒟sinθabsent𝑔𝜃𝒟𝜃\displaystyle=g(-\theta)+\mathcal{D}\sin\theta= italic_g ( - italic_θ ) + caligraphic_D roman_sin italic_θ
×[g(θ)2πsinθ0πsin2θg(θ)dθ],absentdelimited-[]𝑔𝜃2𝜋𝜃superscriptsubscript0𝜋superscript2superscript𝜃𝑔superscript𝜃differential-dsuperscript𝜃\displaystyle\times\left[g(-\theta)-\frac{2}{\pi}\,\sin\theta\int_{0}^{\pi}% \sin^{2}\theta^{\prime}g(-\theta^{\prime})\,{\rm d}\theta^{\prime}\right]\ ,× [ italic_g ( - italic_θ ) - divide start_ARG 2 end_ARG start_ARG italic_π end_ARG roman_sin italic_θ ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_π end_POSTSUPERSCRIPT roman_sin start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_θ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT italic_g ( - italic_θ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ) roman_d italic_θ start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ] , (S4)

where 0<θ<π0𝜃𝜋0<\theta<\pi0 < italic_θ < italic_π are the reflected electrons and π<θ<0𝜋𝜃0-\pi<\theta<0- italic_π < italic_θ < 0 are the incident ones. For the sake of simplicity, we wrote the boundary condition for an edge parallel to θ=0𝜃0\theta=0italic_θ = 0. Here, 𝒟πh2hkF31𝒟𝜋superscript2superscriptsuperscriptsubscript𝑘𝐹3less-than-or-similar-to1\mathcal{D}\equiv\sqrt{\pi}h^{2}h^{\prime}k_{F}^{3}\lesssim 1caligraphic_D ≡ square-root start_ARG italic_π end_ARG italic_h start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_h start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT italic_k start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT ≲ 1 is the dispersion coefficient, with hhitalic_h the edge’s bumps mean height and hsuperscripth^{\prime}italic_h start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT its correlation length. If 𝒟1much-less-than𝒟1\mathcal{D}\ll 1caligraphic_D ≪ 1 the edge is almost specular (see Ref. [28] for details of the theoretical model).

Appendix D Numerical methods

We solve the Boltzmann equation numerically with a conformal Galerkin finite element method [44]. We approximate the solution as

g(𝒓,θ)=n=1Nm=1Mϕn(𝒓)φm(θ)𝑔𝒓𝜃superscriptsubscript𝑛1𝑁superscriptsubscript𝑚1𝑀subscriptitalic-ϕ𝑛𝒓subscript𝜑𝑚𝜃g(\bm{r},\theta)=\sum_{n=1}^{N}\sum_{m=1}^{M}\phi_{n}(\bm{r})\varphi_{m}(\theta)italic_g ( bold_italic_r , italic_θ ) = ∑ start_POSTSUBSCRIPT italic_n = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT ∑ start_POSTSUBSCRIPT italic_m = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_M end_POSTSUPERSCRIPT italic_ϕ start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ( bold_italic_r ) italic_φ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( italic_θ ) (S5)

For the spatial part, {ϕn}n=1Nsuperscriptsubscriptsubscriptitalic-ϕ𝑛𝑛1𝑁\{\phi_{n}\}_{n=1}^{N}{ italic_ϕ start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT } start_POSTSUBSCRIPT italic_n = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT is the set of tent functions and the products of two tent functions defined on a triangular mesh [45] for the antidot geometry. We impose a maximum triangle size h<0.1d0.1𝑑h<0.1ditalic_h < 0.1 italic_d (or h<0.2d0.2𝑑h<0.2ditalic_h < 0.2 italic_d under a magnetic field) for which N2000similar-to𝑁2000N\sim 2000italic_N ∼ 2000 and convergence is ensured. For the angular part, {ϕm}n=1Msuperscriptsubscriptsubscriptitalic-ϕ𝑚𝑛1𝑀\{\phi_{m}\}_{n=1}^{M}{ italic_ϕ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT } start_POSTSUBSCRIPT italic_n = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_M end_POSTSUPERSCRIPT is a set of periodic functions defined on the interval [0,2π)02𝜋[0,2\pi)[ 0 , 2 italic_π ) and we use M=16𝑀16M=16italic_M = 16 and M=32𝑀32M=32italic_M = 32. We write the weak formulation of (1), add an equation to set a uniform density of carriers and solve the resulting linear system iteratively with a least square approximation in Matlab. At the edges, we impose the boundary condition (S4) for reflected electrons. We solve the system on a rectangular cell of size 22d×2d22𝑑2𝑑2\sqrt{2}d\times\sqrt{2}d2 square-root start_ARG 2 end_ARG italic_d × square-root start_ARG 2 end_ARG italic_d, and impose periodic boundary conditions. We set the potential difference between two cells across the longitudinal direction, and determine the Hall potential across the transverse direction by adding an additional equation that imposes no net current across the transverse direction. The hydrodynamic model that we used for the geometry optimization was also solved using finite elements [28]. We used a Runge-Kutta 4 method to solve the electronic trajectories in the streamlines, and numerical integration to find the total current. Given the density of carriers n=0.3×1012cm2𝑛0.3superscript1012superscriptcm2n=0.3\times 10^{12}\,\rm cm^{-2}italic_n = 0.3 × 10 start_POSTSUPERSCRIPT 12 end_POSTSUPERSCRIPT roman_cm start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT which was determined from the quantum Hall effect, the reduced units from the numerical simulation are translated into resistances. The resistances include a geometrical aspect ratio L/W𝐿𝑊L/Witalic_L / italic_W, where L3μm𝐿3𝜇mL\approx 3\,\rm\mu mitalic_L ≈ 3 italic_μ roman_m is the length of the region containing the antidots and W=3μm𝑊3𝜇mW=3\,\rm\mu mitalic_W = 3 italic_μ roman_m its width, so L/W=1𝐿𝑊1L/W=1italic_L / italic_W = 1 .

Appendix E Sample characterization

We need to determine the mean free path lesubscript𝑙𝑒l_{e}italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT to characterize the electrical properties of the antidot superlattice. In this section, we show the characterization of the sample and determine the values of lesubscript𝑙𝑒l_{e}italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT.

Refer to caption
Figure S4: Sample characterization. (a)-(c) Experimental Dirac peaks at T=1.5K𝑇1.5KT=1.5\,\rm Kitalic_T = 1.5 roman_K. (d)-(f) Estimation of the mean free path for collisions against impurities out of the Dirac peaks.

Figure S4(a-c) shows the Dirac peaks for the three regions (d=100, 200𝑑100200d=100,\,200italic_d = 100 , 200, and 300nm300nm300\,\rm nm300 roman_nm) of the device at T=1.5K𝑇1.5KT=1.5\,\rm Kitalic_T = 1.5 roman_K. Notice the peaks broaden due to the nanostructuring of the device, which can add more impurities and tensions near the edges as well as charge inhomogeneity discussed in Appendix F. Since electron-electron collisions are negligible at this temperature, the device resistance mainly depends on the parameter lesubscript𝑙𝑒l_{e}italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT and the particular lattice geometry. Hence, we can use simulations of the Boltzmann transport equation to estimate the mean free paths lesubscript𝑙𝑒l_{e}italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT that give rise to the experimental resistances. Our results are shown in Fig. S4(d-f), where we find maxima near Vg=1.5Vsubscript𝑉𝑔1.5VV_{g}=1.5\,\rm Vitalic_V start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT = 1.5 roman_V. Interestingly, this maximum coincides with the maximum of the ratio R1/R2subscript𝑅1subscript𝑅2R_{1}/R_{2}italic_R start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / italic_R start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT in the quasi-Poiseuille law, see Fig. 4(a).

We notice that the resistance is mainly determined by the antidot geometry and not by scattering against impurities. On the one hand, this causes the estimated value of lesubscript𝑙𝑒l_{e}italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT to be enormously sensitive to R𝑅Ritalic_R and to the experimental inaccuracies on its measurement, with estimations that may differ for each region. On the other hand, it is an advantage for this experiment, since the exact electrical properties are quite robust, not being very sensitive to lesubscript𝑙𝑒l_{e}italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT.

In the low-temperature limit, we assume le=700nmsubscript𝑙𝑒700nml_{e}=700\,\rm nmitalic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT = 700 roman_nm at Vg=3.9Vsubscript𝑉𝑔3.9VV_{g}=3.9\,\rm Vitalic_V start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT = 3.9 roman_V, similar to the one achieved in other graphene nanostructures [30], which also accounts for the additional damage to the structure due to the structuring process. Last, we include phonon scattering, with a Bloch-Grüneisen temperature of 54K54K54\,\rm K54 roman_K to account for the low-temperature non-linear dependence [32]. This assumes that phonons do not affect the electrical properties below 10Kabsent10K\approx 10\,\rm K≈ 10 roman_K and a high-temperature linear dependence [31] of le,ph170μmK/Tsubscript𝑙𝑒𝑝170𝜇mK𝑇l_{e,ph}\approx 170\,{\rm\mu m\,K}/Titalic_l start_POSTSUBSCRIPT italic_e , italic_p italic_h end_POSTSUBSCRIPT ≈ 170 italic_μ roman_m roman_K / italic_T.

Appendix F Charge inhomogeneity

Refer to caption
Figure S5: Charge inhomogeneity in antidot superlattices. (a)-(b) Transverse section of a graphene conductor, encapsulated between two layers of hBN, on top of a back gate at a given potential, for d=100𝑑100d=100italic_d = 100 and 300nm300nm300\,\rm nm300 roman_nm hole sizes. (c)-(d) Induced density of carriers in the flake, normalized to the n0subscript𝑛0n_{0}italic_n start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT density corresponding to an ideal capacitor. (e)-(f) Color maps of the electrical potential inside the sample with the streamlines for the electron fluid simulated with the Boltzmann transport equation. Left (right) panels show results for antidots of diameter d𝑑ditalic_d (1.1d1.1𝑑1.1d1.1 italic_d).

In this experiment, we engineer the geometry of the device in order to bend the electron flow. As a side effect, the back gate induces an inhomogeneous density of carriers n𝑛nitalic_n. Let us quantify the underlying electrostatic effect for a graphene flake encapsulated between two layers of hBN, a dielectric medium with ε3.5similar-to𝜀3.5\varepsilon\sim 3.5italic_ε ∼ 3.5, with the top and bottom thicknesses of 8nm8nm8\,\rm nm8 roman_nm and 60nm60nm60\,\rm nm60 roman_nm respectively. For this purpose, we consider a simplified two dimensional problem by assuming that the holes extend across the direction perpendicular to the figure, where we can solve the Poisson equation with a finite differences numerical method. Figure S5 shows the solution: panels (a) and (b) show the electrostatic potential and (c) and (d) account for the induced density of carriers. Indeed, the density of carriers is inhomogeneous, and the accumulation next to the edges would be more relevant for superlattices of shorter d𝑑ditalic_d.

Charge inhomogeneity flattens the Dirac peak of geometrically engineered samples, shortens the effective mean free paths lesubscript𝑙𝑒l_{e}italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT, and leads to small discrepancies when the same path is used for all the antidot regions. We quantify the effect of charge inhomogeneity by studying the most dramatic scenario where the region adjacent to the holes does not contribute to electrical conduction, for example, when it is depleted of electrons. Therefore, we can simply simulate a system with the same center-to-center distance, but with bigger antidots. The results for archetypal experimental parameters at low temperature and d/le=0.15𝑑subscript𝑙𝑒0.15d/l_{e}=0.15italic_d / italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT = 0.15 are shown in panels (e) and (f). We find that even a major change of 10%percent1010\%10 % in the antidot diameter only results in a change of 17%percent1717\%17 % in the resistance. Although it is quite remarkable, it is not enough to explain the most relevant results of our work Last, we notice that, due to the doping of the graphene flake with holes, the applied voltage that we need to achieve a n𝑛-n- italic_n density of holes is smaller, in modulus, than to achieve the same n𝑛nitalic_n density of electrons. Thus, the inhomogeneity is less noticeable when working with holes. This may explain the better agreement with the simulations when working for holes, see Fig. 2(b) and (c). Last, if the scaling laws were caused by the inhomogeneity, they should be dramatically different for electrons and holes. Consequently, the fact that the scaling R1/R2>2subscript𝑅1subscript𝑅22R_{1}/R_{2}>2italic_R start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / italic_R start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT > 2 prevails both for types of carriers (see Fig. 4(a)) indicates that it is not due to charge inhomogeneity.

Appendix G Charge neutrality

Refer to caption
Figure S6: Resistance measured near the charge neutrality point. (a) Arrhenius law characterizes thermally activated conduction. (b)-(c) Resistance as a function of the temperature for two gate voltages.

Superballistic conduction is often studied far away from the charge neutrality point [9, 16, 46]. Conversely, in this section, we will focus on superballistic conduction near the charge neutrality point. Figure S6(a) shows an Arrhenius plot near the charge neutrality point. On the one hand, the high-temperature region is dominated by the thermal activation of carriers, as it is characteristic of an intrinsic semiconductor or an insulator. On the other hand, the experimental data does not fit to the Arrhenius law for low temperatures T100Kless-than-or-similar-to𝑇100KT\lesssim 100\,\rm Kitalic_T ≲ 100 roman_K [31], meaning that the low-temperature resistance is not dominated by thermal activation. Indeed, Fig. S6(b) shows two distinguishable steps in the resistance: superballistic conduction occurs at T100Kless-than-or-similar-to𝑇100KT\lesssim 100\,\rm Kitalic_T ≲ 100 roman_K and thermal activation for T100Kgreater-than-or-equivalent-to𝑇100KT\gtrsim 100\,\rm Kitalic_T ≳ 100 roman_K. Fig. S6(c) accounts for another density of carriers, and it also shows distinguishable superballistic conduction at T100Kless-than-or-similar-to𝑇100KT\lesssim 100\,\rm Kitalic_T ≲ 100 roman_K. Therefore, we conclude that the Gurzhi ratio 𝒢=1R(50K)/R(1.5K)𝒢1𝑅50K𝑅1.5K\mathcal{G}=1-R(50{\,\rm K})/R(1.5{\,\rm K})caligraphic_G = 1 - italic_R ( 50 roman_K ) / italic_R ( 1.5 roman_K ) effectively decouples superballistic conduction from thermal activation.

Appendix H Boundary scattering and universality

One of the main questions regarding electron hydrodynamics is that of the edge scattering, which determines some electrical properties [37]. The cryo-etching technique gives control of the graphene edge, with bumps of mean height hhitalic_h and correlation lengths hsuperscripth^{\prime}italic_h start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT in the order of the nm [25]. Since kF11nmmuch-greater-thansuperscriptsubscript𝑘𝐹11nmk_{F}^{-1}\gg 1\,\rm nmitalic_k start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ≫ 1 roman_nm unless we are next to the Dirac peak, we assume 𝒟=πh2hkF31𝒟𝜋superscript2superscriptsuperscriptsubscript𝑘𝐹3much-less-than1\mathcal{D}=\sqrt{\pi}h^{2}h^{\prime}k_{F}^{3}\ll 1caligraphic_D = square-root start_ARG italic_π end_ARG italic_h start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_h start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT italic_k start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT ≪ 1 for the dispersion coefficient in the simulations. Indeed, we compute all the results in the article with 𝒟=0.011𝒟0.01much-less-than1\mathcal{D}=0.01\ll 1caligraphic_D = 0.01 ≪ 1, which is a practically perfect slip. In this sense, this is the same as the boundary condition for electrostatically defined edges in GaAs heterostructures [16], providing a mechanism for the control of edge scattering in graphene, whose gaplessness does not allow electrostatically defined edges. This provides results in agreement with the experiment.

Refer to caption
Figure S7: Robustness of the results regardless of boundary scattering. (a) Simulations of the Boltzmann transport equation for the resistance as a function of the electron-electron collision rate d/lee𝑑subscript𝑙𝑒𝑒d/l_{ee}italic_d / italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT, being le=700nmsubscript𝑙𝑒700nml_{e}=700\rm\,nmitalic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT = 700 roman_nm. We show the result for a diffusive (DF) edge and several values of 𝒟𝒟\mathcal{D}caligraphic_D in a partially specular edge. (b) Simulations of the magnetoresistance for some cases shown in (a). (c)-(d) Streamlines and velocity profiles for an specular boundary with 𝒟=0.01𝒟0.01\mathcal{D}=0.01caligraphic_D = 0.01 and a diffusive one, for le=3dsubscript𝑙𝑒3𝑑l_{e}=3ditalic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT = 3 italic_d, leedmuch-greater-thansubscript𝑙𝑒𝑒𝑑l_{ee}\gg ditalic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT ≫ italic_d.

In this section, we investigate the role of edge scattering by performing simulations of the Boltzmann transport equation. Figure S7(a-b) shows the resistance as a function of the electron-electron scattering rate and the magnetic field, for several boundary conditions. The result is almost constant for several values of 𝒟𝒟\mathcal{D}caligraphic_D in a partial slip boundary condition. A perfect slip boundary condition is shown in Fig. S7(c), where the electron streamlines perfectly follow the shape of the holes. The electrical properties are not very different for the diffusive edge, whose velocity profile is shown in Fig. S7, and whose streamlines are slightly separated from the holes.

Consequently, we prove that the results are robust regardless of the edge scattering mechanism. Both the edge control of the cryo-etching technique and the fact that the electron bending is mainly controlled by the antidot geometry and not by edge scattering, ensure an almost universal viscous electron flow in our graphene devices.

Appendix I Reproducibility

Refer to caption
Figure S8: Superballistic conduction only arises in regions with antidots. (a) Optical image of the second device, with a pristine region and antidots of diameters d=100, 200𝑑100200d=100,\,200italic_d = 100 , 200, and 300nm300nm300\,\rm nm300 roman_nm. (b) Dirac peaks in all the samples. (c)-(d) Resistance as a function of the temperature for the densities of carriers n=0.3×1012cm2𝑛0.3superscript1012superscriptcm2n=0.3\times 10^{12}\rm\,cm^{-2}italic_n = 0.3 × 10 start_POSTSUPERSCRIPT 12 end_POSTSUPERSCRIPT roman_cm start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT and n=1.25×1012cm2𝑛1.25superscript1012superscriptcm2n=1.25\times 10^{12}\rm\,cm^{-2}italic_n = 1.25 × 10 start_POSTSUPERSCRIPT 12 end_POSTSUPERSCRIPT roman_cm start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT, respectively.

In order to further support our findings, we fabricated a second device. The new Hall bar includes the same three regions with antidots, as well as a pristine Hall-bar-like region where no antidots were defined. The experimental procedure is described in Appendix B, but now a standard 290-nm-thick SiO2 substrate was used to enable higher gate voltages without dielectric breaking. While the top hBN flake was kept comparable to the one used for the graphite back-gate sample (10similar-toabsent10\sim 10\,∼ 10nm) the bottom one was significantly thinner (30similar-toabsent30\sim 30\,∼ 30nm) since there was no risk of electrical short from the graphene towards an underlying metallic layer (graphite). Fig. S8(a) shows the Dirac peaks of the new sample, with the resistance as a function of the gate voltage, which is already referred to the center of the peaks. Fig. S8(b) and (c) shows the superballistic conduction in the regions with antidots, being qualitatively similar to the ones studied in Fig.2 for different densities of carriers. Most importantly no superballistic effect arise in the pristine region. In conclusion, the values of dR/dT<0d𝑅d𝑇0{\rm d}R/{\rm d}T<0roman_d italic_R / roman_d italic_T < 0 in the d=100, 200𝑑100200d=100,\,200italic_d = 100 , 200, and 300nm300nm300\,\rm nm300 roman_nm regions, further justifies the finding of superballistic conduction in antidot superlattices. Moreover, the absence of the effect in the pristine region with no antidots suggests that the geometrical engineering of the device is responsible for the superballistic conduction.

Appendix J Tomographic and hydrodynamic regimes

Refer to caption
Figure S9: Tomographic and hydrodynamic transport. (a) Simulations of the Boltzmann equation for electron streamlines and current profiles, being d/le=0.25𝑑subscript𝑙𝑒0.25d/l_{e}=0.25italic_d / italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT = 0.25, d/leeeven=0.4𝑑superscriptsubscript𝑙𝑒𝑒even0.4d/l_{ee}^{\rm even}=0.4italic_d / italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_even end_POSTSUPERSCRIPT = 0.4 and d/leeodd=0𝑑superscriptsubscript𝑙𝑒𝑒odd0d/l_{ee}^{\rm odd}=0italic_d / italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_odd end_POSTSUPERSCRIPT = 0 in the tomographic approach followed in this paper. (b) Same with d/leeeven=d/leeodd=0.4𝑑superscriptsubscript𝑙𝑒𝑒even𝑑superscriptsubscript𝑙𝑒𝑒odd0.4d/l_{ee}^{\rm even}=d/l_{ee}^{\rm odd}=0.4italic_d / italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_even end_POSTSUPERSCRIPT = italic_d / italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_odd end_POSTSUPERSCRIPT = 0.4 and d/leeodd=0𝑑superscriptsubscript𝑙𝑒𝑒odd0d/l_{ee}^{\rm odd}=0italic_d / italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_odd end_POSTSUPERSCRIPT = 0 in a fully hydrodynamic description. (c) Magnetoresistance simulations in the fully hydrodynamic approach.

We work under the tomographic approach [35, 36], where leeeven=leesuperscriptsubscript𝑙𝑒𝑒evensubscript𝑙𝑒𝑒l_{ee}^{\mathrm{even}}=l_{ee}italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_even end_POSTSUPERSCRIPT = italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT and leeoddleemuch-greater-thansuperscriptsubscript𝑙𝑒𝑒oddsubscript𝑙𝑒𝑒l_{ee}^{\mathrm{odd}}\gg l_{ee}italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_odd end_POSTSUPERSCRIPT ≫ italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT. Let us compare against the conventional hydrodynamic description leeeven=leeodd=leesuperscriptsubscript𝑙𝑒𝑒evensuperscriptsubscript𝑙𝑒𝑒oddsubscript𝑙𝑒𝑒l_{ee}^{\mathrm{even}}=l_{ee}^{\mathrm{odd}}=l_{ee}italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_even end_POSTSUPERSCRIPT = italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_odd end_POSTSUPERSCRIPT = italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT. This comparison is possible going beyond the approximated anti-Matthiessen rule [15, 9] and solving the Boltzmann equation. Figure S9(a) and (b) shows the current distribution under both approaches. Also, Fig. S9(c) shows the magnetic response under a purely hydrodynamic description, to be compared with Fig. 3(d) for the tomographic regime. As the magnetic field increases, it rotates the velocity of the electrons and makes the distinction between even and odd parity modes less noticeable. However, the behavior is different in the absence of a magnetic field. Therefore, simulations show that there is a difference between the electrical response depending on the microscopic scattering mechanisms.

Appendix K Scaling laws

Hydrodynamic flow features properties that strongly depend on the scale. Let us analyze some relevant scaling laws affecting the device resistance in the hydrodynamic regime, R1/dαproportional-to𝑅1superscript𝑑𝛼R\propto 1/d^{\alpha}italic_R ∝ 1 / italic_d start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT. In the latter, the Boltzmann equation can be reduced to a modified Navier-Stokes equation [28] together with the continuity equation as follows

𝒖=𝒖absent\displaystyle\nabla\cdot{\bm{u}}=∇ ⋅ bold_italic_u = 00\displaystyle 0 (S6)
ν2𝒖+vFle𝒖=𝜈superscript2𝒖subscriptv𝐹subscript𝑙𝑒𝒖absent\displaystyle-\nu\nabla^{2}{\bm{u}}+\frac{{\rm v}_{F}}{l_{e}}{\bm{u}}=- italic_ν ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT bold_italic_u + divide start_ARG roman_v start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT end_ARG start_ARG italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT end_ARG bold_italic_u = emV𝑒𝑚𝑉\displaystyle\frac{e}{m}\nabla Vdivide start_ARG italic_e end_ARG start_ARG italic_m end_ARG ∇ italic_V (S7)

where 𝒖𝒖\bm{u}bold_italic_u is the fluid drift velocity and ν=vF(le1+lee1)1/4𝜈subscriptv𝐹superscriptsuperscriptsubscript𝑙𝑒1superscriptsubscript𝑙𝑒𝑒114\nu={\rm v}_{F}\cdot(l_{e}^{-1}+l_{ee}^{-1})^{-1}/4italic_ν = roman_v start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT ⋅ ( italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT + italic_l start_POSTSUBSCRIPT italic_e italic_e end_POSTSUBSCRIPT start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT / 4 is the viscosity. For simplicity, we write the equation in the absence of a magnetic field. Equation S7 includes a dissipative term that accounts for collisions against impurities and phonons. These equations are solved for a particular boundary condition, that in our experimental setup describes a perfectly specular edge. This imposes the no-trespassing u=0subscript𝑢perpendicular-to0u_{\perp}=0italic_u start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT = 0 and the perfect slip u=0subscriptperpendicular-tosubscript𝑢parallel-to0\partial_{\perp}u_{\parallel}=0∂ start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT italic_u start_POSTSUBSCRIPT ∥ end_POSTSUBSCRIPT = 0 conditions, where usubscript𝑢parallel-tou_{\parallel}italic_u start_POSTSUBSCRIPT ∥ end_POSTSUBSCRIPT and usubscript𝑢perpendicular-tou_{\perp}italic_u start_POSTSUBSCRIPT ⟂ end_POSTSUBSCRIPT are the components of the velocity parallel and perpendicular to the edge.

Let us first study the physical situation such that the viscous term dominates

ν2𝒖V=emVV.𝜈superscript2subscript𝒖𝑉𝑒𝑚subscript𝑉𝑉-\nu\nabla^{2}{\bm{u}_{V}}=\frac{e}{m}\nabla V_{V}\ .- italic_ν ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT bold_italic_u start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT = divide start_ARG italic_e end_ARG start_ARG italic_m end_ARG ∇ italic_V start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT . (S8)

In this case the solution for any geometrical size d𝑑ditalic_d reads as 𝒖V(𝒓)=𝒖~𝑽(𝒓/d)subscript𝒖𝑉𝒓subscriptbold-~𝒖𝑽𝒓𝑑{\bm{u}_{V}}({\bm{r}})={\bm{\tilde{u}_{V}}}({\bm{r}}/d)bold_italic_u start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT ( bold_italic_r ) = overbold_~ start_ARG bold_italic_u end_ARG start_POSTSUBSCRIPT bold_italic_V end_POSTSUBSCRIPT ( bold_italic_r / italic_d ) and VV(𝒓)=d1V~V(𝒓/d)subscript𝑉𝑉𝒓superscript𝑑1subscript~𝑉𝑉𝒓𝑑{V_{V}}({\bm{r}})=d^{-1}\tilde{V}_{V}({\bm{r}}/d)italic_V start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT ( bold_italic_r ) = italic_d start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT over~ start_ARG italic_V end_ARG start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT ( bold_italic_r / italic_d ), where 𝒖~Vsubscript~𝒖𝑉\tilde{\bm{u}}_{V}over~ start_ARG bold_italic_u end_ARG start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT and V~Vsubscript~𝑉𝑉\tilde{V}_{V}over~ start_ARG italic_V end_ARG start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT are functions that do not depend on d𝑑ditalic_d. Notice that these functions fulfill the perfect slip boundary condition. Furthermore, although this reasoning would not be valid for partial slip boundary conditions, it also applies to the n- slip u=0subscript𝑢parallel-to0u_{\parallel}=0italic_u start_POSTSUBSCRIPT ∥ end_POSTSUBSCRIPT = 0 condition commonly used for the derivation of the Poiseuille law in conventional fluids. As a consequence of the scale dependence, VV(𝒓)1/dproportional-tosubscript𝑉𝑉𝒓1𝑑V_{V}({\bm{r}})\propto 1/ditalic_V start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT ( bold_italic_r ) ∝ 1 / italic_d is the voltage drop in a region of length d𝑑ditalic_d, and the resistance for devices of the same geometry and different sizes scales as RV1/d2proportional-tosubscript𝑅𝑉1superscript𝑑2R_{V}\propto 1/d^{2}italic_R start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT ∝ 1 / italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, i.e. α=2𝛼2\alpha=2italic_α = 2.

On the contrary, if the diffusive term dominates Eq. S7 reduces to

vFle𝒖D=emVD,subscriptv𝐹subscript𝑙𝑒subscript𝒖𝐷𝑒𝑚subscript𝑉𝐷\frac{{\rm v}_{F}}{l_{e}}{\bm{u}_{D}}=\frac{e}{m}\nabla V_{D}\ ,divide start_ARG roman_v start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT end_ARG start_ARG italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT end_ARG bold_italic_u start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT = divide start_ARG italic_e end_ARG start_ARG italic_m end_ARG ∇ italic_V start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT , (S9)

which is solved by the family of solutions 𝒖D(𝒓)=𝒖~D(𝒓/d)subscript𝒖𝐷𝒓subscriptbold-~𝒖𝐷𝒓𝑑{\bm{u}_{D}}({\bm{r}})={\bm{\tilde{u}}_{D}}({\bm{r}}/d)bold_italic_u start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT ( bold_italic_r ) = overbold_~ start_ARG bold_italic_u end_ARG start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT ( bold_italic_r / italic_d ) and VD(𝒓)=dV~D(𝒓/d)subscript𝑉𝐷𝒓𝑑subscript~𝑉𝐷𝒓𝑑{V_{D}}({\bm{r}})=d\leavevmode\nobreak\ \tilde{V}_{D}({\bm{r}}/d)italic_V start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT ( bold_italic_r ) = italic_d over~ start_ARG italic_V end_ARG start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT ( bold_italic_r / italic_d ). Independently of the boundary condition, this results in a constant resistance RDsubscript𝑅𝐷R_{D}italic_R start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT, i.e. α=0𝛼0\alpha=0italic_α = 0.

In the general scenario considered in Eq. S7, there is no trivial expression for α𝛼\alphaitalic_α. However, we propose to make the following ansatz: 𝒖V(𝒓)=𝒖D(𝒓)=𝒖(𝒓)subscript𝒖𝑉𝒓subscript𝒖𝐷𝒓𝒖𝒓{\bm{u}_{V}}({\bm{r}})={\bm{u}_{D}}({\bm{r}})={\bm{u}}({\bm{r}})bold_italic_u start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT ( bold_italic_r ) = bold_italic_u start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT ( bold_italic_r ) = bold_italic_u ( bold_italic_r ). Namely, there is a single velocity field 𝒖(𝒓)𝒖𝒓{\bm{u}}({\bm{r}})bold_italic_u ( bold_italic_r ) that solves both the equation with the viscous term and the dissipative term. The ansatz seems reasonable if we look at Fig. 4(d) and (e), whose streamlines are quite similar. However, the associated potentials and resistances may not be the same, so we define the bending coefficient β=RV/RD𝛽subscript𝑅𝑉subscript𝑅𝐷\beta=R_{V}/R_{D}italic_β = italic_R start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT / italic_R start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT. The latter only depends on the geometry such that a higher β𝛽\betaitalic_β corresponds to a more irregular electron flow. In order to estimate β𝛽\betaitalic_β we consider the total resistance in Eq. S7 as proportional to the sum of two terms (associated with the limiting cases considered in Eq. S8 and Eq. S9) that are independently simulated:

RβνvFd2+1leα21+vFd2βνle.proportional-to𝑅𝛽𝜈subscriptv𝐹superscript𝑑21subscript𝑙𝑒𝛼similar-to21subscriptv𝐹superscript𝑑2𝛽𝜈subscript𝑙𝑒R\propto\frac{\beta\nu}{{\rm v}_{F}d^{2}}+\frac{1}{l_{e}}\ \Rightarrow\ \alpha% \sim\frac{2}{1+\frac{{\rm v}_{F}d^{2}}{\beta\nu l_{e}}}\ .italic_R ∝ divide start_ARG italic_β italic_ν end_ARG start_ARG roman_v start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + divide start_ARG 1 end_ARG start_ARG italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT end_ARG ⇒ italic_α ∼ divide start_ARG 2 end_ARG start_ARG 1 + divide start_ARG roman_v start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG italic_β italic_ν italic_l start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT end_ARG end_ARG . (S10)

As a result, we estimate β4.4similar-to-or-equals𝛽4.4\beta\simeq 4.4italic_β ≃ 4.4 for the particular geometry considered in the experiments. Last, notice that this expression is consistent with an exponent 0<α<20𝛼20<\alpha<20 < italic_α < 2.

References

  • Gurzhi [1963] R. N. Gurzhi, Minimum of resistance in impurity-free conductors, J. Exp. Theor. Phys. 17, 521 (1963).
  • Gurzhi [1968] R. N. Gurzhi, Hydrodynamic effects in solids at low temperature, Sov. Phys. Usp 11, 255 (1968).
  • Polini and Geim [2020] M. Polini and A. K. Geim, Viscous electron fluids, Phys. Today 73, 28 (2020).
  • Narozhny [2022] B. N. Narozhny, Hydrodynamic approach to two-dimensional electron systems, Riv. del Nuovo Cim. 45, 661 (2022).
  • Varnavides et al. [2023] G. Varnavides, A. Yacoby, C. Felser, and P. Narang, Charge transport and hydrodynamics in materials, Nat. Rev. Mater. 8, 726 (2023).
  • Bandurin et al. [2016] D. A. Bandurin, I. Torre, R. K. Kumar, M. Ben Shalom, A. Tomadin, A. Principi, G. H. Auton, E. Khestanova, K. S. Novoselov, I. V. Grigorieva, L. A. Ponomarenko, A. K. Geim, and M. Polini, Negative local resistance caused by viscous electron backflow in graphene, Science 351, 1055 (2016).
  • Bandurin et al. [2018] D. A. Bandurin, A. V. Shytov, L. S. Levitov, R. K. Kumar, A. I. Berdyugin, M. Ben Shalom, I. V. Grigorieva, A. K. Geim, and G. Falkovich, Fluidity onset in graphene, Nat. Commun. 9, 4533 (2018).
  • Sulpizio et al. [2019] J. A. Sulpizio, L. Ella, A. Rozen, J. Birkbeck, D. J. Perello, D. Dutta, M. Ben-Shalom, T. Taniguchi, K. Watanabe, T. Holder, et al., Visualizing poiseuille flow of hydrodynamic electrons, Nature 576, 75 (2019).
  • Krishna Kumar et al. [2017] R. Krishna Kumar, D. A. Bandurin, F. M. D. Pellegrino, Y. Cao, A. Principi, H. Guo, G. H. Auton, M. Ben Shalom, L. A. Ponomarenko, G. Falkovich, K. Watanabe, T. Taniguchi, I. Grigorieva, L. S. Levitov, M. Polini, and A. Geim, Superballistic flow of viscous electron fluid through graphene constrictions, Nat. Phys. 13, 1182 (2017).
  • Kumar et al. [2022] C. Kumar, J. Birkbeck, J. A. Sulpizio, D. Perello, T. Taniguchi, K. Watanabe, O. Reuven, T. Scaffidi, A. Stern, A. K. Geim, et al., Imaging hydrodynamic electrons flowing without landauer–sharvin resistance, Nature 609, 276 (2022).
  • Stern et al. [2022] A. Stern, T. Scaffidi, O. Reuven, C. Kumar, J. Birkbeck, and S. Ilani, How electron hydrodynamics can eliminate the landauer-sharvin resistance, Phys. Rev. Lett. 129, 157701 (2022).
  • Nagaev and Ayvazyan [2008] K. Nagaev and O. Ayvazyan, Effects of electron-electron scattering in wide ballistic microcontacts, Phys. Rev. Lett. 101, 216807 (2008).
  • Nagaev and Kostyuchenko [2010] K. Nagaev and T. Kostyuchenko, Electron-electron scattering and magnetoresistance of ballistic microcontacts, Phys. Rev. B 81, 125316 (2010).
  • Renard et al. [2008] V. Renard, O. Tkachenko, V. Tkachenko, T. Ota, N. Kumada, J.-C. Portal, and Y. Hirayama, Boundary-mediated electron-electron interactions in quantum point contacts, Phys. Rev. Lett. 100, 186801 (2008).
  • Guo et al. [2017] H. Guo, E. Ilseven, G. Falkovich, and L. S. Levitov, Higher-than-ballistic conduction of viscous electron flows, PNAS 114, 3068 (2017).
  • Keser et al. [2021] A. C. Keser, D. Q. Wang, O. Klochan, D. Y. Ho, O. A. Tkachenko, V. A. Tkachenko, D. Culcer, S. Adam, I. Farrer, D. A. Ritchie, et al., Geometric control of universal hydrodynamic flow in a two-dimensional electron fluid, Phys. Rev. X 11, 031030 (2021).
  • Poiseuille [1844] J. L. Poiseuille, Recherches expérimentales sur le mouvement des liquides dans les tubes de très-petits diamètres (Imprimerie Royale, 1844).
  • Batchelor [2000] G. Batchelor, An introduction to fluid dynamics (2000).
  • Landau and Lifshitz [2013] L. D. Landau and E. M. Lifshitz, Fluid mechanics: Landau And Lifshitz: course of theoretical physics, Volume 6, Vol. 6 (Elsevier, 2013).
  • Gooth et al. [2018] J. Gooth, F. Menges, N. Kumar, V. Süβ𝛽\betaitalic_β, C. Shekhar, Y. Sun, U. Drechsler, R. Zierold, C. Felser, and B. Gotsmann, Thermal and electrical signatures of a hydrodynamic electron fluid in tungsten diphosphide, Nat. Commun. 9, 4093 (2018).
  • Baker et al. [2024] G. Baker, T. W. Branch, J. Bobowski, J. Day, D. Valentinis, M. Oudah, P. McGuinness, S. Khim, P. Surówka, Y. Maeno, et al., Nonlocal electrodynamics in ultrapure pdcoo 2, Phys. Rev. X 14, 011018 (2024).
  • Aharon-Steinberg et al. [2022] A. Aharon-Steinberg, T. Völkl, A. Kaplan, A. K. Pariari, I. Roy, T. Holder, Y. Wolf, A. Y. Meltzer, Y. Myasoedov, M. E. Huber, et al., Direct observation of vortices in an electron fluid, Nature 607, 74 (2022).
  • Dean et al. [2010] C. R. Dean, A. F. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, T. Taniguchi, P. Kim, K. L. Shepard, et al., Boron nitride substrates for high-quality graphene electronics, Nature Nanotech. 5, 722 (2010).
  • Palm et al. [2024] M. L. Palm, C. Ding, W. S. Huxter, T. Taniguchi, K. Watanabe, and C. L. Degen, Observation of current whirlpools in graphene at room temperature, Science 384, 465 (2024).
  • Clericò et al. [2019] V. Clericò, J. A. Delgado Notario, M. Saiz-Bretín, A. V. Malyshev, Y. M. Meziani, P. Hidalgo, B. Méndez, M. Amado, F. Domínguez-Adame, and E. Diez, Quantum nanoconstrictions fabricated by cryo-etching in encapsulated graphene, Sci. Rep. 9, 13572 (2019).
  • Rhodes et al. [2019] D. Rhodes, S. H. Chae, R. Ribeiro-Palau, and J. Hone, Disorder in van der waals heterostructures of 2d materials, Nature Materials 18, 541 (2019).
  • Clericò et al. [2018] V. Clericò, J. A. Delgado Notario, M. Saiz-Bretín, C. Hernández Fuentevilla, A. V. Malyshev, J. D. Lejarreta, E. Diez, and F. Domínguez-Adame, Quantized electron transport through graphene nanoconstrictions, Phys. Status Solidi A 215, 1701065 (2018).
  • Estrada-Alvarez et al. [2024] J. Estrada-Alvarez, F. Dominguez-Adame, and E. Díaz, Alternative routes to electron hydrodynamics, Comm. Phys. 7 (2024).
  • Holder et al. [2019] T. Holder, R. Queiroz, T. Scaffidi, N. Silberstein, A. Rozen, J. A. Sulpizio, L. Ella, S. Ilani, and A. Stern, Ballistic and hydrodynamic magnetotransport in narrow channels, Phys. Rev. B 100, 245305 (2019).
  • Sandner et al. [2015] A. Sandner, T. Preis, C. Schell, P. Giudici, K. Watanabe, T. Taniguchi, D. Weiss, and J. Eroms, Ballistic transport in graphene antidot lattices, Nano Lett. 15, 8402 (2015).
  • Vaquero et al. [2023] D. Vaquero, V. Clericò, M. Schmitz, J. A. Delgado-Notario, A. Martín-Ramos, J. Salvador-Sánchez, C. S. Müller, K. Rubi, K. Watanabe, T. Taniguchi, et al., Phonon-mediated room-temperature quantum hall transport in graphene, Nat. Comm. 14, 318 (2023).
  • Hwang and Sarma [2008] E. Hwang and S. D. Sarma, Acoustic phonon scattering limited carrier mobility in two-dimensional extrinsic graphene, Phys. Rev. B 77, 115449 (2008).
  • Giuliani and Vignale [2008] G. Giuliani and G. Vignale, Quantum theory of the electron liquid (Cambridge university press, 2008).
  • Callaway [1959] J. Callaway, Model for lattice thermal conductivity at low temperatures, Phys. Rev. 113, 1046 (1959).
  • Ledwith et al. [2019] P. Ledwith, H. Guo, A. Shytov, and L. Levitov, Tomographic dynamics and scale-dependent viscosity in 2d electron systems, Phys. Rev. Lett. 123, 116601 (2019).
  • Kryhin et al. [2023] S. Kryhin, Q. Hong, and L. Levitov, t𝑡titalic_t-linear conductance in electron hydrodynamics, arXiv preprint arXiv:2310.08556  (2023).
  • Kiselev and Schmalian [2019] E. I. Kiselev and J. Schmalian, Boundary conditions of viscous electron flow, Phys. Rev. B 99, 035430 (2019).
  • Pezzini et al. [2012] S. Pezzini, C. Cobaleda, E. Diez, and V. Bellani, Quantum interference corrections to magnetoconductivity in graphene, Phys. Rev. B 85, 165451 (2012).
  • Eroms and Weiss [2009] J. Eroms and D. Weiss, Weak localization and transport gap in graphene antidot lattices, NJP 11, 095021 (2009).
  • Power et al. [2017] S. R. Power, M. R. Thomsen, A.-P. Jauho, and T. G. Pedersen, Electron trajectories and magnetotransport in nanopatterned graphene under commensurability conditions, Phys. Rev. B 96, 075425 (2017).
  • Yagi et al. [2015] R. Yagi, R. Sakakibara, R. Ebisuoka, J. Onishi, K. Watanabe, T. Taniguchi, and Y. Iye, Ballistic transport in graphene antidot lattices, Phys. Rev. B 92, 195406 (2015).
  • Afanasiev et al. [2021] A. Afanasiev, P. Alekseev, A. Greshnov, and M. Semina, Ballistic-hydrodynamic phase transition in flow of two-dimensional electrons, Phys. Rev. B 104, 195415 (2021).
  • Di Ventra [2008] M. Di Ventra, Electrical Transport in Nanoscale Systems (Cambridge University Press, 2008).
  • Whiteley [2017] J. Whiteley, Finite Element Methods, Mathematical Engineering (Springer International Publishing AG, 2017).
  • Engwirda and Ivers [2016] D. Engwirda and D. Ivers, Off-centre steiner points for delaunay-refinement on curved surfaces, Computer-Aided Design 72, 157 (2016).
  • Ginzburg et al. [2023] L. V. Ginzburg, Y. Wu, M. P. Röösli, P. R. Gomez, R. Garreis, C. Tong, V. Stará, C. Gold, K. Nazaryan, S. Kryhin, et al., Long distance electron-electron scattering detected with point contacts, Phys. Rev. Res. 5, 043088 (2023).