Guilherme Ribeiro

Lorenzo Paulatto

Félix Mouhat

Tommaso Gorni

Renero-Lecuna Carlos

Perez-Villa Andrea

José Alexandre

Francesca Gualdi

Giulia Miozzi

Guillaume Consuma

David Morard

Ornella Santos-Cottin

Yannick Vaccarelli

Matteo Gillet

Benoit Giantomassi

Van Troeye

Sadia

Vinciane Viet-Anh David

Rebeh Zitouni

Younes Rachedi

Mabarek Belkadi

Khiredine Kadri

Hamide Houasse

Billel Khalef

Amine, Raphaelle, Anile, Edoardo Khiredine Mabrouk

I want to thank all my colleagues from the IMPMC and IMCN laboratories

 2. Composite bands are a set of bands isolated from the rest of the spectrum by a gap. The valence bands in an insulator is an example.

Introduction

When a medium is illuminated with light, many phenomena can be observed both at the surface and in the bulk (see Fig. A). At the surface, the incident radiation may be either reflected or transmitted to the bulk. Once inside the bulk, the transmitted light can be scattered, absorbed, or simply transmitted out of the material without any change. The absorbed radiation will be dissipated as heat for example. The propagating radiation inside the medium may be scattered either elastically or inelastically depending on whether or not the energy of the incident photon is conserved. At this level, the inelastically scattered light by lattice acoustic waves is usually referred to as Brillouin scattering (this phenomenon was independently discovered by Mandelstam). Scattering of light by other elementary excitations, such as optical phonons or plasmons, is known as Raman scattering (Sir C.V. Raman received the Nobel Prize in Physics in 1930 for discovering, in 1928 in Calcutta, India, the effect named after him). In this thesis, only phonons are considered. Within this assumption, a Raman scattering is said to be first, second, third.., n th order, depending on the number of emitted (Stokes processes) or absorbed (anti-Stokes processes) phonons [START_REF] Cardona | Light scattering in solids II: basic concepts and instrumentation[END_REF].

Figure A -Schematic diagram showing the linear optical processes that occur at the surface and in the interior of a medium. Adapted from Ref. [START_REF] Peter | Fundamentals of semiconductors: physics and materials properties[END_REF].

iii A Raman experiment consists in measuring the intensity of the scattered light as a function of the difference in frequency between the incoming and the scattered radiations (usually called the Raman shift). The energy of the laser light usually used in Raman experiments is of the order of 2 eV. In molecules and large gap insulators (∆ > 2 eV) the laser light is unable to excite the electrons over the gap. In diamond for example, the gap is 5.5 eV and a laser light of 2 eV will always excite the electrons into virtual states. These electrons will absorb or emit a certain number of phonons then return to the valence band by emitting the scattered radiation. Any other excitation energy that is much smaller than the material band gap, will always promote the electrons into virtual states and the same process of Raman scattering will take place. Since virtual states are not real states of the material, the Raman response is completely independent of the excitation energy. In this kind of situations, where the laser excitation energy is much smaller than the material band gap, Raman scattering is said to be non-resonant. On the contrary, in small gap insulators (PbTe, SnTe) or metals and semimetals (graphene, graphite, copper), the laser light can always promote electrons from the valence to the conduction band and the laser light is resonant with electron excitations in the material system. If we change the laser energy, the involved excited states change accordingly. The Raman response, now, depends on the excitation energy because the involved excited states are real states of the material system. In this kind of situations, where the incident laser energy matches an electronic transition in the material, Raman scattering is said to be resonant. Experimentally, this phenomenon of resonance is characterized by the appearance of sharp peaks in the Raman spectrum.

Resonance effects (which constitute a central point of the present thesis) are important phenomena to look for in Raman experiments. The shape, width, and intensity of the observed resonant peaks carry important informations about fundamental properties of the studied material (e.g. phonon frequencies, electron-phonon coupling ..etc.). This offers a very efficient way to quantify them using Raman spectroscopy. However, this quantification (or characterization) can not be achieved without a theoretical understanding of the resonant peaks (origin, shape, width.. etc.), in other words, without performing Raman simulations. The general aspects of the present thesis fall within this category of research, i.e. the theoretical understanding of resonant Raman spectra. The materials that are considered in this work are graphene, few-layer graphene, and graphite where resonance phenomena are present due to the semimetallic character of these systems. Before going into more details about the present thesis, we would like to give an overview about the different theoretical approaches used to simulate Raman spectra. Generally speaking, one can group these approaches into two main categories depending on the type of Raman scattering to be simulated whether it is resonant or not.

For non-resonant Raman scattering, the dielectric approach [START_REF] Lazzeri | First-principles calculation of vibrational raman spectra in large systems: Signature of small rings in crystalline SiO 2[END_REF][START_REF] Veithen | Nonlinear optical susceptibilities, raman efficiencies, and electro-optic tensors from first-principles density functional perturbation theory[END_REF], based on the concept of polarizability, has played a very important historical role in the theory of Raman scattering. Within this approach, one formulates the scattering cross sections in terms of an expansion of the polarizability tensor in power series of the normal coordinates. The thermal fluctuations induced by atomic vibrations in the medium (i.e. phonons) change the polarizability and consequently the frequency of the scattered light [START_REF] Knoll | Raman scattering of atomic vibrations in anharmonic potentials[END_REF]. Efficient first-principles computational methods have been implemented to evaluate v non-resonant Raman intensities for first order Raman scattering within the dielectric approach [START_REF] Lazzeri | First-principles calculation of vibrational raman spectra in large systems: Signature of small rings in crystalline SiO 2[END_REF][START_REF] Veithen | Nonlinear optical susceptibilities, raman efficiencies, and electro-optic tensors from first-principles density functional perturbation theory[END_REF]. The polarizability tensor (or equivalently, the dielectric function) can be computed from first principle within density functional perturbation theory using the 2n+1 theorem [START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF]. This fact has simplified enormously the computation of non-resonant Raman spectra. The most important application of this DFT-dielectric approach was the realization of the WURM project [START_REF] Bobocioiu | The wurm project -a freely available web-based repository of computed physical data for minerals[END_REF] which contains over 400 Raman spectra of different semiconductors and insulators. It should be highlighted, however, that near the resonance (when the laser energy is comparable to the material band gap) the dielectric approach calculations based on the standard approximation used within DFT (e.g. the local density approximation), become less reliable. In this situation, one needs more sophisticated approaches for the evaluation of Raman spectra. It is worth mentioning that, recently, a new approach to include excitons effects in Raman simulations within the dielectric approach has been established (see Refs. [START_REF] Gillet | First-principles study of excitonic effects in raman intensities[END_REF][START_REF] Del Corro | Atypical exciton-phonon interactions in WS 2 and WSe 2 monolayers revealed by resonance raman spectroscopy[END_REF] for more details).

For resonant-Raman scattering, Many Body Perturbation Theory (MBPT) [START_REF] Hedin | New method for calculating the one-particle green's function with application to the electron-gas problem[END_REF][START_REF] Hybertsen | Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies[END_REF] can also be used to evaluate Raman cross sections with resonance effects automatically included. In this approach, we express the Raman cross section as a function of electron-phonon and electron-photon interactions as well as any other coupling to elementary excitations one might be interested in. Such an approach has been formulated by Ganguly and Birman in 1967 [START_REF] Ganguly | Theory of lattice raman scattering in insulators[END_REF] for arbitrary electronic excitations which may include strongly correlated electron-hole pairs (strong excitons effects). The MBPT approach provides an accurate description of elementary excitations in materials (excitons, plasmons, magnons ..etc.) as well as the interaction between them. Hence, it provides us with a precise description of optical properties in general. Unfortunately, the computational load of a full many body calculation of Raman cross sections is extremely expensive and beyond the computing power afforded by High Performance Computing centers nowadays. Resonant Raman simulations based on MBPT calculations are still challenging and constitute and intense field of research.

Let us now turn to the discussion of more specific aspects of the present thesis. The three systems that are studied in the present work are : graphene, few-layer graphene (FLG), and graphite. Several theoretical and experimental studies have shown that the electronic an optical properties of these materials change drastically according to the number of layers and the way these layers are put on top of each others i.e. the stacking order. While counting the number of layers in few-layer graphene has been enormously simplified using Raman spectroscopy (the Raman 2D peak shape and width) after the pioneer work of A.C. Ferrari et al. [START_REF] Ferrari | Raman spectrum of graphene and graphene layers[END_REF], stacking order determination still a challenging problem in the field of few-layer graphene and graphite.

Stacking in FLG and graphite, which is the main subject of the present thesis, can be either of Bernal type (AB-stacking) or rhombohedral type (ABC-stacking) as shown in Fig. B. The Bernal stacking [START_REF] Bernal | The structure of graphite[END_REF] is the most stable form of FLG and graphite systems. Recently, however, rhombohedral stacked multi-layers graphene (RMG), attracted an increasing attention as theoretical calculations suggest the occurrence of a dispersionless electronic band (bandwith smaller than 2 meV) at the Fermi level [START_REF] Kopnina | Surface superconductivity in rhombohedral graphite[END_REF][START_REF] Xiao | Density functional investigation of rhombohedral stacks of graphene: Topological surface states, nonlinear dielectric response, and bulk limit[END_REF]. This flat band with extremely large effective mass indicates the existence of highly correlated system of electrons and make RMG a promising candidate for the appearance of highly correlated states of matter such as magnetism [START_REF] Pamuk | Magnetic gap opening in rhombohedral-stacked multilayer graphene from first principles[END_REF] or room-temperature su- perconductivity [START_REF] Precker | Identification of a possible superconducting transition above room temperature in natural graphite crystals[END_REF]. Thus, a powerful non-destructive characterization tool is highly needed in order to investigate the presence of the rhombohedral stacking in FLG and graphite. This fact constitutes our main motivation for studying these systems. The main objective of the present thesis is the understanding of the 2 nd order Raman peak in FLG and graphite (the so-called 2D peak) for all possible stacking and the extraction (from first-principle calculations) of the Raman fingerprints of each type of stacking. These fingerprints will be used as prerequisite to experimental nondestructive identification and synthesis of RMG and rhombohedral graphite using Raman spectroscopy. Given that we are interested in simulations of a resonant two phonons Raman peak, two main difficulties will be encountered :

-First, the fact that the 2D Raman peak is a two phonons process imposes a challenging difficulty that it is worth to explain here. In first order Raman scattering, only phonon from the Brillouin zone (BZ) center can contribute to the scattering cross section. This is the so-called fundamental selection rule which is due to momentum conservation. However, in 2 nd order Raman (i.e. two phonons peaks) the quantity that needs to belong to BZ center is the sum of the two phonons momenta. This does not impose any restriction on the phonon wave vector which can span, now, the entire BZ. Thus, for two phonons Raman scattering (using the dielectric approach or the many body approach) the response of the crystal must be calculated from the whole BZ.

-Second, the resonant character of the 2D peak implies the use of a many body approach to evaluate the Raman cross section. As said before, performing this kind of computation from first-principles is challenging. As an alternative, one can use tight-binding approaches to perform Raman simulations. These approaches are known to give remarkable results at least in the case of graphene and few-layer graphene (e.g. the work of P. Venezuela et al. [START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF] and the work V. Popov [START_REF] Popov | Two-phonon raman bands of bilayer graphene: revisited[END_REF]). Unfortunately, in this kind of approaches, the electron-phonon coupling which is a key ingredient in Raman simulations, is not screened. This may affect our vii simulations which is required to be as precise as possible in order to detect the tiny modification of the 2D peak that might be induced by a change in the stacking order. For that reason, our approach will be based on a an ab-initio, parameter-free, many body approach.

Our methodology to cope with the previous two difficulties will be essentially based on the use of Wannier functions to interpolate (with an ab-initio precision) the matrix elements needed in Raman simulations (e.g. electron-phonon and electron-light couplings). In particular, the electron-phonon coupling will be first calculated within density functional theory, then interpolated using Wannier functions. After interpolating all the required matrix elements, the calculation of Raman intensities will be carried out using the generalized Fermi golden rule which is the same procedure used in tight-binding approaches. In this way, we will have both the efficiency (the rapidity) of a tight-binding approach and the precision of an ab-initio calculation. The present thesis is organized as follows. In the first chapter, an overview about Raman spectroscopy and its applications in FLG will be presented with a special focus on graphene. The main motivations of the present work will be introduced at the end of this chapter. In the second chapter, we review the theory of Raman scattering and introduce the main challenges concerning the heavy computational load that characterizes Raman simulations of double resonant peaks. In the third chapter, Our methodology to cope with these challenges is presented and discussed. In the fourth and fifth chapters, we present the application of our methodology to calculate resonant Raman spectra in few-layer graphene (trilayer and tetralayer) and bulk graphites (Bernal, rhombohedral, and ABCB-graphites).

Chapter 1

Raman spectroscopy in graphene and few-layer graphene

In this introductory chapter we will present an overview of the use of Raman spectroscopy in graphene and few-layer graphene systems. The power of Raman spectroscopy as a structural and characterization tool will be presented and discussed. The chapter is organized as follows. First, general definitions and properties of carbonbased materials will be presented. Second, electronic and vibrational properties of graphene and few-layer graphene will be discussed within the framework of density functional theory. Third, the Raman effect will be introduced and the use of Raman spectroscopy in graphene and few-layer graphene will be discussed. Finally, the main subject and motivations of the present thesis which concern stacking order determination in few-layer graphene and graphite will be introduced and discussed. Although the present thesis has a theoretical background, the mathematical content of this chapter will be kept to a minimum in order to provide a clear view of the physical interests of the present work.

Carbon-based Materials

Carbon-based materials is a family of compounds formed by any stable arrangement of carbon atoms at different levels of hybridization (sp, sp 2 or sp 3 ). Diamond is an example of a stable sp 3 arrangement where the rigidity is simply due to the presence of the σ bonds that exist as a consequence of the sp 3 hybridization. Graphite (3 dimensional), graphene (2 dimensional), carbon nanotubes (1 dimensional ), and fullerene (0 dimensional) are well known examples of stable sp 2 arrangements. The sp arrangements are mostly found within the vast hydrocarbon family of carbyne. The physical and chemical properties of carbon-based materials are of great interest to the physical and industrial communities. In this chapter, we will be interested in sp 2 bonded carbon crystals (see Fig. 1.1) which have attracted a lot of attention within the condensed matter community. This interest was started by Bernal's discovery of the hexagonal symmetry of natural graphite and its layered structure (Bernal, 1924 [14]) during x-ray diffraction experiments (the studied graphite has since been named Bernal graphite). This family includes, in particular, graphite, graphene, few-layer rings can be seen (adapted from Ref. [START_REF] Odom | Atomic structure and electronic properties of single-walled carbon nanotubes[END_REF]). (g) Atomic structure of carbon nanotubes with the three possible arrangements of the carbon rings (armchair, zigzag and chiral).

graphene and carbon nanotubes. In the next paragraphs we will discuss some properties of these materials, before looking at how these properties can be studied by Raman spectroscopy.

Graphite is the hexagonal close-packed arrangement of carbon atoms. It belongs to the P 6 3 /mmc space group with four atoms per unit cell and a six-fold rotation axis along the stacking direction (see Figs. 1.1(b) and (c)). The cohesion of the whole crystal is ensured by the sp 2 bonds in the plane directions and by Van Der Waals interaction [START_REF] Bernal | The structure of graphite[END_REF][START_REF] Wallace | The band theory of graphite[END_REF] in the out-of-plane direction (c-axis direction). This gives rise to the so-called layered structure which is at the origin of the exceptional electronic and mechanical properties of graphite. From a thermal point of view, graphite has a high melting points (due to σ bonding) and a good thermal conductivity, making it suitable for heat transfer applications (e.g. nuclear reactors [24]). From a structural point of view, the possibility of intercalation in bulk graphite opened a vast field of research and gave rise to an independent class of materials called graphite intercalation compounds. The physical and chemical properties of these materials are tuned by controlling the nature and the amount of the intercalated element [25,[START_REF] Ohzuku | Formation of lithium-graphite intercalation compounds in nonaqueous electrolytes and their application as a negative electrode for a lithium ion (shuttlecock) cell[END_REF][START_REF] Inoshita | Electronic structure of potassiumgraphite intercalation compound: C 8 K[END_REF].

Another interesting aspect of graphite is the possibility to exfoliate a single or few-layer from the graphite bulk material. Over the last ten years, and following the discovery of graphene in 2004 [START_REF] Novoselov | Electric field effect in atomically thin carbon films[END_REF] (Nobel prize in 2010), this property has opened a vast field of research on the so-called 2-dimensional (2D) materials. The discovery of graphene has shown that 2D materials can exist under normal conditions of temperature and pressure. Since then, new types of 2D materials, such as phosphorene and transition metal dichalcogenides, have been isolated and studied.

Properties of graphene and few-layer graphene systems

Graphene is a two-dimensional arrangement of sp 2 bonded carbon atoms that follow a honeycomb pattern (see Fig. 1.1 (d)). It belongs to the P 6/mmm space group with two atoms per unit cell and a six-fold rotation axis in the out-of-plane direction. Fig. 1.1 (d) is a high resolution TEM (Transmission Electron Microscopy) image of oneatom-thick exfoliated sample, in which the hexagonal arrangement of carbon atoms is clearly observable. Graphene has two important properties that make it one of the most interesting materials nowadays. First, its high electronic mobility (15000-200000 cm 2 /V.s) exceeds any mobility of all the known semi-conductors. Second, its rigidity (that is due to the sp 2 bonding) makes it one of the hardest materials. These two properties are relevant features in a wide domain of applications such as optoelectronics, transistors, sensors, solar cells and many other areas [START_REF] Geim | The rise of graphene[END_REF][START_REF] Novoselov | Electric field effect in atomically thin carbon films[END_REF]30,31].

Few-layer graphene systems are systematically found during the exfoliation procedure from bulk graphite. Exfoliated samples may exhibit different thicknesses that are multiples of the inter-plane distance in graphite, as shown in Atomic Force Microscopy (AFM) measurements [32]. These systems are called multi-layer graphene or few-layer graphene. The electronic and optical properties of FLG have been shown to be different from those of graphene and depend of the number of layers (composing the FLG sample) and the stacking mode sequence [START_REF] Nilsson | Electronic properties of bilayer and multilayer graphene[END_REF]. Fig. 1.1 (e) shows a typical high resolution Scanning Tunneling Microscope (STM) image of a few-layer graphene sample.

Another type of sp 2 bonded carbon crystals is a one dimensional arrangement of carbon atoms to form the so-called carbon nanotubes. These systems are one dimensional materials in the sense that their wave function decays exponentially in the directions perpendicular to the nanotube. Thus, the charge carriers are confined to just one dimension. This suggests that 1D materials have many interesting features one can only find in one dimension (e.g. Van Hove singularities). A SWNT (single wall carbon nanotube) can be described as a single layer of a graphene crystal that is rolled up into a seamless cylinder. Fig. 1.1 (f) is typical STM image which reveals the atomic arrangement of the sp 2 rings within one single carbon nanotube. SWNTs are predicted to be metallic or semiconductors [START_REF] Mintmire | Are fullerene tubules metallic?[END_REF][START_REF] Saito | Electronic structure of graphene tubules based on C 60[END_REF] depending on their diameter and the helicity of the sp 2 rings arrangement in their walls (zigzag, armchair or chiral, see Fig. 1.1 (g)). This dependence on the diameter and the helicity, which was confirmed in many experimental works [START_REF] Odom | Atomic structure and electronic properties of single-walled carbon nanotubes[END_REF], has opened a large area of possible applications in electronics and optics (for further reading, see the book of R. Saito, G. Dresselhaus and M.S. Dresselhaus [36] on carbon nanotubes).

In the rest of this chapter, we will focus on graphene and multilayer graphene systems, presenting firstly their electronic and vibrational properties.

Properties of graphene and few-layer graphene systems

In this section, the electronic and vibrational properties of graphene and FLG systems are discussed using Density Functional Theory (DFT) and Density Functional perturbation Theory (DFPT). For the sake of simplicity, only fundamental properties are presented in this chapter, leaving more detailed study (such as the effect of electrons correlation) to be presented in the subsequent chapters. The understanding of a key Chapter 1. Raman spectroscopy property, called Kohn anomaly, is essential to the broader understanding of the Raman spectra in FLG systems and will be presented in an independent section.

The adiabatic approximation

Let us consider a crystal formed by electrons and nuclei interacting via the Coulomb law. The time scale associated with the motion of nuclei is usually much smaller than the time scale associated with electrons due to the high mass ratio M proton /m e ∼ 1836. This difference in time scale suggests that the quantum treatments of electrons and nuclei are different. Indeed, this observation is taken into account within the adiabatic approximation [START_REF] Born | Zur quantentheorie der molekeln[END_REF], where we assume that the nuclei motion is decoupled from the electrons system. In the adiabatic approximation, the electrons are assumed to be capable of instantaneously adjusting their motion to the relatively slow motion of the ions. This can be cast in a formal mathematical way by assuming that the total wave function Ψ of the system can be written as a product of two wave functions containing the electronic and nuclear degrees of freedom separately [38] :

Ψ(r, R) = Ψ e (r, [R]) Ψ i (R) (1.1)
where Ψ(r, R) is an eigenstate of the total Hamiltonian H given by :

H = i -2 2m ∇ r i + I -2 2M I ∇ R I + e 2 2 i =j 1 |r i -r j | + e 2 2 I =J Z I Z J |R i -R j | - iI e 2 Z I |r i -R I | (1.2)
with the eigenvalue E such that:

H Ψ(r, R) = E Ψ(r, R) (1.
3)

The substitution of Eq. 1.1 into Eq. 1.3 gives rise to two decoupled electronic and nuclear Schrödinger equations :

H e [R] Ψ e (r, [R]) = E e Ψ e (r, [R]) (1.4) H p (E e ) Ψ i (R) = E i Ψ i (R) (1.5)
In these expressions r and R are the electron and ion position vectors respectively, namely : R = {R 1 , .., R I , .., R N i } and r = {r 1 , .., r i , .., r Ne }, where N i and N e are the number of ions and electrons in the system.the quantities m and M I are the electron and nuclei masses respectively. The dependence of the electronic Hamiltonian H e [R] and the wave function Ψ e (r, [R]) on the atomic positions is just parametric and it is indicated by square brackets. A full proof of the adiabatic Approximation which justifies the validity of the previous system of equations is beyond the scope of the present thesis and will not be presented here (it can be found in Ref. [START_REF] Born | Dynamical theory of crystal lattices[END_REF] for example).

In the next chapter, we give some hints of the full proof. This approximation was successfully applied to the majority of materials including insulators metals and semimetals. However, the full proof does not give any guarantee of the validity of the adiabatic approximation on metallic systems. In a metal (or a semimetal), the energy difference between excited electronic states in Eq. 1.4 is zero and the electrons may, in principle, jump between excited states while being perturbed by the nuclear motion. This means that the two systems can no longer be considered independent. Nowadays, the adiabatic approximation is always used to study the electronic and vibrational properties of metallic and semimetallic materials with an accuracy of a few percent and it will be used in our work to study graphene, few-layer graphene and graphite.

Electronic properties

Fig. 1.2 (a) shows a schematic for the honeycomb structure of graphene in which the unit cell has two atoms labeled A and B. The whole lattice can be regarded as the superposition of two sub-lattices (A and B sub-lattices) arising from translation symmetry operations acting on the A and B atoms in the unit cell. The overlap between the hybridized orbitals of the carbon atoms during the formation of graphene gives rise to a σ bonding in the graphene plane while the remaining π orbitals overlap gives rise to a π bonding. The σ bonds are mainly responsible for the mechanical properties of graphene and they ensure the in-plane cohesion of the whole crystal, while the π bonds are mostly responsible for optical and electronic properties of graphene.

Fig. 1.2 (c) shows the electronic structure of graphene calculated with density functional theory (DFT) in the LDA approximation 1 . The Fermi surface in graphene is reduced to two points : K and K , edges of the first Brillouin zone (BZ) (see Fig. 1.2 (b)). Due to the hexagonal symmetry of the graphene honeycomb lattice with two atoms per unit cell, the dispersion of the electronic bands in graphene near the Fermi level is found to be linear. This implies that the charge carriers in graphene are described by a Dirac-like Hamiltonian, with a Fermi velocity of v f = 10 6 m/s ∼ 0.003 c, where c is the light velocity. This characteristic makes graphene very interesting from a fundamental point of view as it offers a solid bridge between theoretical electrodynamics of massless fermions in high-energy physics and the practical low-energy condensed matter physics.

In few-layer graphene systems, in addition to the in-plan σ and π-bondings, we have the interaction between the π orbitals from different graphene sheets. This interaction makes the electronic bands dispersion in FLG different from that in graphene, as there are more conduction and valence bands in the vicinity of the Fermi level (see the case of bilayer in Fig. 1.2 (d)). Although the linear dispersion is lost in FLG, the semimetallic character is preserved and the Fermi surface is always reduced to the points K and K (See Fig. 1.2 (d) for the bilayer case). After the isolation of FLG (2004) it was soon realized that new optical and electronic features can be found within FLG. For two layers, gated bilayer graphene was shown to possess a widely tunable band gap (up to 250 meV) [START_REF] Zhang | Direct observation of a widely tunable bandgap in bilayer graphene[END_REF]. When there is more layers (N layers > 2), the possibility of having different stacking sequences is shown to be of crucial importance [START_REF] Bao | Stacking-dependent band gap and quantum transport in trilayer graphene[END_REF][START_REF] Bao | Stacking-dependent band gap and quantum transport in trilayer graphene[END_REF]. We will study the stacking order in FLG more rigorously when we deal with the application of Raman spectroscopy on stacking order determination.

Vibrational properties

The motion of the nuclei of given system determines its vibrational properties. The quantum treatment of nuclei is carried out, in general, within the adiabatic approximation. Let us consider a crystal structure belonging to a given Bravais lattice with N a atoms per unit cell. The positions of the atoms may be labeled by an index I; this contains the unit cell index l to which a given atom belongs and the position of the atom within that unit cell indexed by s; I ≡ {l, s}. The position of the I th atom is thus :

R I = R l + τ s + u I (1.6)
where R l is the position of the l th unit cell in the Bravais lattice; τ s is the equilibrium position of the atom in the unit cell, and u I = u l,s indicates a deviation from the equilibrium position. In order to access the nuclei eigenstates Ψ i (R), one has to solve the Schrödinger equation in Eq. 1.5 , namely :

H p (E e ) Ψ i (R) = E i Ψ i (R) (1.7)
where the nuclear Hamiltonian, within the adiabatic approximation, is given by :

H p = I -2 2M I ∇ 2 R I + E e (R) (1.8)
The quantity E e is the total energy of the electrons system calculated at nuclear configuration R. The resolution of Eq. 1.7 is usually performed within the harmonic approximation. In this approximation, the energy E e (R) is written as a Taylor series in u I up to the second order; which simplifies the Hamiltonian in Eq. 1.8 into2 (∇ R I = ∇ u I ):

H p = I -2 2M I ∇ 2 u I + 1 2 I,J ∂E e ∂u α I ∂u β J u α,β I,J =0 u α I u β J (1.9)
where u α I is the α th Cartesian component of u I . Usually, the thermal wavelength of nuclei in solids is small compared to their distance. This allows us, in general, to treat nuclei as classical particles which results in the so-called classical nuclei approximation. In this approximation, one solves the Newtonian equations of motion for the ions and not the full Schrödinger equation in Eq. 1.7. To do that, one shall use the classical version of the previous Hamiltonian, which describes a set of decoupled harmonic oscillators. The vibrational frequencies ω of these oscillators are determined as the solution of the following linear system :

det 1 √ M I M J ∂E e ∂u α I ∂u β J -ω 2 .1 = 0 (1.10)
where 1 is the identity matrix. The derivative of the energy with respect to atomic displacement in the previous equation is called the inter-atomic force constants (IFC), usually denoted by C α,β s,t (R l , R m ). Due to translation invariance, the IFC is a function of R n = R l -R m . The linear system in Eq. 1.10 is infinite and its resolution is usually carried out by the Fourier transform technique [START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF]. Let us define the mass scaled Fourier transform of the IFC :

D α,β s,t (q) = Rn 1 √ M s M t e -iq.Rn C α,β s,t (R n ) (1.11)
where q is a vector belonging the first BZ in reciprocal space. The quantities ω 2 are now defined for each wave vector q as the solution of the following linear system : The subscripts s, t in Eq. 1.12 run over the atoms within the unit cell. The quantity

det D αβ st (q) -ω 2 (q).1 = 0 (1.
D αβ st (q) is called the dynamical matrix. It is a second rank tensor of dimension (3 × N a ) × (3 × N a ) (N a
is the number of atoms per unit cell) with eigenvalues ω 2 qν and eigenvector e qν representing the squared-frequency and the polarization of the decoupled harmonic oscillators of the system in Eq. 1.9 now called phonons. The displacement u I of an atom in the unit cell R l at a position τ s is now given by :

u α l,s = 1 √ M s q e s,α qν e iq.R l (1.13)
where M s is the mass of the considered atom. Within DFT, one can show that the dynamical matrix (or the IFC) is a functional of the first derivative of the electronic density at the equilibrium nuclear positions (see Ref. [START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF]). This fundamental result, first stated by De Cicco and Johnson (1969) and by Pick, Cohen, and Martin (1970) [START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF], transform the nuclear problem in Eq. 1.10 into a simpler one where we just need to calculate (self-consistently) the change of the electronic density with respect to a phonon perturbation. This type of calculation, is usually carried out using density functional perturbation theory (DFPT) [START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF] which constitutes our method to study vibrational properties.

An atomic displacement generated by a given phonon mode may or not be invariant under the action of a crystal symmetry. The set of symmetry operations under which an atomic displacement is invariant, is called the symmetry group of the phonon mode that generate this atomic displacement. Clearly, this symmetry group is a sub-group of the crystal space group. It is important to determine the symmetry properties of phonons in a crystal, since they determine whether or not a phonon mode is Raman active. This kind of analysis is carried out using group theory. A full discussion of phonon symmetries is beyond the scope of the present work, but the reader can find good references in the domain (see for example Ref. [43]). Table 1.1 recapitulates commonly used notations in the literature to label the phonon modes in crystals with respect to their symmetry properties. 3) . The "i" and "o" prefixes in labeling the phonon modes stand for "in-plane" and "out-of-plane" vibrations respectively. (c) Schematics for the Γ-point phonon displacement pattern for graphene and bilayer graphene, where each phonon mode in graphene gives to two phonon modes in bilayer. The Raman (R) and infrared (IR) activity is also indicated. Figure (c) is adapted from Ref. [START_REF] Reich | Raman spectroscopy of graphite[END_REF].

Properties of graphene and few-layer graphene systems

We will now introduce now some useful definitions regarding the type of phonon looking at it independently from its symmetry properties. A phonon is called longitudinal (transverse) if the corresponding displacement of the atoms u I is parallel (perpendicular) to the phonon wave vector q. Acoustic (optical) phonons are those for which the displacement of atoms in the unit cell is in-phase (out of phase). If the atomic displacement of the atoms following a given phonon mode is in the plane (out of plane) in a layered material, the phonon is called an in-plane (out of plane) vibration.

Applying the above theory to the case of graphene gives the phonon bands dispersion presented in Fig. 1.3 (a). Graphene has two atoms per unit cell, so we have in total 6 phonon modes, namely 3 optical and 3 acoustic modes. At the Γ point, they are reduced to 4 independent modes due to symmetry. These modes are depicted in Fig. 1.3 (c) with their symmetry properties.

Usually, in crystals with a center of symmetry (e.g. graphene), phonons that have the index "g" (from the German "gerade") are symmetric with respect to inversion and are Raman active (see Fig. 1.3 (c)). In graphene, at the Γ point, there are only two phonon modes that are Raman active. The first Raman mode is the E 2g (see Fig. 1.3 (c)) mode which is a doubly degenerate in-plane optical vibration. The second mode is the B mode (Fig. 1.3 (c)). Although the B mode is Raman active, this mode is undetected in Raman scattering experiments because of its weak electron-phonon coupling which can be put down to its out-of-plane vibration character (see Sec. 1.3).

Few-layer graphene have almost the same phonon band dispersion as graphene (see the case of bilayer in Fig. 1.3). Although phonons symmetries are different, we still have the E 2g vibration at the Γ point in the BZ which is Raman active. Furthermore, in FLG, a set of new phonon modes which are not present in the case of single layer graphene can be found. These modes are related to layers displacement parallel or perpendicular to the stacking direction. Namely, the Layer Breathing Modes (LBM) and the shear mode (see Fig. 1.3 (c)). These phonons have weak vibrational frequencies because of the weak Van Der Waals interaction between the graphene plan and this makes them hard to detect in Raman or neutron scattering experiments.

Kohn anomalies

Kohn anomalies (KA) [START_REF] Kohn | Image of the fermi surface in the vibration spectrum of a metal[END_REF] in metallic and semimetallic systems are abrupt softening of the phonon frequencies for phonons with wave vector q which connect two points on the Fermi surface. To explain why these anomalies occur at these particular phonon wave vectors, let us first remind the expression of the dressed phonon frequency according to perturbation theory when electron-phonon 3 interaction is taken into account. We have [START_REF] Grimvall | The electron-phonon interaction in metals[END_REF] :

