Electronic transport properties and quantum
localization effects monitored by selective
functionalization in Bernal bilayer graphene

Jouda Jemaa Khabthani, Ahmed Missaoui, Didier Mayou, Guy Trambly de Laissardière

► To cite this version:

Jouda Jemaa Khabthani, Ahmed Missaoui, Didier Mayou, Guy Trambly de Laissardière. Electronic transport properties and quantum localization effects monitored by selective functionalization in Bernal bilayer graphene. Physical Review B, 2021, 104 (24), pp.245125. 10.1103/Phys-RevB.104.245125. hal-03488199

HAL Id: hal-03488199

https://hal.science/hal-03488199

Submitted on 25 Aug 2023

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Electronic transport properties and quantum localization effects monitored by selective functionalization in Bernal bilayer graphene

Jouda Jemaa Khabthani,<sup>1,*</sup> Ahmed Missaoui<sup>1,2,†</sup> Didier Mayou<sup>1,3,4,‡</sup> and Guy Trambly de Laissardière<sup>1,2,§</sup>

<sup>1</sup>Laboratoire de Physique de la Matière Condensée, Département de Physique, Faculté des Sciences de Tunis,
Université de Tunis El Manar, Campus universitaire 1060 Tunis, Tunisia

<sup>2</sup>Laboratoire de Physique Théorique et Modélisation, CY Cergy Paris Université, CNRS, 95302 Cergy-Pontoise, France

<sup>3</sup>Université Grenoble Alpes, Inst. NEEL, F-38042 Grenoble, France

<sup>4</sup>CNRS, Inst. NEEL, F-38042 Grenoble, France

<figure><img src="image_2.png" /></figure>

(Received 6 May 2021; revised 24 November 2021; accepted 29 November 2021; published 15 December 2021)

Monitoring electronic properties of two-dimensional (2D) materials is an essential step to open a way for applications such as electronic devices and sensors. From this perspective, Bernal bilayer graphene (BLG) is a fairly simple system that offers great possibilities for tuning electronic gap and charge carriers' mobility by selective functionalization (adsorptions of atoms or molecules). Here we present a detailed numerical study of BLG electronic properties when two types of adsorption site are present simultaneously. We focus on realistic cases that could be realized experimentally with adsorbate concentration *c* varying from 0.25% to 5%. For a given value of *c*, when the electronic doping is lower than *c*, we show that quantum effects, which are ignored in usual semiclassical calculations, strongly affect the electronic structure and the transport properties. A wide range of behaviors is indeed found, such as gap opening, metallic behavior, or abnormal conductivity, which depend on the adsorbate positions, the *c* value, the doping, and eventually the coupling between midgap states which can create a midgap band. These behaviors are understood by simple arguments based on the fact that BLG lattice is bipartite. We also analyze the conductivity at low temperature, where multiple scattering effects cannot be ignored. Moreover, when the Fermi energy lies in the band of midgap states, the average velocity of charge carriers cancels but conduction is still possible thanks to quantum fluctuations of the velocity.

DOI: 10.1103/PhysRevB.104.245125

I. INTRODUCTION

Monolayer graphene (MLG) is a two-dimensional carbon layer that has been of increasing interest for the scientific community since its first experimental realization in 2004 [1–3]. Indeed, its chirality and linear dispersion at low energies are responsible for its fascinating properties [4] such as Klein tunneling [5], quantum Hall effect [6], and their potential applications in electronic devices, graphene-based nanocomposites, or chemical sensors [7–12]. However, these applications are severely limited by the absence of a gap. Hence, the band-gap opening and the control of graphene bilayer become essential for applications in various electronic devices. One way to create a gap in graphene is the selective functionalization, which has been used, for example, with hydrogen adsorption on a moiré of graphene-Ir(111) [13]. A functionalization by an adatom (or admolecule) covalently bounded to a carbon atom is a resonant scatterer for conduction states which strongly affect electronic structure and transport properties [14–19]. Since graphene is a zero-gap material with a bipartite lattice, such functionalization states

create so-called *midgap states* at the Dirac energy *E*<sub>D</sub>. Bernal bilayer graphene (BLG) is a system formed by two layers of MLG translated from one to the other. One of its advantages is the control of its gap by applying an external gate voltage [20–22], which opens the way to multiple applications for nanodevices [23–25]. On the other hand, the BLG devices can be based on changes in their electrical conductivity, which can be performed with using the influence of substrate [26], vacancies, adatoms, or admolecules adsorbed on the surface of BLG [27–34]. Recently, it has been shown that single- and double-sided fluorination affect strongly conductivity, exhibiting insulating and conducting behavior, respectively [34]. From a theoretical point of view, the study of transport by semiclassical methods has been well done (see for instance Refs. [4,22]). This approach is valid when *E*<sub>F</sub> is far enough from the Dirac energy. But for *E*<sub>F</sub> close to Dirac energy, abnormal transport due to quantum localization has been predicted for a random distribution of absorbates [30] and some very specific cases of selective functionalization [29,31]. These effects are important when the resonant scatterer concentration (defect concentration) is large with respect to the charge carrier concentrations; indeed, each resonant scatterer creates one midgap state at Dirac energy *E*<sub>D</sub>. Since these quantum effects, beyond the semiclassical behavior, are extremely dependent on the type of functionalization, a more systematic theoretical study is still needed to understand current experimental results and stimulate new experimental studies.

<footnote>*jouda.khabthani@fst.utm.tn</footnote>

<footnote>†ahmed.missaoui@cyu.fr</footnote>

<footnote>‡didier.mayou@neel.cnrs.fr</footnote>

<footnote>§guy.trambly@cyu.fr</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Sketch of the crystal structure of AB stacked (Bernal) bilayer graphene (BLG). Atoms A₁ and B₁ are on the lower layer (layer 1); A₂ and B₂ on the upper layer (layer 2).</figcaption></figure>

The unit cell of Bernal BLG contains four carbon atoms, A₁, B₁ in layer 1 and A₂, B₂ in layer 2 (Fig. 1). Atoms A have three B first neighbors in the same layer and one A neighbor in the other layer, while atoms B have only three A first neighbors in the same layer. Thus, the local environment of A and B atoms is different, and the probability that an atom or molecule will stick to an atom A or an atom B should be different. It is thus reasonable to think that the functionalization of B atoms is favored. This simple argument has been confirmed by DFT calculations [35] showing that H adsorption energy difference between A site and B site is about ΔE = 85 meV in favor of B site, when the number of adsorbates is very low. For a larger number of adsorbates, one can therefore expect competition between two contradictory effects: on the one hand preferential adsorption on the B sites of the bilayer, and on the other hand adsorption on different sublattices of the same layer as expected in MLG [35,36]. Indeed in MLG, it exists an interaction between defects states that favors configurations with adsorbates on different sublattices. Such asymmetric adsorption properties between sublattice A and sublattice B have been recently suggested by experimental measurements [32], where the distribution of hydrogen adsorbates on the sublattices is adequately controlled. Overall, BLG lattice is a bipartite lattice of the two sublattices α {A₁, B₂} and β {A₂, B₁}, from which one expects very specific electronic properties produced by selective functionalization. Since BLG is metallic, an isolated functionalization creates an isolated state that is a kind of “midband” state, so-called midgap states by analogy with MLG. In a previous paper [31] we have considered the limiting cases where adsorbates are randomly distributed only on A sublattice or B sublattice of layer 1 while layer 2 remains pristine. On one hand, such a selective functionalization leads to the creation of a gap when sublattice B₁ is functionalized. This gap is a fraction of 1 eV of at least 0.5 eV for a concentration c of adsorbates larger than 1% of the total number of atoms. On the other hand, functionalization of sublattice A₁ decreases the effective coupling between layers, and thus the conductivity increases when c increases, since the pristine layer is less perturbed by the disordered layer when c increases. These two types of selective functionalization exhibit very different and unusual behaviors. This opens the way to the control of electronic properties through selective functionalization, which is experimentally feasible

