Jorge Linares
email: author.jorge.linares@uvsq.fr

Catherine Cazelles

Pierre-Richard Dahoo

Devan Sohier

Thomas Dufaud

] ] Kamel Boukheddaden

Shape, size, pressure and matrix effects on 2D Spin Crossover nanomaterials studied using density of states obtained by dynamic programming

Keywords: Phase transitions, Spin crossover, nanomaterials, dynamic programming

In the present work, numerical simulations based on a new algorithm specific for 2D configurational topology of spin crossover nanoparticles embedded in a matrix are presented and discussed in the framework of the Ising-like model taking into account for short-(J) and longrange (G) interactions as for surface effects (L). The new algorithm is applied to calculate the density of states for each macro-state, which is then used to calculate exactly the thermal behavior of spin-crossover nanoparticles under an applied pressure. We find that the pressure plays the role of a conjugate parameter of the temperature. Thus, increasing pressure is somehow equivalent to reducing the temperature

Introduction

Spin-crossover (SCO) compounds which belong to the field of molecular magnetism are typical examples of a first-order phase transition with thermal hysteresis [1][START_REF] Gütlich | Chem. Spin Croosover in Transition Metal Compounds I-III[END_REF][START_REF] Gütlich | [END_REF][4][5][6][7][8][9] . Fe(II) based SCO complexes are among the most studied switchable molecular materials in which the spin transition takes place between two spin states, namely the high spin (HS) state, stable at high temperature, and the low-spin (LS) state, which is stable at lower temperature. It is quite well known that in Fe(II) d 6 spin-crossover complexes, the electronic structure of the HS state is with a spin S=2, whereas that of the LS state is with a spin S=0. Because of the antibonding character of the orbitals, the strength of the bond between the Fe(II) and the ligand is weaker (resp. stronger) in the HS (resp. LS) state with two (resp. zero) electrons in the orbital. Consequently, the bond length between the ligand and the Fe(II) metal centre is longer in the HS than that in the HS by about [START_REF] Gütlich | Chem. Spin Croosover in Transition Metal Compounds I-III[END_REF] .From the elastic rigidity point of view, the SCO materials are then soft and distortable in the HS state and rigid in the LS state, which has a higher bulk modulus than the HS state. It is important to mention for the non-specialist reader that the SCO solids are paramagnetic in the HS and do not show magnetic ordering, while they are diamagnetic in the LS state. The reason of the absence of the magnetic ordering can be attributed to their molecular structures, which leads to large distances ( nm) between the iron sites. The competition between the spin states and the elastic interactions between the molecular units is therefore responsible for the resulting "bi-stable" character of the spin-crossover solids. Indeed, due to the local volume change accompanying each spin transition of the molecules, long-and short-range range elastic interactions [10][11][12][13] take place between the SCO units causing their cooperative hysteretic first-order thermal transitions, driven by the large entropy (originating from spin and vibrations) of the HS state compared to that of the LS. In this respect, experimental and theoretical studies of solid state cooperative effects have become a very popular topic and as a result the investigation on the potential industrial applications offered by SCO compounds is worldwide.

One of the quest in many research groups concerns new materials for information storage and these compounds as bistable materials are the focus of research work. Indeed, spin transition materials have been studied for a long time for their thermally [START_REF] Gütlich | Chem. Spin Croosover in Transition Metal Compounds I-III[END_REF][START_REF] Gütlich | [END_REF] , pressure [14][15][16][17] , electric [5] , magnetic-field [18] , and light [19] induced bistabilities. They have potential applications as sensors of pressure and temperature [7] but also as displays and actuators [8] and recently they have shown abilities as materials for molecular spintronics [1] due to the interplay between the conducting properties and their spin states. Since ten years now, nanoparticles of spincrossover materials (SCO) have been synthesized and studied with the aim to understand how the reduction of the size influences their switching properties. In addition to the study of the size dependence of ensemble of free nanoparticles, other investigations related to nanoparticles embedded into polymeric matrices have also been developed in order to clarify the role of the environment or surface nanoparticles on their thermal behaviors [20][21][22] . Several experimental results demonstrated that the physical properties of the nanoparticles crucially depend on their size, shape and environment and a large number of models have been proposed to mimic the experimental facts.

In this work, we will focus solely on the modeling of thermal and pressure effects in 2D SCO lattices. From the theoretical side, the modelling of SCO phenomena helps to get new insights on the physical properties of SCO nanoparticles, particularly those displaying phase transitions, whose behavior as a function of their size and their interaction parameters play an important role in the control of the switching features of these complex materials, with the ultimate goal to help to design novel and specific systems with tailored properties.

From the general point of view, the modeling of phase transition allows a better understanding of the microscopic physical processes from which originates the cooperativity of the materials, as well as the nature of the interaction's mechanisms (magnetic, elastic, electrostatic, electronic …) involved during their transformation. In spin-crossover materials, the elastic interactions have been identified as playing a key role in their macroscopic thermo-, photo-and piezo-transformations between LS and the HS states. These elastic interactions are mainly due to the molecular volume change accompanying the spin transition at the microscopic scale, which deploys in the long-range way over the lattice through acoustic phonons. The long-range elastic interactions which stabilize homogeneous phases in most of the cases, compete with the ligand field which has a local and electronic nature, following which the thermal-dependence of the HS fraction curve may exhibit several shapes, among which, one can quote: (i) first order phase transitions with thermal hysteresis, (ii) incomplete spin transitions characterized by the presence of residual HS fractions at low-temperature, (iii) gradual and stepwise one with two or three steps, as well as re-entrant transitions. In the special case of SCO nanoparticles, these behaviors can be obtained by changing the size of the nanoparticles, which influences the balance of interactions between the bulk and surface contributions as well as other effects related to external stimuli (light, pressure etc.) The manuscript is organized as follows: in Sec. 2 the model and the principles of calculations are explained. Sec. 3 summarizes the numerical simulations' results on 2D lattices and a detailed discussion of the behavior of the high-spin fraction under various stimuli. Finally, in Sec. 4 we conclude and draw some possible extensions of this work.

Model and principle of calculation

Various physical techniques and methods, such as magnetic, optical spectroscopy, X-ray and neutron diffraction and heat capacity measurements are used to study the switching properties of SCO materials. The switching process from LS to HS triggered by an external stimuli, have been simulated using different types of models such as, the Ising-like model [23- 26] , the atom-phonon coupling [27] , mechano-elastic [13,28] , electro-elastic [29,30], vibronic [31,32], elastic [33].

Below, we first present the Ising-like model, with includes short-(J) and long-(G) range interactions together with the interaction (L) between surface molecules with the surrounding matrix. The density of states, needed to perform our numerical simulation is obtained with a new algorithm that will introduced and described in details.

Extended Ising Model with short-and long-range interactions

Forty years ago, Wajnflasz and Pick [23] proposed the use of the Ising model to describe the spin transition behavior in which only "short-range interaction" between the spincrossover sites were considered. Later, Bousseksou et al. [24] who made a thorough analysis of this model in the mean field approximation, reproduced, by introducing ferro-magnetic like and short-range "antiferromagnetic-like" interaction, the double-step character of some experimentally observed spin transition solids. Furthermore, Linares et al. [25] , by the introduction of a long range interaction in the Hamiltonian besides the short range interaction, were able to reproduce hysteresis in 1D compounds, as well.

With these contributions, the Hamiltonian now includes the short-and long-range interactions to which we added an energetic contribution "L" [11-34-35] , «to the ligand-field of edge molecules of a SCO nanoparticle ». This term allows us to explicitly take into account for the interactions between molecules on the edge and the environment in close contact, which weaken the molecular field.

In addition to these energetic terms, the present study considers also the contribution of an external isotropic pressure, P, whose effect is to renormalize the gap energy Δ (energy difference between HS and LS states), the effective expression of which is given by 𝛥 + 𝑃 × 𝛥𝑉 , where 𝑃 is the applied external pressure and 𝛥𝑉 = 𝑉 𝐻𝑆 -𝑉 𝐿𝑆 is the volume change of the material between the HS and the LS states.

The behaviour of SCO nanoparticles under the external constraints of pressure or surface effects as experimentally studied by one of the authors [15][16] is simulated using an Hamiltonian operator based on an extended Ising model. Simulation is performed by applying the algorithm described in section 2.2 to the Hamiltonian of the extended model described in a previous work by some of the present authors [11,36] and which is given in equation 1 for selfconsistency and clarity.

𝐻 = ∆ + 𝑃 × ∆𝑉 -𝑘 𝐵 𝑇 𝑙𝑛 𝑔 𝜎 𝑖 𝑁 𝑖=1 -𝐺 𝜎 𝑖 𝑁 𝑖=1 < 𝜎 > -𝐽 𝜎 𝑖 <𝑖,𝑗> 𝜎 𝑗 -𝐿 𝜎 𝑘 𝑀 𝑘=1 ( )
In (1), 𝜎 is a fictitious operator, which has two eigenvalues +1 and -1, respectively associated with the HS and LS states. 𝑇 is the temperature, 𝑔 is the degeneracy ratio between the HS and LS states and 𝑘 𝐵 is the Boltzmann constant, 𝜎 𝑘 represents the fictitious spin state of the molecules on the edge.

The way Hamiltonian (1) is expressed allows separation of the coupling terms pertaining to the interactions of the environmental matrix at the core and the edge as illustrated in figure 1, in terms of black dots and red dots respectively. It is then straight forward to discriminate the effects of these terms on the simulated behavior of the material. The first right hand side term in (1), containing the effective ligand field renormalized by the pressure effects, is a one site contribution which concerns all (bulk and surface) lattice atoms.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
There, 𝑔 = 𝑔 𝐻𝑆 𝑔 𝐿𝑆 represents the ratio of the degeneracies between the HS and LS states. The volume change accompanying the spin transitions of the molecules is accounted for by the second term of equation ( 1) which models the long range elastic interactions, which spreads throughout the lattice. This effect is taken into account here through the uniform and infinite long-range interactions, written as a mean-field term, in which < 𝜎 > is the average magnetization per site. The third contribution expresses short-range interactions between the SCO species, while the last term pertains to the coupling of surface molecules with their specific environment. Finally, parameters N and M describe the total number of molecules in the lattice and at the surface or edge respectively.

For self-consistency we recall the definition of the macroscopic parameters [36] m ("total magnetization"), s (pairwise lattice correlations) and c ("surface magnetization"), used to describe thermodynamic properties and which are implemented as such in the algorithm of section 2.2.

𝑚 = 𝜎 𝑖 𝑁 𝑖=1 = 𝑁 < 𝜎 > 𝑠 = 𝜎 𝑖 <𝑖,𝑗> 𝜎 𝑗 𝑐 = 𝜎′ 𝑘 𝑀 𝑘=1 (2)
Equation 1 can then be expressed in terms of m, s and c as in equation 3:

𝐻 = . ∆+∆𝑉.𝑃-𝑘 𝐵 𝑇 𝑙𝑛 𝑔 2 -𝐺 < 𝜎 >/ 𝑚 -𝐽𝑠 -𝐿𝑐 (3)
The number of configurations corresponding to the same m, s and c values leads to the calculation of the density of macrostates, 𝑑(𝑚, 𝑠, 𝑐) from which the thermodynamic properties of the system of nanoparticles can be determined. If the density of macrostates is denoted by 𝑑(𝑚 𝑖 𝑠 𝑖 𝑐 𝑖 ), then < 𝜎 >the average value of the operator 𝜎 can be calculated from equation 4 [31] .

< 𝜎 >= 𝑚 𝑖 𝑁 𝑁𝐿 𝑖=1 𝑑(𝑚 𝑖 𝑠 𝑖 𝑐 𝑖 )𝑒𝑥𝑝 -𝛽(-ℎ𝑚 𝑖 -𝐽𝑠 𝑖 -𝐿𝑐 𝑖 ) 𝑑(𝑚 𝑖 𝑠 𝑖 𝑐 𝑖 )𝑒𝑥𝑝 -𝛽(-ℎ𝑚 𝑖 -𝐽𝑠 𝑖 -𝐿𝑐 𝑖 ) 𝑁𝐿 𝑖=1 (4)
Where, 𝛽 = and NL is the number of possible *𝑚, 𝑠, 𝑐+ configurations. The field-like term, h, appearing in Eq. ( 4) is expressed as