ω 2 qν = Ω 2 qν -2 Ω qν Re[Π(q, ω qν )] (1.14)
where ω qν and Ω qν are the dressed and the bare phonon frequencies respectively. Re is the real part and Π is the phonon self-energy. The phonon self-energy at the lowest order in the perturbation series is given by :

Π(q, ω qν ) = 1 N k,n,m |g ν k n,k+q m | 2 (f k n -f k+q m ) k n -k+q m + ω qν + iδ (1.15)
where N is the number of k-points used to evaluate the previous sum and g is the electron-phonon coupling (to be introduced in the next chapter). The factors f refer to the Fermi-Dirac distribution function. In the denominator, k n denotes the electronic state energies (labeled by n, m) and δ is the phonon life-time due to electron-phonon interaction. As we can see, the above expression has to be solved in a self-consistent way because of the appearance of the phonon frequency ω qν in the left and right hand sides. However, the phonon energies are often one order of magnitude lower than the electronic ones and may be neglected. In general, two cases can be distinguished:

3. The electron-phonon interaction is presented in the next chapter. Phonon frequency (cm -1

) q = ( ,0,0) q = (0, ,0) GE a exp GI a exp q = ( ,0,0) Phonon Wave Vector (2 /a units) Figure 1.4 -(Color online) Upper panel: Lines are the phonon dispersion of graphene (GE), calculated at the experimental and equilibrium lattice spacings (a exp = 2.46Å and a th = 2.248Å). Experimental data from Ref. [START_REF] Maultzsch | Phonon dispersion in graphite[END_REF]. The red straight lines at Γ and K are obtained from evaluating the phonon self-energy within DFT. The two lower panels correspond to the dotted window in the upper panel. Here, graphite (GI) computed frequencies are also shown. The points are theoretical frequencies obtained by direct calculation. A single GE band corresponds to two almost-degenerate GIbands. Adapted from Ref. [START_REF] Piscanec | Kohn anomalies and electron-phonon interactions in graphite[END_REF].

Static case

If the phonon frequencies are neglected in the denominators in Eq. 1.15, the effect of the electron-phonon interaction on the phonon frequencies is said to be static. This approach is called static in the sense that the obtained dressed frequencies can be recast in a perturbative approach where the phonon is treated as a time independent perturbation. In this case, the energy denominators go to zero for wave vector q such that k,n -k+q,m = 0 on the Fermi surface. This singular behavior gives rise to KA and a softening of the phonon frequency at wave vector q that connect two points on the Fermi surface.

In graphene, graphite and FLG, this condition is satisfied for q = Γ and q = K. The phonon dispersions of graphene and graphite, as measured experimentally, exhibit an anomalous behavior at these particular wave vectors (see Fig. 1.4). This anomalous behavior has been explained by the presence of KA in Ref. [START_REF] Piscanec | Kohn anomalies and electron-phonon interactions in graphite[END_REF]. In Fig. 1.4 (upper panel), the two sharp kinks in the phonon dispersion at Γ and K points indicate the two Kohn anomalies at K and Γ. Due to the presence of Kohn anomalies, the phonon frequencies in graphene, graphite and FLG undergo a softening as high as 80 cm -1 [START_REF] Piscanec | Kohn anomalies and electron-phonon interactions in graphite[END_REF].

Experimentally, as shown in Fig. 1.4, we observe that the Kohn anomaly in graphene and graphite is visible only for the E 2g phonon mode at Γ and the A 1 phonon mode at K. The other branches do not soften, eventhough, their phonon wave vector satisfies the Kohn anomaly condition. In fact, in Eq. 1.15 we see that the electron-phonon coupling enhances the KA. Since the E 2g mode at Γ or the A 1 mode at K possess the highest electron-phonon couplings, Kohn anomalies are visible only for these two branches. The other branches have an electron-phonon coupling which is one order of magnitude smaller than the A 1 phonon at K [START_REF] Piscanec | Kohn anomalies and electron-phonon interactions in graphite[END_REF] and do not show any softening.

Taking into account the existence of Kohn anomalies while studying the Raman spectrum of graphene, FLG, carbon nanotubes and graphite is of great importance. In fact, a redshift of the peaks that arise from phonons affected by the Kohn anomaly, is expected. In carbon nanotubes (see Sec. 1.1), while metallic nanotubes possess Kohn anomalies, semiconducting nanotubes do not. Thus, Raman spectroscopy can be used to distinguish between metallic and semiconducting carbon nanotubes [START_REF] Piscanec | Kohn anomalies and electron-phonon interactions in graphite[END_REF] simply by observing the occurrence the peak's redshift.

Dynamical case

If phonon frequencies are not neglected in the denominators of Eq. 1.15, the treatment of the electron-phonon coupling is said to be dynamic. We notice that dynamical effects move the Kohn anomaly with respect to the static case since the singularity in the denominators is now shifted due to the presence of the ω term. Doped graphene, is a spectacular system where the static approximation fails miserably [START_REF] Pisana | Breakdown of the adiabatic born-oppenheimer approximation in graphene[END_REF] 4 . In fact, the static phonon frequency of the E 2g mode in graphene at the BZ center does not show any dependence on the doping, in contradiction with experimental measurements (see Ref. [START_REF] Pisana | Breakdown of the adiabatic born-oppenheimer approximation in graphene[END_REF] and references therein). As revealed in Ref. [START_REF] Pisana | Breakdown of the adiabatic born-oppenheimer approximation in graphene[END_REF], the observed dependence on doping as well as the stiffening of the G Raman peak in doped graphene is a consequence of moving the KA away from the q = 0 case due to dynamical effects. Indeed, by evaluating the integral (or the sum over k, n, m) in the expression given the self-energy with linearized electronic bands (graphene) we find the following expression for the dressed E 2g phonon frequency at the Γ-point [START_REF] Pisana | Breakdown of the adiabatic born-oppenheimer approximation in graphene[END_REF]:

ω Γ,E 2g = Ω Γ,E 2g + α| f | (1.16)
where α is a constant : α = 4.39 × 10 -3 and f is the Fermi level shift due to doping. Thus, while the static approximation works very well for most metallic systems, dynamical effects are very important while studying doped graphene. We will return to discuss these effects and their implications on the Raman spectrum of doped graphene in Sec. 1.4.2.

Raman scattering

Inelastic light scattering is one of the most powerful tools used to study fundamental properties of materials. Raman scattering, a well known example, is associated with inelastic scattering of light by elementary excitations of the materials. Historically, Raman scattering has been associated with phonons, although other elementary 4. Carbon nanotubes also have important dynamical effects, see for example Ref. [START_REF] Piscanec | Optical phonons in carbon nanotubes: Kohn anomalies, peierls distortions, and dynamic effects[END_REF]. excitations can be involved (plasmons for example) which may give rise to interesting phenomena [START_REF] Peter | Fundamentals of semiconductors: physics and materials properties[END_REF]. In this work, only scattering by electrons and phonons will be considered. In Raman scattering, the energy difference between the incoming and outgoing light is equal to the sum of the absorbed or emitted phonon energies. This energy difference is usually called the Raman shift. A Raman experiment consists of measuring Raman spectra which give the intensity of the scattered light as a function of the Raman shift.

The electromagnetic field interacting with the material system, during the Raman scattering, can be seen as quanta of energy called photon. These photons are characterized by a wave vector k and frequency ω L as dictated by standard electrodynamics. The incoming radiation ω i L will scatter either elastically (

ω i L = ω o L ) or inelastically (ω i L = ω o L ),
where ω o L is the frequency of the scattered light. In the latter case, the material system is promoted to an excited state which, in our case, corresponds to a phonon. The phonon (of frequency ω qν ) may either be emitted (Stokes processes) or absorbed (anti-Stokes processes) by the material system. Energy conservation law implies : ω qν = ω i L -ω o L for Stokes processes and ω qν = ω o L -ω i L for anti-Stokes processes, see Fig. 1.5. If the laser energy ω i L is much smaller than the material band gap, Raman scattering is said to be non-resonant. In this regime Raman scattering can be formulated in terms of the polarizability tensor [START_REF] Knoll | Raman scattering of atomic vibrations in anharmonic potentials[END_REF][START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF][START_REF] Veithen | Nonlinear optical susceptibilities, raman efficiencies, and electro-optic tensors from first-principles density functional perturbation theory[END_REF]. This is the so-called macroscopic view of Raman scattering (or equivalently, the dielectric approach to Raman scattering). If, on the other hand, the laser energy is comparable to the material band gap, Raman scattering is said to be resonant. In this regime, the macroscopic view is no longer adequate and Raman scattering has to be formulated from a micro-Chapter 1. Raman spectroscopy scopic view in which Raman cross sections are calculated directly from perturbation theory [START_REF] Ganguly | Theory of lattice raman scattering in insulators[END_REF]. In the next two paragraphs a brief discussion of the two formulations is given. Since graphene, FLG, and graphite are semimetallic systems, Raman scattering in these materials is always resonant with any external laser energy. Hence, this work will be based on a microscopic view of Raman scattering.

Macroscopic view

In this section, Raman scattering is formulated within a simple classical macroscopic view which is relevant to introduce some important quantities and definitions such as the Raman tensor. Suppose a given dielectric material is interacting with an oscillating electromagnetic field :

E(r, t) = E(k i , ω i ) cos(k i .r -ω i t)
where k i and ω i are the momentum and frequency of the incident electromagnetic wave. This interaction induces a polarization P(r, t) = P(k i , ω i ) cos(k i .r -ω i .t) in the dielectric with the same space and time periodicity. The relation between the amplitudes of the incident radiation and the induced polarization is by definition given through the polarizability tensor. Namely, in the linear regime we have :

P α (k i , ω i ) = χ α,β (k i , ω i ) E β (k, ω i ) (1.17)
where α and β indicate the Cartesian directions. If the material system is at a finite temperature, the polarizability χ is subject to thermal fluctuations. In the presence of a phonon with polarization vector e qν and frequency ω qν (one single monochromatic phonon perturbation), the collective atomic displacement Q(r, t) associated with this phonon can be expressed as a plane wave:

Q(r, t) = Q(q, ω qν ) cos(q.r -ω qν t) (1.18)
where the vector Q(q, ω qν ) is the amplitude of the atomic vibrations. The components of this vector are the vibrational amplitudes of each atom in the unit cell. Normally, the atomic displacements are small compared to the lattice spacing so one can develop the polarizability tensor as Taylor expansion in Q(r, t). At the first order of this expansion one has :

χ(k, ω i ) = χ 0 (k, ω i ) + ν ∂χ ∂Q ν Q=0 Q ν (r, t) (1.19)
where χ 0 is the polarisability of the medium without fluctuations and Q ν are the components of the the vector Q. By replacing the previous equation in Eq. 1.17, the polarization reads :

P α = P 0 α + P ind α (1.20)
where P 0 α is the non-fluctuating part arising from χ 0 . The fluctuation dependent part P ind reads :

P ind α (r, t) = ν ∂χ α,β ∂Q ν Q=0 Q ν (q, ω ν ) E β (k i , ω i ) × cos(k i .r -ω i .t) × cos(q.r -ω qν t) (1.21)

Raman scattering

Using the trigonometric identity :

cos(a) × cos(b) = 1 2 [cos(a + b) + cos(a -b)] (1.22)
The fluctuating part of the polarisability is now given by [START_REF] Peter | Fundamentals of semiconductors: physics and materials properties[END_REF]:

P ind α (r, t) = ν 1 2 ∂χ α,β ∂Q ν Q=0 Q ν (q, ω qν )E β (k i , ω i ) × cos[(k i -q).r -(ω i -ω qν )t] + cos[(k i + q).r -(ω i + ω qν )t] (1.23)
The first term between brackets in the previous equation corresponds to Stokes scattering, while the second corresponds to anti-Stokes scattering. Within this classical description, Raman scattering is simply related to the dependence of the polarizability tensor on the ionic degrees of freedom. As well known, the classical description of Raman scattering does not provide the correct ratio between Stokes and anti-Stokes processes. In order to describe correctly this ratio, one needs to pass to the quantized version of the previous theory where the atomic displacements are treated as quantum operators.

Momentum conservation implies k s ± q = k i where k s is the momentum of the scattered light and the ± sign is for Stokes and anti-Stokes processes. Since |k i,s | is much smaller than the typical size of the BZ in crystalline systems, the phonon wave vector involved in the one phonon processes described above is restricted to zone center phonons. This constitutes the fundamental selection rule for Raman scattering. If the Taylor expansion in Eq. 1.19 is pushed to higher orders, Raman scattering will involve multi-phonon emission or absorption and the fundamental selection will impose that the sum of all the emitted or absorbed phonon wave vectors is zero. Let us denote the incident (scattered) light polarization by e i (e s ). According to the dipole radiation formula at long distances, the scattered intensity is given by :

I s ∝ e s . ∂ 2 P ind (t) ∂t 2 2 (1.24)
which can be rewritten as :

I s ∝ |e s .R.e i | 2 (1.25)
where :

R = ν ∂χ ∂Q ν Q=0 × Q ν (q, ω qν ) |Q(q, ω qν )| (1.26)
is called the Raman tensor.

More selection rules can be derived from symmetry properties of the Raman tensor which can be deduced from the combined symmetry of the involved phonon and the polarizability tensor. The simplest one would be for centrosymmetric systems where only symmetric phonons at BZ center are Raman active (the modes that have a subscript "g" as discussed in Sec. 1.2.3). The inversion symmetry in these systems implies that the Raman tensor must be invariant under inversion. Given that the tensor ∂χ/∂Q is symmetric under inversion, anti-symmetric phonons give rise to an anti-symmetric Raman tensor which is in contradiction with the inversion symmetry requirement.

Electron Phonon

Microscopic view

Let us assume the total system formed by the material and the electromagnetic field to be described by the total Hamiltonian H tot = H + H l + H el , where H is the Hamiltonian describing the material system, H l is the Hamiltonian describing the radiation field and the term H el represents the interaction between the two systems. The Hamiltonian H can be decomposed, in a first approximation, into an ionic and an electronic part according to the adiabatic approximation (see Sec. 1.2.1), H = H e + H p . Furthermore, the interaction between electrons and phonons is described by the Hamiltonian H ep (see the next chapter for an exact description of the electronphonon interaction). Under these conventions, Raman scattering with phonons can be see as a three steps process (see Fig. 1.6) : (a) first, an electron hole pair is created after the annihilation of an incoming photon via H el (b) second, the electron or the hole is diffused by emitting or absorbing a certain number of phonons n with a momentum q i=1,n such that n i=1 q i = 0 (c) the electron and hole meet again and recombine to emit the scattered photon via the Hamiltonian H el . Perturbation theory up to the n + 2 order (where the Hamiltonians H el and H el are treated as perturbations) can be used to describe such a process. In this section, the formula that gives the scattered Raman intensity as a function of electron and phonon properties is presented. Its detailed derivation is provided in the next chapter. The scattered intensity of n-order Raman scattering is given by the generalized Fermi golden rule [START_REF] Peter | Fundamentals of semiconductors: physics and materials properties[END_REF] :

I ∼ f s 0 ,...,sn i|H el |s 0 s 0 |H ep |s 1 . . . s n-1 |H ep |s n s n |H el |f (E i -E 0 + iγ)(E i -E 1 + iγ) . . . (E i -E n + iγ) 2 δ(E i -E f ) (1.27)
where i and f stand for initial and final states and E i , E f the corresponding total energies (material system and the radiation field). The kets |s i represent the total system intermediate states with energies E 0,..,n and γ is the sum of electron and hole lifetimes. In this work, we will be interested in Raman scattering such that the initial and final electronic states of the material system are assumed to be the electronic ground state. A more general Formula can be derived for Raman scattering with defects by including the defect perturbation matrix elements in Eq. 1.27 (See Ref. [START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF]).

In this thesis, we will focus on two-phonon Stokes Raman scattering where the ionic part of the material final state contains two phonons. The processes described by Eq. 1.27 are generally associated with lines which are much weaker than firstorder Raman lines, since they arise from higher order processes in the perturbation series. Graphene, FLG and graphite are notable exceptions. During the intermediate virtual transition the energy is not necessarily conserved and the three denominators of Eq.1.27 are generally different from zero. However, in graphene, FLG and graphite two or more of the denominators of Eq. 1.27 can be equal to zero simultaneously and the intermediate electronic states are real excited states of the material. In the literature, this is called resonance condition and double-resonant in the case of two phonons Raman scattering. Raman peaks associated to resonant processes have an intensity comparable to that of lowest-order processes and are called resonant peaks and double resonant in the case of two phonons scattering.

Raman spectroscopy in graphene

The example of graphene is very suitable to introduce the reader to the field of Raman spectroscopy in carbon-based materials in general and FLG and graphite in particular. It is also a good example where Raman spectroscopy has been proven to be very useful. In this section we will discuss the most important Raman peaks in graphene and how they are used to probe its electronic and vibrational properties.

The G peak at 1580 cm -1 (Fig. 1.7 (a) and (b)) is a first order peak, and is related to the bond stretching phonon (E 2g ) at the Γ points in the BZ. The Raman process giving rise to the G mode (see Fig. 1.7 (b)) can be described as follows. First, an electron is exited from the conduction band to an empty state in the valence band. Second, the exited electron will be scattered by a phonon from the BZ center within the same Dirac cone. Finally, the recombination process takes place with the emission of the scattered photon. Although it is a simple first order peak, it has been shown (see Ref. [START_REF] Basko | Calculation of the raman G peak intensity in monolayer graphene: role of ward identities[END_REF]) that the structure of the G peak involves electronic transitions from a wide region in the BZ. Thus, even if it is a first order peak, ab-initio simulations of this peak is quite challenging. Only recently, the G peak was finally reproduced from first-principles in Ref. [START_REF] Reichardt | Ab initio[END_REF].

The 2D Raman peak, which will be extensively studied in this work, is situated around 2700 cm -1 (Fig. 1.7 (a) and (c)). It is a two phonons Raman peak and it is double resonant. Now we will discuss the exact meaning of this double resonance character. For graphene (and graphite), if one tries to carry out a simple numerical evaluation of the Raman cross section in Eq. 1.27 by neglecting all the matrix elements in the numerators for a given laser energy, one finds that for a certain phonon wave vector q near the K point in the BZ, the scattered intensity is clearly enhanced. This idea was explored by C. Thomson and S. Reich in Ref. [START_REF] Thomsen | Double resonant raman scattering in graphite[END_REF] where they showed that this enhancement is a result of the double resonant character of two phonons peaks in graphene as discussed in the next paragraph.

The Raman process responsible for the 2D peak consists of (see Fig. 1.7 (c)): (a) the creation of an electron-hole pair near the Dirac point K at wave vector K + k i (b) the scattering by a phonon with momentum q of the electron or a the hole (c) the scattering by a phonon with momentum -q of the electron or a the hole (-q because of momentum conservation) (d) the recombination of the electron-hole pair and the emission of the scattered photon. If the phonon that is involved in the previous process is around the K point in the BZ, namely q = K + d where d is a small vector, from momentum conservation, the scattered electron will have a momentum of

K + k i + K + d ∼ 2K + k i ∼ K + k i .
Let us neglect the phonon energy for simplicity. The intermediate electronic states before the scattering (at momentum K + k i ) and after the scattering (at momentum K + k i ) are real excited states of the material (Fig. 1.7 (c)). In Eq. 1.27, this means that, at least, two denominators in the expression of I are simultaneously equal to zero. This kind of process is called double resonant Raman scattering. Due to this resonance, the intensity of the 2D peak is enhanced and is comparable to that of the first order G peak. According to this view, the double resonant condition acts as a selection rule on the phonon wave vector and only phonons around K (and K ) contribute to the Raman cross section of the 2D peak. The D peak (respectively D peak, see Fig. 1.7 (b)) at around 1350 cm -1 (1600 cm -1 resp.) (Fig. 1.7 (a) and (c)), characterizes the presence of defects. These peaks are absent in pristine graphene. The Raman process responsible for the D peak is the same as the one responsible for the 2D peak except that the electron (or the hole) is scattered by a defect and by a phonon (Fig. 1.7 (c)). Following the same argument as the one presented in the previous paragraph, one finds that the D peak is also double resonant and is assisted by phonons around the K and K points in the BZ.

Uniaxial strain

Uniaxial strain in graphene consists of compressing or stretching the whole crystal out of its equilibrium shape along a given direction. The modification of the carboncarbon bond due to strain has a direct consequence on the G Raman peak, as this latter is related to the E 2g phonon mode at Γ. While in unstrained graphene at the Γ-point, the E 2g phonon (1581 cm -1 ) is two-fold degenerate, strain introduces a lifting of this degeneracy, giving rise to two independent phonon modes [START_REF] Mohr | Two-dimensional electronic and vibrational band structure of uniaxially strained graphene from ab initio calculations[END_REF]. The Raman The splitting of the 2D Raman mode in strained graphene in the zigzag direction. The red and blue curves are for light polarizations parallel and perpendicular to strain direction respectively. Adapted from Refs. [START_REF] Mohiuddin | Uniaxial strain in graphene by raman spectroscopy: G peak splitting, grüneisen parameters, and sample orientation[END_REF][START_REF] Huang | Phonon softening and crystallographic orientation of strained graphene studied by raman spectroscopy[END_REF] signature of this splitting is a splitting of the G peak in graphene into two sub-peaks called G + and G -. Each sub-peak is related to one of the two E 2g strain-induced modes. This was observed and confirmed through the experiments reported in Refs. [START_REF] Mohiuddin | Uniaxial strain in graphene by raman spectroscopy: G peak splitting, grüneisen parameters, and sample orientation[END_REF][START_REF] Huang | Phonon softening and crystallographic orientation of strained graphene studied by raman spectroscopy[END_REF]. Fig. 1.8 (a) shows the red-shift and the splitting of the G mode under uniaxial strain in graphene. An important observation was made by authors of Refs. [START_REF] Mohiuddin | Uniaxial strain in graphene by raman spectroscopy: G peak splitting, grüneisen parameters, and sample orientation[END_REF][START_REF] Huang | Phonon softening and crystallographic orientation of strained graphene studied by raman spectroscopy[END_REF] concerning the polarization behavior of the G + and G -sub-peaks. As these two peaks are assisted by two orthogonal phonon modes in graphene, their polarization behavior is predicted to be different. Polarization study 5 confirms this prediction. Furthermore, the angle between the strain direction and the phonon polarization vectors can be easily extracted from the polarization behavior of the G + and G -peaks. On the other hand, the two phonon polarization vectors arising from the splitting of the E 2g mode are directly related to the crystallographic orientation of the sample. Therefore, the relative position between the strain direction and the phonon polarization vectors of the two strain-induced E 2g modes determine completely the orientation of the sample. Thus, Raman spectroscopy offers a simple and efficient way to determine the crystallographic orientation of the graphene sample which turns out to be crucial in a wide domain of applications such as in nanoribbons [START_REF] Han | Energy band-gap engineering of graphene nanoribbons[END_REF] and quantum dots [START_REF] Ponomarenko | Chaotic dirac billiard in graphene quantum dots[END_REF] technologies.

In the previous paragraph, we discussed the effect of strain on the G Raman mode of graphene and its potential applications. A more recent study (Ref. [START_REF] Huang | Probing strain-induced electronic structure change in graphene by raman spectroscopy[END_REF]) shows that the 2D peak is also sensitive to strain and gives a distinct approach to probe strained-graphene electronic structure. As reported in Ref. [START_REF] Huang | Probing strain-induced electronic structure change in graphene by raman spectroscopy[END_REF] excitation and polarization studies 6 of the Raman 2D mode in single layer graphene under uniaxial strain 5. Polarization study consists of tracing the behavior of the scattered light intensity with a given polarization as function of the incident light intensity and polarization.

6. Laser excitation study consists of tracing the behavior of the a Raman peak's position (and eventually its intensity) as a function of the incident laser energy.

provide a powerful tool to probe the two types of modifications of the low-energy electronic structure of strained graphene : the deformation of the Dirac cones and their displacement away from the K point. Fig. 1.8 (b) shows the splitting of the 2D mode of graphene under strain. This splitting is directly related to the displacement of the Dirac cones and has a simple explanation. In the the case of the armchair configuration, for example, the induced strain moves the Dirac cones either close to or away from each other. Raman scattering between close Dirac cones gives rise to a blue-shifted 2D peak (called 2D + ) while scattering between Dirac cones shifted apart gives rise to a red-shifted 2D peak (called 2D -). Thus, the frequency difference between the 2D + and 2D -peaks is directly related to how much the Dirac cones are shifted.

A more complicated aspect of the electronic structure of strained graphene is the deformation of the Dirac cones. Due to the fact that the uniaxial strain changes the length of the carbon-carbon bond in the strain direction, whilst leaving it unchanged in the other one, we would expect a lower Fermi velocity in the strained direction. This effect can also be probed by the excitation energy dependence of the 2D + and 2D -peaks. In Ref. [START_REF] Huang | Probing strain-induced electronic structure change in graphene by raman spectroscopy[END_REF], the authors show the different behavior of the 2D + and 2D -peaks under laser excitation study. The different dispersions of the 2D + and 2D - peaks position as a function of the excitation energy directly reflect the anisotropy of the Fermi velocity. As expected, the greater shift rate of the 2D + mode indicates a lower Fermi velocity for electrons with a wave vector parallel to the strain. This offers a very efficient way to quantify the ellipsoidal character of the Dirac cones in strained graphene using Raman spectroscopy.

Doping

Doping in graphene consists of shifting the Fermi level away from the Dirac points either by electric field techniques, electrochemical doping [START_REF] Das | Monitoring dopants by raman scattering in an electrochemically top-gated graphene transistor[END_REF] or more recently the socalled space charge doping [START_REF] Paradisi | Space charge induced electrostatic doping of two-dimensional materials: graphene as a case study[END_REF] where high doping levels ( charge excess > 10 14 cm -2 ) can be reached. Both the G peak and the 2D peaks are sensitive to the Fermi level shift in graphene. Essentially, doping has two major effects on the phonon frequencies. First, doping moves the Kohn anomaly in graphene away from q = 0 (zone center) phonon which implies a stiffening of the E 2g phonon at Γ [START_REF] Ferrari | Raman spectroscopy of graphene and graphite: disorder, electronphonon coupling, doping and nonadiabatic effects[END_REF][START_REF] Lazzeri | Nonadiabatic kohn anomaly in a doped graphene monolayer[END_REF]. As explained in Sec. 1.2.4, due to the failure of the static-adiabatic approximation in doped graphene [START_REF] Pisana | Breakdown of the adiabatic born-oppenheimer approximation in graphene[END_REF], this stiffening is sensitive to doping [START_REF] Pisana | Breakdown of the adiabatic born-oppenheimer approximation in graphene[END_REF]. Second, the charge excess, which results from doping, changes the lattice parameters of the graphene sheet and, by consequence, influences the phonon frequencies at BZ borders.

The G peak sensitivity to doping comes as a result of the first effect. Raman measurements in doped graphene show, indeed, that the G peak always stiffens and sharpens for both electron and hole doping (see Fig. 1.9 (a)). A much more subtle effect in the behavior of the G peak as a function of doping is the variation of the G peak intensity. This variation can be explained within the non-resonant model for the G peak proposed by Basko [START_REF] Basko | Calculation of the raman G peak intensity in monolayer graphene: role of ward identities[END_REF]. Due to the fact that the G peak involves electronic contributions from states at the vicinity of the Fermi level, the intensity of the G peak is sensitive to doping simply because of Pauli blocking 7 . In fact, by blocking 7. Pauli blocking refers to a situation where an electronic transition is forbidden because the final state is not empty and that electrons being fermions can not occupy the same state with the same (c) The I(2D)/I(G) peaks ratio as a function of doping. Notice the high sensitivity of this ratio to doping. The solid blue lines are the predicted non-adiabatic trends from Ref. [START_REF] Pisana | Breakdown of the adiabatic born-oppenheimer approximation in graphene[END_REF]. Adapted from Ref. [START_REF] Das | Monitoring dopants by raman scattering in an electrochemically top-gated graphene transistor[END_REF].

the electronic transitions at the vicinity of the Fermi level, the contribution of the remaining transitions changes because of the absence of interference with the blocked ones (for more details see Ref. [START_REF] Ferrari | Raman spectroscopy as a versatile tool for studying the properties of graphene[END_REF]).

The 2D peak is also sensitive to doping. However, it is not straightforward to understand how doping can affect the 2D peak since it is assisted by phonons in the BZ where non-adiabatic effects are expected to be negligible 8 . Nevertheless, the change in the cell parameters, due to doping, affects the frequency of the phonons responsible for spin.

8. Recall that the non-adiabatic effects carry the main influence of doping, see Sec. 1.2.4. These effects are important near the Γ and K points only. the 2D peak in doped graphene. As revealed experimentally and explained theoretically in Refs. [START_REF] Lazzeri | Nonadiabatic kohn anomaly in a doped graphene monolayer[END_REF][START_REF] Das | Monitoring dopants by raman scattering in an electrochemically top-gated graphene transistor[END_REF], high hole doping yields to a blue-shifted 2D peak, while high electron doping yields to a red-shifted one ( see Fig. 1.9 (b)).

The different behaviors of the G and 2D peaks in graphene, with respect to doping, allow us to fully monitor doping in graphene. The 2D peak can be used to discriminate between electron and hole doping and the intensity ratio I 2D /I G (see Fig. 1.9 (c)) can be used as an indicator of the amount of doping in the sample [START_REF] Ferrari | Raman spectroscopy of graphene and graphite: disorder, electronphonon coupling, doping and nonadiabatic effects[END_REF].

Edges, defects and disorder

The study of defects in graphene was simplified enormously by the use of Raman spectroscopy. there exist Raman peaks characteristic of defects exist (they are present only if the sample has defects) : the D peak around 1350 cm -1 and the D peak around 1600 cm -1 . This allows for the characterization of defected graphene through the use of Raman spectroscopy.

Several type of defects can be found in graphene. They can be regrouped in two big families, namely, point defects and extended defects. Examples of point defects include : vacancies, divacancies, Stones-Wales or charged impurities. On the other hand, extended defects are delocalized along a 1D region within the graphene sheet. Edges and grain boundaries are the most studied extended defects.

Point defects have been studied since the early days of Raman spectroscopy in graphite. Tuinstra and Koening [START_REF] Tuinstra | Raman spectrum of graphite[END_REF] proposed a formula that gives the dependence of I(D), the intensity of the D peak, on the defect concentration :

I(D) ∼ 1/L 2
D where L D is the average defect-defect distance [START_REF] Tuinstra | Raman spectrum of graphite[END_REF]. This relation was further studied and confirmed in the work of Cançado et. al. [START_REF] Cançado | Quantifying defects in graphene via raman spectroscopy at different excitation energies[END_REF]. However, in Ref. [START_REF] Cançado | Quantifying defects in graphene via raman spectroscopy at different excitation energies[END_REF], the authors show the existence of two distinct regimes. For low defect concentration (L D > 3 nm) the Tuinstra and Koening formula holds. This is completely understandable within an independent defect picture where the contributions arising from all the defects add up constructively and I(D) is simply proportional to defect concentration (1/L 2 D ). For high defect concentration (L D < 3 nm, see Fig. 1.10 (a)) the intensity of the D peak is shown to be inversely proportional to defect concentration and the Tuinstra and Koening formula no longer holds. This can be understood in a simple way. When the defect-defect distance reaches the limit of the electron-hole free path length before a phonon scattering (which is around ω D ∼ 3 nm), the Raman process that gives rise to the D peak will be less probable with an increasing defect concentration. While the scattering with a defect becomes more probable, the scattering with a phonon (which is essential for the appearance of the D peak) becomes less probable. This is due to the fact that, at high defect concentration, the electron-hole pair will generally recombine before the scattering with a phonon takes place [START_REF] Ferrari | Raman spectroscopy as a versatile tool for studying the properties of graphene[END_REF].

The ratio I(D)/I(G) was proposed as a good estimator for defect concentration in defected graphene [START_REF] Cançado | Quantifying defects in graphene via raman spectroscopy at different excitation energies[END_REF]. The existence of the previously discussed two regimes seems, at first, to be a limitation for the use of Raman spectroscopy, since for each I(D)/I(G) ratio, two values for defect concentration exist (see Fig. 1.10 (a)). However, the width dependence of the G peak with respect to defect concentration (see Fig. 1.10 (b)) makes it to possible to distinguish between the two regimes. More precisely, the G peak width is a monotone function of the defect concentration. It becomes larger for higher defect concentration.

Edges can be seen as extended defects because they break the translation symmetry (at least perpendicular to the edges orientation). Perfect edges can exist within two configurations: either armchair or zigzag. While the armchair edges activate the D and the D peak in the Raman spectrum of graphene, the zigzag edges have no D peak but always a D peak. Thus, Raman spectroscopy (more precisely the I(D)/I(G) ratio) can be used to determine the type of perfect edges.

Counting the number of layers

In this section, we will show how Raman spectroscopy can be used to determine the number of layers in few-layer graphene. The study presented here is relatively close to the main subject of the present thesis which concerns stacking order determination of graphene layers by mean of Raman spectroscopy.

