Weyl–Wigner description of massless Dirac plasmas: ab initio quantum plasmonics for monolayer graphene

We derive, from first principles and using the Weyl–Wigner formalism, a fully quantum kinetic model describing the dynamics in phase space of Dirac electrons in single-layer graphene. In the limit ℏ → 0, we recover the well-known semiclassical Boltzmann equation, widely used in graphene plasmonics. The polarizability function is calculated and, as a benchmark, we retrieve the result based on the random-phase approximation. By keeping all orders in ℏ, we use the newly derived kinetic equation to construct a fluid model for macroscopic variables written in the pseudospin space. As we show, the novel ℏ-dependent terms can be written as corrections to the average current and pressure tensor. Upon linearization of the fluid equations, we obtain a quantum correction to the plasmon dispersion relation, of order ℏ 2, akin to the Bohm term of quantum plasmas. In addition, the average variables provide a way to examine the value of the effective hydrodynamic mass of the carriers. For the latter, we find a relation in which Drude’s mass is multiplied by the square of a velocity-dependent, Lorentz-like factor, with the speed of light replaced by the Fermi velocity, a feature stemming from the quasi-relativistic nature of the Dirac fermions.


Introduction
Since its experimental realization in 2004 [1], graphene has gather a considerable amount of attention due to its outstanding optical, electronic and mechanical properties, being the main topic of research in the field of two-dimensional (2D) materials. Since then, many different approaches have been applied to study a variety of phenomena in graphene, ranging from electronic transport and hydrodynamics [2][3][4][5], to quantum Hall effect and topology [6][7][8][9]. On the top of all its remarking properties is the fact that the low-energy electrons in graphene behave like massless Dirac fermions, traveling at a constant (Fermi) velocity of v F 10 6 m s −1 . Such behavior is a consequence of graphene's 2D crystal structure, made out of carbon atoms arranged in a hexagonal lattice, which comprises two interpenetrating triangular Bravais lattices [10,11]. The two atoms per unit cell give rise to two independent solutions, labelled by a band index s (often called pseudospin) representing the conduction (s = +1) and valence (s = −1) bands. This additional degree of freedom grants the carriers a chiral nature and provides an additional playground to test quantum electrodynamics [12].
Its unique properties turn graphene into a prominent candidate for constructing transparent electronic devices, ultra-sensitive photodetectors, and other high-performance optoelectronic structures [13][14][15]. Furthermore, the collective oscillations of the electron fluid lead to the formation of plasmons (or plasma waves) [16]. In fact, plasmonic of 2D materials, such as graphene, transition metal dichalcogenides, and hexagonal boron nitride [17,18], is nowadays a very celebrated field of research [19]. Graphene-based plasmonics finds a variety of applications, as the versatility of graphene enables the manufacture of optical devices working in different frequency ranges, namely in the terahertz (THz) and the infra-red domains [20]. Thus, the investigation of graphene dynamics deserves special attention not only because of the theoretical interest that lies in the chiral and massless nature of its carriers, but also on account of its tremendous potential to produce novel electronics.
From the theoretical point of view, a variety of techniques have been developed to establish the dynamics of Dirac fermions in graphene, alternating from semiclassical hydrodynamical models [21][22][23][24], to quantum formulations that include collective Green's functions [25], time-dependent Hartree-Fock models [26,27], or time-dependent density functional theory [28,29]. Whilst it is often the case that cumbersome equations, of very reduced utility, crop up when going for a complete quantum description, it is also true that the semiclassical approach, based on the Boltzmann equation, is inadequate in the low temperature and high density regimes [30]. As it has been shown [31], the behavior of graphene carriers has a purely quantum character, so a quantum approach is clearly the most indicated one to start with. Although the semiclassical models treat the collisional term in the kinetic equation from a quantum mechanical point of view, the collisionless part of the said equation is inherited from the classical Boltzmann equation, after mere replacement of the parabolic velocity by the Dirac velocity. As we show below, such a procedure is not enough to correctly describe the collection of Dirac particles, and extra contributions do arise in the very collisionless part of the kinetic equation when the latter is properly derived from quantum principles.
The intention of this paper is to fill the gap between the previous semiclassical models and the more advanced quantum calculations, which is done by establishing a kinetic formalism based on the Weyl-Wigner formulation of quantum mechanics [32,33]. We chose this technique for two main reasons: firstly, it takes into account rigorous quantum-mechanical considerations; secondly, it provides a clear path to recover the results based on the more standard distribution theory, upon which the semiclassical models are constructed. We stress that, although previous works have already been focused on the Wigner representation of graphene [34,35], the approach followed by their authors concerns solely the interaction with an external field, without treating the electron-electron interaction, thus making them unsuitable for graphene plasmonics. More importantly, the results therein do not provide a clear connection with the semiclassical models that have been put forward and worked out extensively, something our own use of the Weyl-Wigner approach does allow for. Note that a similar method has already been applied to bilayer graphene as well [36], but this constitutes a significantly simpler system, as it is comprised of particles with parabolic dispersion relation, in which case the collisionless part of the kinetic equation indeed reduces to that of the classical Boltzmann equation.
As for the results, we obtain the plasmon dispersion relation after recovering the usual result based on the collisionless random-phase approximation (RPA), thus showing that our kinetic equation correctly predicts the polarizability function for graphene. In addition, after establishing the relevant macroscopic variables, hydrodynamical equations are obtained allowing for a fluid description of the system. We show that an infinite number of purely quantum contributions arise as corrections to the current density and pressure terms, due to the linear dispersion of the electrons. These corrections would be absent in the case of massive particles, which reinforces the need for an ab initio quantum formulation for the case of massless fermions. Then, we show that the classical limit applied to the fluid model appropriately recovers the semiclassical equations present in the literature. Further ahead, we are able to derive as well the effective mass of the fluid element and relate it to Drude's mass, thus contributing to the understanding of a still-open question in graphene hydrodynamics.
The paper is organized as follows: in section 2, we review the basic properties of graphene carriers and corresponding low-energy Hamiltonian. In section 3, we derive the quantum kinetic equation based on the Weyl-Wigner approach for the Dirac plasma. In section 4, we obtain the fluid equations, valid at the macroscopic scale, by taking the moments out of the kinetic equation, and discriminate its classical and quantum contributions. Finally, in section 5, some conclusions are drawn and future perspectives are outlined.

