A BGK model for charge transport in graphene

TL;DR

This paper introduces and analyzes a BGK kinetic model tailored for charge transport phenomena in graphene, aiming to simplify the complex Boltzmann equation while capturing essential microscopic behaviors.

Contribution

The paper proposes a novel BGK model specifically designed for charge transport in graphene, extending kinetic theory applications to this material.

Findings

The model effectively approximates charge transport in graphene.

Analytical and numerical analysis demonstrates the model's accuracy.

The approach simplifies computations compared to the full Boltzmann equation.

Abstract

The Boltzmann equation describes the detailed microscopic behaviour of a dilute gas, and represents the basis of the kinetic theory of gases. In order to reduce the difficulties in solving the Boltzmann equation, simple expressions of a collision operator have been proposed to replace the true Boltzmann integral term. These new equations are called kinetic models. The most popular and widely used kinetic model is the Bhatnagar-Gross-Krook (BGK) model. In this work we propose and analyse a BGK model for charge transport in graphene.

Equations

$$\frac{\partial f}{\partial t}+\frac{1}{\hslash }{\nabla }_{\mathbf{k}}\epsilon \cdot {\nabla }_{\mathbf{x}}f-\frac{e}{\hslash }\mathbf{E}\cdot {\nabla }_{\mathbf{k}}f={\left(\frac{df}{dt}\right)}_{coll},$$

$${\int }_{{\mathbb{R}}^2}S({\mathbf{k}}^{\prime },\mathbf{k})f(t,\mathbf{x},{\mathbf{k}}^{\prime })\left[1-f(t,\mathbf{x},\mathbf{k})\right]d{\mathbf{k}}^{\prime }-{\int }_{{\mathbb{R}}^2}S(\mathbf{k},{\mathbf{k}}^{\prime })f(t,\mathbf{x},\mathbf{k})\left[1-f(t,\mathbf{x},{\mathbf{k}}^{\prime })\right]d{\mathbf{k}}^{\prime },$$

$$\mathbb{F}(\mathbf{k},\mu )=\frac{1}{1+\exp \left(\frac{\epsilon (\mathbf{k})-\mu }{k_{B}T}\right)},$$

$$Q(f)(t,\mathbf{x},\mathbf{k})={\int }_{{\mathbb{R}}^2}S({\mathbf{k}}^{\prime },\mathbf{k})\mathbb{F}({\mathbf{k}}^{\prime },\mu )\left[1-f(t,\mathbf{x},\mathbf{k})\right]d{\mathbf{k}}^{\prime }-{\int }_{{\mathbb{R}}^2}S(\mathbf{k},{\mathbf{k}}^{\prime })f(t,\mathbf{x},\mathbf{k})\left[1-\mathbb{F}({\mathbf{k}}^{\prime },\mu )\right]d{\mathbf{k}}^{\prime }.$$

$${\int }_{{\mathbb{R}}^2}Q(f)(t,\mathbf{x},\mathbf{k})d\mathbf{k}=0.$$

$${\Phi }_0(\mathbf{k},\mu )={\int }_{{\mathbb{R}}^2}S({\mathbf{k}}^{\prime },\mathbf{k})\mathbb{F}({\mathbf{k}}^{\prime },\mu )d{\mathbf{k}}^{\prime }\text{mbox}and{\Phi }_1(\mathbf{k},\mu )={\int }_{{\mathbb{R}}^2}S(\mathbf{k},{\mathbf{k}}^{\prime })\left[1-\mathbb{F}({\mathbf{k}}^{\prime },\mu )\right]d{\mathbf{k}}^{\prime },$$

$$Q(f)(t,\mathbf{x},\mathbf{k})={\Phi }_0(\mathbf{k},\mu )\left[1-f(t,\mathbf{x},\mathbf{k})\right]-{\Phi }_1(\mathbf{k},\mu )f(t,\mathbf{x},\mathbf{k}).$$

$${\Phi }_0(\mathbf{k},\mu )\left[1-\mathbb{F}(\mathbf{k},\mu )\right]-{\Phi }_1(\mathbf{k},\mu )\mathbb{F}(\mathbf{k},\mu )$$

