Numerical analysis of electronic conductivity in graphene with resonant adsorbates: comparison of monolayer and Bernal bilayer
Ahmed Missaoui Jouda Jemaa Khabthani Nejm-Eddine Jaidane Didier Mayou Guy Trambly de Laissardière





Numerical analysis of conductivity in bilayer graphene
A. Missaoui et al.

                     


Laboratoire de Spectroscopie Atomique Moléculaire et Applications, 
Département de Physique, Faculté des Sciences de Tunis, Université Tunis El Manar,  Campus Universitaire 1060 Tunis, Tunisia
Laboratoire de Physique de la matière condensée,
Département de Physique, Faculté des Sciences de Tunis, Université Tunis El Manar,  Campus Universitaire 1060 Tunis, Tunisia
Université Grenoble Alpes, Inst NEEL, 38042 Grenoble, France
CNRS, Inst NEEL, 38042 Grenoble, France
Laboratoire de Physique théorique et Modélisation, CNRS and 
Université de Cergy-Pontoise, 
95302 Cergy-Pontoise, France





 
We describe the electronic conductivity, as a function of the Fermi energy, in the Bernal bilayer graphene (BLG) in presence of a random distribution of vacancies that simulate resonant adsorbates. We compare it to monolayer (MLG) with the same defect concentrations. These transport properties are related  to the values of  fundamental length scales such as the elastic mean free path L_e,  the localization length ξ  and the inelastic mean free path L_i. 
Usually the later, which reflect the effect of inelastic scattering by phonons, strongly depends on temperature T. In  BLG an additional characteristic distance l_1 exists which is the typical traveling distance  between two interlayer hopping events. We find that when the concentration of defects is smaller than 1%–2%, one has l_1 ≤ L_e ≪ξ and the BLG has transport properties that differ from those of the  MLG  independently of  L_i(T). Whereas for larger concentration of defects L_e < l_1 ≪ξ, and depending on L_i(T), the transport in the BLG can be equivalent (or not)  to that of two decoupled MLG. 

We compare two tight-binding model Hamiltonians with and without hopping beyond the nearest neighbors.






72.15.Lh 
72.15.Rn 
73.20.Hb 
72.80.Vp 



73.23.-b  

      
 


Numerical analysis of electronic conductivity in graphene with resonant adsorbates: 
comparison of monolayer and Bernal bilayer 

    Ahmed Missaoui1,
 Jouda Jemaa Khabthani2,
  Nejm-Eddine Jaidane1,
  Didier Mayou3,4 Guy Trambly de Laissardière5 



    January 27, 2017
==================================================================================================================================




§ INTRODUCTION


Graphene consists of a monolayer (MLG) carbon atoms, with sp^2 hybridization, forming a 2D honeycomb lattice, with two equivalent atoms –atom A and atom B– in a unit cell. 

Linear dispersion relation of the p_z electron states close to the Fermi energy induces many fascinating transport properties which give rise to potential device applications <cit.>. 
Few-layer graphene also present unusual properties. In particular the Bernal bilayer graphene (BLG) with AB stacking, as in graphite, 
breaks the atom A / atom B symmetry and leads to quadratic dispersion relation <cit.>.

Electronic transport is sensitive to static defects which are for example screened charged impurities, or local defects like vacancies or adsorbates, 

(hydrogen, adatoms or admolecules,
chemically bound to one carbon atom of the surface of graphene layer).
 
Theoretical studies of the effects introduced by the adsorbates on the conductivity  has been done for MLG 
(Refs. <cit.> and Refs. therein),
and  for BLG <cit.>.

Most of them consider a standard Hamiltonian that takes only into account the  hopping between the nearest neighbors orbitals.

Yet some studies show the importance of hopping beyond nearest neighbors on electronic structure and transport properties <cit.>.

In this paper, we present electronic properties of MLG and BLG obtained by two tight-binding (TB) models: 
the standard model with nearest neighbor only (TB1) and a TB model including the effect of the hopping beyond nearest neighbors (TB2).
This second model, which is more realistic, predicts some differences in the transport at energies close to the resonant energy of scatters. 


We consider local defects, such as adsorbates or vacancies, that are resonant scatters.  
Local defects tend to scatter electrons in an isotropic way for each valley and lead also to strong intervalley scattering. The T matrix of a local defect usually depends strongly on the energy. In the case of simple vacancies or adsorbates (atoms or molecules) that create one covalent bound with a carbon atoms of MLG (BLG), the T matrix diverges at the energy E_MG 
(with TB1 model E_MG=E_D=0).
For this reason, theses scatters are called resonant scatters.