[32]. However, these extreme cases (A₁ or B₁ functionalization only) seem too simple to correspond to the experimental sample. Indeed, the complexity of the bipartite BLG lattice requires further theoretical studies of other selective adsorbate distributions. This is why it is necessary to study a combined functionalization of several sublattices. In particular, we have to consider cases where midgap states are coupled to each other and thus form a midgap band, leading to new diffusivity properties that are not a simple combination of the extreme situations studied in Ref. [31], in which midgap states are not coupled together.

In this paper we present a detailed study of the electronic structure and quantum transport in BLG with adsorbates (defects) located on two different sublattices among the four sublattices A₁, A₂, B₁, and B₂. We analyze how the symmetry is broken between sublattices under this selectivity, which may cause either a gap or abnormal behavior of the conductivity. We will pay particular attention to cases where B atoms are preferentially functionalized, since these cases should be energetically favorable. For instance, under some specific conditions (adsorbates on B₁ and B₂ sublattices), a spectacular increase of diffusivity of charge carrier of midgap states band edge is obtained when the concentration c of adsorbates increases. The study of conductivity—taking into account all the effects of quantum interference—requires a distinction between several cases, depending on the value of the inelastic mean free path Lᵢ, mainly due to temperature. At high temperatures (typically room temperature) we calculate the microscopic conductivity σ<sub>M</sub>; then we will analyze the quantum corrections at low temperature (very large Lᵢ values), i.e., at the localization regime. In the latter regime we also study how localized states due to defects (midgap states) are at the origin of a particular quantum conductivity that cannot be explained by the Boltzmann’s transport theory, and which is similar to the one found in quasicrystals [37,38], twisted bilayer graphene [39], and recently graphene with defects inducing flat bands [40,41].

The remainder of this paper is organized as follows. Section II introduces the model and the formalism to compute the density of states (DOS) and the conductivity. Sections III and IV focus on selective distributions of vacancies distributed in layer 1 only, and the two layers, respectively. Localization effects on conductivity are discussed in Sec. V. Finally, Sec. VI provides a summary and conclusions.

## II. ELECTRONIC STRUCTURE AND NUMERICAL METHODS

### A. TB Hamiltonian

The tight-binding (TB) Hamiltonian model for BLG with the p<sub>z</sub> orbitals only is given by

$$H = \sum_{(i,j)} t_{ij} |i\rangle\langle j| \quad (1)$$

where i is the index of p<sub>z</sub> orbitals, the sum runs over neighbor sites i, j, and t<sub>ij</sub> is the hopping element matrix between site i and site j. In this paper we consider only the coupling between first neighbors orbitals. There are thus two types of coupling (Fig. 1): for an intralayer coupling term between first
neighboring orbitals $A_1$ and $B_1$ ($A_2$ and $B_2$), $t_{ij} = -\gamma_0 = -2.7$ eV; and for an interlayer coupling term between first neighboring orbitals $A_1$ and $A_2$, $t_{ij} = \gamma_1 = 0.34$ eV [30]. For this kind of calculation, a more realistic TB model with coupling terms above first neighbors leads qualitatively to similar results [30,31]. We have also checked that such a TB model leading to the results presented here are similar, but the first neighbors TB model allows us to better analyze and discuss the physical mechanisms involved as it preserves the electron-hole symmetry. In the Hamiltonian [Eq. (1)], the on-site energies are taken equal to zero so that the Dirac energy $E_D$ is therefore equal to zero.

## B. Adsorbate simulation

We consider that resonant adsorbates are simple atoms or molecules—such as H, OH, $CH_3$—that create a covalent bond with the carbon atom of the BLG. To simulate this covalent bond, we assume that the $p_z$ orbital of carbon, just below the adsorbate, is removed [42–44]. In our calculations the vacancies are randomly distributed in two of the four sublattices $A_1$, $A_2$, $B_1$, and $B_2$, with finite concentration $c$ with respect to the total number of atoms. Here we study all possible cases of the double type of vacancies:

(1) $A_1B_1$-Va: Vacancies randomly distributed on sublattices $A_1$ and $B_1$. An asymmetric distribution, $A_1^x B_1^{1-x}$-Va, where $x$ is the proportion of vacancies in the sublattice $A_1$, is also considered.

(2) $A_1A_2$-Va: Vacancies randomly distributed on sublattices $A_1$ and $A_2$.

(3) $A_1B_2$-Va: Vacancies randomly distributed on sublattices $A_1$ and $B_2$.

(4) $B_1B_2$-Va: Vacancies randomly distributed on sublattices $B_1$ and $B_2$.

In the following, we call X-midgap states the states created by a random distribution of vacant atoms on the X sublattice, with $X = A_1, A_2, B_1, B_2, A_1B_1, A_1A_2, A_1B_2$, or $B_1B_2$.

## C. Quantum transport calculation

We used the real space Kubo-Greenwood (RSKG) method [45–49] which has already been used to study quantum transport in disordered graphene, chemically doped graphene and bilayer (see for instance [15–18,30,31]), functionalized carbon nanotubes [50–52], and many other systems (see for instance the recent review Ref. [53] and references therein). This numerical method connects the DC conductivity $\sigma$, $\sigma = e^2nD$, to the density of states $n$ and the diffusion coefficient

$$D(E, t) = \frac{\Delta X^2(E, t)}{t}, \qquad (2)$$

where the average square spreading $\Delta X^2$ is calculated at every energy $E$ and time $t$ by using the polynomial expansion method [45–49],

$$\Delta X^2(E, t) = \frac{\mathrm{Tr}\{[X, U(t)]^\dagger \delta(E - H) [X, U(t)]\}}{\mathrm{Tr} \delta(E - H)}, \quad (3)$$

where $U(t)$ is the evolution operator at time $t$, $\delta$ is the Dirac function, and Tr is the trace. This numerical approach has the advantage of using efficiently the method in real space. It takes

into account all quantum effects due to a random distribution of static scatterers in a very large supercell containing more than $10^7$ orbitals. Here all calculations are done in a supercell of $1500 \times 1500$ cells of Bernal bilayer (four atoms), with periodic boundary conditions. Considering such a huge cell, it is possible to evaluate the traces TrA in Eq. (3) by the average $\langle A \rangle$ on a random phase state [49]. Such a calculation may be done by the recursion method (Lanczos algorithm) where the Hamiltonian is written as a tridiagonal matrix in real space [54] of dimension $N_r$. Here we use $N_r = 1500$ and we checked that presented results do not change significantly when $N_r$ increases. Lanczos method, which has been used in our previous papers [30,31,55], leads to a convolution of the DOS by a Lorentzian function which a small width $\epsilon$. The DOS is thus obtained by a Lorentzian broadening of the spectrum and $\epsilon$ is a kind of energy resolution of the calculation. But for systems with a gap, to avoid the tail expansion of the Lorentzian function in the gap, it is more efficient to diagonalize the tridiagonal Hamiltonian of dimension $N_r \times N_r$ and to compute the DOS by Gaussian broadening of the spectrum [56]. In the present work, a Gaussian broadening is used with the Gaussian standard deviation of 5 meV. Note that for energies that are not close to the gap the two methods give almost the same results, except for small oscillations that look like regular beatings. These oscillations are numerical artifacts depending on convergence parameters that we used (see the Supplemental Material [57] Sec. S1). They have no effect on the physics discussed here.

The Hamiltonian $H$ [Eq. (1)], written in a supercell, takes into account the effects of elastic collisions (static defects, here vacancies). Therefore, in the framework of a tight-binding model, all quantum effects—including all multiple-scattering effects—are taken into account to calculate the average square spreading $\Delta X^2$ and the diffusive coefficient [Eq. (2)] without inelastic scattering, i.e., at zero temperature. At finite temperature $T$, the inelastic scattering caused by the electron-phonon interactions are implanted by using the approximation of relaxation time approximation (RTA). For details of the implementation of the RTA see the Appendix of Ref. [18]. The conductivity in the $x$ direction is thus given by

