Some fluid-dynamic models for quantum electron transport in graphene via entropy minimization

We derive some fluid-dynamic models for electron transport near a Dirac point in graphene. We start from a kinetic model constituted by a set of spinorial Wigner equations, we make suitable scalings (hydrodynamic or diffusive) of the model and we build moment equations, which we close through a minimum entropy principle. In order to do this we make some assumptions: the usual semiclassical approximation ($\hbar\ll 1$), and two further hypothesis, namely Low Scaled Fermi Speed (LSFS) and Strongly Mixed State (SMS), which allow us to explicitly compute the closure.


Graphene electronic properties.
Graphene is a single layer of carbon atoms disposed as an honeycomb lattice, that is, a single sheet of graphite. This new material has attracted the attention of many physicists and engineers thanks to its remarkable electronic properties, which make it an ideal candidate for the construction of new electronic devices with strongly increased performances with respect to the usual silicon semiconductors [1,3,8,12]. The great interest around graphene is attested by the Nobel prize gained in 2010 by Geim and Novoselov for its discovery.
Physically speaking, graphene is a zero-gap semiconductor: in the energy spectrum the valence band intersects the conduction band in some isolated points, named Dirac points; moreover, around such points the electron energy is approximately linear with respect to the modulus of momentum. The system Hamiltonian can be approximated, for low energies and in absence of potential, by the following Dirac-like operator 1 : where v F ≈ 10 6 m/s is the Fermi speed (at zero temperature) and, as usual, denotes the reduced Planck constant; however, in this paper we are not going to use (1) as the system Hamiltonian, because a rigorous discussion of a fluid-dynamic model involving (1) 1991 Mathematics Subject Classification. Primary: 35Q40, 76Y05; Secondary: 82D37. Key words and phrases. Quantum hydrodynamics, electron transport, graphene, quantum entropy principle. 1 We recall the Pauli matrices where, for later convenience, we added the identity matrix σ 0 .
1 would require considering the Fermi-Dirac entropy instead of the Maxwell-Boltzmann, due to the lower unboundedness of the energy spectrum of (1) (we defer such a discussion to a future work). Indeed, in the rest of the paper we will make the hypothesis that the system Hamiltonian is well approximated by the following functional: 2 (2) H = Op with m > 0 parameter (with the dimensions of a mass), whose energy spectrum is bounded from below.
1.2. Quantum fluid-dynamic models, generalities. Recently, mathematical models of fluid-dynamic type has been developed in order to describe quantum transport in semiconductors [2, 4, 5, 6, 9, 10, 11,15,18]. Such models rely on a kinetic formulation of quantum mechanics (QM) by means of Wigner-type equations, and are built by taking suitable moments of these latter; the resulting equations involve the chosen fluid-dynamic moments and usually additional expressions (referred to as not closed terms) which cannot be written as functions of the previous moments without further hypothesis. In order to solve the so-called closure problem, that is, to compute the not closed terms from the known moments, many techniques are employed, e.g. the pure-state hypothesis (which allows to obtain, for a scalar Hamiltonian of typeĤ = − 2 2m ∆ + V (x), the so-called Madelung equations, see [14]), the ad-hoc ansatz (like the Gardner's equilibrium distribution, see [11]), and a strategy of entropy minimization (which will be followed in this paper, in analogy to the method employed in the closure of classical fluid-dynamic systems derived from the Boltzmann transport equation in the classical statistical mechanics, see [4, 5, 6, 13]).
The main advantages of fluid-dynamic models with respect to "basic" tools like Schrödinger, Von Neumann, Wigner equations, are basically two. The first advantage is about physical interpretation: fluid-dynamic models contain already the most physically interesting quantities (like particle, momentum and spin densities), while other models usually involve more "abstract" objects (such as wavefunctions, density operators, Wigner functions), which do not have an immediate physical interpretation; in this latter case, further computations have to be made in order to obtain the quantities of real physical interest from the solution of the model. The second and most important advantage is about numerical computation: fluid-dynamic models for quantum systems with d degrees of freedom are sets of PDEs in d space variables and 1 time variable, while other models usually have more complicated structures (for example, Wigner equations are sets of PDEs in 2d space variables and 1 time variable); so fluid-dynamic models are usually more easily and quickly solvable via numerical computation than other models. 2 The Weyl quantization of a symbol γ is the functional Op (γ) such that, for all ψ test functions, x + y 2 , p ψ(y)e i(x−y)·p/ dydp (for more details see Ref. [7]).
In the following we will present two types of models for quantum electron transport in graphene: • Quantum Diffusion Equations (QDE); • Quantum Hydrodynamic Equations (QHE). Both classes of models are obtained through the following strategy. First, we consider the Wigner equations for the system and we perform on them a diffusive or hydrodynamic scaling; then we take moments of the scaled Wigner system, so that we obtain suitable sets of fluid-dynamic equations containing not closed terms; finally we close these latter by choosing the Wigner function equal to an equilibrium distribution, defined as a minimizer of a suitable entropy functional under the constraints of given moments. We will consider here a kinetic model for quantum transport in graphene, derived from the one-particle Hamiltonian (2), and we will build four (semiclassical) fluid-dynamic models for the system under consideration: two QDE models and two QHE models. We will perform different hypothesis and approximations during the derivation of the fluid-dynamic models in order to overcome the not trivial computational difficulties arising from the process, putting in evidence (when possible) their physical meaning.
The paper is organized as follows. In section 2 we present the kinetic model for quantum transport in graphene that we are going to employ in the paper. In section 3 we build a first diffusive model; in section 4 we build a second diffusive model; in section 5 we build a first hydrodynamic model; in section 6 we build a second hydrodynamic model; in section 7 we present our conclusions and intentions for future work.

A kinetic model for graphene
Let w = w(r, p, t) the system Wigner function, defined for (r, p, t) ∈ R 2 × R 2 × (0, ∞). Notice that, due to the presence of the spin, w is a complex hermitian matrix-valued function instead of a real scalar function; so we can write w = 3 s=0 w s σ s with w s Pauli components of w. 3 Moreover let: The Wigner equations for quantum transport in graphene, associated with the one-particle Hamiltonian H + V , with H defined by (2), are: 3 Given a complex hermitian 2 × 2 matrix a, its Pauli components are real numbers given by: with the pseudo-differential operator (Θ (V )w)(r, p) is defined by: We refer to [17,18] for details about the derivation of (3) from the Von Neumann equation. The terms on the right side of (3) are relaxation terms of BGK type, with g the local thermal equilibrium Wigner distribution, which will be defined later, and τ c is the mean free time (the mean time interval between two subsequent collisions experienced by a particle).

A first diffusive model for graphene.
In the next two sections we will build two diffusive models for quantum electron transport in graphene starting from the kinetic model (3).

3.1.
Wigner equations in diffusive scaling. In order to construct diffusive models for graphene, we make the following diffusive scaling of the Wigner equations: withr,t,p,V satisfying: moreover let us define the semiclassical parameter ǫ, the diffusive parameter τ and the scaled Fermi speed c as: Notice that, if we choose as m the electron mass m e , then c 2 = E F /E cl is the ratio between the Fermi energy E F = m e v 2 F and the classical energy E cl = k B T of the electrons. We perform two main approximations here: the well-known semiclassical hypothesis ǫ ≪ 1, and the following assumption, which we call Low Scaled Fermi Speed (LSFS): By performing the scaling (4)-(6) on the equations (3) under the previous hypothesis, we obtain the following scaled Wigner system: where: 4 (9) and: for the hypothesis (7), ∂ s = ∂ ∂r s for s = 1, 2, 3, and η sjk denotes the only antisymmetric 3 × 3 tensor which is invariant for cyclic permutations of indexes and such that 3.2. Equilibrium Wigner distribution in the diffusive case. We are going to build two diffusive models for graphene by taking moments of Eqs. (8) and closing the resulting fluid-dynamic equations by choosing the Wigner distribution w equal to the equilibrium distribution g, defined as a minimizer of a suitable quantum entropy functional under the constraint of given moments. The moments we choose are the following two: (11) n ± (r) = w 0 (r, p) ± p | p| · w(r, p) dp .
The n ± are the so-called band densities, that is, the partial trace (w.r.t. p) of the quantum operators band projections Π ± : 5 n ± (r) =Tr(Π ± S|r) = tr(P ± (p)w(r, p)) dp ; here S = Op (w) is the density operator which represents the state of the system, and the matrices P ± (p) are the projection operators into the eigenspaces of the classical symbol h 4 From now on we adopt the Einstein summation convention. 5 Given a classical, complex hermitian matrix-valued symbol a(r, p), the expectation value of the quantum observable A = Op (a) when the system is in the state S = Op (w) is: if Tr(SA) is finite. of the quantum Hamiltonian H, that is: and E ± (p) are the eigenvalues of h, that is, the energy bands related to the Hamiltonian H: We define now the quantum entropy functional with two equivalent formulations.
• Wigner function formulation: • Density operators formulation: where Log (w) = Op −1 log Op (w) is the so-called quantum logarithm 6 of w, k B is the Boltzmann constant, T > 0 is the system temperature (which we assume constant) and H = Op (h) is the Hamiltonian (2).
We define now the equilibrium Wigner distribution through the minimum entropy principle. Let W the set of all Wigner functions w defined in R 2 × R 2 with complex hermitian matrix values; for all w ∈ W we set: Moreover, let f = f dp for all f (p) ∈ L 1 (R 2 ), and let n + , n − fluid-dynamic moments. We define the Wigner distribution at local thermal equilibrium related to moments n + , n − as the solution g[n + , n − ] of the problem: We are going to solve the problem (15) through the Lagrange multipliers technique. Let us define the Lagrangian functional in the following equivalent ways: • with the Wigner functions formalism: for w ∈ W and ξ ± (r) real functions; • with the density operators formalism: for S density operator and ξ ± (r) real functions.
The Lagrange multipliers theory tells us (see [4] for details) that a necessary condition for the g[n + , n − ] to solve the problem (15) is that the Gâteaux derivative of L with respect to S must vanish at S = Op (g[n + , n − ]), ξ + = ξ * + , ξ − = ξ * − , for a suitable choice of ξ * + (r), ξ * − (r) (Lagrange multipliers): It can be proved (see [5]) that the solution of (16) has the following form: where )/2 have to be determined in such a way that (recall (14)): 3.3. Formal closure of the fluid equations in the diffusive case. Let n τ + , n τ − the moments of a solution w = w τ of (8) with g given by (17), (18), and let: for w ∈ W smooth enough function. We claim that: Indeed, it is immediate to verify that Eq. (20) is satisfied if g 0 [n τ + , n τ − ] is an even function of p and g[n τ + , n τ − ] is an odd function of p; as a matter of fact, g[n τ + , n τ − ] has this property, because of (17). The proof of this claim is quite similar to the proof of proposition 5.1 at page 293 in [2]: one only has to consider the operator T given by (9), (19) instead of that one used in the paper and consider C as the set of all the p−dependent 2 × 2 matrices with the parity structure: (even, odd, odd, odd) instead of: (even, even, odd, even) .
The following (formal) result holds: Theorem 1. Let us suppose that: n τ ± → n ± as τ → 0 ; then n + , n − satisfy: We refer to [2] for the proof (the operator T is slightly different there, but the proof is still valid in our case).
The system (21) is closed because we already defined (at least formally) the equilibrium distribution g[n τ + , n τ − ]. Nevertheless it is very implicit, as the quantum exponential which appears in Eqs. (21) through Eq. (17) is very difficult to handle with analytically and numerically. As anticipated, in the following we will search for an approximated but more explicit version of Eqs. (21).
3.4. First diffusive model: Explicit construction. In order to write (21) in an explicit way, we will exploit the approximations we have done, that is, the semiclassical and the LSFS hypothesis (given by (7)). We will expand the equilibrium distribution g[n + , n − ] at the first order in ǫ, neglecting O(ǫ 2 ) terms; to do this we start by approximating the quantum exponential of an arbitrary classical symbol with linear ǫ-dependence. Let: for all f (r, p), g(r, p) scalar smooth functions. We apply here the general strategy for computing the semiclassical expansion of the quantum exponential (see [10,11] for details). Let a = a 0 σ 0 + a · σ, b = b 0 σ 0 + b · σ be arbitrary matrix hermitian-valued classical symbols, and let us consider the function: Let us recall the definition of the so-called Moyal product: Op ǫ (f 2 )) between arbitrary classical symbols f 1 , f 2 . It is known [6] that the Moyal product has a semiclassical expansion: and the first three terms of this expansion (the only terms needed in this work) are: (24) Now let us differentiate with respect to β the function g ǫ (β) given by (22). By using the definition (23) of the Moyal product we obtain: and g ǫ (0) = σ 0 . So by expanding the expressions in (25) in powers of ǫ we find: with the initial conditions: The equations (26), (27) with the initial conditions (28) can be exactly solved in this order to obtain the O(ǫ 2 )−approximation of Exp ǫ (a) = g ǫ (1): in fact, each equation is a linear ODE with constant coefficients. It is easy to find the leading term in the expansion of g ǫ (β): (29) g (0) (β) = exp(βa) = e βa 0 cosh(β| a|)σ 0 + sinh(β| a|) a | a| · σ .
We now have to explicitly compute the first order correction of g ǫ (β) from (27); to this aim, it is useful to employ some properties of the Pauli matrices. It is easy to verify that, for a, b arbitrary hermitian matrix-valued classical symbols holds: where we defined: The relations (30) allow us to reduce the calculus of the matrix g (1) (β) to that of its Pauli components; if fact, due to (30), (27) becomes: In order to solve (31), let us consider the homogeneous problem: The problem (32) can be solved with elementary techniques, finding that the vector X(β) = [x 0 (β), x(β)] is given by: X(β) = S a (β)X(0) β > 0 , with the semigroup operator S a (β) defined by: S a (β) = e βa 0 cosh(β| a|) sinh(β| a|) α T sinh(β| a|) α (cosh(β| a|) − 1) α ⊗ α + I 3×3 with α ≡ a/| a|, and α T denotes the transpose of the vector α. Now the semigroup theory allows us to write the solution of (31): finally, from (29), (34), (35), (36) we obtain: So we have explicitly computed the first-order semiclassical expansion of g ǫ (β) = Exp ǫ (β(a+ ǫb)). We point out that in the scalar case the odd order terms in the semiclassical expansion of the quantum exponential are zero, while this does not happen in the spinorial case, due to the noncommutativity of the matrix product, which increases much the complexity in computation with respect to the scalar case. The next step is the computation of g[n + , n − ] given by (17), (18). To achieve this goal, we substitute in (29), (37), (38) β = 1 and a + ǫb = −h ξ ; this means, due to (17), (10): By making all the straightforward computations needed, we finally find: (39) with: and: (41) n 0 = 1 2 (n + + n − ) particle density, n σ = 1 2 (n + − n − ) pseudo-spin polarization.