The adsorbate is simulated by a simple vacancy in the plane of p_z orbital 
as usually done <cit.>. 
Indeed the covalent bonding between the adsorbate and the carbon atom of graphene to which it is linked, 
eliminates the p_z orbital from the relevant energy window. The scatterers are distributed randomly in both planes and with the same concentration in both planes. 

We consider  here that the up and down spins are degenerate i.e. we deal with  a paramagnetic state. Indeed the existence and the effect of a magnetic state for various adsorbates or vacancies is still debated <cit.>.  Let us emphasize that  in the case of a magnetic state  the up and down spins  give two different contributions to the conductivity but the individual contribution of each spin can be analyzed from the results discussed here. 



We first determined the density of states (DOS) in disordered MLG and BLG in presence of static scatterers (vacancies) with various concentrations c from 0.5% to 10%. 

Elastic mean free path L_e, which depends on the distribution of scatters and on energy E, is also computed.
From diffusive properties of wave packet in the structure, the electrical conductivity σ 

is computed versus E and the inelastic mean free path L_i. L_i due to electron-phonons interaction or magnetic fied.
Roughly speaking, large L_i correspond to low temperature limit and small L_i to room temperature.
The numerical method used takes into account all quantum effects.
We show that difference between MLG and BLG is explained by considering the average distance l_1 over which a charge carrier travels in  a layer between two interlayer hoppings <cit.>. 
As explained in this paper, l_1  ≃ 2nm and it is almost independent on the concentration of defects. 
Indeed BLG has some similar properties to MLG for small L_i and large c, i.e. when L_e ≃ L_i < l_1. 
But for small c values, L_e ≃ L_i ≥ l_1,  the effects of interlayer hopping affect the electronic properties of BLG with respect to MLG case. 
In this case, the conductivity of the BLG varies with c like in usual metals where DOS is finite whereas MLG behaves like a semi-metals with Dirac electrons.
Moreover, different regimes of transport in BLG are found depending on the values of the energy, 
like in MLG <cit.>.



Finally the localization length ξ is computed and we study the localization regime (which can be observed experimentally in a very low temperature regime, i.e. large L_i). 

For the concentrations studied here, l_1 < ξ, which means that in this localization regime the coupling between the two layers plays always a significant role. Thus behavior of the BLG and MLG are different and  localization length is larger in the BLG.


In Sect. <ref>, the two tight binding (TB) models used are described and the corresponding density of states are discussed. 
Transport properties are presented in Sect. <ref>. We first describe rapidly the computational method and the relevant lengths to analyze conductivity in BLG (Sect. <ref>). Then elastic scattering length L_e (Sect. <ref>) and microscopic conductivity σ_M (Sect. <ref>) are presented.
Finally quantum interference corrections to the conductivity and localization length ξ  are analyzed (Sect. <ref>).  
Sect. 4 is devoted to concluding remarks. 



§ ELECTRONIC STRUCTURE




 §.§ Tight binding Hamiltonian models



BLG can be considered as two coupled MLGs with the top layer shifted by a carbon bond from the bottom layer. Consequently, BLG consists of four atoms in its unit cell, two carbons A_1, B_1 from the unit cell in the bottom layer and A_2, B_2 in the the top layer where A_2 sits at the top of A_1 (Fig. <ref>). 
We used a tight binding scheme (TB) <cit.>.
Only p_z orbitals are taken into account since we studied the electronic properties around the Fermi energy  level. Interlayer interactions are not restricted to the ppσ terms but ppπ terms have also to be introduced. The Hamiltonian has the form:

    Ĥ =∑_iϵ_i|i⟩⟨ i| + ∑_(i,j)t_ij|i⟩⟨ j| ,

where i is the orbital located at r_i and 
the sum runs over all neighboring i, j sites. 
The energy on the site is  taken equal to zero for first nearest neighbor model. t_ij is the hopping element matrix between site i and site j, computed from the Slater-Koster parameters,

    t_ij
     =    n_c^2 V_ppσ(r_ij)  +  (1 - n_c^2) V_ppπ(r_ij),

where V_ppσ and V_ppπ depend from the distance r_ij and n_c is the cosines direction along the Oz axis. It is either equal to zero or to a constant because the two graphene layers have been kept flat in our model. 
The Slater Koster parameters are exponentially decaying function of the distance: 

    V_ppπ(r_ij)       =   -γ_0 exp( q_π(1-r_ij/a) ) , 
    
    V_ppσ(r_ij)    =   γ_1 exp( q_σ(1-r_ij/a_1) ) .