In the early days of graphene, the number of layers in isolated graphene flakes was determined by either Atomic Force Microscopy (AFM) or optical absorption in selected substrates [START_REF] Novoselov | Two-dimensional atomic crystals[END_REF]. Until 2006, these two techniques (accompanied with STM and HRSTM) have been the only techniques used to identify monolayer and FLG. The complexity of the equipment required for these techniques was a major limitation in the research filed of FLG. On the other hand, Raman spectroscopy was already known to be one of the most available non-destructive characterization tools and it was worthy to investigate the Raman response from graphene samples with different number of layers. Indeed, this investigation was carried out by Ferrari et. al. who showed that the shape and width of the 2D Raman peak in graphene and FLG drastically change according to the number of layers [START_REF] Ferrari | Raman spectrum of graphene and graphene layers[END_REF].

In the case of graphene, the simple structure of the electronic bands (two Dirac cones at K and K ) implies the existence of one single resonant phonon mode according to the double resonant model [START_REF] Thomsen | Double resonant raman scattering in graphite[END_REF] (see Fig. 1.11 (c)) and consequently gives rise to a single Lorentzian 2D peak (Fig. 1.10 (a)). In bilayer graphene the electronic structure is more complicated, with two valence and two conductions bands near the Fermi level. This "double bands" structure gives rise to four resonant phonon modes and, as shown in Ref. [START_REF] Herziger | Two-dimensional analysis of the double-resonant 2D raman mode in bilayer graphene[END_REF], to a more complicated 2D Raman peak with four sub-peaks (see Fig. 1.11 (a) and (d)). For N layers > 3, the situation is more or less the same, with more complicated electronic bands near the Fermi level and more resonant phonons modes. The resulting 2D peaks are larger with more sub-peaks contributions (see Fig. 1.11(a)) [START_REF] Ferrari | Raman spectrum of graphene and graphene layers[END_REF].

Other Raman peaks can be used to probe the number of layers in FLG : the C shear mode peaks around 40 cm -1 , as shown in Fig. 1.11 (b), and the layer breathing mode peaks (LBM) (1700-1800 cm -1 ). However, these peaks are less intense than the 2D peak and hard to detect. The shear mode peak (see Sec. 1.2.3) is related to layer displacements parallel to each other. A direct evidence for single layer graphene would be the absence of the shear mode in the Raman spectrum. Though, we should not use the absence of a peak as a characterization tool (because one can never be sure why a peak is absent [START_REF] Ferrari | Raman spectroscopy as a versatile tool for studying the properties of graphene[END_REF]). For N layers > 2, the frequency of this mode evolves with the number of layers and can also be used to probe N layers . The LBM (see Sec. 1.2.3) corresponds to the displacement of layers along the stacking direction. These phonon vibrations are detected in the Raman spectrum between 1650 and 1800 cm -1 . Contrary to the 2D peak which is due to in plane vibrations and sensitive to the electronic structure change with the number of layers, these last two signatures (the C and the LBM peaks) are related to phonons that are sensitive to interlayer coupling; thus, they are a direct probe of stacking in FLG.

To summarize the last four sections, we can say that : Raman spectroscopy is a powerful characterization tool that can be used to probe electrons, phonons, defects and the structure of few-layer graphene. The reason for such a wide domain of applications is that Raman spectroscopy is sensitive to almost all of the fundamental properties of FLG systems.

We turn now to discuss a more recent domain of application of Raman spectroscopy, namely, stacking study in FLG. Whilst it is relatively easy to count the number of layers in FLG, stacking determination is far more difficult to study. In the next section, we will present a state-of-the-art of methods for studying stacking in FLG. Following that, the main subject of the present thesis will be introduced.

Stacking order determination in few-layer graphene

Following the discovery of graphene, one of the most exciting fields of research has been to explore the electronic and optical properties of graphene layers. It was soon realized that FLG systems possess completely different electronic structures and new features can be found. For more than two layers, three possible stacking orders can occur in FLG : the Bernal, rhombohedral and mixed stacking (see Fig. 1.12 (a)). In the Bernal stacking mode (also called AB stacking), half of the atoms in a given layer are positioned perpendicularly between neighboring atoms from the plans on top and below, while the other half are between the centers of the hexagons. On the contrary, in rhombohedral stacking (also called ABC stacking), half of the atoms in each layer are positioned directly above of their neighbors on the layer below, and below the centers of the hexagons from the layer above, the opposite occurs for the other half. Mixed stacking is a mixture of rhombohedral and Bernal sequences.

In the case of trilayer graphene, while Bernal stacked trilayer is always a semimetal, transport measurements on suspended ABC-trilayer show the occurrence of a 42 meV gap at the Fermi level [START_REF] Lee | Competition between spontaneous symmetry breaking and single-particle gaps in trilayer graphene[END_REF] that has been attributed to the occurrence of a magnetic ground states (non vanishing spin density) [START_REF] Lee | Competition between spontaneous symmetry breaking and single-particle gaps in trilayer graphene[END_REF][START_REF] Pamuk | Magnetic gap opening in rhombohedral-stacked multilayer graphene from first principles[END_REF]. This finding was the first indication of the importance of stacking in few-layer graphene. Another aspect of rhombohedralstacked multilayer graphene (RMG) is the existence of an extremely flat band in the electronic structure. This feature implies the existence of a highly correlated electrons system and the possible appearance of exotic states of matter such as magnetic order or room temperature superconductivity. The flat band in RMG is a surface state with an increasing extent in reciprocal space with an increasing number of layers [START_REF] Xiao | Density functional investigation of rhombohedral stacks of graphene: Topological surface states, nonlinear dielectric response, and bulk limit[END_REF][START_REF] Kopnina | Surface superconductivity in rhombohedral graphite[END_REF]. Evidences for this flat band were reported in flakes with 5 layers in Ref. [START_REF] Pierucci | Evidence for flat bands near the fermi level in epitaxial rhombohedral multilayer graphene[END_REF] through scanning tunneling spectroscopy and angle-resolved photoemission spectroscopy. More recently, Henni et. al. [START_REF] Henni | Rhombohedral multilayer graphene: A magneto-raman scattering study[END_REF] were able to isolate multilayer graphene flakes with ABC sequences exceeding 17 graphene sheets. These samples were tentatively attributed to ABC stacking via magneto-Raman measurements.

In view of the previous findings, stacking in FLG is now understood to be of crucial importance. In order to identify the stacking mode in FLG, the method of reference, has been, so far, infrared absorption [START_REF] Mak | The evolution of electronic structure in few-layer graphene revealed by optical spectroscopy[END_REF][START_REF] Mak | Electronic structure of few-layer graphene: experimental demonstration of strong dependence on stacking sequence[END_REF][START_REF] Li | Structure-dependent fano resonances in the infrared spectra of phonons in few-layer graphene[END_REF] 9 . Raman spectroscopy has been also reported to be a reliable method for distinguishing layer stacking in trilayer and tetralayer graphene [START_REF] Lui | Imaging stacking order in few-layer graphene[END_REF]. However, the reported study of Ref. [START_REF] Lui | Imaging stacking order in few-layer graphene[END_REF] were first based on infrared absorption as a reference method. This means that, in a first stage, light absorption experiments were performed in order to identify the stacking order by comparing the obtained absorption spectra to simulation results of the optical conductivity. Raman spectra were later attributed to Bernal and rhombohedral stackings accordingly.

Infrared absorption technique has several limitations such as the low spatial resolution and the high sensitive instrumentation that is required for the detection of the infrared (IR) signal [START_REF] Lui | Imaging stacking order in few-layer graphene[END_REF]. Moreover, for stacking determination, it is reliable only in the case of trilayer graphene where there are exactly two possible stacking modes : the Bernal and rhombohedral stackings (see Fig. 1.12 (b)). In this case, the light absorption can be identified easily as originating from of a Bernal or a rhombohedral stacking by a simple comparison to optical conductivity simulations. All other possible spectra can be seen as arising from large domains with a mixed stacking. This is possible since the infrared signal is not very localized in space (compared to laser sources used in Raman experiments). For tetralayer graphene the situation is different. Three possible stackings can be found for tetralayer graphene, namely, the ABAB, ABCA and the ABCB (which is equivalent to the ABAC stacking). The ABCB stacking is predicted by to be less stable [START_REF] Aoki | Dependence of band structures on stacking and field in layered graphene[END_REF]. However, it can exist under certain (T ,P ) conditions, or produced by means of other technique then exfoliation which was used in Ref. [START_REF] Aoki | Dependence of band structures on stacking and field in layered graphene[END_REF]. In this case one can not be sure, for example, whether the absorption spectra of the infrared light arise from an ABCB stacking or from a mixture of ABAB and ABCA stacking signals.

In view of the previous discussions, we understand that a powerful technique to determine the stacking in FLG is highly required. In the next section, we will show how Raman spectroscopy can be used to fulfill this requirement.

The role of Raman simulations in stacking order determination

A Raman spectroscopy simulation consists of providing the theoretical Raman spectra for a given system based on the knowledge of its electron and phonon properties.

9. Not to be confused with infrared spectroscopy. Here, we mean simply light absorption in the infrared domain. It can be carried out within either parametric models (e.g. tight-binding) or from first-principles. The theoretical frame work for this type of computation is perturbation theory. Raman spectroscopy, contrary to the IR absorption technique, is a simple and powerful tool as its spectral sensitivity reaches the theoretical Density Functional Theory (DFT) limit on the phonon frequencies (few cm -1 ). Furthermore, it can reach high spatial resolution nowadays and, as discussed in previous sections, many authors have reported its reliability in characterizing graphene and few-layer graphene system.

The purpose of the present thesis is to investigate theoretically the use Raman spectroscopy to determine the type of stacking in FLG and graphite systems. Raman spectra will be calculated from first principles for all possible stacking and then compared to experimental results in order to identify the stacking sequence. We expect this contribution to widen the use of Raman spectroscopy to include stacking determination in few-layer graphene and bulk graphites.

Of all the peaks present in the Raman spectrum of FLG (and also graphite), the 2D peak is expected to be more sensitive to stacking. This is the case at least for trilayer graphene, where the 2D peak has been shown to be useful in distinguishing between ABA and ABC stacking orders while the G peak is almost the same for the two stackings [START_REF] Lui | Imaging stacking order in few-layer graphene[END_REF]. The sensitivity of the double resonant 2D peak to stacking is analogous to its sensitivity to the number of layers. This similarity, which justify our use of the 2D peak to differentiate between different stackings, can be understood as follows. The double resonant condition verified by the 2D peak in FLG and graphite imposes a selection rule on the phonon wave vector [START_REF] Thomsen | Double resonant raman scattering in graphite[END_REF] and only phonons around K and K contribute to the Raman cross section. From a geometrical point of view, as shown in Fig. 1.13, if one tries to connect points in the Dirac cones near then K and K points using a given phonon wave vector, only few electronic states will be connected (these states are indicated by red lines in Fig. 1.13). This means that the electronic region contributing to the 2D peak is very narrow. Thus, it is possible, in principle, to probe tiny modifications in the electronic structure of few-layer graphene such as the ones induced by a change in the number of layers or the stacking order. If the electronic regions contributing to the Raman cross section were large, the 2D peak would not be very sensitive to small deformation of the bands.

In view of the previous discussions, simulation of the Raman 2D peak will be at the heart of this work. We will aim to give the theoretical Raman signature for all possible stacking in FLG and graphite.

Chapter 2

Raman simulations and challenges

In this chapter we will discuss the theoretical background of Raman scattering. First, the key ingredients : electron-phonon and electron-photon coupling will be discussed, and then a derivation of Raman intensities from perturbation theory will be presented. Second, the implications of electron correlations beyond the local density approximation (LDA) in graphene, few-layer graphene, and graphite will be discussed. Finally, challenges that are present in Raman cross section calculation will be studied. Our methodology to cope with these challenges, namely, the Wannier interpolation technique and Brillouin zone (BZ) reduction method will be introduced briefly. The full discussion of these two methods will be presented in the next chapter.

Theory of Raman scattering : An overview

The response of a given material to an external light excitation is of great importance since it allows us to trace back its fundamental properties. One of the most important theoretical developments in this field came in the 80's when Density Functional Perturbation Theory (DFPT) was being developed [START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF][START_REF] Veithen | Nonlinear optical susceptibilities, raman efficiencies, and electro-optic tensors from first-principles density functional perturbation theory[END_REF]. However, in what concerns Raman scattering, all these developments were carried out for the so-called non-resonant regimes i.e. for light excitations smaller than the material band gap. In such a regime, the inclusion of excited states is not very important since the electron will almost never visit these states during the scattering as they are far away from the valence band. On the other hand, if the light excitation energy matches an electronic transition (the so-called resonant regime) a precise description of excited states is required. This fundamental difference between resonant and non-resonant regimes is the main reason for which Raman simulations have been historically separated into resonant and non-resonant types (for a review on that see P. Cardona [START_REF] Peter | Fundamentals of semiconductors: physics and materials properties[END_REF] and Refs. [START_REF] Merlin | Overview of phonon raman scattering in solids[END_REF][START_REF] Profeta | Theory of resonant raman scattering of tetrahedral amorphous carbon[END_REF]). In the non-resonant regime, DFT approaches used to evaluate Raman cross sections [START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF][START_REF] Veithen | Nonlinear optical susceptibilities, raman efficiencies, and electro-optic tensors from first-principles density functional perturbation theory[END_REF] are mainly based on derivatives of the dielectric function with respect to phonon displacements (dielectric approach) [START_REF] Knoll | Raman scattering of atomic vibrations in anharmonic potentials[END_REF]. These derivatives are accessible, within DFT, via the 2n + 1 theorem [START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF][START_REF] Veithen | Nonlinear optical susceptibilities, raman efficiencies, and electro-optic tensors from first-principles density functional perturbation theory[END_REF]. The dielectric approach led to important contributions for a wide range of semiconductors and insulators (e.g. the WURM database [START_REF] Bobocioiu | The wurm project -a freely available web-based repository of computed physical data for minerals[END_REF][START_REF] Knoll | Raman scattering of atomic vibrations in anharmonic potentials[END_REF]). However, it is not suitable for graphene, FLG and graphite systems which 29 Chapter 2. Raman simulations and challenges are semimetallic and resonant with any external laser excitation. For this reason, our approach to calculate Raman intensities is based on a direct evaluation of Raman cross sections from perturbation theory and not on the dielectric approach.

The direct perturbative approach (which is presented below) and the dielectric approach are shown to be equivalent in the non-resonant regime through the socalled Placzek approximation [START_REF] Placzek | Handbuch der radiologie[END_REF][START_REF] Brüesch | Phonons: Theory and Experiments II: Experiments and Interpretation of Experimental Results[END_REF]. The Placzek approximation states that, within the Born-Oppenheimer approximation, and in the off-resonant regime with a nondegenerate ground state, the Raman cross section of n-order Raman scattering as derived from perturbation theory is proportional to the n th derivative of the dielectric function with respect to phonon displacements. For a modern review on the link between the two approaches see Ref. [START_REF] Gillet | Ab initio study of Raman and optical spectra of crystalline materials and their temperature dependence[END_REF]. In the next four sections of this chapter, we will present a derivation of the Raman cross section from perturbation theory (PT) following Refs. [START_REF] Basko | Calculation of the raman G peak intensity in monolayer graphene: role of ward identities[END_REF][START_REF] Reichardt | Raman spectroscopy of graphene[END_REF][START_REF] Ganguly | Theory of lattice raman scattering in insulators[END_REF]. The main steps are to :

1. Pass to second quantization for electrons and phonons (Sec. 2.1.1).

Write down the expression of the electron-phonon coupling (vertex) in second

quantization (Sec. 2.1.2).

3. Write down the expression of the electron-photon coupling (vertex) in second quantization (Sec. 2.1.3).

4. Write down the PT expression of the Raman cross section (Sec. 2.1.4).

Second quantization for electrons and phonons

Our starting point will be the general expression of the Hamiltonian H describing the material system considered to be a crystal :

H = i -2 2m ∇ r i + I -2 2M I ∇ R I + e 2 2 i =j 1 |r i -r j | + e 2 2 I =J Z I Z J |R i -R j | - iI e 2 Z I |r i -R I | (2.1)
where r and R are the electrons and ions position vector operators respectively, namely R = {R 1 , .., R I , .., R N i } and r = {r 1 , .., r i , .., r Ne }, where N i and N e are, respectively, the number of ions and electrons in the system. The last three terms in Eq. 2.1 describe the electron-electron, ion-ion, and electron-ion interactions and will be denoted by U ee , U ii and U ei respectively. The Schrödinger equation describing the material systems is

H Ψ(r, R) = E Ψ(r, R) (2.2)
In the adiabatic approximation, the state Ψ(R, r) is written as a product of electronic and ionic wave functions :

Ψ(r, R) = Ψ e (r, [R]) Ψ i (R) (2.
3)

The occurrence of square brackets in the previous equations indicate a parametric dependence according to the adiabatic approximation (see first chapter). For a given nuclear configuration R, the electronic wave function is defined as the solution of the electronic systems described by the Schrödinger equation :

H e [R] Ψ e (r, [R]) = E e ([R]) Ψ e (r, [R]) (2.4)
where

H e [R] = - 2 2m i ∇ r i + U ee + U ei [R] + U ii [R] (2.5)
In order to satisfy the full Schrödinger equation in Eq. 2.2, The ionic part must satisfy

I -2 2M I ∇ R I + E e ([R]) + ∆H(r, R) Ψ i (R) = E Ψ i (R) (2.6)
where ∆H is a term that depends on the operators r and R. It is possible to define a decoupled system of ions if the term ∆H in the previous equation is neglected. Indeed, in the adiabatic approximation, this term is proven to be small compared to H [START_REF] Grimvall | The electron-phonon interaction in metals[END_REF] and the system of ions is described by the following Schrödinger equation :

H p (E e ) Ψ i (R) = E Ψ i (R) (2.7)
where

H p (E e ) = I -2 2M I ∇ R I + E e ([R]) (2.8)
Upon the resolution of Eqs. 2.4 and 2.7, one can construct the full spectrum H from the two decoupled electronic and ionic systems1 . In the next paragraphs, we will give a second quantization formulation for electron and ions systems starting from Eqs. 2.4 and 2.7. This formulation is relevant for perturbation theory derivation of the Raman cross sections.

Electrons

In our approach, the electronic problem described by the Schrödinger equation in Eq. 2.4 is treated within DFT as formulated by Kohn and Sham [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF]. Within this formulation, instead of solving the complicated many-body problem in Eq. 2.4, we solve the following independent particles problem [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF]:

Ψ(r) = i ψ i (r i ) (2.9) n(r) = i |ψ i (r i )| 2 (2.10) H KS [n(r)] ψ i (r) = E i ψ i (r i ) (2.11)
where Ψ is the total wave function of the independent particles system. The functions ψ i (r i ) are the one particle states (the so-called Kohn-Sham orbitals) and n(r) is the associated density. The Kohn-Sham Hamiltonian in Eq. 2.11 is given by [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF] :

H KS [n(r)] = - 2 2m i ∇ r i + V scf KS [n(r)] (2.12)
where the effective interaction (called the self-consisted Kohn-Sham potential)

V scf KS [n(r)
] is given by2 :

V scf KS [n(r)] = e 2 2 ˆn(r ) |r -r | dr + U ei + U ii + H xc [n(r)] (2.13)
In the previous equation, the first term in the right hand side is the classical electrostatic self-interaction of the electron charge-density distribution (Hartree term). The term H xc is called the exchange-correlation term. This term is introduced into the Kohn-Sham Hamiltonian to enforces the independent particles system to have the same ground state density n(r) as the whole interacting system described by the full Schrödinger equation in Eq. 2.4. This effective interaction is unknown, however, within Kohn-Sham DFT (see the original paper of Kohn and Sham [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF]) it can be expressed in a good approximation as the effective interaction of a homogeneous electrons gas if we assume that, locally, the density n(r) is a slowly varying function. This is the socalled the Local Density Approximation (LDA) which will be used in our work. Upon the resolution of the independent particles system in Eq. 2.11, with an appropriate exchange correlation term (LDA in our case) one can calculate the ground state density and all the ground state properties of the full interacting system.

In a crystal, translation invariance implies that the one particle states of the electron (which are the Kohn-Sham orbitals ψ(r i ) within DFT) are the so-called Bloch states :

ψ(r i ) ≡ ψ e nkσ (r i ) = e ik.r i u nkσ (r i ) (2.14)
where u(r i ) is a crystal periodic function 3 . The electron crystal-momentum k belongs to the first BZ and n, σ are indexes for electronic bands and spin. It is important to stress that the Kohn-Sham energies E i ≡ E e knσ do not have a direct physical meaning, thought they do provide a good starting points for more sophisticated approaches to the electronic problem (many body approaches, see Sec. 2.2). As DFT is a single particle theory, following standard second quantization rules for fermions, the electronic Hamiltonian in Eq. 2.4 can be expressed as [START_REF] Grimvall | The electron-phonon interaction in metals[END_REF][START_REF] Jishi | Feynman diagram techniques in condensed matter physics[END_REF] :

H e = knσ E e knσ c † knσ c knσ (2.15)
where c † knσ and c knσ are the electron creation and annihilation operators and E e knσ are the Kohn-Sham energies.

Ions

Lattice vibrations were already treated in the first chapter (Sec. 1.2.3). Here we present a second quantization formulation, while keeping the same notations. Let us consider a crystal structure belonging to a given Bravais lattice with N unit cell. The positions of the atoms are labeled by an index I, that contains the unit cell index l to which a given atom belongs and the index s of its positions within that unit cell, namely, I ≡ {l, s}. The position of the I th atom is thus :

R I = R l + τ s + u l,s (2.16)
where R l is the position of the l th unit cell in the Bravais lattice, τ s is the equilibrium position of the atom in the unit cell, and u I ≡ u l,s indicates the deviation from the equilibrium position. Suppose that we solved the nuclear problem as discussed in the first chapter. Let ω qν and e qν be the frequency and polarization vector of the vibrational normal mode of ions. The Hamiltonian in Eq. 2.7 in the harmonic approximation can be written as4 :

H p = I P 2 u I 2M I + 1 2 Iα,Jβ C α,β I,J u α I u β J (2.17)
where P u I = -i ∇ u I and C α,β I,J is the inter-atomic force constants matrix. M I is the mass of the atom I and α, β indicate Cartesian directions. Following Ref. [START_REF] Ashcroft | Solid State Physics[END_REF], App. L, we define the operators :

b q,ν = 1 √ N i l,s,α e -iq.R l e s,α qν   M s ω qν 2 u α l,s + i 1 2 M s ω qν P u α l,s   b † q,ν = 1 √ N i l,s,α e -iq.R l e s,α qν   M s ω qν 2 u α l,s -i 1 2 M s ω qν P u α l,s  
One can show that b and b † have the usual properties of creation and annihilation operators for bosons. The ionic Hamiltonian H p in the new operators reads :

H p = qν ω qν 1 2 + b † qν b qν (2.18)
which represents a sum of independent harmonic oscillators (called phonons) with a frequency spectrum ω qν (1/2 + n), n being {0, 1, 2..} and q in the BZ. The action of the operators b and b † in the n th quantum state of these harmonic oscillators is given by :

b † |n = (n + 1) 1/2 |n + 1 (2.19) b |n = (n) 1/2 |n -1 (2.20)
Finally, the atomic displacement u l,s of an atom at unit cell l and position τ s , can be expressed in terms of creation and annihilation operators, namely one has [START_REF] Ashcroft | Solid State Physics[END_REF] :

u α l,s = q,ν 1 2N M τ ω qν 1/2 e iq.R l e s,α qν (b q,ν + b † -q,ν ) (2.21)

Electron-phonon coupling

Electron-phonon coupling is a key ingredient used to describe many phenomena in condensed matter physics such as : superconductivity, electronic transport, or Raman scattering. It arises from the possible coupling between electronic states via the effective interaction felt by the electron during the atomic vibrations. In order to describe this interaction, let us first remind ourselves of the so-called Bloch relation. Consider the matrix elements :

n e , n p | H n e , n p (2.22)
in which |n e , n p are the material states taken in the adiabatic approximation and labeled by n e and n p that are sets of quantum numbers which determine the electron and phonon states respectively. If we make the adiabatic approximation also for H i.e. if we neglect completely the ∆H term defined in Eq. 2.6, then the non-diagonal matrix elements in the previous equation are zero. However, if we retain the ∆H term in the action of H, assuming that it is a small contribution which does not break the adiabatic approximation, then 2 nd order perturbation theory tells us that the nondiagonal matrix elements n e , n p | ∆H n e , n p describe electronic scattering processes between different electronic states (n e ↔ n e ) accompanied with a change in the ionic states (n p ↔ n p ) (i.e. lattice vibration). Thus, the term ∆H describes the electronlattice interaction. If the lattice is described within the harmonic approximation, then ∆H represents the electron-phonon interaction i.e. ∆H ≡ H ep . The evaluation of the matrix elements of H ep can be simplified using the separability of eigenstates of H into electronic and ionic parts according to the adiabatic approximation. Up to the first order in u (the atomic deviation from equilibrium position), this yields to the so-called Bloch relation

n e , n p | H ep n e , n p = I n p | u I n p • n e | W (r, R 0 I ) n e (2.23)
where R 0 is the equilibrium position of the atoms. Thus, in the decoupled electron and phonon states, H ep takes the form :

H ep (r, R) = I u I • W (r, R 0 I ) (2.24)
where the dots indicates a scalar product. In these expressions, the quantity W defines the action of H ep on the electronic part of the total wave function |n e , n p . The sum is extended over all the atoms in the crystal indexed by I and u I is a small deviation from the equilibrium position of the atom I which is denoted by R 0 I . The interaction W is given by (see Ref. [START_REF] Grimvall | The electron-phonon interaction in metals[END_REF] for a detailed derivation) :

u I • W (r, R 0 I ) = U e (r, R 0 I + u I ) -U e (r, R 0 I ) (2.25)
where U e (r, R 0 I ) is the total effective interaction felt by an electron at position r during the atomic vibration around the equilibrium positions R 0 I . The phonon part in the Bloch relation is immediately obtained from the relation in Eq. 2.21 which gives the atomic displacement in terms of phonon creation and annihilation operators. All the difficulty in the electron-phonon calculation lies in the electronic part (i.e. W ). In order to evaluate these matrix elements, we will describe two methods : the rigid ions approximation and the first principle self-consistent methods.

Rigid ions approximation

In this approximation, we assume that the interaction U e follows rigidly the motion of the ions. For small displacements u I , one can write :

U e (r, R 0 I + u I ) -U e (r, R 0 I ) = U e (r -R 0 I -u I ) -U e (r -R 0 I ) = u I .∇ R I U e (r -R I ) = -u I .∇ r U e (r -R I ) (2.26)
In the rigid approximation, the term U e (r -R I ) in the previous equations is simply replaced by the bare electron-ion potential V e-i :

U e (r, R 0 I + u I ) -U e (r, R 0 I ) = -u I .∇ r V e-i (r -R I ) (2.27)
A comparison between the previous equation and Eq. 2.25 leads to the following definition of the interaction W within the rigid ion approximation:

W (r, R 0 I ) = -∇ r V e-i (r -R 0 I ) (2.28)
This approximation was widely used in the early days of electron-phonon interaction calculation due to its simplicity. However, this method suffered from the divergence at q = 0 in the Fourier transform of V e-i . Besides this problem, it was already known that the application of the rigid ions approximation is very limited since it neglects the so-called electronic screening of the electron-ion interaction as it will be discussed in the next paragraph.

First principle self-consistent method

The major disadvantage of the rigid ions approximation is to assume that the response of the electron systems to an atomic vibration is only carried out via the bare electron-ion interaction. To see this disadvantage, within DFT for example, the electrons may also respond to any atomic vibration via the electron-electron interaction term because of its dependence on the electron density which is affected by the atomic vibration. This effect is called the electronic screening. As a remedy to the lack in the rigid-ion description, one must include the electronic screening explicitly .i.e. by calculating the self-consistent potential felt by the electron due to the atomic vibration. This approach gave rise to the so-called self-consistent method (used for the first time by Dacorogna [START_REF] Dacorogna | Self-consistent calculation of the q dependence of the electron-phonon coupling in aluminum[END_REF] in 1985 within a frozen phonon approach). Within DFT, this is done by calculating the derivative with respect to the electron position r of the whole

Chapter 2. Raman simulations and challenges

Khon-Sham potential and not only the electron-ion potential as done in Eq. 2.28, namely :

W (r, R 0 I ) = -∇ r V scf KS (r -R 0 I ) (2.29)
where V scf KS is the Kohn-Sham self-consistent potential already introduced in Eq. 2.13. The arrival of Density Functional Perturbation Theory (DFPT [START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF][START_REF] Veithen | Nonlinear optical susceptibilities, raman efficiencies, and electro-optic tensors from first-principles density functional perturbation theory[END_REF]) simplified further the calculation of the electron-phonon coupling. Nowadays, first principle selfconsistent calculation within DFPT is the standard technique used to calculate electronphonon coupling (hence the name first-principles self-consistent method). In our work, we will make use of this method as implemented in the Quantum ESPRESSO package to calculate the electron-phonon coupling. Let us now present a second quantization version of the first-principles self-consistent method calculation of the electronphonon coupling following Ref. [START_REF] Giustino | Electron-phonon interactions from first principles[END_REF].

It is easy to construct a second quantization version of the electron-phonon coupling from Eq. 2.24 by treating the operators u and W as one particle operators (for phonon and electron respectively). Hence, combining Eq. 2.24 with Eqs. 2.29 and Eq. 2.21, gives the following form for the electron-phonon Hamiltonian H ep [START_REF] Giustino | Electron-phonon interactions from first principles[END_REF] :

H ep = 1 √ N k,q,mnν g ν mn (k, q) c † m,k+q c n,k (b † -qν + b qν ) (2.30)
Here, k and q are the electron and phonon momentum respectively and n, m are electronic band indexes. The index ν stands for the phonon branch, c, c † (b, b † ) are the electron (phonon) creation/destruction operators previously introduced, and N is the number of unit cells in the crystal. The constants g are called electron-phonon matrix elements which, in DFPT, are given by [START_REF] Giustino | Electron-phonon interactions from first principles[END_REF] :

g ν mn (k, q) = u m,k+q ∆v scf,ν KS (q) u n,k (2.31)
where u n,k is the periodic part of the electron Bloch function and ∆v scf,ν KS (q) is given by :

∆v scf,ν KS (q) = s,α 2M s ω q,ν e ν s,α (q) ∂ s,α,q v scf KS ∂ τ,α,q v scf KS = l e iq•(r-R l ) ∂V scf KS ∂u α l,s R 0 (2.32)
Here V scf KS is the self-consistent Kohn-Sham potential and e ν (q) is the eigenvector of the dynamical matrix at phonon momentum q. The index l is a unit cell index and τ runs over the possible position vectors of atoms within the unit cell. Finally, α denotes a Cartesian direction.

Electron-photon coupling

The electromagnetic field interacting with the material system is quantized following standard rules of electrodynamics. We will use the so-called radiation gauge or the Coulomb gauge in which the scalar potential is set to zero and the vector potential A obeys ∇.A = 0. Thus, the electromagnetic field is entirely described by the vector potential A. The Hamiltonian of the electromagnetic field H l in the second quantized form is given by [START_REF] Jishi | Feynman diagram techniques in condensed matter physics[END_REF] :

H l = Q,λ Ω λ (Q) 1 2 + a † λ (Q)a λ (Q) (2.33)
where a † λ (Q) and a λ (Q) are photon creation and annihilation operators. The index λ runs over the 2 possible polarizations of the photon and the quantities Q and Ω λ (Q) define its momentum and frequency. The Hamiltonian H l also represents a set of independent harmonic oscillators analogous to the phonon case described previously, with a different dispersion relation. Namely, for light we have : Ω λ (Q) = c |Q| where c is the light velocity. The eigenstates of these harmonic oscillators (photons) are labeled by a positive integer m, with eigenenergies Ω λ (Q)(1/2 + m). The action of the operators a and a † is given by :

a † |m = (m + 1) 1/2 |m + 1 (2.34) a |m = (m) 1/2 |m -1 (2.35)
The vector potential can be expressed in terms of creation and annihilation operators as:

A(r) = 4π c 2 2Ω(Q)V 1/2 ˆd3 Q (2π 3 ) λ Q,λ a(Q, λ) + -Q,λ a † (-Q, λ) e iQ.r (2.36)
where V is the volume in which we quantize the electromagnetic field and are the polarization vectors of the photon. In order to express the electron-photon coupling in quantum mechanics, the minimal coupling prescription must be used starting from the classical expression. Namely, the Hamiltonian describing the electron moving in an electromagnetic field is given by 5 [85]