A second diffusive model for graphene
We are going to build another diffusive model for quantum transport in graphene, starting again from the Wigner equations in diffusive scaling (8), considering the same fluiddynamic moments n ± of the Wigner distribution w(r, p, t) and taking again as the equilibrium distribution the one given in (17), (18); however, we will make stronger assumptions than (7), which will allow us to consider also O(ǫ 2 )−terms in the fluid equations. 4.1. Assumptions. We make the semiclassical hypothesis ǫ ≪ 1 like in the previous model, and another hypothesis, stronger than (7), which we call Strongly Mixed State (SMS): (Recall the definitions (6), (17) of c and B, respectively). These further assumptions are necessary to overcome the computational difficulties arising from the spinorial nature of the problem: without these hypothesis, it would be hard to compute the second order expansion of the equilibrium distribution.
We will see that the two approximations will result in the fact: This means, from a physical point of view, that the charge carriers have approximately the same probability of being found in the conduction band or in the valence band of the energy spectrum, or equivalently, there is little difference (with respect to the total charge density) between the electron density and the hole density. Analytically, the main consequence of the SMS hypothesis (44) will be the decoupling of the modified Hamiltonian in a scalar part of order 1 and a spinorial perturbation of order ǫ; this fact will be very useful in computations.
4.2. Second diffusive model: Explicit construction. For the sake of simplicity let us redefine B → ǫB and consider B = O(1). Under our hypothesis, the classical symbol of the modified Hamiltonian becomes: that is, as anticipated, the modified Hamiltonian h ξ decouples in a scalar part of order O(1) and a spinorial part of order O(ǫ), so that h ξ can be seen as a small perturbation of a scalar Hamiltonian. We are going to see that this fact leads to remarkable simplifications in computations.
In order to explicitly compute the equilibrium distribution, let us recall (29), (31). In this case we have: so (29) becomes the scalar function: (47) g (0) (β) = e βa σ 0 , and (31) takes the form: then the solution is: Now we can compute also the second-order correction to g ǫ (β), thanks to the approximations we have done. From (25) we obtain: so one finds that: 8 So from (46), (47), (48), (49) we find: whereg(β) is the second-order approximation of the quantum exponential of the scalar symbol βa: It is convenient to work from now on with the moments n 0 , n σ given by (41) instead of n ± ; for this reason we will write g[n 0 , n σ ] instead of g[n + , n − ], being no possibility of misunderstanding in this. The constraint (18) can so be rewritten in terms of g[n 0 , n σ ] in the following way 9 : (52) g 0 [n 0 , n σ ] = n 0 , g σ [n 0 , n σ ] = n σ ; By imposing the constraint on the moments and making all the straightforward computations needed, we obtain the following result: whereg[n] is the O(ǫ 3 )−approximation of the scalar quantum Maxwellian with g[n] = n, that is:g for an arbitrary positive scalar function n(r) and I is the 3 × 3 identity matrix. We point out that, during the computations, we find (see the Appendix): (54) so the constraint (45) is satisfied. Finally, by using Eqs. (53) to compute the terms in (21) and making all the long but straightforward computations needed, we obtain the second semiclassical quantum diffusive model we were looking for: where Γ,γ are given by (43), V B is the so-called Bohm potential : and we defined the following constants: We will dedicate the remaining part of this paper to the construction of hydrodynamic models for graphene.