It allows, according to the value of q_π, to take into account both first neighbors or second neighbors. 
a is the nearest neighbor distance within a layer, a=1.418 Å, and  a_1 is the interlayer distance, a_1=3.349 Å. First neighbors interaction in a plane is characterized by the commonly used value γ_0=2.7 eV and the second neighbors interaction γ_0^' is set to 0.1γ_0. 
The ratio q_π/a in Eq. <ref> is fixed by the value of the γ_0^'. The interlayer coupling between two p_z orbitals in 
π configuration is γ_1,  γ_1=0.48 eV, and it is fixed to obtain a good fit with  ab initio calculation around Dirac energy in AA stacking and AB bernal stacking <cit.>. 
It is worth to note that we choose the same decay coefficient  for V_ppπ and V_ppσ:

    q_σ/a_1 = q_π/a = log(γ_0/γ_0' )/a' - a =  2.218Å^-1,

with a'=2.456 Å the distance between second neighbors in a plane. All p_z orbitals have the same on-site energy ϵ_i in two planes. 
In order to obtain a Dirac energy E_D equal to zero in MLG, one fixes ϵ_i equal to -0.78 eV for TB model with hopping beyond first neighbors. 
This is necessary because hopping beyond first neighbors breaks the electron/hole symmetry and then shifts E_ D value 
<cit.>.
We have used this model Hamiltonian in our previous work <cit.> to study the electronic structure rotated graphene bilayer. 
In order to analyze the effect of hopping beyond the first-neighbor distances, we consider also the simplest TB model (TB1) with first-neighbor hopping only and the same parameters than the complete TB2 model described above.

As explained in the introduction we consider that resonant adsorbates –simple atoms or molecules such as H, OH, CH_3–  create a covalent bond with some carbon atoms of the BLG. To simulate this covalent bond, we assume that the p_z orbital of the carbon, that is just below the adsorbate, is removed. In our calculations the mono-vacancies are distributed at random between the both planes with a finite concentration c. 




 §.§ Density of states







Figure <ref>(b) shows the total density of states (total DOS) n(E) in BLG for different concentrations c = 0.5% to 10% of defects randomly distributed.

For comparison, the DOSs of MLG with same TB Hamiltonian 

<cit.> 

are  also shown (Fig. <ref>(a)).
Without defects, c=0, BLG DOS and MLG DOS differ for energies E such as 
-γ_1 < E < γ_1 <cit.>. For small c concentrations, c< ∼ 1%, this distinction is still observed. 
But remarkably, MLG and BLG have very similar DOS, for concentration of defects c larger than ∼ 1%. 

With TB1 model (first-neighbor hopping only), states occur at energy E_ MG=0. This is reminiscent of the midgap state produced by asymmetry between the number of atoms A and  B in monolayer graphene <cit.>. 

In agreement with previous findings for MLG <cit.> and BLG <cit.>,

for large values of energies  E the DOS is weakly affected by the presence of disorder. Finally near the Dirac energy there is an intermediate regime where the pseudo-gap is filled (“plateau”). 

With the TB2 model including all neighbors, the midgap state is no longer at E=0, but it becomes a broad peak at negative energy E_ MG. 
E_ MG value varies from E_  MG≃-0.2 eV to E_ MG≃-0.3 eV when defect concentration increases from c=0.5% to 10%, like in monolayer graphene (Fig. <ref>(a)) 

<cit.>. 

In the following we distinguish  three cases according to energy values:


  (i) Sufficiently large energies, BLG DOS and MLG DOS are similar and they are not strongly modified by the presence of resonant defects. 

  (ii) Energies in the “plateau” due to vacancies but not in the midgap states E ≠ E_ MG.

  (iii) Energies in the midgap states, i.e. E=E_ MG=0 for TB1 (with first-neighbor hoppings only) and E≃ E_ MG for TB2 (with hopping beyond first neighbors). 

These three cases correspond to different transport 
regimes.





§ TRANSPORT PROPERTIES





 §.§ Computational method and relevant lengths



In the framework of the Kubo-Greenwood formula for electronic transport properties, 
the quantum diffusion coefficient D (diffusivity) and conductivity σ are computed by using the polynomial expansion method, developed by Mayou, Khanna, Roche and Triozon <cit.>.  
This numerical approach
allows very efficient calculations by recursion (Lanczos algorithm) in real-space which take into account all quantum effects.
It has been used to study quantum transport in disordered graphene, chemically doped graphene and bilayer (see Refs. <cit.> and Refs. therein). 
Our calculations are performed on sample containing up to a few 10^7 carbon atoms, which corresponds to typical sizes of about one micrometer square and allows to study systems with elastic mean-free length of the order of few hundred nanometers.


Elastic scattering events are taken into account in the Hamiltonian, but effects of inelastic scattering by phonons at temperature T are not included in the Hamiltonian. To consider the inelastic processes, we introduce an inelastic scattering time τ_i(T) beyond which the propagation becomes diffusive due to the destruction of coherence by these inelastic processes. The effect of a magnetic field on the electron propagation is not included  directly in the TB model, but a magnetic field B can have also a similar 

incoherent

dephasing effect.
This dephasing effect  occurs on a length L_i(B) such that the flux of the magnetic field enclosed in the disk of radius L_i(B) is equal to the flux quantum h/e, i.e. L_i(B)≃√(h/eB). 
We treat these two dephasing effects in a phenomenological way through a Relaxation Time Approximation (RTA) as described here after.
In the RTA, the conductivity  along the x-axis is given by, <cit.>

    σ(E_F,τ_i)    =    e^2n(E_F)D(E_F,τ_i) , 
    
    D(E_F,τ_i)    =   L_i^2(E_F,τ_i)/2τ_i ,

where E_F is the Fermi energy, n(E_F) is the density of states (DOS)  and L_i is the inelastic mean-free path conductivity  along the x-axis. L_i(E_F,τ_i) is the typical distance of propagation during the time interval τ_i for electrons at energy E.

We compute the distance L_i, the diffusivity D and the conductivity σ at all inelastic scattering times τ_i and all energies E for model Hamiltonian that includes inelastic scatters distributed randomly in the super-cell.
At short times τ_i –i.e. τ_i lower than elastic scattering time τ_e– the propagation is ballistic and the conductivity σ increases when τ_i increases (Fig. <ref>), 

    L_i(E,τ_i) ≃ V_0(E) τ_i   when τ_i ≪τ_e,

where V_0(E) is a velocity  at the energy E and short time t.
In crystals,  V_0 ≥ V_B where V_B is the Boltzmann velocity (intra-band velocity) <cit.>.
In BLG and MLG, V_0 and V_B have the same order of magnitude: 
V_0( BLG) = V_0( MLG) = √(2)V_B( MLG) <cit.>.
According to the renormalization theory <cit.> in 2D systems with static defects, diffusivity D always goes to zero at very large τ_i. 
At each energy, the microscopic diffusivity D_M (microscopic conductivity σ_M) is defined as the maximum value of D(τ_i) (σ(τ_i)).
We compute also the elastic mean-free path L_e along the x-axis, from the relation <cit.>,

    L_e(E) = 1/V_0(E) Max_τ_i{L_i^2(E,τ_i)/τ_i} = 2 D_M(E)/V_0(E).
 
L_e is the average distance between two elastic scattering events. At each energy, the elastic scattering times τ_e is deduced from L_e by L_e(E) = V_0(E) τ_e(E). 


In our calculations τ_i and L_i  are considered as adjustable parameters.
Roughly speaking, when L_i ≪ L_e (τ_i ≪τ_e) the inelastic disorder dominates; it should correspond to very high temperatures. 
When L_i ≃ L_e (τ_i ≃τ_e), the conductivity is equal to microscopic conductivity, which should correspond to high temperature cases, typically room temperature. 
And when L_i ≫ L_e (τ_i ≫τ_e), quantum localization will dominates transport properties; this is the  low temperature limit. 
In the following we therefore discuss the two important cases: L_i ≃ L_e (τ_i ≃τ_e) and L_i ≪ L_e (τ_i ≪τ_e).




Figure <ref> shows the variation of the conductivity σ and the inelastic mean free path L_i versus τ_i for energies corresponding to the three previous cases (see Sect. <ref>).
The first case (i) (E=1.5 eV in Fig. <ref>)

corresponding to a  Boltzmann transport: for large values of τ_i, the conductivity σ is almost constant as expected in a diffusive regime. This regime corresponds to energies for which the DOS is weakly affected by scatters.
The third case (iii) (E=E_ MG in Fig. <ref>) is determined by the transport of the midgap states which are localized states.  
The latter case (ii) is an intermediary regime between the previous two: 
for τ_i closed to the elastic scattering time τ_e, there is a diffusive behavior where the σ(τ_i) reaches a maximum, σ_M; for larger values of τ_i, τ_i ≫τ_e, σ(τ_i) decreases progressively as expected in localization regime due to Anderson localization in 2D <cit.>.

In BLG another relevant time is the average traveling time t_1 between two interlayer hoppings of the charge carriers, which is associated to an
average traveling distance l_1 in a layer between two interlayer hoppings  <cit.>.
In perfect BLG typical values of t_1 and l_1 can be easily estimated : t_1 = ħ / Γ_1 ≃ 2 × 10^-15 s, where Γ_1 ≃ 0.4eV is the 
interlayer hopping parameters of the Hamiltonian, and  

l_1 ≃ V_m t_1 ≃ 2  nm 

<cit.>

where V_m is the velocity in MLG, 
V_m ≃  V_0 ≃ 10^6 ms^-1. 
When there is elastic disorder such that τ_e < t_1 the value of t_1 can be modified. A simple argument may be given as follows: A Bloch state  of the  MLG is still coupled to Bloch states of the other layer by the same intensity, typically Γ_1, but these states are no longer eigenstates and have a typical lifetime τ_e. Because of that they have a spectral width W≃ħ /τ_e. From the Fermi Golden rule the typical time needed to jump from one layer to the other will be such that  ħ/t'_1 ≃Γ_1^2/W. 
Therefore the new value of the interlayer hopping time t'_1 will be larger and will be such that  t'_1≃ t_1 (ħ/Γ_1 τ_e). Since the propagation is diffusive on the timescale t>τ_e with diffusion coefficient D≃ V_0^2 τ_e, the new length l'_1 in presence of defects is obtained from the relation, l'_1^2/t_1' ≃ V_0^2 τ_e,
and thus l_1' ≃ l_1 ≃ 2 nm depends weakly on disorder (i.e. l_1' almost independent  on τ_e). 