$$\sigma(E_F, \tau_i) = e^2 n(E_F) D(E_F, \tau_i), \quad (4)$$

$$D(E_F, \tau_i) = \frac{L_i^2(E_F, \tau_i)}{2\tau_i}, \qquad (5)$$

$$L_i^2(E_F, \tau_i) = \frac{1}{\tau_i} \int_0^\infty \Delta X^2(E_F, t) e^{-t/\tau_i} dt, \quad (6)$$

where $E_F$ is the Fermi energy, $\tau_i$ is the inelastic scattering time, $n(E) = \mathrm{Tr} \, \delta(E - H)$ is the total density of states (total DOS), $D$ the diffusivity along the $x$ axis, and $L_i$ is the inelastic mean free path. $L_i(E_F, \tau_i)$ is the typical distance of propagation during the time interval $\tau_i$ for electrons at energy $E$. $\tau_i$ is the time beyond which the velocity autocorrelation function goes exponentially to zero [18].

$L_i$ is the distance beyond which a wave packet loses its phase coherence due to inelastic scattering processes, whereas elastic scattering events do not destroy the phase coherence. We know that $L_i$ decreases when the temperature $T$ increases,
however, the exact function of $L_i$ versus $T$ is unknown. This is why we consider different cases according to different possible values of $L_i$. Indeed, three different transport regimes may exist depending on $L_i$ value with respect to the elastic mean free path $L_e$, which is the average distance between two elastic scattering events. When $L_i \gg L_e$, multiple scattering effects (such as weak or strong localization) strongly affect the transport and the conductivity is “macroscopic” in the sense that it is established over large sample sizes. This happens at sufficiently low temperature $T$, and then $\sigma$ decreases when $L_i$ increases (i.e., $T$ decreases). For smaller $L_i$ values, since $L_i > \sim L_e$, i.e., larger temperature, $\sigma(L_i)$ reaches a conductivity plateau close to the maximum $\sigma$ value $\sigma_M$, as shown in Sec. V. This regime is called the diffusive regime, where $\sigma(L_i)$ is almost independent on $L_i$ over a large $L_i$ range depending on the energy $E_F$. Examples presented in Sec. V show that the conductivity plateau corresponds to $L_i$ values from few nm to few 10 nm, which may correspond to high temperature and room temperature, respectively. In this case, the conductivity of a sample depends only on the quantum scattering over small distances which are typically of the order of magnitude of the distances between static defects ($L_e$); this is the reason why we call $\sigma_M$ the “microscopic” conductivity. The situation $L_i < L_e$ is an extreme case that one should not often reach in real materials. This corresponds to the case of very pure materials with very few static defects. The conductivity is independent of static defects, and thus $\sigma(L_i)$ increases when $L_i$ increases.

At each energy, the microscopic diffusivity $D_M$ and microscopic conductivity $\sigma_M$ are defined as the maximum value of $D(\tau_i)$ and $\sigma(\tau_i)$, respectively. It is also interesting to have an estimate of the $L_e$ values, and the $L_i$ values corresponding to the diffusive regime, i.e., $\sigma(L_i) \simeq \sigma_M$. We compute the elastic mean free path $L_e$ along the $x$ axis, from the usual phenomenological formula [18],

$$L_e(E) = \frac{1}{V_0(E)} \text{Max}_{\tau_i} \left\{ \frac{L_i^2(E, \tau_i)}{\tau_i} \right\} = \frac{2D_M(E)}{V_0(E)}, \quad (7)$$

where the velocity $V_0$ is the slope of $L_i(\tau_i)$ at very small $\tau_i$. It is important to note that such a definition of $L_e$ is not very accurate, and this calculation can only give an order of magnitude of the average distance between two elastic scattering events. Indeed, the formula (7) is not always valid when the electronic structure is strongly modified by static defects. Moreover, $V_0$ is overestimated since the numerical calculations include not only the intraband terms but also the interband terms. In the case of graphene monolayer, we have shown [39] that these latter increase $V_0$ by a factor of $\sqrt{2}$ which leads to an underestimation of the $L_e$. However, roughly speaking, $L_e$ is the $L_i$ value above which conductivity curve $\sigma(L_i)$ reaches the plateau of a diffusive regime due to elastic scattering. To better define the $L_i$ values corresponding to the diffusive regime, we define the lengths $L_{i1}$ and $L_{i2}$ such as: $\forall L_i \in [L_{i1}; L_{i2}]$, $\sigma(L_i) > 0.9\sigma_M$. We also determine the value $L_{im}$ such as $\sigma$ is maximum, i.e., $\sigma(L_{im}) = \sigma_M$. The values of $L_e, L_{i1}, L_{im}$, and $L_{i2}$ are shown in Fig. S4 in the Supplemental Material [57] for different concentrations of the four types of vacancies studied. The results show that $L_e \ll L_{i1}$ with the same order of magnitude, and the ratio $L_{i2}/L_{i1}$ varies from

5–10 to very large values, depending on the type of defects and their concentrations.

Microscopic conductivity, which corresponds to the situation where $\sigma(L_i) \simeq \sigma_M$, i.e., large (or room) temperature limit, is analyzed in Secs. III and IV. The $L_i \gg L_e$ limit, i.e., $\sigma(L_i) < \sigma_M$, which corresponds to the localization regime at low temperature, is analyzed in Sec. V.

### III. VACANCIES IN ONE LAYER ONLY

In this section we are focusing on the impact of the vacancies distributed on one layer (layer 1) of BLG. It should simulate adsorbates or defects that come from the preparation process [19] or induced by the substrate [58]. For example, in epitaxial graphene on Pt(111) [58], the authors have shown the appearance of covalent bonds between the carbon atoms of graphene and the atoms of Pt. Since the $B_1$ atoms of layer 1 do not have a first neighbor in layer 2, it is likely that their functionalization is favored, but the experimental results [32] do not show functionalization only on B atoms. It is thus important to study an asymmetric functionalization of $B_1$ or $A_1$ sublattice. We first consider a majority functionalization of the $B_1$ atoms ($A_1$ atoms), and we analyze the effect of defect concentrations on a symmetric distribution of vacancies.

#### A. $A_1B_1$-Va asymmetrically distributed

We consider an asymmetric distribution of vacancies: $A_1^x B_1^{1-x}$-Va, where $x(1-x)$ is the proportion of vacancies on sublattice $A_1$ ($B_1$). Considering the cases with a total number of vacancies corresponding to a concentration $c = 3\%$ with respect to the total number of atoms, the density of states $n(E)$ and the microscopic conductivity $\sigma_M(E)$ are shown in Fig. 2 for different $x$ values. As presented in Fig. S5 of the Supplemental Material [57], the results for $c = 0.5\%$ show very similar behaviors.

The different disorder distributions, i.e., value of $x$ between $x = 0$ ($B_1$ vacancies only) and $x = 1$ ($A_1$ vacancies only), affect strongly the regime around the Dirac energy $E_D$. Midgap states at $E_D$ always appeared in both layers. Indeed, each $A_1$ missing orbital of layer 1 produces a $A_1$-midgap state at Dirac energy $E_D$ that spread on $B_1$ sublattice (layer 1) only, and $B_1$ missing orbital produces a $B_2$-midgap state that spread on $A_1$ (layer 1) and $B_2$ (layer 2) sublattices [31]. $A_1$-midgap states and $B_1$-midgap states are coupled by the Hamiltonian and form a band of midgap states with specific transport properties. In the extreme cases of vacancies distributed over a single sublattice $B_1$ ($x = 0$), we have shown [31] that a gap around the Dirac energy $E_D$ is created. This gap is a consequence of the reduction of the average number of neighbors of atoms in a sublattice. For intermediate $x$ values, the gap disappears under the effect of the interactions between midgap states. Depending on $x$ values, two scenarios emerge:

(i) For $x \in [0; 0.3]$ and $x \in [0.7; 1]$, the number of $A_1$-midgap states and $B_1$-midgap states are rather different, and many of those states are not coupled to each other and remain isolated with energy $E_D$. The small number of mixed midgap states leads to a small DOS at intermediates energies [Fig. 2(a)].
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 2. BLG with $A_1^x B_1^{1-x}$-Va for different distributions $x$ of vacancies between $A_1$ and $B_1$ sites: (a) and (b) $x \in [0; 0.5]$ (mainly $B_1$-Va) and (c) and (d) $x \in [0.5; 1]$ (mainly $A_1$-Va). (a)–(c) Density of states $n(E)$, the integrated density of states is represented in the left inset while the density of states around the Dirac energy $E_D$ is in the right inset. (b)–(d) Microscopic conductivity $\sigma_M(E)$ for the same disorder configurations. The total concentration of vacancies is 3%. $G_0 = 2e^2/h$.</figcaption></figure>

Concerning the conductivity, two different behaviors are obtained according to the dominant concentration of $B_1$ vacancies ($x \in [0; 0.3]$) or $A_1$ vacancies ($x \in [0.7; 1]$). The behavior of $\sigma_M(E)$ around Dirac energy for $x \in [0; 0.3]$ is determined mainly by the effects of the $B_1$ vacancies. For energies $E$ in the intermediate regime with $E \le \gamma_1 = 0.34$ eV, $\sigma_M$ increases when the coupling between midgap states increases, i.e., when $A_1$ and $B_1$ vacancy concentrations are close to each other. For $x \in [0.7; 1]$, results are very sensitive to the concentration of $A_1$ vacancies. $\sigma_M$ increases when $x$ increases. This effect of $A_1$ vacancies affects the microscopic conductivity on a range of energy that does not exceed 1 eV as it is shown in Fig. 2(b). In the extreme case $x = 1$, a gap appears in the average DOS for the layer with defects (layer 1) [31]. It is proportional to the concentration $c$ of vacancies and layer 2 behaves more and more like a pristine MLG which gives the ballistic behavior. When $x$ is close to 1, $x \lesssim 1$, the gap in layer 1 disappears, and thus the microscopic conductivity increases when $x$ (close to 1) increases.

(ii) The interactions between midgap states are important for $x \in [0.4; 0.6]$, and it is maximum for $x = 0.5$. Therefore $n(E)$ is larger for energy $E \ne E_D$ [right inset of Fig. 2(a)].

The conductivity behavior is similar to that found in the following section for $x = 0.5$.

## B. $A_1B_1$-Va symmetrically distributed

We now study a random distribution of defects equally distributed in sublattices $A_1$ and $B_1$, labeled $A_1B_1$-Va. Total DOS $n(E)$, LDOS, and microscopic conductivity $\sigma_M(E)$ are shown in Fig. 3 for several values of vacancy concentrations $c$ with respect to the total number of atoms. Since the electron transport through the BLG is mainly determined by the electrons which have energy close to the Dirac point, the conductivity is displayed within a small energy region around the charge neutrality energy $E_D = 0$. By inspecting Figs. 3(a)–3(c), one can identify several important features. (i) For all concentrations $c$ and energy around $E_D$, $0.02 < |E - E_D| < 0.1$ eV, $\sigma_M$ presents a minimum plateau at conductivity $\sigma_M \approx 1.2 G_0$, with $G_0 = 2e^2/h$. Thus $\sigma_M \approx 2\sigma_M^{\text{mono}}$, where $\sigma_M^{\text{mono}} \approx 0.6 G_0$ is the monolayer graphene (MLG) microscopic conductivity [18,30,59,60]. This shows that the defects affect both planes similarly, although one of the two planes is defect-free. Moreover, the presence of a plateau, almost independently of the
<figure><img src="image_8.png" /><figcaption>FIG. 3. Electronic properties in BLG with A₁B₁ vacant atoms (A₁B₁-Va), with equal distribution of vacancies between A₁ and B₁ sublattices: (a) total DOS (dashed line is the total DOS without vacancies), (b) average local DOS on A₁, B₁, A₂, B₂ atoms for c = 0.25% (dashed line and dot line are LDOS on A and B atom without vacancies), and (c) microscopic conductivity σ<sub>M</sub>(E). c is the concentration of vacancies with respect to the total number of atoms in BLG. G₀ = 2e²/h.</figcaption></figure>

concentration, shows that the microscopic quantities in the
BLG are not affected directly by interlayer coupling terms,
which gives them a behavior similar to MLG. This behavior
is understandable since the elastic mean free path L<sub>e</sub> (see the
Supplemental Material [57] Figs. S3 and S4) is smaller than
the traveling distance l<sub>1</sub> in a layer between two interlayer
hoppings, l<sub>1</sub> ≃ 1–2 nm [30]. (ii) For energies far from E<sub>D</sub>,
|E − E<sub>D</sub>| > 0.1 eV, two behaviors of the conductivity is ob-

served: for c ≤ 2%, σ<sub>M</sub> ≃ σ<sub>B</sub>, where σ<sub>B</sub> is calculated with the
Bloch-Boltzmann approach [22,61], and then conductivity is
proportional to 1/c. While for c ≥ 2%, σ<sub>M</sub> seems to depend
less on c, and even slightly increases when c increases, such
as for A<sub>1</sub> vacancies alone or B<sub>1</sub> vacancies alone [31].

IV. VACANCIES IN BOTH LAYERS

In this section we study the combined effect of vacan-
cies distributed in two sublattices that do not belong to the
same layer. The case B₁B₂-Va, which seems to be the most
favored case for functionalization, is considered first. These
midgap states are coupled to each other and form a midgap
band characterized by a very unusual quantum diffusion of
charge carriers. After, we study the cases of A₁A₂-Va and
A₁B₂-Va that both produce uncoupled midgap states at energy
E = E<sub>D</sub> = 0.

A. B₁B₂-Va cases

B₁- and B₂-midgap states are distributed over all the struc-
ture with different weights on the atoms A₁, A₂, B₁, and B₂
[Fig. 4(b)]. They form a band since B₁-Va midgap states and
B₂-Va midgap states are coupled by the Hamiltonian. Their
electronic properties are thus very different from those of B₁
vacancies in BLG for which a gap proportional to c is formed
around E<sub>D</sub> [31]; while with B₁B₂-Va, the B₁- and B₂-midgap
states are coupled, and thus the gap is filled or partially filled
by a midgap states band. Several regimes are present depend-
ing on both energy E and vacancy concentration c.

For small concentrations c, typically c ⩽ 1%, there is no
gap in the DOS [Fig. 4(a)] and states around E<sub>D</sub> form a
narrow midgap states band. The corresponding microscopic
conductivity σ<sub>M</sub> presents a plateau [see the inset Fig. 4(c)] at
a value independent on c, σ<sub>M</sub> ≃ 2σ<sub>M</sub><sup>mono</sup>.