$$={\int }_{{\mathbb{R}}^2}S({\mathbf{k}}^{\prime },\mathbf{k})\mathbb{F}({\mathbf{k}}^{\prime },\mu )\left[1-\mathbb{F}(\mathbf{k},\mu )\right]d{\mathbf{k}}^{\prime }-{\int }_{{\mathbb{R}}^2}S(\mathbf{k},{\mathbf{k}}^{\prime })\mathbb{F}(\mathbf{k},\mu )\left[1-\mathbb{F}({\mathbf{k}}^{\prime },\mu )\right]d{\mathbf{k}}^{\prime }=0.$$

$$Q(f)(t,\mathbf{x},\mathbf{k})$$

$$\kappa (\mathbf{k},\mu )=\frac{{\Phi }_0(\mathbf{k},\mu )}{\mathbb{F}(\mathbf{k},\mu )}={\Phi }_0(\mathbf{k},\mu )\left[1+\exp \left(\frac{\epsilon (\mathbf{k})-\mu }{k_{B}T}\right)\right],$$

$$\frac{\partial f}{\partial t}(t,\mathbf{x},\mathbf{k})+\frac{1}{\hslash }{\nabla }_{\mathbf{k}}\epsilon \cdot {\nabla }_{\mathbf{x}}f(t,\mathbf{x},\mathbf{k})-\frac{e}{\hslash }\mathbf{E}\cdot {\nabla }_{\mathbf{k}}f(t,\mathbf{x},\mathbf{k})=\kappa (\mathbf{k},\mu )\left[\mathbb{F}(\mathbf{k},\mu )-f(t,\mathbf{x},\mathbf{k})\right],$$

$${\int }_{{\mathbb{R}}^2}\kappa (\mathbf{k},\mu )\left[\mathbb{F}(\mathbf{k},\mu )-f(t,\mathbf{x},\mathbf{k})\right]d\mathbf{k}=0.$$

$$\frac{\partial f}{\partial t}(t,\mathbf{k})=\kappa (\mathbf{k},\mu )\left[\mathbb{F}(\mathbf{k},\mu )-f(t,\mathbf{k})\right]\text{mbox}and{\int }_{{\mathbb{R}}^2}\kappa (\mathbf{k},\mu )\left[\mathbb{F}(\mathbf{k},\mu )-f(t,\mathbf{k})\right]d\mathbf{k}=0,$$

$$K(\mathbf{k},\mu )\text{mbox}suchthat\frac{\partial K}{\partial t}=\kappa (\mathbf{k},\mu )\forall t\ge 0.$$

$$\exp (K(\mathbf{k},\mu ))\frac{\partial f}{\partial t}=\exp (K(\mathbf{k},\mu ))\kappa (\mathbf{k},\mu )\mathbb{F}(\mathbf{k},\mu )-\exp (K(\mathbf{k},\mu ))\kappa (\mathbf{k},\mu )f(t,\mathbf{k})$$

$$\frac{\partial \text{mbox}}{\partial t}\left[\exp (K(\mathbf{k},\mu ))f(t,\mathbf{k})\right]=\exp (K(\mathbf{k},\mu ))\kappa (\mathbf{k},\mu )\mathbb{F}(\mathbf{k},\mu )\ge 0.$$