ℎ = -. ∆+∆𝑉.𝑃-𝑘 𝐵 𝑇 𝑙𝑛 𝑔 2 -𝐺 < 𝜎 >/ (5)
The HS fraction, 𝑁 𝐻𝑆 , that is the probability to occupy the HS state is given by the equation,

𝑁 𝐻𝑆 = +<𝜎> (6)
𝑁 𝐻𝑆 is determined from equation 4 using values calculated for 𝑑(𝑚, 𝑠, 𝑐), and the different corresponding *𝑚, 𝑠, 𝑐+ configurations both by Monte Carlo Entropic Sampling [26] and by the algorithm described in section 2.2

Under the assumption that the transition temperature of the system is the result of a null total effective ligand-field, the transition temperature, 𝑇 𝑒𝑞 of a square lattice can be calculated from:

∆+∆𝑉.𝑃-𝑘 𝐵 𝑇 𝑒𝑞 𝑙𝑛 𝑔 2 × (𝑁 𝑥 -) 2 + ∆+∆𝑉.𝑃-2𝐿-𝑘 𝐵 𝑇 𝑒𝑞 𝑙𝑛 𝑔 2 × (𝑁 𝑥 -) = (7)
Where the number of core and edge sites are given by (𝑁 𝑥 -) and (𝑁 𝑥 -) respectively when the length of one side of the square is given by 𝑁 𝑥 .

As shown by Muraoka et al. [9] , the transition temperature 𝑇 𝑒𝑞 can be expressed in terms of 𝑇 𝑒𝑞 𝑠𝑢𝑟𝑓 and 𝑇 𝑒𝑞 𝑏𝑢𝑙𝑘 , the surface and bulk transition temperatures respectively as :

𝑇 𝑒𝑞 = 𝑁 𝑐 𝑁 𝑇 𝑒𝑞 𝑏𝑢𝑙𝑘 + 𝑁 𝑠 𝑁 𝑇 𝑒𝑞 𝑠𝑢𝑟𝑓 (8)
where

𝑇 𝑒𝑞 𝑏𝑢𝑙𝑘 = ∆+∆𝑉.𝑃 𝑘 𝐵 𝑙𝑛 𝑔 𝑎𝑛𝑑𝑇 𝑒𝑞 𝑠𝑢𝑟𝑓 = ∆+∆𝑉.𝑃-𝐿 𝑘 𝐵 𝑙𝑛 𝑔 (9)
To study the effect of the interaction between the nanoparticles at the surface and the matrix environment through the coupling term L, different cases have been considered. Keeping the temperature at a fixed value and varying the pressure (isothermal conditions) and keeping the pressure at a fixed value and varying the temperature (isobaric conditions) for L=0 (absence of surface effects) and 𝐿 ≠ (including surface effects).

Derivation of the density of states d[m,s,c] from dynamic programming algorithm

For a number of atoms 𝑛 , the computation of 𝑑(𝑚, 𝑠, 𝑐) by the naive enumeration algorithm is of complexity 𝑂( 𝑛 ) , which is considered intractable even at nanoscale. Stochastic approach with Monte Carlo Metropolis method or Monte Carlo Entropic sampling is normally used [26] .

Nevertheless, considering only 𝑛, it is possible to reduce the complexity and perform exact computation based on the Transfer Matrix method as it has been proposed to compute 𝑑(𝑠) in the ferro-magnetic field [37] .

We propose here a general technique based on a dynamic programming approach [38] that, in a sense, extends the Transfer Matrix method [39] and enables to compute 𝑑(𝑚, 𝑠, 𝑐) for any geometry. It leads to a parallel algorithm with complexity 𝑂 .𝑛 5 2 𝑛 / in time and space. This exact density of states has been used for the study of SCO compounds firmly expected for the design of a three-state electronic storage system [40] .

In this section, we propose a general algorithm for the exact computation of the density of states 𝑑(𝑚, 𝑠, 𝑐) which is denoted by 𝑧(𝑚, 𝑠, 𝑐) in the following and then explain the particular case of the 2D grid.

Problem statement

We consider a set 𝐴 of atoms, each of spin ± . Note 𝜎 𝑎 the spin of atom 𝑎 ∈ 𝐴. Pairs of atoms may be linked, note 𝐿 ⊂ 𝐴 × 𝐴 the set of pairs of adjacent atoms. The considered set of atoms has a border 𝐵 through which it interacts with the outer world.

To compute the Hamiltonian associated to this situation, one needs to compute the number 𝑧(𝑚, 𝑠, 𝑐) of configurations characterized by the three already defined quantities, which are recalled in the following :

-𝑚 = 𝑎∈𝐴 𝜎 𝑎 , is the sum of all spins; -𝑐 = 𝑎∈𝐵 𝜎 𝑎 , is the sum of the spins at the edge of the grid; -𝑠 = (𝑎,𝑏)∈𝐿 𝜎 𝑎 𝜎 𝑏 , is the sum of the product of spins of adjacent atoms. Thus, the knowledge of 𝑧 ′ is sufficient to compute 𝑧. From the equalities above, one sees that the parity of 𝑚, 𝑠 and 𝑐 only depends on that of |𝐴|, |𝐿| and |𝐵|. In other words, for a given particle topology, indexing the 𝑧 array with 𝑚, 𝑠 and 𝑐 imposes that seven elements in eight are zero; the 𝑧 ′ array is much denser and thus allows a better exploitation of the memory available. In the following of this section, we compute 𝑧 ′ and, for the sake of simplicity, drop the primes.

Let 𝜎 𝑎 ′ = +𝜎 𝑎 , which is 1 if 𝜎 𝑎 = ,

Dynamic programming algorithm

When attaching a new atom 𝑎 to the set 𝐴, 𝑚 is incremented if the new atom has spin 1, 𝑐 is incremented if the new atom is an edge atom and has spin 1, and 𝑠 increases by the number of atoms with a different spin to which the new atom is attached. Thus, computing the new (𝑚, 𝑠, 𝑐) requires to know its former value but also, in order to compute 𝑠, the spins of the atoms to which 𝑎 is attached.

Generalizing this observation, we can compute the array 𝑧 𝐴 (𝑚, 𝑠, 𝑐) for a set of atoms 𝐴 by adding the atoms in 𝐴 one after the other, and remembering the number of configurations with a particular value of (𝑚, 𝑠, 𝑐) and of each atom not all of the neighbours of which have yet been attached. Thus, note 𝑧 𝐴 ′ (𝑚, 𝑠, 𝑐, 𝜎) the number of configurations with given 𝑚, 𝑠, 𝑐 for a given set of atoms 𝐴 ′ , and ∀𝑎 ∈ 𝐴 ′′ , 𝜎(𝑎) , 𝐹 being the set of nodes in 𝐴 ′ that have neighbours in 𝐴 ∖ 𝐴 ′ . Thus, with 𝐹 the set atoms to which 𝑎 is attached, and to which no other nodes will be attached, if further atoms will be attached to 𝑎:

Configuration (𝑚, 𝑠, 𝑐, 𝜎) arises from any configuration with a valuation 𝜏 of atoms in 𝐹 in which the number of atoms with spin 1 was 1 less if the new atom is of spin 𝜎(𝑎) = , or the same if the atom is of spin 𝜎(𝑎) = .

The number of unbalanced edges is increased by the number of atoms (either in 𝐹, in which case 𝜎 is not defined but 𝜏 is, or not, in which case 𝜎 is); the number of edge atoms with spin 1 is less by one than that of the target configuration if 𝑎 is an edge atom ( 𝐵 (𝑎) is 1 if 𝑎 ∈ 𝐵 and 0 else) and the spin of 𝑎 is 1; last, the valuation function is the same than that of 𝜎, except for 𝑎 for which it is undefined, and it is defined on 𝐹 by 𝜏.

If no further atoms will be attached to 𝑎:

The configurations in the first sum are those to which attaching 𝑎 with spin -1 (𝜎(𝑎) = ) leads to the target configuration (𝑚, 𝑠, 𝑐, 𝜎); the second sum lists the configurations leading to the target configurations if 𝑎 is of spin 1.

The initialization of this algorithm is 𝑧 ∅ ( , , , ⊥) = , 𝑧 ∅ is 0 for all others configurations.

Complexity

The resulting algorithm consists in a loop on the atom 𝑎 ∈ 𝐴, the body of which consists in updating all possible values for (𝑚, 𝑠, 𝑐, 𝜎), the update consisting in summing over all partial valuations 𝜏. The outer loop has thus complexity |𝐴|. 𝑚 can take values from 0 to |𝐴|, 𝑠 from 0 to |𝐿|, and 𝑐 from 0 to |𝐵|. 𝜎 and 𝜏, together, build a partial function on 𝐴, and thus can take at most 𝑙 different values, with 𝑙 the maximal number of atoms in 𝐹, which depend on the topology of the system and the order with which atoms are added in 𝐴. The overall complexity is thus: 𝛩(|𝐴| 2 |𝐵||𝐿| 𝑙 ). The memory occupation is defined by the size of the array 𝑧: 𝛩(|𝐴||𝐵||𝐿| 𝑙 ).

The 2D cartesian grid case

Algorithm

Let us consider a nanoparticle composed of 𝑛 × 𝑛 atoms on a 2D Cartesian grid as illustrated in the next figure.

In the following we consider both the geometry of the particle and the histogram which stores the densities of states 𝑧(𝑚, 𝑠, 𝑐). In the particular case of a 𝑛 × 𝑛 grid, the operations of the above algorithm, can be scheduled by adding the atoms according to the topology of the grid, row by row. Then, at step (𝑖, 𝑗), 𝑖 rows are complete, and the 𝑖 + st row contains 𝑗 atoms. The atoms to which further atoms will be attached are the 𝑛 atoms constituted by the 𝑗 atoms of the row currently built, and the 𝑛 -𝑗 last atoms of the previous row. Adding a new atom discards the atom below, to which no new node will be attached; this creates two new edges (except if 𝑗 = or 𝑖 = ) that can be balanced or not. 𝜎 can thus be coded as a vector of 𝑛 binary values (at most: for the first row, 𝜎 is incomplete).