For high concentrations c, the density of states [Fig. 4(a)]
around E<sub>D</sub> increases significantly, and as a direct conse-
quence, the plateau of conductivity increases σ<sub>M</sub> > 2σ<sub>M</sub><sup>mono</sup>.
As explained above (Sec. IV B), in each layer the gap due
to B-Va increases when c increases, therefore for large c the
midgap states bandwidth becomes smaller than the gap, and
the midgap states band becomes isolated from other states
by small gaps at |E| ≳ γ<sub>1</sub> [Fig. 4(a)]. The width of this
isolated band is Δw ≃ 2γ<sub>1</sub>, i.e., E ∈ [−γ<sub>1</sub>, γ<sub>1</sub>]. For large con-
centrations c, the edge states (E ≃ ±γ<sub>1</sub>) have a very exotic
conductivity σ<sub>M</sub> which strongly increases when c increases,
whereas DOS does not change too much. Roughly speaking
this spectacular behavior can be explained by considering the
coupling between the B<sub>1</sub>-Va monolayer midgap states and the
B<sub>2</sub>-Va monolayer midgap states. In monolayer, B-Va midgap
states are located on the A sublattice around the B vacancy.
B-Va midgap states of each layer are not coupled to each other.
But, since each A orbital are coupled with an A orbital of
the other layer, a B<sub>1</sub>-Va midgap state is coupled with a B<sub>2</sub>-Va
midgap state, with a hopping term γ<sub>B<sub>1</sub>-B<sub>2</sub></sub>. γ<sub>B<sub>1</sub>-B<sub>2</sub></sub> ≃ γ<sub>1</sub>, for the
smallest d<sub>B<sub>1</sub>-B<sub>2</sub></sub> distance between the B<sub>1</sub>-Va and the B<sub>2</sub>-Va
(typically first neighbor B<sub>1</sub>-B<sub>2</sub>), and γ<sub>B<sub>1</sub>-B<sub>2</sub></sub> decreases when
d<sub>B<sub>1</sub>-B<sub>2</sub></sub> increases. When c increases, the average d<sub>B<sub>1</sub>-B<sub>2</sub></sub> distance
decreases and thus the average γ<sub>B<sub>1</sub>-B<sub>2</sub></sub> value increases. As a
result, by a kind of percolation between monolayer B-midgap
<figure><img src="image_9.png" /><figcaption>FIG. 4. Electronic properties in BLG with $B_1B_2$ vacant atoms: (a) Total DOS (dashed line is the total DOS without vacancies), and (b) average local DOS on $A_1$, $B_1$, $A_2$, $B_2$ atoms for $c = 0.25\%$ (dashed line and dot line are LDOS on A and B atom without vacancies). The average local DOS on $A_2$, $B_2$ atoms is obtained by a symmetry with relative atoms $A_1$, $B_1$, respectively. (c) Microscopic conductivity $\sigma_M(E)$. $c$ is the concentration of vacancies with respect to the total number of atoms in BLG. $G_0 = 2e^2/h$.</figcaption></figure>

states of the two layers, the conductivity through the BLG increases strongly when $c$ increases.

Finally, the presence of the conductivity plateau for all concentrations [inset Fig. 4(c)] can be understood considering the elastic mean free path $L_e$ shown in the Supplemental Material [57] (Figs. S3 and S4). Around $E_D$ energy ($E \in [-0.2; 0.2]$ eV), $L_e < l_1$, where $l_1 \simeq 1-2$ nm is the traveling distance between two interlayer hopping events [30]. Thus, the diffusion of the charge carriers is not affected by the

interlayer coupling. The diffusive regime is reached in each layer independently, and it takes the MLG minimum value in each layer. Note that like for other types of vacancies, for energy away from Dirac energy, $|E - E_D| \gg \gamma_1$, Boltzmann behavior is always found.

## B. $A_1A_2$-Va and $A_1B_2$-Va cases

The double-type vacancies: $A_1A_2$-Va (vacancies randomly distributed on $A_1$ and $A_2$ sublattices) and $A_1B_2$-Va (vacancies randomly distributed on $A_1$ and $B_2$ sublattices) are characterized by the absence of coupling between midgap states and thus all midgap states remain at energy $E_D = 0$. Indeed, in the case of $A_1A_2$-Va, $N$ vacancies on atoms $A_1$ ($A_2$) sublattice produce a set of $N$ uncoupled midgap states at Dirac energy $E_D = 0$ that are located on the orbitals $B_1$ ($B_2$) of the same layer [31]. As $B_1$ orbitals and $B_2$ orbitals are not directly coupled by the Hamiltonian, midgap states located on $B_1$ and $B_2$ sublattices are not coupled together. In the case $A_1B_2$-Va, vacancies are vacant atoms of the same sublattice $\alpha$ of the BLG lattice. Corresponding midgap states are thus uncoupled states at $E_D$, located on the $\beta$ sublattice with a greater weight on the $B_1$ atoms. For clarity these isolated states at $E_D = 0$ are not shown in the DOSs drawn Figs. 5 and 6 (see the Supplemental Material [57] Sec. S1).

In the $A_1A_2$-Va case, $A_1$ vacancies and $A_2$ vacancies act on both layers symmetrically and independently because the midgap states of a layer are not coupled with midgap states of the other layer. Thus, the result is simply the sum of two independent MLG. In MLG, vacancies in sublattice A (B) produce midgap states at $E_D$ that are located in sublattice B (A). As shown in our previous paper [31] by an analysis of the spectrum of bipartite Hamiltonian, when the concentration $c$ of vacancies increases, a gap increases around the Dirac energy. This gap is a consequence of the reduction of the average number of neighboring atoms of sublattice's atoms which do not contain vacancies. Thus, $A_1$ vacancies ($A_2$ vacancies) create a gap in layer 1 (layer 2) as it is clearly shown on the local DOS of atoms $A_1$ and $B_1$ [Fig. 5(b)]. The total DOS has a gap proportional to the concentration of vacancies $c$ around the Dirac energy $E_D$ [Fig. 5(a)].

$A_1B_2$-Va create also a gap because they are distributed randomly on the same sublattice $\alpha$ {$A_1B_2$} of BLG. Total and local DOSs [Figs. 6(a) and 6(b)] confirm the presence of a gap around Dirac energy $E_D$.

The microscopic conductivity $\sigma_M(E)$ for both types of vacancies $A_1A_2$-Va and $A_1B_2$-Va are shown in Figs. 5(c) and 6(c), respectively. The midgap states at energies $E = E_D$ do not contribute to the conductivity $\sigma_M$ since they are isolated localized states around each vacancy. Beyond the gap, $\sigma_M$ decreases when $c$ increases, following a typical Boltzmann behavior [61].

## V. CONDUCTIVITY VERSUS INELASTIC SCATTERING

In the two previous sections we have studied the microscopy conductivity $\sigma_M$ which is equal to the maximum value of $\sigma_M(\tau_i)$ (Sec. IIC). We now consider $\sigma$ versus the inelastic scattering time $\tau_i$ or the inelastic scattering length $L_i$. Indeed, the inelastic scattering events, which depend on the temperature, can lead to new behaviors at low temperature
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 5. Electronic properties in BLG with A₁A₂ vacant atoms: (a) Total DOS (dashed line is the total DOS without vacancy), (b) average local DOS on A₁, B₁ atoms for *c* = 0.25% (dashed line and dotted line are LDOS on A and B atom without vacancy). The average local DOS on A₂, B₂ atoms is obtained by a symmetry with relative atoms A₁, B₁, respectively. (c) Microscopic conductivity σ<sub>M</sub>(E). *c* is the concentration of vacancies with respect to the total number of atoms in BLG. For clarity the midgap states at *E*<sub>D</sub> = 0 are not shown (see text). *G*₀ = 2e²/h.</figcaption></figure>

due to the multiple scattering, i.e., when $L_e \ll L_i$. This reveals
new quantum effects such as the Anderson localization and the
universal conductivity of the midgap states.

A. Anderson localization

In the framework of the relaxation time approximation
(RTA) (Sec. II C), we compute the inelastic mean free path

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 6. Electronic properties in BLG with A₁B₂ vacant atoms: (a) Total DOS (dashed line is the total DOS without vacancy), (b) average local DOS on A₁, B₁, A₂, B₂ atoms for *c* = 0.25% (dashed line and dotted line are LDOS on A and B atom without vacancy), and (c) microscopic conductivity σ<sub>M</sub>(E). *c* is the concentration of vacancies with respect to the total number of atoms in BLG. For clarity the midgap states at *E*<sub>D</sub> = 0 are not shown (see text). *G*₀ = 2e²/h.</figcaption></figure>