As shown in the following, electronic properties of disordered BLG depend on the values of the length L_e, ξ and l_1 which are characteristic of the BLG and of the amount of elastic scatters:


  * For low concentration of defects, c < c_l = 1%-2%, and for E  E_MG, we have, ł_1 ≤ L_e ≪ξ, and thus electronic  properties of BLG are influenced by  interlayer hoppings for every L_i ≥ L_e values. 

  * For larger concentration of defects, c > c_l, one obtains, Ł_e < l_1 ≪ξ.  
When the effect of quantum interferences on conductivity is small, i.e.  when L_e ≃ L_i < l_1, BLG  behaves as two decoupled MLG. 
For Ł_e < l_1 ≤ L_i ≪ξ, the coupling between the two planes can affect the propagation of charge carriers before inelastic scattering makes the propagation diffusive (i.e. on the length scale L_i). In this regime, quantum corrections to transport are not the same in the BLG and in two decoupled MLG. This influences strongly the localization regime as we discuss in Sect. <ref>. 












 §.§ Elastic mean-free path




The elastic mean-free path L_e (Eq. (<ref>)) along the x-axis 
as a function of the E is shown figure <ref> for different values of vacancy concentrations c in MLG and BLG with both TB models.
It depends on the energy even in the intermediate regime and it takes a finite and non-zero value for E=E_D but stays comparable to the distance d between adsorbates (vacancies) defined by, d ≃1/√(n_a), 
where n_a is the adsorbates density.
Numerical results (Fig. <ref>) show that L_e values in BLG and MLG 

