1 Introduction

In its pristine form graphene is a monolayer sheet of carbon atoms assembled in a honeycomb structure and its electronic properties are due to the behavior of the electric charges in proximity of the maximum of the valence band and of the minimum of the conduction band, that correspond to the six vertices of the hexagonal Brillouin zone. In these points the two bands are in contact making graphene a gapless conductor. The charge properties are studied in only two equivalent points of the Brillouin zone, indicated by K and \(K'\) (Dirac points) and the dispersion relation is

$$\begin{aligned} \varepsilon (\mathbf{k}) = \pm \hbar v_F |\mathbf{k}|, \end{aligned}$$
(1)

(\(+\) is for the conduction band and − is for the valence band), with \(\hbar \) the Dirac constant, \(v_F\) (constant) the Fermi velocity and \(|\mathbf{k}|\) the modulus of the wave number.

This peculiar band structure makes graphene a semimetal with high mobility and good conductivity and, as a consequence, one of the best candidates for the realization of new nanodevices [1, 2]. On the other hand the absence of the band gap does not permit to switch off the current and even at equilibrium conditions the conductivity has a minimum value of the order of \(\sim 4 e^2 / \hbar \) [2]. This makes graphene in its pristine form not suitable for the realization of Fields Effect Transistors (FETs).

A possible solution to overcome this problem is the use of graphene nanoribbons, narrow strips of graphene that are semiconductors because the quantum confinement along with one direction creates a bandgap tuned as a function of the width and of the edge structure [3]. In idealized situations the edges have zigzag or armchair structure, however today it is very difficult the production of graphene ribbon with one well defined edge, but they are available with a high degree of edge disorder. This leads to the investigation of the behavior of the charge transport in graphene taking into account the bandgap [4] and the scattering with edge [5, 6] as well as the well known scattering with the lattice [8].

Here we develop a hydrodynamical model for simulating charge transport in graphene nanoribbons by using the maximum entropy principle (hereafter MEP), a theoretical approach already examined theoretically for graphene and other semiconductors in [8,9,10].

The evolution equations for the macroscopic variables like density, energy, velocity and energy-flux are obtained by taking the moments of the transport equations. The constitutive relations, needed to have a closed system of balance laws, are deduced by resorting to MEP in a similar way to the Levermore moment approach [11]. It has been proved in [8, 9] that maximization problem, MEP leads to, is globally solvable in the physically relevant region of the field variables and that in the same region the evolution equations form a hyperbolic system of conservation laws.

Other macroscopic models are also available in the literature. For example in [12, 13] a spinorial Wigner function has been employed to get drift-diffusion like models for charge transport in graphene, including spin transport as well. Quantum effects have been also introduced with a Wigner formalism in [14,15,16,17].

The plan of the article is the following one. In Sect. 2 the semiclassical description is recalled and in Sect. 3 the macroscopic model is developed. In Sect. 4 the closure relations are furnished and discussed while in the subsequent section the production terms are calculated. In penultimate section the results of a homogeneous case are presented and discussed.

2 Kinetic Model

The confinement of the charges in graphene nanoribbons has been widely treated and we refer the reader to [3, 4, 18,19,20] for details. Here we suppose a high degree of edge disorder and assume the model proposed by [4] that, in good agreement with DFT calculation, furnishes the following energy dispersion relation

$$\begin{aligned} \varepsilon _{\pm }(\mathbf{k}) = \pm \hbar v_F \sqrt{k^2_x+k^2_y+\left( \frac{\pi }{L} \right) ^2} \end{aligned}$$
(2)

that unlike the (1) depends on the width L of the ribbon, so there is an energy gap given by

$$\begin{aligned} \varDelta \varepsilon = \frac{2 \hbar v_F \pi }{L}. \end{aligned}$$

By considering the interaction with the optical phonons (see forward for the details) an electron can absorb a phonon and jump from the valence band to the conduction band (inter-band scattering) or jump in the opposite direction by emitting a phonon. In this paper we assume a value of L so that \(\varDelta \varepsilon > \hbar \omega \) (where \(\hbar \omega \) is the highest phonon energy) and a Fermi energy greater than the bottom of the conduction band. In this way the inter-band scattering is avoided and the charge current is given only by electrons in the conduction band.

Under these conditions the charge transport in the framework of the semiclassical models is described by the Boltzmann equation

$$\begin{aligned} \frac{\partial f(\mathbf{r}, \mathbf{k}, t)}{\partial t} + \mathbf{v} \cdot \nabla _\mathbf{r} f (\mathbf{r}, \mathbf{k}, t) - \frac{e}{\hbar } \mathbf{E} \cdot \nabla _\mathbf{k} f (\mathbf{r}, \mathbf{k}, t) =\mathcal {C} (\mathbf{r},\mathbf{k}, t) \end{aligned}$$

where \(f (\mathbf{r}, \mathbf{k}, t)\) is the distribution function for the electrons at position \(\mathbf{r}\), with wavevector \(\mathbf{k}\) and time t. e is the absolute value of the elementary (electron) charge, \(\nabla _\mathbf{r}\) and \(\nabla _\mathbf{k}\) are the gradients with respect to the position and wavevector respectively, \(\mathbf{E}\) is the electric field. \(\mathbf{v}\) is the microscopic velocity whose expression is

$$\begin{aligned} \mathbf{v}=\frac{1}{\hbar } \nabla _\mathbf{k} \varepsilon = \frac{v_F}{\sqrt{k^2_x+k^2_y+ \left( \frac{\pi }{L} \right) ^2 } } \mathbf{k}. \end{aligned}$$

The collisional term has the following generic form (to avoid cumbersome notation we omit the dependence on \(\mathbf{r}\) and t)

$$\begin{aligned}&\mathcal {C}(\mathbf{k})=\frac{1}{(2 \pi )^2} \int _{\mathbb {R}^2} \left[ w(\mathbf{k},\mathbf{k'}) f(\mathbf{k'}) \left( 1-f(\mathbf{k}) \right) - w(\mathbf{k'},\mathbf{k}) f(\mathbf{k}) \left( 1-f(\mathbf{k'}) \right) \right] d^2 \mathbf{k'} . \end{aligned}$$

and describes a single scattering of the electron from a state with wavevector \(\mathbf{k}\) to a state with wavevector \(\mathbf{k}'\) with transition rate \(w(\mathbf{k},\mathbf{k'})\).