$L_i(E, \tau_i)$ at every energy $E$ and inelastic scattering times $\tau_i$ (Sec. II C). Figure 7 shows the conductivity $\sigma$ drawn versus $L_i$ for different types of vacancies and different energies close to $E_D$. The microscopic conductivity $\sigma_M(E)$ discussed in previous sections (Secs. III and IV) is the maximum value of the curves $\sigma(L_i)$ at the corresponding energy $E$. Each curve $\sigma(L_i)$ has three parts. (1) For small $L_i$, typically $L_i \ll L_e$, the static defects have no direct effect and $\sigma \propto L_i$. This regime is possible at finite temperature only when the defect concen-
<figure><img src="image_16.png" /><figcaption>FIG. 7. Conductivity $\sigma$ as a function of inelastic scattering length $L_i$ for $c = 1\%$. (a) Vacancies randomly distributed on atoms A₁ and B₁, (b) vacancies randomly distributed on atoms B₁ and B₂, (c) vacancies randomly distributed on atoms A₁ and A₂, and (d) vacancies randomly distributed on atoms A₁ and B₂. $G_0 = 2e^2/h$. The vertical dashed lines show the value of $L_e$ calculated by Eq. (7) for each energy value. In (a) and (b), the black dashed lines are the extrapolation of $\sigma(L_i)$ curves, using Eq. (8), to find the localization length $\xi$ at the limit: $\sigma(L_i = \xi) = 0$.</figcaption></figure>

tration is extremely low. (2) For $L_i > \sim L_e$, $\sigma(L_i)$ reaches a plateau at $\sim \sigma_M$. For small defect concentrations $c$, this regime can be found for a wide range of $L_i$ values. (3) For large $L_i$ values, $L_i \gg L_e$, localization regime is reached and $\sigma(L_i)$ decreases when $L_i$ increases. In this regime, the so-called quantum corrections $\Delta\sigma(L_i) = \sigma(L_i) - \sigma_M$ govern the transport properties.

Inelastic scattering collisions are mainly due to electron-phonon interactions, and thus $L_i$ decreases when the temperature $T$ increases. Realistic $L_i$ values are difficult to known, but it is reasonable to consider that at room temperature and higher temperature $L_i$ is such as $\sigma(L_i) \simeq \sigma_M$ (plateau regime) and thus the quantum corrections are negligible. At low temperatures, i.e., $L_i \gg L_e$, quantum interferences dominate transport properties.

In 2D materials, Anderson localization due to quantum interferences leads to a conductivity varying linearly with $\ln L_i$ [62], and can be written [18,63],

$$ \sigma(E, L_i) = \sigma_0(E) - \alpha G_0 \ln \left( \frac{L_i}{L_e(E)} \right), \quad (8) $$

where $G_0 = 2e^2/h$, and $\sigma_0$ values are on the range of $\sigma_M$ values. The second term of the right side of Eq. (8) is the quantum correction of the conductivity. The linear behavior of $\sigma(L_i)$ is clearly seen for cases A₁B₁-Va and B₁B₂-Va [Figs. 7(a) and 7(b)].

The fit of the $\sigma(L_i)$ curve for large $L_i$ gives the $\alpha$ value $\alpha \simeq 0.34$. This value is close to the result found in MLG [18], BLG with random vacancies [30], twisted bilayer graphene [55], and close to the prediction of perturbation theory of 2D Anderson localization [62], for which $\alpha = 1/\pi$. The localization length $\xi$ can be extracted from the expression (8) by extrapolation of $\sigma(L_i)$ curves [Figs. 7(a) and 7(b)] when $\sigma(L_i = \xi) = 0$, giving the following expression:

$$ \xi(E) = L_e(E) \exp \left( \frac{\sigma_0(E)}{\alpha G_0} \right). \quad (9) $$

Since $\alpha$ is a constant, this leads to a simple relationship between $\xi$ and $L_e$, $\xi \simeq 50L_e$, which is between the monolayer graphene value with random vacancy distributions ($13L_e$) [18] and that of BLG in the same case ($13^2L_e$) [30].

For A₁A₂-Va and A₁B₂-Va cases, at energies around the edge of the gap [Figs. 7(c) and 7(d)], the decrease of $\sigma(L_i)$ does not follow Eq. (8). This behavior is more similar to what is generally expected for the conduction by midgap states of graphene [18], which are very localized states with abnormal diffusion behavior.

## B. Universal conductivity of the midgap states

It is also interesting to focus on the conduction by flat band midgap state themselves, i.e., here midgap states at energy
$E_D = 0$ that are not coupled to each other by the Hamiltonian ($A_1A_2$-Va and $A_1B_2$-Va cases). In these midgap states, the average velocity is zero but conduction is possible due to the inherent quantum fluctuations of the velocity which are due to the interband contributions of the velocity correlation function [37,39–41]. Indeed, in the presence of inelastic scattering these fluctuations are modified [41] and do not cancel completely at large times which allows electronic diffusion. It results in non-Boltzmann conductivity, similar to the one found in quasicrystals [37,38], twisted bilayer graphene at the magic angle [39], and graphene with particular defects inducing flat bands [40,41]. In $A_1A_2$-Va and $A_1B_2$-Va, microscopic conductivity, i.e., small inelastic mean free time $\tau_i$, at midgap-states energy is negligible. But at large $\tau_i$ (large $L_i$), the Kubo-Greenwood conductivity of midgap states is [41]

$$\sigma(E, \tau_i) = e^2 n_i(E, \tau_i) D(E, \tau_i), \quad (10)$$

where $n_i$ and $D$ are the DOS and the diffusivity [Eq. (5)] in the presence of inelastic scattering. Since midgap states are nondispersive states at $E = 0$, isolated by gaps (cases $A_1A_2$-Va and $A_1B_2$-Va), $n_i$ is the broadening of the Delta function $c\delta(E)$ by a Lorentzian with a width at half-maximum $\eta$, $\eta = \hbar/\tau_i$. Thus at Dirac energy $E_D = 0$,

$$\sigma(E = 0, \tau_i) = \frac{16}{S} G_0 c \tau_i D(E = 0, \tau_i), \quad (11)$$

where $S$ is the surface of the unit cell. As shown in Fig. 8, for large $\tau_i$, $\sigma(E = 0, \tau_i)$ reaches a constant universal value, independent of the defect concentration $c$, which is twice that of graphene [41]: $\sigma(E = 0) \simeq 1.3 G_0$. As shown in Sec. S5 of the Supplemental Material [57], similar behavior is also seen for the midgap states of $A_1$-Va only and $B_1$-Va only (Fig. S6).

## VI. CONCLUSION

We have studied numerically the effects on the electronic properties of selective functionalization distributed over different sublattices of the Bernal bilayer graphene (BLG). We consider the covalent adsorptions of atoms or molecules. For Fermi energy $E_F$ far from Dirac energy, typically corresponding to a charge carrier concentration greater than the defect concentration $c$, the adsorbates act as weak scatterers, and the usual semiclassical transport calculations are possible. But for smaller doping, typically when the doping is smaller than $c$, $E_F$ is close to Dirac energy the quantum effects—such as midgap states or midgap-band gap, unusual localization—dominate transport properties. Our numerical approach includes all these quantum effects.

We prove theoretically that controlled functionalization can be an excellent way to tune BLG conductivity. This is in agreement with recent experimental results [32,34] showing that it is possible to control the functionalization with an adsorbate rate of the order of 1% of the total number of atoms and to fabricate single and double side adsorbed bilayer graphene. We find a wide variety of original behaviors and have classified them according to the functionalized sublattices, the adsorbate concentration $c$, and the energy. For example, we give the conditions for opening a mobility gap of several 100 meV. Experimentally, and according to Ref. [32], the hydrogen adsorption on the B atoms in one

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>FIG. 8. Conductivity $\sigma(E = E_D = 0)$ as a function of inelastic scattering time $\tau_i$. (a) Vacancies randomly distributed on atoms $A_1$ and $A_2$. (b) Vacancies randomly distributed on atoms $A_1$ and $B_2$. In both cases midgap states are uncoupled states at $E_D = 0$ isolated by gaps. $G_0 = 2e^2/h$.</figcaption></figure>

layer is energetically favored. For this reason, the study of the specific cases of $B_1B_2$ adsorbates is very interesting. An isolated midgap states band is predicted. Spectacularly, for $c \sim 1\%$, its edge states have a high electrical conductivity due to the large diffusivity of charge carriers, which deserves further investigation. As the functionalization of atoms can be performed experimentally, one can imagine that those of the $B_1B_2$ adsorbates can also be carried out, which makes it possible to control the conductivity.

