The Pauli principle and the Monte Carlo Method for charge transport in graphene

The attempt to include the Pauli principle in the Monte Carlo method by acting also on the free flight step and not only at the end of each collision is investigated. The charge transport in suspended monolayer graphene is considered as test case. The results are compared with those obtained in the standard Ensemble Monte Carlo technique and in the new Direct Simulation Monte Carlo algorithm which is able to correctly handle with Pauli's principle. The physical aspects of the investigated approach are analyzed as well.


INTRODUCTION
The Monte Carlo approach is nowadays of widespread use in the simulations of charge transport in semiconductor devices. Most of the methods employed in this field are based on the Ensemble Monte Carlo (EMC) procedure of Lebwohl and Price [1], and on the approach developed by Bosi, Jacoboni and Reggiani [2]- [3]. When the simulations involve degenerate materials, the inclusion of the Pauli principle becomes essential; to this aim, Lugli and Ferry [4] improved the EMC method and included the Pauli principle by means of a rejection technique at the end of each scattering process. Unfortunately, with this procedure the charge distribution can exceed the maximum value of one, leading to unphysical results. Some attempts were made in the literature at overcoming this difficulty: in [5,6] ad hoc scattering out terms to force the distribution function to be smaller than one were introduced; in [7,8] some approximations of the distribution of the final states were used. These efforts improved the simulation results but did not lead to a correct reconstruction of the distribution function, that is of fundamental importance, since it enters the collisional terms of the Boltzmann equation and influences the determination of the scattering probabilities.
An important issue, which we address in this paper, arises when degenerate materials are considered. In this case, there is an ongoing debate in the literature on whether the Pauli principle should be applied to the free flight step or not. This is an important theoretical question, also for computational purposes. An alternative procedure was proposed in [9] where the rejection technique was adopted not only at the end of each scattering event but also at the end of each free flight. This method, although having been cited several times in the literature (see for example [5,6,8,10,11]), to the best of our knowledge it has never been used. In the literature an analysis of this approach is missing (it has received only a brief comment in [5]); it is the aim of this work to fill this gap and to present a coherent discussion on the inclusion of the Pauli principle in a Monte Carlo procedure. We do so by comparing the numerical results for a suspended monolayer graphene obtained with the Monte Carlo Method presented in [9] with those obtained by using the standard EMC in [4] and the new Direct Simulation Monte Carlo (DSMC) in [12], that are by now well established in the semiconductors field and cross-validated with deterministic solutions, for example those based on the discontinuous Galerkin method ( [12], [13]- [18]) or on WENO schemes [19].
The new DSMC scheme in [12] is able to correctly include the Pauli principle; the streaming term of the Boltzmann equation is treated deterministically by means of a splitting procedure, resulting in a rigid translation of the distribution function as a whole. The scattering events are then simulated and the rejection technique is applied at the end of each collision.
The results obtained by using the EMC and the new DSMC are already studied in [12] and allows us to quantify the correctness of the method proposed in [9]. Besides, graphene is a material whose peculiar energy bands make the degeneracy effects relevant, so representing a useful choice as test case.
The paper is organized as follows. In section the semiclassical mathematical model for spatially homogeneous graphene is presented and the simulation procedures are introduced; in sections we present and discuss the results of our simulations; section contains our conclusions.