are close to each other.

Moreover,
L_e(E) < l_1 ≃ 2 nm (i.e. τ_e < t_1)
for c > c_l ≃ 1%–2%; 
whereas  for smaller c, L_e(E) ≥ l_1 (i.e. τ_e ≥ t_1)
for c < c_l.






 §.§ Microscopic conductivity








As shown Fig. <ref>, the static scattering events perturb strongly the wave packet propagation and a maximum value of the conductivity σ(τ_i), called microscopic conductivity σ_M, σ_M(E) = e^2n(E)D_M(E), is reached. 
σ_M calculated from both TB models are shown in Figs. <ref> and <ref> for different  concentrations c of vacancies  in MLG and BLG.
According to the renormalization theory <cit.>, this value is obtained when the inelastic mean free path L_i and the elastic mean free path L_e are comparable, L_e ≃ L_i, which corresponds to τ_i ≃τ_e. 
As L_i and τ_i decrease when the temperature T increases, the microscopic conductivity is a good estimation of the high temperature conductivity (or room temperature conductivity).  

For energies corresponding the to Boltzmann regime, i.e. regime (i) described in Sec. <ref>, σ_M≃σ_B, where σ_B is calculated with the Bloch-Boltzmann approach <cit.>. In this regime the conductivity decreases with the concentration of defects. 

In the intermediate energy values regime (ii), the semi-classical approach fails and the behavior depends on c. 
From 

Fig. <ref>(a),

for small c values, typically 
c< c_l ≃ 1%–2% (i.e. L_i ≃ L_e ≥ l_1), 
σ_M seems to reach a constant minimum value (“plateau”), but this minimum  σ_M value decreases as c  increases. 

This concentration dependence is specific to BLG and is not observed in MLG.  

For larger c 

(Fig. <ref>(b)), 

c>c_l (i.e. L_i ≃ L_e < l_1),
σ_M reaches a minimum value independent on c value: 
σ_M≃ 1.2 G_0 where G_0 = 2 e^2/h. 
This values for BLG is <cit.> two times the universal value of the conductivity, ∼ 4e^2/(π h),
expected 