H = H[(p → p + e c A)] = H + e m c A(r).p + e 2 2 mc 2 [A(r)] 2 = H + H el (2.37)
where p is the electron momentum and e, m are the electron charge and mass respectively. The quantum mechanic version of the electron-photon coupling can be obtained now by promoting p and A in the expression above into operators. Usually, the term proportional to A 2 in the expression giving H el involves two-photon scattering processes. For weak fields, this term is generally far less important than the term proportional to A which involves single-photon processes; henceforth, the A 2 term will be ignored. Giving the expression of A(r) in Eq. 2.36, one can use Eq. 2.37 to express the electron-photon coupling Hamiltonian in second quantization [START_REF] Jishi | Feynman diagram techniques in condensed matter physics[END_REF] :

H el = i e cm A(r i ).p i = knm Q,λ γ nm,λ k,k+Q c mk+Q c † n,k (a Q,λ + a † -Q,λ ) (2.38)
5. We work with Gaussian (or cgs) units. In SI units, the factor e/c is to be replaced by e.

where :

γ nm,λ k,k+Q = 2 4π c 2 2Ω(Q)V 1/2 mk + Q| e iQ•r p. λ (Q) |nk (2.39)
where the factor 2 accounts for spin. a and a † are the usual electron creation and annihilation operator and k n, m are the electron momentum and band indexes respectively. The γ coefficients in the expression above are called the elctron-photon coupling matrix elements. They play the same role as the coupling coefficient g in the case of phonons (see Sec. 2.1.2).

The expression that gives the γ coefficients in Eq. 2.40 is often simplified further by the dipole approximation in which we replace the exponential in Eq. 2.40 by 1, namely : exp(iQ • r) ≈ 1. This approximation is justified by the fact that the wave length of the light (∼300-800 nm for light in the visible spectrum) is much larger than the typical crystal length scale given by the lattice constant (2.46 Å for graphene), i.e. Q • r ≈ 0. In other words, for the description of Raman spectra, we can treat the external photons as having a zero wave vector (we set Q = 0) yet a finite frequency. Thus, only vertical transitions are allowed, and the electron-photon coupling takes the form :

γ nm,λ k = 2 4π c 2 2Ω(Q)V 1/2 λ . mk| p |nk (2.40)

Second order Raman scattering

In this section, we assume that the reader is familiar with perturbation theory technique. The total Hamiltonian describing the material system and the radiation field is written as :

H tot = H + H l + H el (2.41)
where H is the Hamiltonian describing the material system assumed to be a crystal, H l is the radiation field Hamiltonian and H el is the Hamiltonian describing the interaction between the material system and the electromagnetic field. Furthermore, the Hamiltonian H can be decomposed into an ionic and electronic contribution within the adiabatic approximation (H p and H e respectively). We also must add a term which represents the interaction between the electrons and lattice vibrations. Namely, we have :

H = H e + H p + H ep (2.42)
In our approach, we assume that the treatment of the ionic part is done within the harmonic approximation, which maps the ionic problem onto a system of independent phonons. The interaction term H ep represents the electron phonon interaction which was treated in Sec. 2.1.2. Within these conventions one can rewrite H tot in the following form :

H tot = (H e + H p + H l ) + (H el + H ep ) = H 0 + H 1 (2.43)
where

H 0 = H e + H p + H l (2.44
)

H 1 = H el + H ep (2.45)

Theory of Raman scattering : An overview 39

The reason for this decomposition is to separate the contribution to H tot which will be treated as perturbation for the derivation of the Raman cross section. This contribution is the sum of electron-phonon and electron-light interactions (H 1 ). Second order Raman scattering which is discussed here corresponds to a scattering where, in the initial state, the material system is assumed to be in its ground state.

In the final state, the electronic part of the material system is always assumed to be the electronic ground states, whereas the phonon part contains 2 emitted phonons. Within these assumptions, one can carry out the calculation of second order Stokes Raman scattering via the generalized Fermi golden rule applied to the perturbation H 1 . Namely, the probability that an incident photon with the polarization λ in and frequency Ω in is scattered into an element of the solid angle dΩ s with the polarization λ out and with the frequency Ω out after the emission of 2 phonons ω qν and ω -qµ (-q for momentum conservation) is given by [START_REF] Reichardt | Ab initio[END_REF][START_REF] Basko | Theory of resonant multiphonon raman scattering in graphene[END_REF][START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF]]

dI dΩ s ∝ Ω 2 out (2π) 2 c 4 f |M if | 2 × δ Ω in -Ω out -(ω qν + ω -qµ ) (2.46)
where the i and f indexes stand for initial and final states. The quantity M if is the scattering matrix which describes the transition rate of the total system from the initial state to the final state. The evaluation of the scattering matrix elements M is relatively simple for first order Raman scattering. For second and higher order Raman scattering, the evaluation of the M matrix is tedious and usually carried out using Feynman diagrams technique.

For second order Raman scattering, one finds 3 independent Feynman diagrams 6 which are represented in Fig. 2.1. In literature, most standard textbooks [START_REF] Peter | Fundamentals of semiconductors: physics and materials properties[END_REF] prefer to work with the so-called Goldstone diagrams instead of Feynman diagrams. Goldstone diagrams are Feynman diagrams that are time ordered .i.e. the time ordering of the vertexes matters in counting the total number of diagrams. This means that two topologically equivalent Feynman diagrams which have different time ordering of the vertexes are considered as two independent Goldstone diagrams. Expressing the matrix M in terms of Feynman or Goldstone diagrams is completely equivalent. In our approach, we choose to work with Goldstone diagrams. For second order Raman scattering, the 3 possible Feynman diagrams lead to 12 Goldstone diagrams where only 8 can be resonant 7 . The non-resonant diagrams give negligible contributions to the Raman cross section compared to the resonant ones, thus, only resonant Goldstone diagrams will be retained in our theory. These diagrams are represented in Fig. 2.1 . Now, we will make the sum over all possible Goldstone diagrams and intermediate states explicit in the expression of M. Let us drop the index i in M if since it always indicates the initial ground state. If the final state contains two phonons with 6. If we take into account higher order terms in the electron-phonon and electron-photon interactions (e.g. 2-phonon vertexes) one finds that there are 11 Feynman diagrams, see Ref. [START_REF] Reichardt | Raman spectroscopy of graphene[END_REF].

7. The number of Goldstone diagrams a Feynman diagram can generate is equal to the number of external lines attached to the internal loop of the diagram in the case where two lines do not meet in the same vertex, see Ref. [START_REF] Reichardt | Raman spectroscopy of graphene[END_REF].

Light:

Electron:

Hole: Phonon:

Feynamn diagrams momentum q and -q from branches µ and ν (f ≡ {µ, ν, q}), then :

M µ,ν (q) = 8 α=1 M α µ,ν (q) (2.47) M α µ,ν (q) = 1 N k ijlm k M α,ijlm µ,ν (k, q) (2.48)
The first equation represents the sum over the eight Goldstone diagrams while the second one gives the contribution of each diagram after integrating over the internal variables k and i, j, l, m which are the electronic momentum and the electronic bands index of the intermediate states respectively. N k is the number of electronic states over which the sum over the momentum k is performed. The explicit expressions of M α,ijlm µ,ν (k, q) for a given diagram α is given in App. B. The generic from of M α,ijlm µ,ν (k, q) reads :

M α,ijlm µ,ν (k, q) = (-1) S(P α ) ×P α      f m f l f j f i × γ λ in k i,m g µ k±q m,l g ν k±q l,j γ λout k j,i ( m -l -iδ)( l -j -iδ)( j -i -iδ)     
where P α is a permutation of the electron-phonon vertexes g and the electron-light vertexes γ that gives rise to the Goldstone diagram α and (-1) S(P α ) its signature. In the denominators, refers to the total (electron+phonon+photon) energies of the states involved during the scattering and δ is the electronic broadening (electron+hole life times). Finally, the prefactor f i,j,l,m corresponds to the statistical occupation number of the electronic states i, j, l, m.

If one performs the summations in Eqs. 2.47 and 2.48, the next step to obtain the Raman cross section is to perform the sum over all final states that have an outgoing photon in a given direction with a given Raman shift ω = Ω in -Ω out . This sum is performed over all the possible two phonon states (weighted by the Bose-Einstein distribution) with momentum conservation. If we denote this sum by I(ω) then Eq. 2.46 reads :

I(ω) ∝ 1 N q q,µ,ν |M µ,ν (q)| 2 δ(ω -ω q,µ -ω -q,ν ) × [n(ω q,µ ) + 1] × [n(ω -q,ν ) + 1] (2.49)
where M µ,ν (q) is given by Eqs. 2.47 and 2.48. In the previous equation, N q is the total number of the final states (phonon states) over which the summation over the variable q is performed. This finalizes our derivation of the Raman cross section for second order Raman scattering.

Electron correlations effects beyond LDA

Besides the effects of the electrons interaction that are already included in the exchange and correlation term at the DFT-LDA level, additional effects of this interaction may arise. This is mainly due to the long range nature of the coulomb interaction which is not fully described by LDA within DFT. One remedy to this lack of description, is the use of Many Body Perturbation Theory (MBPT) to describe the so-called quasi-particle energies. In this section, we will examine the implication of electron correlations beyond LDA on electrons, phonons and electron-phonon coupling within the framework of MBPT at the GW level. A full discussion of the GW methods is not presented here, but the reader can find well established references about these approaches; see for example Refs. [START_REF] Aryasetiawan | The gw method[END_REF][START_REF] Hybertsen | Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies[END_REF][START_REF] Hedin | New method for calculating the one-particle green's function with application to the electron-gas problem[END_REF]. Usually, DFT-LDA results are used as bare quantities so that they serve as a starting point for more sophisticated many body calculations.

Effects on electron energies

In graphene, few-layer graphene (FLG), and graphite, Angle Resolved Photoemission Spectroscopy (ARPES) shows that LDA underestimates the slope of the electronic bands and the trigonal warping 8 [START_REF] Grüneis | Tight-binding description of the quasiparticle dispersion of graphite and few-layer graphene[END_REF] as shown in Figs. 2.2 and 2.3 for the case of 8. The trigonal warping refers to the deviation from the isotropic dispersion of the electron energies as described by the low energy Hamiltonian near the Fermi level in graphene, FLG and graphite. It can already be seen in a tight-binding description of these systems if higher nearest neighboring coupling is included, see. Ref. [START_REF] Grüneis | Tight-binding description of the quasiparticle dispersion of graphite and few-layer graphene[END_REF] (a) (a) (b) (c) graphite. The experimental measured Fermi velocity is shown to be ∼ + 20% higher than the LDA value [START_REF] Grüneis | Electronelectron correlation in graphite: a combined angle-resolved photoemission and first-principles study[END_REF]. If one includes electron correlations at the GW level (G 0 W 0 approximation) [START_REF] Grüneis | Tight-binding description of the quasiparticle dispersion of graphite and few-layer graphene[END_REF], one finds that the electronic energies are in better agreement with the experimental results as shown in Fig. 2.2 and the Fermi velocity is indeed renormalized by + ∼ 18%. Namely, one finds that, within GW [START_REF] Grüneis | Electronelectron correlation in graphite: a combined angle-resolved photoemission and first-principles study[END_REF][START_REF] Grüneis | Tight-binding description of the quasiparticle dispersion of graphite and few-layer graphene[END_REF] :

GW n,k ( DF T n,k -f ) × 1.18 + f (2.50)
where GW,DF T n,k are the (GW, DFT-LDA) electronic energies of the band n at electron momentum k and f is the Fermi level. Consequently, the trigonal warping is also shown to be affected by electron correlations. This can be seen if one tries to plot constant energy contours around the K -H high symmetry line in graphite for example [START_REF] Grüneis | Tight-binding description of the quasiparticle dispersion of graphite and few-layer graphene[END_REF] as shown in Fig. 2.3.

Important remark.

For the sake of clarity, we provide in this paragraph some important details about the GW correction for electrons. In our approach no real GW calculation was performed and our calculations of the electron energies in FLG and graphite were entirely based on the results reported in the two references [START_REF] Grüneis | Electronelectron correlation in graphite: a combined angle-resolved photoemission and first-principles study[END_REF] and [START_REF] Grüneis | Tight-binding description of the quasiparticle dispersion of graphite and few-layer graphene[END_REF]. For graphite, that was treated in Ref. [START_REF] Grüneis | Electronelectron correlation in graphite: a combined angle-resolved photoemission and first-principles study[END_REF], the authors have calculated the electronic Kohn-Sham energies within DFT-LDA and corrected them within GW. Then, a tightbinding Hamiltonian for graphite was constructed upon these GW-corrected electronic bands in order to avoid the heavy computational load of a full GW approach to describe the quasi-particle energies for graphite in the entire BZ. For graphene and FLG that were treated in Ref. [START_REF] Grüneis | Tight-binding description of the quasiparticle dispersion of graphite and few-layer graphene[END_REF], the authors used the GW-tight-binding parameters that were found in their previous work on graphite (i.e. in Ref. [START_REF] Grüneis | Electronelectron correlation in graphite: a combined angle-resolved photoemission and first-principles study[END_REF]) to construct new tight-binding Hamiltonians for graphene and few layer graphene. These GW-based tight-binding Hamiltonians were then used to describe the quasi-particles energies in FLG. The important finding was the fact that the electronic bands in graphite and FLG undergo the same rescaling law which is given by Eq. 2.50 ( see Ref. [START_REF] Grüneis | Tight-binding description of the quasiparticle dispersion of graphite and few-layer graphene[END_REF] for more details).

Effects on phonon frequencies and the electron-phonon coupling

The inclusion of electron correlations beyond DFT-LDA may affect the phonon frequencies because of their dependence on the electronic band energies and the electronphonon coupling. We already saw in the previous paragraph how the electronic energies are renormalized (with respect to DFT-LDA ones). Now we will examine how the electron-phonon coupling changes within GW.

The electron-phonon coupling can be computed with linear response as explained in Sec. 2.1.2. However, nowadays, the use of this technique within GW is not feasible. Alternatively, the EPC associated with a phonon mode can be determined by the variation of the electronic band energies by displacing the atoms according to the considered mode (the so-called frozen phonon approach). Indeed, this method was applied in Ref. [START_REF] Lazzeri | Impact of the electronelectron correlation on phonon dispersion: Failure of LDA and GGA DFT functionals in graphene and graphite[END_REF] for graphene and graphite (around the Γ and K points in the BZ). Following Ref. [START_REF] Lazzeri | Impact of the electronelectron correlation on phonon dispersion: Failure of LDA and GGA DFT functionals in graphene and graphite[END_REF] one can show that for the TO branch (highest optical branch) near K we have 9 :

D 2 K = lim d→0 1 8 ∆E K d (2.51)
where D 2 K is the deformation potential that is proportional the square of the electronphonon coupling and ∆E K is the π-bands splitting at the K points due to the phonon distortion calculated within DFT+GW. The parameter d measures the atomic displacement from equilibrium due to the phonon distortion. The previous equation was obtained in Ref. [START_REF] Lazzeri | Impact of the electronelectron correlation on phonon dispersion: Failure of LDA and GGA DFT functionals in graphene and graphite[END_REF] using symmetry consideration of the low energy Hamiltonian in graphene near the K-points. The authors of Ref. [START_REF] Lazzeri | Impact of the electronelectron correlation on phonon dispersion: Failure of LDA and GGA DFT functionals in graphene and graphite[END_REF] show that the GW-corrected Figure 2.4 -Dispersion of the highest optical phonon in graphite near K. Calculations are from DFT, or corrected to include GW renormalization of the electron-phonon coupling dots, triangles and squares ares experimental data. For more details about the calculation of DFT and DFT+GW phonon dispersion, see Ref. [START_REF] Lazzeri | Impact of the electronelectron correlation on phonon dispersion: Failure of LDA and GGA DFT functionals in graphene and graphite[END_REF]. Adapted from Ref. [START_REF] Lazzeri | Impact of the electronelectron correlation on phonon dispersion: Failure of LDA and GGA DFT functionals in graphene and graphite[END_REF].

electron-phonon coupling is enhanced by almost +60 % due to electron correlations. More precisely near K we have :

D 2 K GW D 2 K DF T ∼ 1.6 (2.52)
Now, since we know how the electronic band energies and the electron-phonon coupling are renormalized due to electron correlations, we turn to discuss the implications of electron correlations on phonon frequencies in graphene and graphite. Let us first remind ourselves the expression of the dynamical matrix within linear response (DFPT) ( [START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF][START_REF] Lazzeri | Impact of the electronelectron correlation on phonon dispersion: Failure of LDA and GGA DFT functionals in graphene and graphite[END_REF]) :

D q = Dq + Π q (2.53) Π q = 4 N k knm |D k+qm,kn | 2 × (f kn -f k+qm ) kπ -k+qπ * (2.54)
where f is the Fermi-Dirac distribution function and Dq and Π q are the analytic and non-analytic contributions respectively (see Ref. [START_REF] Calandra | Adiabatic and nonadiabatic phonon dispersion in a wannier function approach[END_REF] for explicit expression) 10 . k and q are the electron and phonon momentum. The sum in the equation above is performed over the whole BZ. In graphite, the experimental phonon frequencies of the highest optical mode around K in the BZ are lower than the DFT prediction by almost 100 cm -1 [START_REF] Lazzeri | Impact of the electronelectron correlation on phonon dispersion: Failure of LDA and GGA DFT functionals in graphene and graphite[END_REF] (see Fig. 2.4). This huge discrepancy can be taken into account within GW which affects the Π term in Eq. 2.54. The denominators in the expression of Π are larger in GW than in DFT-LDA (see previous section). Thus, inclusion 10. The non analytic term is defined as the contribution where singularities of the phonon frequencies are expected of the GW correction in the electronic bands alone results in a larger denominator providing a higher phonon frequency than the DFT-LDA one and a worse agreement with experiments. The overestimation of the TO phonon frequency near K (hence the underestimation of the TO phonon slope near K) in DFT is, thus, due to the electron-phonon coupling. Indeed, the calculated phonon dispersion with the GWcorrected electron-phonon coupling is in a good agreement with experiments as shown in Fig. 2.4 (dotted line). In view of the discussion outlined in the last two sections, we can conclude that the effects of electron correlations beyond DFT-LDA are of great importance in Raman simulations. Thus, the GW correction will be fully taken into account in our approach. In Chap. 4 and 5 we will see that the GW-dressed phonon frequencies induce a redshift in peak positions while the GW-renormalized electron-phonon coupling changes the dispersion of the peak positions as a function of the laser energy.

Important remark. The previous discussion concerns only graphene and graphite. FLG systems, on the other hand, were not considered. Except for the bilayer case, neither experimental measurements nor theoretical calculations show how electronelectron correlations may affect the phonon frequencies in FLG systems. In our approach, we assume that the in-plane vibrations in FLG, and more precisely, the HBO branches near the K point undergo the same correction as the one for graphene and graphite. This means that, the formula in Eq. 2.52 will be used to correct the electronphonon coupling in all the systems that are treated in our work. The validity of this assumption in mainly due to the in-plane character of the corrected phonon branches (TO branches).

A real calculation of the phonon frequencies within GW would require a prohibitively expensive GW calculation of the total energy and its derivatives with respect to atomic displacements. This is unfortunately not affordable by actual computational facilities and justify the frozen-phonon approach described previously. In the rest of this thesis, we will refer to the phonon frequencies that are calculated within this GW-frozen-phonon approach as the GW-corrected phonon frequencies.

Difficulties of ab-initio resonant Raman calculation

One can distinguish two main categories of difficulties in resonant Raman cross section calculations. The first one concerns the problem of excited states. Though this has already been discussed (in the overview of Sec. 2.1), here, we give a reminder of its origin. In resonant Raman scattering, the excited states visited by the electron have to be evaluated, in principle, within a full many body approach. Unfortunately, this type of calculations is, nowadays, limited to coarse grids due to its heavy computational load. The phonon states should also be computed on the excited Born-Oppenheimer energy surfaces, because the corresponding electronic states (in the intermediate transitions) are not the electronic ground state. These two observations represent fundamental difficulties if one tempts to evaluate resonant Raman intensities (and optical responses in general) in an ab-initio (parameter-free) approach [START_REF] Merlin | Overview of phonon raman scattering in solids[END_REF][START_REF] Profeta | Theory of resonant raman scattering of tetrahedral amorphous carbon[END_REF]. In the off-resonance regime (external light energy less than the material band gap), however, calculations are not concerned by this difficulty and the dielectric approach to Raman intensities give remarkable results 11 (e.g. the WURM project [START_REF] Bobocioiu | The wurm project -a freely available web-based repository of computed physical data for minerals[END_REF][START_REF] Knoll | Raman scattering of atomic vibrations in anharmonic potentials[END_REF]). Nowadays, one can say that Raman intensities calculation in the resonant regime still a challenging problem. In our approach, we believe that the inclusion of many-body effects at the GW level as presented in the last two sections provide us with an accurate description of electron and phonon energies. This description is reliable since in the electrons case, the energy windows spanned by the laser energy usually used in Raman experiments fall within an energy widows for which one has a simple formulas to correct the electronic energies when many-body effects are taking into account (the +18 % factor). For phonons, the situation is the same. Only the TO phonon branch near the K point in the BZ is concerned by the GW correction. For this phonon branch a simple GW-frozen phonon approach (that was used in the work of Lazzeri et al. [START_REF] Lazzeri | Impact of the electronelectron correlation on phonon dispersion: Failure of LDA and GGA DFT functionals in graphene and graphite[END_REF]) has been shown to give a good description the observed experimental TO phonon dispersion.

The second category of difficulties concerns the high computational time needed to evaluate the matrix elements required in Raman cross section calculations. This difficulty is independent from the resonant character of Raman scattering and applies to resonant and non-resonant calculations of Raman intensities. Contrary to first order Raman for which the phonon wave vector is restricted to zone center phonons only (fundamental selection rule), for second order Raman scattering, the phonon wave vector q in Eq. 2.49 can be any vector in the BZ. This implies that one needs to calculate the scattering by phonons (i.e. electron-phonon coupling) from the whole BZ. Furthermore, it is very well known that the spectral resolution of Raman experiments is very high and can reach the standard DFT-LDA precision on the phonon frequencies (few cm -1 ). For example, a standard Raman experiment in a FLG sample can detect peak shifts of the order of ∼ 2 cm -1 [START_REF] Lui | Imaging stacking order in few-layer graphene[END_REF]. This implies that the electron and phonon properties needed for the Raman cross section calculation have to be evaluated in ultradense sampling of the BZ which requires a high computational cost. In this thesis, it is this second category of difficulties which is treated. For an examination of the the previous difficulty in a greater detail, let us first establish the scaling law of the computational time (t) in a given Raman simulation within our approach. We already demonstrated that the Raman cross section is given by :

I(ω) ∝ 1 N q q,µ,ν |M µ,ν (q)| 2 δ(ω -ω q,µ -ω -q,ν ) × [n(ω q,µ + 1)] × [n(ω -q,ν + 1)] (2.55)
where :

M µ,ν (q) = 1 N k 8 α=1 ijlm k M α,ijlm µ,ν (k, q) (2.56)
By looking at the loops entering the previous formula, one can easily establish the following scaling law :

t ∝ 8 × N q × N k × N 4 bands × N 2 modes (2.57)
11. It is worthy to stress, however, that if the laser energy starts to approach the material band gap without matching any electronic transition, many body effects can be seen. This situation was investigated in the case of silicon recently in Ref. [START_REF] Gillet | First-principles study of excitonic effects in raman intensities[END_REF] where the exciton effects have been examined.

where N bands is the number of bands included in the Raman simulation, usually taken to be the the number of states near the Fermi level. In the case of graphene, FLG and graphite N bands is simply the number of atoms per unit cell N atoms . N mode is the number of phonon modes : N modes = 3 × N atoms . This results in N 6 atoms in total which is very high. Concerning N k and N q , typical ultra-dense phonons and electrons momentum grids of 480 × 480 were used in Refs. [START_REF] Herziger | Two-dimensional analysis of the double-resonant 2D raman mode in bilayer graphene[END_REF][START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF] to converge the Raman cross section in single and double layer graphene. More recently, grids up to 920 × 920 were used to converge the Raman cross section of the G peak in graphene [START_REF] Reichardt | Ab initio[END_REF]. Recall that in each point of these grids one must know the following matrix elements :

-Electron bands energy -Phonon frequencies.

-Electron-Phonon coupling -Electron-photon coupling It is clear that the available DFT implementations are not suitable for such kinds of heavy calculations in these ultra-dense grids. The most time-consuming steps are clearly the calculations of phonon frequencies and the electron-phonon coupling on ultra-dense grids. The reason for that is the fact that these matrix elements have to be evaluated self-consistently (e.g. taking into account the screening effects as discussed in Sec. 2.1.2 ). Moreover, one usually uses grids such as 12×12 to calculate phonon frequencies and electron-phonon coupling using standard DFT implementations (PWscf, abinit..). Previous attempts to overcome the electron-phonon calculation difficulty is the use of tight-binding approaches [START_REF] Popov | Two-phonon raman bands of bilayer graphene: revisited[END_REF][START_REF] Popov | Electron-phonon and electron-photon interactions and resonant raman scattering from the radial-breathing mode of single-walled carbon nanotubes[END_REF][START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF] which gave a remarkable agreement with experiments in the case of graphene, bilayer an carbon nanotubes. However, in these approaches, the electron-phonon is not screened which may affect the obtained spectra. For that reason, we preferred not to follow a tight-binding approach. The calculations of electron energies and electron-photon coupling are less time-consuming steps since they can be calculated non-self-consistently. Note that the time scaling low has N 6 atoms dependence. This imposes a limit on the maximum number of layers for which one can calculate the Raman cross section since, in FLG, we have N atoms = 2 × N layer .

In our approach, we have adopted two strategies in order to overcome the previous computational problems. The first strategy is to use the Wannier interpolation scheme to calculate the previous matrix elements in ultra-dense grids with an ab-initio precision. This interpolation technique reduces the computational load required by Raman simulations. The second strategy consists of an automatic reduction technique which allows us to reduce the number of points in the BZ for which one has to calculate the Raman cross section. This latter strategy is mainly based on the resonant character of 2D Raman peaks in graphene, FLG and graphite. In the next chapter we will give a detailed presentation of Wannier interpolation schemes and our reduction method and test their validity in the case of single and double layer graphene.

Maximally localized Wannier functions

In a crystal, a first-principles calculation yields a set of Bloch eigenfunctions labeled by a band index n and a reciprocal wave vector k in the BZ. Wannier functions are a set of functions w nR (r) labeled by a band like index n and a real Bravais lattice vector R such that [START_REF] Marzari | Maximally localized generalized wannier functions for composite energy bands[END_REF][START_REF] Calandra | Adiabatic and nonadiabatic phonon dispersion in a wannier function approach[END_REF][START_REF] Souza | Maximally localized wannier functions for entangled energy bands[END_REF]:

w nR (r) = 1 √ N k k e -ik.R ψ nk (r) (3.1)
where N k is the number of electron momenta for which the Bloch functions have been calculated. It is always possible to choose a normalization for the Bloch states such that the generated Wannier functions satisfy the orthogonality condition

w nR | w n R = δ n,n δ R,R
. Unlike Bloch functions, that are delocalized over the whole crystal, Wannier functions are localized functions by construction 1 [START_REF] Marzari | Maximally localized generalized wannier functions for composite energy bands[END_REF]. Hence, Wannier functions constitute a natural orthonormal and localized basis set. Historically, this is the main reason for which Wannier functions were used [START_REF] Wannier | The structure of electronic excitation levels in insulating crystals[END_REF]. As the phase of Bloch states can be chosen arbitrary, Wannier functions are not unique. In fact, the following definition :

w nR (r) = 1 √ N k k e iφ(k) e -ik.R ψ nk (r) (3.2)
where φ(k) is any real function that is periodic in reciprocal space, provides also a good definition of Wannier functions. This definition would not change the physical description of the system, however, it does change the localization properties of the generated Wannier functions as a function of r. More generally, for a composite set of N bands 2 that are isolated from the rest of the spectrum (composite bands), any definition of of the from :

w nR (r) = 1 √ N k km e -ik.R U nm (k) ψ mk (r) (3.3)
where U (k) is unitary matrix is a good definition for Wannier functions. This corresponds to a mixing (rotating) of Bloch states ψ nk before the sum over k is performed. Different choices of the matrix U (k) leads to different localization properties of the generated Wannier functions. Marzari and Vanderbilt [START_REF] Marzari | Maximally localized generalized wannier functions for composite energy bands[END_REF], have exploited this degree of freedom in choosing the U (k) matrix (called gauge freedom) to construct Wannier functions that are maximally localized. The strategy of Ref. [START_REF] Marzari | Maximally localized generalized wannier functions for composite energy bands[END_REF] consists in minimizing the sum of the quadratic spreads of the Wannier probability distributions |w nR (r)| 2 defined by :

Ω = N n r 2 n -r 2 n (3.4)
where the sum is over the chosen group of bands and r α n = ´rα |w nR (r)| 2 dr. The resulting Wannier functions upon the minimization of the spread Ω over all the possible choices of the matrix U are called Maximally Localized Wannier Functions (MLWFs). Nowadays it is possible to determine MLWFs using the Wannier90 package (Mostofi et al. [START_REF] Mostofi | wannier90: A tool for obtaining maximally-localised wannier functions[END_REF], 2008). Moreover, the original algorithm of Marzari and Vanderbilt was also extended to deal with situations where a composite set of bands cannot be identified (entangled bands). This happens notably when the bands of interest can not be isolated (e.g. metals near the Fermi level). For these cases, Souza et al. [START_REF] Souza | Maximally localized wannier functions for entangled energy bands[END_REF] in 2001 developed a band disentanglement procedure, which extracts a subset of composite bands out of a larger set of states.

It turns out that the spread minimization criterion for the determination of MLWFs is of great fundamental and practical application. MLWFs are usually comparable in size to atomic orbitals. This makes them ideally suited for interpolation since their localization enforce the localization of operators needed to be interpolated (for a review see [START_REF] Marzari | Maximally localized wannier functions: Theory and applications[END_REF]). Other localized basis set (e.g. atomic orbitals) can also be used to interpolate ab-initio quantities. The main advantage of MLWFs over other localized basis set is that no fit procedure is needed to reproduce the ab-initio electronic bands which are by construction, as we will see in the next paragraph, equal to the calculated ab-initio bands. Furthermore the number of Wannier functions needed to obtain a good description of the electronic structure is smaller.

Exponential localization of MLWFs was demonstrated in one spatial dimension by Kohn [START_REF] Kohn | Analytic properties of bloch waves and wannier functions[END_REF] and by Brouder et al. [START_REF] Brouder | Exponential localization of wannier functions in insulators[END_REF] for system with time-reversal symmetry in two and three dimensions. The more complex problem of exponential localization of ML-WFs for composite bands in 2D and 3D was finally proven by Panati and Pisante [START_REF] Panati | Bloch bundles, marzari-vanderbilt functional and maximally localized wannier functions[END_REF]. In the case of entangled bands, no exponential localization is guaranteed. However, the disentanglement procedure of Souza et al. [START_REF] Souza | Maximally localized wannier functions for entangled energy bands[END_REF] always yields highly localized Wannier functions.

Wannier Interpolation schemes

Suppose that the minimization of the spread Ω in Eq. 3.4 was performed (wannierization step). Our task will be to interpolate the matrix elements of a given operator which is known on a coarse momentum grid N k 1 ×N k 2 ×N k 3 (called the wannierization grid) into a denser grid N k1 × N k2 × N k3 (called the interpolation grid). We start with the simple case of operators of the form O(k) where k is the crystal momentum of the electron. Examples of such operators are the Hamiltonian and the dipole. Operators of the form O(k, q) which depend on two points in the BZ need some special treatment and will be studied later in the case of electron-phonon vertex g(k, q).

Interpolation of the Hamiltonian and the dipole matrix elements

The minimization of Ω yields the unitary matrix (denoted by U) which rotates the initial Bloch states in order to maximize the localization of the associated Wannier functions. The rotated Bloch states leading the MLWFs reads :

ψ (W ) nk = N m=1 U nm (k) |ψ nk (3.5)
where |ψ nk are the original Bloch states (from ab-initio calculation) and N is the number of bands of interest 3 . Note that, if the bands of interest can not be isolated from the rest of the spectrum (entangled bands), one needs a disentanglement step first to extract the relevant bands from the rest of the spectrum and obtains the matrix U (see Ref. [START_REF] Yates | Spectral and fermi surface properties from wannier interpolation[END_REF]). The type of operators O to be interpolated here are Hamiltonian and dipole like operators which are diagonal in k-space i.e. they do not couple Bloch functions with different crystal momentum. This means :

ψ mq | O |ψ nk = δ q,k ψ mk | O |ψ nk ≡ δ q,k ψ mk | O(k) |ψ nk (3.6)
The matrix elements of O in the rotated Bloch functions basis reads :

O (W ) nm (k) = ψ (W ) mk O(k) ψ (W ) nk = ψ mk | U † (k)O(k)U(k) |ψ nk (3.7)
and its Fourier transform reads :

O (W ) nm (R) = 1 N k k e -ik.R O (W ) nm (k) (3.8)
This operation is done once and for each of the N k lattice vectors R lying in a supercell conjugate to the wannierization grid. One can show that the previous matrix elements

O (W )
nm (R) are the matrix elements of the operator O between Wannier functions [START_REF] Yates | Spectral and fermi surface properties from wannier interpolation[END_REF]:

O (W ) nm (R) = w n0 | O |w mR (3.9)
Once the matrix elements O

O (W ) nm ( k) = R e i k.R O (W ) nm (R) (3.10)
where k lies on the interpolation grid. If the interpolation grid is equivalent or included in the wannierization grid, then :

O (W ) nm ( k) = O (W )
nm (k). This means that the interpolated quantities are, by construction, equal to the ab-initio ones on the wannierization grid. For a general interpolation grid, the inverse Fourier transform in Eq. 3.10 leads to accurate interpolated quantities only if the overlap between Wannier functions at different supercell is negligible as shown in Fig. 3 nm ( k) in Eq. 3.10 are expressed on the rotated basis. To find the matrix elements between the original Block functions we need the rotation matrix U( k) at the specific momentum k. Recall that the U matrices are known on the wannierization grid only. In order to find these matrices on the interpolation grid we proceed as follows. We start by noting that, the matrices U are completely independent of the type of operator O. If one repeats the interpolation algorithm described previously not for any operator but for the Hamiltonian H then, arriving to Eq. 3.10, we use Eq. 3.7 to derive the following relation :

H (W ) nm ( k) = ψ (W ) mk H( k) ψ (W ) nk = U † ( k)H( k)U( k) nm (3.11)
where H( k) is the Hamiltonian at wave vector k. Inverting the previous relation leads to :

H nm ( k) ≡ ψ mk | H( k) |ψ nk = U( k)H (W ) ( k)U † ( k) nm = kn δ n,m (3.12)
The last equality is due to the fact that H( k) is diagonal in the original Bloch basis. Thus, the matrix U † ( k) is nothing but the matrix that diagonalizes H (W ) ( k). The computational time to find U from the previous equation (i.e. by diagonalize H (W ) ), is negligible since the matrices H (W ) are small matrices in the Wannier representation (H (W ) ( k) is of dimension N × N where N is the number of bands of interest). Once, the matrices U( k) have been calculated, one can use them to express the operator O (W ) nm ( k) (Eq. 3.10) in the original Bloch states basis. This finalizes the Wannier interpolation algorithm. In summary, the important thing in the Wannier interpolation scheme is the localization of the Wannier functions which justifies the use of denser reciprocal grids in Eq.3.10. As discussed above, The need of the inverse rotation matrix at the interpolation momentum k, implies that, before interpolating any operator, we need to perform the interpolation of the Hamiltonian and diagonalize it in order to extract the inversion matrices U( k).

Interpolation of electron phonon coupling and phonon frequencies

Let us first rewrite the expression of the electron-phonon coupling (EPC), introduced in the second chapter, as :

g ν nm (k + q, k) = s e s qν • d s mn (k + q, k)/ 2M s ω qν (3.13)
where

d s mn (k + q, k) = ψ k+qm δV SCF δu qs |ψ nk (3.14)
is called the deformation potential [START_REF] Calandra | Adiabatic and nonadiabatic phonon dispersion in a wannier function approach[END_REF][START_REF] Giustino | Electron-phonon interaction using wannier functions[END_REF]. In these expressions e qν and ω 2 qν are the eigenvectors and eigenvalues of the dynamical matrix at phonon momentum q. The phonon branch index is denoted by ν. The quantity u qs is the Fourier transform of the atomic displacement u ls where the index s runs over the atoms within the l th unit cell. The integration the previous equation is understood to be on the unit cell. The quantity V SCF is the static self-consistent potential from ab-initio calculations. Before presenting the EPC and phonon frequency interpolation, we first present the Wannier interpolation of the deformation potential.

Interpolation of deformation potential. Let us define the quantity :

d s mn (R e , R p ) = w m0 | ∂V SCF ∂u ps |w nRe (3.15)
which is the representation in the Wannier functions basis of the deformation potential The quantities R e and R p are real space lattice vectors. The term between bracket is the variation of the self-consistent potential with respect to a displacement of the atom at position R p + τ s . The bra and ket are the Wannier functions already introduced. By inserting the definition of Wannier functions in terms of Bloch states (Eq. 3.3) in the previous equation, one can show that the deformation potential defined in Eq. 3.14 can be rewritten as [START_REF] Giustino | Electron-phonon interaction using wannier functions[END_REF][START_REF] Calandra | Adiabatic and nonadiabatic phonon dispersion in a wannier function approach[END_REF]

d s mn (k + q, k) = 1 N k ReRp e i(k.Re+q.Rp) × m n U mm (k + q) d s m n (R e , R p ) U † nn (k) (3.16)
which define the Wannier-to-Bloch transform of the deformation potential. In the previous expression, N k is the number of k-points in the Wannierization grid. One can also define the inverse Fourier transform (Bloch-to-Wannier transform) to express d s mn (R e , R p ) as a function of d s mn (k + q, k). Namely, we have [START_REF] Giustino | Electron-phonon interaction using wannier functions[END_REF][START_REF] Calandra | Adiabatic and nonadiabatic phonon dispersion in a wannier function approach[END_REF]:

d s mn (R e , R p ) = 1 N q ReRp e -i(k.Re+q.Rp) × m n U † mm (k + q) d s n m (k + q, k) U nn (k) (3.17)
where N q is the number of phonon momenta used in the calculation of the deformation potential (i.e. the number of phonon momenta where the dynamical matrices have been calculated). The last two equations define the Wannier interpolation of the deformation potentials. Namely, in the Wannier-to-Bloch step in Eq. 3.16 , one can use denser interpolation grids if the following two requirements are satisfied [START_REF] Giustino | Electron-phonon interaction using wannier functions[END_REF] :

Wannier interpolation

55

-The wannierization grid must be commensurate with the phonon grid over which the ab-initio calculation of g ν mn (k, q) has been carried out. In other words, one can always find a reciprocal vector G such that : k + q = k + G, where k belongs to the wannierization grid.

-The quantity d s mn (R e , R p ) must decay rapidly as a function of |R e | and |R p |. The first requirement is due to the fact that the Wannier interpolation requires only one single minimization of the spread functional to determine the MLWFs. In the Bloch-to-Wannier step (i.e. Eq. 3.17) we need both the matrices U(k) and U(k + q). Thus, if k + q does not belong to the wannierization grid, in order to obtain the matrix U(k + q), one has to perform a new wannierization starting from a new set of Bloch functions where the state at k + q is included. This yields to a set of matrices U(k) completely different from the ones obtained during the first wannierization. If one tries to use the matrices U(k) from the two wannierization procedures, the localization properties of the matrix elements to be interpolated are completely lost. In the purpose of this thesis, the wannierization grid has been always taken to be the same as the grid over which we have calculated the dynamical matrices.

The second requirement is essential to have well interpolated matrix elements. Generally speaking, the decay of d s mn (R e , R p ) as a function of |R e | is driven by the localization of the Wannier functions, while the decay of d s mn (R e , R p ) as a function of |R p | is controlled by the decay rate of the deformation potential ∂V SCF /∂u ps in real space [START_REF] Giustino | Electron-phonon interaction using wannier functions[END_REF]. We mainly distinguish two cases. In non polar semiconductors this requirement is proven to be satisfied (see Ref. [START_REF] Giustino | Electron-phonon interaction using wannier functions[END_REF] and references therein). In metals the situation is more complicated due to the presence of Kohn anomalies. A full discussion of the localization of the deformation potential is beyond the scope of this thesis and will not be presented. The reader is referred to Refs. [START_REF] Calandra | Adiabatic and nonadiabatic phonon dispersion in a wannier function approach[END_REF][START_REF] Giustino | Electron-phonon interaction using wannier functions[END_REF] and references therein for more elaborated studies. In our work, the localization of d s mn (R e , R p ) will be verified numerically at each step for graphene, FLG, and graphite.

Interpolation of phonon frequencies.

It is worthy to stress that the Wannier interpolation of electron-phonon coupling requires the knowledge of phonon frequencies and eigenmodes on the interpolation grid. Actually, in Eq 3.13 one needs ω qν and e qν in order to obtain g ν nm (k, k + q) from the interpolated deformation potential. Thus, we must first interpolate the dynamical matrices. To do that, two cases can be distinguished :

-The phonon dispersion is smooth over the whole BZ.

-The phonon dispersion suffers from the presence of Kohn anomalies and the dispersion is not smooth.

In the first case, a Fourier interpolation of the dynamical matrix is sufficiently safe generate e qν and ω qν in the interpolation grid. This is allowed since, in this case, the interatomic force constants matrix (IFC) is short ranged. Conversely, in the second scenario, the presence of Kohn anomalies indicates a long rang IFC and a direct Fourier interpolation is not safe. Unfortunately, this is the case for graphene, FLG, and graphite as discussed in the first chapter. This second case need a special treatments.

Following Calandra et al. [START_REF] Calandra | Adiabatic and nonadiabatic phonon dispersion in a wannier function approach[END_REF], if Kohn anomalies are present, we first separate the dynamical matrix D(q) into two terms. Namely we define [START_REF] Calandra | Adiabatic and nonadiabatic phonon dispersion in a wannier function approach[END_REF] :

D(q) = D(q) + Π(q) (3.18)
where D(q) and Π(q) are the analytic and non-analytic terms arising from the short and long range parts of the IFC respectively. The non-analytic terms 4 (phonon selfenergy) is given by 5 [START_REF] Calandra | Adiabatic and nonadiabatic phonon dispersion in a wannier function approach[END_REF]:

Π sr (q) = 2 N s kij f ki -f k+qj ki -k+qj × d s ij (k + q, k) × d r ij (k, k + q) (3.19)
where f is the Fermi distribution function and ki are the Kohn-Sham energies. N s is the number of electron momenta used to converge the non-analytic term at the considered smearing. After this decomposition, since the analytic term is smooth (short range IFC), it can be Fourier interpolated safely. On the other hand, the interpolation of the non-analytic term must be carried out using Wannier functions [START_REF] Calandra | Adiabatic and nonadiabatic phonon dispersion in a wannier function approach[END_REF].

As one can see from Eq. 3.19, in order to interpolate the non-analytic term, one needs the interpolated deformation potential and the interpolated Kohn-Sham energies first.

After obtaining these matrix elements, one can obtain the non-analytic term and sum it together with the Fourier interpolated analytic term and construct the interpolated dynamical matrix. The last step will be to diagonalize this matrix and obtain the interpolated eigenmodes and phonon frequencies.

Interpolation of electron phonon coupling

Only after obtaining the interpolated deformation potential and the interpolated phonon frequencies and eigenmodes, one can use the relation in Eq. 3.13 to obtain the interpolated EPC.

Summary of Wannier interpolation

In summary, the practical implementation of the above theoretical formulation for Wannier interpolation schemes consists of :

1. Perform a first-principles self-consistent calculation for the ground state, then a non-self-consistent calculation to generate the Kohn-Sham states |ψ nk and energies kn on the wannierization grid.

2. Generate the Maximally localized Wannier functions from the Kohn Sham basis obtain in step 1. The only information needed from this step is the rotation matrices U(k). We do not need the explicit expressions of the generated MLWFs.

3. Perform a standard linear response calculation to obtain the dynamical matrices D(q) and the deformation potential d(k + q, k) on a grid which is taken to be the same as the wannierization grid used in step 2.

4. In practical implementation of phonon frequency interpolation, the expression of the analytic term is not relevant since it can be obtained simply by D(q) = D(q) -Π(q) where D(q) are the standard output of DFPT calculations and Π(q) is given by Eq. 3.19.

5. Contrary to Ref. [START_REF] Calandra | Adiabatic and nonadiabatic phonon dispersion in a wannier function approach[END_REF], here we work in the adiabatic regime only.

4. Calculate the dipole matrix elements γ λ nm (k) on the wannierization grid for all polarizations λ. As discussed in the previous chapter, the dipole matrix elements (electron-light coupling in the dipole approximation) are proportional to momentum operator matrix elements.

Perform a Wannier interpolation of the Hamiltonian to obtain the Kohn-Sham

energies on the interpolation grid kn . We also need the unitary matrix that diagonalize the Hamiltonian at each interpolated electron momenta k. This is done following the procedure in Sec. 1.

6. Perform a Wannier interpolation of the dipole matrix elements to obtain γ λ nm ( k) in the interpolation grid. This is done using the matrix elements γ λ nm (k) from step 4 and the unitary matrices U from steps 2 and 5 following the procedure in Sec. 1.

7. Wannier Interpolate the deformation potential following the procedure in paragraph 2.

8. Perform a Wannier interpolation of the Dynamical matrices as explained in paragraph 2.

9. Obtain the Wannier interpolation EPC from the interpolated quantities obtained in the last two steps.

Case study : single and double layer graphene

Application of the presented interpolation schemes to the case of single and double layer graphene is presented in this section.

Localization

Fig. 3.2 shows the localization of the Hamiltonian, the dipole and the deformation potential matrix elements in the Wannier basis for both compounds for a wannierization grid of 6 × 6. In practical calculations, one has to test many wannierization grids and choose the smallest one which reproduces most accurately the ab-initio electronic bands. In the case of graphene an bilayer graphene grids from 3 × 3 to 9 × 9 were tested. It appears that, the wannierization grid 6 × 6 is the smallest grid to reproduce accurately the ab-initio results.

Interpolation of electron energies

Interpolation of phonon frequencies

Recall that the 2D peak is the overtone of D peak, meaning that, it is assisted by the TO phonons near K [START_REF] Thomsen | Double resonant raman scattering in graphite[END_REF]. Thus, an accurate description of the TO branch dispersion near K is crucial in order to correctly simulate the 2D Raman peak. The interpolated phonon frequencies with the inclusion electrons correlation at the LDA+GW level is presented in Fig. 3.4 for graphene and bilayer. The comparison with LDA phonon frequencies (Black curves) shows the importance of electron correlations.

0 2 4 6 8 | R | (Å) 1 Arb. unite | H 11 (R) | 0 2 4 6 8 | R | (Å) 1 | γ 12 (R) | 0 2 4 6 8 | R | (Å) 1 | d 1 11 (R,0) | | d 1 11 (0,R) | 0 2 4 6 8 | R | (Å) 1 Arb. unite | H 11 (R) | 0 2 4 6 8 | R | (Å) 1 | γ 12 (R) | 0 2 4 6 8 | R | (Å) 1 | d 1 11 (R,0) | | d 1 11 (0,R) |

Reduction method

After presenting the Wannier interpolation, now we are able to calculate the Raman cross section for moderate electron and phonon grid that are denser than ab-initio grids but not dense enough to converge Raman cross section calculations. In order to achieve this convergence on ultra-dense grids, one needs an other computational time reduction step. In this section we present the BZ reduction method.

General formulation

The double resonant Raman intensity for two phonons scattering is given by (see previous chapter Eq. 2.49) :

I(ω) ∝ 1 N q qµν |M µν (q)| 2 δ(ω L -ω qµ -ω -qν ) × [n(ω qµ + 1)] × [n(ω -qν + 1)] (3.20)
where ω qµ and n(ω qµ ) are the phonon frequencies and the Bose distribution for mode µ, respectively. The laser energy is denoted by ω L . The sum is performed over a uniform phonon-momentum grid in the Brillouin zone (BZ) where N q is the total number of points. The scattering probability for a specific couple of phonon modes µ, ν is :

M µ,ν (q) = 1 N k 8 α=1 ijlm k M α,ijlm µ,ν (k, q) (3.21)
where k here stands for the electron wave vector and M α,ijlm µ,ν (k + q) with α = 1, .., 8 are the scattering probabilities (α denotes a specific Goldstone diagram). Explicit expressions for these probabilities are given in App. B. The sum in Eq. 3.21 is performed over a uniform electron-momentum grid in the BZ (called electron grid) where N k is the total number of points. Ultra-dense phonons and electrons momentum grids are needed to compute the Raman cross section (summations in Eqs. 3.20 and 3.21). Typically, momentum grids of 480 × 480 × 1 were used in Refs. [START_REF] Herziger | Two-dimensional analysis of the double-resonant 2D raman mode in bilayer graphene[END_REF][START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF] for monolayer and bilayer graphene. However, only a small portion of the grids actually satisfies resonant conditions and contributes to I(ω). Here below we explain how we manage to condense points in the relevant part of the grids, without calculating the cross section for the portion of the BZ not satisfying the resonant condition.

Phonon Grid

For the phonon grid (sum over q in Eq. 3.20), we proceed as follows: 1. We first calculate I(ω) in Eq. 3.20 for uniform and coarse electron and phonon grids. 2. We calculate the intensity I(ω) on the coarse phonon grid and select the set of points that gives a substantial contribution to the cross section, namely those satisfying the condition:

I pp qµν > I max S ph (3.22)
Where S ph is a parameter (typical values ranges from 10 to 200). This subset of q-points constitute the so-called resonant region for phonon momenta. We also verify that Eq. 3.20 restricted to the resonant region leads to the same I(ω) for the coarse grid.

Electron grid

A similar procedure is adopted for identifying the resonant region in for electronmomenta. To do that, we have to examine the exact form of the kernel M α,ijlm µ,ν (k + q) (see App. B) The denominators in the expression giving M α,ijlm µ,ν (k + q) have the following form (for the sake of simplicity we neglect the phonon energy):

(ω L -π * k + π k -iδ) × (ω L -π * k±q + π k -iδ) × (ω L -π * k±q + π k±q -iδ) (3.23)
Where π * k±q ( π k±q ) are the energies of the electronic π * and π states, δ is the electronic+hole lifetimes and ω L is the laser energy. According to the double resonant model [START_REF] Thomsen | Double resonant raman scattering in graphite[END_REF][START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF] the first term (creation of the electron-hole pair : ω L -π * k + π k -iδ) and the third term (annihilation of the electron-hole pair : ω L -π * k±q + π k±q -iδ) in Eq. 3.23 are responsible for the double resonance during the Raman scattering [START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF]. The non-resonant regions for electron-momenta are the set of k-points of the grid where, at least, one of these two denominators is off-resonance. In order to avoid these regions we adopted the following procedure :

1. We first calculate I(ω) in Eq. 3.20 for a coarse electron grid. 2. We determine the resonant region for electrons from the following conditions that have to be simultaneously verified for each phonon momentum q determined in the previous subsection:

| π * k -π k | < ω L + S el × δ (3.24) | π * k±q -π k±q | < ω L + S el × δ (3.25)

General procedure for Raman simulations 61

Where S el is a parameter ranging typically from 5 to 15.

3. We recalculate the I(ω) and verify that it does not differ significantly from the one calculated with the full coarse grid.

Once the resonant regions of the electron and phonon grid are obtained on coarse grids, we proceed as follows :

-Fist, we converge the sum over k in Eq. 3.21 keeping the coarse phonon grid, namely we consider the electron-momenta belonging to the resonant region and we shift them randomly of an amount given by the distance between the k-points times a random number between 0 and 1 6 . This generates many randomly displaced resonant regions. We calculate the cross section summing over all electron grids and calculate the spectrum until we are at convergence.

-Second, we generate a new randomly displaced resonant region for phonon momenta as we did for the electrons. For this new grid we converge again the sum over electron-momenta as done previously. We then calculate the cross section for each randomly displaced resonant region for phonon momenta. If n-randomly displaced resonant regions are generated, we calculate the cross section and obtain the n-th approximant I n (ω). We then calculate

I n (ω) = 1 n n 1 I i (ω) (3.26)
When I n (ω) is converged, we stop.

Case study : single and double layer graphene

The application of the reduction procedure described above for phonon grids in the case of single and double layer graphene gives the results depicted in Fig. 3.5. The results concern the 2D peak only. This means that I max introduced in Eq. 3.22 is the maximum intensity of the 2D peak. The reduced grids shown in Fig. 3.5 were obtained with S ph = 50 for both compounds. As one can see in Fig. 3.5, the set of points that satisfy the condition in Eq. 3.22 falls within a corona around the K and K points in the BZ. This result is in a complete agreement with the double resonant model where the important phonon momenta for the 2D Raman peak are predicted to be located around the Dirac points [START_REF] Thomsen | Double resonant raman scattering in graphite[END_REF].

General procedure for Raman simulations

In this section, we summarize the whole methodology of calculating the Raman spectrum using Wannier functions interpolation schemes and the BZ reduction method. First, we perform a first-principles calculation of the Kohn-Sham states and generate the corresponding MLWFs. This is to be followed by a linear response calculation 6. Generate a grid randomly means that each time a grid is generated, the position of the points inside this grid is shifted by a small random vector. This is just a computational procedure to generate dense grids for the calculation of Raman spectra avoiding symmetry equivalent points.

Graphene Bilayer

SCF + n-SCF Linear response Wannierization

Storage Storage Storage

Figure 3.6 -Schematic that shows : the first-principles calculation (scf+n-scf step), the wannierization, and the linear response steps. The obtained quantities have to be stored for later use. The notations are as follows : kn are the Kohn-Sham energies, D(q) are the dynamical matrices, g ν nm (k + q, k) are the electron-phonon coupling matrix elements, and γ nm,λ k are the dipole matrix elements. The matrices U(k) are the unitary matrices that lead to MLWFs. of the matrix elements of : H, γ and g on the wannierization grid where we have calculated the Kohn-Sham states (see Fig. 3.6). Second, we perform the BZ reduction method on coarse electron and phonon grids. At the end of this step, we will have the resonant regions for electrons and phonons. Finally, we use these resonant regions to generate dense electron and phonon grids (as dense as required by the Raman cross section to converge) and calculate the converged Raman spectrum (see Fig. 3.7).

Wannier interpolation

Raman cross section calculation R = ---------------------- The label "Theory" refers to our results. The experimental data are from Ref. [START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF] for graphene and from [START_REF] Herziger | Two-dimensional analysis of the double-resonant 2D raman mode in bilayer graphene[END_REF] for bilayer.

Tests of Raman simulations in single and double layer graphene.

In this section, we present the results of tests of our methodology in the case of single and double layer graphene. For the sake of comparison with previous theoretical works reported on bilayer graphene [START_REF] Herziger | Two-dimensional analysis of the double-resonant 2D raman mode in bilayer graphene[END_REF], the spectra that are presented here were obtained with an electronic broadening of 84 meV and a phonon broadening of 8 cm -1 .

For graphene, Fig. 3.8 (a) shows the excellent agreement between the experimental and the calculated 2D Raman peak. The peaks position, width, and shape are well described, in particular, the symmetric Lorentzian shape of the 2D peak is well reproduced. The same good agreement is observed for the D+D peak at around 2450 cm -1 and the asymmetric Lorentzian shape of this peak is also well reproduced. Furthermore the intensity ratio between the 2D and D+D peaks is in a good agreement with experiments.

For bilayer graphene the agreement is good with previous theoretical and experimental works [START_REF] Herziger | Two-dimensional analysis of the double-resonant 2D raman mode in bilayer graphene[END_REF]. The shape of the 2D peak and all its sub-peaks are well reproduced. However, we observe a blueshift of the theoretical spectra with respect to the experimental ones. Also, the theoretical spectra are larger than the experimental ones. This discrepancy was already encountered in the work of F. Herziger et al. [START_REF] Herziger | Two-dimensional analysis of the double-resonant 2D raman mode in bilayer graphene[END_REF] in which the calculation of resonant Raman intensities was carried out using an approach which equivalent to ours but do not include the BZ reduction algorithm. Thus, this discrepancy can not be attributed to the reduction of the BZ.In our work, we confirm that the more broadened theoretical 2D peak is not due to the phonon lifetime. Calculations at lower phonon broadening yield a 2D peak which is always larger than the experimental one.

The experimental Raman spectra of bilayer graphene presented in Fig. 3.8 (b) were Raman shift (cm -1 ) Raman intensity (a.u.)

0 V 0.8 V 1.6 V 2.4 V 3.2 V

eV

Figure 3.9 -The Raman 2D peak in bilayer graphene as a function of doping (adapted from Ref. [START_REF] Das | Phonon renormalization in doped bilayer graphene[END_REF] where doping was achieved by the application of an external electric field). measured in freestanding samples where doping due to the presence of the substrate can occur. It is well known that, in bilayer graphene the 2D Raman peak sharpens with increasing doping as we can see from Fig. 3.9. Thus, the discrepancy between our theoretical calculations and the experimental measurements, in the case of bilayer graphene, can be attributed to the fact that the Raman measurements were performed on doped samples.

For the case of graphene, our BZ reduction procedure allows to gain a factor of 300 on the computational effort with the respect to the brute force calculation (no use of the resonant regions). This factor is less for bilayer. This is expected since the resonance in the case of graphene is very sharp (a sharp 2D peak) so that the resonant regions (for electrons and phonons) in the BZ are sharper than for graphene for bilayer graphene (see Fig. 3.5).

Chapter 4 Application for 2D systems : trilayer and tetralayer graphene.

This chapter is organized as follows. First, a state-of-the-art about Raman simulations in few-layer graphene (FLG) is given. Second, Wannier interpolated electron and phonon properties in trilayer and tetralayer graphene are presented. Third, Raman spectra that were calculated using the methodology developed in the last two chapters are studied and the Raman signatures of different stackings in FLG is established. Fourth, we present a detailed analysis of the obtained Raman spectra in ABA and ABC trilayer graphene. We also show how Raman spectroscopy can be used to extract the dispersion of the TO phonon modes in ABA-trilayer graphene. Finally, we explain the difference between the 2D Raman modes in ABA and ABC trilayer graphene and reveal its origin.

State-of-the-art of Raman simulation in few-layer graphene

For single layer graphene, previous work done by P. Venezuela et al. [START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF] shows that a tight-binding approach gives a good agreement with experiments. The bilayer case was also investigated more recently by V. Popov [START_REF] Popov | Two-phonon raman bands of bilayer graphene: revisited[END_REF] within a tight-binding approach and a good agreement between theory and experiments has been obtained. Bilayer graphene was also investigated by F. Herziger et al. [START_REF] Herziger | Two-dimensional analysis of the double-resonant 2D raman mode in bilayer graphene[END_REF] using the same methodology as the one developed in this work except for the reduction method that was not used. Their approach was also shown to give a good agreement with experiments.

Up to now, Raman simulations for trilayer and tetralayer graphene were discussed in term of symmetry-allowed processes [START_REF] Malard | Group-theory analysis of electrons and phonons in n-layer graphene systems[END_REF][START_REF] Malard | Raman spectroscopy in graphene[END_REF][START_REF] Dresselhaus | Characterizing graphene, graphite, and carbon nanotubes by raman spectroscopy[END_REF] (15 in the case of trilayer graphene) and, to our knowledge, no explicit calculation of the corresponding Raman spectra has been reported. Contrary to the bilayer case, for trilayer and tetralayer graphene, Raman simulations are extremely expensive without the reduction method.

In the remaining of this chapter we will present the application of the methodology developed in the previous two chapters for the calculation of Raman spectra in trilayer and tetralayer graphene for all possible stacking.

Electron bands dispersion

Different stackings can occur during the formation of trilayer graphene (see Fig. 4.1). The most stable and common ones are : the Bernal (ABA) and the rhombohedral (ABC) stackings. Other possibilities which involve a sequence of type "AA" (two equivalent graphene layers) are predicted to be less stable [START_REF] Charlier | First-principles study of the stacking effect on the electronic properties of graphite(s)[END_REF] 1 . Bernal and rhombohedral trilayer graphene have been already observed and differentiated experimentally using optical absorption [START_REF] Lui | Imaging stacking order in few-layer graphene[END_REF][START_REF] Cong | Raman characterization of aba-and abc-stacked trilayer graphene[END_REF]. Although the rhombohedral stacking is less stable, 15 % of the experimental samples were identified to be in this stacking mode [START_REF] Lui | Imaging stacking order in few-layer graphene[END_REF].

In both stackings, the π orbitals overlap gives rise to 3 valence bands, and 3 conduction bands labeled π 1..6 from low to high energy near the Fermi level at the K points. Figs 4.2a and 4.2b show the Wannier interpolated electronic bands compared to the DFT-LDA ones for ABA and ABC stackings. The results depicted in these 1. Although unstable, stackings of type AA can be found in defected few-layer graphene samples with grain boundaries. For example, in bilayer graphene (AB stacking) if the grain boundary is well adapted, one can have a large area of AA stacking due to the reorientation of the second layer [START_REF] Brown | Twinning and twisting of tri-and bilayer graphene[END_REF]. In tetralayer graphene the number of non-equivalent stacking modes (excluding stackings that contain "AA" sequences) are the three following possibilities : ABAB (Bernal stacking), ABCA (rhombohedral stacking) and ABAC (mixed stacking which is equivalent to ABCB). In all the these stackings, the π orbitals overlap gives rise to 4 valence bands and 4 conduction bands (4 atoms/unit cell). These three polytypes have close electron and phonon bands structure. Using the same wannierization grid (6 × 6) as for the trilayer case, and 4 Wannier functions, we obtain an excellent agreement between the Wannier interpolated and ab-initio bands for tetralayer graphene. The inclusion of electron correlations effects on the electron energies is shown in Figs. 4.4 for the three stacking modes in tetralayer graphene.

Phonon frequencies

Remark. In all the figures that are presented in the this chapter (and the next one as well) we have adopted the following convention for the legends that concerns both electrons and phonon bands:

-The legend LDA refers to the ab-initio LDA bands.

-The legend Wan. LDA and Wan. LDA+GW refer to the Wannier interpolated LDA and LDA+GW results respectively.

-The legend LDA+GW refers to the GW corrected ab-initio LDA bands. ABC trilayer graphene suggests that the difference between their Raman 2D mode can not be attributed to the difference of their phonons bands structure. This difference at the half the frequency of the 2D mode (around 1350 cm -1 ) is of few cm -1 and can not be the main responsible for the observed difference between the Raman 2D modes in ABA and ABC trilayer graphene. In Sec. 4.6 we investigate more rigorously this observation and we show that the difference between the 2D Raman peaks in ABA and ABC stacked trilayer graphene is, indeed, not due to the phonon dispersions variation, but has an electronic origin. Recall that, in both trilayer and tetralayer graphene, the 2D Raman peak is assisted by the TO phonon branch near the K points [START_REF] Thomsen | Double resonant raman scattering in graphite[END_REF]. Hence, an accurate description of this branch is crucial. The TO splitting, which occurs between the TO sub-branches (3 in trilayer graphene and 4 in tetralayer graphene) near the K point, is larger in LDA+GW compared to LDA calculations. For trilayer graphene (see Fig. 4.5), the maximum splitting of the TO phonons near the K point is 6 cm -1 in LDA, while in GW it is 10 cm -1 . The same enhancement of the TO splitting is observed for tetralayer graphene for all possible stackings (see Fig. 4.5). The TO splitting is of great importance in the calculation of Raman cross sections since it determines the position of the sub-contributions to the 2D Raman mode of different symmetry-allowed processes, as it will be discussed in the analysis section of this chapter.

Phonon frequencies

Technical details. first-principles calculations were performed using the Quantum ESPRESSO package [START_REF] Giannozzi | Quantum espresso: a modular and open-source software project for quantum simulations of materials[END_REF]. The ground state calculation was performed using Perdew-Zunger (LDA) norm-conserving pseudopotential with a frozen 2s 2 electrons. and an energy cutoff of 70 Ry. Electronic ground state integration was performed with a smearing of 0.001 Ry and a 64 × 64 sampling of the BZ. The hexagonal unit cell parameter was fixed to a 0 = 2.46 Å and the interlayer spacing was fixed to 3.35 Å. A vacuum of 10 Å along the z axis was used to avoid periodic image effects. The phonon frequencies calculation was performed on a 6 × 6 × 1 grid and the electronphonon coupling was calculated using the Wannier information encoded in the ".chk" and ".eig" files from the wannierization step. This calculation was performed using the ph package of Quantum ESPRESSO . The same package (ph) was also used to carry out the calculation of the electron-photon coupling (we have actually calculated the momentum matrix elements which are proportional to the dipole matrix elements; see Chap. 2).

Calculated spectra

In this section, we report a complete comparative study (in trilayer and tetralayer graphene) between the Raman 2D modes of all possible stackings (ABA and ABC for trilayer and ABAB, ABCA, and ABAC for tetralayer graphene). The calculated spectra for trilayer are depicted in Fig. 4.7 for several laser energies compared to experimental data. The calculated spectra for tetralayer graphene are presented in Fig. 4.8 for the three possible stackings at laser energy 2.41 eV and 2.54 eV. In the case of ABCB tetralayer graphene, no experimental spectra have been reported for this stacking to our knowledge, that is why we have plotted the theoretical results only.

In both trilayer and tetralayer graphene, the overall agreement with experiments is good. Especially, all the sub-peaks of the 2D Raman modes and their relative intensities are well reproduced for all stackings. Furthermore, our calculations reproduce correctly the experimental blueshift of the 2D Raman modes position with an increasing laser energy. This fact is due to the correct description of the linear behavior of the TO branch near the K as explained in the next paragraph.