Graphene preliminaries
The electronic dynamics can be captured starting with a general form for the free Hamiltonian, written in a second quantized fashion as where R and R run over the graphene lattice, and l and l over the two sublattices (A and B). Additionally, u † l (R) and u l (R) denote the creation and annihilation operators, respectively, for each lattice point R and sublattice l. Moreover, the tight-binding approximation can be settled with a proper restriction on the matrix elements u l , R| H| u l , R . By allowing hopping only between nearest neighbors, we can set all matrix elements to zero with the exception of the cases of u l , R| H| u l , R + η i = −t, with t 2.97 eV the hopping integral and η i the nearest-neighbor vectors [10]. Resorting to the relation u l (R) = k u lk e ik·R / √ N, where N is the total number of carbon atoms, equation (1) reduces to where ϕ k = ( u Ak , u Bk ) T and Δ = i e −ik·η i . By working out the single-particle approximation, we can concentrate on a fixed wave-vector in the above summation. As such, the effective free Hamiltonian that we adopt comes after expanding equation (2) around the Dirac point K = (4π/3 √ 3d, 0), keeping only the linear terms and project it onto a fixed k. We obtain with k = ( k x , k y ) the 2D wave-vector operator, σ = (σ x , σ y ) the 2D vector of Pauli matrices, v F = 3ta/(2 ) 10 6 m s −1 the Fermi velocity, and a 2.6 Å the lattice parameter. It is important to note that equation (3) is appropriate to study the transport properties in graphene, specially in the semiclassical limit, but does not capture the full collective behavior. However, for the purposes of the current work, it suffices to remain within this level of approximation. After diagonalizing equation (3), we find its two eigen-pseudospinors where s = ±1 labels the conduction (+1) and valence (−1) bands and φ k = arctan(k y /k x ). We have where V(t) is a diagonal potential on pseudospin space. We determine the potential by imposing its position representation, V(r, t) . = r| V(t)|r , to verify the Poisson equation, where n is the total electronic density, e is the elementary charge and ε is the medium permittivity. A solution is found to be where r is the position operator. This particular form for V(t) is consistent with the initial single-particle approximation as it corresponds to the Hartree potential, thus discarding the contribution of the Fock and correlation terms. It is valid for small values of the coupling parameter r s ∼ α s /ε r (ratio of the average potential energy to the average kinetic energy), with α s = e 2 /(4πε 0 v F ) 2.2 the graphene structure constant, ε 0 the vacuum permittivity and ε r = ε/ε 0 . As a result, the Hartree approximation remains reliable as long as ε r 2.2.
To proceed, we perform our calculations with the density matrix ρ, given that it has a connection with the Wigner function, to be defined in the next section. The density matrix evolves in time according to the quantum Liouville equation [37], which reads with the right-hand side (rhs) representing the collision operator. The density in equation (7) can then be written as the trace and we recast equation (8) in a more convenient form, where we have set O ss