in presence of resonant scatters in MLG <cit.>.


Results with TB2 model (including hopping beyond nearest neighbors) 

(fig. <ref>)

show that a plateau of the microscopic  conductivity near the Dirac energy exists 

in MLG and BLG,

but is not symmetric due to the symmetry breaking electron-hole. 
Nevertheless in this case, 

for c>c_l, 

σ_M values are still close to the universal conductivity plateau;

whereas for c<c_l, the minimum value of σ_M decreases as c  increases like with TB1 model.


For energies in the midgap states (iii) with TB1 model (first neighbor hopping only),
an anomalous behavior of the conductivity is obtained and there is a peak of σ_M at E=E_ MG=0. 
With TB2 (including hopping beyond nearest neighbors), this anomalous behavior is still slightly present at 
E ≃ E_  MG, but the change in the conductivity is rather small. 
Thus, as in monolayer graphene <cit.>, conduction by the midgap states is very specific to TB1 model. 




 §.§ Quantum localization regime







In this section we consider the case of large inelastic mean-free path, 
L_i ≫ L_e (i.e. τ_i ≫τ_e), 
which should corresponds to the  low temperature limit.
The conductivity in the intermediate energy case (ii) (see Sect. <ref>) presents localization effects that  are a consequence of quantum interferences. 
In that case, for all vacancy concentration values c, L_i ≫ l_1 (i.e. τ_i ≫ t_1). 
Therefore BLG should have different localization behavior than MLG.
As shown in Fig. <ref>, σ(L_i) follows the linear variation with the logarithm of the inelastic mean free path L_i, 
like in the case of monolayer graphene <cit.>,

    σ(E,L_i) = σ_0(E) - α G_0 log(L_i/L_e(E)),

where G_0 = 2 e^2/h, L_e(E) is defined by Eq. (<ref>), and σ_0 values are in the range of σ_M values. 
For low concentration of defects c=0.5% and 1% one can estimate α≃ 0.26, and for larger c, c = 3% and 5%,
α≃ 0.32.
These values are close to the result found in monolayer graphene with same computational method <cit.> and close to the prediction of the perturbation theory of 2D Anderson localization for which α = 1/π <cit.>. 
As in the monolayer case, this linear variation of σ with log L_i is found for both models, with nearest neighbor hopping only and with hopping beyond first nearest neighbors.
 
We finally define the localization ξ length from the expression (<ref>) by extrapolation of σ(L_i) curves  (Fig. <ref>) when : σ(L_i = ξ)=0, and then 

    ξ(E) = L_e(E) exp( σ_0(E)/α G_0).

The ξ values for energies in the plateau of σ_M (i.e. case (ii) described in Sect. <ref>)  are shown figure <ref>. 

For large concentrations of defects (c>c_l), ξ in MLG and BLG is almost independent on the energy. 
Moreover ξ(BLG) is always larger than ξ(MLG). 

This difference results from the fact that σ_0, BLG≃ 2 σ_0, MLG (Sect. <ref>) and then, for c>c_l,

    ( ξ(E)/L_e(E) )_BLG≃ ( ξ(E)/L_e(E) )_MLG^2 ,

where L_e are similar in MLG and BLG (Sec. <ref>). 
It is thus a multilayer effect on quantum interferences that modifies the 2D behavior with respect to MLG cases.  
For low defects concentration (c<c_l), interlayer coupling modifies also quantum interferences.

Therefore,  for every resonant scatterer concentration,  quantum interference corrections to the conductivity in BLG and in MLG  are not similar. 




§ CONCLUSION

To conclude we have studied numerically the quantum diffusion of charge carriers in monolayer graphene and Bernal bilayer graphene in the presence of local defects. These defects are simulated by simple vacancies  randomly distributed in the structure. Among the fundamental length scales in the MLG and BLG  there are the elastic mean free path  L_e, the localization length ξ and the  inelastic mean-free path L_i which in real systems depends on temperature or on magnetic field. For the bilayer, there is an additional length scale which is the  typical distance l_1≃ 2 nm over which an electron travels in a plane before hopping to the other plane. We have compared the Bilayer and monolayer transport properties for identical concentrations of vacancies. We show that these properties can be either similar or different depending on the  comparison between l_1 and the other three length scales L_e, L_i, ξ. This relation explains essentially the numerical results detailed in this paper.

Our results show that for strong concentration of defects, c > c_l ≃ 1%–2%, the bilayer graphene could be equivalent to two independent disordered monolayers of graphene, because the elastic mean free path L_e is smaller than the average distance l_1. 
Therefore for c > c_l, the universal aspects of the conductivity  are present in bilayer, as in monolayer graphene, with (TB1) or without (TB2) the hopping beyond nearest neighbors.
 