With |𝐴| = 𝑛 1 × 𝑛 2 , |𝐵| = (𝑛 1 + 𝑛 2 -) , |𝐿| = 𝑛 1 (𝑛 2 -) + 𝑛 2 (𝑛 1 -) , and 𝑙 = 𝑛 2 , the complexity of this algorithm is 𝛩((𝑛 1 𝑛 2 ) 3 (𝑛 1 + 𝑛 2 ) 𝑛 2 ) = 𝑂(𝑛 compute_2dCart_msc(z,z,i,j,m,s,c). Those three functions are called in the inner loop c.

At each step it gets the data reference in histogram by calling index_2D(𝑚, 𝑠, 𝑐). It updates one line of index K1D, of the array 𝑧 considering the 8 dependencies. Indices of the dependencies are stored in the array K1D_deps and correspond to the triplet (𝑚, 𝑠 -, 𝑐), (𝑚, 𝑠 -, 𝑐), (𝑚 -, 𝑠, 𝑐), (𝑚 -, 𝑠 -, 𝑐), (𝑚 -, 𝑠 -, 𝑐), (𝑚 -, 𝑠 -, 𝑐 -).

Complexity and experimental platform

In the t-loop of compute_2dCart_msc , the update consists in multiple reads in 𝑍,𝐾 𝐷_𝑑𝑒𝑝𝑠array, and 𝑛 independent write in 𝑍,𝐾 𝐷-. We can exploit the later for a parallel loop implementation on a shared memory architecture.

The dynamic programming in the 2D case decreases dramatically the time complexity, but multiply the memory complexity by 𝑛 according to the study above. Nanoparticles of our interests corresponds to a range from 𝑛 = × to 𝑛 = × . Let us assume that the frequency of a processor unit is about 2.6 GHz and a memory on a node is about 32 GB, we can illustrate this complexity given order of human perception.

The table below shows that, under those assumptions, we can now consider nanoparticle simulation from hours to days, and few weeks in the largest case with our approach.

The figure 3 shows the memory complexity function of the size of the problem for the dynamic programming approach. The complexity has increased and exhibit the needs of distributed memory architecture.

In this article simulation for grids of atoms up to × where performed on a single shared memory node with 16 cores.

Time complexity 𝑶 .𝒏

Numerical results and discussion

In a first stage, we look for the size effects on the thermal and the pressure behaviour of the HS fraction of a SCO nanoparticle. Next, we study the shape effects on the pressure behaviour of the thermal-dependence of the HS fraction. In each stage, we consider the two following alternatives: calculations with ligand-field 𝐿 and calculations without ligand-field.

In the present study, we will analyze the size and shape effects on the thermal and pressure dependence of the HS fraction for two typical spin-crossovers compounds characterized by a set of thermodynamic parameters (enthalpy and entropy changes). The first one (system 1) has /k B =1300K, ln(g) =6.01 leading to the transition temperature Teq=216 K and the second one (system 2) is characterized by /k B =1978,6 K, ln(g) =6.9 (transition temperature Teq=286 K).

2D size effects -The case L≠0

Figure 4 gives the temperature effect on the HS fraction for a given applied isotropic pressure when the interaction at the boundary of the SCO system with the surrounding medium cannot be neglected (𝐿 ≠ ). The analysis of the results obtained from 3×3, 4×4, 5×5 and 12×12 square lattices shows the presence of a thermal hysteresis with the smaller sizes, 3×3 and 4×4 squares compared to the larger 5×5, 12×12 which exhibit a gradual spin transition. In fact, when edge effects increase as a result of a higher perimeter to surface ratio, which is the case with the smaller squares compared to the larger ones, the transition temperature shifts to the reference edge temperature, 𝑇 𝑒𝑞 𝑠𝑢𝑟𝑓 = (𝛥 -𝐿)/𝑘 𝐵 𝑙𝑛 𝑔. Because hysteresis phenomenon is the fingerprint of first order transition, this result means that 𝑇 𝑒𝑞 must be lower than the critical temperature, 𝑇 𝐶 , at which order-disorder transition of the corresponding Ising model occurs.

Parameter values are chosen from typical experimental data of standard SCO compounds and are fixed at: /k B =1300K, G/k B =172.2K, J/k B = 13.5 K, L/k B =130 K and ln(g) =6.01.

Using these values, we obtain for the curie temperature, T C = (2.269J +G)/k B =203 K, when short range and long range interactions J and G are present and for the surface and bulk transition temperatures, 𝑇 𝑒𝑞 𝑠𝑢𝑟𝑓 = 7 K and 𝑇 𝑒𝑞 𝑏𝑢𝑙𝑘 = K, respectively. Compared to the curie temperature, the surface and bulk transition temperatures fulfill the following inequalities, 𝑇 𝑒𝑞 𝑠𝑢𝑟𝑓 < 𝑇 𝐶 , whereas 𝑇 𝑒𝑞 𝑏𝑢𝑙𝑘 > 𝑇 𝐶 . A gradual transition is obtained with the latter relation, while the former implies a thermal hysteresis for surface lattice. It means that decreasing the size of the nanoparticles favors the surface effects for which the conditions of appearance of a first-order transitions are fulfilled. Concerning the effect of the pressure on the thermal properties of the SCO nanoparticules, we can see from Eq. ( 8) that the hydrostatic pressure appears as a conjugate parameter of temperature [41] . Increasing pressure stabilizes the LS state and it is equivalent to reducing the temperature.

Figure 5 summarizes the pressure dependence of the HS fraction under isothermal conditions at T=300 K, for the several lattice sizes. It is noted, as expected, that strong (resp. weak) pressure values stabilize the LS (resp. HS) state. As a result, when hydrostatic pressure is increased the molecules changes from the HS to the LS state. Moreover, and in contrast with the thermal behaviors of Fig. 4, the increase of the system's size leads to the increase of the width of the pressure hysteresis loop (12x12). It is observed that this isothermal transition depends on temperature, whose variation affects not only the width of the pressure hysteresis but also the character (first-order or gradual) of the transition.

-The case L=0

We now compare the thermal and pressure-dependences of the HS fraction when the interactions between the edges molecules and their surrounding environment are not taken into account. First, we start with the case of thermal hysteresis of Fig. 4, for which we fix 𝐿 = . Hence, the effective ligand field becomes the same for surface and bulk atoms leading to a unique surface and bulk transition temperature, 𝑇 = K, as clearly depicted in Fig. 6. The thermal behavior for the case L=0 is almost identical in the three numerical simulations: when the size increases the hysteresis width is increased and the equilibrium temperature is the same for the three cases. Indeed, according to Eqs. ( 7) and ( 8) the equilibrium temperature remains independent of the lattice size then L=0 (the system does not experience any surface effects). In contrast the order-disorder transition T OD increases with lattice size and then the thermal hysteresis becomes wider (see Fig 3) when the lattice size increases.

When the surface effects are absent (L=0), when one applies an isotropic pressure, under isothermal conditions at T=300K, the results obtained of the pressure dependence of the HS fraction (see Fig. 7) are significantly different from those of Fig. 5. If one considers the lattices 5x5 and 12x12, we clearly see that their hysteresis widths are much wider than those of corresponding lattices of Fig. 5. For example, the lattice 12x12, exhibits in Fig. 7 a pressure hysteresis width given by ΔP ~ 13 MPa compared to 7 MPa for Fig. 5 and a pressure transition P 1/2 bulk = (k B Tln(g)-Δ)/ΔV ~ 52 MPa, common to all other curves, which is smaller than that of P 1/2 =90 MPa case, which is clearly affected by the surface effects.

2D Shape effects under pressure

In a previous study, Guerroudj et al [35] analyzed the thermal-dependence of the HS fraction for different lattice shapes of SCO in nano configuration topology for a fixed total number of atoms ( 144). Here, we analyze the pressure-dependence of the HS fraction for different lattice shapes comprising 144 SCO molecules. We start with a 12x12 square lattice and through successive transformations from a square to a rectangle, one length is increased at the expense of the width sweeping all the different possible shapes until an elongated rectangular 2x72 lattice is obtained. For a × lattice, we define the aspect ratio parameter, = ( + -)/

, as the ratio between surface and volume numbers of molecules. Thus, for a 12x12 square-shape, we calculate = . , while for rectangular-shaped lattices of size 9x16, 8x18, 6x24, 4x36, 3x48 and 2x72, we found the respective values, 0.32, 0.33, 0.39, 0.53, 0.68 and 1. Figure 8 shows two different shapes for a system of 144 molecules.

-The case L≠0

As shown in Fig. 9a, the 2x72 (t=1) and 3x48 (t=0.68) lattices exhibit a gradual thermal spin transition. A first-order transition with a hysteresis loop appears when the "t" value decreases. We note that the width of the loop increases when t decreases and seems to saturate for the lattices 8x18, 9x16 and 12x12. In addition, Fig. 9b reporting the behavior under pressure indicates that the transition pressure, defined as (P up +P down )/2 which is approximately the center of the hysteresis, follows an increasing linear plot (see Fig. 10) with the ratio "t". -The case L=0

Figure 11 displays the pressure-dependence of the HS fraction for several lattice shapes in isothermal conditions (T=300 K). We observe that, except for the 2x72 (t=1) lattice which exhibits a gradual spin transition between the HS and the LS state, all other systems present a hysteresis loop, typical of a first-order transition. The width of this hysteresis loop increases when t decreases. This behavior is in line with that of Fig. 9b, obtained for L=0. Moreover, one can also notice that the pressure transition, 𝑃 , in Fig. 10 remains invariant (as in Fig. 10) whatever the shape of the material.

Summary and conclusion

In this paper, the computation of the Hamiltonian is based on the exact density of states, obtained by the algorithm presented in Section 2.2, and we analyze the behavior of 2D spincrossover compounds under temperature and pressure. We focus on the effect of the interactions between nanoparticles at the edges and the matrix in their specific local environment, modelled with the coupling parameter, 𝐿. When 𝐿 = , the temperature and the pressure play as conjugate parameters, and both of them lead to thermal hysteresis as the size of the nanoparticle is increased. In contrast, when 𝐿 ≠ , opposite effects for the pressure and temperature are observed on the thermal hysteresis as function of the nanoparticle size. The main effect underlying these differences in the simulated behaviours is solely related to the relationship between the equilibrium temperature and the order-disorder temperature when 𝑔 = (unit degeneracy ratio) and 𝛥 = and 𝐿 =0 (absence of ligand field), i.e., mainly a temperature dependence. Indeed, the pressure effect is independent of L at all sizes except in the limit of very small lattices. This property is due to a concomitance of the unequal effect of pressure in the bulk and at the surface because of a weaker ligand field at the surface compared to that in the bulk and the surface/volume ratio effect which is parametrized by 𝑡.

On the other hand, one can easily imagine that the effective short-and long-range interactions,J and 𝐺 , which have an elastic origin, may depend on the applied pressure through the change of the lattice parameter and the bulk modulus under pressure. In addition, since the elastic and phonon properties of the surface and the bulk are different, additional degrees of freedom may also emerge. For example, the effect of an expansion of the interaction parameters in terms of a reduced pressure parameter relative to a maximum pressure or in terms of a pressure law may be investigated in this respect for a deeper understanding of the SCO properties under various stimuli. This will benefit to the necessary reliability of new technological applications based on the switchable properties of these materials. Nevertheless, the results reported in this work are consistent with available experimental data as far as the effect of the internal interaction strength, the local environment at the surface or the system's size on the thermal and pressure induced transitions are concerned.

Figure 1 :

 1 Figure 1: Schematic view of the SCO nanoparticle. Black dots are bulk sites, while red dots are surface atoms that have specific interaction with their immediate environment.

  and 0 else. Then, let -𝑚 ′ = 𝑎∈𝐴 𝜎 𝑎 ′ = |𝐴|+𝑚 , the number of atoms with spin 1; -𝑐 ′ = 𝑎∈𝐵 𝜎 𝑎 ′ = |𝐵|+𝑐 , the number of edge atoms with spin 1; -𝑠 ′ = (𝑎,𝑏)∈𝐿 (𝜎 𝑎 ′ -𝜎 𝑏 ′ ) = |𝐿|+𝑠 (because 𝜎 𝑎 𝜎 𝑏 = -(𝜎 𝑎 -𝜎 𝑏 ) = -(𝜎 𝑎 ′ -𝜎 𝑏 ′ ) ), the number of « unbalanced » edges, the extremities of which have different spins;-𝑧 ′ (𝑚 ′ , 𝑠 ′ , 𝑐 ′ ) the number of configurations with given 𝑚 ′ , 𝑠 ′ and 𝑐 ′ .From these definitions, we have the following relations: 𝑚 = 𝑚 ′ -|𝐴|, 𝑐 = 𝑐 ′ -|𝐵| and 𝑠 = 𝑠 ′ -|𝐿|, and 𝑧(𝑚, 𝑠, 𝑐) = 𝑧 ′ .

  𝑧 𝐴 ′ ∪*𝑎+ (𝑚, 𝑠, 𝑐, 𝜎) = 𝜏:𝐹↦*0,1+ 𝑧 𝐴 ′ (𝑚 -𝜎(𝑎), 𝑠 -𝑏∈𝑁 𝑎 ∩𝐴 ′ (𝜎(𝑎) -𝜎(𝑏)) 2 ⬚ -𝑏∈𝑁 𝑎 ∩𝐹 (𝜎(𝑎) -𝜏(𝑏)) 2 , 𝑐 -𝐵 (𝑎) × 𝜎(𝑎), (𝜎 |𝐴 ′ ) ∨ 𝜏)

⬚

  𝑧 𝐴 ′ ∪*𝑎+ (𝑚, 𝑠, 𝑐, 𝜎) = 𝜏:𝐹↦*0,1+ 𝑧 𝐴 ′ (𝑚, 𝑠, 𝑐, 𝜎 ∨ 𝜏) ), 𝑐 -𝑎∈𝐵 , 𝜎 ∨ 𝜏

Figure 2 .

 2 Figure 2. (a) schematic 2D Cartesian particle grid. (b) Histogram for a 6x6 lattice.

Figure 3 .

 3 Figure 3. Memory complexity of the computation of the density of state for a 2D network; for n=6x6 to n=20x20

Figure 4 .

 4 Figure 4. Thermal evolution of the HS fraction for different lattice sizes: 3x3 (green left triangle), 4x4 (blue up triangle), 5x5 (black squares), 12x12 (red stars). The computational parameters are: Δ/k B = 1300 K, G/k B = 172.2 K, J/k B = 13.5 K, L/k B = 130 K, ln(g) = 6.01.