The present study contributes to the understanding the electronic properties of localized states—flat bands—due to the combined effect of quantum interferences and geometrical properties (here bipartite lattice). This physics of flat bands is currently a major one in condensed matter, either for field topological insulators or for remarkable electronics (correlation effect, superconductivity) of the moiré flat bands in twisted bilayer graphene at magic angle [64,65].

## ACKNOWLEDGMENTS

The authors wish to thank G. Bouzerar, L. Magaud, P. Mallet, G. Jemaï, and J.-Y. Veuillen for fruitful discussions. Calculations have been performed at the Centre de Calculs (CDC), CY Cergy Paris Université, and using HPC resources from GENCI-IDRIS (Grant No. 910784). We thank Y. Costes
and B. Mary, CDC, for computing assistance. This work was supported by the ANR project J2D (ANR-15-CE24-0017) and

the Paris//Seine excellence initiative (Grant No. 2019-055-C01-A0).

[1] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Electric field effect in atomically thin carbon films, *Science* **306**, 666 (2004).

[2] C. Berger, Z. Song, T. Li, X. Li, A. Y. Ogbazghi, R. Feng, Z. Dai, A. N. Marchenkov, E. H. Conrad, P. N. First, and W. A. de Heer, Ultrathin epitaxial graphite: 2D electron gas properties and a route toward graphene-based nanoelectronics, *J. Phys. Chem. B* **108**, 19912 (2004).

[3] A. Hashimoto, K. Suenaga, A. Gloter, K. Urita, and S. Iijima, Direct evidence for atomic defects in graphene layers, *Nature* (London) **430**, 870 (2004).

[4] S. Das Sarma, S. Adam, E. H. Hwang, and E. Rossi, Electronic transport in two-dimensional graphene, *Rev. Mod. Phys.* **83**, 407 (2011).

[5] M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, Chiral tunnelling and the Klein paradox in graphene, *Nat. Phys.* **2**, 620 (2006).

[6] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Two-dimensional gas of massless Dirac fermions in graphene, *Nature* (London) **438**, 197 (2005).

[7] F. Schedin, A. K. Geim, S. V. Morozov, E. W. Hill, P. Blake, M. I. Katsnelson, and K. S. Novoselov, Detection of individual gas molecules adsorbed on graphene, *Nat. Mater.* **6**, 652 (2007).

[8] J. Wu, H. A. Becerril, Z. Bao, Z. Liu, Y. Chen, and P. Peumans, Organic solar cells with solution-processed graphene transparent electrodes, *Appl. Phys. Lett.* **92**, 263302 (2008).

[9] M. A. Rafiee, J. Rafiee, Z. Wang, H. Song, Z.-Z. Yu, and N. Koratkar, Enhanced mechanical properties of nanocomposites at low graphene content, *ACS Nano* **3**, 3884 (2009).

[10] S. Stankovich, D. A. Dikin, G. H. B. Dommett, K. M. Kohlhaas, E. J. Zimney, E. A. Stach, R. D. Piner, S. T. Nguyen, and R. S. Ruoff, Graphene-based composite materials, *Nature* (London) **442**, 282 (2006).

[11] X.-M. Huang, L.-Z. Liu, S. Zhou, and J.-J. Zhao, Physical properties and device applications of graphene oxide, *Front. Phys.* **15**, 33301 (2020).

[12] L. Chen, F. Ouyang, S. Ma, T.-F. Fang, A.-M. Guo, and Q.-F. Sun, Enhancement of electron transport and band gap opening in graphene induced by adsorbates, *Phys. Rev. B* **101**, 115417 (2020).

[13] J. H. Jørgensen, A. G. Čabo, R. Balog, L. Kyhl, M. N. Groves, A. M. Cassidy, A. Bruix, M. Bianchi, M. Dendzik, M. A. Arman, L. Lammich, J. I. Pascual, J. Knudsen, B. Hammer, P. Hofmann, and L. Hornekaer, Symmetry-driven band gap engineering in hydrogen functionalized graphene, *ACS Nano* **10**, 10798 (2016).

[14] N. Leconte, D. Soriano, S. Roche, P. Ordejon, J.-C. Charlier, and J. J. Palacios, Magnetism-dependent transport phenomena in hydrogenated graphene: From spin-splitting to localization effects, *ACS Nano* **5**, 3987 (2011).

[15] A. Lherbier, S. M.-M. Dubois, X. Declerck, Y.-M. Niquet, S. Roche, and J.-C. Charlier, Transport properties of graphene containing structural defects, *Phys. Rev. B* **86**, 075402 (2012).

[16] S. Roche, N. Leconte, F. Ortmann, A. Lherbier, D. Soriano, and J.-C. Charlier, Quantum transport in disordered graphene: A theoretical perspective, *Solid State Commun.* **152**, 1404 (2012).

[17] A. Cresti, F. Ortmann, T. Louvet, D. Van Tuan, and S. Roche, Broken Symmetries, Zero-Energy Modes, and Quantum Transport in Disordered Graphene: From Supermetallic to Insulating Regimes, *Phys. Rev. Lett.* **110**, 196601 (2013).

[18] G. Trambly de Laissardière and D. Mayou, Conductivity of Graphene with Resonant and Nonresonant Adsorbates, *Phys. Rev. Lett.* **111**, 146601 (2013).

[19] G. Yang, L. Li, W. B. Lee, and M. Cheung Ng, Structure of graphene and its disorders: A review, *Sci. Technol. Adv. Mater.* **19**, 613 (2018).

[20] 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, and A. H. Castro Neto, Biased Bilayer Graphene: Semiconductor with a Gap Tunable by the Electric Field Effect, *Phys. Rev. Lett.* **99**, 216802 (2007).

[21] E. McCann and V. I. Fal'ko, Landau-Level Degeneracy and Quantum Hall Effect in a Graphite Bilayer, *Phys. Rev. Lett.* **96**, 086805 (2006).

[22] E. McCann and M. M. Koshino, The electronic properties of bilayer graphene, *Rep. Prog. Phys.* **76**, 056503 (2013).

[23] Y. Zhang, T.-T. Tang, C. Girit, Z. Hao, M. C. Martin, A. Zettl, M. F. Crommie, Y. R. Shen, and F. Wang, Direct observation of a widely tunable bandgap in bilayer graphene, *Nature* (London) **459**, 820 (2009).

[24] H. Overweg, A. Knothe, T. Fabian, L. Linhart, P. Rickhaus, L. Wernli, K. Watanabe, T. Taniguchi, D. Sánchez, J. Burgdörfer, F. Libisch, V. I. Fal'ko, K. Ensslin, and T. Ihn, Topologically Nontrivial Valley States in Bilayer Graphene Quantum Point Contacts, *Phys. Rev. Lett.* **121**, 257702 (2018).

[25] A. Kurzmann, M. Eich, H. Overweg, M. Mangold, F. Herman, P. Rickhaus, R. Pisoni, Y. Lee, R. Garreis, C. Tong, K. Watanabe, T. Taniguchi, K. Ensslin, and T. Ihn, Excited States in Bilayer Graphene Quantum Dots, *Phys. Rev. Lett.* **123**, 026803 (2019).

[26] S. Y. Zhou, G.-H. Gweon, A. V. Fedorov, P. N. First, W. A. de Heer, D.-H. Lee, F. Guinea, A. H. Castro Neto, and A. Lanzara, Substrate-induced bandgap opening in epitaxial graphene, *Nat. Mater.* **6**, 770 (2007).

[27] O. Leenaerts, B. Partoens, and F. M. Peeters, Hydrogenation of bilayer graphene and the formation of bilayer graphane from first principles, *Phys. Rev. B* **80**, 245422 (2009).

[28] R. E. Mapasha, A. M. Ukpong, and N. Chetty, *Ab initio* studies of hydrogen adatoms on bilayer graphene, *Phys. Rev. B* **85**, 205402 (2012).

