{'fraction_non_alphanumeric': 0.03173320402995528, 'fraction_numerical': 0.011098539949356176, 'mean_word_length': 4.54089552238806, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical evidence of electron hydrodynamic whirlpools in graphene samples

We present an extension of recent relativistic Lattice Boltzmann methods based on Gaussian quadratures for the study of fluids in (2+1) dimensions. The new method is applied to the analysis of electron flow in graphene samples subject to electrostatic drive; we show that the flow displays hydro-electronic whirlpools in accordance with recent analytical calculations as well as experimental results.


Introduction
Relativistic fluid dynamics has so far been generally confined to the study of astrophysical phenomena, and mostly for ideal non-viscous fluids. However, it has been remarked recently that a fluid dynamics approach may be able to capture interesting aspects of the behavior of systems at much smaller scales. A pioneering example is the study of the behavior of the ultra-relativistic quarkgluon plasma formed in the collision of high-energy heavy ions in particle accelerators [1,2,3]. More recently, it has been suggested that relativistic fluid dynamics is relevant to understand the behavior of quantum states that can now be studied in several condensed matter experimental setups [4], within the context of so-called AdS-CFT holographic fluids [5,6,7,8].
Recent experimental studies have shown that certain features of the flow of electrons in graphene can be explained through a pseudo-relativistic hydrodynamic approach [9], confirming earlier theoretical predictions [10,11].
In this paper, we present preliminary results on the development and validation of computationally efficient numerical approaches aimed at capturing details of the electron behavior in these systems, addressing the specific case of graphene.
In brief, electrons in graphene follow an "ultra-relativistic" dispersion relation, so they can be considered as a fluid of massless (quasi-)particles whose energy depends on the momentum as E = v f p, with v f ∼ 10 6 m/s the Fermi speed, mimicking the role of the speed of light in true relativistic systems.
The observation of hydrodynamic regimes is predicted to be simpler in doped graphene sheets [12], which are characterized by large viscosities.
The Lattice Boltzmann Method (LBM) is a class of computational fluid dynamics solver which has attracted much interest in the past three decades for the solution of the Navier-Stokes equations. The method stems from the kinetic level and consists of a quadrature-based discretization of the Boltzmann equation, allowing the use of a regular grid of points to exactly match the moments of an equilibrium distribution function up to a desired order. Rigorous mathematical analyses based on the Chapman Enskog expansion are typically employed to bridge between the kinetic and the macroscopic level for near equilibrium fluids. Recently, new formulations have been introduced to extend the method to the study of relativistic hydrodynamics. The relativistic Lattice Boltzmann Method (RLBM) [13,14] offers an appealing solution for the study of dissipative relativistic hydrodynamics, since viscosity is naturally included in its formulation and, furthermore, it does not involve second order derivatives in space, preserving relativistic invariance and causality by construction [15].
The interest towards the numerical study of electrons flow in graphene has motivated the development of two-dimensional RLBM solvers [16,17,18,19].
Most of these numerical methods are based on a second order expansion of an equilibrium distribution function following the Fermi-Dirac statistics, and they have been applied to study e.g. low-viscosity pre-turbulent regimes.
In [20], working in three dimensions, we have shown that third order expansions of the equilibrium distribution function are the minimum requirement to correctly handle dissipative effects in simulations of the relativistic regime.
In this paper, we present a new RLBM in two dimensions, based on a third order expansion of the equilibrium Maxwell-Jüttner distribution. Quantum effects are not described in this model, a choice which simplifies the algorithmic derivation allowing us to retain one of the main LBM features, namely perfect streaming. This could be regarded as the first step in the derivation of a truly accurate model for the fluid dynamics descriptions of electrons in graphene, but we also expect that quantum effects should have a limited impact on the averages involved in hydrodynamical bulk observables. As a result, we expect that the present model should be able to provide new useful insights into the physics of relativistic electron flow in graphene devices. We provide a first validation test of our approach simulating a doped single layer graphene sheet in the so-called "vicinity-geometry", which was considered in a series of papers [21,12,22] to outline phenomena such as negative nonlocal resistance and current whirlpools. The numerical method is tested in a steady-state regime, for which semi-analytical solutions are available, showing satisfactory agreement with previous works; we wish to emphasize that the present numerical method allows to describe time-dependent, nonlinear flows which escape analytical treatments. Hence, the future plan is to address electron flows of experimental interest. This paper is organized as follows. In Section 2 we describe our new RLBM model, summarizing the general procedure used in its derivation and providing details of the quadrature and of the external forcing scheme. In Section 3 we carry out simulations of the "vicinity-geometry" replicating the formation of current whirlpools and providing quantitative comparisons of the electrochemical potential against analytical approximations available in the literature.

Model Description
This section briefly describes our RLBM model, which is a "downsizing" to two spatial dimensions of a similar model handling relativistic flows in 3D; for a more detailed description, the reader is referred to [23,20].

Relativistic Boltzmann equation
We consider an ideal non-degenerate relativistic fluid, consisting at the kinetic level of a system of interacting particles of rest mass m. The particle distribution function f ((x α ), (p α )), depending on space-time coordinates (x α ) = (ct, x) and momenta (p α ) = p 0 , p = p 2 + m 2 , p (c is the speed of light, x, p ∈ R 2 ), describes the probability of finding a particle with momentum p at a given time t and position x. We adopt Einstein's summation convention over repeated indexes, and use Greek indexes to denote (2 + 1) space-time coordinates and Latin indexes for 2 dimensional spatial coordinates. The particle distribution function obeys the relativistic Boltzmann equation, here taken in the Anderson-Witting [24,25] relaxation-time approximation: with τ the relaxation (proper-)time, (U α ) = γ · (c, u) the macroscopic (2 + 1)- , K α the external forces acting on the system (for simplicity we assume they do not depend on the momentum), and f eq the local equilibrium. In this work f eq will follow a Maxwell-Jüttner distribution: where A is a normalization constant and k B the Boltzmann constant.
At the macroscopic level the Anderson-Witting model correctly reproduce the conservation equations, i.e. ∂ α N α = 0 and ∂ β T αβ = 0, with N α the particle (2+1)-flow and T αβ the energy-momentum tensor. At equilibrium N α and T αβ can be described by the moments of the equilibrium distribution function: T αβ E = f eq p α p β dp where n is the particle number-density, the energy density, P the pressure and η αβ the Minkowski metric tensor. To be noted that the normalization constant in Eq. 2 has to be chosen in such a way to satisfy the relation with the fluid particle density in Eq. 3. In the following we will use η αβ = diag(1, −1, −1), and adopt natural units for which c = k B = 1 .

Lattice discretization
In this section we revise the general procedure, used in the derivation of previous non-relativistic [26,27,28,29] and relativistic LBMs [30,31,23], for the discretization of the Boltzmann equation on a lattice.
We start from an expansion of the equilibrium distribution function f eq in a basis of polynomials, orthogonal with respect to a weighting function ω corresponding to f eq in the fluid rest frame (where u = 0). It is simple to verify that in the rest frame Eq. 2 reduces to where the normalization factor N R is taken such that ω(p 0 ) dp/p 0 = 1. Starting from the basis V = {1, p α , p α p β , . . . } one derives the set of orthogonal polynomials {J (i) , i = 1, 2 . . . } by following a Gram-Schmidt procedure, with the inner product defined using the weighting function in Eq. 5. In Appendix A we provide an example of polynomials up to the third order for m = 0. The polynomials are then used to build the expansion: where a (k) are the projection coefficients defined as Observe that by construction the coefficients a (k) coincide with the moments of the distribution function; this is a crucial aspect since it follows that f eq N , obtained truncating the summation in Eq. 6 such to include only the terms of order up to N , correctly preserves the moments of the distribution up to the N − th order. The next step consists in determining a Gauss-type quadrature on a Cartesian grid, with the aim of i) ensuring exact streaming by requiring that all quadrature points lie on lattice sites ii) preserving the moments of the distribution up to a desired order N . The discretized version of the equilibrium distribution can be then written as follows: where w i and p µ i are the weights and the nodes of the quadrature, respectively. The analytic expression of f eq iN for N = 3 and m = 0 is given in Appendix C. At this stage it is possible to formulate the discrete Boltzmann equation, which in the relativistic case reads as A detailed description of the algorithmic derivation for the 3-dimensional case is given in [23]. The algebraic complexities in the calculation of the polynomials and the expansion of the equilibrium distribution significantly simplify in 2-D. The full details will be described at length in a future expanded version of this work.