Figure 5 .

 5 Figure 5. Pressure-dependence of the HS fraction in isothermal conditions (T=300 K) for different lattice sizes: 3x3 (green left triangle), 4x4 (blue up triangle), 5x5 (black squares), 12x12 (red stars). The computational parameters are: Δ/k B = 1978.6 K, G/k B = 172.2 K, J/k B = 48 K, L/k B = 130 K, ln(g) = 6.906, ΔV = 100 Å 3 .

Figure 6 .

 6 Figure 6. Thermal evolution of the HS fraction for different lattice sizes: 4x4 (blue up triangle), 5x5 (black squares), 12x12 (red stars). The computational parameters are: Δ/k B = 1300 K, G/k B = 172.2 K, J/k B = 16.8 K, L/k B = 0 K, ln(g) = 6.01.

Figure 7 .

 7 Figure 7. Pressure-dependence of the HS fraction for T=300 K. for different lattice sizes: 4x4 (blue up triangle), 5x5 (black squares), 12x12 (red stars). The computational parameters are: Δ/k B = 1978.6 K, G/ k B = 172.2 K, J/ k B = 48, L/k B = 0 K and ln(g) = 6.906, ΔV = 100 Å 3 .

Figure 8 .

 8 Figure 8. (a) Square lattice with size 12x12 with a number core atoms 𝑁 𝐶 = and atoms at the edge, 𝑁 𝑒 = . (b) Rectangle-shaped lattice with size 4x36, containing 𝑁 𝐶 = 8 core atoms and 𝑁 𝑒 = 7 edge sites. Black and red dots represent respectively the core (N c ) and the edge (N e ) lattice sites.

Figure 9 .

 9 Figure 9. (a) Evolution of the high spin versus the pressure for several shapes : 2x72 (black square), 3x48 (red circle), 4x36 (green up triangle), 6x24 (blue down triangle), 8x18 (purple left triangle), 9x16 (dark yellow right triangle), 12x12 (orange star); the computational parameters are : Δ/k B = 1978.6 K, G/ k B = 172.2 K, J/k B = 49 K, L/k B = 120 K, T=300 K, ln(g)=6.906, ΔV=100 Å 3 ; (b) Zoom around the hysteresis-loop for the system shapes 8x18 (purple left triangle), 9x16 (dark yellow right triangle) and 12x12 (orange star).

Figure 10 .

 10 Figure 10. The pressure transition 𝑃 versus the surface/volume parameter t.