Weyl-Wigner formalism
A very handy way to treat the electronic system of equation (10) is to use the Weyl-Wigner picture of quantum mechanics [32,33], which allows for description in the phase-space, in close analogy with the classical case. In the classical limit, the Wigner function denotes the probability density of finding a particle in a given infinitesimal phase-space volume dr dk centred in (r, k). In the quantum case, due to the commutation relation between r and k, the Heisenberg uncertainty principle prevents particles to localise in a specific phase-space point, and a proper distribution function (i.e. non-negative everywhere) is not possible [38]. However, we can still construct the much renowned Wigner function W(r, k, t), whose properties are similar to those of a classical distribution function, i.e. it is a function of both a spatial coordinate r and a wave-vector k. Although r and k are classical-like variables and not conjugate quantum mechanical operators, they both give information about the spatial and momentum distributions of the system. Consequently, we may have W(r, k, t) < 0, with the phase-space regions where it takes negative values being purely quantum and having no classical analogue. For this reason, the Wigner function is often referred to as a quasiprobability function. The mathematical definition of the Wigner function is given by the Weyl transform of the density operator [32], The expectation value of operators can be computed by integrating the corresponding Weyl transforms multiplied by W(r, k, t) over the phase space, very much like in the classical case [39]. Introducing the completeness relation for Dirac states |ks , the Wigner function becomes from which we immediately get the Fourier transform We stress the fact that W(q, k, t) carries the dynamical information contained in the matrix elements of ρ, plus the lattice information contained in the overlap of the pseudospinors. In the case of pure plane waves (e.g. a single-band massive plasma, E k ∼ k 2 ), we would have equation (13) with v ss k,k = 1. By defining the matrix element W ss (q, k, t) When written in this manner, it becomes evident that the Wigner function correctly incorporates the dynamics of each band plus the contribution of the mixing terms, thus accounting for the intra-and inter-band excitations. It is indeed known that the introduction of the spin (or pseudospin) makes us move to a tensorial formulation [40]. One can show that the total density, following equation (9), has a simple form in terms of the Wigner function n(r, t) where n ss (r, t) . = A −1 k W ss (r, k, t). As discussed before, the diagonal elements correspond to the density of electrons in the conduction (s = +1) and valence (s = −1) bands, whereas the off-diagonal elements represent the coherence between both populations, i.e. the fraction of the electronic population that is constantly changing between the conduction and valence bands.
Next, we derive an equation of motion for the Wigner elements in momentum space. Upon time differentiating W ss and using equation (10), we are lead to where ΔE ss = Av ss k,k L ss k,k is the collision integral components, and the kernel C{W ss } is defined as Moreover, we have used with V(k, t) the Fourier transform of V(r, t). By Fourier transforming equation (7), we find accounting for the (2D) time-independent Coulomb potential.