Since the 2D mode is double resonant, if we change the laser energy, the phonon wave vector which satisfies the resonance condition changes accordingly. More specifically, the resonant wave vector as calculated from the K point is larger with higher laser energy. This can be seen by a simple geometrical construction as shown in Fig. 4.11 or, more rigorously, from Fig. 4.14 in the next sections. As a consequence, the 2D peak position also changes, giving rise to the observed blue shift with increasing laser energy. The inclusion of the GW-correction for the TO phonon branch renormalizes its phonon slope near the K point as explained in Chap. 2. The renormalized phonon slope (Figs. 4.5, red curves) controls the change of the resonant phonon wave vector and consequently, the dispersion of the 2D Raman peak position as a function of the laser energy. Thus, a good description of the phonon slope, as provided by the inclusion of the GW correction, leads to the correct blueshift of the 2D peak.

The most important consequence of the above good agreement between theory and experiment is that : Raman simulations give for each stacking (in trilayer and tetralayer graphene) its own Raman signature that is the shape of its 2D Raman peak. This signature, which we have calculated for each stacking, can be used to identify this stacking directly. To be more quantitative, in Fig. 4.9 we show for each laser energy the comparison between the theoretical and the experimental difference of the ABA and the ABC trilayer spectra (after area normalization). As one can notice, the theoretical curves follow the experimental ones reproducing their maxima and minima. This results allow us to say that the main objective of the present thesis has been achieved in the case of trilayer and tetralayer graphene.

Technical details. Double resonant cross sections were calculated completely from first-principles. The electron-phonon and dipole matrix elements were first calculated within DFT-LDA and then interpolated using Wannier functions. The reduction method was used with the converged parameters S el = 9 and S ph = 70 (see the previous chapter, Sec. 3.2 for the definition of S el and S ph ). The electronic broadening δ was calculated as in Ref. [START_REF] Herziger | Two-dimensional analysis of the double-resonant 2D raman mode in bilayer graphene[END_REF] for double layer graphene, namely we have defined :

δ = 0.081832 × ω L 2 -0.1645 (meV) (4.1)
where ω L (in eV) is the laser excitation energy. The phonon broadening (FWHM) was fixed to 5.3 cm -1 following L. Paulatto et al. [START_REF] Paulatto | Anharmonic properties from a generalized third-order ab initio approach: Theory and applications to graphite and graphene[END_REF] where phonon-phonon interaction was studied to determine the phonon lifetime from first-principles.

The double resonant Raman cross section was calculated on ultra-dense phonon and electron grids for reciprocal space integration, namely, grids of 300 × 300 and 256 × 256 were used for electrons and phonons respectively.

Intensity

Analysis of trilayer graphene Raman spectra

In this section we present several types of analysis that were performed on the obtained spectra for trilayer graphene in order to understand the structure of the 2D peak. The performed analysis is quite similar to the one reported in Refs. [START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF] and [START_REF] Herziger | Two-dimensional analysis of the double-resonant 2D raman mode in bilayer graphene[END_REF] for single and double layer graphene. The results presented in next paragraphs were obtained at laser energies of 1.96 eV and 2.71 eV. An analysis for other laser energies can be found in App. D.

Let us first remind the mechanism responsible for the Raman scattering that gives rise to the 2D peak (see Fig. 4.10b). First, an electron-hole pair is created due to annihilation of the incoming photon. Second, the electron (or the hole) is scattered by a phonon with momentum q. Third, the electron (or the hole) is scattered by a phonon with momentum -q (-q from momentum conservation). Finally, the electron-hole pair recombines and the outgoing photon is emitted. According this description, many processes can take place depending on the electronic states between which the creation (recombination) or the scattering by the phonon occurs. However, dipole selection rules and phonon symmetries reduce the total number of possible combinations between initial and intermediate electronic states, to a set of symmetry allowed processes. As an example, graphene has only one process and bilayer graphene has four. Here, we examine the case of trilayer graphene.

In ABA trilayer graphene the symmetry of electronic bands is as follows [START_REF] Malard | Group-theory analysis of electrons and phonons in n-layer graphene systems[END_REF] (see Fig. 4.10a)

π 1 → T -π 2 → T + π 3 → T - and π 4 → T -π 4 → T + π 6 → T -
where T ± are appropriate representations of the symmetry group of electronic bands near the K point. The incoming light belongs to the T + representation which connects bands with the same symmetry only. This leaves us with 5 allowed optical transitions that are presented in Fig. 4.10a.

TO phonon symmetries have also their implications on Raman scattering. Recall that trilayer graphene has 3 TO branches near the K point which can be either symmetric or antisymmetric. In the second step of the Raman process described previously, a symmetric phonon will scatter the electron (or the hole) between electronic states of the same symmetry (T -→ T -or T + → T + ). This gives rise to the so-called symmetric processes as shown in Fig. 4.10b. Conversely, an antisymmetric phonon will scatter the electron (or the hole) between electronic states of different symmetry (T -→ T + or T + → T -). This gives rise to the so-called antisymmetric processes as shown in Fig. 4.10b. By counting all the possible combinations between the allowed optical transitions and the allowed phonon scattering, we find that trilayer graphene has 15 = 5 × 3 allowed Raman processes that contribute to the 2D Raman peak. We define the symmetric (anti-symmetric) contribution to the Raman cross section as the sum of all the symmetric (anti-symmetric) processes.

In the case of ABC trilayer graphene, due to the lack of group theory analysis, we have explored the allowed Raman process numerically. By calculating the dipole matrix elements between electronic bands near the Fermi level, we found that ABC trilayer graphene has the same allowed optical transitions as ABA trilayer (5 transitions). This allowed us to assign the symmetry properties of the electronic bands in ABA trilayer to ABC ones. In order to classify the allowed Raman processes in ABC trilayer graphene into symmetric and anti-symmetric processes, we proceed as follows 2 . After an electron is excited into an empty state from the conduction band near the K point by absorbing the incoming photon, it will scatter to a state near the K point. This state may have the same or the opposite symmetry as the state near K (see Fig. 4.10). Now, we define a symmetric Raman process in ABC trilayer graphene as a process where the electron scatters between electronic states with the same symmetry. Conversely, an anti-symmetric process is defined as the process where the electron scatters between electronic bands with opposite symmetry. In this way, we will have the same Raman allowed processes for ABA and ABC trilayer graphene which is very useful for the sake of comparison between the Raman 2D peaks of the two systems.

In the case of tetralayer graphene, the number of allowed Raman processes is very large. These processes give many contributions which interfere to form the 2D structure and it is a hard task to isolate and study every single one of them. For this reason, we did not perform any analysis for tetralayer graphene.

Polarization study

Fig. 4.12 (a) shows, for ABA and ABC trilayer graphene, light polarization analysis of the Raman cross section at 1.96 eV (the same analysis at 2.71 eV is presented in Fig. 4.13 (a)). We define the xx (resp. xy) contribution as the contribution where the incident light is polarized along the x direction and the scattered light is polarized along the x (resp. y) direction. The contributions xy and yx are defined in a similar way. As one can see in Fig. 4.12 (a) the symmetry x ≡ y is preserved in our calculation3 . This is not surprising since FLG have a hexagonal symmetry which implies that the Raman tensor should be symmetric. If we define I to be the parallel contribution to the Raman cross section where the scattered light is parallel to the incident one and, conversely, the perpendicular contribution I ⊥ , to be the contribution where the scattered light is perpendicular to the incident one, then trilayer graphene have the same ratio I I ⊥ ≈ 2.7 as for graphene [START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF]. This ratio is found to be constant and independent of the laser excitation energy.

Inner-outer and e-h contributions

Fig. 4.12 (b) show, for ABA and ABC trilayer graphene, the electron-hole (e-h) contribution to the Raman cross section for both stackings as well as the inner-outer decomposition at 1.96 eV (the same analysis at 2.71 eV is presented in Fig. 4.13 (b)). Let us first give some definitions. The e-h contribution is defined by the sum of all Raman processes where both the hole and the electron are scattered. Conversely, electron-electron : e-e (hole-hole : h-h) contributions are defined by the sum of Raman processes were only the electron (hole) is scattered twice. A process is called inner if the resonant phonon wave vector stems from a sector of ±30 • next to the K -Γ direction with respect to the K point in the BZ. These processes involve electrons that are near the K -M line as shown in Fig. 4.11 (a) and (c). Conversely, outer processes have phonon wave vector that stems from ±30 • next to the K -M direction and involve electrons that are near the K -Γ line as shown in Fig. 4.11 (a) and (b) [START_REF] Mafra | A study of inner process double-resonance raman scattering in bilayer graphene[END_REF]. As in single layer [START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF] and bilayer graphene [START_REF] Herziger | Two-dimensional analysis of the double-resonant 2D raman mode in bilayer graphene[END_REF], e-h contribution is dominant compared to the e-e and h-h ones. Furthermore, we notice that the inner contribution is dominant over the outer one for ω L = 1.96 eV. By increasing the laser energy, the outer contribution increases and becomes almost comparable to the inner one for laser energy ω L = 2.71 eV. This behavior is different from what happens in single layer [START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF] and double layer graphene [START_REF] Herziger | Two-dimensional analysis of the double-resonant 2D raman mode in bilayer graphene[END_REF] where the dominant processes are the inner ones independently from the laser excitation energy. It is well known in the literature that the matrix elements describing the coupling of the electron to photons and phonons is the main responsible for the dominance of the inner processes over the outer ones (see Ref. [START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF] for example). An outer contribution which is comparable to the inner one could be explained by the homogeneity between the coupling of the electron to phonon and photon in both K -M (outer) and K -Γ (inner) directions.

Processes decomposition

In Fig. 4.12 (c) we show the decomposition over the 15 allowed processes of the 2D mode in trilayer graphene for the ABA and ABC stacking. These processes are labeled following the electron transitions between filled and empty states during the Raman scattering. A Process P jl i refers to a scattering where the electron is initially in the filled state π i . It will be excited to the empty state π j and then scattered by a phonon to the empty state π l . We notice that the interference between these symmetry-allowed processes shown in in Fig. 4.12 (c) is totally additive and that their sum gives the e-h scattering spectrum. This additive interference plays an important role in the 2D structure as it enhances the sum of the sub-contributions and leads to the observed shape of the 2D peak. Although there are theoretically 15 allowed processes for trilayer graphene, we observe in Fig. 4.12 (c), for ABA and ABC trilayer, that 6 of them have a negligible contribution (the processes that are not present in these figures have negligible contribution to the 2D peak). Some of the remaining 9 processes give almost the same contribution; they are called degenerate processes (represented by the same color with solid and dashed lines). These processes are of the from: P jl i and P lj k , namely we have the 3 following pair of degenerate processes : (P 56 2 , P 65 1 ), (P 46 3 , P 64 1 ) and (P 45 3 , P 54 2 ). The remaining 3 processes, namely, the processes: P 55 1 , P 44 3 and the P 66 1 give independent contributions. In view of the previous observations, one can say that theoretically, the 2D peak in ABA and ABC trilayer graphene is composed by 6 different contributions or sub-peaks. This was already observed experimentally in Refs. [START_REF] Cong | Raman characterization of aba-and abc-stacked trilayer graphene[END_REF][START_REF] Lui | Imaging stacking order in few-layer graphene[END_REF] where the obtained 2D peaks are shown to be reproducible using at least a 6 Gaussian functions fitting. However, the authors of Refs. [START_REF] Cong | Raman characterization of aba-and abc-stacked trilayer graphene[END_REF][START_REF] Lui | Imaging stacking order in few-layer graphene[END_REF] were not able to explain why only 6 contributions are observed even if symmetry-allowed processes are 15 in total. In this section, we have given the theoretical explanation of this fact.

Symmetric-antisymmetric decomposition

Fig. 4.12 (d) shows, for ABA and ABC trilayer graphene, the decomposition of the 2D peak upon symmetric and antisymmetric processes at 1.96 eV. For this laser energy the symmetric contribution is more important then the antisymmetric one. Furthermore, the main sub-peak (the one that marks the maximum intensity at around 2650 cm -1 ) arises from symmetric processes. Our analysis at increasing laser energy shows that, the symmetric contribution decreases with increasing laser energy and it becomes comparable to the antisymmetric one at 2.71 eV as shown in Fig. 4.13 (d).

q-resolved Raman cross section (contour plot)

In Figs. 4.14 we present the contour plot of the Raman cross section for ABA and ABC trilayer graphene at different excitation energies. Contour plots show the dependence of the Raman cross section on the phonon wave vector q(q x , q y ) in the BZ, Chapter 4. Application for 2D systems more precisely, in Fig. 4.14 we plot the function I(q x , q y ) defined by :

I(q x , q y ) = Nm µ,ν M µν (q x , q y ) (4.2)
where M µν (q x , q y ) is the Raman cross section for modes µ and ν (See Chap. 2 Eq. 2.47) and N m is the number of modes (N m = 3×N atoms where N atoms is the number of atoms in the unit cell). This kind of analysis allows us to visualize the resonant regions in the BZ. As one can see in Fig. 4.12 (d), there are 6 resonant regions, for both stackings, which manifest as a set of 6 rings around the K point. By calculating the contour plot of every allowed Raman process independently, and compare it to total contour plot in Figs. [START_REF] Veithen | Nonlinear optical susceptibilities, raman efficiencies, and electro-optic tensors from first-principles density functional perturbation theory[END_REF] . This confirms previous experimental finding that the 2D mode in trilayer graphene can be decomposed in term of 6 different contributions [START_REF] Lui | Imaging stacking order in few-layer graphene[END_REF][START_REF] Cong | Raman characterization of aba-and abc-stacked trilayer graphene[END_REF]. Furthermore, this kind of analysis allows us to extract accurately the resonant wave vector for each process. The possibility to extract the exact resonant wave vector is of great usefulness as it allows for the determination of the TO phonon dispersion and the estimation of the TO phonons splitting in trilayer graphene as explained in the next section.

ABA-trilayer

TO branch dispersions from Raman spectra in trilayer graphene

In ABA trilayer graphene, one can calculate easily the contributions to the subpeaks P 55 2 which arise from electron and hole scattering by the same phonon. It appears that, the main contribution to the P 55 2 process comes from the 1 st TO mode 4 . As one can see in Fig. 4.15a, the position of this peak is the same as the main subpeak of the 2D Raman mode. Thus, half of the frequency of the main peak in the 2D Raman mode 5 gives the frequency of the 1 st TO phonon in ABA trilayer graphene. Furthermore, from the P 55 2 process contour plot, one can extract the resonant wave vector of the involved phonon. By repeating these two steps for all excitation energies used in this work, we were able to extract the dispersion of the 1 st TO phonon mode in ABA trilayer graphene from Raman spectra.

A similar analysis was carried out for the P 44 3 process. This latter has a double peak structure, where the right sub-peak is mainly due to the 3 rd TO phonon near K (see Fig. 4.15a). Using the same procedure as the one employed for the P 55 2 process, half of the experimental frequency of the last sub-peak in the 2D Raman mode gives the frequency of the 3 rd TO phonon mode. From the contour plot of the P 44 3 process, we can extract the exact resonant phonon wave vector involved in the P 44 3 Raman scattering. By repeating the previous procedure for all laser excitations, we were able to extract the experimental dispersion of the 3 rd TO mode in ABA trilayer graphene from Raman spectra. Fig. 4.15b shows our extracted phonon dispersion compared to LDA and LDA+GW calculations. As expected, the TO phonon frequencies are overestimated in LDA, be- cause of the wrong description of Kohn anomaly near K. The TO splitting is also underestimated, namely, the maximum experimental TO splitting (10 cm -1 ) is almost double of the one calculated within LDA (6 cm -1 ). The importance of this splitting in trilayer graphene arises from the fact that the inner contribution, which is the dominant contribution to the Raman scattering, is (by definition) assisted by phonons near to Γ -K high symmetry line where the TO splitting is not negligible (contrary to the K-M direction). This TO splitting controls the relative positions of the sub-peaks (or Raman processes) and, thus, determines the shape of the 2D Raman peak in trilayer graphene as shown in Fig. 4.12 (c).

The same kind of analysis was not carried out for ABC trilayer graphene. This is because the processes P 55 2 and P 44 3 do not have a main contribution that arises from just one phonon branch. Both of them have contributions from the 3 TO modes. Consequently, the frequency of the observed experimental sub-peaks of the 2D Raman mode could not be assigned to any single phonon branch.

Origin of the difference between ABA and ABC stackings in trilayer graphene

We will now turn to the discussion of the origin of the difference between the 2D Raman peaks arising from ABA and ABC trilayer graphene. By analyzing the graphs in Fig. 4.12 (b) for ABA and ABC trilayer graphene, it appears that the difference between the two stacking modes comes essentially from the inner contribution (see App. D to examine the results at higher laser energies). The outer contribution is always almost the same (between the two stackings at a given laser energy). Recall that inner processes involve electrons near the K -M line and phonons near the K -Γ line. By analyzing the electronic band dispersions of ABA and ABC trilayer graphene shown in Fig. 4.16 (a), we observe that the two stackings have the same dispersions along the Γ -K direction, while the dispersions along the K -M direction is different. During inner Raman scattering, this difference leads to an emitted phonon with momentum q which is different in the ABC stacking than in the ABA one : |q ABC inner | = |q ABA inner |. By comparing the contour plots of the Raman intensity for the two stackings, already presented in Fig. 4.14, we observe that the rings that show the resonant wave vectors in ABA and ABC trilayer graphene, are more deformed in the inner directions than in the outer ones. In Fig. 4.16 (b) we took the example of the P 55 2 processes (which is responsible of the main peak of the 2D mode) for ABA and ABC trilayer graphene at 1.96 eV. As one can see, while along the outer direction the rings superpose, for the inner direction, the ring of the ABA stacking is situated inside the ABC one indicating a resonant phonon wave vector with a module that is smaller for the ABC case than the ABA one (the module is calculated from the Γ point).

This deformation as small as it is, gives rise to a phonon frequency shift of the order of 5 cm -1 at least. This corresponds to 10 cm -1 shift of the Raman two phonon peaks. Due to this shift, the relative position of the ABA and the ABC inner sub-peaks changes and, consequently, their 2D peak structure. This idea was first developed in Ref. [START_REF] Lui | Imaging stacking order in few-layer graphene[END_REF], explored more and verified in this work.

Conclusion

In summary, we have shown that, for trilayer and tetralayer graphene, the 2D Raman peak can be used to distinguish between different stacking orders. In trilayer graphene, we have explained why the experimentally observed 2D Raman mode can be decomposed in 6 different contributions. For each single contribution, we have assigned a particular Raman process. This allowed us to fully understand the structure of the 2D Raman mode in ABA and ABC trilayer graphene. In the case of ABA-stacked trilayer graphene, a full analysis of the P 55 2 and P 44 3 processes, that were assigned to the main and last sub-peaks of the 2D Raman mode, allowed us to extract the TO phonon dispersion from experimental Raman spectra. The extracted TO-dispersion shows that the TO branches splitting is underestimated in LDA by almost a factor 2. Finally, we have shown that, the origin of the difference between the 2D Raman mode in ABA and ABC trilayer graphene has mainly an electronic origin. More precisely, while the phonon bands structure is almost the same between the two stackings, the electronic bands structure along the K -M direction is different. Raman scattering processes along this direction (inner processes) are responsible for the difference between the 2D Raman modes of ABA and ABC trilayer graphene.

Chapter 5 Application for 3D systems : Bernal and rhombohedral graphites

In this chapter, Raman simulations for the three systems : Bernal, rhombohedral, and ABCB graphites are presented. It is organized as follows. First we present the Wannier interpolated electron band energies and phonon frequencies. Second, a comparative study between the Raman response of the three types of graphite is provided. The results of this comparison yields the Raman fingerprint of what we call short and long range rhombohedral stacking orders. Third, the shape of the 2D Raman peak in Bernal and rhombohedral graphite is explained. Fourth, the analysis of the resulting spectra from Bernal graphite is presented and the structure of their 2D Raman peak is studied. Finally, the dispersion of the symmetric transversal optical (TO) phonon branch in Bernal graphite is calculated from Raman spectra and compared to LDA and LDA+GW calculations.

State-of-the-art of Raman simulation in graphite

Bernal graphite [START_REF] Bernal | The structure of graphite[END_REF] with AB stacked graphene layers is the most stable form of graphite. Recently, however, rhombohedral stacked multilayers graphene (RMG) attracted an increasing attention as theoretical calculations suggest the occurrence of a dispersionless electronic band (bandwidth smaller than 2 meV) at the Fermi level [START_REF] Kopnina | Surface superconductivity in rhombohedral graphite[END_REF][START_REF] Xiao | Density functional investigation of rhombohedral stacks of graphene: Topological surface states, nonlinear dielectric response, and bulk limit[END_REF]. This flat band with extremely large effective mass constitutes a very promising candidate for highly correlated states of matter such as magnetism [START_REF] Pamuk | Magnetic gap opening in rhombohedral-stacked multilayer graphene from first principles[END_REF] or room-temperature superconductivity [START_REF] Precker | Identification of a possible superconducting transition above room temperature in natural graphite crystals[END_REF].

As ABC-stacked graphite (rhombohedral graphite) is metastable [START_REF] Charlier | First-principles study of the stacking effect on the electronic properties of graphite(s)[END_REF], the synthesis of long sequences of ABC graphene layers is a real challenge. For a random sequence of N graphene layers stacked along the c axis, a purely statistical argument states that the probability to obtain N layers with ABC order is 1/2 N -1 1 . In reality the probability is even lower as all stackings are not equally probable as energetics favors the Bernal one with respect to the others. This explains why three and four layer graphene flakes with ABC-stacking are systematically found [START_REF] Mak | Electronic structure of few-layer graphene: Experimental demonstration of strong dependence on stacking sequence[END_REF][START_REF] Lui | Imaging stacking order in few-layer graphene[END_REF][START_REF] Cong | Raman characterization of aba-and abc-stacked trilayer graphene[END_REF], while it is highly im probable to obtain long range ABC-stacking order. Recently, it has been suggested that pentalayers graphene with rhombohedral stacking can be grown epitaxially on 3C-SiC [START_REF] Brown | Twinning and twisting of tri-and bilayer graphene[END_REF] [START_REF] Pierucci | Evidence for flat bands near the fermi level in epitaxial rhombohedral multilayer graphene[END_REF]. Finally, Henni et al. [START_REF] Henni | Rhombohedral multilayer graphene: A magneto-raman scattering study[END_REF] were able to isolate multilayer graphene flakes with ABC sequences exceeding 17 graphene sheets. However, while for three and four ABC stacked graphene layers an optical signature exists [START_REF] Mak | Electronic structure of few-layer graphene: Experimental demonstration of strong dependence on stacking sequence[END_REF][START_REF] Lui | Imaging stacking order in few-layer graphene[END_REF], a clear fingerprint of long-range rhombohedral order is lacking.