Quadrature with prescribed nodes
As discussed in the introduction, here we focus our attention on solving Eq. 10 using polynomials up to the third order.
The lattice discretization of the Boltzmann equation can be reduced to a quadrature problem. In practice, one needs to find the weights and the abscissas of a quadrature able to satisfy the orthonormal conditions up to a desired order [32]: here p µ i are the discrete quadri-momentum vectors. A convenient parametrization of p µ i was given in [23] and writes as follows: where n i ∈ Z 2 are the vectors forming the stencil G = {n i | i = 1, 2, . . . , i max } defined by the (on-lattice) quadrature points, v 0 is a free parameter that can be freely chosen such that v i = v 0 ||n i || ≤ 1, ∀i, and p 0 i is defined as In order to determine a quadrature we proceed as follows: i) select a value for the rest massm = m/T 0 (with T 0 a reference temperature on the lattice), ii) choose a set of velocity vectors G, formed by a sufficient number of elements such that the left hand side of Eq. 10 is a full ranked matrix, iii) look for a solution of Eq. 10 formed by non-negative weights (w i ≥ 0, ∀i).
We point out that the parametrization in Eq. 11 is general and can be used to determine quadratures for wide ranges of values ofm.
As an example we show in Figure 1a a set of vectors that can be used to build a quadrature form = 5. In the remainder of this paper we are interested in particular in the case of massless particles, all traveling at the same speed v i = c = 1, ∀i. Since for m = 0 Eq. 12 is not well defined, we let p 0 i be free parameters (as already suggested in [31]) to be determined such as to satisfy Eq. 10. We can have several energy shells associated to each vector and therefore we add a second index to Eq. 11:

Forcing Scheme
The definition of force in relativity is subject to a certain degree of arbitrariness due to the lack of certain general properties such as, for example, Newton's third law [33]. In the following we will use the definition of the Minkowski force: subject to the condition and To introduce a forcing term in our numerical scheme we make the following two assumptions: i) the force does not depend on the momentum three vector ( ∂K α ∂p α = 0 ) ii) the distribution function in not far from the equilibrium, such that we can approximate the term K α ∂f ∂p α in Eq. 1 with an expansion that uses the same polynomials used for the equilibrium distribution function: with the projection coefficients defined as

Numerical Tests
We now apply the model described in the previous section to the simulation of the (pseudo)-relativistic dynamics of electrons in graphene sheets; as already remarked, in this case the Fermi velocity v f of the simulated system plays the role of the speed of light. We consider an experimental setup consisting of an ultraclean single layer graphene encapsulated between boron nitride crystals in which it has been shown that electrons exhibit a hydrodynamic flow [12]. This setup has been used in a series of works [21,12,22] to highlight peculiar properties such as negative nonlocal resistance and current whirlpools. The geometry is sketched in Figure 2. The total force acting on the system is given by the vector sum of the force due to the electric field F E and the force due to the pressure gradient F P . While F P is naturally described by our RLBM solver, F E is included in the form of an external force. Therefore, in the simulations the external force F (Eq. 16) is given by a self-consistent electric field E = −ρ e ∇φ, with ρ e = ne being the electron charge density. For our initial validation tests on this specific setup, we follow [34] and do not solve explicitly the Poisson equation for the electric potential, but rather use a local capacitance approximation defined as: where C g is the capacitance per unit area.
Using this setup, we simulate a system similar to the one considered in [21,12], where analytical results are obtained in the approximation of an infinitely long channel; we use a lattice with an aspect ratio L/W = 4, that we simulate on a lattice of 2000 × 500 grid points. The translation between physical units and adimensional lattice units is based on the definition of a length-unit on the lattice such that the width of the channel corresponds to the physical value and on an energy unit that we chose as the Fermi-energy of the simulated system.
In Figure 3 we show a snapshot of a simulation, using a constant initial density and a large value for the shear viscosity. As we can see, results are qualitatively comparable with those presented in [21,12]. In particular one can appreciate the (symmetric) formation of electron back-flows in the proximity of the gates, so called current whirlpools. For a more quantitative comparison, we take into consideration the electrochemical potential in the proximity of the injector, for which the following approximate analytic expression was derived in [21]: I is the driving current at the inlet, η is the shear viscosity,n is the equilibrium density, e is the electron charge, r and θ are used to parametrize in polar coordinates the space in the proximity of the inlet. In Figure 4 we compare the prediction of Eq. 20 with the results of our simulations by plotting the electrochemical potential as a function of the polar angle for several lattice points at several distances r from the center of the injector. In particular we show that for different setups, the quantity r 2 φ(r, θ) does not depend on r as predicted by Equation 20: to a good approximation, all curves collapse on the top of each other, as expected.
As a further benchmark we evaluate how the steady state solution reported in Figure 3 varies when tuning the magnitude of the driving forces F E and F P .
To this purpose, we perform simulations with different values of the parameter C g (see Eq. 19) to evaluate the role of the electric potential. Following [35]  one would not expect to observe the effect of Coulomb interactions for static flows. In Figure 5 we show that this is indeed the case; in fact varying C g over several different orders of magnitude does not yield any appreciable effect on the solution. Moreover the results are the same even in the case when F E is neglected (C G = ∞), proving that the model gives a self-consistent description of hydrodynamic theory on long length scales. On the other hand, the electric potential is expected to play a major role on the dynamics of non-linear, timedependent flows, which will make the object of forthcoming studies.