Collective oscillations and the RPA
We now consider small perturbations around an equilibrium configuration, keeping the lowest order contributions to the Wigner components. As such, we expand the density matrix as T the absolute temperature, μ the chemical potential and β = (k B T) −1 . Besides, g s g v accounts for the spin (g s = 2) and valley (g v = 2) degeneracy. The former results from the degeneracy of the spin populations in each energy band, which we have neglected in our treatment so far, and the latter should be incorporated to consistently include the two minima in the first Brillouin zone [41]. The perturbation δρ(t) is assumed to be small, δρ(t) ∼ O( 4 ), such that higher powers in the coupling parameter are discarded. Consequently, the Wigner elements give with δW ss (q, k, t) a small perturbation. Similarly, the density and potential are perturbed as where is the equilibrium density and δn(q, t) a small perturbation verifying |δn(q, t)| 1. Introducing equations (19)-(21) into equation (16), and Fourier transforming the latter in the time domain, yields where F ss = |v ss k,k | 2 is the chirality factor, and a second-order term has been neglected, as well as the collision term. Equation (22) is formally equivalent to Kubo's formula for the linear response of a many-body system [42], and reproduces the features contained in the RPA for the collisionless limit S ss k,k → 0 [31,43]. This formalism is specially advantageous to describe the dynamics of electrons that are far from equilibrium, such as the case of plasma instabilities, with the configuration being solely defined by the equilibrium distributions. Upon summing both sides of equation (22) over k, s and s , we find the plasmon dispersion to be given by (q, ω) = 0, with (q, ω) = 1 + U(q)Π(q, ω) the dielectric function and Π(q, ω) the polarizability, The formal result of equation (23) has been found before, in the context of the RPA [44]. It applies in the weak-coupling limit r s → 0, where the plasmon frequency can be computed using the noninteracting irreducible polarizability.
In what follows, we consider the case of negatively doped graphene, with the conduction band filled up to the Fermi level E F = v F k F μ 0. The Fermi level defines the Fermi wave-number k F , and the latter is related to the doping density n 0 by Typical experimental values of n 0 in the range 10 9 -5 × 10 12 cm −2 are achievable in graphene [45]. In the case E F k B T, the presence of holes (i.e. vacancies of valence electrons) is negligible, which allows us to set = 0 for all q and calculate the zeros of equation (23) explicitly (notice that the off-diagonal terms s = s in the polarizability are also not important since we are interested in the long-wavelength limit, for which we expand F ss k,k+q = δ s,s + O(q 2 ); simultaneously, the zeroth order in q of the remaining expression vanishes). For the sake of simplicity, we use the ultra-cold limit (T → 0) of the Fermi-Dirac distribution, with Θ(x) the Heaviside step function. After expanding equation (23) around q = 0 and keeping only terms up to q 2 , the plasmon dispersion relation is obtained where ω p is the characteristic plasmon frequency, For the characteristic experimental values ε r = 2.5 and n 0 within 5 × 10 9 -10 12 cm −2 , ω p lies in the THz region, between 2.6-37.3 THz. The first term ω ∼ √ q describes the long wavelength signature of plasmons in 2D electron gases [43,46,47]. The most notable difference, when compared to the characteristic plasmon frequency in the 3D parabolic case, ω 3D p = e 2 n 0 /(εm), is the appearance of in leading order, revealing its pure quantum nature. Therefore, no classical counterpart exists for the 2D Dirac plasma.
The same kinetic approach can be used to describe graphene electrons in a field effect transistor structure, i.e. placed between two metallic contacts, source and drain, and controlled by a gate. The gate voltage is related to the carrier density by [48] U(r, t) = en(r, t) where C g and C q are, respectively, the gate and quantum capacitance. The gate capacitance is given by C g = εd 0 , where d 0 is the gate separation. The quantum capacitance C q reflects the change in the potential with the band occupancy and is defined as is the density of states. For typical carrier densities we have C g C q , and so the second term in equation (29) can be neglected. To include this effect, we replace the effective Coulomb potential by −eU(r, t), which amounts to the substitution U(q) → e 2 d 0 /ε. The dispersion relation becomes and S the sound velocity of the electron fluid [3]. It is patent that the gate induces the linear behavior ω ∼ q, as a consequence of the (static) screening produced by the gate potential. Equation (30) is plotted in figure 1, alongside with equation (27) for comparison.