A first hydrodynamic model for graphene
In the following of this paper we will build two hydrodynamic models for quantum electron transport in graphene following a strategy similar to that one employed in the construction of the diffusive models (42), (55). 5.1. Wigner equations in hydrodynamic scaling. Let us reconsider the Wigner equations (3) and perform the following hydrodynamic scaling: withr,t,p,V satisfying: moreover, let us define the semiclassical parameter ǫ and the hydrodynamic parameter τ as: Let c (scaled Fermi speed) be given again by (6). We make the same assumptions as in the first diffusive model, that is, the semiclassical hypothesis ǫ ≪ 1 and the LSFS hypothesis (7). Let γ be defined as in (10). If we perform the scaling (57) -(59) on (3) under the assumptions we have made, we obtain the following scaled Wigner system: g is the quantum thermal equilibrium Wigner distribution, which will be different from the diffusive case because we are going to choose a different set of moments.
5.2. Equilibrium Wigner distribution in the hydrodynamic case. The moments we choose are the following six: n s = w s dp s = 0, 1, 2, 3 , J k = p k w 0 dp k = 1, 2 . n 0 is the particle density, n = (n 1 , n 2 , n 3 ) is the spin vector, J = (J 1 , J 2 , 0) is the flow vector.
Note that the flow vector has only two components because graphene is a two-dimensional object.
We choose (12), (13) as quantum entropy like in the diffusive case, and we define the equilibrium distribution g[n 0 , n, J] again through the MEP. Let n 0 , n, J denote given moments. We define the Wigner distribution at local thermal equilibrium related to moments n 0 , n, J the solution g[n 0 , n, J] of the problem: A(g[n 0 , n, J]) = min A(w) : w ∈ W , w 0 = n 0 , w = n , pw 0 = J .
The proof of theorem 2 is analogue to the proof of Theorem 4.2 in [11], pages 38 − 40: one must consider (60) as Wigner equations and use the weight functions k(p) = {1, p} for the first equation in (60) and k(p) = 1 for the second equation.
Eqs. (65) are a closed system of hydrodynamic equations, indeed g[n 0 , n, J] is a function of the moments n 0 , n, J only; however, the system (65) is a very implicit model; in the next subsection we will build an approximated but more explicit version of (65) by exploiting the hypothesis we have done (semiclassical and (7)). 5.4. First hydrodynamic model: Explicit construction. It is possible to write the first-order approximation of g[n 0 , n, J] given by (63) by following a strategy similar to that one employed to compute the approximation of the equilibrium distribution for the first diffusive model. More precisely, we consider (29), (37), (38) with: then we impose the constraints (64). Skipping the long but straightforward computations needed, we obtain: By using the O(ǫ 3 )-expansion of g[n 0 , n, J] to explicitly calculate the terms inside Eqs. (65) up to O(ǫ 3 ) we find the following hydrodynamic model: where V B is the Bohm potential, defined by (56).

Conclusions and perspectives for the future
In this paper we have proposed four fluid-dynamic models for quantum electron transport in graphene: the first two ones are diffusive-type models, the last two ones are hydrodynamic-type models. All models have been built by using a statistical closure of the moment equations derived from the Wigner system (3) based on the minimization of the quantum entropy functional (12), (13). For each class of models (diffusive or hydrodynamic) we have have presented two strategies of construction: a first approach, consisting in making weaker hypothesis (semiclassical and Low Scaled Fermi Speed) but neglecting O(ǫ 2 ) terms in the fluid equations, which leads to the models (42), (68); and a second approach, consisting in making stronger hypothesis (semiclassical and Strongly Mixed State) but neglecting O(ǫ 3 ) terms, which leads to the models (55), (73).
In the next future, we are going to perform numerical simulations with the models (55) and (73), attempting to reproduce the features of charge transport in graphene, with particular attention to the so-called "Klein paradox" (unimpeded penetration of electrons through arbitrary high potential barriers, see [3,12]). Moreover we plan to build new fluid-dynamic models for quantum transport in graphene, again by an entropy minimization strategy, but using Fermi-Dirac entropy instead of Maxwell-Boltzmann entropy, the former being more suitable for describing electrons than the latter.