Raman spectroscopy, and in particular the 2D double resonant Raman peak, has been proven to be a very powerful technique to investigate structural and physical properties of graphene flakes [START_REF] Ferrari | Raman spectroscopy as a versatile tool for studying the properties of graphene[END_REF]. However, despite its crucial importance, the theoretical understanding of the 2D double resonant Raman spectrum has been obtained only for few-layer graphene [START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: Intensity and line shape of defect-induced and two-phonon bands[END_REF][START_REF] Basko | Theory of resonant multiphonon raman scattering in graphene[END_REF][START_REF] Popov | Theoretical 2D raman band of strained graphene[END_REF][START_REF] Herziger | Two-dimensional analysis of the double-resonant 2D raman mode in bilayer graphene[END_REF]. Even the basic case of bulk Bernal graphite is not completely understood.

The quest to understand the Raman spectrum of graphite started over 40 years ago [START_REF] Tuinstra | Raman spectrum of graphite[END_REF]. Previous theoretical works on the study Raman scattering in graphite used simplified approaches such as the constant matrix elements approximation [START_REF] Thomsen | Double resonant raman scattering in graphite[END_REF][START_REF] Reich | Raman spectroscopy of graphite[END_REF][START_REF] Cançado | Geometrical approach for the study of G' band in the raman spectrum of monolayer graphene, bilayer graphene, and bulk graphite[END_REF][START_REF] Narula | Double resonant raman spectra in graphene and graphite: A two-dimensional explanation of the raman amplitude[END_REF] where one simply assumes that the matrix elements entering the cross section are constants. These approaches allowed one to understand the origin of the Raman peaks and to calculate the phonon dispersion of graphite from experimental Raman spectra [START_REF] Saito | Probing phonon dispersion relations of graphite by double resonance raman scattering[END_REF]. However, the previously calculated spectra do not reproduce the exact shape and width of the Raman peaks, especially the 2D one.

To our knowledge, no previous studies were carried out to explore the presence of rhombohedral stacking in natural graphite using Raman spectroscopy. It is worthy nothing that the first identification of the rhombohedral phase of graphite was reported by H. Lipson and A. R. Stokes in 1942 using X-ray diffraction on natural graphite [START_REF] Lipson | The structure of graphite[END_REF]. In Ref. [START_REF] Lipson | The structure of graphite[END_REF], the authors show that, in addition to the expected lines from the hexagonal symmetry (P 6 3 /mmc space group), additional lines arising from a rhombohedral symmetry (R 3m space group) can be clearly seen. More Interestingly, their study, based on X-ray diffraction, gives the same amount rhombohedral stacking which present in natural graphite (15 %) than the one found in FLG (see Ref. [START_REF] Lui | Imaging stacking order in few-layer graphene[END_REF]).

In this chapter, we will show how Raman spectroscopy can be used to identify what we call long and short rang rhombohedral (or ABC) stacking sequences.

Electron bands and phonon frequencies

The hexagonal stacking of graphene layers gives rise to Bernal graphite (or ABgraphite) with a nearest neighboring distance of 1.42 Å and an interlayer spacing of 3.35 Å. In this stacking mode, half of the atoms in each layer are directly on top of their neighbors on the layer below and below the center of the hexagons from the layer on top, the inverse occurs for the other half (see Fig. 5.1). Bernal Graphite belongs to the P 6 3 /mmc space group, with four atoms per unit cell occupying the Wyckoff positions 2c and 2d. For graphite, as for the bilayer case, the π orbitals overlap gives rise to 2 valence and 2 conduction bands that we label by π 1...4 from low to high energy near the K point. Graphite is a semimetal with a small overlap between the valence and conduction bands near the K point in the BZ (∼ 30 meV).

The infinite stacking of ABC-trilayer graphene gives rise to rhombohedral graphite (or ABC-graphite). This stacking belongs to the the R 3m space group with two atoms in the unit cell occupying the Wyckoff positions ± (u, u, u) with u = 1 6 . The nearest neighboring distance is 1.42 Å and the interlayer spacing is taken to be the same as for Bernal graphite : 3.35 Å2 . In this work, the rhombohedral unit cell was not used to perform Raman spectra calculation in ABC-graphite. Alternatively, we have used a hexagonal supercell cell with 6 atoms. This choice is more adapted for Raman simulations since, in this case, we do not need a dense sampling of the BZ in the outof-plane direction (z direction). The price to pay for choosing this supercell is the use of 6 Wannier functions instead of 2 that can be used if we choose to work with the rhombohedral unit cell 3 . It turns out that the calculation are easily feasible with the hexagonal supercell, while they are more expensive using the rhombohedral unit cell.

In order to detect signatures of different kinds of stackings, we also perform calculations for the case ABCB bulk graphite. ABCB bulk graphite is interesting since it corresponds to a sequence of .. Electron bands. Figs. 5.2a, 5.2c, and 5.2e show the electron bands dispersion of AB, ABC, and ABCB graphites calculated within DFT-LDA and compared to the Wannier interpolated ones. The Wannier bands were obtained using a 6 × 6 × 3 wannierization grid and 2, 6, and 8 Wannier functions for AB, ABC, and ABCB graphites respectively. The obtained agreement between LDA and the Wannier interpolated bands is excellent over a large window of energy. We should stress the importance of the wannierization grid 6 × 6 × 3. Actually, before converging toward this grid, calculations of Raman spectra for Bernal graphite were performed with the wannierization grid 6 × 6 × 2. The obtained spectra were completely wrong and the shape of the calculated 2D peaks did not match the experimental data. After several tests, we realized that the disagreement was due to the bad description of the out-of-plan bands dispersion. The use of only 2 points in the wannierization grid 6×6×2 to sample the out-of-plane direction generates the wrong dispersion between the K and H points in the BZ and consequently yields the wrong spectra. This was the first hint about the importance of the out of plan electronic bands dispersion in calculating Raman spectra for bulk graphites. Later in this chapter, we will discuss the importance of this dispersion in details.

Figs. 5.2b, 5.2d, and 5.2f show a comparison between electronic bands calculated within LDA and LDA+GW for the three types of graphite. The same correction was maintained for the three compounds and taken to be the same as the one used in FLG (see previous chapter).

Phonon frequencies. In AB, ABC, and ABCB graphite, the number of transverse optical (TO) branches near the K points is equal to 2, 3, and 4 respectively. The same correction as the one used in FLG systems (see App. A) was applied for the three types of graphite. As in FLG, the splitting between the TO branches is underestimated in LDA. Fig. 5.3 shows the results for bulk AB-graphite (Bernal graphite). A careful analysis shows that, for AB-graphite, while the TO branches are almost degenerate within LDA, a maximum splitting of 12 cm -1 is observed near the K point after including the GW correction.

Bulk AB

M K H L ---> A Γ <--- Wave vector

Calculated spectra : the Raman fingerprint of rhombohedral graphite

We first consider bulk AB graphite (Bernal graphite) for which several experimental data are available. We calculate the spectra for different laser energies finding an excellent agreement with experimental data (see Fig. 5.4). We then consider in more details the spectra at ω L = 1.96 eV. The 2D peak is composed of a main sub-peak at ∼ 2683 cm -1 and a shoulder around 2640 cm -1 , as shown in Fig. 5.4. Both features are well described by the calculation. The shape and intensities of the D+D overtone structure at ∼ 2456 cm -1 , although at slightly lower energy in the calculation, are also very well reproduced. The blueshift of the theoretical spectra with increasing laser energy is also in excellent agreement with the experimental one.

For the sake of comparison, we performed calculations for the case of ABC and ABCB graphites. As we said before, ABCB bulk graphite is interesting as it corresponds to a sequence ...[ABC](BAB)[CBA](BCB)... that is an equal mixing of trilayer with rhombohedral and Bernal stackings. Thus, the differences between bulk AB and ABCB stackings can be seen as fingerprint of short (or local) rhombohedricity (i.e. few ABC layers) while the differences between bulk ABCB and bulk ABC are signatures of long range rhombohedral order. The results are depicted in Fig. 5.5, where they are compared with the spectrum of Ref. [START_REF] Henni | Rhombohedral multilayer graphene: A magneto-raman scattering study[END_REF] that has been tentatively attributed to 17 layers ABC-stacked graphene. Both the 2D peaks spectra of bulk ABC and ABCB graphite are substantially broader than the one of Bernal graphite. Thus, the increased width of the 2D peak at ω L = 1.96 eV is a fingerprint of short range ABC sequences.

Even if both bulk ABC and bulk ABCB spectra seem similar, they differ by the presence of a feature at ∼ 2576 cm -1 in the ABC case that is completely missing in the ABCB stacking. As the bulk ABC structure differs from the bulk ABCB one by the occurrence of long range rhombohedral order, the feature at ∼ 2576 cm -1 at ω L = 1.96 eV can be seen as a fingerprint of long range rhombohedral order. The good agreement

Figure 5.5 -Theoretical spectra for bulk AB, ABC and ABCB stacked graphite against experiments at 1.96 eV. In (a) the signature of the long rang ABC stacking at around ∼ 2576 cm -1 is shown. In (b) we demonstrate the absence of this signature for short rang ABC-stacking in ABCB graphite. The experimental data are from samples composed of 17 layers of ABC-stacked graphene [START_REF] Henni | Rhombohedral multilayer graphene: A magneto-raman scattering study[END_REF].

between our theoretical calculation and the experimental spectrum in Ref. [START_REF] Henni | Rhombohedral multilayer graphene: A magneto-raman scattering study[END_REF], both from what concerns the 2D peak width and shape as well as the presence of the feature at ∼ 2576 cm -1 suggests that the samples in Ref. [START_REF] Henni | Rhombohedral multilayer graphene: A magneto-raman scattering study[END_REF] contain long range sequences of rhombohedral stacked multilayer graphene. The color bar indicates the normalized q-resolved Raman cross section, where q is the phonon momentum.

Width of the 2D peak in AB and ABC-stacked graphite

The larger width of the 2D peak in bulk ABC graphite with respect to the bulk AB graphite seems counterintuitive. AB-stacked graphite has two couples of π, π * (4 atoms/unit cell) electronic bands, which give rise to four allowed Raman processes while ABC-stacked graphite has only one couple of π, π * (2 atoms/unit cell) and one Raman process. So one could naively think that AB-stacked graphite should have more allowed dipolar transitions and that, consequently, the 2D peak should be broader. However, this is without taking into account the electronic k z dispersion that, as shown in Fig. 5.7a, is substantially different along K -M -K in the two case. The different k z electronic band dispersion implies that different electron and phonon momenta contribute to the Raman cross section. This is clearly seen in Fig. 5.6 where the resonant phonon momenta contributing to the 2D peak cross sections are highlighted in a contour plot (see also Ref. [START_REF] Herziger | Two-dimensional analysis of the double-resonant 2D raman mode in bilayer graphene[END_REF] for more technical explanations). While the Bernal case includes very sharp resonances in phonon momenta, the resonance is much broader in the rhombohedral case due to the different band dispersion along z. This explain the larger width of the 2D peak in the rhombohedral case. In order to explicitly verify the effect of the electronic band dispersion on the 2D peak, we calculate the full Raman cross section with two approximations. First, we restrict the calculation to scattering with transverse optical phonons only (both for electrons and holes), as shown in Fig. 5.7b. The details the method used to calculate the TO contribution to the Raman cross section is presented in App. C. This approximation leads to indistinguishable spectra for the 2D peak from the full calculation, meaning that the 2D feature is only due to transverse optical (TO) modes. On the contrary, under this approximation the D+D peak disappears as it is in part due to scattering to longitudinal acoustic phonons, as it was shown in the case of graphene [START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF]. Second, we calculate the spectrum restricted to TO modes and we set the electron-phonon matrix elements of these phonons to unity, regardless of the electronic band index. The results of this approximation are depicted in Fig. 5.7c. As it can be seen, for both AB and ABC-stacked graphite, all the peaks occurring in the complete calculation are still present. However the intensities are different. Under this last approximation the intensities are simply related to the number of optically-active dipolar-processes contributing to the Raman cross section as they are not anymore weighted by the electron-phonon matrix element. As a result of its different k z dispersion, rhombohedral graphite has a substantially larger number of optically active dipolar processes and this explains the origin of its broader spectrum despite the reduced number of bands in the BZ.

Technical details. Double resonant spectra were calculated from first-principles using the method developed in Chap. 2 and 3. The electron and phonon bands where first calculated using the Quantum ESPRESSO [START_REF] Giannozzi | Quantum espresso: a modular and open-source software project for quantum simulations of materials[END_REF] code/ The ground state was calculated using Perdew-Zunger (LDA) norm-conserving pseudopotential with a frozen 2s 2 electrons and an energy cutoff of 70 Ry. Ground state electronic integration was performed on k-point grids of 64 × 64 × 4 for AB, ABC, and ABCB graphite. The dynamical matrices and the electron-phonon coupling were first calculated in linear response on sparse phonon momentum grid of 6 × 6 × 3 and then both were Wannier interpolated through the Brillouin zone. The electronic bands and phonon frequencies were corrected for the electrons correlations as explained in App. A. The Raman cross sections were calculated on ultra-dense electron and phonon grids for reciprocal space integration. Namely, phonon grids up to 300×300×16 and electron grids as large as 128 × 128 × 16 were reduced using the reduction method described in the Chap. 3, then used to evaluate the Raman cross sections. The electron and phonon lifetimes were chosen as in FLG and kept the same for all calculations.

Spectra analysis of AB-stacked graphite

In this section, we present the analysis that has been performed on the obtained spectra for AB-graphite. The analysis which will be presented, was performed at 1.96 eV excitation energy. The analysis at 2.41 eV laser energies is presented in App. D Before going into the details of the analysis, let us first give some useful definitions. A Raman scattering process P lm ij refers to a scattering where an electron, initially at the filled state π i , is excited to the empty state π j by the annihilation of the incoming photon, then scattered by a phonon to the empty state π l . The hole, on the other hand, is scattered by a phonon from the state π i to the state π j . Finally the recombination process takes place between the π l and π j states. The definitions of inner, outer, symmetric, antisymmetric, and e-h contributions to the Raman cross section are the same as the one employed for FLG. In AB-graphite, dipole selection rules do not impose any restriction on the optical transitions between valence and conduction bands i.e. all the transitions are allowed. However, due to the optical anisotropy in graphite Ref. [START_REF] Malard | Group-theory analysis of electrons and phonons in n-layer graphene systems[END_REF], the only allowed transitions are π 1 -→ π 4 and π 2 -→ π 3 . The electronic scattering by a phonon can be assisted either by a symmetric or an anti-symmetric mode (from the two TO branches). This gives rise to four possible processes which can be indexed by the labels of the involved four electronic bands without ambiguity as : P 44 11 , P 33 22 which are symmetric processes and P 43 12 , P 41 23 which are antisymmetric 4 .

4. Contrary to the case of trilayer graphene, the Raman allowed processes for Bernal graphite are indexed by four integers P lm ij . The j index here (which is absent in the indexation of the trilayer Raman processes) indicates the final electronic band occupied by the electron (or the hole) after the occurrence of the recombination process. I ⊥ is about 2.7 and the symmetry x ≡ y is always preserved. shows that the e-h scattering is the dominant contribution to the 2D-mode over the e-e and h-h processes. However, while in FLG e-e and h-h processes are negligible, in Bernal graphite they have a small but not negligible contribution. Analysis at different laser energies shows that, as for the case of trilayer and tetralayer graphene, the inner contribution is dominant at low energy. The outer contribution on the other hand, is shown to increase with an increasing laser energy.

Inner-outer and e-h contributions

Processes and Symmetric-antisymmetric decompositions

Processes analysis is very important since it allows us to understand the shape of Raman peaks. Fig. 5.8 (c) shows processes decomposition (at 1.96 eV) of the 2D Raman cross section upon the four allowed Raman processes in Bernal graphite discussed in the beginning of this section. As we can see, the main sub-peak of the 2D structure (the peak around ∼ 2700 cm -1 at 1.96 eV) is mainly due to the process P 33 22 , while the other processes, namely, P 44 11 , P 42 13 and P 31 24 form the shoulder of the 2D structure. If we combine these results with the inner outer decomposition from Fig. 5.8 (b) and the symmetric-antisymmetric decomposition from Fig. 5.8 (d), we can say, in a more detailed way, that the main sub-peak of the 2D Raman peak is mainly due to the inner contribution of the symmetric process P 33 22 . The outer-antisymmetric contribution forms the shoulder of the 2D structure.

TO phonon bands dispersion from Raman spectra

We have already shown in Sec. 5.4 that the 2D Raman peak in AB-graphite is due to the TO phonon modes. Moreover, in the previous section, we have also seen that the main sub-peak of the 2D peak in AB-graphite is due to the inner-symmetric contribution of the P 33 22 Raman process. The combination of these two results leads to the conclusion that the main sub-peak of the 2D Raman mode is assisted by the innersymmetric TO phonon branch. Recall that the "inner" characteristic essentially means that the wave vector of the involved phonon stands along the Γ -K direction and that the "symmetric" characteristic means that the involved phonon is the symmetric one. All these findings allow us, in principle, to retrieve the dispersion of the symmetric TO phonon in AB-graphite along the Γ -K Direction. This can be achieved using the same methodology as in ABA-trilayer presented in the previous chapter. Namely, the symmetric TO phonon frequency can be easily determined from the position of main sub-peak of the 2D mode. The only missing information is the phonon wave vector. We have seen in the previous chapter that this can be determined by calculating the contour plot of the considered Raman process (see Sec. 4.6 and Fig. 4.16). However, in the case of bulk graphite, the situation is slightly different because of the third dimension. In fact, the intensity of the 2D peak depends on the three Cartesian components of the phonon momentum q(q x , q y , q z ) :

I(ω) = q I(q x , q y , q z ) I(q x , q y , q z ) = µν M µ,ν (q x , q y , q z ) δ(ω -ω qµ -ω -q,ν ) (5.1)
where ω is the Raman shift and M µ,ν (q x , q y , q z ) is given by Eq. 2.47. Let us Assume that the dispersion of the phonon branches is negligible along q z direction (out-ofplane direction).This can be justified since the Van Der Waals interaction between the graphene layers is relatively weak. In this case, the projected contour plot of the process P 33 22 on the plan containing the Γ point will give us the resonant wave vector. To summarize, the practical method to be used in order to extract the symmetric TO phonon dispersion in AB-graphite from Raman spectra reads :

1. Perform Raman simulations of the 2D mode at different laser energies 2. For each laser energy, extract the frequency of the main sub-peak of the 2D mode.

Half this frequency represents the frequency of the involved phonon mode in the scattering.

3. For each laser energy, calculate the projected contour plot of the P 33 22 Raman process and extract the resonant momentum of the involved phonon.

4. Draw the phonon dispersion using the results from step 2 and 3.

The application of the above procedure yields the experimental phonon dispersion shown in Fig. 5.9 compared to LDA and LDA+GW calculations.

Conclusion

In this chapter we have reported a first-principles calculation of the two-phonon resonant 2D and D+D peaks in AB, ABC, and ABCB graphite. Our calculations, carried out for several laser energies, are in good agreement with experimental data available for AB-graphite. In the case of ABC-stacked graphite, we have compared our calculation with recently synthesized flakes from Ref. [START_REF] Henni | Rhombohedral multilayer graphene: A magneto-raman scattering study[END_REF] that were tentatively attributed to 17 layers with ABC sequences. Our theoretical ABC-graphite spectrum confirm this attribution. Furthermore, we have shown how to distinguish between short and long range rhombohedral order. Furthermore, we have demonstrated that the 2D peak in AB-graphite is mainly assisted by the inner-symmetric TO phonon near the K point in the BZ. This allowed us to extract the dispersion of the TO symmetric phonon in AB-graphite from experimental Raman spectra.

Conclusions and perspectives

The central subject of present thesis was the study of stacking orders in few layer graphene (FLG) and graphite systems by mean of ab-initio Raman simulations. Two main objectives have been achieved. The first one was the understanding of the 2D Raman peak's origin and structure in trilayer and tetralayer graphene for all possible stackings. This understanding allowed us to determine a Raman fingerprint for each stacking which is the shape of its 2D Raman peak. The second main objective was achieved for 3D bulk graphite systems (Bernal, rhombohedral, and ABCB-graphites) where we have shown that, the 2D Raman peak of bulk graphites carry important information about the occurrence of long and short range rhombohedral stacking (rhombohedricity). More specifically, we have demonstrated that, the eventual observed larger width of the 2D Raman peak in natural graphite samples with respect to the Bernal graphite case is a signature of short range rhombohedricity. This short rhombohedricity reflects the possible occurrence of some rhombohedral-stacked FLG inside bulk graphite. Long range rhombohedricity, on the other hand, has a distinctive Raman signature represented by the appearance of an extra-peak at 2576 cm -1 with laser excitation of 1.96 eV in the Raman spectrum of rhombohedral graphite. This extrapeak, which is absent in the Raman spectrum of Bernal graphite, has been observed experimentally in a recent work by Henni et al. [START_REF] Henni | Rhombohedral multilayer graphene: A magneto-raman scattering study[END_REF] on few layer graphene systems exceeding 17 layers. The studied samples in Ref. [START_REF] Henni | Rhombohedral multilayer graphene: A magneto-raman scattering study[END_REF] were tentatively attributed to a rhombohedral stacking of graphene layers using magneto-Raman spectroscopy. In this thesis we confirm this finding and we provide a prerequisite to experimental non-destructive identification and synthesis of rhombohedral multi-layer graphene and graphite.

Double resonant Raman scattering cross sections were calculated completely from first-principles and including electron-electron interaction at the GW level. Since the computational load of such an approach is relatively high, we have used the Wannier interpolation technique which has provided us with interpolated quantities with an ab-initio precision at low computational load. The need of an other time reduction method was highly required, especially, for bulk graphites systems where calculations of Raman spectra using the Wannier interpolation technique still extremely expensive. To do that, we have developed an automatic procedure, based on the double resonant character of the studied Raman spectra, in order to reduce the computational cost. The developed procedure consists of avoiding the regions from the Brillouin zone that do not contribute to the Raman cross section. This reduces the number of points needed to calculate the Raman cross sections and, consequently, lower computational load. The combination of the Wannier interpolation technique with the Brillouin zone 102 Chapter 5. Conclusions and perspectives reduction method allowed us to evaluate Raman cross sections in few layer graphene an bulk graphites with an ab-initio precision at at low computational load.

We believe that the achieved goals and the established methodology in this thesis will improve our understanding of Raman spectra in graphene layers and bulk graphites systems. We also expect this contribution to widen the use of Raman spectroscopy to include stacking order determination. Further studies are definitely required to go beyond the remaining limitations of the present work. In the next few paragraphs we will provide some perspectives and extensions of the present thesis.

The methodology developed in this work to simulate Raman spectra is not restricted to graphene layers and graphite only, on the contrary, it can be extended to study other systems. The Wannier interpolation can be used to interpolate the matrix elements needed in Raman simulations for any system provided that one can find a well localized Wannier functions basis set that reproduces accurately the ab-initio electronic bands. The reduction method, on the other hand, requires the existence of the resonant regions at first. A sufficient but not necessary condition for that will be the resonant character of the Raman peaks to be simulated. These are the two requirements for the extension of our methodology to other systems. Transition metal dichalcogenide (TMDs) are probably the most exciting route for new simulations using our approach. Raman spectroscopy in these materials has been already used and proven to be useful. In these materials also, a theoretical understanding of the resonant peaks is crucial in order to trace back the materials electronic and vibrational properties. However, if we are interested in characterizing semiconducting TMDs, which are more relevant for electronic applications, special care must be given for exciton effects that are usually found within this class of materials and which our approach do not include (since we are dealing with semimetals). This requires the incorporation of the Bethe-Salpeter Equation (BSE) in order to describe accurately the electron-hole pair. Furthermore, electron correlations effects on the bare electron an phonon properties (i.e. DFT and DFPT results with approximated exchange-correlation functionals) must also be included due to their importance in Raman scattering. This latter point may generate some computational issues. Recall that in our work, these effects are very well described by first-principle calculations at the GW level in the energy widows (for both electrons and phonons) relevant for Raman simulations. The same accurate description is required for TMDs if one tempts to simulate the Raman response of these materials.

An other exciting perspective of the present work would be the application of our approach to simulate ultraviolet-Raman spectra of different solids and molecules. The particularity of ultraviolet-Raman (UV-Raman) spectroscopy is the fact that it probes deep electronic levels (like the sigma bands of graphene) in the valence bands (i.e. bands that are responsible for the cohesion of the whole systems). This may provide us with important informations that can be used to investigate phase transitions for example. Our approach can be used, in principle, directly for UV-Raman simulations. However, one needs to take care of the inclusion of electron correlations since a large energy windows (many electronic bands) will contribute to the Raman cross section in this case.

A last but not least perspective, would be the extension of the present formalism to include more elementary excitations (plasmons, polarons..etc.). This perspective needs some additional theoretical development for light coupling to plasmons and other excitations. Once such formalisms are established, this route of simulation will not just allow the characterization of materials by mean Raman spectroscopy, but also the prediction of novel properties related to the existence of the described elementary excitations. Furthermore, it will also allow for the description of electronic Raman scattering (ERS) which still a hot field of research.

ABA-trilayer

Résumé en Français Introduction

Lorsqu'un matériau est éclairé par la lumière, de nombreux phénomènes peuvent être observés à la fois en surface et a l'intérieure du matériau (voir la Fig. R.A). À la surface, la lumière incidente peut être réfléchie ou transmise à l'intérieure du matériau. Une fois transmise, elle peut être diffusée, absorbée. A son tour, le rayonnement qui se propage dans le matériau peut être diffusé élastiquement ou inélastiquement selon que l'impulsion du photon incident est conservée ou pas. À ce stade là, la lumière diffusée inélastiquement par les ondes acoustiques est généralement appelée diffusion de Brillouin (ce phénomène a été découvert indépendamment par Mandelstam). d'autre part, La diffusion par les autres excitations élémentaires du matériau (phonons optiques, plasmons, ..etc.) est appelée la diffusion Raman (Sir C.V. Raman a reçu le Prix Nobel de Physique en 1930 pour la découverte, en 1928 à Calcutta, en Inde, de l'effet qui porte son nom aujourd'hui) [START_REF] Cardona | Light scattering in solids II: basic concepts and instrumentation[END_REF]. Dans cette thèse, seuls les phonons seront considérés. Par conséquent, un processus de diffusion Raman est dit du 1 er , 2 ème , .., n ème ordre suivant le nombre de phonons qui ont été émis (diffusion dite de Stokes) ou absorbés (diffusion dite de anti-Stokes). La spectroscopie Raman consiste à mesurer l'intensité de la lumière diffusée en fonction de la différence de fréquence entre les rayonnements entrants et sortants. Dans ce genre d'expérience, la source de la lumière incidente est un laser délivrant des excitations monochromatiques dont l'énergie est de l'ordre de 2 eV. Dans les molécules et les isolants à grande bande interdite (∆ > 2 eV), ce type d'excitation est incapable de promouvoir les électrons au-delà de la bande interdite (gap). Dans le diamant par exemple, où la valeur du gap est de 5.5 eV, un laser de 2 eV excitera toujours les électrons dans des états virtuels. Durant le processus de diffusion Raman, Ces électrons absorbent ou émettent un certain nombre de phonons puis retournent à la bande de valence en émettant le rayonnement diffusé. Si nous utilisons un laser avec une énergie différente qui restera bien inférieure au gap, les électrons seront toujours excités dans des états virtuels et le même processus de diffusion aura lieu. Puisque les états virtuels ne sont pas des états réels du matériau, la réponse Raman du matériaux est complètement indépendante de l'énergie de l'excitation laser. Dans ce genre de situations, où l'énergie de l'excitation laser est bien inférieure au gap, la diffusion Raman est dite non-résonnante. Dans le cas contraire (les isolants à faible gap (PbTe, SnTe) ou les métaux et les semimétaux (graphène, graphite, cuivre)) un laser à 2 eV peut toujours exciter les électrons de la bande de valence à la bande de conduction et la lumière est dite résonnante avec les excitations électroniques. Si l'énergie du laser est modifiée, les états excités impliqués seront différents, ce qui fait que, la réponse Raman est, maintenant, dépendante de l'énergie d'excitation par ce que les états électroniques impliqués sont des états réels du matériau. Dans ce genre de situation, où l'énergie du laser est en résonance avec une transition électronique, la diffusion Raman est dite résonnante. Expérimentalement, ce phénomène de résonance se manifeste par l'apparition de pics fins dans le spectre Raman.

Les phénomènes de résonance (qui constituent un point centrale de la présente thèse) sont des phénomènes très importants. En effet, la forme, l'élargissement et l'intensité des pics résonnants observés sont reliés a des propriétés fondamentales du matériau étudié (fréquences des phonons, le couplage électron-phonon ...etc). Cela offre un moyen efficace de les étudier en utilisant la spectroscopie Raman. Cependant, 123 cette étude (ou caractérisation) ne peut être achevée sans une compréhension théorique des pics résonnants et de leur lien avec les propriétés fondamentales du matériau étudié,c'est-à-dire, sans effectuer des simulations théorique de la diffusion Raman. Les aspects généraux de la présente thèse s'insèrent dans cette catégorie de recherche, à savoir la compréhension théorique des spectres Raman résonnants. Les trois systèmes qui seront étudiés dans ce travail sont : le graphène, les multicouches de graphène (few layer graphène ou FLG) et le graphite. Dans ces systèmes la diffusion Raman est toujours résonante à cause de leur caractère semi-métallique.

Plusieurs études théoriques et expérimentales ont montré que les propriétés électroniques et optiques des FLG (et leur version 3D, à savoir, le graphite) changent en fonction du nombre de couches de graphène qui les constituent et de la façon dont ces couches sont placées l'une sur l'autre, c'est-à-dire, du type d'empilement. L'empilement dans les FLG et dans le graphite, qui constitue le sujet principale de cette thèse, peut être du type Bernal (AB-stacking) ou du type rhomboédrique (ABC-stacking), comme le montre la Fig. R.B. L'empilement Bernal, étant énergétiquement plus favorisé, représente la forme la plus stable des FLG et du graphite. Cependant, dans les dernières années, les systèmes de multicouches de graphène en empilement rhomboédrique (RMG) ont attiré plus d'attention car des prédictions théoriques suggèrent l'apparition (à la surface de ces matériaux) d'une bande électronique non-dispersive (dispersion inférieure à 2 méV) à proximité du niveau de Fermi [START_REF] Kopnina | Surface superconductivity in rhombohedral graphite[END_REF][START_REF] Xiao | Density functional investigation of rhombohedral stacks of graphene: Topological surface states, nonlinear dielectric response, and bulk limit[END_REF]. Cette bande plate, dont la masse effective est très large, constitue un candidat très prometteur pour les états électroniques fortement corrélés tels que le magnétisme [START_REF] Pamuk | Magnetic gap opening in rhombohedral-stacked multilayer graphene from first principles[END_REF] ou la supraconductivité à température ambiante [START_REF] Precker | Identification of a possible superconducting transition above room temperature in natural graphite crystals[END_REF]. Ainsi, un outil de caractérisation qui soit à la fois précis et non-destructif est essentiel pour étudier la présence de l'empilement rhomboédrique dans les FLG et dans le graphite.

L'objectif principal de la présente thèse est la compréhension des spectres Raman résonants à deux phonons dans les FLG et dans le graphite pour tous les empilements possibles, ainsi que, l'extraction (à partir des calculs ab-initio) des signatures Raman de chaque type d'empilement. Ces signatures servent pour l'identification non-destructive et à la synthèse expérimentale des RMG et du graphite rhomboédrique à l'aide de la spectroscopie Raman.

Méthodologie

Afin d'extraire la signature Raman de chaque type d'empilement dans les FLG et dans le graphite, on s'est proposé de formuler un cadre théorique dans lequel on peut évaluer les sections efficaces de la diffusion Raman en tenant compte des effets de résonance. Notre approche sera formulée dans le cade de la théorie de perturbation à plusieurs corps (Many-Body Perturbation Theory ou MBPT). La seule difficulté de ce type d'approche, est que le temps de calcul nécessaire pour qu'une simulation d'un spectre Raman converge, est très long. En effet, Cette difficulté a deux origines différents :

Le premier, c'est le fait que les pics d'intérêt dans notre travail (les pics 2D) sont des pics à deux phonons. Tandis que dans la diffusion Raman du 1 er ordre, seul les phonons du centre de la zone de Brillouin (BZ) peuvent contribuer à la section efficace 124 Chapitre 5. Résumé en Français Raman, dans la diffusion à deux phonons, toute la BZ est sensée contribuer. La raison, est que dans le Raman du 2 ème ordre, la quantité qui est restreinte d'appartenir au centre de la BZ est la somme de l'impulsion des deux phonons. Cela, par contre, n'impose aucune restriction sur l'impulsion de chaque phonon (pris séparément) qui balaye maintenant toute la BZ.

La Deuxième origine, c'est le fait que le pic 2D est un pic résonant. Cela implique automatiquement l'utilisant de la MBPT pour l'évaluation des intensités Raman. Ce type d'approche, qui est considéré comme étant le plus précis pour la description des propriétés optiques des matériaux, demande, cependant, beaucoup plus de temps pour l'évaluation précise des états excités implqiués dans la diffusion Raman.

Plusieurs approches peuvent être utilisées pour faire face aux difficultées mentionnées. Cependant, tous ces approches sont ou bien très efficaces (rapides) mais alors avec une pauvre précision, ou bien, très précis mais avec simulation qui des demandent un temps infini pour converger. Notre approaches, par contre, est conçu dans le but d'avoir à la fois une bonne efficacité et une haute précision. La méthodologie qui a été développées dans ce travail de thèse pour contourner les difficultés mentionnées précédemment repose essentiellement sur l'utilisation des fonctions Wannier. Le calcule des intensités Raman en utilisant les fonctions de Wannier consiste, essentiellement, en l'interpolation de Wannier des éléments de matrices nécessaires pour l'évaluation de la section efficace Raman. Le grand avantage de l'utilisation de l'interpolation de Wannier à la place des autres méthodes d'interpolation déjà existantes est la précision ab-initio sur les quantités interpolées qui ne peut être atteinte par les anciennes méthodes d'interpolation. Il s'avère que le temps de calcul nécessaire pour interpoler les éléments de matrice avec les fonctions de Wannier est très réduit par rapport a un calcul pure ab-initio (c'est-à-dire sans aucune étape d'interpolation). Après avoir interpolé tous les éléments de matrice nécessaires, le calcul des intensités Raman sera effectué en utilisant la règle d'or de Fermi généralisée qui est, en effet, la même procédure utilisée dans les approches dit de liaison forte. De cette façon, nous aurons à la fois l'efficacité (la rapidité) d'un approche liaison-forte et la précision d'un calcul ab-initio.

En plus de l'utilisation de l'interpolation de Wannier, nous avons mis en place une deuxième méthode de réduction de temps de calcul. En fait, même avec l'utilisation de l'interpolation de Wannier, la simulation des spectres Raman pour les systèmes à 3 et 4 couches de graphène nécessitent toujours beaucoup de temps pour être évaluée. La situation est encore plus compliquée pour le graphite à cause du fait que ce dernier soit en 3 dimension au lieu de 2 comme c'est le cas des multicouches de graphène. La deuxième méthode de réduction sert, essentiellement, à éviter d'évaluer la section efficace Raman sur des régions de la BZ qui donne une contribution négligeable.

En fin, La formule qui sera utilisée pour l'évaluation des sections efficaces Raman, appelées I(ω), est donnée par : I(ω) ∝ 1 N q q,µ,ν |M µ,ν (q)| 2 δ(ω -ω q,µ -ω -q,ν )

× [n(ω q,µ ) + 1] × [n(ω -q,ν ) + 1] (R.1)

Où ω qµ et n(ω qµ ) sont (respectivement) les fréquences des phonons et la distribution de

Bose-Einstein pour le phonon d'impulsion q et appartenant à la branche µ . L'énergie du laser est désignée par ω L . La somme est effectuée sur une grille des impulsion des phonons (appelée grille des phonons) qui est uniforme dans la zone de Brillouin (BZ) où N q est le nombre total de points. Pour finir, les entiers positives i, j, l, m sont les indices des états électroniques (initiaux, intermédiaires et finaux) qui participent à la diffusion Raman. La probabilité de diffusion pour un couple spécifique de modes de phonon µ, ν est donnée par :

M µ,ν (q) = 1 N k 8 α=1 ijlm k M α,ijlm µ,ν (k, q) (R.2)
Où k représente le vecteur d'onde de l'électron et M α,ijlm mu,ν (k, q) avec α = 1, ..., 8 sont les probabilités de diffusion pour chaque diagramme de Goldstone (indexé par α) contribuant à la diffusion Raman. Les expressions explicites de ces probabilités sont données dans l'annexe B. La somme dans l'équation précédente est effectuée sur une grille électronique uniforme dans la BZ (appelée grille d'électrons) où N k est le nombre total de points.

Les grilles d'électrons et des phonons qui sont utilisées dans les deux équations précédentes sont les grilles réduites, c'est-à-dire, les grilles dans lesquelles les régions de la BZ qui ne contribuent pas à la section efficace Raman ont été éliminées.

Résultats et discussion

Application pour les tri-couches et quadri-couches de graphène. La conséquence la plus importante de l'accord obtenu pour les multicouches de graphène est que : les simulations Raman donnent pour chaque empilement sa propre signature Raman encodée dans la forme de son pic 2D. Cette signature, que nous avons calculé dans une approche ab-initio, peut être utilisée pour identifier chaque type empilement via son spectre Raman.

Application dans les graphites

Ayant validé notre calcul théorique dans le cas des tri-et quadri-couches, on passe maintenant au cas des graphites. Nous avons d'abord considéré le graphite-AB pour lequel plusieurs données expérimentales sont disponibles. Afin de trouver les signatures des différents types d'empilements dans le cas des graphites, nous avons effectué des calculs pour le cas du graphite-ABC et du graphite-ABCB. D'une part, le graphite-ABC est intéressant car il représente un empilement infini (3D) de la tri-couches de graphène en empilement rhomboédrique. D'autre part, le graphite-ABCB est aussi intéressant car il correspond à une séquence ... [ -D'autre part, même si les spectres obtenus pour le graphite-ABC et le graphite-ABCB semblent très similaires, ils diffèrent par la présence d'un pic à ∼ 2560 cm -1 qui est complètement absent dans l'empilement-ABCB. Comme la structure ABC diffère de celle en emipelemnt ABCB par l'apparition de longues séquences de couches de graphène en empilement ABC, la caractéristique à ≈ 2576 cm -1 peut être considérée comme la signature Raman d'une rhomboédricité à longue portée.

Le bon accord entre nos calculs théoriques et le spectre expérimental du Ref. [START_REF] Henni | Rhombohedral multilayer graphene: A magneto-raman scattering study[END_REF], qui concerne, à la fois, la largeur du pic 2D ainsi que la présence du nouveau pic à ≈ 2576 cm -1 suggère que, les échantillons du Ref. [START_REF] Henni | Rhombohedral multilayer graphene: A magneto-raman scattering study[END_REF] contiennent bien des séquences à longue portée de graphène en empilement rhombohedrique.

Conclusion

Dans ce travail de thèse, nous avons montré que la spectroscopie Raman est un outil très efficace pour la détermination du type d'empilement dans les multicouches de graphène ainsi que dans les graphites. En particulier, pour les multicouches de graphène, la forme du pic Raman 2D est une signature directe du type d'empilement étudié. Pour les graphites, nous avons pu distingué deux types de situations. La première concerne ce qu'on appelle la rhomboédricité a courte portée dans laquelle une courte séquence de couches de graphène en empilement rhomboédrique est présente à l'intérieure du graphite. Dans cette situation, nous avons montré que la largeur du pic 2D est un indicateur direct de ce type d'empilement. Dans le cas contraire, ou une longue séquence de couches de graphène en empilement rhomboédrique est présente, c'est-à-dire, une rhomboédricité à longue portée, un extra-pic apparaît autour de 2560 cm -1 comme une signature directe de ce type d'empilement. Ce pic est complètement absent dans le cas d'une rhomboédcirité a courte portée.

Abderrezak

  Figure B -Crystal structure of Bernal (ABA) and rhombohedral (ABC) stacked multilayer graphene. Red lines are guide for the eyes.

2 ChapterFigure 1 . 1 -

 211 Figure 1.1 -(a) The sp 2 hybridization of carbon atoms and formation of the σ bonding. (b) A piece of natural graphite. (c) Atomic structure of Bernal graphite showing the ABA stacking sequence and the hexagonal symmetry. (d) and (e) Typical STM images of graphene and few-layer graphene respectively (adapted from Ref. [21]). (f) High resolution STM image of one single carbon nanotube where the structure of the sp 2rings can be seen (adapted from Ref.[START_REF] Odom | Atomic structure and electronic properties of single-walled carbon nanotubes[END_REF]). (g) Atomic structure of carbon nanotubes with the three possible arrangements of the carbon rings (armchair, zigzag and chiral).

Figure 1 . 2 -

 12 Figure 1.2 -(a) Atomic structure of graphene showing the unit cell with two atoms and the first neighbor distance: 1.42 Å (b) First Brillouin zone in graphene. The points K and K , edges of the BZ are called Dirac points (c) Ab-initio (DFT-LDA) electronic bands structure of graphene showing the linear dispersion. (d) Ab-initio (DFT-LDA) electronic bands dispersion in bilayer graphene. (a) and (b) are adapted from [40].

Figure 1 . 3 -

 13 Figure 1.3 -(a) and (b) Ab-initio phonon dispersion of graphene and bilayer graphene respectively (technical details are given in Chap.3) . The "i" and "o" prefixes in labeling the phonon modes stand for "in-plane" and "out-of-plane" vibrations respectively. (c) Schematics for the Γ-point phonon displacement pattern for graphene and bilayer graphene, where each phonon mode in graphene gives to two phonon modes in bilayer. The Raman (R) and infrared (IR) activity is also indicated. Figure(c) is adapted from Ref.[START_REF] Reich | Raman spectroscopy of graphite[END_REF].

1. 3 . Raman scattering 13 Figure 1 . 5 -

 31315 Figure 1.5 -Schematic for a standard Raman spectrum obtained with Stokes and anti-Stokes processes.

Figure 1 . 6 -

 16 Figure 1.6 -Simple diagrammatic representation of Raman scattering according to perturbation theory.

Figure 1 . 7 -

 17 Figure 1.7 -(a) Raman spectrum of pristine and defected graphene showing all the Raman peak between 1000 cm -1 and 3300cm -1 . (b) double-Resonant Raman processes responsible of the G, D and 2D peaks in graphene. The black lines are the Dirac cones, red and blue arrows are optical transitions and dashed lines represent the scattering by a defect or a phonon. Adapted from [53].

Figure 1 . 8 -

 18 Figure 1.8 -(a) Splitting of the G Raman mode in graphene under strain.The splitting and the red shift of the G + and G -peaks increase with increasing strain ( ). (b) The splitting of the 2D Raman mode in strained graphene in the zigzag direction. The red and blue curves are for light polarizations parallel and perpendicular to strain direction respectively. Adapted from Refs.[START_REF] Mohiuddin | Uniaxial strain in graphene by raman spectroscopy: G peak splitting, grüneisen parameters, and sample orientation[END_REF][START_REF] Huang | Phonon softening and crystallographic orientation of strained graphene studied by raman spectroscopy[END_REF]

Figure 1 . 9 -

 19 Figure 1.9 -(a) variation of G peak position with respect to doping in graphene. Notice the stiffening of the G mode for both electron and hole doping. (b) The 2D position as a function of doping. The 2D peak responds differently to electron and hole doping. (c) The I(2D)/I(G) peaks ratio as a function of doping. Notice the high sensitivity of this ratio to doping. The solid blue lines are the predicted non-adiabatic trends from Ref. [49]. Adapted from Ref. [61].

Figure 1 .

 1 Figure 1.10 -(a) Dependence of the I D /I G ratio on the defect-defect distance L D for different Laser energy. We can see clearly the existence of two regimes : before and after L D ∼ 3 nm. (b) The width of the G peak in defected graphene as a function of defect concentration. Adapted from Ref. [65]

Figure 1 . 11 -

 111 Figure 1.11 -Variation of (a) the Raman 2D modes and (b) the shear mode in fewlayer graphene with the number of layers. Adapted from Ref. [53]. (c) and (d) are schematic representations for 2D Raman processes in graphene and bilayer graphene respectively.

Figure 1 .

 1 Figure 1.12 -(a) Bernal (AB) and rhombohedral (ABC) stackings. (b) The Raman 2D peak of Bernal and rhombohedral trilayer graphene at 1.96 eV. (b) is adapted from Ref. [69].

Figure 1 . 13 -

 113 Figure 1.13 -The resonant scattering process responsible for the 2D Raman peak in graphene.

Figure 2 . 1 -

 21 Figure 2.1 -Feynman diagrams and resonant Goldstone diagrams responsible for the two phonon Raman scattering. The text in front of each Goldstone diagram is a label to describe if the electron (e) or the hole (h) is the first to be scattered. Adapted from Ref. [19].

Figure 2 . 2 -

 22 Figure 2.2 -ARPES experimental electronic band dispersion of graphite compared to LDA and LDA+GW calculations at (a) H point (b) point between H and K points and (c) at the K point. Adapted from Ref. [91].

Figure 2 . 3 -

 23 Figure 2.3 -Equienergy contours of the photoemission intensity (logarithmic scale) around the K -H axis for k z = 0.22Å corresponding to Fig. 2.2 (c). Experimental data are denoted by crosses. Lines depict TB-GW (∆) and LDA ( ) calculations, respectively. Trigonal warping and the strong asymmetry in the intensity can clearly be seen. Adapted from Ref. [91].

  have been tabulated, it is possible to evaluate the matrix elements O (W ) nm ( k) at an arbitrary point k in the BZ by inverting Eq. 3.8 :

. 1 .

 1 In the Wannier interpolation presented above, this overlap is negligible due to the localization of Wannier functions which enforces the localization of the quantity O (W ) nm (R) and, hence, higher the accuracy of the interpolated quantity O

3 .Figure 3 . 1 -

 331 Figure 3.1 -Schematic overview of the Wannier interpolation procedure. The left panel shows the BZ mesh k used in the first-principless calculation, where the quantity of interest O(k) is explicitly calculated. The Wannier-transformed quantity O(R) is strongly localized near the origin of the equivalent supercell, shown in the middle panel with Wannier functions at the lattice points. The right panel shows a dense grid of interpolation points k in the BZ, where the quantity O(k ) is evaluated at low cost starting from O(R). Adapted from [100].

Fig. 3 .

 3 Fig. 3.3 shows a comparison between the Wannier interpolated LDA and LDA+GW electron bands for single and double layer graphene.

Figure 3 . 2 -Figure 3 . 3 -

 3233 Figure 3.2 -Upper panel (from left to the right) : Localization of the Hamiltonian, dipole , and deformation potential matrix elements as a function of the distance from the origin in real space in the case of graphene. Lower panel : The same study for bilayer graphene. In these plots we show the localization properties of H, γ and d for some band indexes. Matrix elements corresponding to different band indexes have the same localization behavior.

Figure 3 . 4 -

 34 Figure 3.4 -LDA versus LDA+GW Wannier interpolated phonon frequencies of graphene (left panel) and bilayer (right panel).

Figure 3 . 5 -

 35 Figure 3.5 -Resonant phonon regions (black) which contribute the most to the Raman cross section of the 2D mode in graphene (left panel) and bilayer graphene (right panel). In yellow, the full BZ sampling and k x,y are the Cartesian components of the crystal momentum.

Figure 3 . 7 -Figure 3 . 8 -

 3738 Figure 3.7 -Schematic for the general procedure of calculating the Raman spectrum.

3. 4 .

 4 Tests in single and double layer graphene. 65 2620 2640 2660 2680 2700 2720 2740 2760

Figure 4 . 1 -

 41 Figure 4.1 -Crystal structure of Bernal (ABA) and rhombohedral (ABC) stacked multilayer graphene. Red lines are guide for the eyes.

Figure 4 . 2 -

 42 Figure 4.2 -Wannier versus DFT-LDA electronic bands for (a) ABA and (b) ABC trilayer graphene. The horizontal line indicates the Fermi level.

69 Γ<

 69

Figure 4 . 3 -

 43 Figure 4.3 -LDA versus LDA+GW Wannier interpolated electronic bands for (a) ABA and (b) ABC trilayer graphene.

Figs. 4 .Figure 4 . 4 -Figure 4 . 5 -

 44445 Figs.4.5 and 4.6 show the GW-corrected phonon frequencies compared to LDA calculations for all possible stackings in trilayer and tetralayer graphene respectively. A careful analysis of the difference between the corrected phonon branches in ABA and

Figure 4 . 6 -

 46 Figure 4.6 -Comparison between LDA and LDA+GW-corrected phonon dispersion for the three polytypes of tetralayer graphene, namely : ABAB, ABCA and ABCB stacking around the K point in the BZ.

Figure 4 . 9 -

 49 Figure 4.9 -Experimental (black) and theoretical (red) difference between ABA and ABC 2D Raman peaks at different laser energies. The area of the theoretical and experimental spectra was normalized to 1. The difference in intensity ( vertical axis) is reported with respect the maximum intensity of the experimental spectra at 2.41 eV of the ABA-stacked trilayer and expressed in term of percentage.

Figure 4 .

 4 Figure 4.10 -(a) The 5 allowed optical transitions in trilayer graphene near the K point (figures adapted from Ref. [40]). (b) Cartoon of symmetric, anti-symmetric, inner (red) and outer (blue) double resonant Raman processes in trilayer graphene.

Figure 4 .

 4 Figure 4.11 -(a) Zoom of the BZ around the K point. (b) Outer process involving electrons on the K-Γ line and phonons on the K -M line. (c) Inner process involving electrons on the K -M line and phonons on the K -Γ line. (a) is adapted from Ref. [68]. (b) and (c) are adopted from Ref. [115].

Figure 4 . 12 -

 412 Figure 4.12 -Analysis of the 2D Raman mode in ABA (upper panel) and ABC (lower panel) trilayer graphene at 1.96 eV. (a) polarization analysis, (b) e-h contribution and inner outer decomposition, (c) process analysis: dashed and solid lines with the same color indicate degenerate processes. (d) symmetric and anti-symmetric decomposition. See Sec. 4.5 for relevant definitions.

Figure 4 . 13 -

 413 Figure 4.13 -Analysis of the 2D Raman mode in ABA (upper panel) and ABC (lower panel) trilayer graphene at 2.71 eV. (a) polarization analysis, (b) e-h contribution and inner outer decomposition, (c) process analysis (d) symmetric and anti-symmetric decomposition. See Sec. 4.5 for relevant definitions.

Figure 4 .

 4 Figure 4.15 -(a) Contributions of the P 55 2 and P 44 3 processes the 2D peak of ABA trilayer graphene at 1.96 eV and their decomposition upon the TO phonon branches. (b) Comparison between experimental and theoretical (LDA) 1 st and 3 rd TO phonon dispersions in ABA trilayer graphene. The blue and green connected dots are the extracted phonon frequencies from the experimental Raman spectra used in this work.

Figure 4 .

 4 Figure 4.16 -(a) Bands structure for trilayer graphene in ABA and ABC stackings. (b) 2D contour plot of the resonant phonon wave vector for process P 55 2 . The external ring is for ABC trilayer graphene and the internal one (which coincides with the external one for the outer direction) is for ABA stacking. In (b) the dashed lines indicate the inner directions.

Figure 5 . 1 -

 51 Figure 5.1 -Crystal structure of Bernal (ABA) and rhombohedral (ABC) stacked multilayer graphene. Red lines are drawn as a guide for the eyes.

  .[ABC](BAB)[CBA](BCB)... that is an equal mixing of trilayer with rhombohedral (labeled [ABC] or [CBA]) and Bernal (labeled (BAB) or (BCB)) stackings.

Figure 5 . 2 -Figure 5 . 3 -

 5253 Figure 5.2 -Left panel : Ab-initio (LDA) versus Wannier interpolated electronic bands for (a) AB and (b) ABC, and (c) ABCB graphites. Right panel : LDA versus LDA+GW interpolated electronic bands for (a) AB and (b) ABC, and (c) ABCB graphites. The Fermi level is set to 0.

Figure 5 . 6 -

 56 Figure 5.6 -Phonon momenta contributing to the Raman cross section in bulk AB (a) and bulk ABC (b) graphite around the K point in the BZ at 1.96 eV. The solid and dashed white lines denote the K -M and K -Γ high-symmetry lines respectively.The color bar indicates the normalized q-resolved Raman cross section, where q is the phonon momentum.

Figure 5 . 7 -

 57 Figure 5.7 -(a) Bernal (bulk AB) versus rhombohedral (bulk ABC) electronic bands dispersion. To obtain the gray (red) area, we have calculated the electronic bands dispersion for bulk AB (bulk ABC) graphite along the high symmetry lines (0, 0,k z ) → (1/3, 1/3, k z ) → (1/2, 1/2, k z ) (crystal coordinates) for different values of k z .The obtained band structures were then put on top of each other in the same plot. This gives the grey area for AB graphite and the red area for ABC graphite. (b) Contribution of the traversal optical (TO) modes (dashed lines) to the 2D Raman mode for bulk AB and bulk ABC graphites (solid lines) (c) The same as in panel (b) except that the TO contribution is now calculated with constant electron-phonon matrix elements.

Figure 5 . 8 -Fig. 5 .

 585 Figure 5.8 -Analysis of the 2D Raman mode of bulk AB-stacked graphite (Bernal graphite) at 1.96 eV. (a) polarization analysis, (b) e -h contribution and inner outer decomposition, (c) process analysis, (d) symmetric and anti-symmetric decomposition.

Fig. 5 .

 5 Fig.5.8 (b) shows that the e-h scattering is the dominant contribution to the 2D-mode over the e-e and h-h processes. However, while in FLG e-e and h-h processes are negligible, in Bernal graphite they have a small but not negligible contribution. Analysis at different laser energies shows that, as for the case of trilayer and tetralayer graphene, the inner contribution is dominant at low energy. The outer contribution

5 . 5 K

 55 The definition of these contributions has been already introduced in the previous chapter, Sec. 4.Γ<

Figure 5 . 9 -

 59 Figure5.9 -The extracted symmetric TO phonon dispersion in AB-graphite from experimental Raman spectra compared to LDA and LDA+GW calculations. We have used three laser energies only : 1.96 eV, 2.33 eV, and 2.41 eV.

Figure D. 4 -

 4 Figure D.4 -Analysis of the 2D Raman mode in ABA (upper panel) and ABC (lower panel) trilayer graphene at 2.54 eV. (a) polarization analysis, (b) e-h contribution and inner outer decomposition, (c) process analysis: dashed and solid lines with the same color indicate degenerate processes. (d) symmetric and anti-symmetric decomposition. See Sec. 4.5 for relevant definitions.

Figure D. 5 -

 5 Figure D.5 -Analysis of the 2D Raman mode of bulk AB-stacked graphite (Bernal graphite) at 2.41 eV. (a) polarization analysis, (b) e -h contribution and inner outer decomposition, (c) process analysis, (d) symmetric and anti-symmetric decomposition.

Figure 121 Figure

 121 Figure R.A -Schéma simplifié montrant les processus optiques linéaires à l'intérieur d'un matériau. D'après Ref [2].

  Tous d'abord, nous avons calculé les spectres Raman résonants (le pic 2D principalement) pour trois et quatre couches de graphène avec tous les empilements possibles. Pour ces systèmes, une détermination de la séquence d'empilement a été obtenue auparavant par des mesures optiques[START_REF] Lui | Imaging stacking order in few-layer graphene[END_REF]. Les résultats de nos calcules sont présentés dans la Fig. R.C pour les tri-couches et dans Fig. R.D pour les quadri-couches en utilisant plusieurs énergies du laser et en comparaison avec les données expérimentals. L'accord obtenu entre nos calculs ab-initio et les données expérimentales déjà existantes est très bon. En particulier, nous avons reproduit toutes les caractéristiques spectrales de position, largeur et intensité ainsi que la dépendance en fonction de l'énergie du laser des pics observés expérimentalement.

  Le graphite-AB (ou le graphite Bernal) est l'empilement inifini de type AB (voir Fig. R.B) des couches de graphène. Nous avons calculé les spectres Raman résonants (toujours le pic 2D) pour différentes énergies du laser et nous avons trouvé un excellent accord avec les données expérimentales (voir Fig. R.E). Ensuite, nous avons analysé le spectre à ω L = 1.96 eV de façons détaillée. Le pic 2D du graphite se compose d'un pic principal à ≈ 2683 cm -1 et une épaule autour de 2640 cm -1 , comme le montre la Fig. R.E. Les deux caractéristiques sont bien décrites par nos calculs. De plus, la forme et l'intensité de la structure D + D à ≈ 2456 cm -1 , sont aussi bien décrites. Cependant, pour ce pic un faible décalage en énergie existe entre nos calculs et les données expérimentales.

Figure

  Figure R.F -Les spectres Raman théoriques du graphite AB, ABC et ABCB comparés aux données expérimentales des échantilons du Réf.[START_REF] Henni | Rhombohedral multilayer graphene: A magneto-raman scattering study[END_REF] à 1,96 eV. Dans (a) la signature de l'empilement ABC à longue portée à environ ∼ 2576 cm -1 est montrée. Dans (b), nous montrons l'absence de cette signature pour la rhombohédricité a courte porté dans le graphite-ABCB. Les données expérimentales proviennent d'échantillons composés de 17 couches de graphène en empilement ABC[START_REF] Henni | Rhombohedral multilayer graphene: A magneto-raman scattering study[END_REF].

Table 1 .

 1

	12)

1 -Useful Mulliken notations for labeling phonon modes in crystals according to their symmetry properties.

  (a.u.) Measured versus calculated Raman spectra for the three possible stacking in tetralayer graphene. Experimental data are from Ref.[START_REF] Lui | Imaging stacking order in few-layer graphene[END_REF].

	4.4. Calculated spectra			75
									Tetralayer
		2.41 eV					Trilayer ABA 2.41 eV	2.41 eV
		Exp. Theory 1.96 eV Intensity (a.u.) ABAB Exp. Theory				2.09 eV ABCA	2.33 eV ABCB
		2550	2700		2850	2550	2700	2850	2550	2700	2850
			Raman shift (cm	-1 )	2500 2600 2700 2800 Raman shift (cm -1 )	2500 2600 2700 2800 Raman shift (cm	-1 )
	Intensity (a.u.)	2.41 eV 2.54 eV Exp. Theory Raman intensity (a.u.) ABAB				2.54 eV 2.54 eV ABCA	2.71 eV 2.54 eV ABCB
		2500 2600 2700 2800 2900	2500 2600 2700 2800 2900	2500 2600 2700 2800 2900
		2500 2600 2700 2800 Raman shift (cm -1 )	2500 2600 2700 2800 Raman shift (cm -1 )	2500 2600 2700 2800 Raman shift (cm -1 )
		Raman shift (cm -1 )			Raman shift (cm -1 )	Raman shift (cm	-1 )
	Intensity (a.u.) Figure 4.8 --28% Exp. Theory 1.96 eV + 28 % 0 Intensity (a.u.)			Trilayer ABC 2.09 eV 1.96 eV exp theo	2.33 eV	2.33 eV	2.54 eV
		2500 2600 2700 2800 + 28 %	2500 2600 2700 2800	2500 2600 2700 2800
	Intensity (a.u.)	2.41 eV	-28% 0 Intensity (a.u.)	2.54 eV 2.09 eV 2600 2700 2800	2.41 eV Raman shift (cm ) -1 2600 2700 2800	2.71 eV 2600 2700 2800 2.71 eV
		2500 2600 2700 2800	2500 2600 2700 2800	2500 2600 2700 2800
		Raman shift (cm -1 )			Raman shift (cm -1 )	Raman shift (cm	-1 )

Figure 4.7 -Measured versus calculated Raman spectra of ABA and ABC trilayer for different laser energies. Experimental data are from Refs. [112, 69].

  .14, one can attribute each resonant region (i.e. each ring) from outside to inside the rings to the following processes : P44 3

	, P 45 3 + P 54 2 , P 46 3 + P 64 1 , P 55 2 , P 56 2 + P 65 2 , P 66 1

  Figure 5.4 -Comparison between the calculated 2D Raman mode for bulk AB graphite (Bernal graphite) and the experimental Raman spectra obtained from HOPG graphite at different laser energies. Experimental data are from Ref.[START_REF] Zhang | Review on the raman spectroscopy of different types of layered materials[END_REF] for 1.96 eV and 2.33 eV and from Ref.[START_REF] Ferrari | Raman spectrum of graphene and graphene layers[END_REF] for 2.41 eV.

	5.3. The Raman fingerprint of rhombohedral graphite	93
						Bulk AB
	Intensity (a.u.)	Exp. Theory 1.96 eV			2.33 eV				2.41 eV
		2400	2600	2800	2400	2600	2800	2400	2600	2800
		Raman shift (cm	-1 )	Raman shift (cm	-1 )	Raman shift (cm	-1 )
	Raman intensity (a.u.)	Exp. 17 layers Theory bulk ABC Theory bulk AB ω L = 1.96 eV			Exp. 17 layers Theory bulk ABCB ω L = 1.96 eV
		2400		2550	2700	2850		2400	2550	2700	2850
			Raman shift (cm -1 )					Raman shift (cm -1 )

  Figure R.C -Les spectres Raman mesurés et calculés des tricouches ABA et ABC pour différentes énergies du laser. Les données expérimentales proviennent des Refs. [112, 69].Figure R.D -Comparaison entre les spectres Raman mesurés et calculés pour les trois empilements possibles dans la quadri-couches de graphène. Les données expérimentales proviennent de la Réf. [69]. Pour l'empilement ABCB, on n'a pu trouvé des données expérimentales de sont spectre Raman, c'est pour cela, seule les spectres théoriques ont été montrés.

										127
						Tetralayer
		2.41 eV				Trilayer ABA 2.41 eV			2.41 eV
	Intensity (a.u.)	Exp. Theory 1.96 eV Intensity (a.u.) ABAB Exp. Theory			2.09 eV ABCA			2.33 eV ABCB
		2550	2700	2850	2550	2700	2850	2550	2700	2850
		2500 2600 2700 2800 Raman shift (cm -1 )	2500 2600 2700 2800 Raman shift (cm -1 )	2500 2600 2700 2800 Raman shift (cm	-1 )
	Intensity (a.u.)	2.41 eV 2.54 eV Exp. Theory Raman intensity (a.u.) ABAB			2.54 eV 2.54 eV ABCA			2.71 eV 2.54 eV ABCB
		2500 2600 2700 2800 2900	2500 2600 2700 2800 2900	2500 2600 2700 2800 2900
		2500 2600 2700 2800 Raman shift (cm -1 )	2500 2600 2700 2800 Raman shift (cm -1 )	2500 2600 2700 2800 Raman shift (cm -1 )
		Raman shift (cm -1 )			Raman shift (cm -1 )			Raman shift (cm	-1 )
						Trilayer ABC
		1.96 eV			2.09 eV			2.33 eV
	Intensity (a.u.)	Exp. Theory
		2500 2600 2700 2800	2500 2600 2700 2800		2500 2600 2700 2800
		2.41 eV			2.54 eV			2.71 eV
	Intensity (a.u.)
		2500 2600 2700 2800	2500 2600 2700 2800		2500 2600 2700 2800
		Raman shift (cm -1 )			Raman shift (cm -1 )			Raman shift (cm	-1 )

  ABC] Figure R.E -Comparaison entre les pics Raman 2D calculés et les données expérimentles pour le graphite-AB (Bernal) à différentes énergies du laser. Les données expérimentales proviennent de la Réf. [124] pour 1.96 eV et 2.33 eV, et de la Réf. [13] pour 2.41 eV.

						Bulk AB
	Intensity (a.u.)	Exp. Theory 1.96 eV			2.33 eV				2.41 eV
		2400	2600	2800	2400	2600	2800	2400	2600	2800
		Raman shift (cm	-1 )	Raman shift (cm	-1 )	Raman shift (cm	-1 )
	Raman intensity (a.u.)	Exp. 17 layers Theory bulk ABC Theory bulk AB ω L = 1.96 eV			Exp. 17 layers Theory bulk ABCB ω L = 1.96 eV
		2400		2550	2700	2850		2400	2550	2700	2850
			Raman shift (cm -1 )					Raman shift (cm -1 )

TORCHE 05/10/2017 Sujet : Simulations ab-initio des spectres Raman résonants dans le graphène, les multicouches de graphène et le graphite Résumé

  : Les multicouches de graphène en empilement rhomboédrique sont considérés comme une phase prometteuse du carbone. Cela est due à la particularité de cette phase de pouvoir exhiber des états à forte corrélation électronique comme le magnétisme ou la supraconductivité à haute température critique. Ce qui est due, a son tour, à l'occurrence d'un état de surface avec une dispersion d'énergie électroniques quasi-nulle à proximité du niveau de Fermi. Malgré que le graphite Bernal soit la forme la plus stable du graphite, des échantillons a trois et quatre couches de graphène en empilement rhomboédrique ont pu être synthétisés. Plus récemment, des flocons d'épaisseur dépassant les 17 couches ont été isolés et provisoirement attribués à des séquences d'empilement rhomboédrique. Cette attribution à été faite via des expériences de spectroscopie Raman sous champ magnétique, bien que l'empreinte Raman des multicouche de graphène en empilement rhomboédrique est actuellement inconnue. Même le cas simple du spectre Raman résonnant à deux phonons (le pic 2D) du graphite Bernal n'est pas totalement compris. Dans ce travail de thèse, nous fournissons une description ab-initio complète du pic Raman 2D dans les systèmes de graphène à trois et quatre couches pour tous les empilements possibles, ainsi que pour le graphite Bernal, rhomboédrique et une alternance de graphite Bernal et rhomboédrique. Les calculs ont été effectués pour plusieurs énergie du laser et nous donnons à la fin des prescriptions pratiques pour l'identification des multicouches de graphène susceptible d'avoir un fort empilement rhomboédrique.

Mots clés : Multicouches de graphène, empilement rhomboédrique, empilement Bernal, spectroscopie Raman.

The first order in the Taylor series being neglected assuming that the atoms vibrate around their equilibrium positions.

Decoupled in the adiabatic sense, as discussed here.

scf stands for self-consistent field.

In the formulation above, we dropped the parametric dependence on the nuclear positions for simplicity.

We neglect the first order term in the harmonic approximation assuming that the nuclei vibrate around their equilibrium positions.

This formula was derived for graphene. However, since we are interested in phonons which are an in-plane vibration, the formula is also valid for graphite.

This procedure do not require the knowledge of phonon symmetries in ABC trilayer.

The symmetry x ≡ y means that the Raman cross section is invariant under the permutation x → y and y → x. This is equivalent to say that the Raman tensor is a symmetric tensor. This latter property can be deduced from the invariance of the Raman tensor under 120 • rotations from the hexagonal symmetry.

We label the TO sub-branches by 1 st , 2 nd , . . . from low to high energy near the K point in the BZ.

We take half of the frequency because the position of the 2D peak is the sum of frequencies of the two involved phonons during the scattering.

This is true only if AA stacking is ignored

Some authors[START_REF] Charlier | First-principles study of the stacking effect on the electronic properties of graphite(s)[END_REF] suggest a slightly different interplaner distance of

3.34 Å for rhombohedral graphite. For the sake of comparison in this work, we have taken it the same as for Bernal graphite.3. Recall that the number of Wannier functions in FLG and graphite systems is equal to the number of atoms per unit cell. This is because of the pure π character of the bands near the Fermi level.

Acknowledgments

I would first like to thank my supervisors : M. Matteo Calandra, M. Jean-Christophe Charlier, and M. Francesco Mauri for their patience and support of my Ph.D, for sharing their knowledge and experience with me, and for their efforts to make this collaboration a nice and fruitful experience. I also would like to thank Pr. Pedro Venezuela for his insightful comments and encouragement.

Chapter 3

Calculation of Raman intensities using the Wannier interpolation scheme and the BZ reduction method

This chapter is organized as follows. First, the theoretical formalism of maximally localized Wannier functions is presented. Second, The Wannier interpolation scheme and it applications to matrix elements needed in Raman simulations is given. Single and double layer graphene are taken as a case study. Third, the Brillouin zone (BZ) reduction method is presented and applied to graphene and bilayer graphene. Finally, tests of our approach to calculate Raman intensities in single and double layer graphene are presented.

Wannier interpolation

The need of powerful interpolation schemes is highly requested in many computational areas. This is mainly due to the fact that first-principles calculations of some fundamental properties, like phonon frequencies or the electron-phonon coupling, are very expensive and mostly carried out in a coarse sampling of the Brillouin zone (BZ). However, the convergence of many physical quantities and observables (e.g. Raman cross sections, phonon self-energy, etc ...) requires the use of much finer grids, sometimes as many as 10 6 wave vector in the BZ are needed [START_REF] Giustino | Electron-phonon interaction using wannier functions[END_REF].

The Wannier interpolation technique can provide interpolated quantities with an ab-initio precision.

The key property of Wannier functions that allows the interpolation of the ab-initio quantities with essentially no loss of accuracy is localization. This property constitutes the central property studied in this chapter. It turns out that, this localization is of great usefulness as it controls the precision of the interpolated ab-initio quantities. Here, we explore this property for the purpose of interpolating the matrix elements needed in Raman cross section calculations [START_REF] Marzari | Maximally localized generalized wannier functions for composite energy bands[END_REF][START_REF] Souza | Maximally localized wannier functions for entangled energy bands[END_REF].

ABC ABA

Chapter A

Inclusion of electron correlations effects.

In Chap. 2 we have discussed the importance of electron correlation in calculating Raman intensities. Here we present the practical way we included them in our calculation. The case of electrons is straightforward, since the renormalization of the Fermi velocity (+18%) is relatively simple. After the Wannier interpolation of the Hamiltonian operator, one re-scales the the Kohn-Sham energies as follows :

where k is the interpolated electron momentum and f is the Fermi level.

In the case of phonons, only the non-analytic term (phonon self-energy) is renormalized due to electron correlations. Calculations done in Ref. [START_REF] Lazzeri | Impact of the electronelectron correlation on phonon dispersion: Failure of LDA and GGA DFT functionals in graphene and graphite[END_REF] for graphene and graphite suggest a rescaling factor of r = 1.61 at the vicinity of the K point in the BZ. Namely, near the K point we have :

In Ref. [START_REF] Lazzeri | Impact of the electronelectron correlation on phonon dispersion: Failure of LDA and GGA DFT functionals in graphene and graphite[END_REF] this rescaling factor is a constant and the phonons are studied just in the neighborhood of K. Here, following Refs. [START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF][START_REF] Herziger | Two-dimensional analysis of the double-resonant 2D raman mode in bilayer graphene[END_REF] and in order to obtain a phonon dispersion all over the BZ, the rescaling factor r must depend on q such that r(q) = r near K and smoothly drops to one elsewhere [START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF][START_REF] Herziger | Two-dimensional analysis of the double-resonant 2D raman mode in bilayer graphene[END_REF] :

where "erfc" is the error function, a 0 = 2.46 Å the graphene lattice constant, and K n is the nearest vector to q among those equivalent to K. This formula simply expands the effect of electron correlations around the K points. The above equation was tested in the case of graphene [START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF] and bilayer graphene [START_REF] Herziger | Two-dimensional analysis of the double-resonant 2D raman mode in bilayer graphene[END_REF] and gave the correct phonon dispersion compared to experiments.

Chapter B Goldstone diagrams for double resonant Raman scattering

In this Appendix we give the explicit expressions of the eight Goldstone diagrams that have been used to evaluate Raman cross sections and already introduced in chap.

(see Fig. B.1). The e-e (h-h) label refers to diagrams where only the electron (hole)

is scattered twice (the so called a-a processes [START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF]). Conversely, e-h (h-e) refers to diagrams where the electron then the hole (the hole then the electron) are scattered (the so-called a-b processes [START_REF] Venezuela | Theory of double-resonant raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands[END_REF]).

In the expressions below : k and q are the electron and phonon wave vectors respectively. The states |kπ i,j,l,m where i, j, l, m are positive integers, refer to electronic states that are occupied by the electron (or the hole) during the Raman scattering. H in el and H out el are the incoming and the outgoing dipole operators while H ep,qν is the electron-phonon coupling operator describing a scattering with a phonon of momentum q from the ν branch. In the denominators, the i,j,l,m k are the electron band energies corresponding to the states |π i,j,l,m and δ is the sum of the electron and hole lifetimes. The laser energy is indicated by L . For the sake of simplicity, the electronic band and the phonon branch indexes as well as the phonon and electron momentum were dropped from the right hand side of the equations below.

Diagram ee1: the electron is first scattered by the qν phonon and then by the -qµ one,

Diagram ee2: the electron is first scattered by the -qµ phonon and then by the qν one, one,

.

Diagram hh2: the hole is first scattered by the -qµ phonon and then by the qν one,

first the electron is scattered by the qν phonon and then the hole by the -qµ one,

: first the electron is scattered by the -qµ phonon and then the hole by the qν one,

Diagram he1: first the hole is scattered by the qν phonon and then the electron by the -qµ one,

Diagram he2: first the hole is scattered by the -qµ phonon and then the electron by the qν one,

Chapter C

TO-character calculation

In order to reveal which phonon branches are responsible for the 2D and D+D Raman peaks in AB and ABC-graphites, we have calculated the TO phonon contribution to the Raman cross section. In order to perform this kind analysis, we firstly have to define, for a given phonon, the exact meaning of its TO-character. To do that, one could simply think about ordering all eigenmodes of the dynamical matrix near the K point in the BZ from low to high energy and taking the highest optical branches to be the TO phonons. This is correct, since we know that at the K point the highest optical branch in AB and ABC-graphites is a TO vibration [START_REF] Lazzeri | Impact of the electronelectron correlation on phonon dispersion: Failure of LDA and GGA DFT functionals in graphene and graphite[END_REF]. However, due to the occurrence of several crossing between the phonons branches in the whole BZ, this definition is only correct at the vicinity of the K point. Alternatively, if we examine the behavior of the non-analytic (NA) term of the dynamical matrix, we see that the eigenmodes (eigenvectors) of this term do not cross in a wide region around the K point, meaning that, the order of a TO phonons do not change around K (see Refs. [START_REF] Piscanec | Kohn anomalies and electron-phonon interactions in graphite[END_REF][START_REF] Lazzeri | Impact of the electronelectron correlation on phonon dispersion: Failure of LDA and GGA DFT functionals in graphene and graphite[END_REF] to examine the NA term behavior around K). This allow us to identify the TO eigenmodes of the NA term in a wide region around K. We find that, these phonons are : the 2 highest (3 highest) optical modes in AB-stacked (ABC-stacked) graphite. Once the TO eigenmode for the NA term has been identified, we proceed as follows. A phonon which corresponds to an eigenmode e qν of the dynamical matrix at wave vector q has a TO-character C TO qν defined by :

where the dot indicates a complex scalar product. The index m is the number of the TO modes of the NA term (m = 2 for AB-graphite and m = 3 for ABC-graphite) and the quantity e NA qµ is the eigenmode of the NA term corresponding to the TO phonon mode µ. In the expression above, we have assumed the eigenmode of the dynamical matrix and the NA term are normalized to 1. We defined the TO-contribution to the Raman cross section I TO (ω) by the sum of all Raman processes weighted by the TO character of the involved two phonons. Namely, we have :

where M µν (q) is the contribution to the scattering by the two phonons from the branches µ and ν at wave vector q (See Eq. 2.47). The delta function is for energy conservation and ω, and ω qν are, respectively, the Raman shift and the phonon frequencies Chapter D

Spectra analysis at higher excitation energies

In this Appendix, we report the analysis that has been performed for ABA and ABC-trilayer graphene at laser energies : 2.33 eV, 2.09 eV, 2.41 eV, and 2.54 eV. We also show the analysis performed at 2.41 eV for AB-graphite (Bernal graphite) Raman spectrum. Here we provide a complete first principles description of the 2D Raman peak in three and four layer graphene for all possible stackings, as well as for bulk Bernal, rhombohedral and an alternation of Bernal and rhombohedral graphite, that can be seen as a periodic sequence of ABA and ABC trilayers. Calculations for several laser energies are performed and we give practical prescriptions are proposed to identify long range sequences of ABC multi-layer graphene flakes. Our work is a prerequisite to experimental nondestructive identification and synthesis of rhombohedral multi-layer graphene and graphite.

ABA-trilayer

Keywords : Few layer graphene, rhombohedral stacking, Bernal stacking, Raman spectroscopy.