[29] D. Van Tuan and S. Roche, Anomalous ballistic transport in disordered bilayer graphene: A Dirac semimetal induced by dimer vacancies, *Phys. Rev. B* **93**, 041403(R) (2016).

[30] A. Missaoui, J. J. Khabthani, N.-E. Jaidane, D. Mayou, and G. Trambly de Laissardière, Numerical analysis of electronic conductivity in graphene with resonant adsorbates: Comparison of monolayer and Bernal bilayer, *Eur. Phys. J. B* **90**, 75 (2017).
[31] A. Missaoui, J. J. Khabthani, N.-E. Jaidane, D. Mayou, and G. Trambly de Laissardière, Mobility gap and quantum transport in a functionalized graphene bilayer, *J. Phys.: Condens. Matter* **30**, 195701 (2018).

[32] J. Katoch, T. Zhu, D. Kochan, S. Singh, J. Fabian, and R. K. Kawakami, Transport Spectroscopy of Sublattice-Resolved Resonant Scattering in Hydrogen-Doped Bilayer Graphene, *Phys. Rev. Lett.* **121**, 136801 (2018).

[33] A. K. M. Pinto, N. F. Frazão, D. L. Azevedo, and F. Moraes, Evidence for flat zero-energy bands in bilayer graphene with a periodic defect lattice, *Physica E* **119**, 113987 (2020).

[34] J. Son, H. Ryu, J. Kwon, S. Huang, J. Yu, J. Xu, K. Watanabe, T. Taniguchi, E. Ji, S. Lee, Y. Shin, J. H. Kim, K. Kim, A. M. van der Zande, and G.-H. Lee, Tailoring single- and double-sided fluorination of bilayer graphene via substrate interactions, *Nano Lett.* **21**, 891 (2021).

[35] M. Moaied, J. V. Alvarez, and J. J. Palacios, Hydrogenation-induced ferromagnetism on graphite surfaces, *Phys. Rev. B* **90**, 115441 (2014).

[36] D. W. Boukhvalov, M. I. Katsnelson, and A. I. Lichtenstein, Hydrogen on graphene: Electronic structure, total energy, structural distortions and magnetism from first-principles calculations, *Phys. Rev. B* **77**, 035427 (2008).

[37] G. Trambly de Laissardière, J.-P. Julien, and D. Mayou, Quantum Transport of Slow Charge Carriers in Quasicrystals and Correlated Systems, *Phys. Rev. Lett.* **97**, 026601 (2006).

[38] G. Trambly de Laissardière, C. Oguey, and D. Mayou, Subdiffusive electronic states in octagonal tiling, *J. Phys.: Conf. Ser.* **809**, 012020 (2017).

[39] G. Trambly de Laissardière, O. F. Namarvar, D. Mayou, and L. Magaud, Electronic properties of asymmetrically doped twisted graphene bilayers, *Phys. Rev. B* **93**, 235135 (2016).

[40] G. Bouzerar and D. Mayou, Quantum transport in self-similar graphene carpets, *Phys. Rev. Res.* **2**, 033063 (2020).

[41] G. Bouzerar and D. Mayou, Quantum transport in flat bands and supermetallicity, *Phys. Rev. B* **103**, 075415 (2021).

[42] V. M. Pereira, J. M. B. Lopes dos Santos, and A. H. Castro Neto, Modeling disorder in graphene, *Phys. Rev. B* **77**, 115109 (2008).

[43] J. P. Robinson, H. Schomerus, L. Oroszlány, and V. I. Fal'ko, Adsorbate-Limited Conductivity of Graphene, *Phys. Rev. Lett.* **101**, 196803 (2008).

[44] T. O. Wehling, S. Yuan, A. I. Lichtenstein, A. K. Geim, and M. I. Katsnelson, Resonant Scattering by Realistic Impurities in Graphene, *Phys. Rev. Lett.* **105**, 056802 (2010).

[45] D. Mayou, Calculation of the conductivity in the short-mean-free-path regime, *Europhys. Lett.* **6**, 549 (1988).

[46] D. Mayou and S. N. Khanna, A real-space approach to electronic transport, *J. Phys. I France* **5**, 1199 (1995).

[47] S. Roche and D. Mayou, Conductivity of Quasiperiodic Systems: A Numerical Study, *Phys. Rev. Lett.* **79**, 2518 (1997).

[48] S. Roche and D. Mayou, Formalism for the computation of the RKKY interaction in aperiodic systems, *Phys. Rev. B* **60**, 322 (1999).

[49] F. Triozon, J. Vidal, R. Mosseri, and D. Mayou, Quantum dynamics in two- and three-dimensional quasiperiodic tilings, *Phys. Rev. B* **65**, 220202(R) (2002).

[50] S. Latil, S. Roche, D. Mayou, and J.-C. Charlier, Mesoscopic Transport in Chemically Doped Carbon Nanotubes, *Phys. Rev. Lett.* **92**, 256805 (2004).

[51] H. Ishii, S. Roche, N. Kobayashi, and K. Hirose, Inelastic Transport in Vibrating Disordered Carbon Nanotubes: Scattering Times and Temperature Dependent Decoherence Effects, *Phys. Rev. Lett.* **104**, 116801 (2010).

[52] G. Jemaï, J. J. Khabthani, G. Trambly de Laissardière, and D. Mayou, Quantum localization and electronic transport in covalently functionalized carbon nanotubes, *J. Phys.: Condens. Matter* **32**, 115301 (2019).

[53] Z. Fan, J. H. Garcia, A. W. Cummings, J. E. Barrios-Vargas, M. Panhans, A. Harju, F. Ortmann, and S. Roche, Linear scaling quantum transport methodologies, *Phys. Rep.* **903**, 1 (2021).

[54] D. G. Pettifor and D. L. Weaire, eds., *The Recursion Method and its Applications*, Springer Series in Solid-State Sciences 58 (Springer, Berlin, 1987).

[55] O. F. Namarvar, A. Missaoui, L. Magaud, D. Mayou, and G. Trambly de Laissardière, Electronic structure and quantum transport in twisted bilayer graphene with resonant scatterers, *Phys. Rev. B* **101**, 245407 (2020).

[56] A. Lacroix, G. Trambly de Laissardière, P. Quémerais, J.-P. Julien, and D. Mayou, Modeling of Electronic Mobilities in Halide Perovskites: Adiabatic Quantum Localization Scenario, *Phys. Rev. Lett.* **124**, 196601 (2020).

[57] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevB.104.245125 for numerical treatment of isolated midgap states, elastic mean free path results, DOS, and conductivity for a concentration $c = 0.5\%$ of $A_1B_1$-Va, and the conductivity by $B_1$-midgap states and $A_1$-midgap states.

[58] G. Otero, C. González, A. L. Pinardi, P. Merino, S. Gardonio, S. Lizzit, M. Blanco-Rey, K. Van de Ruit, C. F. J. Flipse, J. Méndez, P. L. de Andrés, and J. A. Martín-Gago, Ordered Vacancy Network Induced by the Growth of Epitaxial Graphene on Pt(111), *Phys. Rev. Lett.* **105**, 216102 (2010).

[59] S. Yuan, H. De Raedt, and M. I. Katsnelson, Electronic transport in disordered bilayer and trilayer graphene, *Phys. Rev. B* **82**, 235409 (2010).

[60] J. W. González, H. Santos, M. Pacheco, L. Chico, and L. Brey, Electronic transport through bilayer graphene flakes, *Phys. Rev. B* **81**, 195406 (2010).

[61] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, *Rev. Mod. Phys.* **81**, 109 (2009).

[62] P. A. Lee and T. V. Ramakrishnan, Disordered electronic systems, *Rev. Mod. Phys.* **57**, 287 (1985).

[63] G. Trambly de Laissardière and D. Mayou, Electronic transport in graphene: Quantum effects and role of local defects, *Mod. Phys. Lett. B* **25**, 1019 (2011).

[64] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, *Nature (London)* **556**, 80 (2018).

[65] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, *Nature (London)* **556**, 43 (2018).