\(\mathcal {C}(\mathbf{k})\) is the sum of the contributions of the scattering of the electrons with acoustic and optical phonons and with the edge of the ribbon.

For the evaluation of production terms in this paper we consider a graphene ribbon of width L along axis x so that the edges are located at \(y=0\) and \(y=L\) and the transport direction is along x direction (see Fig. (1)).

Fig. 1
figure 1

Schematic representation of a graphene nanoribbon, the transport direction is along x axis

Full size image

In the case of acoustic phonons, usually one considers the elastic approximation. Therefore the collision is only intraband and can be written as

$$\begin{aligned} {\mathcal {C}} (\mathbf{k}) = \frac{2}{(2 \pi )^2} \int _{\mathbb {R}^2} w^{(ac)}(\mathbf{k'}, \mathbf{k}) \left( f(\mathbf{k'}) - f (\mathbf{k}) \right) \delta (\varepsilon (\mathbf{k'}) - \varepsilon (\mathbf{k})) \, d^2 \mathbf{k'}, \end{aligned}$$
(3)

where \(w^{(ac)}(\mathbf{k'}, \mathbf{k}) = A^{(ac)} (1 + \cos \theta '')\) with \(A^{(ac)} = \displaystyle {\frac{\pi D_{ac}^2 k_B T}{\rho \hbar v_{ac}}}\). Here \(\theta ''\) is the convex angle between \(\mathbf{k}\) and \(\mathbf{k'}\), \(D_{ac}^2\) is the acoustic phonon coupling constant, \(v_{ac}\) is the sound speed in graphene, \(\rho \) the graphene areal density and T the graphene lattice temperature, which in this paper will be kept constant (the reader interested in thermal effects in graphene can refer to [8, 21]). There are three relevant optical phonon scatterings: the longitudinal optical (LO), the transversal optical (TO) and the K phonons. We label them by the index s. The general expression of the collision term for a generic optical phonon interaction is given by

$$\begin{aligned}&\!\!\!\mathcal {C}_{\alpha \, \alpha '}^{(s)} (\mathbf{k})=\int _{\mathbb {R}^2} \left[ w_{\alpha ' \, \alpha }^{(s)}(\mathbf{k'}, \mathbf{k}) f^{\alpha '} (\mathbf{k'}) \left( 1 - f^{\alpha } (\mathbf{k}) \right) - w_{\alpha \, \alpha '}^{(s)}(\mathbf{k}, \mathbf{k'}) f^{\alpha } (\mathbf{k}) \left( 1 - f^{\alpha '} (\mathbf{k'}) \right) \right] \, d^2 \mathbf{k'} \nonumber \!\!\!\!\!\!\!\\ \end{aligned}$$
(4)

where \(w_{\alpha \, \alpha '}^{(s)}(\mathbf{k}, \mathbf{k'}) = w_{\alpha \, \alpha '}^{(s, +)}(\mathbf{k}, \mathbf{k'}) + w_{\alpha \alpha '}^{(s, -)}(\mathbf{k}, \mathbf{k'})\), \(s=LO, TO,K\). If \(\alpha = \alpha '\) one has an intraband scattering otherwise an inter-band scattering, here we consider only intraband scattering and omit the \(\alpha \) index. In the case \(s = LO, TO\) the expression of \(w_{\alpha \, \alpha '}^{(s, \pm )}(\mathbf{k}, \mathbf{k'})=w^{(s, \pm )}(\mathbf{k}, \mathbf{k'})\) is

$$\begin{aligned} w^{(s, \pm )}(\mathbf{k}, \mathbf{k'}) = A^{(op)} D_{\varGamma }^2 \left[ 1 - \eta _s \alpha \alpha ' \cos (\theta + \theta ') \right] \left( N_B^s + \frac{1}{2} \pm \frac{1}{2} \right) \delta (\varepsilon '_{\alpha '} - \varepsilon _{\alpha } \pm {\hbar {\omega }}) \end{aligned}$$

where \(A^{(op)} = \pi / \rho \omega _s\), \(D_{\varGamma }^2\) is the optical phonon coupling constant and \(N_B^s\) is the Bose-Einstein distribution

$$\begin{aligned} N_B^s = \displaystyle {\frac{1}{e^{{\hbar {\omega }}_s/k_B T} - 1}} \end{aligned}$$

with \({\hbar {\omega }}_s\) phonon energy. \(\theta \) and \(\theta '\) are the convex angles between \(\mathbf{k}\) and \(\mathbf{k'} - \mathbf{k}\) and between \(\mathbf{k'}\) and \(\mathbf{k'} - \mathbf{k}\). \(\eta _s\) is a parameter which takes the following values

$$\begin{aligned} \eta _s = \left\{ \begin{array}{rl} 1 &{} \text{ if } \quad s = \text{ LO }\\ -1 &{} \text{ if } \quad s = \text{ TO }\\ \end{array}\right. \end{aligned}$$

Here we assume for LO and TO phonon energy \(\hbar \omega _{LO}=\hbar \omega _{TO}\) and hereafter it will be indicated by \(\hbar \omega _{opt} \). Moreover, under the previous assumptions, \(w^{(s)}(\mathbf{k},\mathbf{k}') =w^{(s,+)}(\mathbf{k},\mathbf{k}')+w^{(s,-)}(\mathbf{k},\mathbf{k}')\) reads

$$\begin{aligned} w^{(s)}(\mathbf{k}, \mathbf{k'}) = A^{(op)} D_{\varGamma }^2 \left( N_B^s + \frac{1}{2} \pm \frac{1}{2} \right) \delta (\varepsilon '_{\alpha '} - \varepsilon _{\alpha } \pm {\hbar {\omega }}) \end{aligned}$$

If \(s = K\) the transition rate is given by

$$\begin{aligned} w^{(s, \pm )}(\mathbf{k}, \mathbf{k'}) = A^{(K)} D_{K}^2 \left( 1 - \alpha \alpha ' \cos \theta '' \right) \left( N_B^s + \frac{1}{2} \pm \frac{1}{2} \right) \delta (\varepsilon _{\alpha '} - \varepsilon _{\alpha } \pm {\hbar {\omega }}_s) \end{aligned}$$

where \(A^{(K)} = \pi / \rho \omega _s\) and \(D_{K}^2\) is the K-phonon coupling constant. The physical parameters for electron-phonon scatterings are reported in Table 1.

Table 1 Physical parameters for the electron-phonon collision term
Full size table

To include the edge scattering we suppose disordered edges and consider an additional term as proposed by [6] whose expression is

$$\begin{aligned} \mathcal {C}^{(edge)} (\mathbf{k}) =\frac{N_i V_0^2}{2 \pi \hbar W} \int _{\mathbb {R}^2} e^{-2(k_x-k'_x)^2a^2} \left( f(\mathbf{k'})- f(\mathbf{k}) \right) \delta ( \varepsilon (\mathbf{k}) - \varepsilon (\mathbf{k'}) ) d^2 \mathbf{k'} \end{aligned}$$

where \(N_i\) is the linear density of the scatterer (defects) along with the graphene edge, and the quantity

$$\begin{aligned} V_0 e^{-a^2(k'_x-k_x)^2} \end{aligned}$$

is the matrix element due to the potential of a single scatterer at the edge, being \(V_0\) a constant and a a characteristic range, along with the edge direction, so that electron scattering with rather strong \(k_x\)-momentum transfer, \(|k_x-k'_x|>1/a\), is effectively suppressed. For the edge scattering we adopt the parameter reported in Table 2.

Direct simulations of the kinetic models [22, 23] require a huge computational effort (similar difficulties arise also for other semiconductors materials). This has prompted the formulation of hydrodynamical models given by balance equations for the macroscopic quantities of interest like density, energy, current.

Table 2 Physical parameters for the electron-edge collision term
Full size table

3 Moment Equations

Macroscopic quantities can be defined as moments of the distribution function with respect to some suitable weight functions, assuming a sufficient regularity for the existence of the involved integrals. In particular we propose as a set of moment equations for the conduction band the balance equations for average densities \(\rho \), velocities V, energies W and energy fluxes S defined as

$$\begin{aligned} \rho= & {} \frac{2}{(2 \, \pi )^2} \int _{\mathbb {R}^2} f (\mathbf{r},\mathbf{k},t) \,d^2 \mathbf{k},\\ \rho \mathbf{V}= & {} \frac{2}{(2 \, \pi )^2} \int _{\mathbb {R}^2} f (\mathbf{r},\mathbf{k},t) \,\mathbf{v} \, d^2 \mathbf{k},\\ \rho W= & {} \frac{2}{(2 \, \pi )^2} \int _{\mathbb {R}^2} f (\mathbf{r},\mathbf{k},t) \,\varepsilon \, d^2 \mathbf{k},\\ \rho \mathbf{S}= & {} \frac{2}{(2 \, \pi )^2} \int _{\mathbb {R}^2} f (\mathbf{r},\mathbf{k},t) \,\varepsilon \mathbf{v} \, d^2 \mathbf{k}. \end{aligned}$$

By integrating the Boltzmann equation with respect to \(\mathbf{k}\), one has the balance equations for the macroscopic quantities defined above

$$\begin{aligned}&\frac{\partial }{\partial t} \rho + \nabla _\mathbf{r} \cdot \left( \rho \, \mathbf{V} \right) = \rho \, C_{\rho } \nonumber \\&\quad \frac{\partial }{\partial t} \left( \rho \, W \right) + \nabla _\mathbf{r} \cdot \left( \rho \, \mathbf{S} \right) - \, e \rho \mathbf{E} \cdot \mathbf{V} = \rho \, C_{W} \nonumber \\&\quad \frac{\partial }{\partial t} \left( \rho \, \mathbf{V} \right) + \nabla _\mathbf{r} \cdot \left( \rho \, \mathbf{F}^{(0)} \right) -\, e \, \rho \mathbf{G}^{(0)} \cdot \mathbf{E} = \rho \, C_{\mathbf{V}} \nonumber \\&\quad \frac{\partial }{\partial t} \left( \rho \, \mathbf{S} \right) + \nabla _\mathbf{r} \cdot \left( \rho \, \mathbf{F}^{(1)} \right) \, e \, \rho \mathbf{G}^{(1)} \cdot \mathbf{E} = \rho \, C_{\mathbf{S}}. \end{aligned}$$
(5)

Of course other choices are possible, but the accuracy in the choice of these moments has been well tested in [14]. In particular it is crucial to include the average energy and energy-flux because the thermal effects are supposed to be very relevant in graphene, as confirmed at a kinetic level. We note that, besides the average densities, velocities, energies and energy fluxes, additional quantities appear (omitting for the sake of brevity the dependence on space and time)

$$\begin{aligned} \rho \, C_{\rho }= & {} \frac{2}{(2 \, \pi )^2} \int _{\mathbb {R}^2} {\mathcal {C}} (\mathbf{k}) \, d^2 \, \mathbf{k}, \, \rho \, C_{\mathbf{V}} = \frac{2}{(2 \, \pi )^2} \int _{\mathbb {R}^2} \mathbf{v} \, {\mathcal {C}} (\mathbf{k}) \, d^2 \, \mathbf{k}, \\ \rho \, C_{W}= & {} \frac{2}{(2 \, \pi )^2} \int _{\mathbb {R}^2} \varepsilon \, {\mathcal {C}} (\mathbf{k}) \, d^2 \, \mathbf{k}, \, \rho \, C_{\mathbf{S}} = \frac{2}{(2 \, \pi )^2} \int _{\mathbb {R}^2} \varepsilon \, \mathbf{v} {\mathcal {C}} (\mathbf{k}) \, d^2 \, \mathbf{k}\\ \rho \, \left( \begin{array}{c} \mathbf{F}^{(0)} \\ \mathbf{F}^{(1)} \end{array} \right)= & {} \frac{2}{(2 \, \pi )^2} \int _{\mathbb {R}^2} \left( \begin{array}{c} 1\\ \varepsilon \end{array} \right) \mathbf{v} \otimes \mathbf{v} f (\mathbf{r},\mathbf{k},t) \,d^2 \mathbf{k}, \\ \rho \, \left( \begin{array}{c} \mathbf{G}^{(0)} \\ \mathbf{G}^{(1)} \end{array} \right)= & {} \frac{2}{\hbar \, (2 \, \pi )^2} \int _{\mathbb {R}^2} f (\mathbf{r},\mathbf{k},t) \nabla _\mathbf{k} \left( \begin{array}{c} \mathbf{v} \\ \varepsilon \, \mathbf{v} \end{array} \right) \,d^2 \mathbf{k}, \qquad \end{aligned}$$

that must be expressed as function of the basic variables \(\rho \), W, \(\mathbf{V}\), \(\mathbf{S}\) (closure problem).

4 Closure Relations

A well theoretically motivated way to get the desired closure relations is to resort to the maximum entropy principle according to which the electron distribution functions is estimated with the distributions \(f_{MEP}\) obtained solving the following constrained optimization problem [8, 10]. For fixed \(\mathbf{r}\) and t,

$$\begin{aligned} f_{MEP}= & {} {\mathop {{{\,\mathrm{argmax}\,}}}\limits _{f \in \mathscr {F}_{\varPsi } }} S[f]\quad \text{ under } \text{ the } \text{ constraints: } \nonumber \\ \rho= & {} \frac{2}{(2 \, \pi )^2} \int _{\mathbb {R}^2} f(\mathbf{r},\mathbf{k},t) \,d^2 \mathbf{k},\quad \rho W = \frac{2}{(2 \, \pi )^2} \int _{\mathbb {R}^2} f (\mathbf{r},\mathbf{k},t) \,\varepsilon \, d^2 \mathbf{k}, \end{aligned}$$
(6)
$$\begin{aligned} \rho \mathbf{V}= & {} \frac{2}{(2 \, \pi )^2} \int _{\mathbb {R}^2} f (\mathbf{r},\mathbf{k},t) \,\mathbf{v} \, d^2 \mathbf{k}, \quad \rho \mathbf{S} = \frac{2}{(2 \, \pi )^2} \int _{\mathbb {R}^2} f (\mathbf{r},\mathbf{k},t) \,\varepsilon \mathbf{v} \, d^2 \mathbf{k}, \end{aligned}$$
(7)

where S[f] is the entropy of the system, which in the semiclassical approximation reads

$$\begin{aligned} S[f] = - \frac{2 k_B}{(2 \pi )^2}\, \int _{\mathbb {R}^2} \left[ f \ln f + \left( 1 - f \right) \ln \left( 1 - f \right) \right] \, d^2 \, \mathbf{k} \end{aligned}$$

and \(\mathscr {F}_{\varPsi }\) is the space of the distribution functions such that the expectation values with respect to the considered weight functions there exist [8]. In order to take into account the constraints let us introduce the Lagrange multipliers \(\lambda \), \(\lambda _{W}\), \(\lambda _{\mathbf{V}}\), \(\lambda _{\mathbf{S}}\) and the Legendre transform of S

$$\begin{aligned} S'= & {} S + \left[ \lambda \left( \rho - \frac{2}{(2 \, \pi )^2} \int _{\mathbb {R}^2} f (\mathbf{r},\mathbf{k},t) \,d^2 \mathbf{k} \right) \right. \\&\left. + \lambda _{\mathbf{V}} \left( \rho \mathbf{V} - \frac{2}{(2 \, \pi )^2} \int _{\mathbb {R}^2} f (\mathbf{r},\mathbf{k},t) \,\mathbf{v} \, d^2 \mathbf{k} \right) \right. \\&\left. +\lambda _{W} \left( \rho W - \frac{2}{(2 \, \pi )^2} \int _{\mathbb {R}^2} f (\mathbf{r},\mathbf{k},t) \,\varepsilon \, d^2 \mathbf{k}, \right) \right. \\&\left. + \lambda _{\mathbf{S}} \left( \rho \mathbf{S} - \frac{2}{(2 \, \pi )^2} \int _{\mathbb {R}^2} f (\mathbf{r},\mathbf{k},t) \,\varepsilon \mathbf{v} \, d^2 \mathbf{k} \right) \right] . \end{aligned}$$

The optimality condition \(\delta S' = 0\) gives

$$\begin{aligned} f_{MEP} (\mathbf{r},\mathbf{k},t) = \frac{1}{1+ \exp \left[ \lambda +\lambda _{W} \varepsilon + \left( \lambda _{\mathbf{V}}+ \varepsilon \lambda _{\mathbf{S}} \right) \cdot \mathbf{v} \right] }. \end{aligned}$$

The multiplicative constant \(k_B/ 2 \pi ^2\) has been included into the multipliers for the sake of simplifying the notation. In order to invert the relations (6)–(7) and express the Lagrangian multipliers as functions of the basic variables, we linearize the above distributions around the isotropic part by neglecting the anisotropic part of higher order

$$\begin{aligned} f_{MEP} \approx f^{(i)}({\varepsilon })+f^{(a)}({\varepsilon }) \end{aligned}$$
(8)

where

$$\begin{aligned} f^{(i)}({\varepsilon }) = \frac{1}{1+ e^{\lambda +\lambda _{W} \varepsilon }}, \quad f^{(a)}({\varepsilon }) = - \frac{ e^{\lambda +\lambda _{W} \varepsilon }}{\left( 1+ e^{\lambda + \lambda _{W} \varepsilon }\right) ^2 } \, (\lambda _{\mathbf{V}} + \varepsilon \lambda _{\mathbf{S}} ) \cdot \mathbf{v} \end{aligned}$$
(9)

The comparison between linear and non-linear models can be found in [24]. A comparison between both models and numerical solutions of the kinetic equations is presented in [25].

We remark that the expectation values of any tensorial function of the type \(\sum _{i=1}^m \otimes v_i\varepsilon ^n\) with respect to the density (8) there exists for each m and n non negative integers provided \(\lambda _{W} > 0\) which constitutes the realizability condition.

Since at equilibrium \(\lambda _{W} = 1/k_B T >0\), the equilibria are interior points of the realizability region. It is simple to prove the following property.

Theorem

The distribution (8) satisfies the constraints (6)–(7) if the Lagrangian multipliers are given by

$$\begin{aligned} \rho= & {} \frac{1}{\pi (\hbar v_F)^2} I_1 (\lambda , \lambda _{W}), \quad \rho W= \frac{1}{\pi (\hbar v_F)^2} I_2 (\lambda , \lambda _{W}) \end{aligned}$$
(10)
$$\begin{aligned} \left( \begin{array}{c} \lambda _{ \mathbf{V}} \\ \lambda _{\mathbf{S}} \end{array} \right)= & {} \left( \begin{array}{cc} a_{11} &{} a_{12}\\ a_{21} &{} a_{22} \end{array}\right) \left( \begin{array}{c} \rho \mathbf{V} \\ \rho \mathbf{S} \end{array} \right) := {\mathcal {A}} \left( \begin{array}{c} \rho \mathbf{V} \\ \rho \mathbf{S} \end{array} \right) \end{aligned}$$
(11)

where

$$\begin{aligned} {\mathcal {A}}= & {} \frac{2 \pi \hbar ^2}{D e^{\lambda }} \left( \begin{array}{cc} -J_3 (\lambda , \lambda _{W}) &{} J_2 (\lambda , \lambda _{W})\\ J_2 (\lambda , \lambda _{W}) &{} - J_1 (\lambda , \lambda _{W}) \end{array} \right) \end{aligned}$$
(12)

with

$$\begin{aligned} I_n (\lambda , \lambda _W)= & {} \displaystyle {\int _{\frac{\pi \hbar v_F}{L}}^{+\infty } \frac{\varepsilon ^n}{1 + e^{\lambda +\lambda _{W} \varepsilon } } d\, \varepsilon }, \end{aligned}$$
(13)
$$\begin{aligned} J_n (\lambda , \lambda _{W})= & {} \displaystyle {\int _{\frac{\pi \hbar v_F}{L}}^{+\infty }} \frac{\varepsilon ^n e^{\lambda _{W} \varepsilon }}{\left( 1 + e^{\lambda +\lambda _{W} \varepsilon }\right) ^2 } \left[ 1- \left( \frac{\pi \hbar v_F}{L \varepsilon } \right) ^2 \right] d\, \varepsilon \end{aligned}$$
(14)
$$\begin{aligned} D= & {} J_1(\lambda , \lambda _{W}) J_3(\lambda , \lambda _{W})-J_2^2(\lambda , \lambda _{W}). \end{aligned}$$
(15)

For the proof the reader is referred to [8, 9].

It is also possible to prove that the relationships (10) and (11) are locally and globally invertible in the region \(\lambda \in \mathbb {R}\), \(\lambda _{W} >0\).

Also for these propositions the reader can find the proof in [8, 9].

Moreover we observe that inverse relation of (11) is

$$\begin{aligned}&\begin{pmatrix} \rho \mathbf{V}\\ \rho \mathbf{S} \end{pmatrix} =-\frac{e^\lambda }{2 \pi \hbar ^2} \begin{pmatrix} J_1 &{} J_2 \\ J_2 &{} J_3 \end{pmatrix} \begin{pmatrix} \lambda _\mathbf{V} \\ \lambda _\mathbf{S} \end{pmatrix} :=\mathcal {J} \begin{pmatrix} \lambda _\mathbf{V} \\ \lambda _\mathbf{S} \end{pmatrix} \end{aligned}$$

and is simple to prove that \(\det {\mathcal {J}} \ne 0\) and \({\mathcal {J}}^{-1}=\mathcal {A}\). \(\square \)

Once \(f_{MEP}\) has been expressed in terms of the basic variables, from the relative definition we are able to determine the closure relations for fluxes and production terms. In particular for fluxes one obtains

$$\begin{aligned} \rho \, \left( \begin{array}{c} \mathbf{F}^{(0)} \\ \mathbf{F}^{(1)} \end{array} \right)= & {} \frac{1}{2 \pi \hbar ^2} \left( \begin{array}{c} I_1 (\lambda , \lambda _{W})\\ I_2 (\lambda , \lambda _{W}) \end{array} \right) \, \mathbf{I}, \end{aligned}$$
(16)
$$\begin{aligned} \rho \, \left( \begin{array}{c} \mathbf{G}^{(0)} \\ \mathbf{G}^{(1)} \end{array} \right)= & {} \frac{1}{2 \pi \hbar ^2} \left( \begin{array}{c} \mathcal {I}_0 (\lambda , \lambda _{W})\\ 2 I_1 (\lambda , \lambda _{W}) \end{array} \right) \, \mathbf{I} \end{aligned}$$
(17)

where \(\mathbf{I}\) is the identity matrix and

$$\begin{aligned} \mathcal {I}_0(\lambda , \lambda _W)= \int _{\frac{\pi \hbar v_F}{L}}^{+ \infty } \frac{1}{1+e^{\lambda + \lambda _W \varepsilon }} \left[ 1 + \left( \frac{\pi \hbar v_F}{L \varepsilon } \right) ^2 \right] . \end{aligned}$$

The moment equations form a hyperbolic system of balance laws in the realizability region \(\lambda _{W_A} >0\) [8, 9].

Regarding the production terms, they are given by a sum of contributions arising from the different types of phonon scattering

$$\begin{aligned} C_{\psi } = C_{\psi }^{ac} + \sum _{s = LO,TO, K} C_{\psi }^{s} + C_{\psi }^{(edge)}, \end{aligned}$$

with \(\psi = \rho , \mathbf{V}, W, \mathbf{S}\). In the next sections all the productions terms will be explicitly calculated completing the formulation of the model. The constitutive relations involve some integral functions which can be easily evaluated numerically.

5 Production Terms

5.1 Production Terms of Acoustic Phonon Scattering

The acoustic phonons scattering conserves density and energy, so

$$\begin{aligned} C_{\rho _A}^{(ac)} = C_{W_{A}}^{(ac)} = 0 \end{aligned}$$
(18)

while velocity and energy-flux production terms are given by

$$\begin{aligned} \begin{pmatrix} \rho C_\mathbf{V} \\ \rho C_\mathbf{S} \end{pmatrix} =-\frac{A^{(ac)}}{4 D \pi \hbar ^2 v_F^2} \begin{pmatrix} J_2 &{} J_3 \\ J_3 &{} J_4 \end{pmatrix} \begin{pmatrix} -J_3 &{} J_2 \\ J_2 &{} -J_1 \end{pmatrix} \begin{pmatrix} \rho \mathbf{V} \\ \rho \mathbf{S} \end{pmatrix}. \end{aligned}$$
(19)

5.2 Production Terms of LO and TO Optical Phonon Scattering

The production terms due to optical phonon scattering are the sum of two contributions, LO and TO. Each of them gives an intravalley and an intervalley contribution. Here we consider only intravalley contribution. The intraband scattering conserves the particle number

$$\begin{aligned} \rho _A C_{\rho _A}^{(opt)}=0 . \end{aligned}$$
(20)

The energy production term is given by

$$\begin{aligned} \rho C_W^{(opt)}=\frac{A^{(opt)} D^2_{\varGamma } N_B^{(opt)} }{\pi ^2 (\hbar v_F)^4} \hbar \omega _{opt} \left( e^{\lambda _W \hbar \omega _{opt}} - e^{\beta \hbar \omega _{opt}} \right) I_W(\lambda , \lambda _W) \end{aligned}$$
(21)

with

$$\begin{aligned} I_W(\lambda , \lambda _W)=\int _{\frac{\pi \hbar v_F}{L}}^{+\infty } \frac{\varepsilon (\varepsilon +\hbar \omega _{opt})}{1+e^{\lambda + \lambda _w \varepsilon }} \frac{e^{\lambda +\lambda _W \varepsilon }}{1+e^{\lambda +\lambda _W (\varepsilon + \hbar \omega _{opt})}} d \varepsilon . \end{aligned}$$

At equilibrium \(\lambda _W= \beta \) and therefore \(C_{W_A}^{(opt)}\) vanishes.

The other production terms read

$$\begin{aligned} \begin{pmatrix} \rho C_\mathbf{V}^{(opt)}\\ \rho C_\mathbf{S}^{(opt)} \end{pmatrix} =\frac{A^{(opt)} D^2\varGamma N_B^{(opt)}}{\pi (v_F \hbar )^2 D e^\lambda } \begin{pmatrix} Z_1+V_1 &{} Z_2+V_2 \\ Z_2+V_2 &{} Z_3+V_3 \end{pmatrix} \begin{pmatrix} -J_3 &{} J_2 \\ J_2 &{} -J_1 \end{pmatrix} \begin{pmatrix} \rho \mathbf{V}\\ \rho \mathbf{S} \end{pmatrix} \end{aligned}$$
(22)

with

$$\begin{aligned} Z_n=\int _{\frac{\pi \hbar v_F}{L}}^{+\infty } \frac{e^{\beta \hbar \omega _{opt}}+e^{\lambda +\lambda _W (\varepsilon + \hbar \omega _{opt})}}{1+e^{\lambda +\lambda _W (\varepsilon + \hbar \omega _{opt})}} \frac{e^{\lambda +\lambda _W \varepsilon }}{(1+e^{\lambda +\lambda _W \varepsilon })^2} \varepsilon ^n (\varepsilon + \hbar \omega _{opt}) \left[ 1-\left( \frac{\pi \hbar v_F}{L})^2 \right) \right] d \varepsilon \end{aligned}$$

and

$$\begin{aligned} V_n= & {} \int _{\frac{\pi \hbar v_F}{L}}^{+\infty } \frac{1+e^{\beta \hbar \omega _{opt}+\lambda +\lambda _W \varepsilon }}{1+e^{\lambda +\lambda _W \varepsilon }} \frac{e^{\lambda +\lambda _W (\varepsilon +\hbar \omega _{opt})}}{\left( 1+e^{\lambda +\lambda _W (\varepsilon + \hbar \omega _{s})}\right) ^2} \varepsilon (\varepsilon + \hbar \omega _{opt})^n \\&\sqrt{1-\left( \frac{\pi \hbar v_F}{L(\varepsilon +\hbar \omega )}\right) ^2} \sqrt{1-\left( \frac{\pi \hbar v_F}{L\varepsilon }\right) ^2} d \varepsilon . \end{aligned}$$

5.3 Production Terms of K-Phonon Scattering

For the production terms arising from the K-phonon scattering the same considerations as regarding the optical phonon hold. Again we consider only intraband scattering. One has

$$\begin{aligned}&\!\!\! \rho _A C_\rho ^{(K)} = 0\\&\!\!\! \rho C_W^{(K)}=\frac{A^{(K) D^2_K N_B^{(K)} }}{4 \pi v_F^2 \hbar ^2 D e^\lambda } \hbar \omega _K \left( e^{\lambda _W \hbar \omega _K }- e^{\beta \hbar \omega _K} \right) \int _{\frac{\pi \hbar v_F}{L}}^{+\infty } \frac{\varepsilon (\varepsilon +\hbar \omega _K) e^{\lambda +\lambda _W \varepsilon }}{\left( 1+e^{\lambda +\lambda _W \varepsilon } \right) \left( 1+e^{\lambda +\lambda _W (\varepsilon +\hbar \omega _K)} \right) } d \varepsilon \\&\!\!\! \begin{pmatrix} \rho C_\mathbf{V}^{(K)}\\ \rho C_\mathbf{S}^{(K)} \end{pmatrix} =\frac{A^{(K)} D^2_K N_B^{(K)} }{4 \pi v_F^2 \hbar ^2 D e^{\lambda }} \begin{pmatrix} D_{1\,1}(0,0)+ {\bar{D}}_{1\,1}(0,0) &{} D_{1\,2}(0,0)+ {\bar{D}}_{2\,1}(0,0) \\ D_{1\,1}(1,1)+ {\bar{D}}_{1\,1}(1,1) &{} D_{2\,1}(1,1)+ {\bar{D}}_{2\,1}(1,1) \end{pmatrix} \begin{pmatrix} -J_3 &{} J_2\\ J_2 &{} -J_1 \end{pmatrix} \begin{pmatrix} \rho \mathbf{V}\\ \rho \mathbf{S} \end{pmatrix} \end{aligned}$$

with

$$\begin{aligned}&D_{m\,n}(p,q,\lambda ,\lambda _W)=\int _{\frac{\pi v_F \hbar }{L}}^{+\infty } \varepsilon ^m (\varepsilon + \hbar \omega _K)^n \frac{1+e^{\beta \hbar \omega _K + \lambda +\lambda \varepsilon }}{1+e^{\lambda +\lambda _W \varepsilon }} \frac{e^{\lambda +\lambda _W(\varepsilon + \hbar \omega _K)}}{\left( 1+e^{\lambda +\lambda _W (\varepsilon +\hbar \omega _K)} \right) ^2}\\&\quad \left[ 2(\varepsilon +\hbar \omega _K)^p \left( 1-\left( \frac{\pi \hbar v_F}{L(\varepsilon +\hbar \omega _K)} \right) ^2 \right) + \varepsilon ^q \sqrt{1-\left( \frac{\pi \hbar v_F}{L(\varepsilon +\hbar \omega _K)} \right) ^2 } \sqrt{1-\left( \frac{\pi \hbar v_F}{L\varepsilon } \right) ^2 } \right] d \varepsilon ,\\&\quad {\bar{D}}_{m\,n}(p,q,\lambda ,\lambda _W)=\int _{\frac{\pi v_F \hbar }{L}}^{+\infty } \varepsilon ^m (\varepsilon + \hbar \omega _K)^n \frac{e^{\lambda +\lambda _W(\varepsilon +\hbar \omega _K)} +e^{\beta \hbar \omega _K} }{1+e^{\lambda +\lambda _W(\varepsilon +\hbar \omega _K)}} \frac{e^{\lambda +\lambda _W \varepsilon }}{\left( 1+e^{\lambda +\lambda _W \varepsilon } \right) ^2}\\&\quad \left[ 2 \varepsilon ^p \left( 1-\left( \frac{\pi \hbar v_F}{L(\varepsilon +\hbar \omega _K)} \right) ^2 \right) + (\varepsilon +\hbar \omega _K)^q \sqrt{1-\left( \frac{\pi \hbar v_F}{L(\varepsilon +\hbar \omega _K)} \right) ^2 } \sqrt{1-\left( \frac{\pi \hbar v_F}{L\varepsilon } \right) ^2 } \right] d \varepsilon . \end{aligned}$$

5.4 Production Terms of the Edge Scattering

The edge scattering conserves the particle numbers and the energy, so

$$\begin{aligned} \rho C_\rho ^{(edge)}= & {} 0\\ \rho C_W^{(edge)}= & {} 0 \end{aligned}$$

while the expressions of the production terms for the velocity \(\mathbf{V}\) and for the energy flux \(\mathbf{S}\) are

$$\begin{aligned} \begin{pmatrix} \rho C_\mathbf{V}^{(edge)}\\ \rho C_\mathbf{S}^{(edge)} \end{pmatrix} =\frac{N_i V_0^2}{\hbar ^3 v_F^2 L D e^{\lambda }} \begin{pmatrix} ML_0(\lambda , \lambda _W) &{} ML_1(\lambda , \lambda _W) \\ ML_1(\lambda , \lambda _W) &{} ML_2(\lambda , \lambda _W) \end{pmatrix} \begin{pmatrix} -J_3 &{} J_2\\ J_2 &{} -J_1 \end{pmatrix} \begin{pmatrix} \rho \mathbf{V}\\ \rho \mathbf{S} \end{pmatrix} \end{aligned}$$

with

$$\begin{aligned}&ML_n(\lambda ,\lambda _W)\\&\quad =\int _0^{2 \pi } d \theta ' \int _0^{2 \pi } d \theta \int _{\frac{\pi \hbar v_F}{L}}^{+\infty } d \varepsilon \, \varepsilon ^n \left( \varepsilon ^2 - \left( \frac{\pi \hbar v_F}{L} \right) ^2 \right) e^{-\left( \frac{a}{\hbar v_F} \right) ^2 \left( \varepsilon ^2 - \left( \frac{\pi \hbar v_F}{L} \right) ^2 \right) (\cos \theta -\cos \theta ')^2 } \\&(\cos \theta - \cos \theta ') cos \theta \frac{e^{\lambda + \lambda _W \varepsilon }}{\left( 1+e^{\lambda +\lambda _W \varepsilon } \right) ^2} \end{aligned}$$

where \(\theta \) and \(\theta '\) are the angles between \(\mathbf{k}\) and the edge and \(\mathbf{k}'\) and edge respectively.

6 Simulation of Charge Transport in the Homogeneous Case

To test the model we consider a simple case with a constant electric field \(\mathbf{E}\) along the x direction so the distribution is \(f(\mathbf{k},t)\) and equations (5) became

$$\begin{aligned}&\frac{\partial }{\partial t}&\rho = \rho C_{\rho } \\&\quad \frac{\partial }{\partial t}&(\rho W) -e \rho \mathbf{E} \cdot \mathbf{V}= \rho C_W \\&\quad \frac{\partial }{\partial t}&(\rho \mathbf{V}) - e \rho \mathbf{G}^{(0)} \cdot \mathbf{E} = \rho C_\mathbf{V}\\&\quad \frac{\partial }{\partial t}&(\rho \mathbf{S}) - e \rho \mathbf{G}^{(1)} \cdot \mathbf{E} = \rho S_\mathbf{V} \end{aligned}$$

or in array form

$$\begin{aligned} \frac{\partial }{\partial t} \begin{pmatrix} \rho \\ \rho W \\ \rho \mathbf{V} \\ \rho \mathbf{S} \end{pmatrix} -e \mathbf{E} \begin{pmatrix} 0 \\ \rho \mathbf{V} \\ \rho \mathbf{G}^{(0)} \\ \rho \mathbf{G}^{(1)} \end{pmatrix} = \begin{pmatrix} \rho C _{\rho } \\ \rho C _W \\ \rho C _\mathbf{V} \\ \rho C _\mathbf{V} \end{pmatrix} . \end{aligned}$$

To reduce the numerical cost we have considered as independent variables \(\lambda , \lambda _W, \rho \mathbf{V}, \rho \mathbf{S}\) so the evolution equations are

$$\begin{aligned} \mathbf{M} \frac{d}{d t} \begin{pmatrix} \lambda \\ \lambda _W \\ \rho \mathbf{V} \\ \rho \mathbf{S} \end{pmatrix} =e \mathbf{E} \begin{pmatrix} 0 \\ \rho \mathbf{V} \\ \rho \mathbf{G}^{(0)} \\ \rho \mathbf{G}^{(1)} \end{pmatrix} + \begin{pmatrix} \rho C _{\rho } \\ \rho C _W \\ \rho C _\mathbf{V} \\ \rho C _\mathbf{V} \end{pmatrix} \end{aligned}$$

with

$$\begin{aligned} \mathbf{M}=\begin{pmatrix} -\frac{P_1(\lambda ,\lambda _W)}{\pi (\hbar v_F)^2} &{} -\frac{P_2(\lambda ,\lambda _W)}{\pi (\hbar v_F)^2} &{} 0 &{} 0\\ -\frac{P_2(\lambda ,\lambda _W)}{\pi (\hbar v_F)^2} &{} -\frac{P_3(\lambda ,\lambda _W)}{\pi (\hbar v_F)^2} &{} 0 &{} 0\\ 0 &{} 0 &{} 1 &{} 0\\ 0 &{} 0 &{} 0 &{} 1 \end{pmatrix} \quad \text {and} \quad P_n=\int _{\frac{\pi \hbar v_F}{L}}^{+\infty } \frac{\varepsilon ^n e^{\lambda +\lambda _W \varepsilon }}{(1+e^{\lambda +\lambda _W \varepsilon })^2} d \varepsilon . \end{aligned}$$

In the simulation here presented we use an applied electric field of \(0.4 \, V/\mu m\), a Fermi energy of \(0.4 \, eV\) and several values of the width of the ribbon. The results are shown in Fig. 2 and in Fig. 3.

From the top to the bottom of Fig. 2 the mean energy (left) and the mean velocity (right) are plotted for increasing values of the width of the ribbon. The plots show that the edge scattering reduces the mean velocity and the mean energy of the charges and that this degradation disappears by increasing the width of the ribbon. This can be ascribed to the boundary effects, being more relevant if the ribbon is narrow. As consequence, the electron mobility is reduced, leading to worse material performances, despite to the presence of an energy gap, leading to a better behavior of electron device. This is a key point in the design of new field effect transistors made of graphene where the efficiency is related to the electron mobility and the absence of a gap is a drawback because there exists an inversion current due to the holes. The stationary state is reached after 1.0–1.5 ps.

Figure 3 shows the increase of the density current for unit length by the increase of the width of the ribbon. For larger ribbons the density current is higher because the channel is larger.

Fig. 2
figure 2

Mean energy value (left) and mean velocity (right) with an applied electric field of \(0.4 \, V/\mu m\) and a Fermi energy of \(0.4 \, eV\) for different width of the nanoribbon (greater to the bottom): 7 nm, 10 nm, 100 nm, 1000 nm. It is remarkable the degradation of the mean velocity and of the mean energy due to the edge scattering. The effect decreases for higher values of the width and disappears in the limit case of graphene sheet

Full size image
Fig. 3
figure 3

Current density for unit length versus time. The applied electric field is \(0.4 \, V/\mu m\) and the Fermi energy \(0.4 \, eV\)

Full size image

7 Conclusions

The hydrodynamical model for charge transport in graphene nanoribbons here formulated takes into account the energy gap due to the charge confinement and the scattering of the electrons with the lattice as well as the scattering with the edge of the ribbon by supposing an high degree of edge disorder and assuming the energy band model proposed in [4]. The evolution equations have been obtained by taking suitable moments of the semiclassical transport equations and getting the needed closure relations by employing the maximum entropy principle. The model is free from adjustable parameters, apart from those already present at kinetic level.

The simulation has shown that the model describes qualitatively the macroscopic behavior of the charge transport and is able to differentiate the graphene ribbons from the graphene sheets. Moreover the obtained results are comparable with that ones obtained by solving numerically the Boltzmann equation [22, 23, 25, 26] but with a remarkable reduction of the computational time.

The code for the simulation is available upon request.