THE MATHEMATICAL MODEL AND THE MONTE CARLO TECHNIQUES
The carrier population in graphene is made of four components: the electrons of the valence and of the conduction band, which can occupy the states around either of the Dirac points, K and K ′ , of each band. In our simulations, we consider a homogeneous graphene strip of infinite extension in the direction transversal to the applied electric field and only the electrons of the conduction band belonging to the valley around the K point, by considering all valleys as equivalent. We recall that the graphene Brillouin zone B has hexagonal shape, and we shall choose the reference frame of the k-space with the origin to coincide with the K point. The electric field is directed along the x-axes.
With good approximation [20], the dispersion relation for the band energy ε around the equivalent Dirac points is given by where v F is the Fermi velocity and k ℓ is the position of the Dirac point ℓ. We will use Eq.
(1) as dispersion relation because for the electric field strengths usually considered in the applications the charge transport involves almost exclusively the electrons around the K and K ′ points.
Under these conditions and by using the semiclassical approximation, the Boltzmann equation for the charge carriers is where f (t, k) is the electron distribution function of the charge carriers at time t, k = (k x , k y ) is the wave-vector and ∇ k is the gradient with respect to k. The r.h.s. of Eq. 2 is the collision operator which describes the interactions of the carriers with the phonons.
The appropriate initial condition for Eq. 2 in the degenerate case is the Fermi-Dirac where ε F is the Fermi level and T the room temperature, related to the charge distribution by means of where only the spin degeneracy is considered. The Fermi level will be chosen high enough in order to produce a strong degeneracy, according to the situation we want to investigate. This is equivalent to introduce a high n-doping in a traditional semiconductor.
The electron mean energy and velocity are defined as where ε(t, k) and v(t, k) are the particle energy and velocity, respectively, and ρ(t) the time dependent electron density The The general form of the collision term can be written as ( [19,20]) where S(k ′ , k) is the total transition rate and is given by the sum over the types of scatterings described above where the index A runs over the phonon modes. The G (A) (k ′ , k) 2 's are the electron-phonon coupling matrix elements, which describe the interaction mechanism of an electron with an Ath phonon, from the state of wave-vector k ′ to the state of wave-vector k. The symbol δ denotes the Dirac delta function and g A (q) is the phonon distribution for the A-type phonons. In (9), g ± A = g A (q ± ), where q ± = ± (k ′ − k), stemming from the momentum conservation.
For both longitudinal and transversal acoustic phonons, we consider the elastic approximation according to which the transition rate is given by [26] where D ac is the acoustic phonon coupling constant (also called acoustic phonon deformation potential), σ m is the graphene areal density, v A the sound speed of the Ath acoustical phonon branch in graphene and ϑ k ,k ′ is the convex angle between k and k ′ .
The electron-phonon coupling matrix elements of the longitudinal optical (LO), transversal optical (T O) and K phonons are [19] where D O is the optical phonon coupling constant, ω O the optical phonon frequency, D K is the K phonon coupling constant and ω K the K phonon frequency. The angles ϑ k ,k ′ −k and ϑ k ′ ,k ′ −k denote the convex angles between k and k ′ − k and between k ′ and k ′ − k, respectively.
In the Monte Carlo procedure, the scattering rate of each type of scattering plays a fundamental role; the scattering rate Γ A of the Ath type of scattering is defined as For the sake of completeness, we provide here the scattering rates of all types of scatterings used in this paper. For the acoustic phonon scattering we get for the longitudinal and transversal optical phonons where the upper and the lower signs refer to the LO and T O phonons, respectively. H is the Heaviside function and Similarly, the expression of the scattering rate for the K phonon scattering reads These scattering rates are shown in Fig. 1 as functions of the energy.
In this work, we compare the results for charge transport in graphene obtained with three different Monte Carlo procedures: the standard Ensemble MC procedure, a new Ensemble MC procedure and the Free Flight based one.
In general, in a Monte Carlo simulation, the motion of each particle is given by the solution of the semi-classical equation of motion k = −eE, followed by a collisional event after a time ∆t. The total scattering rate is given bỹ and it varies with the energy. A scattering rate Γ ss due to fictitious collision events, called self-scatterings and which do not change the particle state, is also introduced and a new constant total scattering rate is defined as Γ tot =Γ tot + Γ ss . ∆t is then calculated as the ratio (see for example [25]) where η is a random number uniformly distributed in [0, 1]. In general, the value of Γ tot is determined as αΓ max , being Γ max = max (Γ LA + Γ T A + Γ LO + Γ T O + Γ K ), with α > 1 a tuning parameter. Since the previous scattering rates can differ even by two order of magnitude, using the same Γ tot for each time step leads to a very large number of self-scatterings, making the computational cost considerably great. Therefore, in our simulations, we use a variable Γ tot depending on the state of the particle at the current time t: with α = 1.1.
In our simulations, the following three Monte Carlo procedures will be used: After each collision, the new wave-vector k ′ is determined and if the final state is available the initial state k is updated. The availability of the final states has to respect the Pauli principle and this is checked by means of a rejection procedure: a random number ζ, uniformly distributed in [0, 1], is generated and the final state k ′ is available if the condition ζ < 1 − f (k ′ ) holds. This scheme is repeated for each particle.
It consists of two main steps. First, the distribution function is translated as a whole according to the semi-classical equation of motion, and all the particles experience the same free flight; then, for each particle a sequence of collisional events is simulated. The final state after each scattering mechanism is checked by using the rejection technique described above.
In this method, the rejection technique to check the availability of the final states is used not only at the end of each collision but also at the end of each free flight; if the state reached after the free flight governed by the semi-classical equation of motion is not accepted, the particle goes back and nothing happens.

SIMULATION RESULTS
In this section, the results of the three previous procedures to include the Pauli principle in the Monte Carlo method are shown, compared and discussed.
In the simulations, n P = 10 4 (super)-particles are used, the time step is set equal to ∆t = 2.5 fs; the wave-vector space is discretized by means of a uniform square grid , with k xmax = k ymax = 24 nm −1 , and 642 × 642 cells. The physical parameters proposed in [27,28] and reported in Table I are adopted. The NEMC procedure is able to properly take into account the Pauli exclusion principle and the distribution function is correctly reconstructed and bounded between zero and one ( Fig. 2 b), unlike the SEMC method ( Fig. 2 a), with which the distribution takes values larger than one. The section of the distributions along the x direction in Fig. 3 highlights the difference between the results in the two cases. In Fig. 4    In the following subsections, we discuss the main results, in particular for the charge distribution and the mean energy and velocity, when both an initial Fermi-Dirac and a Maxwell-Boltzmann distribution are considered.