In the high temperature limit, i.e when inelastic scattering length L_i is small, L_i ≃ L_e, the conductivity in bilayer  

is almost equal to two times the universal minimum plateau
of microscopic conductivity in monolayer graphene 

(except for the Dirac energy with TB1 model that takes only into account nearest neighbor hopping). 

For smaller c, c < c_l, the BLG should be considered like a usual metal:  with static defects the minimum of microscopic conductivity of BLG increases when the c values decreases.

For the parameters studied here, the localization length ξ is larger than the traveling distance l_1 between two interlayer hopping; and therefore the BLG and MLG have different localization lengths at the same concentration, even if they have similar elastic-mean-free paths. The localization length is the largest in the BLG. In the limit, i.e L_i ≫ L_e, (which is relevant at low temperature) and for all c values,
the conductivity follows a linear variation with the logarithm of L_i in MLG and BLG and for both TB1 models (nearest neighbor hopping only) and TB2 (with hopping beyond nearest neighbors), excepted for energy in the midgap states for TB1. This is in good agreement with two-dimensional Anderson localization and consistent with the expected universal behavior of conductivity of a two dimensional disordered system <cit.>. 




§ ACKNOWLEDGMENT

The authors wish to thank C. Berger, W. A. de Heer, L. Magaud, P. Mallet and J.-Y. Veuillen
for fruitful discussions.
The numerical calculations have been performed 
at  Institut Néel, Grenoble,
and at the Centre de Calculs (CDC),
Université de Cergy-Pontoise.

This work was supported by the Tunisian French Cooperation Project (Grant No. CMCU 15G1306)
and the project ANR-15-CE24-0017.





§ AUTHOR CONTRIBUTION STATEMENT

Numerical calculations have been performed by A. Missaoui and G. Trambly de Laissardière. 
All authors took part in drafting of the manuscript as well as in analysis and interpretation of the
results.


















Wallace47 

P. R. Wallace,  Phys. Rev. 71, 622 (1947)

Berger06
C. Berger, Z. Song, X. Li, X. Wu, N. Brown, C. Naud, D. Mayou, T. Li, J. Hass, A. N. Marchenkov, E. H. Conrad, P. N. First, W. A. de Heer, Science 312, 1191 (2006)

Katsnelson06
M. I. Katsnelson, K. S. Novoselov, A. K. Geim, Nature Phys 2, 9 (2006)

Castro09_RevModPhys
A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, A. K. Geim, 
Rev. Mod. Phys 81, 109 (2009) 

Latil06
S. Latil, L. Henrard, Phys. Rev. Lett 97, 036803 (2006)

Castro07
E. V. Castro, K. S. Novoselov, S. V. Morozov, N. M. R. Peres, J. M. B. Lopes dos Santos, J. Nilsson,
F. Guinea, A. K. Geim, A. H. Castro Neto,

Phys. Rev. Lett 99, 216802 (2007)

MacdoABC
F. Zhang, B. Sahu, H. Min, A. H. MacDonald, Phys. Rev. B 82, 035409  (2010) 

Varchon_prb08
F. Varchon, P. Mallet, J.-Y. Veuillen, L. Magaud, Phys. Rev. B 77, 235412 (2008)

Ohta06
T. Ohta, A. Bostwick, T. Seyller, K. Horn, E. Rotenberg, Science 313, 951 (2006)

Brihuega08
I. Brihuega, P. Mallet, C. Bena, S. Bose, C. Michaelis, L. Vitali, F. Varchon, L. Magaud, K. Kern, J. Y. Veuillen, Phys. Rev. Lett. 101, 206802 (2008)

DasSarma11

S. Das Sarma, S. Adam, E. H. Hwang, E. Rossi
Rev. Mod. Phys 83, 407 (2011) 


McCann13

E. McCann and M. Koshino,
Rep. Prog. Phys. 76 056503 (2013)


Ulstrup14
S. Ulstrup, J. C. Johannsen, F. Cilento, J. A. Miwa, A. Crepaldi, M. Zacchigna,
C. Cacho, R. Chapman, E. Springate, S. Mammadov, F. Fromm, C. Raidel, T. Seyller, F. Parmigiani, M. Grioni, P. D. C. King,  P. Hofmann1,

Phys. Rev. Lett. 112, 257401 (2014)

V.M.Pereira
V. M. Pereira, J. M. B. L. dos Santos, A. H. C. Neto, Phys. Rev. B 77, (2008) 115109


Robinson08
J. P. Robinson, H. Schomerus, L. Oroszlány, and V. I. Fal’ko,