Classical limit and the Boltzmann equation
In this section, we adopt an approximation for equation (16) that simplifies the ensuing calculations. By realizing that V(q, t) ∼ 1/q, it is convenient to evaluate the overlap factors on the rhs at q 0, since the integrand vanishes for very large q . Thus, equation (16) will be approximated by For the diagonal components W s .
= W ss (r, k, t), the real-space version of equation (31) reads where W s ± .
= W s (r, k ± q/2, t), S ss (r, k, t) is the Fourier transform of S ss (t) and represents the quantum kinetic operator. For parabolic systems, with dispersion relation E(k) = 2 k 2 /(2m 2 * ), we would have the local kinetic operator K{W} = ( k/m * ) · ∇W, which would lead to the well-known Wigner-Moyal equation describing, e.g. quantum plasmas in metallic or semi-conductor systems [49][50][51][52][53][54][55][56]. Equation (32) reads as one of the important results of this paper, for it is the quantum analogue of the well-known Boltzmann equation for a Dirac system. We are led to the conclusion that, in order to construct the correct quantum kinetic equation for graphene, it is not enough to infuse, say, some 'quantumness' into the collision integral, either by using a quantum cross-section or taking into account proper statistics. In fact, the very structure of the collisionless, i.e. lhs of equation (32) is fully modified [60], which is not taken into account in the former works based on the semiclassical limit of that same equation (32) [21][22][23][24]. This is yet another consequence of the unconventional nature of Dirac particles.
In the classical limit → 0, we indeed find with v s k .
= sv F k/k the single-particle velocity. Following the same limit, we have with p = k the momentum, which implies Under these approximations, we finally recover the Boltzmann equation, consistently with most 'semiclassical' equations present in the literature (whose lhs is, in fact, classical and only the rhs is treated quantumly). We draw attention to the fact that the Weyl-Wigner formulation of quantum mechanics could also be derived by applying quantum deformation theory to the classical Poisson algebra [60]. Equation (32), despite its algebraical complexity, is possible to treat numerically, by adapting the approach of references [57][58][59], which include recipes for integrating the finite differences in the Wigner function. As for the non-local kinetic term of equation (33), a standard Fourier transform method can be applied.

Fluid model
One of the major advantages of the present description is the possibility of computing the average values of operators, which are naturally given in terms of the Wigner function. Their time evolution can be computed using equation (31), and a set of coupled equations can be put forward. As an example, we calculate the average momentum density P(r, t) = tr P ρ(t) , where is the momentum density operator, p = k, and {, } denotes symmetrisation. It follows that with the tensor P ss containing both the contribution of the particles in each band (s = s ) and of the particles which are in mixing states (s = s ). It is natural to define the average momentum p ss through the relation n ss p ss = P ss , thus setting the relevant macroscopic variables to be and p ss (r, t) = 1 n ss dk (2π) 2 k W ss (r, k, t), (42) after replacing the sums with integrals. Henceforth, and in order to simplify our analysis, we neglect the off-diagonal contributions, in as much as we are interested in the semiclassical limit of negatively doped graphene, for which case they play a minimal role. We are, thus, allowed to drop the double index to facilitate the notation (as, e.g. W s . = W ss ). The time evolution equations for the diagonal quantities in equations (41) and (42) provide and where j s = n s v s + j s Q is the density current, comprising the usual classical term plus a quantum ( -dependent) part stemming from the non-vanishing p-derivatives of the Dirac velocity, v s p .
= sv F p/p. Above, the arrows indicate the terms in which the derivatives act. The pressure tensor takes the form P s = P s cl + P s Q , where is the 2D (classical) pressure, and P s Q introduces quantum-mechanical contributions, Above, [v s p ] m refers to the mth component of v s p . By taking the limit → 0, it follows that j s → n s v s and P s → P s cl , as expected for the classical case. The classical limit of equations (43) and (44) would thus be recovered if one had replaced the Wigner equation by the Boltzmann equation from the start. Note that the new terms arise after taking into account not only the local contribution (∼v s k ), but all (non-local) terms of equation (33). To the best of our knowledge, equations (46) and (48) have not been given anywhere else. In the case of systems with parabolic dispersion, j s Q and P s Q would vanish and the fluid equations derived from the Wigner model would present no advantage over its classical analogue. This can be deduced from equations (46) and (48) by simply assuming the dependence v s k ∼ k. A reason for this particular behavior relies on the purely quantum structure of equation (3), in contrast to the kinetic term for massive particles, which takes the same form irrespective of treating the classical or the quantum case. On the contrary, the Dirac kinetic operator has no classical analogue.
By including the first terms of j s Q and P s Q , we obtain the modified fluid equations and The tensors J s ijl and T s ijl result from the first (∼ 2 ) terms of j s Q and P s Q , and read and Furthermore, G s (r, t) denotes the average value of some generic function G s (p), akin to equation (45). We will provide explicit expressions for these quantities below and discuss how they modify the plasmon dispersion relation.