Initial Maxwell-Boltzmann distribution
To overcome the difficulty related to the frozen dynamics due to the Fermi-Dirac distribution, Tadyszak et al. [9] proposed to use a high temperature Maxwell-Boltzmann distribution as initial condition in place of the Fermi-Dirac one; the distribution temperature is heuris- tically set equal to 80 T , T being the room temperature, for the silicon case. Along these lines, we introduce the following Maxwell-Boltzmann distribution where the free parameters µ and T * are the electro-chemical potential and the temperature,  In the presence of an applied electric field E, the mean energy and velocity calculated  As it is evident from Fig. 9, also by using an initial Maxwell-Boltzmann distribution, the velocity has the same behavior as in Fig. 4 b), with an initial small rising portion followed the end of each free flight is smaller, as noted in [9] as well. Therefore, the absolute value of the velocity is appreciably higher.
The SEMC and NEMC simulations are not affected by the choice of the initial condition; the mean energy and velocity reach the same stationary values; the only difference is an initial overshoot in the velocity that is not present when the initial Maxwell-Boltzmann distribution is used (see Fig. 10).

The charge distribution of the FFMC
The results obtained by using the FFMC procedure, i.e. almost constant mean energies and negative mean velocities, are certainly unphysical and not in accordance with the solution of the Boltzmann equation, as well as with the SEMC and NEMC methods, and deserve a deeper analysis.
The comparison between the charge distributions obtained in the FFMC simulation when the Fermi-Dirac or the Maxwell-Boltzmann distribution are taken as initial conditions is shown in Fig. 11 for ε F = 0.6 eV and E = 20 kV/cm. The initial Fermi-Dirac distribution undergoes only negligible changes during the time evolution (see Figs. 11 a,b); the Maxwell-Boltzmann as initial distribution allows a coarser collocation of the particles with high mean energy because it occupies a portion of the k-space twice larger than that of the initial Fermi-Dirac (see Fig. 11 c). The initial MB distribution does not seem to evolve towards a realistic particle distribution but, as time goes, it becomes more similar to an irregular Maxwell-Boltzmann c) and charge distribution at 5 ps d). ε F = 0.6 eV and E = 20 kV/cm.

Negative mean velocity with the FFMC method
To understand the origin of the negative mean velocity obtained with the FFMC approach, we investigate the main steps of the Monte Carlo procedure in the three methods, i.e. the free flights and the subsequent scattering events. In the NEMC procedure, the distribution function is translated as a whole at each time step ∆t, and the number of the free flights is equal to that of the simulated particles n P = 10 4 . In the SEMC scheme, the particles are followed one by one; this leads to an incorrect reconstruction of the distribution function, even if the number of free flights is of the same order as in the NEMC case, as shown in Fig. 12 (a) for the case of an initial Fermi-Dirac distribution. In the NEMC approach, the Pauli principle is not imposed at the end of the free flight but only after the scatterings, so that all free flights are allowed and take place. If the Pauli principle is imposed also at the end of the free flights, the number of the accepted free flights is strongly reduced as reported in Fig. 12 (b). Also when initially the MB distribution is used, the The ratio between the number of the accepted free flights with respect to the total is around 0.6 until t = 1.5 ps and drops very much afterwards, making the statistics not significant (see Fig. 13). The free flights, even a few ones, give a positive contribution to the mean velocity. Therefore, the negative values seen with the FFMC in Fig. 9 should originate from the scattering events. The percentage of scatterings whose final states are available is shown in Fig. 14 any scattering, is shown. In Fig.17 a)  value with respect to the FD. This is due to the fact that in the early times of the simulation there are more available states and each particle has a free flight event at least once, so that the previous two samples are almost the same and the corresponding mean velocities are almost equal, as it is evident in Fig. 17 (b).
When a particle experiences a change of state due to either free flight or scattering with the crystal lattice phonons, it keeps moving in a free motion, without changing its velocity during the rest of the simulation; when this velocity is negative, it will remain negative until the simulation has finished. This fact is supported by the count of the percentage of particles that do not change their velocity between two consecutive time steps, as reported in Fig. 18; after about 1.5 ps, almost all the particles perform only a free motion.
In Fig. 19, the contribution of each type of scattering to the mean velocity, obtained by averaging only over the particle population which experienced the corresponding collision events, is shown. The prevalence of the backward scatterings with LO and K phonons is