Phys. Rev. Lett. 101, 196803 (2008)
Yuan10b
S. Yuan, H. De Raedt, M. I. Katsnelson,

Phys. Rev. B 82, 115448 (2010)
Yuan12
S. Yuan, T. O. Wehling, A. I. Lichtenstein, M. I. Katsnelson,

Phys. Rev. Lett. 109, 156601 (2012)


Lherbier12
A. Lherbier, S. M.-M. Dubois, X. Declerck, Y.-M. Niquet, S. Roche, J.-C. Charlier, Phys. Rev. B 86, 075402 (2012)

Roche12
S. Roche, N. Leconte, F. Ortmann, A. Lherbier, D. Soriano, and J.-C. Charlier, Solid States Comm. 152, 1404 (2012)

Roche13
A. Cresti,  F. Ortmann, T. Louvet, D. Van Tuan, S. Roche, Phys. Rev. Lett. 110, 196601 (2013) 


Kretinin13
A. Kretinin, G. L. Yu, R. Jalil, Y. Cao, F. Withers, A. Mishchenko, M. I. Katsnelson, K. S. Novoselov, A. K. Geim, F. Guinea,

Phys. Rev. B 88, 165427 (2013)


Trambly11
G. Trambly de Laissardière, D. Mayou, Modern Physics Lett. B 25, 1019 (2011)

Trambly13

G. Trambly de Laissardière, D. Mayou, Phys. Rev. Lett. 111, 146601 (2013)  

Trambly14
G. Trambly de Laissardière, D. Mayou, Adv. Nat. Sci.: Nanosci. Nanotechnol. 5, 015007  (2014)


Zhao15
P.-L. Zhao, S. Yuan, M. I. Katsnelson, H. De Raedt, 

Phys. Rev. B 92, 045437 (2015)


Koshino06
M. Koshino, T. Ando,

Phys. Rev. B 73, 245403 (2006)

Adam08
S.  Adam, S. Das Sarma,

Phys. Rev. B 77, 115436 (2008)

Koshino09
M. Koshino,

New J. Phys. 11,  095010 (2009)


Yuan10
S. Yuan, H. De Raedt,  M. I. Katsnelson,

Phys. Rev. B 82,  235409 (2010)

Gonzalez10
J. W. González, H. Santos, M. Pacheco, L. Chico, L. Brey,

Phys. Rev. B 81, 195406 (2010)

Ferreira11
A. Ferreira, J. Viana-Gomes, J. Nilsson, E. R. Mucciolo, N. M. R. Peres, A. H. Castro Neto,

Phys. Rev. B 83, 165402 (2011)

VanTuan16
D. Van Tuan, S. Roche, 
Phys. Rev. B 93, 041403 (2016) 


Nair12

R. R. Nair, M. Sepioni, I.-L. Tsai, O. Lehtinen, J. Keinonen, 
A. V. Krasheninnikov, T. Thomson, A. K. Geim, and I. V. Grigorieva,

Nature Phys. 8, 199 (2012).
Scopel16
W. L. Scopel, W. S. Paz, J. C. C. Freitas,

Sol. Stat. Comm. 240, 5 
(2016)

Snyman07
I. Snyman, C. W. J. Beenakker,

Phys. Rev. B 75, 045322 (2007)



Trambly12
G. Trambly de Laissardière, D. Mayou, L. Magaud, Phys. Rev. B 86, 125413  (2012) 

Trambly10
G. Trambly de Laissardière, D. Mayou, L. Magaud, Nano Lett. 10, 804 (2010)

Trambly16
G. Trambly de Laissardière, O. F. Namarvar, D. Mayou, L. Magaud, Phys. Rev. B 93, 235135 (2016)

Trambly14CRAS
G. Trambly de Laissardière, D.  Mayou,
C. R. Physique 15, 70 (2014)

Lee
P. A. Lee, T. V. Ramakrishnan, Rev. Mod. Phys. 57, 287 (1985)

P.W.Brouwer 
P. W. Brouwer, E. Racine, A. Furusaki, Y. Hatsugai, Y. Morita, and C. Mudry, Phys. Rev. B 66, 014204 (2002) 

Mayou88
D. Mayou, Europhysics Lett. 6, 549 (1988)

Mayou95
D. Mayou and S.N. Khanna, J. Phys. I France 5, 1199 (1995) 

Roche97
S. Roche and D. Mayou, Phys. Rev. Lett. 79, 2518 (1997) 

Roche99
S. Roche and D. Mayou

Phys. Rev. B 60, 322 (1999)

Triozon02
F. Triozon,  J. Vidal, R. Mosseri, D. Mayou, Phys. Rev. B 65, 220202 (2002) 