$$\exp (K(\mathbf{k},\mu ))f(t,\mathbf{k})\ge {\left[\frac{}{}\exp (K(\mathbf{k},\mu ))f(t,\mathbf{k})\right|}_{t=0}.$$

$$\frac{\partial g}{\partial t}=\kappa (\mathbf{k},\mu )\left[\left(1-\mathbb{F}(\mathbf{k},\mu )\right)-g\right]$$

$$\exp (K(\mathbf{k},\mu ))g(t,\mathbf{k})\ge {\left[\frac{}{}\exp (K(\mathbf{k},\mu ))g(t,\mathbf{k})\right|}_{t=0},$$

$${\left|G^{(\nu )}(\mathbf{k},{\mathbf{k}}^{\prime })\right|}^2\left[\left(n_{\mathbf{q}}^{(\nu )}+1\right)\delta \left(\epsilon ({\mathbf{k}}^{\prime })-\epsilon (\mathbf{k})+\hslash {\omega }_{\mathbf{q}}^{(\nu )}\right)+n_{\mathbf{q}}^{(\nu )}\delta \left(\epsilon ({\mathbf{k}}^{\prime })-\epsilon (\mathbf{k})-\hslash {\omega }_{\mathbf{q}}^{(\nu )}\right)\right].$$

$$n_{\mathbf{q}}^{(\nu )}=\frac{1}{e^{\hslash {\omega }_{\mathbf{q}}^{(\nu )}/k_{B}T}-1}.$$

$$2n_{\mathbf{q}}^{(ac)}{\left|G^{(ac)}({\mathbf{k}}^{\prime },\mathbf{k})\right|}^2=\frac{1}{(2\pi {)}^2}\frac{\pi D_{ac}^2{k}_{B}T}{2\hslash {\sigma }_m{v}_p^2}\left(1+\cos {\vartheta }_{\mathbf{k},{\mathbf{k}}^{\prime }}\right),$$

$${\left|G^{(LO)}({\mathbf{k}}^{\prime },\mathbf{k})\right|}^2$$

$${\left|G^{(TO)}({\mathbf{k}}^{\prime },\mathbf{k})\right|}^2$$

$${\left|G^{(K)}({\mathbf{k}}^{\prime },\mathbf{k})\right|}^2$$

$${\Phi }_0(\mathbf{k},\mu )\left\{1-f(t,\mathbf{x},\mathbf{k})\left[1+\exp \left(\frac{\epsilon (\mathbf{k})-\mu }{k_{B}T}\right)\right]\right\}.$$

$${\int }_{{\mathbb{R}}^2}{\mathcal{G}}^{(\nu )}(\cos {\vartheta }_{\mathbf{k},{\mathbf{k}}^{\prime }})\left[\left(n_{\mathbf{q}}^{(\nu )}+1\right)\delta \left(\epsilon ({\mathbf{k}}^{\prime })-\epsilon (\mathbf{k})-\hslash {\omega }_{\mathbf{q}}^{(\nu )}\right)+n_{\mathbf{q}}^{(\nu )}\delta \left(\epsilon ({\mathbf{k}}^{\prime })-\epsilon (\mathbf{k})+\hslash {\omega }_{\mathbf{q}}^{(\nu )}\right)\right]\mathbb{F}({\mathbf{k}}^{\prime },\mu )d{\mathbf{k}}^{\prime }$$

$$\text{mbox}=\frac{n_{\mathbf{q}}^{(\nu )}+1}{1+\exp \left(\frac{\epsilon (\mathbf{k})+\hslash {\omega }_{\mathbf{q}}^{(\nu )}-\mu }{k_{B}T}\right)}{\int }_{{\mathbb{R}}^2}{\mathcal{G}}^{(\nu )}(\cos {\vartheta }_{\mathbf{k},{\mathbf{k}}^{\prime }})\delta \left(\epsilon ({\mathbf{k}}^{\prime })-\epsilon (\mathbf{k})-\hslash {\omega }_{\mathbf{q}}^{(\nu )}\right)d{\mathbf{k}}^{\prime }$$

$$\text{mbox}+\frac{n_{\mathbf{q}}^{(\nu )}}{1+\exp \left(\frac{{\left[\epsilon (\mathbf{k})-\hslash {\omega }_{\mathbf{q}}^{(\nu )}\right]}_{+}-\mu }{k_{B}T}\right)}{\int }_{{\mathbb{R}}^2}{\mathcal{G}}^{(\nu )}(\cos {\vartheta }_{\mathbf{k},{\mathbf{k}}^{\prime }})\delta \left(\epsilon ({\mathbf{k}}^{\prime })-\epsilon (\mathbf{k})+\hslash {\omega }_{\mathbf{q}}^{(\nu )}\right)d{\mathbf{k}}^{\prime },$$

$${\int }_{{\mathbb{R}}^2}{\mathcal{G}}^{(\nu )}(\cos {\vartheta }_{\mathbf{k},{\mathbf{k}}^{\prime }})\delta \left(\epsilon ({\mathbf{k}}^{\prime })-\epsilon (\mathbf{k})-\hslash {\omega }_{\mathbf{q}}^{(\nu )}\right)d{\mathbf{k}}^{\prime }$$

$$\text{mbox}={\int }_0^{+\infty }d{r}^{\prime }{\int }_0^{2\pi }d{\vartheta }^{\prime }{\mathcal{G}}^{(\nu )}(\cos ({\vartheta }^{\prime }-\vartheta ))\delta \left(\hslash v_F{r}^{\prime }-\hslash v_{F}r-\hslash {\omega }_{\mathbf{q}}^{(\nu )}\right)r^{\prime }$$

$$\text{mbox}=\left[{\int }_0^{2\pi }{\mathcal{G}}^{(\nu )}(\cos ({\vartheta }^{\prime }-\vartheta ))d{\vartheta }^{\prime }\right]\frac{1}{\hslash v_F}{\int }_{\mathbb{R}}\delta \left(r^{\prime }-r-\frac{{\omega }_{\mathbf{q}}^{(\nu )}}{v_F}\right)r^{\prime }H(r^{\prime })d{r}^{\prime }$$

$$\text{mbox}=\left[{\int }_0^{2\pi }{\mathcal{G}}^{(\nu )}(\cos {\vartheta }^{\prime })d{\vartheta }^{\prime }\right]\frac{1}{\hslash v_F}\left(r+\frac{{\omega }_{\mathbf{q}}^{(\nu )}}{v_F}\right)H\left(r+\frac{{\omega }_{\mathbf{q}}^{(\nu )}}{v_F}\right)=\left[{\int }_0^{2\pi }{\mathcal{G}}^{(\nu )}(\cos {\vartheta }^{\prime })d{\vartheta }^{\prime }\right]\frac{\epsilon (\mathbf{k})+\hslash {\omega }_{\mathbf{q}}^{(\nu )}}{(\hslash v_F{)}^2},$$

$${\int }_{{\mathbb{R}}^2}{\mathcal{G}}^{(\nu )}(\cos {\vartheta }_{\mathbf{k},{\mathbf{k}}^{\prime }})\delta \left(\epsilon ({\mathbf{k}}^{\prime })-\epsilon (\mathbf{k})+\hslash {\omega }_{\mathbf{q}}^{(\nu )}\right)d{\mathbf{k}}^{\prime }=\left[{\int }_0^{2\pi }{\mathcal{G}}^{(\nu )}(\cos {\vartheta }^{\prime })d{\vartheta }^{\prime }\right]\frac{{\left[\epsilon (\mathbf{k})-\hslash {\omega }_{\mathbf{q}}^{(\nu )}\right]}_{+}}{(\hslash v_F{)}^2}.$$

Full text

In this work we propose and analyse a BGK model for charge transport in graphene.

MSC-class: 82C40 - 82C70 (Primary) 82D37 (Secondary)

1 Introduction

Graphene is a gapless semiconductor made of a sheet composed of a single layer of carbon atoms arranged into a honeycomb hexagonal lattice [1]. In view of application in graphene-based electron devices, it is crucial to understand the basic transport properties of this material. In a semiclassical kinetic setting, the charge transport in graphene is described by four Boltzmann equations, one for electrons in the valence ($\pi$) band and one for electrons in the conductions ($\pi^{*}$) band, that in turn can belong to the $K$ or $K^{\prime}$ valley. In this paper we study the case of a single distribution function $f$ for electrons belonging to a conduction band. This corresponds to a physical case, where a n-type doping or equivalently a high value of the Fermi potential is considered, and the electrons, belonging to a conduction band, do not move to the valence band. Moreover $K$ and $K^{\prime}$ are considered equivalent. Under these assumptions Boltzmann equation writes

[TABLE]

where $\mathrm{t}$ is the time, $\mathbf{x}$ and $\mathbf{k}\in\mathbb{R}^{2}$ are the position and the wave-vector of a charge particle, respectively. We denote by $\nabla_{{\mathbf{x}}}$ and $\nabla_{\mathbf{k}}$ the gradients with respect to the position and wave-vector, respectively. With a very good approximation [1] the energy bands is given by $\varepsilon(\mathbf{k})=\hbar\,v_{F}\left|\mathbf{k}\right|$, where $v_{F}$ is the (constant) Fermi velocity, and $\hbar$ is the Planck constant divided by $2\,\pi$. The elementary (positive) charge is denoted by $e$, and ${\bf E}$ is the electric field. The collision operator of Equation (1) describes the interaction of electrons with acoustic, optical and $K$ phonons.

If the electric field ${\bf E}$ is not constant, then it must be self-consistently evaluated by coupling the Boltzmann equation (1) with the Poisson equation for the electrostatic potential.

More recently accurate numerical solutions, based on the discontinuous Galerkin method, to the Boltzmann equation have been shown [2], [3]. A comparison with results, obtained by using a direct simulation Monte Carlo approach, shows an excellent agreement, which gives a further validation of the numerical scheme.

In this work the collision operator of Equation (1) is replaced by a relaxation (BGK) collision operator. We propose an operator, which retains the fundamental properties of the Boltzmann equation, such the mass conservation, the same equilibrium distribution functions and it properly deals with Pauli’s exclusion principle in the degenerate case.

2 The BGK model

The collision operator of Equation (1) writes [4]

[TABLE]

where the transition rate $S(\mathbf{k},\mathbf{k}^{\prime})$, related to electron-phonon scatterings, is described in detail in the following.

The collision operator (2) vanishes if $f$ is the Fermi-Dirac distribution

[TABLE]

where $k_{B}$ is the Boltzmann constant, $T$ the constant graphene lattice temperature, and $\mu$ is the chemical potential. If the electric field $\bf E$ is null, then Fermi-Dirac distributions are solutions, which do not depend on time $t$ and space coordinates $\mathbf{x}$, of the Boltzmann equation (1).

In this paper we propose the following BGK collision operator

[TABLE]

It is derived from Equation (2), replacing the distribution $f(t,\mathbf{x},\mathbf{k}^{\prime})$, inside the integrals, with a Fermi-Dirac distribution. Now the chemical potential $\mu$ becomes a new unknown. Therefore we must add an equation to the kinetic model. It is

[TABLE]

This equation guarantees the mass conservation, exactly as the same for the Boltzmann equation.

If we define

[TABLE]

then (3) becomes

[TABLE]

It is evident that, for every $\mathbf{k}$ and $\mu$,

[TABLE]

Using this identity, we have

[TABLE]

By defining the collision frequency of the BGK model

[TABLE]

the kinetic model writes

[TABLE]

A solution of system (9)-(10) consists of the two functions $f(t,\mathbf{x},\mathbf{k})$ and $\mu(t,\mathbf{x})$, which must satisfy suitable regularity conditions. The common definitions of the kernel $S$ guarantees that the function $\kappa(\mathbf{k},\mu)$ is a no negative continuous function in the set $\mathbb{R}^{2}\times\mathbb{R}$.

2.1 The Pauli exclusion principle

The distribution function $f$ must be not negative and less or equal to one, according to the Pauli’s exclusion principle. We prove that, if the electric field is null, then every spatial homogeneous solution of Equations (9)-(10) satisfies these conditions.

Theorem 1

Let be $f(t,\mathbf{k})$ a function defined in $[0,+\infty[\times\mathbb{R}^{2}$, differentiable with respect to $t$ for every $\mathbf{k}\in\mathbb{R}^{2}$, and integrable, with respect to $\mathbf{k}$, over the domain $\mathbb{R}^{2}$ for each time $t$. Moreover the function $\mu(t)$ is continuous in $[0,+\infty[$. If $f$ and $\mu$ satisfy the equations

[TABLE]

and $0\leq f(0,\mathbf{k})\leq 1$, for all $\mathbf{k}\in\mathbb{R}^{2}$, then $0\leq f(t,\mathbf{k})\leq 1$ for $(t,\mathbf{k})\in[0,+\infty[\times\mathbb{R}^{2}$.

Proof. Let be

[TABLE]

Therefore we have

[TABLE]

that is

[TABLE]

and then

[TABLE]

Hence the distribution function $f(t,\mathbf{k})$ is always not negative, because $f(0,\mathbf{k})\geq 0$ for every $\mathbf{k}\in\mathbb{R}^{2}$.

Since $f(0,\mathbf{k})\leq 1$ for every $\mathbf{k}\in\mathbb{R}^{2}$, we define $g(t,\mathbf{k})=1-f(t,\mathbf{k})$, and we note that $0\leq g(0,\mathbf{k})\leq 1$. It is evident that $g$ satisfies the equation

[TABLE]

and then, as before, taking into account that ${\mathbb{F}}(\mathbf{k},\mu)<1$, from

[TABLE]

it follows $g(t,\mathbf{k})\geq 0\Leftrightarrow f(t,\mathbf{k})\leq 1$. $\Box$

2.2 The kernel of the Boltzmann collision operator

The kernel $S(\mathbf{k},\mathbf{k}^{\prime})$ is given by the sum of terms of the kind

[TABLE]

The index $\nu$ labels the $\nu$th phonon mode, $G^{(\nu)}(\mathbf{k}^{\prime},\mathbf{k})$ is the scattering rate, which describes the scattering mechanism between phonons $\nu$ and electrons. The symbol $\delta$ denotes the Dirac distribution function, $\omega^{(\nu)}_{\mathbf{q}}$ the $\nu$th phonon frequency, and $n^{(\nu)}_{\mathbf{q}}$ is the Bose-Einstein distribution for the phonon of type $\nu$

[TABLE]

For acoustic phonons, usually one considers the elastic approximation, and

[TABLE]

where $D_{ac}$ is the acoustic phonon coupling constant, $v_{p}$ is the sound speed in graphene, $\sigma_{m}$ the graphene areal density, and $\vartheta_{\mathbf{k}\,,\mathbf{k}^{\prime}}$ is the convex angle between $\mathbf{k}$ and ${\mathbf{k}^{\prime}}$.

There are three relevant optical phonon scatterings: the longitudinal optical (LO), the transversal optical (TO) and the ${K}$ (K) phonons. The scattering rates are

[TABLE]

where $D_{O}$ is the optical phonon coupling constant, $\omega_{O}$ the optical phonon frequency, $D_{K}$ is the K-phonon coupling constant and $\omega_{K}$ the K-phonon frequency. The angles $\vartheta_{\mathbf{k}\,,\mathbf{k}^{\prime}-\mathbf{k}}$ and $\vartheta_{\mathbf{k}^{\prime}\,,\mathbf{k}^{\prime}-\mathbf{k}}$ denote the convex angles between $\mathbf{k}$ and $\mathbf{k}^{\prime}-\mathbf{k}$ and between $\mathbf{k}^{\prime}$ and $\mathbf{k}^{\prime}-\mathbf{k}$, respectively.

We used the same physical parameters of the Table 1 of Ref. [3].

Since the phonon frequency of longitudinal and transversal optical phonons coincide, then, as we sum the corresponding terms (11), the function $\cos(\vartheta_{\mathbf{k}\,,\mathbf{k}^{\prime}-\mathbf{k}}+\vartheta_{\mathbf{k}^{\prime}\,,\mathbf{k}^{\prime}-\mathbf{k}})$ can be eliminated, easily. Therefore, in this case, $S$ is a sum of terms (11) where now $\left|G^{(\nu)}(\mathbf{k},\mathbf{k}^{\prime})\right|^{2}={\cal G}^{(\nu)}(\cos\vartheta_{\mathbf{k}\,,\mathbf{k}^{\prime}})$.

2.3 The collision operator of the BGK model

Taking into account Equation (9) and the definition (8), the collision operator of the BGK model writes

[TABLE]

This expression can be simplified. To this aim, we consider the function $\Phi_{0}(\mathbf{k},\mu)$, which is given by the sum of the following integrals

[TABLE]

where $[z]_{+}=\max\{z,0\}$.

Introducing polar coordinates $\mathbf{k}=r(\cos\vartheta,\sin\vartheta)$ and $\mathbf{k}^{\prime}=r^{\prime}(\cos\vartheta^{\prime},\sin\vartheta^{\prime})$, we have

[TABLE]

where $H$ is the Heaviside function. Analogously

[TABLE]

Hence

[TABLE]

In order to simplify the notation, we define

[TABLE]

So

[TABLE]

Then the BGK collision operator (12) is a linear combination, with constant positive coefficients, of the terms

[TABLE]

Since $n^{(\nu)}_{\mathbf{q}}\geq 0$ and $\hbar\,\omega^{(\nu)}_{\mathbf{q}}\geq 0$, then we can consider the function $\psi$, defined in (13), only for $a\geq 0$, $b\geq 0$, and $\xi>0$.

Lemma 1

If $0\leq f(t,\mathbf{x},\mathbf{k})\leq 1$, then the collision operator (12) is a strictly increasing function of $\mu$ on $\mathbb{R}$.

Proof. Since Equation (12) is a linear combination, with constant positive coefficients, of the functions (15), it is sufficient to prove that every function (15) is strictly increasing on $\mathbb{R}$ with respect to the variable $\mu$.

Again to simplify the notation, we consider the general application

[TABLE]

where $0\leq\varphi(\mathbf{k})\leq 1$ for all $\mathbf{k}\in\mathbb{R}^{2}$. Equation (16) gives the expression (15), choosing the parameters $a$, $b$, $\xi$ appropriately, and $\varphi(\mathbf{k})=f(t,\mathbf{x},\mathbf{k})$, for fixed $t$ and $\mathbf{x}$. If we define

[TABLE]

then

[TABLE]

where we have taken into account the identity

[TABLE]

Hence we can write

[TABLE]

which is strictly decreasing with respect to the positive variable $\xi$, because $0\leq\varphi(\mathbf{k})\leq 1$.

Now the result follows for $\displaystyle\xi=\exp\left(\dfrac{-\,\mu}{k_{B}\,T}\right)$. $\Box$

It is useful, for the following, to establish some inequalities. From

[TABLE]

we obtain

[TABLE]

and

[TABLE]

We remark that the variable $\xi$ has not been involved in inequalities.

3 The mass conservation

Equation (10), that guarantees the conservation of mass, can be written as

[TABLE]

i.e.

[TABLE]

Theorem 2

If $0\leq f(t,\mathbf{x},\mathbf{k})\leq 1$ and $\varepsilon(\mathbf{k})\,f(t,\mathbf{x},\mathbf{k})$ is integrable with respect to $\mathbf{k}$ over $\mathbb{R}^{2}$, for all $t$ and $\mathbf{x}$, then there exists a unique $\mu$ satisfying Equation (19).

Proof. If there exists a solution $\mu$ of Equation (19), then it must be unique due to the Lemma (1). To prove the theorem, we show that every integrals of Equation (19) is negative for $\mu\rightarrow-\infty$, and positive for $\mu\rightarrow+\infty$. We do not consider the meaningless case $f(t,\mathbf{x},\mathbf{k})=0$ almost everywhere. Taking into account the inequality (17), we have

[TABLE]

where the last integral is finite due to hypotheses of the theorem. If $\mu\rightarrow-\infty$, then $\xi\rightarrow+\infty$, and

[TABLE]

This implies that the l.h.s. of Equation (19) is negative for $\mu\rightarrow-\infty$.

Now we consider the inequality (18). We have

[TABLE]

where both integrals exist. Since, for $\mu>0$, we have

[TABLE]

Therefore

[TABLE]

This implies that the l.h.s. of Equation (19) is positive for $\mu\rightarrow+\infty$, and it concludes the proof. $\Box$ This theorem establishes that the kinetic model (9)-(10) is correctly-set.

4 Conclusions

In this paper we propose a BGK model for charge transport in graphene. The collision operator of this model replaces the true non linear Boltzmann integral operator. A further equation, for a new unknown, is added to the kinetic equation in order to guarantee the mass conservation. The model would furnish good results when the electric field is not strong, so that the distribution function remains near an equilibrium Fermi-Dirac distribution. Moreover the simple collision term allows analytical investigations, which may be prohibitive for the Boltzmann equation. A numerical scheme, based on a discontinuous Galerkin method, for finding approximate solutions to the model, does not seem more simple than the scheme used for solving the Boltzmann equation in Ref. [5], due to the non-linearity of the integral equation for the mass conservation.

Acknowledgements

The author acknowledges the financial support provided by the project Modellistica, simulazione e ottimizzazione del trasporto di cariche in strutture a bassa dimensionalit , (2016-2018) University of Catania.