Calculation of the macroscopic variables
Having set the relevant transport equations, we move now to a more detailed discussion concerning the averaged momentum and velocity fields. Being interested in the case of full valence band and doped conduction band (μ k B T), we concentrate on the conduction electrons only, thus dropping the band index (s = 1). As one can readily observe by comparing equations (42) and (45), a proportionality relation of the form p = mv does not hold for massless particles, with a constant m [4]. However, by allowing a space and time dependence on the mass, we are able to define an effective 'mass tensor', m ij . = m ij (r, t), as The meaning of such field should be clear: although the carriers have no mass, since they are Dirac particles, the fluid velocity and momentum fields can be related via an effective mass, providing a measure of the inertia of a fluid element. A tensorial structure such as that of equation (53) is required since, at a first glance, we should anticipate the components of p in a given direction i to be given by a linear combination of the velocity components in each direction j, with coefficients m ij . Nevertheless, if one is dealing with a rotationally invariant system, the mass tensor becomes diagonal, m ij = δ ij m.
Let us now modify the equilibrium used in the last section in order to accommodate for a (possibly local) average velocity of the electronic fluid, namely, where μ . = μ(r, t), u . = u(r, t) and β . = 1/[k B T(r, t)] are local variables to be calculated by imposing the conservation of number of particles, momentum and energy, respectively [61]. Additionally, this particular choice for the Wigner function cancels the collision term on the Wigner equation [22,62]. Computations with equation (54) are valid under the assumptions that the macroscopic quantities are slow-varying in both space and time, which is the natural requirement to go from a microscopic field theory to a macroscopic fluid model. The fluid description is valid if changes in the macroscopic quantities take place on large spacial and temporal scales, i.e. if q k F and ω ω p , respectively. This requirement is fulfilled if the characteristic time of the kinematic processes, t ∼ 1/ω, is much longer than the inverse collision frequency, 1/ν c , and the typical length L is much greater than the mean free path l ∼ v F /ν c , so that the plasma can be regarded locally as in a quasi-equilibrium configuration. However, due to the quantum nature of the model, we were able to capture the relativistic Dirac structure in a rigorous way. Actually, in the context of 2D quantum plasmas achieved in semiconductor structures, the De Broglie wavelength is replaced by the Thomas-Fermi screening length λ TF = 2π/(g s g v r s k F ). This is the analogue of the classical Debye length in plasmas, and differs from the Fermi wavelength, λ F = 2π/k F (∼10 nm − 10 μm for typical graphene parameters). Consequently, the fluid limit is expected to be valid provided the conditions qλ TF 1, qλ F 1 (55) apply. For graphene, we have λ F /λ TF 8.8/ε r , such that the second condition is the most stringent. Typical values of k F are found between 10 3 and 10 6 cm −1 . At smaller wavelengths, the microscopic structure becomes important and the fluid approximation no longer holds. As expected, equation (45) provides v = u. In its turn, equations (41) and (42) yield and where ξ(r, t) .
= v(r, t)/v F is the reduced velocity, n 0 is the density at rest, and is the local mass. Li n denotes the polylogarithm function of order n. Both m and n 0 admit amenable forms in the limit of strong dopping (μ k B T), for which case we find n 0 πk 2 F and with M = k F /v F being Drude's mass. Note that this result differs from that which can be found in the literature [22]. 2 The classical pressure reduces to P cl (r, t) = P 0 (f 1 with and In the above expressions, Ω = arctan(v y /v x ) is the polar angle of the fluid velocity and is the pressure in the absence of flow. The average value of the quantum tensor J ijl vanishes under the assumption of equation (54), i.e. J ijl = 0 for any {ijl}, which yields the density current to be approximated by its classical counterpart up to third order in , j = nv + O 4 . On the contrary, we find a non-zero contribution from the average value of T ijl = (T x ijl , T y ijl ) = (p x J ijl , p y J ijl ). Given the symmetry properties of equations (51) and (52), the independent components are found to be dispersion relation in the hydrodynamical regime, The first two terms enclose the classical behavior of the plasma waves [63,64]. The third one accounts for a quantum correction (∼ 2 ) known as the Bohm term [65], which is modified by a small numerical factor. The Bohm potential plays the role of a quantum pressure and here results from the contribution of T x xxx . For Dirac particles, we find a negative sign in the quantum correction (i.e. the Bohm pressure softens the plasmon mode), contrary to the what has been reported for the case of parabolic fermions, in which case it takes the form 2 q 4 /(4m 2 ), with m being the usual (constant) mass [50,66]. Bohm dispersion is known to play an important role in the case of generic dense plasmas [67]. In the present case, since the magnitude of this correction is small, it only becomes important for intermediate wave-number values, as can be seen in figure 2.
Going back to equations (43) and (44), it could be of interest to keep more terms than those that lead to the simplified system of equations (77)-(79). A numerical analysis, at this stage, should be convenient, and extensions to the common hydrodynamical codes could be used (see, e.g. reference [68]), with the Winger function taken as that of equation (54) (i.e. assuming local equilibrium). The main difficulties arise from equations (46) and (48), which are expressed as an infinite sum. Therefore, one should truncate the series at some order in , depending on the desired accuracy.