Conclusions
In this work we have described a new solver for the study of (2+1)-dimensional  view we have successfully reproduced the current whirlpools highlighted by recent experimental works [12]. Besides, we have provided a more quantitative validation, with a comparison of the electrochemical potential in the proximity of the current injector against previous analytic predictions [21]. We consider this to be a first step in the derivation of an accurate model for the study of the hydrodynamics behavior of electrons flow in graphene. Future works will deal with more robust comparisons of simulations against experimental data and more detailed simulations of actual experimental setups. This work will allow a proper evaluation of the loss of accuracy due to the neglect of quantum effects, alongside with further parameters that should be taken into account (such as electrons collisions with impurities and phonons) to expand the capabilities of the present model.

B. Third order 2D Orthogonal Projections
In this appendix we provide the analytic expressions of the orthogonal projections a (k) , up to the third order, for the ultra-relativistic case. The notation follows the one introduced in Appendix A for the orthogonal polynomials. All the projections are scaled with the particle density n, thereby ensuring the correct normalization of the equilibrium distribution function (Eq. 2).

C. Third order expansion of the equilibrium distribution function
The third order expansion of the Maxwell-Jüttner distribution in two dimension and for m = 0, which allows to recover the first, the second and the third order moments of Eq. 2, was derived using the polynomials defined in Appendix A and projections in Appendix B. It reads as follow: −1) − 4T u 0 + 2 + (p 0 i − 1)(T u 0 − 1) + p x i p y i T 2 u x u y + p x i T u x + p y i T u y + 1