Figure 11 .

 11 Figure 11. Evolution of the high spin versus the pressure for several shapes in isothermal conditions (T=300 K) for, 2x72 (black square), 4x36 (green up triangle), 8x18 (purple left triangle), 12x12 (orange star) lattices. The computational parameters are : Δ/k B = 1978.6 K, G/ k B = 172.2 K, J/k B = 49 K, L/k B = 120 K, ln(g)=6.906, ΔV = 100 Å 3 .

  𝑛 2 and choosing 𝑛 as the smallest dimension of the grid. The algorithm below consists of loops over the geometry (𝑖, 𝑗) and histogram (𝑚, 𝑠, 𝑐)

							7 2	𝑛 ), taking
	𝑛 = 𝑛 1 which are respectively illustrated in the two figures 2a and 2b.
	We	create	three	functions	index_2D	(	m, s, c	),
	get_global_dependencies_2dCart_msc(k1D_deps,	m,	s,	c),		and

  𝑘 𝐷-, 𝑧,𝑘 𝐷_𝑑𝑒𝑝𝑠-, 𝑖, 𝑗, 𝑚, 𝑠, 𝑐). It consists in a loop over 𝜎.

	This	array	is	updated	at	each	step	by	calling	the	function
	get_global_dependencies_2dCart_msc(𝑘 𝐷_𝑑𝑒𝑝𝑠, 𝑚, 𝑠, 𝑐).
	The	update		of	the		intermediate	state	is	performed	by
	compute_2dCart_msc(𝑧,Algorithm: Dynamic	programming	for	2D	cartesian	grid
	For = : 1 For 𝑗 = : 𝑛 2						{𝑛 < 𝑛 }
		For 𝑚 = 𝑚𝑎𝑥𝑀: -:				{(𝑖 + ) * 𝑛 + 𝑗 + }
			For 𝑠 = 𝑠𝑚𝑎𝑥 EndFor
			EndFor
	EndFor
	EndFor
	EndFor

𝑚 : -: For 𝑐 = : 𝑐 <= 𝑚.and.𝑐 < 𝑛 + * (𝑖 + ) 𝑘 𝐷 = 𝑖𝑛𝑑𝑒𝑥_ 𝐷(𝑚, 𝑠, 𝑐, 𝑛, 𝑠𝑖𝑧𝑒𝐶, 𝑠𝑚𝑎𝑥, 𝑠𝑖𝑧𝑒 𝐷) 𝑔𝑒𝑡_𝑔𝑙𝑜𝑏𝑎𝑙_𝑑𝑒𝑝𝑒𝑛𝑑𝑒𝑛𝑐𝑖𝑒𝑠_ 𝑑𝐶𝑎𝑟𝑡_𝑚𝑠𝑐(𝑘 𝐷_𝑑𝑒𝑝𝑠, 𝑚, 𝑠, 𝑐) 𝑐𝑜𝑚𝑝𝑢𝑡𝑒_ 𝑑𝐶𝑎𝑟𝑡_𝑚𝑠𝑐(𝑧,𝑘 𝐷-, 𝑧,𝑘 𝐷_𝑑𝑒𝑝𝑠-, 𝑖, 𝑗, 𝑚, 𝑠, 𝑐)

Acknowledgments

The authors are grateful to Sallah E. Allal and Camille Harlé who did the preliminary simulations during their master internship. CHAIR Materials Simulation and Engineering of UVSQ-UPSAY, the French "Ministère de la Recherche", the Université de Versailles St. Quentin-en-Yvelines, Université Paris-Saclay, CNRS and ANR BISTA-MAT (ANR-12-BS07-0030-01) are warmly acknowledged for the financial support.


Omid Faizy Namarvar

Ahmed Missaoui

Laurence Magaud
email: laurence.magaud@grenoble.cnrs.fr

Didier Mayou
email: didier.mayou@grenoble.cnrs.fr§guy.trambly@u-cergy.fr

Guy Trambly De Laissardière

Electronic structure and quantum transport in twisted bilayer graphene with resonant scatters

Staking layered materials revealed to be a very powerful method to tailor their electronic properties. It has indeed been theoretically and experimentally shown that twisted bilayers of graphene (tBLG) with a rotation angle θ , forming Moiré pattern, confine electrons in a tunable way as a function of θ . Here, we study electronic structure and transport in tBLG by using tight-binding numerical calculations in commensurate twisted bilayer structures and a pertubative continuous theory, which is valid for not too small angles (θ >∼ 2 • ). This two approaches allow to understand the effect of θ on the local density of states, the electron lifetime due to disorder, the dc-conducitivity and the quantum correction of the conductivity due to multiple scattering effects. We distinguished the cases where disorder is equaly distributed in the two layer or only in one layer. When only one layer is disordered, diffusion properties depends strongly on θ , showing thus the effect of Moiré electronic localisation at intermediate angles, ∼ 2 • < θ <∼ 20 • .

I. INTRODUCTION

Staking layered materials is a very powerful method to tailor their electronic properties. [START_REF] Geim | Van der Waals heterostructures[END_REF] The properties not only depend on the choice of materials to be stacked but also on the details of the relative arrangement of the layers. It has thus been theoretically [START_REF] Lopes Dos Santos | Graphene bilayer with a twist: Electronic structure[END_REF][START_REF] Trambly De Laissardière | Localization of dirac electrons in rotated graphene bilayers[END_REF][START_REF] Suárez Morell | Flat bands in slightly twisted bilayer graphene: Tight-binding calculations[END_REF][START_REF] Bistritzer | Transport between twisted graphene layers[END_REF][START_REF] Bistritzer | Moiré bands in twisted double-layer graphene[END_REF][START_REF] Trambly De Laissardière | Numerical studies of confined states in rotated bilayers of graphene[END_REF] and experimentally [START_REF] Li | Observation of van hove singularities in twisted graphene layers[END_REF][START_REF] Luican | Single-layer behavior and its breakdown in twisted graphene layers[END_REF][START_REF] Brihuega | Unraveling the intrinsic and ro-bust nature of van hove singularities in twisted bilayer graphene by scanning tunneling microscopy and theoretical analysis[END_REF][START_REF] Huder | Electronic spectrum of twisted graphene layers under heterostrain[END_REF] shown that twisted bilayers of graphene, forming Moiré pattern, confine conduction electrons in a tunable way as a function of the angle of rotation of one layer with respect to the other. Recently, it has been experimentally proven that this electronic localization by geometry can induce strong electronic correlations [START_REF] Cao | Correlated insulator behaviour at half-filling in magic-angle graphene superlattices[END_REF] and a superconducting state [START_REF] Cao | Unconventional superconductivity in magic-angle graphene superlattices[END_REF] for certain angles called magic angles. [START_REF] Bistritzer | Moiré bands in twisted double-layer graphene[END_REF] Despite numerous studies of the electronic structure of these systems, [START_REF] Lopes Dos Santos | Graphene bilayer with a twist: Electronic structure[END_REF][START_REF] Trambly De Laissardière | Localization of dirac electrons in rotated graphene bilayers[END_REF][START_REF] Suárez Morell | Flat bands in slightly twisted bilayer graphene: Tight-binding calculations[END_REF][START_REF] Bistritzer | Transport between twisted graphene layers[END_REF][START_REF] Bistritzer | Moiré bands in twisted double-layer graphene[END_REF][START_REF] Trambly De Laissardière | Numerical studies of confined states in rotated bilayers of graphene[END_REF][START_REF] Li | Observation of van hove singularities in twisted graphene layers[END_REF][START_REF] Luican | Single-layer behavior and its breakdown in twisted graphene layers[END_REF][START_REF] Brihuega | Unraveling the intrinsic and ro-bust nature of van hove singularities in twisted bilayer graphene by scanning tunneling microscopy and theoretical analysis[END_REF][START_REF] Huder | Electronic spectrum of twisted graphene layers under heterostrain[END_REF][START_REF] Latil | Massless fermions in multilayer graphitic systems with misoriented layers: Ab initio calculations and experimental fingerprints[END_REF][START_REF] Lopes Dos Santos | Continuum model of the twisted graphene bilayer[END_REF][START_REF] Trambly De Laissardière | Electronic properties of asymmetrically doped twisted graphene bilayers[END_REF][START_REF] Chari | Resistivity of rotated graphite-graphene contacts[END_REF][START_REF] Le | Electronic structure and optical properties of twisted bilayer graphene calculated via time evolution of states in real space[END_REF][START_REF] Chung | Transport measurements in twisted bilayer graphene: Electron-phonon coupling and landau level crossing[END_REF][START_REF] Wu | Boosting carrier mobility of synthetic few layer graphene on sio2 by interlayer rotation and decoupling[END_REF][START_REF] Ribeiro-Palau | Twistable electronics with dynamically rotatable heterostructures[END_REF][START_REF] Andelković | Dc conductivity of twisted bilayer graphene: Angle-dependent transport properties and effects of disorder[END_REF][START_REF] Woo Jeon | Experimental evidence for interlayer decoupling distance of twisted bilayer graphene[END_REF][START_REF] Rickhaus | Transport through a network of topological channels in twisted bilayer graphene[END_REF] the consequences of the electronic localization by a Moiré on electrical transport properties are still very poorly known.

Graphene can be formed in multilayers on SiC [START_REF] Ohta | Controlling the electronic structure of bilayer graphene[END_REF][START_REF] Coletti | Charge neutrality and band-gap tuning of epitaxial graphene on sic by molecular doping[END_REF][START_REF] Brihuega | Quasiparticle chirality in epitaxial graphene probed at the nanometer scale[END_REF][START_REF] Hass | Highly ordered graphene for two dimensional electronics[END_REF][START_REF] Hass | The growth and morphology of epitaxial multilayer graphene[END_REF][START_REF] Emtsev | Interaction, growth, and ordering of epitaxial graphene on sic0001 surfaces: A comparative photoelectron spectroscopy study[END_REF][START_REF] Hass | Why multilayer graphene on 4h-SiC(0001) behaves like a single sheet of graphene[END_REF][START_REF] Sprinkle | First direct observation of a nearly ideal graphene band structure[END_REF][START_REF] Hicks | Symmetry breaking in commensurate graphene rotational stacking: Comparison of theory and experiment[END_REF] but also on metal surfaces such as Ni [START_REF] Luican | Single-layer behavior and its breakdown in twisted graphene layers[END_REF] and in exfoliated flakes, [START_REF] Li | Observation of van hove singularities in twisted graphene layers[END_REF] where interactions between successive layers play a crucial role. While on the Si face of SiC, multilayers have an AB Bernal stacking and do not show graphene properties, [START_REF] Ohta | Controlling the electronic structure of bilayer graphene[END_REF][START_REF] Coletti | Charge neutrality and band-gap tuning of epitaxial graphene on sic by molecular doping[END_REF][START_REF] Brihuega | Quasiparticle chirality in epitaxial graphene probed at the nanometer scale[END_REF][START_REF] Latil | Charge carriers in few-layer graphene films[END_REF][START_REF] Varchon | Ripples in epitaxial graphene on the si-terminated sic(0001) surface[END_REF][START_REF] Zhang | Band structure of abc-stacked graphene trilayers[END_REF][START_REF] Mccann | The electronic properties of bilayer graphene[END_REF][START_REF] Rozhkov | Electronic properties of graphene-based bilayer systems[END_REF] on the C-face multilayers are twisted multilayers of graphene with various angles of rotation between two layers. For large twisted angle θ between two layers, multilayers present graphene like properties even when they involve a large number of graphene planes. Indeed, as shown by ARPES, [START_REF] Emtsev | Interaction, growth, and ordering of epitaxial graphene on sic0001 surfaces: A comparative photoelectron spectroscopy study[END_REF][START_REF] Hass | Why multilayer graphene on 4h-SiC(0001) behaves like a single sheet of graphene[END_REF][START_REF] Sprinkle | First direct observation of a nearly ideal graphene band structure[END_REF][START_REF] Hicks | Symmetry breaking in commensurate graphene rotational stacking: Comparison of theory and experiment[END_REF] STM, [START_REF] Miller | Observing the quantization of zero mass carriers in graphene[END_REF] transport [START_REF] Berger | Electronic confinement and coherence in patterned epitaxial graphene[END_REF] and optical transitions, [START_REF] Sadowski | Landau level spectroscopy of ultrathin graphite layers[END_REF] they show properties characteristic of a linear graphene like dispersion. Therefore, in twisted bilayer graphene (tBLG) interlayer interaction does not systematically destroy graphene like properties, but it can lead to the emergence of very peculiar and new behaviors induced by the Moiré patterns that is accentuated for θ smaller than ∼ 20 • . Theoretical studies have predicted [START_REF] Lopes Dos Santos | Graphene bilayer with a twist: Electronic structure[END_REF][START_REF] Trambly De Laissardière | Localization of dirac electrons in rotated graphene bilayers[END_REF][START_REF] Suárez Morell | Flat bands in slightly twisted bilayer graphene: Tight-binding calculations[END_REF][START_REF] Bistritzer | Transport between twisted graphene layers[END_REF][START_REF] Bistritzer | Moiré bands in twisted double-layer graphene[END_REF][START_REF] Trambly De Laissardière | Numerical studies of confined states in rotated bilayers of graphene[END_REF][START_REF] Lopes Dos Santos | Continuum model of the twisted graphene bilayer[END_REF] the existence of three domains: (1) for large rotation angles (θ > 20 • ) the layers are decoupled and behave as a collection of isolated graphene layers. (2) For intermediate angles ∼ 2 • < θ < 20 • the dispersion, around Fermi energy E F , remains linear but the velocity is renormalized. Consequently, the energies of the two van Hove singularities E - and E + are shifted to Dirac energy E D when θ decreases, as it has been shown experimentally. [START_REF] Luican | Single-layer behavior and its breakdown in twisted graphene layers[END_REF][START_REF] Brihuega | Unraveling the intrinsic and ro-bust nature of van hove singularities in twisted bilayer graphene by scanning tunneling microscopy and theoretical analysis[END_REF][START_REF] Ohta | Evidence for interlayer coupling and moiré periodic potentials in twisted bilayer graphene[END_REF][START_REF] Cherkez | Van hove singularities in doped twisted graphene bilayers studied by scanning tunneling spectroscopy[END_REF] (3) For the lowest θ , θ <∼ 2 • , almost flat bands appear and result in electronic localization: states of similar energies, belonging to the Dirac cones of the two layers interact strongly. In this regime, the velocity of states at Dirac point goes to almost zero for specific angle so-called magic angles. [START_REF] Trambly De Laissardière | Localization of dirac electrons in rotated graphene bilayers[END_REF][START_REF] Bistritzer | Moiré bands in twisted double-layer graphene[END_REF][START_REF] Trambly De Laissardière | Numerical studies of confined states in rotated bilayers of graphene[END_REF] In this paper, we study the consequence of the tunable effective coupling between layers by angle θ with intermediate values, ∼ 2 • < θ <∼ 20 • , on local density of states (LDOS) and transport properties. We combine tight-binding numerical calculations for commensurate tBLG and a perturbative continuous theory (see Appendix) that gives us deeper insight on θ effect. To analyze transport properties numerically in bulk 2D systems, we consider local defects, [START_REF] Lazar | Adsorption of small organic molecules on graphene[END_REF][START_REF] Brihuega | Selective hydrogen adsorption in graphene rotated bilayers[END_REF] 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 adsorbate is simulated by a simple vacancy in the plane of p z orbital as usually done. [START_REF] Castro Neto | The electronic properties of graphene[END_REF][START_REF] Trambly De Laissardière | Conductivity of graphene with resonant and nonresonant adsorbates[END_REF][START_REF] Missaoui | Numerical analysis of electronic conductivity in graphene with resonant adsorbates: comparison of monolayer and bernal bilayer[END_REF] 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. 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. [START_REF] Nair | Spin-half paramagnetism in graphene induced by point defects[END_REF][START_REF] Scopel | Interaction between single vacancies in graphene sheet: An ab initio calculation[END_REF] 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 considerthe case (i) where defects are located in the two layers with respect to the case (ii) where defects are located on one layer (layer 2) only.

In Sec. II, tight-binding (TB) local Density of states (LDOS) in pristine tBLG and the effect of disorder on total DOS (TDOS) are analyzed with respect to our analytical model for commensurate tBLG. The spatial modulation of the DOS shows an increase of the DOS in AA region of the Moiré. This is a precursor of the localization in the AA region for very low angles less than ∼ 2 • . [START_REF] Trambly De Laissardière | Localization of dirac electrons in rotated graphene bilayers[END_REF][START_REF] Trambly De Laissardière | Numerical studies of confined states in rotated bilayers of graphene[END_REF] The electrical dc-conductivity at high temperature (microscopic conductivity) is studied Sec. III A, and quantum corrections of conductivity (low temperature limit) are presented Sec III B. The method to compute de conductivity is given in the appendix A. Numerical resuts of the paper are analyzed using the analytical continous model presented in appendix B and C. This pertubative theory recovers known results for the velocity renormalization, 2,15 but provides new analytical results concerning LDOS and state lifetime versus θ values.

The method to built commensurate tBLG is well known and explained in many articles. Here we use the notations used in our previous papers [START_REF] Trambly De Laissardière | Localization of dirac electrons in rotated graphene bilayers[END_REF][START_REF] Trambly De Laissardière | Numerical studies of confined states in rotated bilayers of graphene[END_REF][START_REF] Trambly De Laissardière | Electronic properties of asymmetrically doped twisted graphene bilayers[END_REF] where each tBLG is built from two index n and m (table I). For |m -n| = 1 the cell of the bilayer contains one Moiré cell, whereas for |m -n| > 1 the cell of the bilayer contains several Moiré cells.

II. DENSITY OF STATES

A. Without defect

We first analyze the local density of states (LDOS) in pristine twisted bilayer graphene (tBLG) computed with the TB Hamiltonian detailed in the Refs 7. It is now well known theoretically [START_REF] Trambly De Laissardière | Numerical studies of confined states in rotated bilayers of graphene[END_REF][START_REF] Brihuega | Unraveling the intrinsic and ro-bust nature of van hove singularities in twisted bilayer graphene by scanning tunneling microscopy and theoretical analysis[END_REF][START_REF] Lopes Dos Santos | Continuum model of the twisted graphene bilayer[END_REF][START_REF] Trambly De Laissardière | Electronic properties of asymmetrically doped twisted graphene bilayers[END_REF] and experimentally [START_REF] Brihuega | Unraveling the intrinsic and ro-bust nature of van hove singularities in twisted bilayer graphene by scanning tunneling microscopy and theoretical analysis[END_REF][START_REF] Cherkez | Van hove singularities in doped twisted graphene bilayers studied by scanning tunneling spectroscopy[END_REF] that the energies E -and E + of Van Hove singularities vary linearly with the angle θ for θ >∼ 2 • . This is clearly seen in the LDOS on p z orbital of atom located at the center of AA area of the Moiré (Fig. 1). Since our TB Hamiltonian includes coupling beyond the first neighboring atoms, the electron / hole symmetry is slightly broken and E -is not strictly equal to -E + .

The LDOS in one layer of the bilayer as a function of position r in the Moiré structure is

ρ(E, r) = r|δ (E -H)| r . ( 1
)
To compare LDOS in bilayer with LDOS in monolayer we compute the relative variation of the LDOS due to interlayer interaction ∆ρ(E, r)/ρ m (E), with ∆ρ(E, r) = ρ(E, r)ρ m (E), where ρ m (E) is the LDOS in monolayer that does not depend on the position r. (12,13) and (6,7) bilayers at the energy E = E D + 0.05 eV. The strong increase of LDOS in AA areas with respect to the AB zone are clearly seen. As expected this difference between LDOS in AA area and AB area decreases as θ increases. Moreover our numerical TB calculation recovers the difference in the LDOS of the two inequivalent atoms in AB area. Indeed in AB area, as in AB Bernal stacking, C atoms lying above an C atom of the other layer have a lower LDOS than C atom that do not lied with a C atom of the other layer. That leads to a triangular contrast in the density map that has been observed in STM images ? in AB Bernard bilayer. Laurence as-tu une ref? je n'ai pas réussi à lire le nom de ce phénomènes dans des notes According to the perturbative analytical model presented in Appendix (Sec. C 5), the relative variation of the LDOS is independant of E for small E and it can be estimated by the simple formula,

∆ρ(E, r) ρ(E) θ 1 θ 2 6 ∑ j=1 cos( G j • r), (2)
where G j are 6 equivalent vectors of the reciprocal space of the Moiré lattice. The constant θ 1 is given by,

θ 1 = √ 2t (hcK D ) , (3)
where K D is the modulus of the wave-vector in Dirac point of graphene. Using the interlayer coupling value t 0.12 eV (Appendix Sec. B 1), one finds that the value of θ 1 is close to θ 1 1 • . Equation (2) does not depend on the type of atom ε (ε = A or B atom), it oscillates with G j as expected. As it is clear the maximum value is obtained for r = 0 which is at the center of AA area, and relative variation of the LDOS varies as θ -2 . As shown in Fig. 2 the overall agreement between TB numerical calculation and TB analytical model is very good. We just note a small triangular contrast in AB zone which is not reproduced by the analytical model (see Appendix for a discussion) [verifier qu'on en parle dans l'annexe]. We observe in particular a reinforcement of the DOS in the AA region and a lowering in the AB regions. This behavior is a precursor of the electronic localization in AA region which is observed in the very low angle limit θ < 2 • . [START_REF] Trambly De Laissardière | Localization of dirac electrons in rotated graphene bilayers[END_REF][START_REF] Trambly De Laissardière | Numerical studies of confined states in rotated bilayers of graphene[END_REF][START_REF] Huder | Electronic spectrum of twisted graphene layers under heterostrain[END_REF]

B. With resonant adsorbates

To study the effect of static defects on the electronic confinement by the Moiré we include atomic vacancies (vacant atoms) that simulate resonnant adsorbates atoms or molecules. [START_REF] Trambly De Laissardière | Conductivity of graphene with resonant and nonresonant adsorbates[END_REF][START_REF] Missaoui | Numerical analysis of electronic conductivity in graphene with resonant adsorbates: comparison of monolayer and bernal bilayer[END_REF][START_REF] Vitor | Disorder induced localized states in graphene[END_REF][START_REF] Vitor | Modeling disorder in graphene[END_REF][START_REF] Wehling | Resonant scattering by realistic impurities in graphene[END_REF][START_REF] Trambly De Laissardière | Conductivity of graphene with resonant adsorbates: beyond the nearest neighbor hopping model[END_REF][START_REF] Fan | Anderson localization in two-dimensional graphene with short-range disorder: One-parameter scaling and finite-size effects[END_REF][START_REF] Missaoui | Mobility gap and quantum transport in a functionalized graphene bilayer[END_REF] For each vacancies concentrations c with respect to the total number of Carbon atoms in tBLG, we consider two cases:

(i) vacancies are randomly distibuted in both layers, (ii) vacancies are randomly distibuted in layer 2 only. Total DOSs (tDOSs) in (12,13) tBLG and (6,7) tBLG are drawn Fig. 3 for different concentrations of vacancies in cases (i) and (ii). For small c values, the Van Hove singularities are still clearly seen but they are enlarged by disorder. This shows that static disorder destroys the confinement by Moiré in AA areas. For c >∼ 0.5 % peaks of the Van Hove singularities are destroyed by vacancy states. With TB Hamiltonian including only first neighbor hopping terms, the vacancy states are midgap states at Dirac energy. [START_REF] Vitor | Disorder induced localized states in graphene[END_REF][START_REF] Vitor | Modeling disorder in graphene[END_REF] But, as in monolayer graphene [START_REF] Trambly De Laissardière | Conductivity of graphene with resonant adsorbates: beyond the nearest neighbor hopping model[END_REF] and Bernal bilayer graphene, [START_REF] Missaoui | Numerical analysis of electronic conductivity in graphene with resonant adsorbates: comparison of monolayer and bernal bilayer[END_REF] the TB hoppings beyond first neighbor enlarges the midgap states and shift it to the negative energies, typically around -0.2 eV. As shown in Fig. 4, when vacancies are located on layer 2 only (case (ii)), the vacancy states is appear only on LDOS p z orbitals of layer 2. Note that average DOS in layer 1 is slightly modified by the vacancies located in layer 2 (Fig. 4). This effect seems similar to modification du to nonresonant scatters. [START_REF] Trambly De Laissardière | Conductivity of graphene with resonant and nonresonant adsorbates[END_REF] Figs. 3 and4 show that, as far as the DOS is concerned, the rotated angle θ does not change the effect of vacancies. As we will see in next section, the effect of θ is more pronounced on wave-packet quantum diffusion and thus on transport properties. D'AUTRES COMMENTAIRES ???? III. QUANTUM TRANSPORT Within the Kubo-Greenwood formalism we compute the conductivity σ (E F ) versus the Fermi energy E F using the real space method developped by Mayou, Khanna, Roche and Triozon, [START_REF] Mayou | Calculation of the conductivity in the short-meanfree-path regime[END_REF][START_REF] Mayou | A real-space approach to electronic transport[END_REF][START_REF] Roche | Conductivity of quasiperiodic systems: A numerical study[END_REF][START_REF] Roche | Formalism for the computation of the rkky interaction in aperiodic systems[END_REF][START_REF] Franc ¸ois Triozon | Quantum dynamics in two-and three-dimensional quasiperiodic tilings[END_REF] in the famework of the Relaxation Time Approximation (RTA) to account [START_REF] Trambly De Laissardière | Conductivity of graphene with resonant and nonresonant adsorbates[END_REF] effects of inelastic scatters due to electron-phonon interactions (see Appendix A). Elastic scattering events due to local defects (vacant atoms) are included in the Hamiltonian itself in a large unit cell containing more than 10 7 atoms with boundary periodic conditions.

A. High temperature conductivity

We first consider the high temperature case (or room temperature case) where the inelastic scattering time τ i is close to the elastic scattering time τ e due to static defects. In that case, the conductivity is called microscopic conductivity, σ M , because it takes into account quantum interference effects accuring during time less or equal to τ e τ i . σ M is close to semiclassical conductivity that does not take into account the quantum corrections due to multiple scattering effects. Typically, this quantity represents a room temperature conductivity when multiple scattering effects are destroyed by dephasing due to the electron-phonon interactions. In Fig. 5, σ M (E) is shown for three tBLG with rotated angle θ equal to 2.656 • , 5.086 • and 9.430 • , respectively, and in Fig. 6, σ M (θ 2 ) is shown for different energy values close to the Dirac energy E D .

For vacancies distribution (i) -ie vacancies randomly distruted in the two layers-, σ M (E) is almost independent of θ value. When vacancies concentration c is large (Fig. 6, c = 1% and 2%) behavior is similar to that of MLG and σ M 2σ M,MLG , where σ M,MLG is the conductivity for MLG. [START_REF] Trambly De Laissardière | Conductivity of graphene with resonant and nonresonant adsorbates[END_REF] σ M,MLG reaches to the well known universal minimum of the conductivity so-called conductivity "plateau" -independent of defect concentration-at energies arroud E D . [START_REF] Castro Neto | The electronic properties of graphene[END_REF] For smaller concentration (Fig. 6, c = 0.5%), σ M increases when the concentration c increases. These two regims, are similar to that found in AB Bernal bilayer graphene. [START_REF] Missaoui | Numerical analysis of electronic conductivity in graphene with resonant adsorbates: comparison of monolayer and bernal bilayer[END_REF] Roughly speaking, for large c values, the elastric mean free path L e in MLG (see Fig. 4(a) in Ref. 48) is smaller than the average traveling distance [START_REF] Missaoui | Numerical analysis of electronic conductivity in graphene with resonant adsorbates: comparison of monolayer and bernal bilayer[END_REF] l 1 in a layer between two interlayer hoppings of the charge carriers, and thus carriers behaves like in MLG. Whereas for small c values, L e > l 1 and thus interlayer hopping are involved in the diffusive regim and BLG conductivity properties are different that MLG ones.

For vacancies distribution (ii) -ie vacancies randomly distributed in layer 2-, and large rotated angle (Fig. 5(c)), conductivity is larger than in the first case (i). Indeed for large θ , typically θ > 10 • , eigenstates are located mainly in one layer ("decoupled" layers) [START_REF] Trambly De Laissardière | Localization of dirac electrons in rotated graphene bilayers[END_REF][START_REF] Trambly De Laissardière | Electronic properties of asymmetrically doped twisted graphene bilayers[END_REF] and thus conductivity of the bilayer is the sum of the conductivity of two almost independent layers,

σ M σ M,1 + σ M,2 , (4)
corresponding to conductivity of layer 1 and 2, respectively. Conductivity of layer with defects is close to MLG conductivity σ M,2 σ M,MLG and conductivity of layer without defects σ M,1 is affected by the presence of defects in layer 2. With increasing θ , the eigenstates are more and more located on one layer, thus layers are more and more decoupled, and the σ M,1 increases as layer 1 becomes more and more like a pristine MLG. Consequently the conductivity of the tBLG increases when θ increases. In theses cases numerical results (Figs. 6) show that σ M increases like θ 2 . For small angles (Fig. 5(a) and Fig. 6), eigenstates are located almost equally on both layer for all energies around Dirac energy; [START_REF] Trambly De Laissardière | Electronic properties of asymmetrically doped twisted graphene bilayers[END_REF] therefore they are affected in a similar way by the two kinds of vacancies distributions (i) and (ii). Conductivity is thus very similar in the two cases.

The analytical model presented in Appendix Sec. C 4, allows to understand why σ M increases as θ 2 increases when defects are located only in layer 2 (cases (ii)). From Einstein conductivity formula, conductivity in layer p, p = 1, 2, is

σ M,p (E) = e 2 n p (E)V 2 τ p , (5)
where ρ p and τ p are the average DOS in layer p and the average elastic scattering time in layer p, respectively. For energy values in the plateau of conductivity arround E D , the layer 2 -with defects-has conductivity close to universal minimum of MLG, [START_REF] Trambly De Laissardière | Conductivity of graphene with resonant and nonresonant adsorbates[END_REF] σ M,2 (E) σ M,MLG , thus from equations ( 4) and ( 5), the conductivity in the bilayer is

σ M (E) σ M,MLG 1 + ρ 1 (E) ρ 2 (E) τ 1 τ 2 , (6)
where the ratio between scattering times can be estimated from the formula (C17) obtained in the Appendix. Thus, with θ 0 related to θ 1 (Appendix equation (C2)),

σ M (E) σ M,MLG 1 + ρ 1 (E) ρ 2 (E) θ 2 θ 2 0 , (7)
θ 0 = √ 3θ 1 , (8)
i.e. θ 0 2 • (Appendix Sec. C 3). Since ρ 1 (E) and ρ 2 (E) are different (Fig. 4) and depend on the energy values and the defect concentration c, the slope of σ M versus θ 2 also depends on E and c (Fig. 6).

B. Low temperature conductivity

In the low temperature limit, inelastic scattering time τ i is larger than elastic scattering time τ e , and multiple scattering effects may reduced the conductivity with respect to microscopic conductivity σ m . The average inelastic length L i thus satisfies L i L e and L i l 1 . τ i and L i increase when temperature increases. To evaluate this effect we compute [START_REF] Trambly De Laissardière | Conductivity of graphene with resonant and nonresonant adsorbates[END_REF][START_REF] Missaoui | Numerical analysis of electronic conductivity in graphene with resonant adsorbates: comparison of monolayer and bernal bilayer[END_REF] the conductivity σ versus L i at every energy E (Fig. 7) for the two vacancy distribution cases ((i) in layer 1 or (ii) in the two layers). As expected in disordered 2D systems, [START_REF] Lee | Disordered electronic systems[END_REF] for large L i , σ (L i ) follows a linear variation with the logarithm of L i , like in the case of monolayer graphene [START_REF] Trambly De Laissardière | Conductivity of graphene with resonant and nonresonant adsorbates[END_REF][START_REF] Trambly De Laissardière | Electronic transport in graphene: Quantum effects and role of loacl defects[END_REF] and Bernal bilayer graphene, 48

σ (E, L i ) = K -αG 0 log(L i ), (9)
where K is a constant depending on σ M and L e , and the slope α is almost independent on energy E, the defect concentration and the repartition of the defects (in on layer or in both layers).

From numerical results one obtains α ??? [Ahmed: valeur de alpha ?] which is close to monolayer value [START_REF] Trambly De Laissardière | Conductivity of graphene with resonant and nonresonant adsorbates[END_REF] and Bernal bilayer value. [START_REF] Missaoui | Numerical analysis of electronic conductivity in graphene with resonant adsorbates: comparison of monolayer and bernal bilayer[END_REF] Localization length ξ can be estimated from the equation σ (L i = ξ ) = 0 and the linear extrapolation of σ versus log L i at large L i (see dashed lines Fig. 7). ξ versus θ for various energies in the plateau of conductivity are shown in Fig. 8. As σ M , ξ is almost independent of θ when defects are located in both layers, but ξ increases strongly when defects are located in one layer only.

Commentaires ? conclusion sur conductivity a basse temperature

CONCLUSION

We have presented a numerical study of the local electronic density of states (LDOS) and the conductivity in pristine and covalently functionalized twisted graphene bilayers (tBLG), with an angle of rotation θ > 2 • . Those results are understood using a perturbative analytical model described in the Appendix. The atomic structure in Moiré induces a strong modulation in the LDOS between AA stacking areas and AB stacking areas, which varies as θ -2 following a simple analytic expression. We show that disorder breaks the interlayer effective coupling due to Moiré pattern. Therefore when defects are randomly distributed in both layer, the conductivity σ M is almost independent of θ , whereas σ M ∼ θ 2 when defects are randomly distributed in one layer only. Such a nonsymmetric distribution of defects may often occur in experimental situation because of the effect of substrate, adatoms or admolecules. Finally the quantum correction to the conductivity are computed and localization length is calculated versus θ .

IV. ACKNOWLEDGMENT

The authors wish to thank C. Berger, W. A. de Heer, P. Mallet and J.-Y. Veuillen, V. Renard, C. Chapelier for fruitful discussions. TB calculations have been performed at the Centre de Calculs (CDC), Université de Cergy-Pontoise. We thank Y. Costes and Baptiste Mary, CDC, for computing assistance. DFT calculations was performed using HPC resources from GENCI-IDRIS (Grants A0060910784 and Laurence numero de ton projet ). This work was supported by the ANR project J2D (ANR-15-CE24-0017) and the Paris//Seine excellence initiative.

Appendix A: Kubo-Greenwood conductivity

In Kubo-Greenwood approach for transport properties, the quantum diffusion D, is computed by using the polynomial expansion of the average square spreading, ∆X 2 , for charge carriers. This method, developed by Mayou, Khanna, Roche and Triozon, [START_REF] Mayou | Calculation of the conductivity in the short-meanfree-path regime[END_REF][START_REF] Mayou | A real-space approach to electronic transport[END_REF][START_REF] Roche | Conductivity of quasiperiodic systems: A numerical study[END_REF][START_REF] Roche | Formalism for the computation of the rkky interaction in aperiodic systems[END_REF][START_REF] Franc ¸ois Triozon | Quantum dynamics in two-and three-dimensional quasiperiodic tilings[END_REF] allows very efficient numerical calculations by recursion in real-space that take into account all quantum effects. Static defects are included directly in the structural modelisation of the system and they are randomly distributed on a supercell containing up to 10 7 Carbon atoms. Inelastic scattering is computed [START_REF] Trambly De Laissardière | Conductivity of graphene with resonant and nonresonant adsorbates[END_REF] within the Relaxation Time Approximation (RTA). An inelastic scattering time τ i beyond which the propagation becomes diffusive due to the destruction of coherence by inelastic process. One finally get the Einstein conductivity formula, [START_REF] Trambly De Laissardière | Conductivity of graphene with resonant and nonresonant adsorbates[END_REF] σ

(E F , τ i ) = e 2 n(E F )D(E F , τ i ), (A1)
where E F is the Fermi level, D(E, τ i ) is the diffusivity (diffusion coefficient at energy E and inelastic scattering time τ i ),

D(E, τ i ) = L 2 i (E, τ i ) 2τ i , (A2)
n(E) is the density of states (DOS) and L i (E, τ i ) is the inelastic mean-free path. L i (E, τ i ) is the typical distance of propagation during the time interval τ i for electrons at energy E,

L 2 i (E, τ i ) = 1 τ i ∞ 0 ∆X 2 (E,t) e -t/τ i . (A3)
Without static defects (static disorder) the L i and D goes to infinity when τ i diverges. With statics defects, at every energy

E, σ (τ i ) reaches a maximum value, σ M (E F , τ i ) = e 2 n(E F ) Max τ i {D(E F , τ i )} , (A4)
called microscopic conductivity. σ M corresponds to the usual semi-classical approximation (semi-classical conductivity). This conductivity is typically the conductivity at room temperature, when inelastic scattering τ i (inelastic mean free path L i ) is close to elastic scattering τ e (elastic mean free path L e ), τ e (E) = L e (E)/v(E) and L e (E) = D M (E)/2v(E), where D M (E) is the maximum value of D(τ i ) at energy E and v(E) the velocity at very small times (slope of ∆X(t)).

For larger τ i and L i , τ e τ i and L e L i , quantum interferences may result in a diffusive state, D(τ i ) D M , or a subdiffusive state where D decreases when τ i and L i increase. For very large L i , L i close to localization length ξ , the conductivity goes to zero.

Appendix B: Tight-Binding Model

Real space couplings

In the tight-binding (TB) scheme only p z orbitals are taken into account since we are interested in electronic states close to the Fermi level. The TB model used in this paper is the same as in our previous work on twisted bilayers graphene [START_REF] Trambly De Laissardière | Localization of dirac electrons in rotated graphene bilayers[END_REF][START_REF] Trambly De Laissardière | Numerical studies of confined states in rotated bilayers of graphene[END_REF][START_REF] Trambly De Laissardière | Electronic properties of asymmetrically doped twisted graphene bilayers[END_REF] and in AB Bernal bilayer graphene. [START_REF] Missaoui | Numerical analysis of electronic conductivity in graphene with resonant adsorbates: comparison of monolayer and bernal bilayer[END_REF][START_REF] Missaoui | Mobility gap and quantum transport in a functionalized graphene bilayer[END_REF] The Hamiltonian has the form,

H = ∑ i ε i |i i| + ∑ (i, j) t i j |i j|, ( B1
)
where i is the p z orbital located at r i the on-site energy ε i , and the sum runs over all neighboring i, j sites. t i j is the hopping element matrix between site i and site j, computed from the usual Slater-Koster parameters as given in Ref 7. Since the layers are rotated, interlayer neighbors are not on top of each other (as is the case of the Bernal AB stacking). Therefore, the interlayer interactions are then not restricted to ppσ terms but ppπ terms have also to be introduced. [START_REF] Trambly De Laissardière | Localization of dirac electrons in rotated graphene bilayers[END_REF][START_REF] Trambly De Laissardière | Numerical studies of confined states in rotated bilayers of graphene[END_REF] Moreover interactions are not restricted to first neighbouring orbitals and they decrease exponentially with the interatomic distance. A cutoff distance r c is introduced which should be as large which must be large enough so that the results do not depend on it.

We have check that r c = 0.6 nm is enough. Note that if r c is to small non-physical small gap may appeared at the Dirac energy as shown in Fig. 9.

The matrix element of the interlayer Hamiltonian H c between one orbital at r in plane 1 and one orbital at r in plane 2 is given by essential in the following is also real and depends only on the modulus of the wave vector. From Fourier transformation we write

r |H c | r = H c (| r -r |). ( B2
H c ( r) = H c ( k) e i k• r d 2 k , (B3)
and

H c ( k) = 1 (2π) 2 H c ( r) e -i k• r d 2 r. (B4)
Here also the coupling H c ( k) decreases when | k| increases. We shall see below that the important value of H c ( k) is for | k| close to the modulus of a Dirac point which is represented by K D in Fig. 10.

Interlayer Coupling between Bloch states

We want to compute the coupling between two Bloch states of layer 1 and layer 2. Each graphene layer are honeycomb lattice with two atoms, atoms A and atoms B, in a unit cell. Let us consider normalized Bloch states made of atomic p z orbitals A or B in layer α, α = 1 or 2,

|A k α = 1 √ N ∑ R e i k• R A |A R α , (B5)
|B k α = 1 √ N ∑ R e i k• R B |B R α , ( B6
)
where N is the number of unit cells of the crystal and the summation is performed on all cells of crystal ( R). In the following A or B are indicated by ε according to the following convention,

ε = A for A atom B for B atom (B7) α = 1 lower layer 2 upper layer (B8)
The positions of the atoms in layer 1 are,

r ε R = R if ε = A r ε R = R + u if ε = B (B9)
and in layer 2,

r ε R = R if ε = A r ε R = R + u if ε = B (B10)
where u and u are vectors connecting the two atoms in the unit cells, i.e. A and B atoms in upper and A and B atoms in down layers respectively. Writing

H c |ε k 2 = ∑ i t(ε i k i , ε k) |ε i k i 1 , (B11)
where t(ε i k i , ε k) ≡ t i is the transfer matrix element. We find a selection rule such that Finally for k i = k + K r = k (mod K r ), we derive formula for coupling matrix after some calculations [START_REF] Faizy | Electronic structure and quantum transport in graphene nanostructures[END_REF] we switch to the following expression of the Hamiltonian, where K r and K r are vectors of reciprocal lattices,

k + K r = k + K r . (
t i ( k + K r ) = 4π 2 S H c ( k + K r ) e i ( k+ K r )•(ε u -ε u+ ∆) . (B13)
S is area of unit cell, ∆ is translation between the two layers. However this translation of the two layers just translate the overall Moiré pattern and can be set to zero without loss of generality. By symmetry of interaction between two orbitals, coupling depends only on the modulus of k + K r i.e H c ( k

+ K r ) H c (| K D + K r |).
The modulus of t i is represented in Fig. 10. One sees that the most important value of |t i | is one corresponding to the smallest possible value of k + K r . By careful examination it can be shown that for electronic states close to the Dirac point this minimum corresponds to the modulus of wave-vector in Dirac point (K D = || K D || 17.2 nm -1 ). From Fig. 10, it is easy to compute coupling value close Dirac is around t 1 0.12 eV. All the other contributions are much smaller and will be neglected here.

Selecting only this contribution means that K r is such that k + K r belongs to one of three equivalent valleys. Therefore a set of two Bloch states with a given wave vector (equations (B5) and (B6)) in one plane will be coupled to three sets of two Bloch states in other plane corresponding to three different wave vectors. This simplifies much the structure of Hamiltonian and the analytical calculations presented here.

In the following we shall count the vectors k and k from their respective Dirac point K D1 and K D1 . K D1 is obtained from K D1 by a rotation of an angle θ around the vector ζ which is perpendicular to the planes 1 and 2. Therefore one has

k = δ k + K D1 , (B14)
k = δ k + K D1 . (B15)
Finally one get for the selection rule

δ k δ k -θ ζ × K Di , (B16)
where the indice i takes the values i = 1, 2, 3. K Di and K Di are the three equivalent Dirac point in plane 1 and 2. K Di is obtained from K Di by a rotation of an angle θ around the vector ζ which is perpendicular to the planes 1 and 2.

or self-energy. From this effective Hamiltonian we shall get the velocity renormalization, the electron lifetime in plane 1 due to disorder in plane 2 and the modulation of the density of states close the charge neutrality point. The theory that is developed here is perturbative and assumes that the rotation angle θ is not too small. In particular we emphasize that the perturbation theory is valid for

z,t, ∆ hcK D θ , (C1)
where K D = || K D ||, z is the energy of calculation, t is the interlayer coupling (t t 1 0.12 eV, Sec. B 1) and ∆ is a possible difference in on-site energy between the two layers. The condition on t implies that θ > θ 1 where

θ 1 = √ 2t (hcK D ) . ( C2
)
The value of θ 1 is close to θ 1 1 • . The condition on z implies that the current energy at which the quantities are calculated is smaller than the typical energy of the Van Hove Singularities (VHS) which depends linearly on θ . The difference in energy ∆ of the two layers must also be smaller than the energy of the VHS.Note that the VHS have been clearly observed with STM experiments on twisted graphene bilayer.

Effective one-plan Hamiltonian

We consider first a Bloch state in plane 1 with wave vector δ k 0 . It can be coupled to a Bloch state δ k in plane 2 then propagates freely in plane 2 and is scattered again to a Bloch state in plane 1 with a wave vector δ k f . Applying the selection rule (B16) to each interplane interaction we find that δ k f and δ k 0 are related by

δ k f δ k 0 -θ ζ × ( K Di -K D j ). (C3)
Therefore the coupling between planes 1 and 2 induces an effective coupling between Bloch states of plane 1 with the selection rule (C3). Note that the indices i and j take the values i = 1, 2, 3. When i = j a Bloch state with δ k 0 is coupled only to the Bloch states with the same wave-vector δ k f = δ k 0 . This process gives a self-energy which renormalizes the energy of the state of the single plane 1 (see below).

When i and j are different then δ k f and δ k 0 are different,

δ k f δ k 0 + G k . (C4) G k = θ ζ × ( K Di -K D j
) is a reciprocal lattice vector of the Moiré lattice, where K Di -K D j is a reciprocal lattice vector of graphene. These vectors takes six possible values named G k in the main text. These are vectors of the reciprocal lattice of the Moiré pattern. As we show below this coupling between Bloch states of different wave vector will create eigenstates with mixing of different oscillating components which leads to oscillations in the DOS with wave-vectors components G k (see below). We note also that the coupling introduces only small spatial frequencies and in particular it does not connect states of the two non equivalent Dirac cones.

Self-energy

We are interested in the self-energy of coupling of states in upper plane due to the coupling with states of lower plane. Indeed the real-part of self-energy ℜσ (z) is associated to modification of dispersion relation and will allow us to discuss velocity renormalization. The imaginary part of self-energy is associated to the electron lifetime. It will allow us to discuss lifetime of the electron in one plane when there is disorder in other plane.

Using matrix notations defined in Appendix B we have

Σ 1 (z) = ∑ K r T + ( K r ) G 2 ( K D + K r ) T -( K r ), (C5)
where K r is the vector of reciprocal lattice which has three values that connect one Dirac point to itself or to the two other equivalent Dirac points. T describes the coupling between two plane and Green operator at wave vector

θ ζ × K dµ is G 2 (z, θ ζ × K dµ ) = 1 z -H -(θ ζ × K dµ ) , (C6)
where K dµ counts the three equivalent Dirac points. And for the Hamiltonian

H 2 θ ζ × K dµ ) = ∆ -γ 0 f (θ ζ × K dµ ) -γ 0 f * (θ ζ × K dµ ) ∆ , ( C7
) where ∆ is potential difference between two layers (plane 1 is in potential 0 and plane 2 is in potential ∆). Note that this matrix is evaluated at θ ζ × K dµ . Indeed for k sufficiently close to Dirac point k, because h

(| k -K d |) γ 0 | f (θ ζ × K dµ )| and
we can neglect the dependence on the k in H 2 , G 2 and Σ 2 (z). This corresponds to the general conditions of validity of the presnet perturbation theory (see above the introduction of appendix C).

So now after some calculations we get for the self-energy

Σ 0 (z) = σ (z) I, (C8) with σ (z) θ 2 0 θ 2 ∆ -z , (C9)
where we have introduced θ 0 ,

θ 0 = 3 √ 2π t γ 0 . (C10)
Using the values of t t 1 0.12 eV (Sec. B 1) and γ 0 2.7 eV one finds that the value of the angle θ 0 is θ 0 1.7 • .

Velocity renormalization

The eigenvalues are the poles of the Green's function. Therefore the energy E( k) is given by Eσ (E) = ±hv| k|.

(C11)

For | k| = 0, we have solution E = E 0 such that

E 0 -σ (E 0 ) = 0. (C12)
For small k, we can write E( k) = E 0 + δ E( k). Eventually we have a nice formula:

δ E = ±hv| k| 1 -σ (E 0 ) . ( C13
)
Finally the renormalized velocity r is Therefore using a well established tight-binding model, we recover velocity renormalization consistent with that of Ref. 2. In addition we find that this velocity renormalization is independent of the difference in potential of two planes. As it is shown in Fig. 11, a systematic study of the renormalization of the velocity close to the Dirac point is done [START_REF] Trambly De Laissardière | Localization of dirac electrons in rotated graphene bilayers[END_REF] , compared to its value in a monolayer graphene, for rotation angles θ varying between 0 • and 60 • (Fig. 11 ). The renormalization of the velocity varies symmetrically around θ = 30 • . Indeed, the two limit cases θ = 0 • (AA stacking) and θ = 60 • (AB stacking) are different, but Moiré patterns when θ → 0 • and when θ → 60 • are similar because a simple translation by a vector transforms an AA zone to an AB zone.

r = 1 1 + θ 2 0 /θ 2 . ( C14
Focusing on angles smaller than 30 • , it is defined 3 three regimes as a function of the rotation angle θ (Fig. 11). For large θ (20

(3 • ≤ θ ≤ 20 • )
the perturbative theory of Lopez dos Santos et al. predicts correctly the velocity renormalization which is also in accordance to the above formula equation(C14). For the small rotation angles (θ < 3 • ) a new regime occurs where the velocity tends to zero and perturbation theory can't be applied.

Electron lifetime

The two planes of the bilayer can have very different amount of disorder due to their different exposure to environment . For example the lower plane will be in contact with a substrate and the upper plane is exposed either to vacuum or to a gas (sensor application). Therefore it is of high interest to consider the limit case where defects are present in one plane and absent from the other plane. In the following we consider that defects are present only in the lower plane. If the two planes were decoupled, defects in one plane would affect electron lifetime in that plane but not in other one. Since the planes are coupled defects in one plane will also affect electronic lifetime in the other plane. In this chapter we discuss how such a repartition of defects impacts the electron lifetime. In chapter III based the present results we shall discuss how electron lifetime affects the overall electronic conduction of the bilayer.

If there is disorder in the lower layer (layer 1) the Bloch states of this layer will have a contribution to their self-energy which is imaginary. This can be represented in the simple possible model by a purey imaginary part of the potential evergy ∆,

∆ = -i h τ 1 , ( C15
)
where τ 1 is the lifetime in the layer 1 due to disorder in the layer 1. Using formula (C9) we see that electrons in the upper plane 2 acquire an imaginary self-energy

ℑσ (z) = - ih τ 2 = - ih τ 1 θ 2 0 θ 2 . ( C16
)
Therefore the lifetimes τ 1 and τ 2 in the layer 1 and layer 2 are related through:

τ 2 τ 1 = θ θ 0 2 , ( C17
)
where θ 0 is given by equation (C10), and is same quantity as in the velocity renormalization expression (C14).

Spatial variation of density of states

As explained above the coupling between Bloch states of different wave-vectors in plane 1 (due to interlayer coupling with plane 2) corresponds to the selection rule

δ k f δ k 0 + G k , ( C18
)
where the vector G k such that G k = θ ζ × ( K Di -K D j ) is a reciprocal lattice vector of the Moiré lattice. The typical difference in energy between Bloch states of δ k f and of δ k 0 is ∆E hcθ |K Di |. This difference is nearly independent of δ k 0 provided that it is sufficiently close to zero. The typical coupling is t e f f t 2 /(hcθ |K Di |).

Then the mixing between states of wave vector close to ( K Di ) and wave vector close to K Di + G k will be of order t e f f /∆E i.e. of order (θ 1 /θ ) 2 . We note that the mixing is independent of δ k 0 provided that it is sufficiently close to zero. Therefore the relative variation of the DOS of a state is independent of this state and then the relative variation of the DOS is constant sufficiently close to the Dirac point and depends only on the position in the Moiré pattern. After a lengthy calculation the precise calculation provides the expression given in the main text (equation ( 2)), ∆ρ(E, r)

ρ(E)

θ 1 θ 2 6 ∑ j=1 cos( G j • r), (C19)
where G j are 6 equivalent vectors of the reciprocal space of the Moiré lattice and where the rotation angle θ 1 is given by

θ 1 = √ 2t (hcK D ) = θ 0 √ 3 . (C20)
Using the interlayer coupling value t 0.12 eV (Appendix Sec. B 1) one finds that θ 1 is close to one degree. We emphasize that the present theory is perturbative in the coupling t. This perturbation theory is valid for sufficiently large values of θ as explained in the introduction of appendix (C).The other assumption is to neglect Fourier components of the interlayer Hamiltonian that couple a Bloch state with other states having wave vectors away from the Dirac cones. This approximation can lead to the under estimation of modulations of the DOS at spatial frequencies high with respect to the Moiré period. This could explain why the DOS modulation on sub lattices A and B can differ by about ±15% as compared to averaged DOS whereas the present perturbative theory does not predict this difference. Note that the average DOS of two neighboring A and B atoms is well reproduced by the analytical model.

FIG. 1 .

 1 FIG. 1. (color online) Local density of states (LDOS) at the center of a AA zone in pristine tBLG listed table I for tBLG with different rotation angle θ [Deg.]. Some LDOS are taken from Ref. 7. E D = 0.

2 ) 5 FIG. 2 .

 252 FIG. 2. (Color online) Relative variation of the LDOS on top layer at energy E = E D + 0.05 eV, close to the Dirac energy E D , in (a) (12,13) tBLG and (b) (6,7) tBLG: (a.1) (b.1) TB results and (a.2) (b.2) analytic results from equation (2). To be compared with analitic results the TB plots are made by a continuous extrapolation of LDOS on atomic orbitals. The same arbitrary unit are used for all the 4 LDOS. E D = 0.The LDOS on each carbon atoms of Moiré has been calculated using TB. So that density map ρ(E, r) where r are the

2 FIG. 3 .FIG. 4 .

 234 FIG. 3. (color online) Total DOS in (a) (12,13) tBLG and (b) (6,7) tBLG, for various concentrations c (%) of vacancies with respect to the total number of atom in tBLG: (Dashed line) with vacancies in both layers and (full line) with vacancies in layer 2. (Dotted line) is the DOS in pristine monolayer graphene (MLG). E D = 0.

FIG. 5 .

 5 FIG. 5. (color online)Microscopic conductivity σ M in (a)(12,13) tBLG, (b) (6,7) tBLG, (c)(3,4) tBLG, for the two cases: (Full line) with vacancies in layer 2 and (dashed line) with vacancies in both layers. c is the concentration of vacancies with respect to the total number of atom in tBLG. Inserts: σ M around the Dirac energy E D = 0.

FIG. 6 .

 6 FIG. 6. (color online) Microscopic conductivity σ M in tBLG versus rotated angle θ 2 for energy values E: (Full line) with vacancies in layer 2, (Dashed line) with vacancies in both layers. c is the concentration of vacancies with respect to the total number of atom in tBLG. E D = 0.

FIG. 7 .

 7 FIG. 7. (color online) Conductivity in bilayer versus inelastic scattering L i , at the energies E, for concentration c = 2% of vacancies with respect to the total number of atom in bilayer: (θ 1 = 2.656 o ) (12,13) tBLG, (θ 2 = 5.086 o ) (6,7) tBLG, (θ 3 = 9.430 o ) (3,4) tBLG.(line) with vacancies in layer 2, (dashed line) with vacancies in both layers. For(3,4) tBLG (θ 3 = 9.430 o ) the localization regim appears at very large times for which very time consuming calculations are necessary; that is why this regime is only roughly estimated by extrapolation.

FIG. 8 .

 8 FIG.8. (color online) Localization length versus angle θ , at the energies E, for concentration c = 2% of vacancies with respect to the total number of atom in bilayer: (line) with vacancies in layer 2, (dashed line) with vacancies in both layers.

)FIG. 9 .

 9 FIG. 9. (color on line) Band dispersion E( k): (red points) DFT calculation using VASP code (for details on the calculation see Ref. 7), and (red lines) TB calculation, for (a) (1,3) tBLG (θ = 32.20 o ), (b) (1,4) tBLG (θ = 38.21 o ), calculated with a large interlayer cutoff distance r c = 0.60 nm, whereas Bleu dashed line is TB bands with a too small r c , r c = 0.34 nm. In the latter case a non-physical gap appears at energy E D = 0. Dot black line is MLG. Insert: Bands arround the K point. E D = 0.

FIG. 10 .

 10 FIG.10. Modulus of the interlayer coupling t versus in-layer distance r and modulus k of the wave-vector, calculated from the Tightbinding model described in Ref.7.

TABLE I .

 I Studied (n, m) bilayer structures. N is the number of atoms, θ the rotation angle.

	(n, m)	θ [deg.]	N
	(12,13)	2.656	1876
	(10,11)	3.150	1324
	(8,9)	3.890	868
	(6,7)	5.086	508
	(5,6)	6.009	364
	(4,5)	7.341	244
	(7,9)	8.256	772
	(10,13)	8.613	1596
	(3,4)	9.430	148
	(8,11)	10.417	1092
	(2,3)	13.174	76
	(5,9)	18.734	604

  )

		1
		0.8							θ = 1.6 Deg θ = 1.8 Deg
	mono	0.6							θ = 2.0 Deg θ = 2.2 Deg
	V bi / V	0.4							θ = 2.4 Deg θ = 2.6 Deg θ = 2.8 Deg
		0.2							θ = 3.0 Deg VASP TB
		0
		0	2	4	6	8	10	12	14	16	18	20
						θ (degres)

FIG. 11. Velocity ratio bi / mono for commensurate (n,m) bilayer cell versus rotation angle θ , computed from equation(C10) with varius θ 0 values. Circle, DFT calculation using VASP code; cross, TB calculations, from Ref. 3.

  • ≤ θ ≤ 30 • ) the Fermi velocity is very close to that of graphene. For intermediate values of θ

Appendix C: Effect of interlayer coupling

We consider a plane 1 coupled to plane 2 which is rotated by an angle θ with respect to plane 1. If one considers the time evolution within plane 1 or more generally the restriction of the total Green's function to plane 1 the coupling to plane 2 amounts to the addition of an effective Hamiltonian