Conclusions and outlook
Using the Weyl-Wigner description, we have derived a quantum fluid model for a Dirac plasma in single-layer graphene. The massless nature of the quasi-particles was captured by the low-energy limit of the many-body Hamiltonian. Then, the equation of motion for the Wigner matrix elements in phase space was established, with the interacting potential being given by the Hartree approximation. In order to capture both the collisionless and the hydrodynamical regimes, we formally introduced a collision integral. We were able to recover the polarizability function as given by RPA for the collisionless limit, which acts as a benchmark of our framework. Additionally, a closed set of generalized quantum fluid equations was obtained by taking the moments of the kinetic equation. We found an infinite number of -dependent terms arising in the both the fluid current and the pressure, as a consequence of the linear dispersion relation of the carriers. As it is shown, in the case of parabolic systems (i.e. for which the single-particle dispersion relation is ∼p 2 ), these contributions vanish. To the best of our knowledge, such corrections have not been discussed so far in the literature. By neglecting the quantum corrections, our model reduces to existing semiclassical models for Dirac fermions [21,24]. Moreover, a nonlinear relation between the fluid velocity and momentum fields has been put forward in order to construct a hydrodynamic mass. For low temperatures, k B T/μ 1, we found a relation of the form p = γ 2 Mv, where M is Drude's mass and γ = (1 − v 2 /v 2 F ) −1/2 is a local Lorentz-like factor. We stress the fact that a quantitative description for the hydrodynamic mass is a major issue in the current research, therefore making our findings on the hydrodynamic mass an important result of the current kinetic description. We have also computed a corrected plasmon dispersion relation by taking into account the first quantum corrections to the fluid equations, which led to a novel quantum term (∼ 2 q 4 ) that we identified as the well-known Bohm contribution.
We expect our formalism to be particularly suited to describe out-of-equilibrium quantum dynamics of Dirac electrons in a more systematic manner. While the linear response can be indeed described within the framework of the RPA (which we recover upon linearization of the collisionless Wigner transport equation), the nonlinear regimes deserve a more complete description, which is captured by the kinetic equation that we have put forward. This is crucial, for example, in the description of saturation in the late stages of dynamical instabilities. In addition, quantum diffusive processes may be captured by adapting the strategy of quasilinear theories accounting for particle-wave interactions [69]; by working out the collision integral at different levels of approximation, we may obtain quantum versions of the Fokker-Planck equation for relativistic-like particles. Furthermore, quantum mechanical corrections to the collision operator could be a matter of future research. In fact, as shown in reference [60], the very own collision integral may undergo quantum deformations, which go beyond the fermionic or bosonic statistics of the carriers.
Another interesting aspect of our work is the contribution to the hydrodynamical modelling of graphene plasmas. With the ab initio establishment of a hydrodynamic mass, we are able to simulate experimentally relevant configurations, such as hydrodynamical instabilities in field-effect transistor THz emitters [3]. Moreover, at a phenomenological level, our hydrodynamical model is particularly suited to accommodate both shear and odd (or Hall) viscosity terms. In the presence of strong magnetic fields, the Chapman-Enskog approach to hydrodynamics yields a term in the fluid equations of the form ∼ ν o ∇ 2 ij v j , where ν o is some function of the magnetic field and temperature and ij is the Levi-Civita symbol [70]. The appearance of ij term makes the viscosity tensor off-diagonal and, thus, dissipation free, contrary to what happens in the case of even viscosity. Odd viscosity has been identified to play a special role in the topology of acoustic waves [71], and we hope that plasmons may also feature such topological transitions, opening the venue for a range of applications.