DFT Density Functional Theory TDDFT Time-Dependent Density Functional Theory EMT Effective Medium Theory IEMT Inverse Effective Medium Theory EEL/EELS Electron Energy Loss Spectroscopy FT Fourier Transform GGA Generalized Gradient Approximation HF Hartree-Fock HSE Heyd-Scuseria-Ernzehof IPA Independent Particle Approximation KS Kohn-Sham (approximation) LDA Local Density Approximation LF Local Field NLF No Local Field LR Long Range PBE Perdew, Burke and Ernzerhof (functional) RPA Random Phase Approximation

SR Short Range SC Super Cell xc exchange-correlation || component, parallel to the slab/surface ⊥ component, perpendicular to the slab/surface p-parameter npwmat_z = 2

Résumé

Le but de la thèse était de simuler le spectre d'absorption de méta-matériaux pour les applications photovoltaïques. Par méta-matériaux, nous entendons une assemblée d'objets de taille nanométrique situés à distance mésoscopique. L'idée sous-jacente est qu'en modifiant la taille du nano-objet et l'arrangement géométrique, on peut ajuster le seuil d'absorption. Pour calculer ces quantités, j'ai utilisé l'état de l'art du formalisme, c'est-à-dire des méthodes ab initio.

La première étape du travail a été dédiée au calcul de l'absorption d'un objet isolé (tranche de silicium, graphène, hBN). Dans le cadre de codes périodiques, on utilise une supercellule avec du vide pour isoler l'objet, et une méthode a été développée précédemment dans le groupe de Spectroscopie Théorique du LSI, pour obtenir des résultats indépendants du vide. Elle est appelée Selected-G, et a été appliquée avec succcès aux surfaces de silicium. Pour une tranche isolée, une expression modifiée du potentiel coulombien dans l'espace réciproque, appelé "slab potential", doit être utilisée. Pour valider l'utilisation du potentiel de slab pour le calcul de la matrice diélectrique microscopique, j'ai simulé les spectres de perte d'énergie d'électrons pour des empilements de quelques plans de graphène, et reproduit avec succès les données expérimentales disponibles. Cela a offert la possibilité d'étudier la dispersion du plasmon d'un plan de graphène, et discuter la nature des excitations électroniques dans ce système (transitions interband ou plasmon 2D).

La second étape a été consacrée à l'étude du spectre d'absorption d'une assemblée de tranches en interaction. Comme il a été mis en évidence que le formalisme de supercellule agit comme une théorie de matériau moyen avec du vide, avec l'effet erroné d'avoir des spectres dépendant de la taille de la supercellule, j'ai renversé la procédure pour extraire le spectre de la tranche en interaction, affranchi du problème du vide. La faisabilité a été démontrée sur les tranches de hBN, dont le caractère semi-conducteur à large bande interdite évite les instabilités numériques. Cela a permis de comprendre la raison pour laquelle l'absorption de la tranche en interaction de silicium apparaît à plus basse énergie que celle du matériau massif : cela vient de la présence des états de surface dans la bande interdite de la structure de bandes du massif. Néanmoins, la différence avec la tranche isolée doit être encore étudiée.

La troisième partie a été dédiée à l'étude de matériaux utilisés, ou candidats, aux applications photovoltaïques comme InP et InSe. J'ai étudié dans un premier temps les structures de bandes des massifs. Pour corriger la sous-estimation de la bande interdite calculée dans l'approximation de la densité locale (LDA), j'ai calculé les corrections GW, et utilisé la fonctionnelle d'échange et corrélation de Heyd-Scuseria-Ernzerhof (HSE). Le spectre d'absorption de InP massif a été calculé en résolvant l'équation de Bethe-Salpeter, qui permet de tenir compte des effets excitoniques. Comme ce calcul est très lourd numériquement, j'ai également comparé avec le calcul beaucoup plus léger de TDDFT avec le kernel à longue portée pour introduire les effets excitoniques. Pour le massif de InSe, j'ai calculé les corrections HSE pour les valeurs propres et obtenus un bon accord avec la bande interdite expérimentale. Les spectres obtenus en TDDFT, avec le kernel à longue portée, donne de bons résultats. J'ai commencé l'étude de tranches de ces deux matériaux. Des couches épaisses de InP et InSe ont été considérées et une reconstruction de surface (2x2) a été réalisée sur InP pour obtenir une surface semi-conductrice. La structure de bande LDA et les spectres d'absorption ont été calculés. Comme des systèmes d'une telle taille sont hors de portée des calculs de corrections

Introduction

In the last two hundreds of years, with developing of manufacturing and industrial revolution, humanity was particularly keenly aware of the need for new sources of energy. With the passing of the time, humanity is trying to convert environment around into energy with increasing efficiency. The first energy source called fire was used about 1 million years BC [START_REF] Berna | Microstratigraphic evidence of in situ fire in the Acheulean strata of Wonderwerk Cave, Northern Cape province, South Africa[END_REF]. Ancestors of people of our time converted wood into fire and heat which gave ability to survive in a cold regions or during glacial periods which have a periodicity of appearance. There was a great way from the first fire sources to the modern energy sources.

Nowadays people increased significantly their knowledge about energy sources and use different materials to produce the energy. The well-known energy sources are coal and oil, last tens of years nuclear energy has been greatly developed. But all these sources are not safe for the environment. Coal burning emits a lot of soot and smoke in the atmosphere, the petrol, the fraction of oil, contains heavy metal ions and this acts very negatively on the environment. From the other side, although coal and oil deposits are huge, they are finite and at certain moment inevitably will be fully consumed. Returning to the nuclear energy, in the recent history, there are two catastrophes namely Chernobyl and Fukushima. All of the spoken above pays attention to the fact that it is necessary to turn to other, less dangerous, less polluting and easily accessible, energy sources.

The solution is renewable energy sources. The well-known examples of a renewable energy are wind power, hydro power, geothermal energy, biomass energy and solar energy. In the thesis the attention will be focused on the solar cells energy and on the discipline called photovoltaics which learn how to transform the energy of the sun into the electricity.

In 1839, French physicist Edmond Becquerel discovered the photoelectric effect which means that illuminated electrolyte can produce electricity [2]. In 1873, Willoughby Smith discovered photo-conductive properties of the selenium. In 1883 Charles Fritts, American inventor, relying on the discover of Willoughby Smith, creates the first solar cell made from selenium, so the origin date of photovoltaics could be considered 1883 year. Theoretical and experimental works on photoelectric effect has been done by Heinrich Hertz [3], by Aleksandr Stoletov [4], who linked photoelectric current with intensity of incident radiation. In 1905, Albert Einstein has described the photoelectric effect [5] and in 1921 was awarded by Nobel Prize for this theory. Einstein described the mechanism of the photoelectric effect which was the basis for the solar energy.

Nowadays three generations of solar cells batteries exist. The first generation of solar cells with efficiency about 6% was produced from crystalline silicon [6]. The next generation tries to reduce the manufacturing costs of the batteries and the amount of toxic waste dur-1

Part I

Background Chapter 1 The Many-Body problem

The description and the understanding of the electronic properties of a material requires to deal with an system of a huge number of electrons, which are interacting: this is called the Many-Body problem. This problem is completely contained in the Schrödinger equation, which takes into account all the possible interactions that electrons experiment in a compound.

The Schrödinger equation

Stationary Schrödinger equation have a form presented below:

ĤΨ = EΨ (1.1)
The Hamiltonian Ĥ in the equation plays the key role since it contains all the information about energy in a quantum system. For the system of electrons and nuclei the Hamiltonian will be:

Ĥ = -2m e i ∇ 2 i - I 2 2M I ∇ 2 I - i,I Z I e 2 |r i -R I | + 1 2 i =j e 2 |r i -r j | + 1 2 I =J Z I Z J e 2 |R I -R J | (1.2)
where m e , e are the mass and charge of the electron, is a Planck constant 1 and M, Z are the mass and charge of the nucleus. Indexes i, j correspond to the electrons and I, J correspond to nucleus. The first and second terms in the equation 1.2 are the kinetic energy of the electrons and the kinetic energy of the nuclei, the three other terms are Coulomb interactions. The third one is the electron-nuclei interaction, the fourth term is the electronelectron interaction and the last term is the nuclei-nuclei interaction.

The complexity of the Hamiltonian rapidly increases with the number of interacting particles and the direct solution of the Schrödinger equation with the Hamiltonian (1.2) is not practically possible for a real systems with lot of interacting particles. Even for a small piece of matter about 1 gram the number of particles is about 10 23 which makes impossible the brute force solving of the Schrödinger equation. The only way out to treat this equation is to make approximations.

Born-Oppenheimer approximation

The first thing we could notice is that the mass of the nuclei is several orders of magnitude larger then the mass of the electrons. The mass of the proton is approximately 1836 times more then the mass of electron. The nuclei moves more slowly and kinetic energy of the nuclei is much smaller then the kinetic energy of the electrons. In this context we can perform the decoupling of the Schrödinger equation into a set of two equations, one will describe the motion of the electrons and another one will encapsulate the motion of the nuclei. This approximation is called Born-Oppenheimer [9] or adiabatic approximation. We make an assumption that the initial wavefunction ψ(r.R) can be rewritten as a product of two wavefunctions which are solutions of the equation of the electrons and equation of nuclei Ψ R (r) and Φ(R) respectively:

ψ(r, R) = Ψ R (r)Φ(R)
The Schrödinger equation for the electrons will have following form

[T e + V ee + V Ie ] Ψ R (r) = E R Ψ R (r) (1.3)
where the T e is a kinetic energy of the electrons, V ee is a electron-electron interaction and V Ie is a nuclei-electron interaction. The index R tell us that the position of the nuclei is fixed and R is the parameter, while the position of the electrons r is a variable. In these equations r stands for (r 1 , r 2 , ...) and R for (R 1 , R 2 , ...).

In the thesis, we will consider the nuclei motion to be very slow compared to the motion of the electrons and we will take into account only the electronic equation 1.3, R being kept fixed at the equilibrium position. However, if one needs to take into account the motion of the nuclei, for example in case of phonons, the nuclear equation will have to be considered.

Hartree and Hartree-Fock approximations

The main problem of solving the Schrödinger equation in a Born-Oppenheimer approximation is that the wavefunction which describes N electrons ψ(r 1 , r 2 , ..., r N ) is a N -variable function. If one tries to solve the problem numerically on a grid and takes, for instance, 10-point grid, for two electrons the only x-coordinate requires a table 10x10, for 3 electrons this is the 10x10x10 cube of information. It is clear that to model one single atom of carbon in three dimensions which is an important component of all the living beings, we need a 9-dimension table with 10 9 elements. And the quality of the solution will be extremely poor.

The first attempt in order to solve the equation was suggested by Hartree [START_REF] Hartree | Self-consistent field, with exchange for beryllium[END_REF]. It is based on the assumption that the wavefunction can be written as a product of one-particle 1.2. HARTREE AND HARTREE-FOCK APPROXIMATIONS wavefuntions called orbitals: ψ(r 1 , r 2 , ..., r N ) = ϕ 1 (r 1 )ϕ 2 (r 2 )...ϕ N (r N ) If the modulus of the wavefunction has a physical meaning, the orbitals doesn't. This approximation has the flaw that the product does not take into account the anti-symmetry principle, coming from the fact that electrons are fermions. This principle says that we don't know which of the electrons we describe with the given orbital, so the orbitals can exchange. And after every exchange the wavefunction has to change the sign. The Hartree approximation does not satisfied the anti-symmetric condition.

The next development has been suggested by Fock [11]. The wavefunction is now approximated by a Slater determinant, constructed with the one-particle orbitals. The general form for a system of N electrons will be 1.4:

ψ HF (ξ 1 , ξ 2 ...ξ N ) = 1 √ N ! ϕ 1 (ξ 1 ) ϕ 2 (ξ 1 ) ... ϕ N (ξ 1 ) ϕ 1 (ξ 2 ) ϕ 2 (ξ 2 ) ... ϕ N (ξ 2 ) ... ... . . . ... ϕ 1 (ξ N ) ϕ 2 (ξ N ) ... ϕ N (ξ N ) (1.4)
with ξ i = (r i , σ i ), r i being the position and σ i the spin.

The electronic density for the Hartree-Fock approximation is expressed as:

ρ(ξ) = N i=1 |ϕ i (ξ)| 2 (1.5)
Spin-orbitals can be found by minimizing the total energy of the system. Thus, we vary the expectation value of the Hamiltonian (1.2) with respect to the single-particle orbitals:

δ δϕ * α (ξ) Φ| Ĥ|Φ - N i=1 i dξ ϕ i (ξ )ϕ * i (ξ ) = 0 (1.6)
After mathematical manipulations, one gets

- 2 2m ∇ 2 + V Ie (ξ) + dξ ρ(ξ )v(ξ, ξ ) ϕ α (ξ) -ρ(ξ, ξ )v(ξ, ξ )ϕ α (ξ )dξ = α ϕ α (ξ) (1.7)
where the electronic density ρ(ξ ) is given by Eq. (1.5) and ρ(ξ, ξ ) can be expressed in terms of the orbitals:

ρ(ξ, ξ ) = N i=1 ϕ * i (ξ )ϕ i (ξ) (1.8)
This is one of the fundamental approximations for solving the Schrödinger equation.

Nevertheless, the Hartree-Fock approximation does not give an exact solution to the Many-Body problem.

Ab initio formalisms aim to give to this problem a in principle exact description. (The approximations will appear later, in the practical resolution since some quantities involved in the formalism are unknown, as it will be see in the following.)

Two ab initio frameworks can be distinguished. One of them proposes to deal with independent electrons evolving in a effective potential mimicking the interactions of all the surrounding electrons. It is the Density Functional Theory (DFT). Since DFT [and its extension to describe electronic excitations Time-Dependent Density Functional Theory (TDDFT)] are the formalism mainly used in this thesis, they will be described in details in further chapters.

The second ab initio framework proposed to solve the fully interacting electrons system is called Many-Body Perturbation Theory (MBPT) [12][START_REF] Mahan | Many Particle Physics[END_REF][START_REF] Schwinger | [END_REF]. Since it has been used in the thesis without developement, it will be briefly summarized here.

Many-Body Perturbation Theory 1.3.1 The GW approximation

The main idea is the concept of quasi-particles [15]. Since electrons are not independent, once the electron has left its orbital, there is a hole in the place of this electron and since there is a positive charge, electrons are attracted and tend to screen this positive coulomb potential (Fig. 1.1, right). It is called quasi-hole. Similarly, the electron which has moved repels the other electrons and feels the positive remaining charge, then the electron appears as surrounded with a positive charge cloud: it is called a quasi-electron. The quasi-electrons interact with a reduced Coulomb interaction, and for this reason it can be treated within perturbation theory.

Ground state

Independent particle

Quasiparticle ħν ħν The differences between the quasiparticles and the bare ones are contained in a nonlocal, non-hermitian and energy dependent operator Σ, called self-energy, which contains all 1.3. MANY-BODY PERTURBATION THEORY the many-body exchange and correlation effects. In this context of Many-Body Perturbation Theory, the quasi-particles energies and wavefunctions are determined by solving the Schrödinger-like equation:

- 1 2 ∇ 2 + V ext (r) + V H (r) ψ(r) QP + dr Σ(r, r ; E QP k )ψ QP (r ) = E k ψ QP (r) (1.9)
A standard and very efficient framework used to calculate the self-energy is the so-called GW approximation. G is the one particle Green's function. It describes the propagation of a particle from a point r 1 at time t 1 , to a point r 2 at time t 2 . W is the screened Coulomb interaction and it is calculated with W = ε -1 v, where ε -1 is the inverse dielectric matrix and v is the bare Coulomb interaction. Many calculations using this approximation were performed successively and made this method the state-of-the-art for the evaluation of the band structures [16,17].

Compare to LDA-approximation GW-method give relatively good result. The image 1.1 is taken from the paper [18]:

Bethe-Salpeter equation

Among the energies and wavefunctions of the quasi-particles, one is interested to calculate absorption spectra. In this kind of spectroscopy, the excited electron (quasi-electron) and the hole (quasi-hole) are simultaneously present in the material, and they can interact, giving rise of an exciton. In order to describe such a phenomenon, one needs 2-particles Green's function G [START_REF] Berna | Microstratigraphic evidence of in situ fire in the Acheulean strata of Wonderwerk Cave, Northern Cape province, South Africa[END_REF]2,3,4), which describes the propagation of an electron from a point r 1 at time t 1 , to a point r 3 at time t 3 , and the propagation of the hole a point r 2 at time t 2 , to a point r 4 at time t 4 .

The 4-point polarizability L(1, 2, 3, 4) is defined according to L(1, 2, 3, 4) = iL 0 (1, 2, 3, 4) -G(1, 2, 3, 4) (1.10) where L 0 (1, 2, 3, 4) = -iG(1, 3)G (2,4) describes the propagation of the electron and the hole separately. This function satisfies a Dyson equation, known as the Bethe-Salpeter equation [19]:

L(1, 2, 3, 4) = L 0 (1, 2, 3, 4) + d(5678) L 0 (1, 2, 5, 6) K (5,6,7,8) L (7,8,3,4) with K(5, 6, 7, 8) = δ(5, 6) δ(7, 8) v(5, 7) + i δΣ (5, 3 δG(8, 4) (

To obtain the absorption spectrum in such a framework, one solves the modified Dyson equation:

L(1, 2, 3, 4) = L 0 (1, 2, 3, 4) + d(5678) L 0 (1, 2, 5, 6) K(5, 6, 7, 8) L (7,8,3,4) with K(5, 6, 7, 8) = δ (5,6) δ (7,8) v(5, 7) + i δΣ (5, 3 δG(8, 4) (1.12)

where v is the Coulomb potential without the long range component.

The absorption spectrum will be given by Im[ε M (ω)] ε M (ω) = 1 -lim q→0 v(q) drdr e iq(r-r ) L(r, r, r , r , ω) (

.

Chapter 2

Density Functional Theory

The starting point for another approach of solving the Schrödinger equation is the electronic density. It is called density functional theory. The basic idea of the theory is that any property of the many-body system can be found through a functional of the density n 0 (r). This method became one of the general methods for calculating of electronic structures, [START_REF] Martin | Interacting Electrons. Theory and Computational approaches[END_REF].

Hohenberg-Kohn theorem

Density functional theory is based on the two theorems which are presented below.

Theorem I: For any system of interacting particles in an external potential V ext (r), the potential V ext (r) is determined uniquely, except for a constant, by the ground state particle density n 0 (r). It means that there is unique correspondence between the ground state density of many-electron system and the external potential [START_REF] Hohenberg | Inhomogeneous Electron Gas[END_REF].

As a consequence, it could be considered that ground-state properties of the manyelectron system can be written as a functional of the density only. If we rewrite any property as an expectation value, we will get

φ 0 | Ô | φ 0 = O [n]
where φ 0 is a many-body ground state wavefunction.

Theorem II: A universal functional for the energy E[n] in terms of the density n(r) can be defined, valid for any external potential V ext (r). For any particular V ext (r), the exact ground state energy of the system is the global minimum value of this functional , and the density n(r) that minimizes the functional is the exact ground state density n 0 (r). [START_REF] Hohenberg | Inhomogeneous Electron Gas[END_REF]

Kohn-Sham auxiliary system

The Thomas-Fermi-Dirac functional [START_REF] Thomas | The calculation of the atomic fields[END_REF][START_REF] Fermi | Un metodo statistico per la determinazione de alcune priorieta dell'atome[END_REF][START_REF] Dirac | Note on exchange phenomena in the Thomas-Fermi atom[END_REF] was let's say a pre-step in the way to DFT, then the theory has been developed by Kohn and Sham in 1965 [25]. They replaced the many-11 body system of fully interacting particles with a fictitious auxiliary system of non-interacting particles, having the same electronic density as the fully interacting one.

The auxiliary system is described by an independent particle Hamiltonian of the form 1 :

Ĥaux (r i )ψ i (r i ) = E i ψ i (r i ) Ĥaux (r i ) = - 1 2 ∇ 2 i + V KS (r i ) (2.1)
where i refers to the i-th electron.The effective potential V KS is given by

V KS (r) = V ext (r) + V Hartree (r) + V xc (r), (2.2)
V ext being the potential due to other nuclei and any others external fields and V Hartree being the Classical Coulomb interaction defined as

V Hartree [n](r) = d 3 r n(r ) |r -r | (2.3)
Here V xc is the exchange-correlation potential, which encapsulates all the many-body effects.

It is an unknown quantity which has to be approximated. The schematic image of the Khon-Sham auxiliary system with one electron which interacts with some mean potential is presented on the figure 2 A key to solution of a Kohn-Sham equation is to find the best approximation for the V xc . Further some well-known approximations of V xc will be shown.

The way of solving the Kohn-Sham equation is an iteration process. Since the set of KS equation is self-consistent, first we guess the electronic density, then substitute into the V KS and solve the KS equation (2.1). The new electronic density relying on Eq. (1.5) is calculated and, if the difference between both densities is significant compared with our level of accuracy, we substitute the new density into the V KS and repeat the procedure again. 1 We will consider atomic units, i.e. m e = e = = 1.

Functionals

Khon-Sham equation (2.1) contains one term which takes into account all the difficulties of quantum nature of the electrons. This is the exchange and correlation potential V xc . All the difficulties are included there and this quantity has to be approximated. In this section, some approximations will be shown.

LDA

The simplest approximation is the local density approximation (LDA). The exchange-correlation energy has the general form

E xc [n] = d 3 r n(r) xc ([n], r) (2.4)
where xc ([n], r) is the energy per electron at the point r. The idea of this approximation is to replace the functional of the density xc ([n], r) by the exchange-correlation energy per electron of the homogeneous electron gas (HEG) (for a detailed review see [START_REF] Martin | Electronic structure: basic theory and practical methods[END_REF]):

E LDA xc [n] = d 3 r n(r) hom xc (n(r)) (2.5)
where xc (n(r)) is a function (not a functional) of the density n(r) at point r. In LDA, the exchange-correlation energy is typically separated into the exchange part and the correlation part: xc = x + c . The exchange part takes the form x = n 1/3 . For homogeneous interacting electron gas, reliable parametrizations of c are available from quantum Monte-Carlo simulations [START_REF] Ceperley | [END_REF]28].

The exchange-correlation potential is given in terms of the energy as

V xc (r) = xc ([n], r) + n(r) δ xc ([n], r) δn(r) (2.6)
The LDA exchange term is known and reads

V LDA xc (r) = 4 3 hom xc (n(r)).
Even if this method is very simple, it can give accurate results for a wide range of insulators.

The LDA has been further extended to spin resolved local density approximation:

xc (n) -→ xc (n ↑ , n ↓ ),
where n ↑ and n ↓ are the up and down spin densities respectively. Nevertheless, the LDA has a tendency to over-estimate the exchange-correlation energy [29]. To correct for this tendency, a new class of functionals has been developed, where E xc becomes a functional of the gradient of the density:

xc (n ↑ , n ↓ ) -→ xc (n ↑ , n ↓ , ∇n ↑ , ∇n ↓ ).
This allows for corrections based on the changes in density away from the coordinate. These expansions are referred to as generalized gradient approximations (GGA) [30][31][32].

One example of such a parameter-free GGA functional is the one developped by Perdew, Burke and Ernzerhof (PBE) [33]. It is known for its general applicability and gives rather accurate results for a wide range of systems.

Hybrid functional

Hybrid functionals are a class of approximations to the exchange-correlation energy functional in DFT that incorporate a portion of exact exchange from Hartree-Fock theory with exchange and correlation from other sources (ab initio or empirical).

For example, the hybrid functional called P BE0 (Perdew -Burke -Ernzerhof) [START_REF] Perdew | Rationale for mixing exact exchange with density functional approximations[END_REF][START_REF] Adamo | Toward reliable density functional methods without adjustable parameters: The PBE0 model[END_REF] has the form:

E P BE0 xc = αE HF x + (1 -α)E P BE x + E P BE c (2.7)
where E P BE x and E P BE c corresponds to the exchange and correlation parts of the PBE GGA functional [30][31][32]. α is the mixing factor.

The Heyd-Scuseria-Ernzerhof (HSE) exchange-correlation functional [36] keeps the correlation part of the P BE0 functional, but applies a screened Coulomb potential to the exchange interaction in order to screen the long-range part of the Hartree-Fock exchange. The Coulomb interaction is cut into a short-range (SR) and a long-range (LR) contributions, with an error function depending of the parameter ω which governs the extent of the short-range interactions. It leads to:

E P BE0 x = αE HF,SR x (ω) + αE HF,LR x (ω) + (1 -α)E P BE,SR x (ω) + E P BE,LR x (ω) -αE P BE,LR x (ω)
Numerical tests with realistic values of ω shows that E HF,LR

x and E P BE,LR x cancel and it results a screened Coulomb potential hybrid density functional of the form:

E wP BEh xc = aE HF,SR x (ω) + (1 -a)E P BE,SR x (ω) + E P BE,LR x (ω) + E P BE c (2.8)

Chapter 3

Time-Dependent DFT

The calculation of absorption spectra cannot be obtained in the framework of static DFT, as it requires the knowledge of excited states. To access these quantities, one can extend the DFT to Time-dependent Density Functional Theory (TDDFT), as introduced by Runge and Gross [37].

Runge-Gross theorem

Following the procedure used for DFT, where the full many-body problem was transformed into a non-interacting particle system, a time-dependent external perturbation is now added to the static potential. The electronic density of the system is denoted by n(r, t) and we consider that the external perturbation is switched-on adiabatically. For such a system, Runge and Gross established the following theorem [37].

Runge-Gross theorem: Two densities n(r, t) and n (r, t), evolving from a common initial many-body state Ψ 0 under the influence of two different potentials v(r, t) and v (r, t) = v(r, t)+c(t) (both assumed to be Taylor-expandable around t 0 ), will start to become different infinitesimally later than t 0 . Therefore, there is a one-to-one correspondence between the time-dependent density and the time-dependent potential, for any fixed initial many-body state. This theorem plays the same role as the Hohenberg-Kohn theorem of DFT.

Relying on this theorem, Runge and Gross suggested the use of an auxiliary system of non-interacting particles, leading to a time-dependent version of the Kohn-Sham equations for the time-dependent electronic density.

Time-dependent Kohn-Sham equations

The exact electronic density can be obtained from the single-particle orbitals, satisfying the time-dependent Kohn-Sham equation.

Under the influence of time-dependent potential V ext (r, t) the time-dependent electronic density becomes:

n(r, t) = N j=1 |φ j (r, t)| 2 (3.1)
where the single-particle orbitals φ(r, t) can be found by solving :

- ∇ 2 2 + V KS ([n], r, t) φ j (r, t) = i ∂ ∂t φ j (r, t) (3.2)
together with the initial conditions

φ j (r, t 0 ) = φ 0 j (r) (3.3)
The time-dependent Khon-Sham potential becomes:

V KS ([n], r, t) = V ext (r, t) + d 3 r n(r , t) |r -r | + V xc ([n], r, t) (3.4)

Perturbation theory

If we perturb the system, the system will react. The quantity which connects the external perturbation and the reaction of the system is called susceptibility and is represented by χ.

If the perturbation is weak, the response of the system can be expanded in power of the perturbation and it is convenient to split the response in two types: linear and nonlinear.

In the linear regime, the response of the system to the perturbation depends linearly on the perturbation. In the framework of TDDFT, we describe the response of the system to the perturbation in terms of the induced electronic density n ind (see Figure 3.1). The relation between n ind and the external potential V ext is given by

✓

n ind (r, t) = χ(r, t; r , t ) V ext (r , t ) dr dt (3.5)
where χ is the density response function or fully-interacting response function, defined also by

χ(r, t; r , t ) = δn(r, t) δV ext (r , t ) (3.6)
The density is "measured at time" t and position r, while the potential has been applied at time t and position r .

In the auxiliary non-interacting system, the induced density is related to the Kohn-Sham potential V KS by

n ind (r, t) = χ 0 (r, t; r t ) V KS (r , t ) dr dt (3.7)
where V KS is given in Eq. 3.4. Using the "chain rule" δF δψ(x)

= dy δF δG(y) δG(y) δψ(x) (3.8) and the functional relation n(x) n(y) = δ(x -y), we have

χ(r, t; r , t ) = dr 2 dt 2 δn(r, t) δV KS (r 2 , t 2 ) χ 0 δV KS (r 2 , t 2 ) δV ext (r , t ) = = dr 2 dt 2 χ 0 (r, t; r 2 , t 2 )   δV ext (r 2 , t 2 ) δV ext (r , t ) + δV H (r 2 , t 2 ) δV ext (r , t ) + δV xc (r 2 , t 2 ) δV ext (r , t )   (3.9)
The first derivative is

δV ext (r 2 , t 2 ) δV ext (r , t ) = δ(r 2 -r ) δ(t 2 -t ) (3.10)
We have for the Let's combine all the results of three terms together. And, in shorter form:

dr 2 dt 2 χ 0 (r, t; r 2 , t 2 ) δ(r 2 -r ) δ(t 2 -t ) + dr 3 dt 3 v(r 2 , r 3 ) χ(r
χ = χ 0 + χ 0 [v + f xc ]χ (3.15)
The exchange-correlation kernel f xc takes into account all the difficulties of exchange and correlation in the dynamical linear regime. The simplest approximation for f xc is a Random Phase Approximation (RPA) when f xc = 0, but more advanced approximations are available, for instance based on the functional derivative of V LDA xc with respect to the electronic density.

The Coulomb potential v is responsible for the so-called local field effects (LF). They account for the fluctuations of the induced electronic density at the microscopic scale, corresponding to the atomic structure of the matter.

Due to the time invariance, the time Fourier transform gives a one frequency dependence. Assuming a periodic crystal and performing space Fourier transform, the equation (3.15) can be rewritten in the reciprocal space. Assume that the reciprocal vector k is equal k = G + q where q is in the first Brillouin zone, the space integrations becomes products of matrices, and we finally get:

χ G,G (q, ω) = χ 0 G,G (q, ω) + G 2 χ 0 G,G 2 (q, ω)v G 2 (q)χ G 2 ,G (q, ω) (3.16)
where the χ 0 G,G (q, ω) is calculated by summing the transitions between Khon-Sham states {|ϕ j k , ε j k }, calculated according to

χ 0 G,G (q, ω) = 1 N k Ω v,c,k (f v k -f c k+q ) ε v k -ε c k+q + ω + iη ϕ v k |e -i(G+q)•r |ϕ c k+q ϕ c k+q |e i(G +q)•r |ϕ v k (3.17)
Here k refers to the k-point grid of the first (irreducible) Brillouin zone. N k is the number of these points. v, c are the states between which a transition is calculated. Ω is the volume of the unit cell in real space.

These quantities describe microscopic phenomena. The term v 0 correspond to the longrange part of the Coulomb potential, and the other components v G for G = 0 allows the inclusion of the local field effects.

The knowledge of χ gives access to the absorption spectra with the following relation:

ε -1 (r, r ) = δ(r -r ) + V dr v(r -r )χ(r , r ) (3.18)
what has in reciprocal space the following form:

ε -1 G,G (q, ω) = δ G,G + v G (q)χ G,G (q, ω) (3.19)
The microscopic inverse dielectric function gives access to the absorption spectrum defined as the imaginary part of the macroscopic dielectric function: Im[ M (ω)]. Within TDDFT, M (ω) is calculated according to [38,39]:

ε M (ω) = lim q-→0 1 ε -1 00 (q) (3.20)
Chapter 4

Supercell formalism

In solid state physics, crystals are considered to be infinite in the three directions of space and periodic. This periodicity defines a Bravais lattice and a basis of primitive vectors (a 1 , a 2 , a 3 ). At each node of the Bravais lattice is attached a motif, which correspond to the different atoms composing the crystal. The repeating pattern is called the unit cell of the structure. Due to this periodic description, a natural framework is to use the reciprocal space, and a planewave basis set. Nevertheless, such an object does not take into account neither the presence of the surface, nor the finite thickness of the sample. To model a surface or a slab in a 3D periodic formalism, we have to construct a supercell, composed of the piece of matter we want to describe and vacuum along the z-direction perpendicular to the surface of the sample.

I will first give a brief description of the periodic boundary conditions, illustrate the difficulties of modeling the surfaces in this context and summarize the so-called "Selected-G" method developed previously [START_REF] Tancogne-Dejean | Ab initio description of second-harmonic generation from crystal surfaces[END_REF].

Periodic structure of crystal, Bloch's Theorem and planewave basis set

Since the crystal is described as an infinite repetition of the unit cell, the ions build a periodic lattice and the interaction between electrons and ions is represented by a periodic potential

V (r + R n ) = V (r), where R n = n 1 a 1 + n 2 a 2 + n 3 a 3 .
According to the Bloch's theorem, the eigenstates for the single-particle Hamiltonian with this periodic potential can be written as

ψ n,k (r) = e ik.r u n,k (r) (4.1)
where n label the electronic level of energy, k is a wavevector and u n,k (r) is a function which has the periodicity of the crystal lattice:

u n,k (r + R) = u n,k (r) (4.2)
For this reason, it can be expended on a basis of planewaves, whose wavevectors are the reciprocal lattice vectors {G}: {e iG.r }

G = n 1 b 1 + n 2 b 2 + n 3 b 3 (n i ∈ Z) (4.3)
where (b 1 , b 2 , b 3 ) are the basis vectors of the reciprocal space (b i .a j = δ ij ).

u n,k (r) = G c n,k (G)e iG.r (4.4)
In the following, the wavevector k will be written as a sum of a vector of the reciprocal space and a vector inside the first Brillouin zone (q): k = G + q,

Modelisation of an isolated object or a surface

To model an isolated object in a framework inherited from the 3D-periodicity, in order to take advantage of the efficiency of reciprocal space and the planewaves-based formalism, we design a supercell, composed of the object itself surrounded by vacuum. We will consider the example of a slab, which is infinite and periodic in the (x, y) plane, and finite in the z-direction, with thickness L mat z . To isolate the slab, one use a supercell as defined on fig 4.1. In that case, the supercell is usually defined according to (a 1 , a 2 , na 3 ), where n is the ratio between the thickness of matter and the height of the supercell.

The effect of using a supercell approach has been studied in Ref [START_REF] Tancogne-Dejean | Ab initio description of second-harmonic generation from crystal surfaces[END_REF][START_REF] Tancogne-Dejean | [END_REF] and I will show some of the results obtained for the silicon (2 x 1) surface. The surface is characterized by asymmetric dimers formed between two atoms at the top and the bottom of the structure. The DFT calculation within local density approximation has been performed to get groundstate electronic properties. This configuration is shown on the figure 4.2.

The absorption spectra, calculated in the Independent-Particle Approximation (IPA) and in the Random Phase Approximation (RPA), that is including the local field effects, are shown for three different amounts of vacuum depicted in figure (4.3) and called void1 (half matter/half vacuum: L SC = 2L mat z ), void2 (2 times more thickness of vacuum than matter: L SC = 3L mat z ) and void3 (3 times more thickness of vacuum than matter: L SC = 4L mat z ). The imaginary part of M for the in-plane and out-of-plane components, using IPA, is presented in Fig. (4.3). One clearly sees that the amplitude of the spectra decreases when increasing the vacuum. This decrease is actually only a scaling factor: when one scales each curve by the factor L SC L mat z , all the curves are superimposed to the black curve labeled "renorm". This effect comes from the fact that in such a calculation, the IPA response function is given by the χ 0 00 , which is proportional to 1 O mega ∝ 1 L SC . The imaginary part of M for the out-of-plane component, including local field effects, is presented in figure (4.4). The inset in this figure shows the spectrum for void1 compared to the bulk one. One sees that the absorption peak for the supercell void1 is located around 12 eV, and the amplitude is strongly reduced. When comparing the spectra for the different thicknesses, the amplitude is not only changing, but the peak position is shifting towards high energy. This result is not physical. It has been shown that the regular TDDFT in the supercell formalism is equivalent to Effective Medium Theory (EMT) with vacuum [START_REF] Tancogne-Dejean | [END_REF]. In this reference, we have distinguished the two limiting cases of minimum screening and maximum screening. We have shown that the minimum screening case was also associated with small local field effects, and was suitable to describe the spectra for in-plane components, leading to a simple scaling factor. It reads:

ε SC M = f + (1 -f ) ε bulk M (4.5)
The maximum screening corresponds to large local field effects and allows us to interpret the out-of-plane component; it can be written as:

1 ε SC M = f + (1 -f ) ε bulk M (4.6)
The label SC means supercell and corresponds to the spectra calculated with the regular TDDFT procedure. f is the proportion of vacuum in the supercell.

For the out-of-plane component, the application of Eq. (4.6), using the macroscopic dielectric function of the bulk system for ε bulk M gives the spectra in dashed lines in Fig. 4.5, which account very well for the behavior of the standard TDDFT spectra (solid lines) with vacuum.

The selected-G method

In the previous section, it was shown that the calculation using a supercell with vacuum leads to nonphysical results, and this spurious effect must be corrected.

It has been done during the PhD-thesis of N. Tancogne-Dejean [START_REF] Tancogne-Dejean | Ab initio description of second-harmonic generation from crystal surfaces[END_REF], where the method called "Selected-G" has been developed. It allows one to solve the problem of the vacuum size dependency for the description of a surface. In the framework of the Selected-G method, we define auxiliary response functions (Independent-particle response functions and susceptibilities). These functions are periodic with the period L mat z along the direction perpendicular to the slab and they have the same values as the corresponding quantities for the isolated slab in the region inside the matter [-L mat z , 0]. Therefore, one has a new set of G-vectors, associated to the periodicity of the auxiliary functions, Gz = 2πn L mat z , which are different from the one having the periodicity of the supercell:

G z = 2πn
Lz . The vectors of the reciprocal lattice are defined according to: G = ( Ḡ, Gz ), where Ḡ stands for the in-plane component.

Since the slab keeps its periodicity for the in-plane directions, we consider that the inplane Fourier transform (FT) has already been done, and we only consider the z, z remaining FT. In such a framework, the FT of the auxiliary response functions has the following form:

χ G, G (q; ω) = 1 L mat z 0 -L mat z dz 0 -L mat z dz e -i(qz+ Gz)z × χ Ḡ, Ḡ (q, z, z ; ω)e i(qz+ G z )z χ Ḡ, Ḡ (q; ω) = 1 L mat z qz Gz, G z e i(qz+ Gz)z × χ G, G (q; ω)e -i(qz+ G z )z (4.7)
These auxiliary response functions χ satisfy a modified Dyson equation:

χ G G (q; ω) = χ 0 G G (q; ω) + G1 G2 χ 0 G G1 (q; ω) Ṽ G1 G2 (q) χ G2 G (q; ω) (4.8)
where Ṽ G1 G2 (q) is defined according to

Ṽ G1 , G2 (q) = 1 L mat z 0 -L mat z dz 1 0 -L mat z dz 2 e -i(qz+ Gz1 )z 1 ×v Ḡ1 (q; z 1 , z 2 ) e -i(qz+ Gz2 )z 2 δ Ḡ1 , Ḡ2 (4.9) with v Ḡ1 (q; z 1 z 2 ) = 1 |q + Ḡ1 | e -|q + Ḡ1 ||z 1 -z 2 | is the 2D Fourier transform of the Coulomb potential. The full expression of Ṽ G1 G2 (q) is [41]. Ṽ G1 , G2 (q) = 4π |q + G1 | 2 δ G1 , G2 + ξ 4π δ Ḡ1 , Ḡ2 |q + G1 | 2 |q + G2 | 2   - e -|q + Ḡ1 |L mat z sin(q z L mat z ) L mat z (2q z + Gz1 + Gz2 ) + e -|q + Ḡ1 |L mat z cos(q z L mat z ) -1 L mat z |q + Ḡ1 | [|q + Ḡ1 | 2 -(q z + Gz1 )(q z + Gz2 )   (4.10)
This potential is be called slab potential. The first term is the standard expression for the 3D-FT of the Coulomb potential, and the second term is a correction due to the finite thickness of the slab, proportional to 1/L mat z and depending on two Gz reciprocal lattice vectors. The second term contains a phase factor ξ, which will be explain latter.

The corresponding microscopic dielectric matrix is

ε -1 G G (q; ω) = δ G G + G1 Ṽ G G1 (q) χ G1 G (q; ω) (4.11)
To describe a surface, one considers the limit L mat z -→ ∞. In that case, the slab potential is reduced to the standard diagonal 3D-Fourier transform of the Coulomb potential. The absorption spectrum is deduced from ε -1 00 (q; ω) with the usual relation : The Selected-G formalism has been successfully applied to the case of surfaces [START_REF] Tancogne-Dejean | [END_REF]42]. The absorption spectra of the (2x1) silicon surface, including local field effects, are presented in figure 4.6. Bulk silicon (black solid line) is presented for reference. The in-plane components xx , yy and the outof-plane component ⊥ are shown with blue, green and red colors respectively.

ε M (ω) = lim q-→0 1 ε -1 00 (q) (4.12)
We can see that the in-plane components have spectra comparable with the bulk, both for the amplitude and energy position. The difference between the in-plane components of the (2x1) surface as compared to the bulk, comes from the presence of surface states in the (2x1) silicon surface, due to the creation of the dimers after the reconstruction of the surface. These states are responsible for the structures at low energy (below 2.5 eV), and explain the anisotropy between the xx and yy in-plane components.

The out-of-plane component ⊥ (dashed red) is now located in an energy range closed to the one of the bulk, so the unphysical shift towards high energy for the out-of-plane component has disappeared. We have also demonstrated that this spectrum is independent of vacuum, and move towards the bulk spectrum when increasing the thickness of the slab [START_REF] Tancogne-Dejean | Ab initio description of second-harmonic generation from crystal surfaces[END_REF].

Part II

Developments and Applications

Chapter 5

Application for finite slab systems

The slab potential

The slab potential has been defined to calculate the dielectric properties of an isolated slab, in a framework of 3D code, based on planewaves and periodic boundary conditions. The standard procedure is to define a supercell, containing the matter and some vacuum, to prevent the slab to interact with artificial replicas through the long-range Coulomb potential. But, we evidenced that the thickness of vacuum strongly influences the absorption spectrum for the excitation perpendicular to the slab, when accounting for the local field effects.

Based on the idea that the response functions are localised in the matter, we considered auxiliary response functions equal to the response functions of the isolated slab, but having the periodicity of the matter L mat z . Then we calculated the Fourier transform of the Dyson equation on the volume corresponding to the matter.

As already stated in the previous chapter, the periodicity of the auxiliary response functions has the consequence that the reciprocal space vectors are multiple of 2π L mat z . They correspond to a subset of the reciprocal space vectors built with the height of the supercell, when the ratio between supercell and matter is an integer. It justifies the name "selected-G" vectors. The Dyson equation is solved only on this reduced basis and has to be modified accordingly [see Eq. (4.8)]:

χ G G (q; ω) = χ 0 G G (q; ω) + G1 G2 χ 0 G G1 (q; ω) Ṽ G1 G2 (q) χ G2 G (q; ω)
where the slab potential is [see Eq. (4.10)]:

Ṽ G1 , G2 (q) = 4π |q + G1 | 2 δ G1 , G2 + ξ 4π δ Ḡ1 , Ḡ2 |q + G1 | 2 |q + G2 | 2   - e -|q + Ḡ1 |L mat z sin(q z L mat z ) L mat z (2q z + Gz1 + Gz2 ) + e -|q + Ḡ1 |L mat z cos(q z L mat z ) -1 L mat z |q + Ḡ1 | [|q + Ḡ1 | 2 -(q z + Gz1 )(q z + Gz2 )  
Ḡ stands for the in-plane component of G.

The first term is diagonal in G. It is the standard expression for the 3D-FT of the Coulomb potential.

CHAPTER 5. APPLICATION FOR FINITE SLAB SYSTEMS

The second term is a correction due to the finite thickness of the slab, proportional to 1/L mat z and depending on two G reciprocal lattice vectors. It is diagonal for the in-plane component Ḡ, which gives for the total potential the following dependence: Ṽ G1 G2 (q) = Ṽ Ḡ1 Ḡ2 ; Gz1 Gz2 (q) δ Ḡ1 , Ḡ2 . The second term also contains a phase factor ξ. It arises from the definition of the cell in the integration. When the integration is done between [-L mat z , 0] or [0, L mat z , 0], ξ = 1, while when the integration is done between [-L mat z /2, L mat z /2], ξ = (-1) n 1 +n 2 , where n i is the reduced coordinate of Gzi = 2πn i L mat z .

To describe the surface, we only used the first term, since we considered the limit L mat z -→ ∞, which has the consequence to suppress the correction.

We will now extend the study to finit slabs and perform the calculations with the full expression of the potential. We will also consider the Electron Energy Loss spectra, since some measurements on slabs of graphene exist, and will be used to validate our formalism. Before presenting the results, we recall the main quantities of interest and the way they are evaluated within this formalism.

Absorption

The absorption spectrum is the imaginary part of the macroscopic dielectric function ε M (ω). It is calculated within TDDFT according to :

ε M (ω) = lim q-→0 ε M (q; ω) (5.1)
where ε M (q; ω) is given, for non local fields (NLF) and local fields (LF) cases respectively, by:

ε N LF M (q; ω) = 1 - 4π q 2 χ 0 00 (q; ω) (5.2) ε LF M (q; ω) = 1 - 4π q 2 χ00 (q; ω) (5.3) with χ G G (q; ω) = χ 0 G G (q; ω) + G1 G2 χ 0 G G1 (q; ω) V G1 G2 (q) χ G2 G (q; ω)
and V defined as V0 G2 = 0, ∀ G2 and V G1 G2 = Ṽ G1 G2 for G1 = 0.

EELS

The EEL spectrum is the imaginary part of the inverse dielectric matrix ε -1 G G(q; ω). For momentum transferred in the first Brillouin zone, this quantity is calculated, for NLF and NLF cases, according to: ε -1,N LF 00 (q; ω) = 1 + Ṽ00 (q) χNLF 00 (q; ω) (5.4) with χNLF 00 (q; ω) = χ 0 00 (q; ω) + χ 0 00 (q; ω) Ṽ00 (q) χNLF 00 (q; ω)

ε -1,LF 00 (q; ω) = 1 + G Ṽ0 G(q) χ G0 (q; ω) (5.5) with χ G G (q; ω) = χ 0 G G (q; ω) + G1 G2
χ 0 G G1 (q; ω) Ṽ G1 G2 (q) χ G2 G (q; ω)

Remarks

One should notice that with the relations Eq. ( 5.3) and ( 5.5),

ε M = lim q-→0 1 ε -1 00 (q) (5.6)
The slab potential, due to the presence of 1/|q|, which diverges at vanishing q, prevents the direct calculation of the component of the absorption perpendicular to the surface of the slab.

Since ε M at vanishing q is a tensor, it is possible to extract this component from the linear combination of spectra calculated for in-plane and (in-plane + out-of-plane) components.

In the case of the EEL spectrum, ε -1 00 does not behaves like a tensor at vanishing q, and it is not possible to extract a q z component from any linear combination. For EELS, we will only consider in-plane spectra.

Moreover, due to the definition of the EELS (Eq. 5.5), the spectrum strongly depends on the value of vanishing q. To calculate the EELS spectra at vanishing q, we used q = 10 -2 a.u..

A detailed behavior of the slab potential is summarized in Appendix A.

The slab potential has been implemented in the DP code, as well as the above expressions, and all the calculations within TDDFT were done with this code [43]. The DFT calculations are done with Abinit [44].

Absorption for a finite slab

To check the behavior of the spectra with the vacuum, I have considered two different systems: a slab of matter of 16 The local fields are included by taken into account the microscopic components associated with G-vectors defined in a rectangular box which size is npwmat_xy × npwmat_z. For the in-plane G xy -vectors, I used 25, and for the out-of-plane G z , I used 33. It corresponds to the following closed shells:

( 0 0 ±G z ) ( 0 1 ±G z ), (0 -1 ±G z ) ( 1 0 ±G z ), (-1 0 ±G z ) ,( 0 2 ±G z ) ( 0 -2 ±G z ) ( 1 1 ±G z ) ( -1 1 ±G z ) ( 1 -1 ±G z ) (-1 -1 ±G z ), (1 2 ±G z ), (-1 2 ±G z ) ( 1 -2 ±G z ) (-1 -2 ±G z ) ( 0 3 ±G z ), ( 0 -3 ±G z ) ( 1 3 ±G z ), (-1 3 ±G z ), ( 1 -3 ±G z ), (-1 -3 ±G z ) ( 2 0 ±G z ), (-2 0 ±G z ), ( 0 4 ±G z ), ( 0 -4 ±G z )
The G z values for 4_4 supercell, due to selected-G procedure are: G z = 0, ±2, ±4, ±6, ±8, ±10, ±12, ±14, ±16, ±18, ±20, ±22, ±24, ±26, ±28, ±30, ±32

The G z values for 4_8 supercell, due to selected-G procedure are: G z = 0, ±3, ±6, ±9, ±12, ±15, ±18, ±21, ±24, ±27, ±30, ±33, ±36, ±39, ±42, ±45, ±48

Figure 5.1 shows the absorption spectra for this slab, calculated using the 2 supercells, with the Selected-G procedure and the slab potential.

The spectra are identical. This conclusion was already obtained in the case of the surface with the 3D Coulomb potential. These results confirm that the Selected-G method, in the presence of the slab potential, allows us to recover spectra independent of the size of the supercell. Moreover, we see that the absorption of the slab of silicon is very similar to the one of the surface. This is explained in appendix (A) where we did a complete analysis of the role of the slab potential.

EEL spectroscopy

We will now consider another kind of spectroscopy commonly used to study finite thickness systems. For this, I studied slabs of few layers graphene, to compare with the measurements of Eberlein et al [45]. This work aims to further validate the formalism.

In electron energy loss spectroscopy (EELS), a sample is illuminated by a beam of fast electrons (energies of about 200keV) and the transmitted electrons are analyzed in terms of their kinetic energy. The energy loss spectrum, defined as the probability to detect a scattered electron having lost a certain amount of energy, gives access to information about the target excitation (plasmon excitations, surface plasmons, core-electron excitations, coreloss fine structure and others applications [46]). Electron energy loss spectroscopy corresponds to the inelastic process where the impiging electron has an energy E 0 and a momentum K 0 , and the outgoing electron has an energy E 1 and a momentum K 1 . E = E 1 -E 0 is the energy loss and K = K 1 -K 0 is the momentum transferred to the sample. Depending on the range of energy loss E, one distinguishes core-loss spectroscopy (E > 50 eV ), where inner-shell electrons are excited, leading to elemental composition analysis, and low-loss spectroscopy (E between 0 and 50 eV ), where valence electrons are excited, giving rise to collective excitations named plasmon resonance [47]. In the following, we will be interested in low-energy excitations and plasmon resonance effect. The history of electron microscopy starts in 1932, guided by E. Ruska and B. von Borries after the invention of the short-focal-length magnetic polepiece lens [48][START_REF] Ruska | Über ein magnetisches Objekt für das Elektronenmikroskop[END_REF][START_REF] Hawkes | The Beginnings of Electron Microscopy[END_REF].

For years, the electron microscopy has been developed both by laboratories and industrial companies. In the 1980s the high resolution images have been obtained by varying the objective lens focus and following reconstruction of the exit-plane wavefunction [START_REF] Urban | [END_REF].

In ten years, aberration-corrected microscopes has become available, leading to the ability of gathering local information about individual atoms, shifts in atomic positions, other defects, stress and other local physical properties. General review of electronic microscopy can be found in [46,[START_REF] Egerton | Physical principles of electron microscopy: an introduction to TEM[END_REF].

When the sample is illuminated with a broad beam, the incident electron can be described by a single plane wave with a momentum K 0 . This is the Angular-Resolved regime. By rotating the sample, one can choose a specific K direction and study the loss spectrum for the corresponding point of the Brillouin zone. Changing the length of K, one can follow the dispersion of the spectra along the chosen directions. However, using electromagnetic lenses, one can focus the beam to study spatially resolved EELS: up to 0.2 nm. Our calculations will concern the angular-resolved EELS for valence electrons. One defines K = G + q, where G is a vector of the reciprocal lattice, and q a vector which belongs to the first Brillouin zone.

EEL(K

; ω) = -Im[ε -1 GG (q; ω)]
(5.7)

EELS for few layers of graphene.

Graphene is a material which possesses semi-metallic properties. Around the metallic point, the bands disperse linearly, corresponding to the massless electrons called Dirac electrons. This phenomenon is at the origin of their huge mobility. Graphene is of great interest for research [START_REF] Novoselov | [END_REF]54] and studies of the graphene properties are developing very intensively. In 2008, Eberlein et al [45] study plasmon spectroscopy in free standing graphene films. The graphene membranes have been prepared by using micromechanical cleavage [START_REF] Novoselov | [END_REF] of natural graphite on top of an oxidized Si wafer. After preparing the samples, the EEL spectra in different places of the sample have been measured in the range 0 -50 eV. The different places chosen correspond to different numbers of layers, namely 1, 2, 5 and more that 10 layers, which is considered to be bulk graphite.

This work is a perfect test case for our slab potential, since it will allow a direct comparison between our calculations and experimental data, and it will provide a validation of our formalism: the Selected-G method has already been successfully used to simulate second harmonic generation with surfaces, but the slab potential was not used.

Vacuum problem for EEL spectra

Before revisiting the work of Eberlein et al [45], I will first show the behavior of EEL spectra with vacuum introduced in the supercell formalism, to allow for the description of an isolated layer in a periodic framework. Figure 5.4 shows the EEL spectra for 1-layer graphene introduced in supercells with increasing vacuum (Fig. 5.3). The size of the cell was chosen to be 12.588, 18.882, 25.176, 31.470, 62.940, 125.880, and 314.700 Bohr. The calculations are done with q = 10 -2 a.u.. In Fig. 5.3, the white empty area corresponds to the amount of vacuum in the supercell and the shaded area corresponds to one graphene layer, with a constant thickness , defined by half of the distance between 2 layers in graphite (6.294 Bohr). The in-plane and out-of-plane components for graphene are presented in Fig. 5.4, including or not the local-field effects. Note that the calculations were done using the standard formalism, i.e. not using the selected-G method. The four panels show clearly that all the spectra (in-and out-of-plane, with and without LF) exhibit a strong dependence with vacuum. For the in-plane component, (left column), the spectra are shifted towards low energy when the vacuum increases. The effect of LF is almost negligible, as expected from a quite homogeneous electronic density in the plane of the layer. For the out-of-plane component (right column), the spectra without LF (top) show a reduction of the amplitude, and a much smaller shift towards low energy, as compared to in-plane component. The LF have a large effect, in agreement with the discontinuity of the electronic density perpendicular to the layer. Indeed, the shape of each spectra is modified, but the influence of the vacuum is less spectacular than for the in-plane case: it is closer to a scaling factor. We emphasize that the behavior of EEL spectra with vacuum is different from the one observed for the absorption. For absorption, the spectra without LF (independently of the in-or out-of-plane component considered) where simply scaled by the ratio between the volume of the matter and the supercell. With LF, the out-of-plane spectrum was strongly shifted; the in-plane spectra, due to the small effects of LF were not affected. For EELS, the in-plane spectra with and without LF are strongly affected: they shift in energy. Since the in-plane LF are rather small, the in-plane spectra with and without LF have the same behavior and are shifted in energy in the same way.

Analysis of the work of Eberlein et al

Coming back to the EEL measurements on samples of few layers of graphene, the results of the experiment for the in-plane component is presented in Fig 5 .5a. Two resonant peaks are clearly visible. The peak located around 5 eV correspond to the π plasmon, and the large structure above 10 eV to π + σ plasmon. The amplitude of the structures increases with increasing the number of layers, the peak position also shifts towards the high-energy region. For thickness larger that 10 layers, the authors consider that the bulk structure have been recovered.

This article [45] also provided numerical results (Fig. 5.5b) using a periodic framework, and a supercell. The numerical simulations were done for 1 layer (single in red), 2 layers (bilayer in green), 3 (trilayer in blue) layers and bulk graphite. One can see the two resonant peaks, which shift when changing the number of layers, in agreement with the experiment. The amplitude of the π plasmon (0 -10 eV) seems to saturate and decreases for graphite.

Despite the fact that the numerical results reproduced the experimental ones quite well, there is still a matter of controversy, that the authors mentioned in the paper. It was quoted that "the precise peak positions depend on the separation of planes and to compare the loss functions for graphite and multiple layers, we choose the separation between periodically repeated multilayers to be 5 times the separation in graphite". In order to compare with their results, I extracted from my calculations the spectra corresponding to the single layer, two and four layers of graphene, with a constant thickness of vacuum. It corresponds to 4 times the separation in graphite (instead of 5 in Ref [45]). The schematic view of the supercell is presented in the figure 5.6a. The white empty area corresponds to the vacuum, the shaded area corresponds to the matter, with also a color correspondence.

Results of the regular TDDFT calculations (without selected-G method and without slab potential) are presented in Fig. 5.6b. The red, green and blue schemes are 1, 2 and 4 layers respectively. I also add the bulk graphite with AB stacking. I got the same behaviour for the spectra, compared to the experimental results and numerical modeling. The peak position is located at a reasonable energy, and shifts towards low energy when decreasing the number of layers; the amplitude of the peaks also decreases. Our results have the same behaviour for the dependence in terms of the number of layers as in the experiment. However, the position of the peaks depends on the size of the vacuum, as it is evidenced on figure 5.4 for the in plane component.

Selected-G method with the slab potential for EELS

We have seen that the vacuum problem was solved for absorption spectrum by applying the Selected-G method, without the slab potential. We will now study the role of the Selected-G method with the slab potential, for the EEL spectra.

I have computed the EEL spectra for supercells containing a fix number of layers of graphene (4 layers) with different thickness of the vacuum, the relation matter-vacuum being then 1:1, 1:2, 1:3 and 1:4 (see Fig. 5.7). The regular TDDFT calculations (without Selected-G method and without slab potential) correspond to the four blue curves and evidence again the dependence with vacuum. The red curves correspond to calculations with the selected-G method only, i.e. without the slab potential (3D potential = surface). There are very close to the curve obtained for graphite (AB stacking), shown in black. The green curves have been obtained with the Selected-G method and with the slab potential. There are all superimposed, showing no dependence in terms of vacuum, and are very different from the surface or bulk spectra. The overall spectrum is in very good agreement with the measured layered spectra [45]. This points out the importance of taking into account the finite nature of the slab in the calculations. The collective oscillation of electrons at plasmon frequency is affect by the thickness of the slab, even for in-plane excitation.

EELS (q || ) Energy, [eV] R2 regular R3 regular R4 regular R5 regular R2 selectedG+V_slab R3 selectedG+V_slab R4 selectedG+V_slab R5 selectedG+V_slab R2 selectedG+V_3D R3 selectedG+V_3D R4 selectedG+V_3D R5 selectedG+V_3D
Graphite AB Before showing the results for different thicknesses, I also want to evidence the influence of the phase factor ξ = (-1)

n 1 +n 2 (if the matter is between [-L mat z /2, L mat z /2]), or ξ = 1 (if the matter is between [-L mat z , 0] or [0, L mat z ]).
= 1 ([-L mat z , 0] or [0, L mat z ] (matcenter=0: R4 solid red, R5: dashed blue) or ξ = (-1) n 1 +n 2 ([-L mat z /2, L mat z /2]
) (matcenter=1: R4 solid green, R5: dashed black), q ΓM = 0.47 Bohr -1 . Figure 5.8 shows that the spectra are strongly dependent on ξ, (red and blue curves calculated with ξ = 1 vs green and black curves calculated with ξ = (-1) n 1 +n 2 ), at least for large values of q. The effect is particularly important for the π + σ plasmon. This means that the way how the integration is treated is important. Nevertheless, the spectra are independent of vacuum whatever the way the integration is done.

The selected-G method with the slab potential (ξ = (-1) n 1 +n 2 ) has been applied to the 1,2,3,4 and 8 layers of graphene. The different systems have been modeled with different amount of vacuum in the supercell, corresponding to an integer factor of the thickness of the matter (this parameter is no more meaningful since we have evidenced that our formalism provides results independent of the vacuum). The size of the supercell and the ratio between matter and vacuum for the different systems are presented in Table 5 Figure 5.9 shows the EEL spectra calculated with the selected-G method using the slab potential. We clearly see that these spectra are in very good agreement with the experimental results. With the selected-G method, together with the slab potential, we are able to get spectra showing the correct behavior in terms of the thickness of the slab of matter, but independent of the amount of vacuum in the supercell.

Conclusion

In this chapter, I have evidenced that the EEL spectra, calculated within standard TDDFT, are strongly dependent on the vacuum introduced in the supercell. Unlike to what was observed for absorption spectra of surfaces, where the size of the vacuum affects mainly the out-of-plane component, with only a normalization factor for in-plane components, the effect is important for the in-plane components for EEL spectra. This was already noticed by the theoreticians, who imposed the vacuum size to obtain a plasmon at the measured energy.

I have shown that within the Selected-G method and using the slab potential, the EEL spectra were located at the measured energy, independently of the vacuum introduced in the supercell. This formalism has allowed me to reproduce the spectra measured by Eberlein et al for increasing stacking of graphene layers.

Moreover, I have shown that once the vacuum problem is removed due to the Selected-G procedure, the use of the slab potential gives an EEL spectrum very different from the surface (where we use the 3D expression for the Coulomb potential in reciprocal space). This result differs from the one obtained for the absorption spectrum: indeed, for the absorption, the surface and the slab calculations are similar. On the contrary, for EEL, the surface spectrum is similar to the one of the graphite (bulk), when the slab spectrum is different and in very good agreement with the the experimental results measured by Eberlein et al for slabs. This confirms that the plasmon, arising form collective excitations, is sensitive to the thickness of the slab, even for the in-plane components.

Chapter 6

Comparison with other cutoffs

Introduction

The slab potential has been defined to treat an isolated slab system in the framework of a 3D codes, based on plane waves periodic boundary conditions. The three ideas are the following:

1) the response functions of the isolated slab are localized inside the matter.

2) we define auxiliary response functions in real space, equal to the response functions of the isolated slab, and periodic with the periodicity L mat z . 3) we calculate the Fourier transform of the Dyson equation on the volume defined by the matter only.

The periodicity of the auxiliary response functions has the consequence that the reciprocal space basis vectors are multiple of 2π

L mat z , instead of 2π
Lz , the so-called "selected-G" vectors: the Dyson equation is solved only on the reduced basis of selected G-vectors.

The limitation of the integration along the direction perpendicular to the slab to the thickness of the matter acts as a cutoff procedure. Another procedure of cutoff has been proposed almost simultaneously some years ago by C. Rozzi et al [55] and S. Ismail-Beigi [56].

In these two cases: 1) the reciprocal space is built on the supercell volume :

k z = 2πnz Lz for integer value of n z ,
2) the Fourier transform for the z component of the Coulomb potential

v 2D c ( k; z, z ) = 2π k e -k|z-z | with k = Ḡ + q (6.1)
is treated as a function of one variable only, namely the distance |z -z |, as it is the case for 3D-infinite systems, leading to a diagonal expression for v 2D c ( k, k z ), ( Ḡ is an in-plane reciprocal lattice vector and q is an in-plane reciprocal space vector spanning the first Brillouin zone),

3) the integration along the direction perpendicular to the slab is limited to [-z c , +z c ], which will be further defined.

The expression proposed by C. Rozzi et al [55] is:

v2D C ( k, k z ) = 4π k 2    1 + e -kzc k z k sin(k z z c ) -cos(k z z c )    (6.2)
The value of z c = L z /2 is suggested to avoid interactions between the artificial replicas. The cases leading to a divergence are calculated explicitly :

v2D C ( k = 0, k z ) = 4π k 2 z 1 -cos(k z z c ) -k z z c sin(k z z c ) + 8 h ln   (α + √ 1 + α 2 )(1 + √ 1 + α 2 ) α   sink z z c k z (6.3)
h and α define the shape of the in-plane Born von Karman cell: h x = αh y = h (much larger than the unit cell). This last term disappears providing that the cutoff procedure is applied in the same way to the ionic potential, ensuring the charge neutrality of the system. The remaining expression is:

v2D C ( k = 0, k z ) = 4π k 2 z 1 -cos(k z z c ) -k z z c sin(k z z c ) (6.4)
The value for k ≡ 0 gives an non-diverging contribution, coming from the limit of the previous expression when k z → 0:

v2D C ( k = 0, k z = 0) = -2πz 2 c . (6.5)
The expression proposed by S. Ismail-Beigi [56] is the same as Eq. 6.2 for k = 0, but in order to avoid the divergence occurring when k → 0, z c is chosen as z c = L z /2. This leads to:

v2D C ( k, k z ) = 4π k 2    1 -(-1) nz e -kzc    (6.6)
where n z is defined with

k z = 2πnz Lz . For k z = 0, it comes: v2D C ( k, k z = 0) = 4π k2 1 -e -kzc (6.7)
and the limit for k → 0 is:

v2D C ( k -→ 0, k z = 0) = 4πz c k (6.8)
It diverges as k → 0, which gives a milder divergence as the untrucated 3D case. The choice of z c has confined the divergence to the single wave vector k = 0. This also corresponds to our limiting case (see Eq. A.9) for the slab potential.

BEHAVIOR OF C. ROZZI ET AL OR S. ISMAIL-BEIGI'S CUTOFFS WITH VACUUM45

Note that the different values for v2D

C ( k = 0, k z = 0) obtained from the expression of C. Rozzi et al and S. Ismail-Beigi come from the order of the limits. For S. Ismail-Beigi, k z is set to zero before k → 0. For C. Rozzi et al, k is set to zero, leading to a new formula, before k z is set to 0.

Finally, during the time we were working on the "Selected-G" formalism leading to the slab potential, D. Novko et al [57] proposed a non-diagonal expression for the Fourier transform of the Coulomb potential:

V G 1 G 2 (q) = 4π |q + G 1 | 2 δ G 1 G 2 (6.9) -p G z1 G z2 4π (1 -e -|q+ Ḡ1 |L ) |q + Ḡ1 |L × |q + Ḡ1 | 2 -G z1 G z2 (|q + Ḡ1 | 2 + G 2 z1 ) (|q + Ḡ2 | 2 + G 2 z2 ) δ Ḡ1 Ḡ2 with p Gz = 1 for k z = 2n z π L z p Gz = -1 for k z = (2n z + 1)π L z , n z ∈ Z.
q is an in-plane vector, contains in the first Brillouin zone.

As in our case for the slab potential, this comes from the finite integration limits in the direction perpendicular to the surface, which imposes to consider the Coulomb potential Eq.(6.1) as a function of two variables (z, z ). However, the difference with our work lies in the definition of the reciprocal lattice vectors, built here according to the volume of the supercell.

In all these cases, except our case based on "selected-G" vectors, the procedure is based on the volume of the supercell, which should lead to an influence of the vacuum on the spectra. To test this effect, we have studied the graphene layered systems with these different cutoff procedures.

Behavior of C. Rozzi et al or S. Ismail-Beigi's cutoffs with vacuum

We have coded the formula of S. Ismail-Beigi, with z c = L z /2, but it is equivalent to C. Rozzi for the general case. The only difference between these two cutoffs comes from the order of limits, as explained before. The choice of this formula is justified by the fact that for the limit k -→ 0, (i) we do not apply any cutoff procedure to the ionic charge density, as it is suggested by Rozzi et al, (ii) since we deal with the in-plane case, we consider the limit k -→ 0, after having set k z ≡ 0 in the calculation, (iii) the expression of this potential for vanishinq q is the same as ours (see A.9).

We have calculated the inverse dielectric function for q in the first Brillouin zone:

ε -1 00 (q) = 1 + 4π q2 χ00 (q) (6.10) with χGG (q) = χ 0 GG (q) + G 1 χ 0 GG 1 (q) v2D C G 1 (q) χGG 1 (q) (6.11)
In such a framework, the potential appearing in the expression for ε -1 00 (q) is the 3D Coulomb potential, since one considers that it describes the external perturbation which does not experience the cutoff.

The absorption spectrum is then deduced from ε -1 00 (q) according to Eq. (3.20):

ε M (q) = 1 ε -1 00 (q)
With such a relationship, one gets:

-Im(ε -1 00 ) = Im( 1 ε M ) = Im(ε M ) Re(ε M ) 2 + Im(ε M ) 2 (6.12)
This expression explains the well-known relation between EEL and absorption spectra observed in bulk material: EEL spectrum present a maximum when Re(ε M ) = 0 and when Im(ε M ) is small. This expression has another consequence when one describes an isolated object in a supercell formalism. When increasing the vacuum, one has Re(ε M ) -→ 1 and Im(ε M ) << 1, which leads to -Im(ε -1 00 ) ≈ Im(ε M ). In such a framework, this equivalence is considered as the criteria to ascertain that the vacuum introduced in the supercell is large enough to produce an isolated object.

Results for graphene

We have calculated the EEL spectra for Graphene (1 layer) for in-plane q with two different lengths: q = 10 -5 Bohr -1 and q = 0.078 Bohr -1 .

To evaluate how the vacuum affects the results, the graphene layer has been introduced in supercells of different thicknesses (R2, R3, R4, R5, R10 with the usual denomination see Table 6 1: EELS for one graphene layer in different supercells (R2, R3, R4, R5, R10) without local field effects (NLF-left) and with local field effects (LF-right) for q = 10 -5 Bohr -1 (top) and for q = 0.078 Bohr -1 (bottom). The calculations are done with the standard cutoff_2D.

the number of bands (nband) and the number of G vectors in-plane and out-of-plane which define the size of the matrix to include local fields (npwmat_xy * npwmat_z). They are summarized in table 6.2. Figure 6.1 presents the EEL spectra without local field (NLF -left) and with local field effects (LF -right) for q = 10 -5 Bohr -1 (top) and for q = 0.078 Bohr -1 (bottom). We observe the same behavior for the two values of q. Without local field effects (left panels), the amplitude of the peaks decreases when increasing the thickness of the vacuum. But it is not a simple scaling factor effect coming from the change of the volume, as the peak positions are shifting towards lower energy. With local field effects, no shift is observed, the amplitude of the spectra decreases when increasing the vacuum, and the decrease corresponds to a scaling factor = L z /L mat z , as it is shown in Fig. ( 6.1 -right) by the scaled curves, which are all superimposed.

We emphasize that the behavior of in-plane EEL spectra with vacuum in the presence of the proposed cutoff, is different from the one observed with standard TDDFT (without cutoff, see Fig. 5.4). Without the cutoff, the EEL spectra with (and without) LF shift towards low energy with increasing vacuum. Here, the shift still occurs for non-LF spectra, while in-plane spectra with LF are simply scaled with the volume.

Moreover, the shape of the EEL spectrum with LF for q = 10 -5 Bohr -1 is very similar to the absorption spectrum of graphene or graphite.

Analysis

To understand the behavior of the spectra, we have compared the associated Real and Imaginary part of ε M (q) to EELS, which are related through Eq.(6.12). The spectra are presented in Fig. 6.2 and Fig. 6.3. . The calculations are done with the cutoff_2D defined in Eq. 6.2 or 6.6. ) for q = 0.078 Bohr -1 (bottom). The calculations are done with the cutoff_2D. defined in Eq. 6.2 or 6.6.

Spectra in figures 6.2 and 6.3 have the same behavior, so the description is done for both. It is based on the Eq. 6.12.

For the NLF spectra (Fig. 6.2 and 6.3 -left column), the π-plasmon (corresponding to the region of energy 0-10 eV) has the expected behavior of a plasmon: its maximum is located at the energy where Re(ε M ) crosses zero from negative to positive values, and where Im(ε M ) is very small. The evolution of Re(ε M ), the values are located around 1 when increasing the vacuum , explains the shift of the peak to low energy. The π + σ -plasmon (between 10 en 14 eV) is associated to a non-vanishing Im(ε M ), and a non-zero crossing Re(ε M ). This part of the spectra arises from interband transitions. Above 14 eV, Im(ε M ) decreases, and Re(ε M ), increases from negative value, crosses zero, and when the vacuum is to large tends to 1. This is at the origin of the shift to low energy of the π + σ peak.

For LF spectra (Fig. 6.2 and 6.3-right column), the structures in the EELS appear at energy where Im(ε M ) is small, and Re(ε M ) is above zero, so the peak does not move anymore in energy.

When the spectra are calculated without LF, one has:

ε N LF M (q) = 1 - 4π q2 χ 0 00 (q) (6.13)
Since χ 0 is normalized with the volume of the supercell (Eq. 3.17), the origin of the scaling factor is obvious. The way Re(ε N LF M ) crosses zero shows clearly the reason why the EELS is shifting with the volume of the supercell, and it is independent of the length of q.

When the spectra include LF, in the limit q → 0, the 2D-Coulomb potential behaves like 4πz c /q (Eq. 6.8). Associated with the fact the χ 0 00 ∝ q2 (Eq. 3.17), one sees from the Dyson equation (Eq. 6.11) that χ00 (q) → χ 0 00 (q). Then

ε -1,LF 00 (q) ≡ 1 + 4π q2 χ00 (q) 1 + 4π q2 χ 0 00 (q) (6.14)
and Im(ε -1,LF 00 (q)) ≡ Im( 4π q2 χ 0 00 (q)) Im(ε N LF M (q)). This is confirmed by the comparison of the spectra in Fig. 6.2 (bottom-right for Im(ε -1,LF

00

) and middle-left for Im(ε N LF M )) and explains why the EELS exhibits only a scaling factor for the amplitude.

For larger q (= 0.078 Bohr -1 ), the EEL spectra with LF are more puzzling. Indeed, the χ00 (q) does not tend to χ 0 00 (q). It is confirmed by the comparison of the spectra Fig. ) (bottom-right). In that case, the full Dyson equation (Eq. 6.11) with the general expression of the potential, which behavior with q is more difficult to handle, must be inverted. Nevertheless, the χ00 (q) is still proportional to the volume of the supercell.

With this use of the cutoff_2d, the behavior of the spectra with and without local fields is very different. This is quite surprising, since the local field effects are expected to be rather small in the plane of the layer. Actually, this difference of behavior comes from the way the EELS is calculated without LF. One first calculated ε M with Eq. (6.13), and we invert this quantity to obtain the EELS. It results from this procedure that no cutoff is apply to the EELS without LF.

On the contrary, when one includes LF, one calculates a microscopic matrix ε = 1 -v2D c χ 0 , we invert this matrix to obtain the EEL: such a procedure is equivalent to solve the Dyson equation (6.11), so the EELS is calculated with the cutoff. The resulting quantity at vanishing q in actually ε -1 00 = 1 + v2D c0 χ00 . Then, one renormalises the response function to obtain (see Eq. 6.10).

If we decide to calculate the EELS without LF directly by solving Eq (6.11) with only v2D c0 and calculate ε -1 00 with this resulting response function in Eq. (6.10), then the EELS spectra without LF are very similar to the ones with local fields. This procedure is the correct one.

We have shown that the cutoff procedure with a cutoff parameter z c = L z /2, where L z is the height of the supercell contains a part of vacuum, and provides EEL spectra which are not independent of the vacuum. Nevertheless, without applying the cutoff, the spectra with local field effect were strongly shifted to high energy (see Fig. 5.4). Here, using this 2D cutoff, the spectra are "cured" from the energy shift, and the effect of the vacuum on the spectra with LF reduces to a scaling factor equal L z /L mat z . Since the resulting spectra seem to be in quite good agreement with experiments, and since the amplitude of the measured spectra is not in absolute unit, this procedure has been considered as satisfying.

Behavior of D. Novko et al cutoff with vacuum

Our slab potential is very similar to the one derived by Novko et al [57] and by switching-off the selected-G procedure, we can reproduce the calculation done by Novko et al.

The EELS is calculated with the same expression as ours (we omit q for simplification):

ε -1,LF 00 = 1 + G V 0G χ G0 with χ GG = χ 0 GG + χ 0 GG 1 V G 1 G 2 χ G 2 G (6.15)
Absorption is deduced from ε -1,LF 00 with the standard relation (Eq. 6.12):

ε M (q) = 1 ε -1 00 (q)
We have calculated the EEL and Absorption spectra for the systems define in the previous section.

Results for graphene

The parameter used in the calculations have been summarized in Table 6.2. For the Selected-G method, we used 19 in-plane vectors, corresponding to 3 closed shells (up to (2 0)), and 7

G z vectors, corresponding to G z = 0 ; G z = ±1 * ratio ; G z = ±2 * ratio ; G z = ±3 * ratio.
The calculations done by D. Novko at al [57] considered only LFE out-of-plane (G x = G y = 0), with 71 G z to reach convergence. Associated to a supercell of heigth L z = 5a = 23.26 Bohr (where a = 4.63 Bohr is the length of the in-plane unit cell of graphite), it corresponds to a length of maximum G z = 9.5 Bohr -1 . A complete analysis of the effect of including only out-of plane or in-and out-of plane LF is done in appendix B.

In order to evidence the evolution of the spectra with the vacuum, we put on the same graph the EELS calculated with Novko at al [57] method with different supercells. Figure 6.4 presents the EEL spectra without local field (NLF -left) and with local field effects (LF -right) for q = 10 -5 Bohr -1 (top) and for q = 0.078 Bohr -1 (bottom).

The behavior for the two ranges of q is now different with vacuum:

• For q = 10 -5 Bohr -1 , the spectra are vacuum independent. The effect of the local fields is very small.

• For q = 0.078 Bohr -1 , the spectra are vacuum-dependent. The effect of the local fields is more important than for vanishing q, and is different for NLF and LF spectra. For NLF spectra (botton-left) the two plasmons shift towards low energy when increasing vacuum. For LF spectra, the amplitude is decreasing with increasing vacuum, but no shift occurs. The apparent effect is a scaling factor, but which is not directly related to the scaling factor coming the volume, as shown by the normalized curves in pale lines. ) for q = 10 -5 Bohr -1 (top) and for q = 0.078 Bohr -1 (bottom). The calculations are done with the cutoff proposed by Novko et al [57] using npwmat_xy=1 and npwmat_z=71 .

Analysis

The spectra are calculating with the same expressions we use with our slab potential:

ε -1,LF 00 = 1 + G V 0G χ G0 with χ GG = χ 0 GG + χ 0 GG 1 V G 1 G 2 χ G 2 G (6.16) ε -1,N LF 00 = 1 + V 00 χ N LF 00 with χ N LF 00 = χ 0 00 + χ 0 00 V 00 χ N LF 00 (6.17)
where the colors refer to Fig. 6.5 and Fig. 6.6.

In order to understand the role of the different terms of the potential, we compare the different response functions: χ 0 00 , as well as χ N LF 00 and χ LF 00 , which enter in the EEL spectra (Fig. 6.5 for vanishing q and Fig. 6.6 for q = 0.078 Bohr -1 ). Left pannels contain χ LF 00 with LF only out-of-plane and right pannels contain χ LF 00 with LF in-and out-of-plane.

The difference between the green and blue curves comes from the V GG for (GG ) = (00). One sees from the two cases (npwmat_xy = 1; npwmat_ z = 57) and (npwmat_xy = 19; npwmat_ z = 25) that when we consider only out-of-plane LF, these terms are extremely small, since the curves are the same. Even when we add in-plane LF, the difference is weak (green vs blue). The difference between the blue curve (χ N LF 00 ) and the red curve (χ 0 00 ) allows us to evaluate the role of the inclusion of V 00 , which tends to zero for vanishing q. Origin of the independence of the spectra with vacuum at vanishing q (Fig. 6. 4

-top):

The EELS without LF (Fig. 6.4 -top -left) is calculated according to Eq. (6.17). It gives:

χ N LF 00 = χ 0 00 1 -V 00 χ 0 00 (6.18) ε -1,N LF 00 = 1 + V 00 χ 0 00 1 -V 00 χ 0 00 (6.19) Since V 00 -→ 2πL SC z q and χ 0 00 ∝ 1 L SC z
one has an exact compensation of the scaling factor coming from the height of the supercell.

At vanishing q, one sees (Fig. 6.5) that the χ LF 00 -→ χ N LF 00 , which means that the V GG is dominated by the V 00 term and we can draw the same conclusion as for the NLF case: one has an exact compensation of the height of the supercell (Fig. 6.4 -top -right).

Origin of the dependence of the spectra with vacuum at non-vanishing q (Fig. 6.4 -bottom):

For NLF EEL spectra (Fig. 6.4 -bottom -left), we still can express ε -1,N LF 00 with Eq. (6.19) but in the case of non-vanishing q, V 00 = 2πL SC z q , which explains the fact that the height of the supercell is not compensated and we still see an influence of the vacuum on the spectra.

In the case of LF, the comparison of the blue and the green curves (Fig. 6.6) shows that the inclusion of the V GG terms for non-zero G is not negligible, so a simple analysis is not possible. Nevertheless, the inclusion of LF has an interesting effect, since the spectra with different vacuums are no more shifted, only the amplitude is affected (see Fig. 6.4 bottom-right). This effect is similar to the one observed for the cutoff2D method, but in that latter case, the scaling factor was the height of the supercell (see Fig. 6.1 bottom-right).

With such a framework, the EEL spectra still suffer from the vacuum problem. One notices that the behavior of the cutoff procedure proposed by Novko et al does not behave like the cutoff_2D proposed by Rozzi et al or Ismail-Beigi. In particular, for vanishing q, the results obtained with the cutoff proposed by Novko et al do not depend on vacuum. For larger q, all the spectra have an amplitude dependent on vacuum, but for the cutoff_2D with local fields, the scaling factor can be deduced from the height of the supercell, contrarily to Novko et al procedure.

On the contrary to the calculations done without cutoff, like did by Eberlein et al, the choice of the size of the supercell is crucial to obtain plasmons peaks at an energy in agreement with experimental data. This is done by Eberlein et al [45] with L SC z = 5d 0 = 31.47 Bohr, where d 0 = 6.294 Bohr is the interlayer distance in graphite. With the two proposed cutoff procedures, on one side Ismail-Beigi et al [56] or Rozzi et al [55] or Novko et al [57], the energy position of the plasmon is correct, but the amplitude depends of vacuum. Since the absorption spectrum is calculated according to Eq. (4.12), it will also suffer from the vacuum problem, as it can bee seen on figure 6.3 middle -right).

Dispersion of the plasmon and nature of electronic excitations

We have studied the dispersion of the plasmon for q ranging from 0.039 up to 0.303 Bohr -1 by step of 0.039 Bohr -1 , along the ΓM direction, using our slab potential, with and without the Select-G method, in order to compare to the results of Novko et al [57].

Calculations with Novko et al procedure

The results corresponding to Novko et al procedure calculated with our code (Slab potential (ξ = (-1) n 1 +n 2 ) but no Selected-G method) are presented in Fig. 6.7. To facilitate the comparison, Fig. 6.8 shows the spectra extracted from the paper of Novko et al [57]. Since we have evidenced that this procedure is dependent on the vacuum introduced in the supercell, the calculations were done with the supercell R4 = 25.173 Bohr = 5 d 0 , (where d 0 = 6.294 Bohr is the distance between two adjacent layers in graphite), to be as close as possible as their results (L SC z = 5a = 23.26 Bohr, where a = 4.63 Bohr is the length of the in-plane unit cell of graphite). We have introduced only out-of-plane LF, with npwmat_z = 71. As already mentioned, this does not affect the position of the π-plasmon, which dispersion will be later followed. q ΓM = 0.313 Bohr -1 Re(ε M ) Im(ε M ) EELS Figure 6.7: Real (red lines) and Imaginary (blue lines) parts of ε M and EELS (green lines) for one graphene layer in R4 supercell for different q ranging from 10 -5 , 0.039, 0.078 for the top and 0.157, 0.235, 0.313 Bohr -1 for the bottom plots from left to right. (The first graph is calculated for vanishing q: Re(ε M ) is not shown and is equal to 1). The local field effects are only along z (npwmat_xy = 1; npwmat_z = 71). The calculations are done with the cutoff proposed by Novko et al [57]. The agreement between spectra is very good, which confirms that the formalism of Novko et al is the same as our, except for the Selected-G part (which is not applied in that case). The major differences come from the k-point mesh: our k-point mesh is 40x40x1 shifted grid, which is much smaller than the one used by Novko et al (200x200x1). This is at the origin of the spurious oscillations of the spectra, which prevent us to see some small structures evidenced by Novko et al in the spectra (specially the peak they attribute to σ -→ σ * 2 transitions). The purpose of our work was not to reproduce their results, but to compare to their work the calculations done with the Selected-G method, and the precision we reach with our parameters is good enough.

Calculations with slab potential and Selected-G method

Since we have demonstrated that the slab potential with the Selected-G allows us to provide results which are independent of the height of the supercell, we applied it to compare to Novko et al's paper. We introduced in-and out-of-plane LF, with (npwmat_xy = 19; npwmat_z = 7) 1 . Results are presented on Fig. 6.9 and can be compared to Fig. 6.7. For the absorption spectra (6.9-top), calculated with ε M = 1 -4π/q 2 χ00 , the difference with Novko et al is more important than for the EEL spectra. In particular, the peak above 10 eV in our absorption spectra (corresponding to the π + σ structure) keeps the same shape for all the range of q; in the calculations of Novko et al, the spectral shape clearly evolves. The first peak of absorption presents the expected behavior coming from interband transitions, diverging at 0 eV for vanishing q due to the metallicity at K point, and dispersing when q increases. This was also the case for Novko et al spectra.

The EELS (Fig. 6.9 -bottom) present mainly the same features as the one calculated with Novko at al. The π plasmon is similar, with the appearance of the structures around 5 eV. The π + σ plasmon is located in the same energy range, the slight modification of the shape can be explained by the introduction of in-plane LF, since we have seen that they were responsible of the shift of the spectral weight to high energy (see Fig. The main difference comes from the amplitude. We recall that the results of Novko et al depend from the vacuum chosen to build the supercell, which has consequences on the amplitude, and on the shape of the absorption spectra. This is not the case of our calculation, which is totally independent of the choice of the supercell.

Dispersion

In order to compare further the behavior of the spectra calculated with the different procedures (Novko et al and our method, which differ from the choice of G-vectors of the microscopic dielectric matrix), we report the dispersion of the main peak of EEL and the corresponding peak in absorption spectra in the 0-10 eV energy range. Figure 6.10 presents data from different origins: black squares are experimental data (from [59]), red dots correspond to dispersion of π plasmon, blue dots follow the dispersion of corresponding structure for Im(ε M ) extracted from Novko et al paper [57]. [59], red dots: dispersion of π plasmon, blue dots: corresponding structure of the absorption extracted from [57]. Cyan (magenta) plain circles + dashed lines corresponds to the π plasmon (absorption) dispersion taken from our calculations, done on a 100x100x1 mesh, with Selectd-G method (see Fig. 6.9).

Considering the results obtained by Novko et al [red (blue) dots on Fig. 6.10], the authors distinguished in the dispersion of the EELS two regimes [57]: for large q, they evidence a linear dispersion, which is in very good agreement with the experimental data of [59], and cannot be extended to vanishing q. For vanishing q, their calculations allows them to conclude that the dispersion is proportional at q2 , similarly to the one of Im(ε M (q)).

Cyan (magenta) plain circles + dashed lines on Fig. 6.10 corresponds respectively to the π plasmon (Im(ε M )) dispersion taken from our calculations done using Selectd-G method (from Fig. 6.9). The agreement between the π plasmon dispersion with both the calculations of Novko et al and the experimental data from Liou et al [59] is very good. Concerning the absorption dispersion, we do not recover the dispersion of Novko et al. Nevertheless, in the absence of experimental mesurements of absorption of graphene, it is difficult to conclude which dispersion would be correct. The dispersion of absorption spectra taken from our calculations is constant for vanishing q, then it has a q2 behavior, then a linear dispersion, parallel to the one of the EELS up to 0.5 a.u. Finally, the curves seem to diverge: the two lines don't converge to each other, as calculated by Novko et al.

Discussion

This discussion about the dispersion of π plasmon of graphene takes place in the context of the nature of electronic excitations at the origin of this spectrum. A quasi-linear dispersion of EELS has been several times measured and explained as a consequence of the Dirac electrons of graphene [58,60,61]. More recently, Nelson el al. [62] claimed that the two features measured in EEL spectrum were not plasmons but single particle interband transitions. Finally, using a very high q resolution TEM, Liou et al [59] evidenced a √ q dispersion of EEL spectra, in agreement with the one expected from the pure 2D plasmon [63].

Novko et al [57] attribute the behaviour of the two dispersions curves for the EEL and absorption spectra [red (blue) dots on 6.10] to the evolution of the screening, calculated as the ratio D(q, ω) = EELS(q, ω)/ Im(ε M (q, ω)). This modification of the screening would be at the origin of the change of nature of the excitations, from single-particle excitations to collective ones when D(q) moves away from 1.

They conclude from their calculated dispersion that for q 0.03 a.u., the spectra are plasmons, with the quasi-linear dispersion linked to the Dirac electrons of graphene [58,60,64]. At vanishing q, they assign their spectra a purely single particle character. They based their conclusion on two arguments: (i) D(q, ω) ≈ 1 and (ii) the calculated dispersion is proportional to q2 similarly to the one of Im(ε M (q)), as a consequence of the band structure of graphene around the M point of the Brillouin zone. This single particle interband transitions nature is in agreement with the conclusion of Nelson et al [62]. Nevertheless, the conclusion of Nelson et al is based on the fact that the Re(ε M ) extracted from their loss measurement never reaches zero. This observation can also be done on Novko et al calculations, even in the cases where they conclude to the collective origin of the EEL spectra.

One also notices that the quadratic dispersion obtained by Novko et al's calculations allows them to reach the two first experimental points, but no intermediate measured point exists in that range to ascertain this dispersion.

Liou et al [59], who did the measurements, do not draw the same conclusion: they show that their data follows a √ q law, expected from the pure 2D plasmon [63] up to 0.8 Å -1 .

Zoom on the dispersion for different ranges of q

We studied the behavior of the dispersion of the spectra calculated with the slab potential and Selected-G method, to compare with the proposed dispersions.

For very small values of q, the calculations are done using k.p theory, since we cannot reach such a thin resolution on the grid. We have checked that for q = 0.016 a.u., the results obtained on the grid were equivalent to the one using k.p theory. In this range of q, the absorption spectra can be superimposed: they do not experience any dispersion. The dispersion of the structure corresponding to the π plasmon is reported with red line in Fig. 6.11. The green line is a fit with a function proportional to q2 . To improve the resolution of the spectra, the calculation has been done with the R5 supercell, using a grid 100x100x1 with a broadening of 0.05 eV. (We recall that our calculations are independent of the the choice of the supercell). The excellent agreement between the two curves confirms that the dispersion of the plasmon is quadratic up to 0.018 a.u.. According to Novko et al, this should be the signature of interband transitions, around M point of the Brillouin zone.

To check this hypothesis, I have compared the response functions at the origin of the EEL and absorption spectra. The plots at the top of Fig. 6.12 compare the EELS spectrum (magenta) with the EEL spectrum truncated to the first term of the summation (cyan) [see Eq. (6.16)] for q = 10 -5 Bohr -1 (left), for q = 0.016 Bohr -1 (right). The perfect superposition of the spectra allows to conclude that only the head of the matrix χ 00 contributes to the spectra in this range of q.

= 1 + G V 0G χ G0 (magenta)
with the EEL spectrum truncated to the first term of the summation 1 + V 00 χ 00 (cyan) [see Eq. (6.16)] for q = 10 -5 Bohr -1 (left), for q = 0.016 Bohr -1 (right). (bottom) The response functions leading to EELS (green) and to absorption (black) for the corresponding q. The spectra have been calculated on the (100x100x1) mesh with npwmat_z=7 and npwmat_xy=19.

The bottom plots show the comparison of χ 00 (green) with the response function at the origin of the interband transitions (black) for the respective q. It is clear that for vanishing q (left), the EEL spectrum is equivalent to interband transitions spectrum, but for q = 0.016 Bohr -1 , which corresponds to one of the largest q still included in the quadratic dispersion, the EELS spectrum is not exactly given by interband transitions. Nevertheless, the energy range of the peaks are very similar. For q 0.02 a.u., I plot the dispersion of the front edge of π transitions for EELS, calculated with the Slab potential and the Selected-G method on R5 supercell, on the grid 100x100x1 with a broadening of 0.05 eV. The calculated points are in cyan, and the data are fitted with √ q function (red line), as well as a linear function (green dashed line). Black squares are extracted from the Liou experiment [59].

EELS slab+selectedG 8√

q + 3.25 10.6 -q + 4.7 Liou-paper-experim-π plasmon For q ranging from 0.02 a.u. up to 0.24 a.u., (see Fig. 6.13), the √ q function (red line) fits very well the calculated points (cyan). This result is on very good agreement with the conclusion drawn by Liou et al [59], who observed a √ q behavior up to √ q ≈ 0.7 given in Å -1 , corresponding to q ≈ 0.26 a.u.. This √ dependence with the transferred momentum suggests that my calculations reproduced the behavior of plasmon with 2D character, in agreement with the high resolved EELS measurements of these authors.

For q 0.25 a. u., I recover the linear dispersion, in agreement with [59], and with others authors in previous works [58][59][60].

Origin of the plasmon in graphene

To further check this hypothesis, I compare ε M and the EEL spectra for this range of q (see Fig. 6.14). Zooming on the real part of ε M , we have added a line at -1, which is the definition of a surface plasmon for a flat thin film in vacuum [65,[START_REF] Raether | Excitation of Plasmons and Interband Transitions by Electrons[END_REF]. for one graphene layer in R5 supercell for different q ranging from 0.016 up to 0.188 Bohr -1 by step of 0.016 Bohr -1 . The calculations are done on the (100x100x1) mesh, using slab potential and Selected-G method with (npwmat_xy = 19; npwmat_z = 7).

DISPERSION OF THE PLASMON AND NATURE OF ELECTRONIC EXCITATIONS63

One can see that for the very first values of q (0.016 up to 0.094 a.u.), the EEL spectrum shifts from an energy similar to the one of Im(ε M ) to an energy corresponding to Re(ε M ) = -1, and keeps this value for the range of q going from 0.1 up to 0.2 Bohr -1 . This means that the loss function calculated in the graphene layer behaves like a surface plasmon. Such a result is consistent with the √ q behavior of the dispersion for q -→ 0.

We did the same analysis for q 0.25 a.u., where the dispersion was found linear (see Fig. 6.15). for one graphene layer in R5 supercell for different q ranging from 0.2 up to 0.7 Bohr -1 . The calculations are done on the (40x40x1) mesh, using slab potential and Selected-G method with (npwmat_xy = 19; npwmat_z = 7).

The first plot shows that the π plasmon can still be associated with a 2D-plasmon. Very rapidly, when increasing q, one sees that Re(ε M ) does not cross anymore neither the -1 nor the 0 line in the energy range 0 -10 eV, and the EEL spectrum has features at energy close to Im(ε M ): the EEL spectrum is still given by inter-band transitions.

One should mention a quite puzzling result: the linear dispersion for large q has been attributed to the presence of Dirac electrons. The reason why it occurs from 4 eV, and not 0 eV, as suggested by the band structure around K point of the Brillouin zone, was explained by the introduction of local field effects, which mix the transitions and transfer the linear dispersion up to 4 eV (see Fig. 3 from [60]). This conclusion was drawn in a context where the isolated graphene layer in supercell formalism was achieved by increasing the size of the supercell, with the criteria that the EELS and the absorption should be the same. For this reason, the EELS spectra without LF were described by the χ 0 00 response function, which follows the dispersion from 0 eV, and the inclusion of the LF was at the origin of the shift of the dispersion from 4 eV. However, the spectra calculated with Novko procedure or with our slab potential and selected-G method still recover the same linear dispersion, but a quick look at the comparison between EELS spectra with and without LF (Fig. 6.16) shows that the EEL spectra without LF (green -dashed) are not described by the χ 0 00 , [Im(ε M ) without LF (blue -dashed)] but are very similar to the ones with LF (green -solid). The solution of the Dyson equation with only the V 00 term χ N LF 00 :

χ N LF 00 = χ 0 00 + χ 0 00 V 00 χ N LF 00 modifies the spectra strongly enough to push it above 4 eV.

In our calculations with the slab potential and selected-G method, the χ 0 00 has not been modified, and has the same features as in [60] (see Fig. 6.14). So it seems that the 2D character of the Coulomb potential is more at the origin of the modification of the spectral weight of the spectrum than the mixing of transitions induced by the local fields. The presence of transitions arising from the K point in the peak at 4 eV should be checked to confirm the origin of the linear dispersion [60]. for one graphene layer in R5 supercell for different q ranging from 0.2 up to 0.7 Bohr -1 . The calculations are done on the (100x100x1) mesh, using slab potential and Selected-G method with (npwmat_xy = 19; npwmat_z = 7). The corresponding dashed lines are the spectra without LF.

Summary

The different behaviors appearing in the dispersion of the main structure of the so-called π plasmon in graphene, for in-plane transferred momenta, are not contradictory, since they seem to correspond to different ranges of q, summarised in table 6.3. range of q (a. u.) q 0.018 0.02 q 0.25 0.25 q dispersion q2 √ q q Table 6.3: Summary of the different dispersions of π plasmon in graphene.

A careful analysis of the q range of the measured spectra, depending on the resolution of the setups, which are regularly improving, indicates that none of the results are in contradiction. The most problematic range is of course the one of the smallest q, since the quadratic dependence arose surely from calculations [both Novko et al [57] and ours (Slab potential with Selected-G method)], but is not really confirmed by experimental measurements due a lack of data. The spectrum at vanishing q can be attributed to inter-band transitions, (as suggested by Novko et al [57] and Nelson et al [62]) but it is quite rapidly no more true.

The nice √ dependence tends to prove the 2D nature of the plasmon, in the range 0.02 q 0.25 a.u. This result was pointed out by Liou et al [59].

Even if the linear dispersion for large q is confirmed by all the formalisms proposed, and very well reproduced both qualitatively and quantitatively the experimental data, the relation between the Dirac electrons of graphene and the inclusion of local fields should be further investigated. With our formalism, it seems to arise from inter-band transitions, with a small influence of local fields. The presence of transitions arising from the K point of the Brillouin zone must still be checked.

Conclusion

Our slab potential has been defined to treat an isolated slab system in a framework of 3D codes, based on plane-waves and periodic boundary conditions. In order to obtain results independent of the vacuum, auxiliary response functions in real space equal to the response functions of the isolated slab, and periodic with the periodicity L mat z has allowed the development of the formalism called "Selected-G".

Due to the periodicity, the integration along the direction perpendicular to the slab is limited to the thickness of the matter. This acts as a cutoff procedure.

In this chapter, we have studied the behaviour of previously proposed cutoff procedures [55][56][57]. We have seen that all these cutoff procedures, which proposed to cut the Coulomb interaction in the middle of the supercell still contains a portion of vacuum, and consequently give spectra depending on the vacuum.

Among the works done with the cutoff procedures, Novko et al [57] derived a non-diagonal expression of the Coulomb potential in reciprocal space similar to ours, except to that the basis of G vectors is the complete one, defined on the height of the supercell, and not the height of the matter, as we do in the Select-G method. They studied the in-plane dispersion of the plasmon in graphene, and compared to the recent measurements of Liou et al [59], and Nelson et al [62]. Their calculations reproduce very well the dispersion measured by Liou et al. Nevertheless, they do not draw the same conclusions as the authors of the measurements.

Liou et al evidenced from their own data a √ q behavior, in agreement with a 2D plasmon, followed by a linear dispersion attributed to the Dirac electrons.

Novko et al identify two regimes in the dispersion curve: one at very small q, proportional to q2 , attributed to inter-band transitions (in agreement with Nelson et al), and one at large q corresponding to the linear dispersion. The √ q dependence (evidenced by Liou et al) is, according to Novko et al, simply the junction between the two regimes. They also conclude from their calculations that the plasmon has an inter-band transitions origin for very small q, as well as for large q, when the dispersion is linear. In between these two regimes, the EEL spectrum originates from collective excitations. This conclusion is done in the light of the comparison with the absorption spectra calculated as the inverse of the microscopic inverse dielectric function, which depends from the height of the supercell.

Using our formalism ("Selected-G" + slab potential), I did the similar calculations and compared the dispersion of the in-plane plasmon in graphene to their results and to the experiments of Liou et al, and Nelson et al [62].

My results show different ranges of dispersion: for very small q (< 0.02 a.u.), the plasmon behaves like interband transitions, with a q2 dependence, in agreement with the conclusion of Novko et al and Nelson et al. Nevertheless, this conclusion arises mainly from calculations, it is not really supported by measurements due too the lack of data in this region of very small q. I also recover the linear dispersion at q > 0.25 a.u., but is arises from collective excitations. Between these two regimes, I evidenced a √ q dispersion, and the comparison with the Re(ε M ) tends to show that it is a 2D plasmon.

Chapter 7

Absorption spectra for interacting slabs

Inverse Effective Medium Theory

We have developed a formalism to calculate the electronic excitations spectra for an isolated slab. Now, we would like to study the absorption spectrum of an array of finite size objects. In that case, the absorption of one object of the array should be modified by the interaction with the surrounded electronic density coming from the presence of the neighbouring objects. Moreover, the neighbouring objects should also experiment an induced electronic density, due to the long range of the Coulomb potential. Such an effect should be accounted for by the local fields, and the expected information should be contained in the TDDFT spectrum of the repeated supercell, namely the regular TDDFT.

But we have already shown that the regular TDDFT in supercell formalism suffers from what we have called the vacuum problem. Indeed, the regular TDDFT is equivalent to Effective Medium Theory (EMT) with vacuum (see Fig. 4.5) [START_REF] Tancogne-Dejean | [END_REF]. Such a spectrum cannot be the expected absorption. In particular for the component perpendicular to the surface of the slab, one cannot imagine that if one slab has an absorption peaked at 4 eV, the absorption of two interacting slabs will be located above 12 eV.

In the paper where we have evidenced that the regular TDDFT in supercell formlism behaves like EMT with vacuum [START_REF] Tancogne-Dejean | [END_REF], we also distinguished the two limiting cases of minimum screening and maximum screening.

We showed that the minimum screening case was also associated with the small local field effects situation, and was suitable to describe the spectra for in-plane components. It reads

ε SC M = f + (1 -f ) ε bulk M (7.1)

69

The maximum screening corresponds to large local field effects and allows us to interpret the out-of-plane component and can be written as:

1 ε SC M = f + (1 -f ) ε bulk M (7.2)
The label SC means supercell and corresponds to the spectra calculated with the regular TDDFT procedure. f is the proportion of vacuum in the supercell and is related to the quantity ratio I used to define the supercell by f = 1 -1/ratio. The results presented in Ref. [START_REF] Tancogne-Dejean | [END_REF] have been obtained using the macroscopic dielectric function of the bulk system for ε bulk M .

The idea is that, among the effective medium theory with vacuum, the regular TDDFT spectra should contain the Coulomb interactions (local field effects) between the neighbouring slabs, which should modify the response as compared to the isolated one. In the ref. [START_REF] Tancogne-Dejean | [END_REF], we used the macroscopic dielectric tensor of bulk silicon ε b ulk M to mimick the regular TDDFT. But, according to the geometry of the repeated supercell, the quantity which is "averaged" with vacuum is the macroscopic dielectric tensor of the "interacting" slab ε I nteract M .

In order to "cure" the spurious influence of the vacuum included in the supercell, we propose to apply a procedure called Inverse Effective Medium Theory (IEMT), to extract the macroscopic dielectric function of the so-called interacting slab, that is a slab which has experienced the local field from the neighbouring ones, but not merged in an artificial vacuum box. For this, we reconsider the equations (7.1) and (7.2) by replacing the bulk macroscopic dielectric function with the corresponding one of the so-called interacting slab. It comes for the minimum screening case:

ε SC M = f + (1 -f ) ε InterSlab M (7.3)
and for the maximum screening case:

1 ε SC M = f + (1 -f ) ε InterSlab M (7.4)
Inverse Effective Medium Theory (IEMT) stands for the inversion of these equations: minimum screening case:

ε InterSlab M = ε SC M -f 1 -f (7.5) maximum screening case: 1 ε InterSlab M = 1 1 -f * ( 1 ε SC M -f ) (7.6)
To test the validity of the model, we have chosen to study hBN slabs as a test case. This choice relies on the fact that the slabs of hBN have an insulating character. We do not need to perform a reconstruction of the surface, leading to the presence of surface states, that could modify the absorption spectrum. They present a large gap, preventing numerical instabilities for small frequencies.

Absorption spectra of hBN slabs 7.2.1 Results with the standard TDDFT

I have calculated the spectra for slabs of 1, 2, 3, 4 and 8 layers in supercells with different ratios (see table 7.1). The different parameters are summarized in tables 7.2, 7.3, 7.4. As already mentioned in ref [START_REF] Tancogne-Dejean | [END_REF], and as it can be easily seen from Eq. (7.5), the absorption spectrum of each slab trivially shows a scaling factor 1/ratio, corresponding to the supercell in which it is introduced. This is directly related to the fact that in the absence of local field (let us remind that in these systems, in-plane local fields are almost negligible), the absorption spectrum is given by the χ 0 00 response function, which is normalized by the volume of the supercell. As a consequence, the position of the peaks in the spectra is unchanged. For the out-of-plane component, we have shown that the regular TDDFT, as a consequence of the local field effects, modifies the spectra. This was particularly dramatic for silicon slabs [START_REF] Tancogne-Dejean | [END_REF], where the position of the peak was pushed towards high energy for more than 10 eV. The effect is less spectacular for hBN systems, since the spectra seem to be "only" scaled, but they are still dependent of vacuum. Moreover, the scaling factor is not related to the ratio of the supercell.

nb-layers R2 R4 R5 R10 1 X X X X 2 X X X X 3 X 4 X X X 8 X

IEMT: convergence effects

To extract the spectrum of the so-called interacting slab, I applied the IEMT for the out-ofplane component of the hBN systems. 7.3 shows the out-of-plane component of Im(ε M ) resulting from IEMT applied to the 2-layer slab into R4 supercell (blue line). I compared it to the bulk spectrum (black line) and to the isolated 2-layer slab (red line). The pre-edge is exactly superimposed to the one of the isolated slab spectrum. Above 17 eV, the IEMT spectrum also reproduces the isolated slab, and the features are very similar to the one of the bulk. Between 8 and 17 eV, the amplitude of the peaks is different from both the bulk and the isolated slab. The peaks are located between these two spectra, in agreement with the idea that interacting slabs have an absorption which "feel" more electronic interactions via the Coulomb potential than the isolated slab, and less than the bulk material.

Nevertheless, the results presented here are extremely sensitive to the number of G z vectors and a careful analysis of the convergence for the local fields has been performed.

Influence of number of G z vectors for local fields

The convergence tests have been done on the 2-layer system for different vacuum sizes. The k-point mesh is very coarse (5x5x1), so the spectra present spurious wiggles, but this does not affect the results. The in-plane local fields are based on npwmat_xy=13.

For each system, the height of the supercell can be associated with a number of "equivalent" layers, defined by Nb_eq_layer * ratio. Then we adjust the value of the maximum G z to be a multiple p of Nb_eq_layer (see table 7.5). This p value aims to represent the frequency of the spatial electronic density fluctuations between the layers building the slab. p(Nb_eq_layer) R2( 4) R4( 8) R5( 10) R10( 20 Figure 7.4 shows the absorption spectra along z-direction after applying the Inverse Effective Medium Theory (focussed on the energy range 4 -13 eV) of the 2-layer slab for the various supercells. The different panels correspond to the different values of the p used to define max G z to include out-of-plane local fields. The black curve is the bulk spectrum and the yellow curve is the one of the isolated slab.

One sees from the top-left panel of Fig. 7.4 (corresponding to p = 3) that we achieved a good agreement for the full range of energy, except the peak located at 9.5 eV. Indeed, the R2, R4 and R5 supercells give a peak whose energy position and amplitude is very similar, when the R10 supercell gives a different result, since the peak is more intense and with a spectral weight at lower energy, that is closer to the bulk hBN than to the isolated slab absorption spectrum. This observation is quite surprising, since one would have expected the R10 supercell to contain less interactions between slabs than the supercell with a smaller height. To check if this observation was physical or spurious, I have increased the number of G z vectors to include local fields in the perpendicular direction. By increasing the p value up to 6 [Fig. 7.4 (bottom-right)], the spectra resulting for IEMT on the R10 supercell converge, and finally all the spectra resulting from IEMT on all the supercells are superposed. The p parameter depends on the ratio of the supercell which is contradictory with the fact that is should represent the fluctuation of electronic density between the layers of the slab. Actually, this is not completely true. For all the supercells, the convergence is better achieved for p = 6, but to be able to do the calculations with in-and out-of-plane local fields, on a k-points grid dense enough in plane to obtain a spectrum without spurious features, we must limit the p value to 3. So, this analysis shows that we cannot use a too large supercell to converge this particular peak, which corresponds to transitions bewteen states very sensitive to small spatial electronic density fluctuations. We also mentioned that the large intense peak is a physical result, and does not come form the coarse k-point mesh used for convergence tests (see Fig. 7.

3).

A comparison of spectra for different supercells [Fig. 7.4 -(last panel)] shows that all the spectra are well superposed, except the R2 one. Since I did the convergence tests on the 2-layer slab, a supercell with ratio = 2 does not contains so much vacuum and we can wonder if the slabs are really isolated at the level of the ground state, that is if their is no more dispersion of the electronic structure along the z-direction, and if a k-point mesh NxNx1 allows to describe properly the system. To check this, I did a convergence study with the number of divisions along the z-direction. Figure 7.5 shows that the spectrum calculated for the N z = 1 k-point mesh is not converged. The vacuum introduced inside the supercell in this R2 case is not large enough to isolate the 2-layer slab. The electronic structure still exhibits a dispersion which is correctly captured with N z = 4.

Influence of k-point mesh

Applying IEMT to the 2-layer R2 supercell calculated on the grid 15x15x4, I do not recovered the same spectrum as the 2-layer R4, and the 2-layer R10 supercells (Fig. 7.6). The slight difference observed is a physical effect, and not a problem of convergence. Indeed, the level of convergence can be seen on the negligible differences at the maximum of the peak at 10 eV for R4, R5, and R10 supercells: all the spectra are converged in amplitude better than 5%. I also checked that for 2-layer R4 was converged with N z = 1. Such a vacuum is large enough to isolate the slab at the level of the ground state. As a consequence, the difference between the R2 system and the other ones should be the signature of the interaction between slabs. As shown by Fig. 7.4 (bottom-right), all the IEMT spectra are the same (except the R2, for the reasons explained just before): no effect due to the height of the supercell is introduced by the TDDFT. The extension of the supercell does not modify the interactions felt by the slabs. This result is quite puzzling. We would have think that by increasing the height, the spectra would have evolved toward the spectrum of the isolated slab. The only case where we see the presence of interactions is when the layers are not isolated even in the ground state.

Inverse Effective Medium Theory with vacuum

Applying the IEMT to the full converged spectra according to equations (7.5) and (7.6), we obtained the spectra on figures 7.7 and 7.8 respectively.

Figure 7.7 presents in-plane absorption spectra for the different hBN slabs. In blue, I plotted the IEMT applied to the regular TDDFT calculation to mimick the interacting slab and I compared to the bulk (red) and isolated slab spectra (green).

The IEMT spectra for in-plane components are almost exactly superposed on the ones of the isolated slab, which are very similar to the one of the bulk material. The only part of the spectrum which exhibits a (very) small difference with the bulk is between 10 and 14 eV. It reduces when the thickness of the slab increases. On in-plane component, this is the only part of the spectrum which involved the part of electronic structure sensitive to the incomplete stacking of hBN layers. In blue, I plotted the IEMT applied to the regular TDDFT calculation to mimick the interacting slab and I compared to the bulk (red) and isolated slab spectra (green).

The following observations can be done:

1) IEMT allows us to recover a spectrum whose absorption is in the energy range of the bulk and of the isolated slab, contrarily to what observed in [START_REF] Tancogne-Dejean | [END_REF].

2) The IEMT applied to the 1-layer slab does not reproduce the pre-edge. This means that this structure comes from states resulting from the overlapping of the p z orbital between the adjacent hBN layers.

3) The pre-edge at 5 eV obtained with IEMT is in the same energy as for the bulk and the isolated slab. The amplitude is in between these two cases for all structures, except the 3-layer. It could be due to the incomplete stacking in case of an odd number of layers. 4) All the supercells give the same results for the peak at 9.5 eV shown to be very sensitive to the small spatial electronic density fluctuations (see Fig. 7.4): its energy position is inbetween the one of the bulk and the isolated slab for thin slabs. Its amplitude is also much more intense. Then it progressively recovers the energy position and the amplitude of the bulk and isolated slab for thick slabs.

5)

The only case where we saw interactions between neighbouring slabs is the one of the 2-layer system with vacuum defined by two times the 2-layer thickness (R2), which reveals too small to isolate the slabs also at the level of ground state. The TDDFT does not add new interactions.

6) The 8-layer slab gives the absorption spectrum the closest to the bulk: the pre-edge structure almost recovers the amplitude of the bulk. As already noticed from the analysis of the 1-later spectrum, it is the signature of the interaction between layers inside the slab: 8-layer slab almost recover the electronic structure of the bulk. 7) All these conclusions are certainly valid for hBN slabs, as well as graphene ones, or more generally systems where the layers are feebly interacting (van der Waals interactions).

For the case of silicon, where the interaction between layers is much larger, since they are due to covalent bonds, one can wonder if the results are the same.

Absorption spectra of silicon slabs

I have studied different slabs of silicon, introduced in different supercells. They correspond to a (1x2) surface reconstruction. The unit block is composed of 4 bilayers of silicon atoms (8 atoms). Then this block is repeated N times to define the slab. The thickness of the slabs considered, as well as the number of atoms, according to their name is summarized in table 7 These slabs have been introduced into supercells of different heights : the second number of the name will correspond to the thickness of vacuum, in the same unit that the thickness of matter. As an example, 4_4 means 16 bilayers of silicon atoms in a supercell which contains the same thickness of vacuum. 4_8 corresponds to 16 bilayers of silicon atoms in a supercell which contains the double thickness of vacuum.

The calculations are done on a 16x8x1 k-point mesh, with the parameters summarized in table 7

IEMT for the slab "4" (thickness 21.72 Å)

I used the "4" slab to illustrate different effects.

As already mentioned, for the in-plane component (Fig. 7.9), the regular TDDFT has the effect to scale the spectrum by a factor L mat z L SC z . This is a consequence of the negligible local fields for in-plane excitations. The isolated slab and the interacting slab (resulting from IEMT) are similar. For the out-of-plane component (Fig. 7.10), the regular TDDFT pushes the spectrum to 12 eV, with a strong reduction of the amplitude. Moreover, the peak is dependent of the size of the supercell as can be seen from the comparison of the two panels of Fig. 7.10.

The inverse effective medium theory applied to the regular TDDFT gives the blue spectrum peaked to 2.5 eV, that is at much smaller energy than the bulk absorption. This effect is independent of the size of the supercell, since the spectrum extracted from the supercell 4_4 gives the same result as the one extracted from the supercell 4_8.

The comparison of the IEMT for the two supercells (embedded the same slab) (Fig. 7.10) exhibits no differences. This confirms that the conclusions drawn for van der Waals interacting layers are also valid for systems with much stronger bonds: we don't see the influence of interacting slab at TDDFT level. The fact that the spectrum extracted from the IEMT applied to the standard TDDFT calculations, is located at much smaller energy than the bulk absorption peak was, at the very begining, a surprise, and our first conclusion was that IEMT could not be applied to extract the absorption spectrum of the so-called interacting slab, namely to correct the spurious effect of the vacuum introduced in the supercell.

The reason why in slab of silicon, we recover an absorption edge at lower energy than for the bulk is actually due to the fact that the silicon slabs need to be build from a reconstructed (2x1) surface, in order to avoid the metallic character of the (1x1) surface. The (2x1) reconstruction, among the semi-conducting behavior, introduces surface states which are located in the gap of the surface projected bulk band structure, and are at the origin of the low energy transitions. Figure 7.11 presents the band structure of the surface projected bulk in blue. To obtain this result, I have used the 3 blocks (24 atoms) built on the (2x1) bulk unit cell, namely without vacuum and without relaxing the atomic positions. By calculating the band structure along the direction Γ -X -S -Y -Γ -S for different cuts along the k z direction, I have obtained the surface projected bulk states. The reason of the use of 3 blocks was simply to increase the density of the lines with a reduced number of cut along k z . The red lines correspond to the band structure calculated for the 4_4 relaxed slab. We clearly see the appearance of bands inside the gap of the blue area. These surface states are at the origin of the pre-edge peak visible for the in-plane component of the 4_4 silicon slab (Fig. 7.9), both for the interacting and isolated slabs.

Nevertheless, even if the low energy position of the absorption peak of the so-called interacting slab can be explained by the transitions arising from the surface states located in the gap of the surface projected bulk states, the result still surprises us. Indeed, the absorption peak of the interacting slab is located at much lower energy than the isolated slab (Fig. 7.10 -(cyan)), where the surface states are also present. This result is quite different from what we have observed for the hBN system. The effect coming from the surface states seems us to be unexpectedly large.

Study of slabs with different thicknesses

I have studied the absorption spectra for three slabs, corresponding to N=3, 4, 6 as defined in table 7.6. Figures 7.12 and 7.13 show the absorption spectra for the three isolated slab as well as for the three interacting slab (IEMT applied to standard TDDFT), for in-plane and out-of-plane components respectively.

Concerning in-plane component (Fig. 7.12), the solid blue, orange and green lines correspond to the absorption spectra of the three isolated slabs (calculated with the Select-G method and the slab potential) : N = 3, N = 4, N = 6 respectively. (The vacuum introduced in the supercell is, for each system, equal to the matter (ratio = 2), but it is meaningless, since the calculations are independent of vacuum). The dotted blue,orange and green curves (Fig. 7.12) correspond to the Inverse Effective Medium Theory applied to the standard TDDFT results, for the three slabs. The isolated and so-called interacting slabs give equivalent spectra. The peak between 0 and 2.5 eV comes from the surface states. The peak between 2.5 eV and 10 eV, coincides in energy with the bulk spectrum (red). One clearly sees that the ratio between the surface and bulk states reduces when the thickness of the slab increases. The amplitude of the bulk peak tends to the bulk amplitude. Bulk "3" Isolated slab (x) "3" IEMT (x) "4" Isolated slab (x) "4" IEMT (x) "6" Isolated slab (x) "6" IEMT (x) Figure 7.12: Comparison of the in-plane absorption of the silicon slabs for isolated and interacting (IEMT) slabs of different thicknesses. "3" (16.29 Å): isolated (blue solid); interacting (blue dotted). "4" (21.72 Å): isolated (orange solid); interacting (orange dotted). "6" (32.59 Å): isolated (green solid); interacting (green dotted). Bulk (red).

In solid blue, green and orange (Fig. 7.13) are plotted the out-of-plane absorption spectra of the three isolated slabs (calculated with the Select-G method and the slab potential) : N = 3, N = 4, N = 6 respectively, (with ratio = 2). As already observed from surface results, larger the thickness, the closer to the bulk the spectrum.

The dotted blue, green and orange curves (Fig. 7.13) correspond to the Inverse Effective Medium Theory applied to the standard TDDFT results, for the three slabs. The smaller the thickness, the lower energy for the absorption peak. This can be explained by the fact that the absorption spectrum arises from transitions from the surface states, and from bulk states. When the slab is thin, the relative ratio of the surface states in the total absorption spectrum is more important: since the spectral weight at low energy comes from surface states transitions, it justifies that the "3" slab of thickness 16.29 Å has its absorption peak at lower energy that the "4" one of thickness 21.72 Å, itself at lower energy than the "6" one, of thickness 32.59 Å. It also explains why the spectra are at lower energy than the bulk.

The absorption spectra of the interacting slabs are very different from the ones of the corresponding isolated systems: the spectra tend to converge to the bulk one when increasing the thickness of the slab, but the isolated systems exhibit a red shift while the so-called interacting ones experiment a blue shift. Even if most of the trends can be explained by the presence of the surface states, the comparison between the so-called interacting slabs and the isolated slabs is still puzzling. Indeed, the surface states are also present in the isolated slab calculations: they are visible at 2.5 eV for the "3" slab (magenta), and decrease with the increasing thickness. The fact that the spectra of interacting slabs exhibit the presence of electronic interactions arising from the electronic density of the surface states, is an expected result, but such a huge effect of the Coulomb interactions between adjacent layers surprises us. A way to understand this effect could be to calculate the spectrum of two interacting layers in real space. If such an effect would be confirmed, it would mean that the important parameter to tune the absorption edge is the thickness of the isolated object, and not the distance between objects, even if neighbouring objects must be present.

Conclusion

In this chapter, we have tested a procedure to calculate the absorption spectrum of an interacting slab.

The idea is to consider that the standard TDDFT contains the correct information from the point of view of the electronic interactions, due to the Coulomb potential between the replicas arising from the periodic repetition of the supercell. However this information is corrupted by the vacuum: the resulting spectrum can be described via an Effective Medium Theory with vacuum, which gives an unphysical behavior of the absorption spectrum with vacuum.

To extract the absorption spectrum of the so-called interacting slab, cured from the spurious effect of vacuum, I have inverted the equations at the origin of the EMT, and named it Inverse Effective Medium Theory.

For hBN slabs, the results are very encouraging, since the IEMT spectra, for all the supercells considered here, behave as one could expect from the absorption of the interacting slab, somewhere in-between the bulk one and the isolated slab. A more surprising result is that when the spectra are well converged (LF along the z-direction), and the slabs isolated at the level of the ground state, the IEMT gives results independent of vacuum. The only case where the absorption of a slab still exhibits a vacuum effect corresponds to a very small amount of vacuum, and the slabs reveal to interact at the DFT level. The TDDFT does not introduce further interactions. One can also notice that for the largest slab (8-layer hBN), we almost recover the absorption of the bulk.

For silicon slabs, the results are more puzzling. I confirm that, as soon as the vacuum is large enough to consider that the slab are isolated at the DFT level, the IEMT spectra are independent of vacuum. The main surprising effect is that the absorption of the interacting slab occurs at much smaller energy than for the bulk one. This can be explained by the surface states, which are located in the gap of the surface projected bulk states (and which were absent in the case of hBN). Nevertheless the difference between the isolated slab and the interacting one needs further investigation. Even if the interaction between slabs arises from the surface states, the effect seems to be huge, and only a real space calculation with two interacting slabs could be a good fingerprint of the validity of the IEMT formalism.

Chapter 8

Application for the photovoltaics

InP

Another prospective material for photovoltaics is Indium phosphide (InP). InP solar cells, with a bandgap of 1.3 -1.4 eV [START_REF] Vurgaftman | [END_REF][START_REF]Landolt-Börnstein Semiconductors[END_REF][START_REF] Lautenschlager | [END_REF][70] and high absorption levels, have a strong potential of reaching the Shockley-Queisser theoretical efficiency level of approximately 33% [71].

It is also expected that coupling this intrinsic properties to specificity of nanotechnologies, like nanowire arrays [72], much better performance could be achieved.

Before reaching the slabs systems of these materials of photovoltaic interest, I studied the bulk counterpart.

InP bulk

The unit cell consist of two atoms: indium (In) with an atomic number 49 and phosphorus with an atomic number 15. The crystalline InP has a face-centered crystal structure (Fig. The electronic configuration of P is [N e] 3s 2 3p 3 . The pseudo-potential used was 5-P.LDA.fhi.

The electronic structure of In is [Kr] 4d 10 5s 2 5p 1 . An important point for these kinds of atom is the presence of the semi-core 4d electrons [75]. In order to test this effect, I first calculated the band structure at the local density approximation (LDA) level, using two different pseudo-potentials without 4d states (49-In.LDA.fhi), and with 4d electrons (49-In-4d.LDA.fhi). The results is shown on Fig. 8.2. The band gap is larger of 0.1 eV without the 4d states. The valence bands are little affected by the presence on 4d electrons (except the presence of the 4d states), the conduction band with 4d states are attracted to lower energy in a non rigid way. It is particularly true for bands above 15 eV around X point. Nevertheless, the part of the spectrum over interest is below 15 eV and the energy level without 4d states can be considered as a rather good approximation, for the case where the system is large, as it will be the case for the slab of InP. For the case of the bulk, where the calculations with the 4d states are tractable, I kept the pseudo-potential with 4d electrons.

Quasiparticles corrections for bulk InP

The well-known underestimation of the band gap in LDA gives a value of 0.361 eV, as compared to the experimental value of 1.3 -1.4 eV. In order to correct the band structure, I calculated the quasiparticle corrections, within GW approximation using Abinit code [44]. The GW value obtained is 0.8 eV, which is still underestimated. I decided to used the formalism of hybrid functionals, to test if it could provide a ground state with a better value for band gap: I chose Heyd-Scuseria-Ernzerhof (HSE) exchange-correlation functional (also using Abinit). The resulting band structure is presented with dots on Fig. 8.3, in dots, superposed to the solid lines of the LDA band structure. The value of the band gap is within this approximation 1.18 eV. It is much closer to the experimental band gap (1.3 -1.4 eV [START_REF] Vurgaftman | [END_REF][START_REF]Landolt-Börnstein Semiconductors[END_REF][START_REF] Lautenschlager | [END_REF][70]) than the GW corrections (0.8 eV), but it is still a bit smaller. GW correction on top of HSE correction would have probably allow us to reach a better value, but it has not been done.

TDDFT-RPA absorption for bulk InP

I calculated the absorption spetrum using Time Dependent Density functional Theory (TDDFT), within the random phase approximation (RPA), i. e. without exchange and correlation kernel to include excitonic effects. The results, both for LDA (black) and using HSE corrections (green) for the energies are presented in Fig. 8.4. The HSE spectrum is moved to higher energy, as expected from the increase of the band gap. It is not an exact shift as compared to the LDA one, but the difference is quite small, and in the perspective of larger calculations, the HSE corrections can be rather good mimicked by a rigid shift. I also add the experimental spectrum (red) from ref [76]). As it can be seen, the agreement between calculations and experiment is quite poor. The calculated spectra are badly converged with the k-point mesh. The grid used is the 8x8x8 with the 4 standard shifts of the Monkorst-Pack scheme. This density of points was correct to reach the convergence of HSE correction but it is to small to obtain a spectrum of good quality. Nevertheless, even if the spurious features of the calculated spectra does not allow an easy comparison, one notices that the absorption edge given by the HSE spectrum is in better agreement with the measurement. One notes that the TDDFT+HSE curve has a spectral weight at too large energy: this is a signature of excitonic effects, which are not taken into account at the TDDFT-RPA level. [76] is in red, experiment [70] is blue.

Excitonic effects for bulk InP

I calculated the excitonic effects using the Bethe-Salpeter equation (BSE). The calculations have been done using EXC code [START_REF]EXC Code[END_REF]. The results are presented on the Fig. 8.5. The experimental curve is in red for comparison. The calculated spectra are still done on the 8x8x8 -4 shifts Monkhorst-Pack mesh, for which I calculated the quasiparticle corrections. The screening has been calculated with LDA eigenvalues. Then I used different quasiparticle energy corrections for the spectrum.

In green, I plotted the spectrum calculated with BSE and using the GW corrections. The blue curve correspond to the BSE with the HSE corrections and the black curve with a scissor (SO) of 1 eV applied on the LDA energies. As expected from the value of the band gap, the edge of the spectrum with GW correction (green -E g = 0.8 eV) is at too low energy, the HSE spectrum (blue -E g = 1.18 eV) is in better agreement, and the scissor spectrum (black -E g = 1.36 eV) is even better. Concerning the overall spectrum, the two main peaks of the HSE and SO spectra are in good agreement with experiment, both from the point of view of the energy position and relative spectral weight. The Bethe-Salpeter equation allows to reproduce the experimental features, and thus correctly accounts for the excitonic effects. In InP, we do not observed a bound exciton, but a excitonic effect located in the continuum. The spectrum obtained with the SO (black) is very similar to the one obtained with HSE correction (blue), meaning that the slight dispersion of the bands that has been evidence in Fig. 8.3), is negligible, and we can safely use the SO of 1 eV to mimick quasiparticule corrections.

Nevertheless, the k-point mesh is not dense enough the give a satisfying spectrum, and we need to increase it a lot, which prevent to calculate brute force the quasiparticle corrections and the screening. Thus, I tested the TDDFT with the long-range kernel f xc = α/q 2 to account for the excitonic effects, when excitons are in the continuum, which is the case for InP. This framework is much less demanding, and I can increase the k-point mesh. Fig. 8.6 shows the result on a Monkhorst-Pack grid of 32x32x32 with the standard 4 shifts. The red curve is the experiment, the black curve is the TDDFT-RPA, and the blue line is the result of TDDFT, with the SO = 1 eV and the long range kernel with α = -0.4. The overall agreement is correct, even if it is not as good as the BSE one. The first peak at 3.2 eV is correctly located, the peak at 5 eV is a bit too high in energy, and finally the transfer of the spectral weight from the high energy peak to the low one is a bit to large.

InP slab

I modeled a slab of InP of thickness 22 Å. This choice was done to be in the same range of value than for silicon slabs, and of the range of experimentally realized setups. Nevertheless, we are still one order of magnitude smaller than the samples proposed in [72]. When we "cut" the bulk along the (1 0 0) direction, the obtained structure is metallic. In order to obtain a semi-conducting compound, with P atoms pending at the external of the surfaces, I needed to model a (2x2) reconstruction surface, which contains 60 atoms. The crystallographic structure of the corresponding slab is shown of Fig. Absorption spectra have been calculated within TDDFT at the RPA level. Figure 8.9 shows the absorption spectra for in-and out-of-plane components with the slab potential, compared with IEMT and the regular TDDFT. As expected from the very small band gap, the spectrum for in-plane component has an onset very close to 0 eV, spectra for interacting and isolated slabs are equivalent. The absorption spectrum for the out-of-plane component of the isolated slab is peaked at 12 eV, the IEMT spectra is shifted towards the low-energy region. Indeed, due the very large parameters required to converge the spectrum, a careful check of the out-of-plane must be drawn. Moreover, we know from the bulk calculation that the DFT-LDA band gap is largely overestimated, but HSE corrections are out of reach of a so huge system, and since we expect that the screening effect will be different in the slab from the one of the bulk, we don't know which value for the scissor could be chosen, and for this reason we do not further investigate this system.

InSe

InSe belongs to the familly of van der Waals (vdW) layered crystals. Their electronic properties depend on the composition, thickness and stacking of the component layers. As one of the typical III -VI semiconductors, InSe has attracted much attention due to its outstanding electronic properties, attractive quantum physics, and large photoresponse, suggesting that it is an excellent material for thin film optoelectronic applications [START_REF] Lei | [END_REF][79][80], and even solar cells [81].

The low interaction between layers allows an easy exfoliation, and the realization of sample of varying thickness, up to the monolayer. Due to the strong quantum confinement effect, few-layer InSe samples exhibit a layer-dependent band gap, spanning the visible and near infrared regions. It evolves from a direct band gap in the bulk (E g = 1.2eV [82]) to an indirect band-gap for few-layer stacking [83][84][85].

InSe bulk

The three-dimentional -polytype of InSe (space group D 1 3h in Schoenflies notation) has been build using the lattice parameters a = 3. The band gap was E g =0.38 eV, since the experimental band gap is 1.2 -1.3 eV [82]. I calculated the quasiparticle corrections using HSE hybrid functional, and the corrected gap is E g = 1.22 eV (Fig. 8.11 -dots), in good agreement with experiment. The spectra have been calculated within TDDFT, with LDA energies (blue) and with the HSE quasiparticle corrections for the energies (green) (Fig. 8.12a and 8.12b for inand out-of-plane components). The calculations have been done on the mesh 24x24x6 with 3989 plane waves, the number of bands is equal to 100 and the number of G-vectors along z-component is 17.

Since InSe is a layered material, it is anisotropic and the absorption spectrum is different for an excitation parallel (x) or perpendicular (z) to the plane of the layers. In both cases, the spectra with quasiparticle energies appear as mainly shifted [green (HSE) to compare with blue (TDDFT-LDA)]. To evaluate to which extend the quasiparticle corrections can be mimicked by a scissor operator (of 0.84 eV, if we target the HSE gap, of 1 eV to reach the experimental one), I calculated the spectra with SO = 0.84 eV (magenta curves). The agreement with the HSE curves very good in the energy region 0 -6 eV, then HSE spectra slightly shift to higher energy. The experimental spectrum taken from Ref. [START_REF] Piacentini | Electronic properties of the III-VI layer compounds GaS, GaSe and InSe[END_REF] is also plotted on Fig. 8.12a (red). Since BSE calculations is very demanding, I have used TDDFT with the long-range kernel to simulated the experiment. Calculations with α = -0.4 is plotted in black on Fig. 8.12a. The pre-edge at 2 eV is in good agreement with the measure. The spectral weight between 4-6 eV is also well reproduced, then the structure above 6 eV is not well located. Nevertheless, the narrow peak at 3 eV is missing: it is probably a bound exciton, which can only be captured within Bethe-Salpeter equation formalism.

InSe slab

Due to the size of the calculations, I decided to concentrate on a quite thin layer: a bilayer of InSe (8 atoms -16.64 Å). This system is of particular interest since is has been shown that in InSe, the gap evolves from a direct band gap at Γ for the bulk to a indirect band gap for few layers slabs [83][84][85]: it is at the origin of a huge activity around this system, as in general around the van der Waals layered crytals.

The band structure at the LDA level is presented on figure 8.13. As already been observed [83][84][85], the band gap evolves from an direct one in the bulk (at Γ) to an indirect one. The value of the direct gap for the 2 layer-InSe is 1.18 eV, since the indirect one is 1.10 eV: it arises between Γ for the conduction band and a point located between Γ and K or Γ and M for the valence band; this is in agreement with the "mexican hat" shape of the last valence band around Γ. These gaps are much larger than the one of the bulk (0.38 eV), at the same level of approximation (LDA).

Since we know from the bulk counterpart that the LDA gap is underestimated, I calculated the HSE corrections for InSe slab [88]: the resulting band structure is plotted with dots on Fig. 8.13. The spectra, within TDDFT, with and without HSE corrections for the energies are presented on Fig. 8.14a and 8.14b. The calculations are done using Selected-G method with slab potential, the IEMT has also been applied. Again, the spectra are mainly shifted by adding HSE corrections. For the in-plane component, the isolated slab has almost the same spectrum as the one with IEMT. For out-ofplane component, the IEMT spectrum has additional peak around 4 eV. These spectra are not the reliable absorption spectra. Based on the experience gained from the bulk calculation, we know that we need to use the Bethe-Salpeter equation. For this we need to calculate the self-energy, which contained screened Coulomb potential. Since we have evidenced that the inverse dielectric matrix was strongly dependent from the vacuum introduced in the supercell, we know that one must calculate the screening using the slab potential.

Conclusion

In this chapter, I have studied slab systems of materials with potential applications for photovoltaics. To defined which level of approximation, or within which formalism, I had to do the calculations, I have first studied the bulk counterpart.

I have evidenced that the local density approximation for the ground state was not a good starting point. I have corrected it by using an hybrid functional (HSE), and evaluate to which extend it could be mimicked by a scissor operator.

The calculation for InP bulk spectrum was done within Bethe-Salpeter equation formalism to correctly account for excitonic effects. Since it is a very demanding calculation, I have also tried the TDDFT formalism with the long-range kernel. The result is reasonably good to provide an absorption spectrum for a slab system, in a material where there is no bound exciton. This is promising for slabs of InP, even if due to the thickness of the slab I have modeled, the calculation was out of reach.

The calculation for InSe bulk was done within TDDFT formalism, but the comparison with experimental spectrum seems to indicated the presence of a bound exciton, which requires the use of the very computationaly demanding BSE calculation. The InSe slab of 2 layers has been studied. The HSE corrections have been calculated. The perspective is now to calculate the screening with the slab potential, in order to calculate the self-energy for the Bethe-Salpeter calculation.

Conclusion

The purpose of the thesis was to simulate the absorption spectrum of meta-materials for photovoltaic applications. By meta-material, we mean an assembly of nanometric size objects at mesoscopic distance. The underlying idea is that by adjusting the size of the nano-object and the geometric arrangement, one could tune the absorption edge. These quantities are calculated using ab initio methods, which are among the state-of-the art in the domain of theoretical spectroscopy. Since these formalisms can only deal with a small number of atoms, we cannot use brut force calculations.

The proposed strategy was to calculate the absorption for an isolated object, which in my case was the thin slab. The second step was to calculate the absorption for the interacting slabs, at nanometric distance in order to be tractable within a TDDFT calculation and to evaluate the influence of the interaction between nano-objects on the absorption edge. These results would allow a further analytic modelisation of these interactions to extend the distance between objects. Different compounds have been studied, depending on the steps of the methology. Graphene or few layers graphene slabs have been used to validate the slab potential, by comparison with EELS experiments. Few layers hBN slabs have been used to study the interacting slabs. This system was chosen because the slabs are semi-conducting, without the need of surface reconstruction, and also because it has a large gap, preventing numerical instabilities at low energy, when both the real and the imaginary parts of the macroscopic dielectric function are required. Silicon slabs were used for their application to photovaltaics, and because we already have a large experience with these slabs and their surface properties. Finally, I studied InP, which is used for photovaltaics in the bulk phase, but is also expected to be a good candidate at the nanometric scale. The last system considered is InSe. It is presently source of an important activity in the context of two-dimensional Van der Waal crystals, and in particular due to its property of having a cross-over from a direct to an indirect band gap when evolving from bulk to thin slabs.

The first four chapters were dedicated to the general description of the theoretical framework used in the thesis.

The ground states properties have been calculated within density functional theory (DFT). In all of the cases dedicated to development, technical tests and feasability, I used the local density approximation (LDA). For some cases, in particular the applications to photovaltaics, I computed the quasiparticule corrections using the GW approximation, or the hybrid functional formalism of Heyd-Scuseria-Ernzerhof.

105

To describe electronic excitations, I used Time-Dependent Density Functional Theory (TDDFT), which is much less computationaly demanding than the Bethe-Salpeter equation (BSE): all the developments have been done using this formalism. Nevertheless, to account for excitonic effects, which can be important in absorption spectra, I also used the code which solves the Bethe-Salpeter equation. It has allowed me to evaluate to which extent the use of the long range kernel in TDDFT could give a correct description of the excitonic effects, at the level of bulk material, before tackle the slab configuration.

In chapter 4, I focussed on the description of isolated objects. In particular, I have summarized the results obtained during the PhD of Nicolas Tancogne-Dejean [START_REF] Tancogne-Dejean | Ab initio description of second-harmonic generation from crystal surfaces[END_REF] in the Theoretical Spectroscopy group of the LSI. In this work, it has been evidenced that, in the framework of 3D periodic codes, the use of a supercell to describe isolated objects gives some non-physical spectra. In particular, when dealing with a surface, the out-of-plane component of the absorption spectrum was affected by local fields with a spurious contribution from the vacuum. To cure the problem, a new method has been proposed, called "Selected-G": the microscopic dielectric matrix is evaluated using a set of reciprocal lattice vectors restricted to the ones of defined by the thickness of the matter, and not the one of the supercell. The Dyson equation has been accordingly modified, with a non-diagonal expression for the reciprocal space Coulomb potential [START_REF] Tancogne-Dejean | [END_REF]. In the limit of a thick slab, which corresponds to the surface case, it reduces to the standard expression of the 3D Coulomb potential, and it has been successfully applied to surfaces [42]. The 2D expression of the Coulomb potential, dedicated to isolated slab, and for this reason called slab potential has not been used. This was the initial framework for my PhD thesis.

In chapter 5, to validate the use of the slab potential, I have studied Electron Energy Loss spectra on few graphene layers, where experimental data were available. I have shown that for EELS, the spurious effect of vacuum, introduced in the supercell to isolate the layers, was important also at the level of in-plane components. I have demontrated that the use of the slab potential allows us to obtain spectra independant of vacuum, and in very good agreement with the experiments on graphene and stacking of few layers graphene of T. Eberlein et al [45].

Since this slab potential acts as a cutoff procedure, I have analysed in chapter 6 the effect of the vacuum with some other cutoff methods which have already been proposed. I have shown that, since all these methods proposed to cut the Coulomb interaction somewhere in the middle of the supercell, they are always dependent of the size of the chosen supercell.

Then, using the slab potential and the Selected-G method, I have studied the in-plane dispersion of the plasmon of a single graphene layer, and compared it with an other similar study, which has developed a similar slab potential, but kept the full basis set for the reciprocal lattice vectors [57]. This work enters in the problematics of the nature of electronic excitations in graphene. My results showed that several ranges can be defined in terms of momentum transfer q, shedding some light on the origin of interband transitions for collective excitations.

One notes that the EEL spectra calculated within the formalism proposed in this thesis is different from the spectra obtained from the corresponding surface. This demonstrates the necessity to use the slab potential to properly describe EEL spectra of isolated slabs.

After having demonstrated the range of application of the slab potential to calculate electronic excitations in an isolated slab, chapter 7 is dedicated to the study of the absorption of interacting slabs. To benchmark the work, I have, in a first step, studied slabs interacting at nanometric distance, to allow the simulation with brut force calculations.

The interactions between slabs must be contained in the spectra resulting from the standard TDDFT, where the inclusion of local fields accounts for the microscopic fluctuations of the induced electronic density. Nevertheless, the spectrum resulting form the standard TDDFT suffers from the vacuum problem. It has been shown during the thesis of Nicolas Tancogne-Dejean [START_REF] Tancogne-Dejean | Ab initio description of second-harmonic generation from crystal surfaces[END_REF] that the spectra resulting from the standard TDDFT was actually equivalent to an effective medium theory (EMT) with vacuum. Thus, I have inverted the equations leading to EMT for small and large local fields. This Inverse Effective Medium Theory (IEMT) allows to extract the absorption of the so-called interacting slab. The validity of the procedure has been evidenced on the absorption spectrum of hBN. For in-plane components, where the local field effects are small, it is equivalent to a scaling factor. For the out-of-plane component, the spectrum is strongly modified. The out-of-plane absorption spectra obtained with IEMT are different from the bulk and the isolated slab ones, but nevertheless located in the same energy range.

I have evidenced that, for a given quantity of matter, the IEMT applied to the different supercells (containing different amount of vacuum), give the same spectra, providing the local fields perpendicular to the slab are well converged. The only case for which the spectrum is different is the one where the vacuum introduced in the supercell is too small to describe an isolated slab in the ground state. When the slabs are isolated in the ground state, the TDDFT step, with the appearance of the induced electronic density, sligthly modifies the spectrum, but in a way which seems to be independent from the distance between slabs. The spectra are dominated by the local fields at the surface of the slab.

The IEMT procedure has been applied to different slabs of silicon. The out-of-plane absorption is at much smaller energy than the isolated slab, and the bulk one. This can be explained by the presence of surface states, arising from the surface reconstruction of the silicon slab. Indeed, it is situated in the same energy range as for the in-plane component, and the ratio between the contribution of the surface states and of the bulk states evolves with the thickness of the slab. It results from the different results that the most important parameter to tune the gap seems to be the thickness of the slab, and not the distance between the slabs.

Nevertheless, the difference between the isolated slab and the so-called interacting one must be further investigated. The resolution of the Dyson equation in real space, with two interacting slabs, could be a first step to confirm these conclusions. These effects must also be studied in the light of excitonic effects, and the modification of the screening as a function of the thickness of the slab.

Finally, in the last part (chapter 8) I studied material commonly used, or expected, for photovoltaic applications: InP and InSe. Before tackling the problem of the slab systems, I have studied the bulk counterpart to evaluate which level of approximation would be required. As these materials are dedicated to photovaltaics, due to their excitonic properties, we should use the state-of-the-art formalism to calculate excitonic spectra, the Bethe-Salpeter equation. ou4nlf.EEL : ε -1,N LF 00 (q; ω) = 1 + Ṽ00 (q) χNLF 00 (q; ω) with χNLF 00 (q; ω) = χ 0 00 (q; ω) + χ 0 00 (q; ω) Ṽ00 (q) χNLF 00 (q; ω) ou4lf.EEL : ε -1,LF 00 ) decreases when q sl decreases (so black and green lines become equivalent). As a consequence, the main effect The fact that χ00 is slightly different from χ 0 00 , and that this difference is independent of q comes from Ṽ G G for G = 0, since it χ00 is the solution of the Dyson equation (A.5):

(q; ω) = 1 + G Ṽ0 G(q) χ G0 (q; ω) with χ G G (q; ω) = χ 0 G G (q; ω) + G1 G2 χ 0 G G1 (q; ω) Ṽ G1 G2 (q) χ G2 G (q; ω) ou1lf.EEL : ε -1,LF 00 (q; ω) = 1 + Ṽ00 (q) χ00 (q; ω) with χ G G (q; ω) = χ 0 G G (q; ω) + G1 G2 χ 0 G G1 (q; ω) Ṽ G1 G2 (q) χ G2 G (q; ω) ou2lf.EEL : G =0 Ṽ0 G(q) χ G0 (q; ω) with χ G G (q; ω) = χ 0 G G (q; ω) + G1 G2 χ 0 G G1 (q; ω) Ṽ G1 G2 (q) χ G2 G (q; ω) 0 2 4
χ G G (q; ω) = χ 0 G G (q; ω) + G1 G2 χ 0 G G1 (q; ω) V G1 G2 (q) χ G2 G (q; ω)
and V defined as V0 G2 = 0, ∀ G2 and V G1 G2 = Ṽ G1 G2 for G1 = 0:

V = 0 0 Ṽ G1 0 Ṽ G1 G2
.

One can also imagine that the contribution from Ṽ G0 will be negligible, since this term should contain a q dependence. The reason why the green ( χLF 00 ) and black ( χ00 ) curves becomes identical when q -→ 0 comes from the fact that the terms in the Dyson equation (A.8) arising from the multiplication with Ṽ0 G, ∀ G tend to zero when q -→ 0. Indeed, Ṽ00 is proportional to 1/q and the χ 0 00 ∝ q2 ; Ṽ0 G ∝ 1/ G and χ 0 G0 ∝ q. For this reason, the χ00 obtained from the solution of the Dyson equation (A.8) gives the same result as the χ00 when q -→ 0. This also confirms that the terms coming from the multiplication with Ṽ G0 are negligible, since due to symmetry reasons Ṽ G0 = Ṽ0 G.

range of q

Figure A.9 presents the Im(ε LF M ) for different length of q vector (q sl ). The black lines are spectra for q x component, the red lines spectra for q y component, the green lines are spectra for q z component (extracted from linear combination of spectra calculated with q x and q x ± q z ), and in blue spectra for q z component (extracted from linear combination of spectra calculated with q y and q y ± q z ). q sl = 10 -2 a.u. ou4MDF (y) q sl = 10 -2 a.u. ou4MDF (z) q sl = 10 -2 a.u. ou4MDF (z) q sl = 10 -3 a.u. ou4MDF (x) q sl = 10 -3 a.u. ou4MDF (y) q sl = 10 -3 a.u. ou4MDF (z) q sl = 10 -3 a.u. ou4MDF (z) q sl = 10 -4 a.u. ou4MDF (x) q sl = 10 -4 a.u. ou4MDF (y) q sl = 10 -4 a.u. ou4MDF (z) q sl = 10 -4 a.u. ou4MDF (z) q sl = 10 -5 a.u. ou4MDF (x) q sl = 10 -5 a.u. ou4MDF (y) q sl = 10 -5 a.u. ou4MDF (z) q sl = 10 -5 a.u. ou4MDF (z)

Figure A.9: Im(ε LF M ) for different q sl (ou4lf.MDF): (black) spectra for q x component; (red) spectra for q y component; (green) spectra for q z component (extracted from linear combination of spectra calculated with q x and q x ±q z ), and (blue) spectra for q z component (extracted from linear combination of spectra calculated with q y and q y ± q z ).

Results are independent of |q|, which allows us to conclude that χ00 is proportional to q 2 . This is a bit different from the silicon case, where the q z term for q sl = 10 -2 a.u. was slightly different from spectra for smaller q sl values. The q z component extracted from the different linear combinations is, as expected, the same. This is an indication of the numerical stability of the procedure, which is particularly delicate for graphitic systems: the metallic point at K leads to numerical divergences. ) [lines without dots (ou4nlf.MDF)] compared with Im(ε LF M ) [lines with dots (ou4lf.MDF)] for q sl = 10-5a.u.: (black) spectra for q x component; (red) spectra for q y component; (green) spectra for q z component (extracted from linear combination of spectra calculated with q x and q x ± q z ), and (blue) spectra for q z component (extracted from linear combination of spectra calculated with q y and q y ± q z ).

Effects of Local Fields

Figure A.10 shows the influence of local fields (LF): LF effects are almost negligible for the region 0 -10 eV corresponding to the π plasmon, for both in-and out-of-plane components. They are more important for 10 -30 eV range. For in-plane components, the effect is a small blue shift, with a small reduction of amplitude. It is much larger for the out-of-plane component, in agreement with the discontinuous electronic density. q sl = 10 -3 a.u. ou4lf.EEL (x) q sl = 10 -4 a.u. ou4lf.EEL (x) q sl = 10 -5 a.u. ou4lf.EEL (x) q sl = 10 -2 a.u. ou4nlf.EEL (x) q sl = 10 -3 a.u. ou4nlf.EEL (x) q sl = 10 -4 a.u. ou4nlf.EEL (x) q sl = 10 q sl = 10 -3 a.u. ou4lf.EEL (x) q sl = 10 -4 a.u. ou4lf.EEL (x) q sl = 10 -5 a.u. ou4lf.EEL (x) q sl = 10 -2 a.u. ou4nlf.EEL (x) q sl = 10 -3 a.u. ou4nlf.EEL (x) q sl = 10 -4 a.u. ou4nlf.EEL (x) q sl = 10 -5 a.u. ou4nlf.EEL (x) * 3.9720368 10 -2 / q sl ) (ou4nlf.EEL) for different length of q vector (q sl ): (black) spectra with LF for q x component; (red) spectra without LF for q x component. Bottom: same spectra normalised with 2 |q|L mat z .

Effects of Local Fields

Local field effects are very small (Fig. A.13). They simply reduced the amplitude, and shift a little bit to higher energy the π + σ plasmon. The effect increases with increasing q sl .

To study the influence of the different contributions coming from the slab potential, I plot for each q sl , the different quantities ( (q; ω) = 1 + Ṽ00 (q) χ00 (q; ω) As already noticed for Si slabs, the influence of the complementary contribution (blue line) is decreasing when q sl is decreasing (so black and green lines become equivalent). For the 8-layer graphene, it is nevertheless almost negligible for q sl < 10 -2 . As a consequence, the main effect (difference between black and red lines) comes from the inclusion of terms of the slab potential among the head of the matrix potential in the matrix inversion. Nevertheless, in the optical limit, all these contributions become small, and it seems that spectra have the shape given by χ 0 00 (furthermore normalized by 4π q 2 or Ṽ00 ). To check this hypothesis, figure A.15 shows the different χ: χ 0 00 (red); the green line is χLF 00 , which enters into the EELS (first term); the black line is χ00 , which enters into the absorption. In fact, neither EELS ( χLF 00 ) nor absorption ( χ00 ) are equal to χ 0 00 . The reason has been explained for the case of silicon.

with χ G G (q; ω) = χ 0 G G (q; ω) + G1 G2 χ 0 G G1 (q; ω) Ṽ G1 G2 (q) χ G2 G (q; ω) ou2lf.EEL : G =0 Ṽ0 G(q) χ G0 (q; ω) with χ G G (q; ω) = χ 0 G G (q; ω) + G1 G2 χ 0 G G1 (q; ω) Ṽ G1 G2 (q) χ G2 G (q; ω) 0 5 10
The EEL spectrum is similar to the absorption for very small q sl . When q increases, the black curve does not change (it is still quite close to χ 0 00 ). The effect seen on the green curve comes from the Ṽ00 term. ) with local fields for the slab (ou4lf.EEL) and the surface (ou5lf.EEL) for different q sl : (black) spectra of slab for q x component; (red) spectra of slab for q y component, (green) spectra of surface for q x component; (blue) spectra of surface for q y component.

Comparison with surface

The EELS of slab is clearly different from the EELS of the surface, even for in-plane component (Fig. A.16). The collective oscillation of electrons at plasmon frequency is affect by the thickness of the slab. The plasmon of the surface is similar to the plasmon of the bulk.

A.2.3 Conclusion

The behavior of the slab potential has been studied in details, on two cases which present different characteristics. The Si slab is cut in a bulk with strong covalent bonds. The graphene slab is part of a bulk with layers slightly linked with van der Waals interactions.

The absorption spectra of the slabs are similar to the surfaces, so the slab potential has almost no influence. For the absorption, the resolution of the Dyson equation with the slab potential does not introduce a dependence with the length of q. The local fields have small effects in-plane, mainly a reduction of the amplitude. The effect is larger for q z component. It is not so spectacular for 8-layers graphene than it is for silicon slab.

For the EELS, the spectra of the slabs are different form the surface ones: the slab potential has a strong influence on the EELS. The collective excitations are sensitive to the thickness of the slab, even for in-plane excitations. The EEL spectrum is dominated by the first term of the summation (Eq. 5.5), the remaining terms increase with increasing q. The resolution of the Dyson equation with the slab potential introduces a large dependence with the length of q, essentially contained in the head of the slab potential matrix. When reducing the length of q, the EELS for 8-layer graphene slab reached quite rapidly the shape of the vanishing q spectrum. It takes longer for the silicon slab.

The EEL spectrum at vanishing q recovers the spectral shape of the absorption, even if it does not result from the standard relation ε M (ω) = lim q-→0 1 ε -1 00 (q; ω)

(A.15)
In the standard procedure, the EELS recovers the absorption at vanishing q when increasing the vacuum, due to Effective Medium theory with vacuum. Indeed,

ε -1 00 = 1 ε M leads to -Im(ε -1 00 ) = Im( 1 ε M ) = Im(ε M ) Re(ε M ) 2 + Im(ε M ) 2
. When the vacuum increases Re(ε M ) -→ 1 and Im(ε M ) -→ 1, so

-Im(ε -1 00 ) = Im(ε M )
Here it is not the case, but we have the same conclusion. For slab systems, the reason why EELS has the same spectral shape as the absorption for vanishing q is due to the dependence in q of the slab potential.

with the standard relation (Eq. 3.20): ε M (q) = 1 ε -1 00 (q) The parameter used in the calculations are summarized in table B The spectra are exactly the same for the three lengths of q vectors (smallest and largest taken into account). For vanishing q, the spectra are equivalent to the no local fields one. Increasing the length of q, local field effects modify the spectra, but even for the largest q considered, we can converge the spectrum with much smaller npwmat_z. We do not understand why Novko and al. [57] need a so huge value. To be comparable with the length of the G z vectors that I have used to reach convergence with the Select-G method, I will further use npwmat_z = 2 * (3 * ratio)+1 in my calculations: it corresponds to npwmat_ z = 25 in the case of a system with ratio = 4. Since npwmat_xy = 19 allows the convergence when including also in-plane, I used this value. For the case of only out-of-plane LF, to be closer to the order of magnitude of the maximum G z vector proposed by Novko et al. [57], I used (npwmat_xy = 1; npwmat_z = 2 * (7 * ratio) + 1) in my calculations which corresponds to (npwmat_xy = 1; npwmat_z = 57) in the case of a system with ratio = 4]. For q = 10 -5 Bohr -1 , the local fied effects are almost negligible. Their inclusion only slightly reduces the amplitude. The two way of including them [only out-of-plane with (mat_xy = 1; npwmat_ z) in-and out-of-plane with (mat_xy = 19; npwmat_z)] does not give exactly the same results. As already mentioned, the only out-of-plane LF (pale lines) is equivalent to the spectra without LF. The effect seems to be independent of the vacuum included in the supercell. One also notice that EEL and absorption spectra are similar. For q = 0.078 Bohr -1 , the effect of local fieds is also to reduce amplitude. For Im(ε M ), one also observes a shift to higher energy. Larger the vacuum, more important this behavior. This effect seems to be less visible for EELS. The two way of including them (only out-of-plane with npwmat_xy = 1 and npwmat_z or in-and out-of-plane with npwmat_xy = 19 and npwmat_z) give exactly the same results for the π plasmon energy range. For the energy range of the π + σ plasmon, and more precisely between 11 and 15 eV, one observes a slight difference between the two LF inclusions: the spectrum with only out-of-plane LF are less depressed as the one with in-and out-of-plane LF. It seems that including also LF in-plane has a (small) influence. Nevertheless, the difference is rather small, and the two ways of including LF will not affect the behavior of the spectra with vacuum. EELS and absorption are different from each other. Title : Numerical design of metamaterials for photovoltaic applications.

Keywords : ab-initio, metamaterials, photovoltaics

Abstract : The purpose of this thesis was to simulate the absorption spectrum of metamaterials for photovoltaic applications. By metamaterial, one means an assembly of nanometric size objects at mesoscopic distance. The underlying idea is that by adjusting the size of the nano-object and the geometric arrangement, one could tune the absorption edge. The stateof-the art formalism, namely ab initio methods, was used. In the framework of periodic codes, one uses a supercell with vacuum to isolate the object. It was evidenced that the supercell formalism acts as an effective medium theory (EMT) with vacuum, with the spurious effect of having spectra dependent on the size of the supercell. A method called "Selected-G" was developed previously in the Theoretical Spectroscopy group at LSI to provide results independent of vacuum. It was successfully applied to the case of surfaces. The first part of the work was dedicated to the calculation of the absorption of an isolated slab. The Selected-G formalism were modified, using the slab potential. The Electron Energy Loss spectra for slabs of few graphene layers were simulated and successfully reproduced available experimental data. This has offered the possibility to study the nature of electronic excitations in graphene. The second step was dedicated to the study of the absorption of an array of interacting slabs. Based on the result that the standard formalism acts as an EMT, the equations of this theory were inverted to obtain the spectrum of the so-called interacting slab, cured from the vacuum problem. In the third part, materials for photovoltaic applications, like InP and InSe, have been studied. The electronic structure of the bulk counterparts, using hybrid HSE functional, have been calculated. The absorption spectra were simulated using the Bethe-Salpeter equation to account for the excitonic effects and using time-dependent density functional theory (TDDFT) with a long range kernel. A slab of InP, based on a (2x2) reconstructed surface, and a slab of 2 layers of InSe were modeled. The absorption spectra were calculated, within TDDFT for the isolated and interacting slabs.

Universit é Paris-Saclay

Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

Figure 1 . 1 :

 11 Figure 1.1: Schematical representation of independent particle approximation (left) and quasiparticle approximation (right).

Figure 1 . 2 :

 12 Figure 1.2: GW compare to LDA.

. 1 :Figure 2 . 1 :

 121 Figure 2.1: Khon-Sham schematic picture representation.

Figure 3 . 1 :

 31 Figure 3.1: Scheme of the perturbation of the ground state.

Figure 4 . 1 :

 41 Figure 4.1: Schematical representation of the supercell.

Figure 4 . 2 :

 42 Figure 4.2: Left: Slab of 32 atoms corresponding to the (2x1) surface reconstruction in a supercell. Right: electronic density obtained by DFT calculation. The red lines shows the limit of matter in the supercell.

Figure 4 . 3 :

 43 Figure 4.3: In-plane component(top) and out-of-plane component (bottom) for absorption of a silicon surface in IPA, for the three different supercells depicted in inset.

Figure 4 . 4 :

 44 Figure 4.4: Absorption spectra along the z-direction for a silicon surface in RPA for the three different supercells.

Figure 4 . 5 :

 45 Figure 4.5: Imaginary part of the out-of-plane component of ε M for the slab of silicon of thickness 21.72 Å introduced in supercells with different vacuum. Solid lines: results of regular TDDFT in supercell formalism. Dashed lines : EMT according to Eq. (4.6).

Figure 4 . 6 :

 46 Figure 4.6: Silicon surface absorption using Selected-G method.

Figure 5 . 1 :

 51 Figure 5.1: Comparison of Im(ε M ) for the 4_4 (black) and of 4_8 (red) supercells. Top: in-plane component. Bottom: out-of-plane component.

Figure 5 . 2 :

 52 Figure 5.2: The schematic representation of an EELS experiment.

Figure 5 . 3 :

 53 Figure 5.3: Schematic view of the supercells R2, R3 and R4 used in the calculations.

Figure 5 . 4 :

 54 Figure 5.4: EEL spectra of 1-layer graphene calculated with regular TDDFT for different supercells (the size of the supercell in Bohr is given in the legend). (top-left): in-plane component without LF, (bottom-left): in-plane component with LF, (top-right): out-ofplane component without LF, (bottom-right): in-plane component without LF.

  (a) Experimental EEL spectra (Eberlein et al[45]) of layers of graphene.(b) Numerical simulation (Eberlein et al[45]) of the EEL spectra of layered graphene. In-plane component.

Figure 5 . 5 :

 55 Figure 5.5: Eberlein EELS experiment (left) and numerical modeling (right).[45]

  EELS of 1, 2, 4 layers of graphene with equivalent vacuum.

Figure 5 . 6 :

 56 Figure 5.6: EEL spectra of regular calculation without slab potential and without selected-G method for in-plane component.

Figure 5 . 7 :

 57 Figure 5.7: EEL spectra of 4 layers graphene with different methods. (blue curves): regular TDDFT calculation within RPA; (red curves): Selected-G method for the surface case (no slab potential) and (green curves) correspond to the Selected-G with a slab potential. The black curve is the bulk graphite with AB stacking,

Figure 5 . 8 :

 58 Figure 5.8: EEL spectra of graphene for two supercells (ratio=4 and ratio = 5) using Selected-G and V slab with ξ = 1 ([-L mat z , 0] or [0, L mat z ] (matcenter=0: R4 solid red, R5: dashed blue) or ξ= (-1) n 1 +n 2 ([-L mat z /2, L mat z /2]) (matcenter=1: R4 solid green, R5: dashed black), q ΓM = 0.47 Bohr -1 .

Figure 5 . 9 :

 59 Figure 5.9: (left) EEL spectra for in-plane component of 1, 2, 3, 4, 8 layers of graphene calculated with selected-G method and slab potential. (right) Experimental data of Eberlein et al (see Fig. 5.5a) for comparison.

Figure 6 .

 6 Figure 6.1: EELS for one graphene layer in different supercells (R2, R3, R4, R5, R10) without local field effects (NLF-left) and with local field effects (LF-right) for q = 10 -5 Bohr -1 (top) and for q = 0.078 Bohr -1 (bottom). The calculations are done with the standard cutoff_2D.

Figure 6 . 2 :

 62 Figure 6.2: Real and Imaginary parts of ε M (top and middle respectively) and EELS (bottom) for one graphene layer in different supercells (R2, R3, R4, R5, R10) without local field effects (NLF-left) and with local field effects (LF-right) for q = 10 -5 Bohr -1 . The calculations are done with the cutoff_2D defined in Eq. 6.2 or 6.6.

6. 2 .Figure 6 . 3 :

 263 Figure 6.3: Real and Imaginary parts of ε M (top and middle respectively) and EELS (bottom) for one graphene layer in different supercells (R2, R3, R4, R5, R10) without local field effects (NLF-left) and with local field effects (LF-right) for q = 0.078 Bohr -1 (bottom). The calculations are done with the cutoff_2D. defined in Eq. 6.2 or 6.6.

  6

  .3: Im(ε N LF M ) (middle-left) and Im(ε -1,LF 00

Figure 6 . 4 :

 64 Figure 6.4: EELS for one graphene layer in different supercells (R3, R4, R5) without local field effects (NLF-left) and with local field effects (LF-right) for q = 10 -5 Bohr -1 (top) and for q = 0.078 Bohr -1 (bottom). The calculations are done with the cutoff proposed by Novko et al[57] using npwmat_xy=1 and npwmat_z=71 .

Figure 6 . 5 :Figure 6 . 6 :

 6566 Figure 6.5: Different response functions implied in EELS for q = 10 -5 Bohr -1 : (red) χ 0 00 , (blue) χ N LF 00 and (green) χ LF 00 , with χ LF 00 calculated with LF only out-of-plane (left) and with LF in-and out-of-plane (right).

Figure 6 . 8 :

 68 Figure 6.8: This figure is extracted from by Novko et al [57].

Figure 6 . 9 :

 69 Figure 6.9: Real (red lines) and Imaginary (blue lines) parts of ε M and EELS (green lines) for one graphene layer in R4 supercell (R4 = 25.173 Bohr) for different q ranging from 10 -5 , 0.039, 0.078 for the top and 0.157, 0.235, 0.313 Bohr -1 for the bottom plots from left to right. (The first graph is calculated for vanishing q). The calculations are done with slab potential and Selected-G method with (npwmat_xy = 19; npwmat_z = 7).

  B.2). The evolution of the spectral weight compares well with the experimental data of Wachsmuth et al[58]. = 25 for the non selected-G case, and allowed the convergence of the spectra (see Fig. B.1).

Figure 6 . 10 :

 610 Figure 6.10: Dispersion of front edge of π transitions for EEL and absorption spectra. Black squares are experimental data[59], red dots: dispersion of π plasmon, blue dots: corresponding structure of the absorption extracted from[57]. Cyan (magenta) plain circles + dashed lines corresponds to the π plasmon (absorption) dispersion taken from our calculations, done on a 100x100x1 mesh, with Selectd-G method (see Fig.6.9).

Figure 6 . 11 :

 611 Figure 6.11: Dispersion of front edge of π transitions for EEL spectra (cyan line) and fit with a q2 function (red line).

Figure 6 .

 6 Figure 6.12: (Top) comparison of the EELS spectrum ε -1,LF 00

Figure 6 . 13 :

 613 Figure 6.13: Dispersion of front edge of π transitions for EEL spectra (cyan line) and fit with a √ q function (red line) and a linear function (green dashed line). Black squares is a Liou et al experiment [59].

Figure 6 . 14 :

 614 Figure 6.14: Real (red lines) and Imaginary (blue lines) parts of ε M and EELS (green lines) for one graphene layer in R5 supercell for different q ranging from 0.016 up to 0.188 Bohr -1 by step of 0.016 Bohr -1 . The calculations are done on the (100x100x1) mesh, using slab potential and Selected-G method with (npwmat_xy = 19; npwmat_z = 7).

Figure 6 . 15 :

 615 Figure 6.15: Real (red lines) and Imaginary (blue lines) parts of ε M and EELS (green lines) for one graphene layer in R5 supercell for different q ranging from 0.2 up to 0.7 Bohr -1 . The calculations are done on the (40x40x1) mesh, using slab potential and Selected-G method with (npwmat_xy = 19; npwmat_z = 7).

-

Figure 6 . 16 :

 616 Figure 6.16: Real (red lines) and Imaginary (blue lines) parts of ε M and EELS (green lines) for one graphene layer in R5 supercell for different q ranging from 0.2 up to 0.7 Bohr -1 . The calculations are done on the (100x100x1) mesh, using slab potential and Selected-G method with (npwmat_xy = 19; npwmat_z = 7). The corresponding dashed lines are the spectra without LF.

Figure 7 . 1 :

 71 Figure 7.1: In-plane components of the dielectric tensor Im(ε M ), for the different hBN slabs (see each panel for the thickness). The blue lines are the results of standard TDDFT for different suprcells. The green dashed line is the isolated slab, and the red curve is the bulk.

Figure 7 .Figure 7 . 2 :

 772 Figure 7.2 shows the out-of-plane absorption spectra for the systems listed in table 7.1.

Figure 7 . 3 :

 73 Figure 7.3: out-of-plane component of Im(ε M ) resulting from IEMT on the 2-layer slab into R4 supercell (blue line) compared to the bulk spectrum (black line) and to the isolated slab (red line). k-point mesh: 40x40x1 for Isolated slab and IEMT ; 40x40x15 for the bulk.

Figure

  Figure 7.3 shows the out-of-plane component of Im(ε M ) resulting from IEMT applied to the 2-layer slab into R4 supercell (blue line). I compared it to the bulk spectrum (black line) and to the isolated 2-layer slab (red line). The pre-edge is exactly superimposed to the one of the isolated slab spectrum. Above 17 eV, the IEMT spectrum also reproduces the isolated slab, and the features are very similar to the one of the bulk. Between 8 and 17 eV, the amplitude of the peaks is different from both the bulk and the isolated slab. The peaks are located between these two spectra, in agreement with the idea that interacting slabs have an absorption which "feel" more electronic interactions via the Coulomb potential than the isolated slab, and less than the bulk material.Nevertheless, the results presented here are extremely sensitive to the number of G z vectors and a careful analysis of the convergence for the local fields has been performed.

Figure 7 . 4 :

 74 Figure 7.4: Absorption spectra along z-direction from IEMT between 4 and 13 eV of the 2-layer slab for various supercell size (R2, R3, R4, R5, R10) as a function of the p used to define max G z to include out-of-plane local fields: p = 3 (top-left); p = 4 (top-right); p = 5 (bottom-left); p = 6 (bottom-right). Mesh for the bulk is 40x40x15, for Isolated slab is 40x40x1, for IEMT 5x5x1.

Figure 7 . 5 :

 75 Figure 7.5: Absorption spectra along z-direction of the 2-layer slab for R2 supercell as a function of k-point meshes 15x15xN z for N z = 1 (solid cyan), 2 (dashed magenta), 4 (long dashed-dotted blue), 6 (short dashed-dotted green), 8 (solid red). Calculations are done within regular TDDFT without LF.

Figure 7 . 6 :

 76 Figure 7.6: IEMT applied to regular TDDFT absorption spectra along z-direction of the 2-layer slab for the different supercells. The R2 case has been calculated with a k-point mesh 20x20x4, the other supercells with a grid 20x20x1. The local field effects are included. npwmat_xy and npwmat_z parameters has been taken according to the table 7.2.

Figure 7 . 7 :

 77 Figure 7.7: Im(ε M ) for in-plane q for the different hBN slabs (top-left) 1-layer, (middle-left) 2-layer, (middle-right) 3-layer, (bottom-left) 4-layer (bottom-right) 8-layer. The solid red line corresponds to the bulk. The green dashed line is the isolated slab and the solid + dot blue line is the interacting slab.

Figure 7 . 8 :

 78 Figure 7.8: Im(ε M ) for out-of-plane q for the different hBN slabs (top-left) 1-layer, (middleleft) 2-layer, (middle-right) 3-layer, (bottom-left) 4-layer (bottom-right) 8-layer. The solid red line corresponds to the bulk. The green dashed line is the isolated slab and the solid + dot blue line is the interacting slab.

Figure 7 .

 7 Figure 7.8 presents out-of-plane absorption spectra for the different hBN slabs. In blue, I plotted the IEMT applied to the regular TDDFT calculation to mimick the interacting slab and I compared to the bulk (red) and isolated slab spectra (green).

Figure 7 . 9 :

 79 Figure 7.9: In-plane absorption of the silicon slab 4_4. Green: regular TDDFT, cyan: isolated slab, blue: IEMT, and red: bulk.

Figure 7 . 10 :

 710 Figure 7.10: Out-of-plane absorption of the silicon slab 4_4, (top) and 4_8 (bottom). Green: regular TDDFT, cyan: isolated slab, blue: IEMT, and red: bulk.

1

 1

Figure 7 . 11 :

 711 Figure 7.11: Surface states of the 4_4 relaxed slab (red) superposed to the surface projected bulk band structure (blue).

Figure 7 . 13 :

 713 Figure 7.13: Comparison of the out-of-plane absorption of the silicon slabs for isolated and interacting (IEMT) slabs of different thicknesses. "3" (16.29 Å): isolated (solid blue); interacting (dotted blue). "4" (21.72 Å): isolated (solid orange); interacting (dotted orange). "6" (32.59 Å): isolated (solid green); interacting (dotted green). Bulk (red).

  8.1 [73]).

Figure 8 . 1 :

 81 Figure 8.1: InP bulk crystallographic structure

  13 electrons, gap=0.361 NO 4d, 3 electrons, gap=0.465 InP with and without 4d -level

Figure 8 . 2 :

 82 Figure 8.2: InP bulk, with and without 4d level.

Figure 8 . 3 :

 83 Figure 8.3: Band structure of InP bulk. Solid lines: LDA (E LDA g = 0.36 eV); dots: HSE corrections (E LDA g

Figure 8 . 4 :

 84 Figure 8.4: InP absorption spectrum calculated within TDDFT-LDA (black), TDDFT-LDA+HSE (green). Experiment from[76] is in red, experiment[70] is blue.

Figure 8 . 5 :

 85 Figure 8.5: InP absorption spectrum calculated within BSE formalism. Green: GW corrections, blue: HSE correctiosn, black: scissor operator 1 eV. (For memory , TDDFT-LDA+HSE is in magenta). Experiment from experiment from [70] is in red.

Figure 8 . 6 :

 86 Figure 8.6: InP absorption: red: Experiment from [70]. Blue: TDDFT with SO=1 eV, and α = -0.4; Black: TDDFT-LDA + SO = 1 eV.

8 . 7 .

 87 The relaxation of the atomic position leads to a semi-conducting band structure (Fig.8.8) with E g = 0.38 eV.

Figure 8 . 7 :

 87 Figure 8.7: Crystallographic structure of InP (2x2) slab, side view.

Figure 8 . 8 :

 88 Figure 8.8: InP slab: LDA band structure (left). Brillouin zone (right) of unit cell and a path in a band structure: ΓXSGZ, the last Z-point is a way in out-of plane direction (0 0 1/2).

Figure 8 . 9 :

 89 Figure 8.9: InP slab absorption spectra: (top) x-component, (bottom) z-component. Parameters for the calculation are: npwwfn = 220 000, nbands 500, npwmat_xy = 5, npwmat_z = 15, mesh = 2x2x1.

  953 Å, d InIn = 2.741 Å, and d SeSe = 5.298 Å [86]. -InSe contains 8 atoms per unit cell. The crystal structure is shown on Fig. 8.10a and 8.10b.

  (a) InSe structure, side view (b) InSe structure, top view Following the same procedure as for InP, I have calculated the LDA band structure of the bulk InSe (Fig. 8.11 -solid lines). I used for In and Se norm-conserving Troullier-Martins pseudo-potentials generated by fhi98PP (49-In-4d.LDA.fhi and 34-Se.LDA.fhi for In and Se respectively). Each In atom participates with 13 valence electrons [4d 10 5s 2 5p 1 ], and each Se atom with 6 [4s 2 4p 4 ].

Figure 8 . 11 :

 811 Figure 8.11: InSe bulk: LDA and HSE band structure, E LDA g

  Absorption spectrum for InSe (x-component), in comparison with the experiment[START_REF] Piacentini | Electronic properties of the III-VI layer compounds GaS, GaSe and InSe[END_REF]. Absorption spectrum for InSe (z-component).

Figure 8 . 12 :

 812 Figure 8.12: Absorption spectra for InSe bulk for mesh 24x24x6.

4 (Figure 8 . 13 :

 4813 Figure 8.13: Band structure of 2 layer-InSe. LDA (lines) and HSE corrections (dots).

  (a) InSe 2 layers slab, x-component.

  (b) InSe 2 layers slab, z-component.

Figure 8 . 14 :

 814 Figure 8.14: Isolated bilayer slab of InSe done with mesh 24x24x6.

Figure A. 6 :

 6 Figure A.6: The different contributions which build the EELS for different q sl .

Figure A. 10 :

 10 Figure A.10: Im(ε N LF M

  .u. ou4lf.EEL (x)

  .u. ou4lf.EEL (x)

Figure A. 13 :

 13 Figure A.13: Top: -Im(ε -1,LF 00

  Fig. A.14):

Figure A. 14 :

 14 Figure A.14: The different contributions which build the EELS for different q sl .

Figure A. 15 :

 15 Figure A.15: Comparison of the different χ 00 for different q sl .

  (x) ou4lf.EEL (y) ou5lf.EEL (x) ou5lf.EEL (y) q sl = 10 -2 a.u.

Figure A. 16 :

 16 Figure A.16: Comparison of -Im(ε -1 00) with local fields for the slab (ou4lf.EEL) and the surface (ou5lf.EEL) for different q sl : (black) spectra of slab for q x component; (red) spectra of slab for q y component, (green) spectra of surface for q x component; (blue) spectra of surface for q y component.

For

  Select-G method, I used 19 in-plane vectors, corresponding to 3 closed shells (up to (2 0)), and 7 G z vectors, corresponding to G z = 0 ; G z = ±1 * ratio ; G z = ±2 * ratio ; G z = ±3 * ratio. The calculation done by D. Novko at al. [57] considered only LFE out-of-plane, with 71 G z to reach convergence. Associated to a supercell of heigth L SC z = 5a = 23.26 Bohr, it corresponds to a length of maximum G z = 9.5 Bohr -1 . In order to reproduce their results, I chose the system with L SC z = 4d 0 = 25.173 Bohr, which is quite close to the one of their supercell. I calculated with npwmat_xy = 1 and npwmat_z = 77, (corresponding to G z = 0 ; G z = ±1 ; G z = ±2....G z = ±38 giving 77 G z vectors with maximum G z = 9.5 Bohr -1 ). Since this value is quite large and prevent to introduce in-plane local fields, I have checked if it could be reduced and I have calculated (npwmat_xy = 1; npwmat_z = 71), (npwmat_xy = 1; npwmat_z = 49) and (npwmat_xy = 1; npwmat_z= 25). The corresponding lengths are summarized in Table B.2. The results are presented on fig B.1.

3 :Figure B. 1 :

 31 Figure B.1: left: EEL; right: Im(ε M) for q = 10 -5 Bohr -1 (top), q = 0.078 Bohr -1 (middle) and q = 0.627 Bohr -1 (bottom) for different out-of-plane (only) LF: npwmat_z = 77, npwmat_z = 71, npwmat_z = 49, npwmat_z = 25. Calculation has been done with cutoff proposed by Novko et al.[57].

Figure B. 2 :

 2 Figure B.2: left: EEL; right: Im(ε M ) for q = 10 -5 Bohr -1 (top), q = 0.078 Bohr -1 (middle) and q = 0.627 Bohr -1 (bottom) for different in-and out-of-plane LF: npwmat_z = 25, and npwmat_xy = 7, npwmat_xy = 19, npwmat_xy = 31. Calculation has been done with cutoff proposed by Novko et al.[57]

  The comparison of the spectra without LF, and with the two different inclusion of LF are presented on Fig.(B.3) for q = 10 -5 Bohr -1 and on Fig.(B.4) for q = 0.078 Bohr -1 .

Figure B. 3 :

 3 Figure B.3: Imaginary part of ε M (blue lines) and EELS (green lines) for one graphene layer in different supercells [R2 (top-left), R3 (top-right), R4 (bottom-left), R5 (bottomright)] for q = 10 -5 Bohr -1 : without LF (pale lines), with local field effects only along z (npwmat_xy = 1; npwmat_ z: lines with crosses) and with local field effects in-and out-ofplane (npwmat_xy = 19; npwmat_ z: bright lines). The different values of npwmat_ z are given in the legend: they depend on the vacuum. The calculations are done with the cutoff proposed by Novko et al. [57]. Re(ε M ) is equal to 1.

Figure B. 4 :

 4 Figure B.4: Real (red lines) and Imaginary (blue lines) parts of ε M and EELS (green lines) for one graphene layer in different supercells [R2 (top-left), R3 (top-right), R4 (bottom-left), R5 (bottom-right)] for q = 0.078 Bohr -1 : without LF (pale dashed lines), with local field effects only along z (npwmat_xy = 1; npwmat_ z: bright dashed lines with dots) and with local field effects in-and out-of-plane (npwmat_xy = 19; npwmat_ z: bright continuous lines). The different values of npwmat_ z are given in the legend: they depend on the vacuum. The calculations are done with the cutoff proposed by Novko et al.[57]

Titre:

  Design num érique de m étamat ériaux pour des applications photovoltaïques. Mots cl és : ab-initio, m étamat ériaux, photovoltaïque R ésum é : Le but de cette th èse était de simuler le spectre d'absorption de m étamat ériaux pour les applications photovoltaïques. Par m étamat ériaux, on entend une assembl ée d'objets de taille nanom étrique situ és à distance m ésoscopique. L'id ée sous-jacente est qu'en modifiant la taille du nano-objet et l'arrangement g éom étrique, on peut ajuster le seuil d'absorption. L' état de l'art du formalisme, c'est-à-dire des m éthodes ab initio, a ét é utilis é. Dans les codes p ériodiques, on utilise une supercellule avec du vide pour isoler l'objet. Il a ét é montr é que ce formalisme agit comme une th éorie de mat ériau moyen (EMT) avec le vide, avec l'effet erron é d'avoir des spectres d épendant de la taille de la supercellule. Une m éthode appel ée "Selected-G" a ét é d évelopp ée pr éc édemment dans le groupe de Spectroscopie Th éorique du LSI pour obtenir des r ésultats ind épendants du vide. Elle a ét é appliqu ée avec succ ès aux cas des surfaces. La premi ère étape de ce travail a ét é d édi ée au calcul de l'absorption d'une couche isol ée. Le formalisme Selected-G a ét é modifi é en utilisant le potentiel de slab. Les spectres de perte d' énergie d' électrons pour des empilements de quelques plans de graph ène ont ét é simul és, et reproduisent tr ès bien les donn ées exp érimentales. Cela a offert la possibilit é d' étudier la nature des excitations électroniques dans le graph ène.La seconde étape a ét é consacr ée à l' étude de l'absorption d'une assembl ée de couches en interaction. En s'appuyant sur le fait que le formalisme standard en supercellule agit comme une EMT avec du vide, les équations de cette th éorie ont ét é invers ées, pour obtenir l'absorption de la couche en interaction, en s'affranchissant du probl ème du vide. La troisi ème partie s'est int éress ée à l' étude de mat ériaux pour des applications photovoltaïques comme InP et InSe. La structure électronique des mat ériaux massifs a ét é calcul ée en utilisant la fonctionnelle hybride HSE. Le spectre d'absorption a ét é calcul é avec l' équation de Bethe-Salpeter, qui contient les effets excitoniques, et avec la th éorie de la fonctionnelle de la densit é d épendant du temps (TDDFT) avec le kernel à longue port ée. Une couche de InP utilisant une reconstruction de surface (2x2), et une couche de deux plans de InSe ont ét é mod élis ées. Les spectres d'absorption ont ét é calcul és en TDDFT pour les couches isol ées et en interaction.

  dr 2 dr 3 dt 2 dt 3 χ 0 (r, t; r 2 , t 2 ) v(r 2 , r 3 ) + f xc (r 2 , t 2 ; r 3 , t 3 ) χ(r 3 , t 3 ; r , t 3 )

3 , t 3 ; r , t ) + dr 4 dt 4 f xc (r 2 , t 2 ; r 4 , t 4 )χ(r 4 , t 4 ; r , t )   = χ 0 (r, t; r , t ) (3.13) + dr 2 dr 3 dt 2 dt 3 χ 0 (r, t; r 2 , t 2 ) v(r 2 , r 3 ) + f xc (r 2 , t 2 ; r 3 , t 3 ) χ(r 3 , t 3 ; r , t ) Finally, we will have the Dyson equation: χ(r, t; r , t ) = χ 0 (r, t; r , t ) (3.14) +

Table 5 .

 5 1: Parameters used to calculate the spectra for the 4_4 and 4_8

		bilayers of silicon (thickness L mat z	= 41.052 Bohr) has
	been introduced in a supercell with either L SC z	= 82.104 Bohr (ratio = 2, called 4_4), or
	L SC z	= 123.156 Bohr (ratio = 3, called 4_8).
	The spectra have been calculated with the parameters summarized in table 5.1, on a grid
	of k-point mesh 16x8x1.
			4_4	4_8
		npwwfn	12803 19209
		nbands	450	450
		npwmat_xy		25	25
		npwmat_z		33	33

Table 5 .

 5 .2. The converged parameters for the DP calculations are presented in Table5.3.

	Number of layers Ratio	L supercell L matter z	L supercell [Bohr]
	1	4		25.176
	2	4		50.352
	3	3		56.646
	4	3		75.528
	8	1		100.704
	Table 5.2: The parameters of the 1, 2, 3, 4 and 8 layers of Graphene for calculation with
	slab and selected G method.
		1l R4 2l R4 3l R3 4l R3 8l R2
	npwwfn	1299 2977 3289 4995 5479
	nbands	90	120	120	120	190
	npwmat_xy	19	19	19	19	19
	npwmat_z	25	45	37	45	49
	npwmat_z selG+slab 7	13	19	25	49

3: Parameters for convergence for DP, carbon layers.

Table 6 .

 6 .1).

	Thickness	R1	R2	R3	R4	R5	R10
	L SC z	(Bohr) 6.294 12.588 18.882 25.173 31.470 62.940

1: Heigth of the different supercells in Bohr.

The R1 length corresponds to half of the graphite AB c-unit cell, also called d 0 . The parameters used for convergence are the number of plane waves for the wavefunctions (npwwfn),

Table 6 . 2

 62

: Parameters used for convergence.

Table 7 .

 7 1: Summary of the differents systems of hBN studied.

	nb-layers R1 bulk	R2	R3 R4	R5	R10
	1	1933/80/25	37 2393/80/49	2885/80/61	5783/80/121
	2 1191/50/25 4689/100/49 73 4689/100/97 5783/100/121 11579/100/241
	3	73	109 145	8681/120/181
	4	4689/150/97 145 9391/150/193 11579/150/241
	8	9295/190/193

Table 7 . 2

 72

	nb-layers R1 bulk	R2	R4	R5	R10
	1	20x20x4 20x20x4 20x20x1 20x20x1
	2 40x40x15 20x20x4 20x20x1 20x20x1 20x20x1
	3			20x20x1
	4	20x20x1 20x20x1 20x20x1
	8	20x20x1

: npwwfn/nbands/npwmat_z parameters. npwmat_xy = 13 for all the cases. The grey color is a recommended value of npwmat_z for systems which have not been calculated.

Table 7 .

 7 3: Mesh for hBN calculations.

	Parameters	1l R5	2l R4	3l R3	4l R2	8l R2
	N k	40x40x1 40x40x1 40x40x1 40x40x1 40x40x1
	N b	80	80	120	150	190
	npwwfn	4001	4483	4599	4991	4989
	npwmat_z	25	25	37	25	49

Table 7 .

 7 Figure 7.1 shows the in-plane absorption spectra for the systems listed in table 7.1.

				1 layer
				bulk AA'
				regular R2
				regular R4
				regular R5
				regular R10
	Im(ε M ) ||			Isolated
	0	5	10	15	20	25	30
			Energy, [eV]
				2 layers			14			3 layers
				bulk AA'						bulk AA'
				regular R2 regular R4		12			regular R5 Isolated
				regular R5 regular R10		10
	Im(ε M ) ||			Isolated		Im(ε M ) ||	6 8
								4
								2
								0
	0	5	10	15	20	25	30	0	5	10	15	20	25	30
			Energy, [eV]					Energy, [eV]
				4 layers			12			8 layers
				bulk AA'						bulk AA'
				regular R2 regular R4		10			regular R2 Isolated
				regular R5
				Isolated			8
	Im(ε M ) ||						Im(ε M ) ||	6
								4
								2
								0
	0	5	10	15	20	25	30	0	5	10	15	20	25	30
			Energy, [eV]					Energy, [eV]

4: Parameters for hBN, selected G + V_slab potential. npwmat_xy = 13.

Table 7 .

 7

	)

5: (Top) npwmat_z (Bottom) corresponding values of max G z to achieved convergence for the 2-layer system as a function of the p parameter. The number in parenthesis is the so-called "equivalent" number of layers. npwmat_z is calculated according to npwmat_z = 2 • [Nb_eq_layer • p] + 1.

Table 7 .

 7 6: Definition of the slabs: first line: multiple of the unit block, second line: thickness (Å), third line: number of atoms.

	.6.
	N	3	4	6
	Slab's thickness (Å) 16.29 21.72 32.59
	Number of atoms	24	32	48

Table 7 .

 7 7: Summary of parameters used to calculate slabs of silicon. (nota bene: npwwfn is largely over converged. It is due to technical reasons).

	.7.
		3_3	4_4	4_8	6_6
	npwwfn	15237 12803 19209 12319
	nbands	350	450	450	650
	npwmat_xy	25	25	25	25
	npwmat_z (Selected-G)	25	33	33	49
	npwmat_z (IEMT)	49	65	97	97

Table B .

 B .1. 1: Parameters used for obtaining the spectra.

	Parameters	R2	R3	R4	R5	R10
	N k	40x40x1 40x40x1 40x40x1 40x40x1 40x40x1
	N b	100	100	100	100	100
	npwwfn	997	1489	1795	2487	3999
	npwmat_xy	19	19	19	19	19
	npwmat_z	13	19	25	31	61

  ) z 31 xy 19 LF EEL mat z 31 mat xy 19 LF Im (ε M ) mat z 71 mat xy 01 LF EEL mat z 71 mat xy 01 LF Im (ε M ) NLF EEL NLF

	0.0014				R2	0.0014				R3
	0012				Im (ε M ) mat z 13 mat xy 19 LF EEL mat z 13 mat xy 19 LF Im (ε M ) mat z 29 mat xy 01 LF EEL mat z 29 mat xy 01 LF Im (ε M ) NLF	0.001 0.0012				Im (ε M ) mat z 19 mat xy 19 LF EEL mat z 19 mat xy 19 LF Im (ε M ) mat z 43 mat xy 01 LF EEL mat z 43 mat xy 01 LF Im (ε M ) NLF
					EEL NLF	0.0008				EEL NLF
						0.0006
						0.0004
						0.0002
						0
	0	2	4	6	8 10 12 14 16 18 20 22 24	0	2	4	6	8 10 12 14 16 18 20 22 24
					Energy, [eV]					Energy, [eV]
	0.0014				R4	0.0014				R5
	0.001 0.0012				Im (ε M ) NLF EEL mat z 57 mat xy 01 LF Im (ε M ) mat z 57 mat xy 01 LF EEL mat z 25 mat xy 19 LF Im (ε M ) z 25 xy 19 LF	0.001 0.0012				Im (ε
	0.0008				EEL NLF	0.0008
	0.0006					0.0006
	0.0004					0.0004
	0.0002					0.0002
	0					0
	0	2	4	6	8 10 12 14 16 18 20 22 24	0	2	4	6	8 10 12 14 16 18 20 22 24
					Energy, [eV]					Energy, [eV]

M

  140APPENDIX B. CONVERGENCE OF LOCAL FIELDS WITH D. NOVKO ET AL. CUTOFF ) mat z 13 mat xy 19 LF Im (ε M ) mat z 13 mat xy 19 LF EEL mat z 13 mat xy 19 LF Re (ε M ) mat z 29 mat xy 01 LF Im (ε M ) mat z 29 mat xy 01 LF EEL mat z 29 mat xy 01 LF ) mat z 19 mat xy 19 LF Im (ε M ) mat z 19 mat xy 19 LF EEL mat z 19 mat xy 19 LF Re (ε M ) mat z 43 mat xy 01 LF Im (ε M ) mat z 43 mat xy 01 LF EEL mat z 43 mat xy 01 LF ) mat z 25 mat xy 19 LF Im (ε M ) mat z 25 mat xy 19 LF EEL mat z 25 mat xy 19 LF Re (ε M ) mat z 57 mat xy 01 LF Im (ε M ) mat z 57 mat xy 01 LF EEL mat z 57 mat xy 01 LF ) mat z 31 mat xy 19 LF Im (ε M ) mat z 31 mat xy 19 LF EEL mat z 31 mat xy 19 LF Re (ε M ) mat z 71 mat xy 01 LF Im (ε M ) mat z 71 mat xy 01 LF EEL mat z 71 mat xy 01 LF Re (ε M ) NLF Im (ε M ) NLF EEL NLF

	3 4 5 6 7				R2 Re (ε Re (ε M ) NLF Im (ε M ) NLF EEL NLF	3 4 5 6 7				R3 Re (ε Re (ε M ) NLF Im (ε M ) NLF EEL NLF
	2					2
	1					1
	0					0
	0	2	4	6	8 10 12 14 16 18 20 22 24	0	2	4	6	8 10 12 14 16 18 20 22 24
					Energy, [eV]					Energy, [eV]
	3 4 5 6				R4 Re (ε Re (ε M ) NLF EEL NLF Im (ε M ) NLF	3 4 5 6				Re (ε	R5
	2					2
	1					1
	0					0
	0	2	4	6	8 10 12 14 16 18 20 22 24	0	2	4	6	8 10 12 14 16 18 20 22 24
					Energy, [eV]					Energy, [eV]

M M M M

Further we will consider m e = e = = 1.5

We remind that 7 G z vectors correspond to G z = 0 ; G z = ±1 * ratio ; G z = ±2 * ratio ; G z = ±3 * ratio. Since we are in the case ratio = 4, is gives 12 for the maximum G z , which is equivalent to consider npwmat_z

This was one of the motivation of the study of the absorption of hBN: (i) the slabs are semi-conducting, which avoid the necessity to relax the atomic positions and consequently to introduce surface states. (ii) hBN has a large gap, which prevents from numerical instabilities at low energy. This system has allows us to validate the fact that IEMT can be applied to obtain the correct absorption spectrum.

Acknowledgments

I would like to thank my supervisors Christine Giorgetti and Valérie Veniard for all the guiding during the thesis, for sharing their knowledge and being ready to answer all the difficult questions and not only scientific ones. I am grateful to Lucia Reining, Francesco Sottile, Matteo Gatti, Andrea Cucca and Igor Reshetnyak for productive discussions and strong assistance. I would like to thank all the ETSF group and all the PhD students for creating very friendly and productive environment during all my thesis. And finally I would like to thank my wife Uliana and my parents for supporting me though all the thesis.

However, it is extremely demanding computationaly, and the TDDFT with a kernel is a much lighter alternative which can give very good results. I calculated quasiparticle energy corrections, within the GW formalism, and also using HSE hybrid fonctional. Absorption spectra have been computed with BSE and TDDFT. For InP, the TDDFT with the long range kernel is a satisfying approximation. For InSe, the presence of a bound exciton requires the BSE formalism. Due to the need of a surface reconstruction for InP, the size of the slab prevents me to do further calculation. For InSe, a thin slab has been calculated, within TDDFT with the long range kernel.

A short term perspective is to implement the calculation of the screening with the slab potential. We could also plan to use the EET formalism [89,90]. It would then offer the possibility to calculate the absorption spectrum within TDDFT with the nanoquanta kernel [91], and obtain absorption spectra for isolated slab taking benefit of TDDFT formalism.

Appendix A The slab potential A.1 Introduction

The slab potential has been defined to treat an isolated slab system in a framework of 3D code, based on planewaves and periodic boundary conditions. The three ideas are:

1) the response functions of the isolated slab are localised on the matter.

2) considering auxiliary response functions in real space equal to the response functions of the isolated slab, and periodic with the periodicity L mat z . 3) calculating the Fourier transform of the Dyson equation on the volume defined by the matter.

The periodicity of the auxiliary response functions implies that the reciprocal space vectors are multiple of 2π L mat z , the so-called "selected-G" vectors: the Dyson equation is solved only on the reduced basis of selected G-vectors and has to be modified accordingly:

where the slab potential is

The vectors of the reciprocal lattice are defined according to: G = ( Ḡ, Gz ), where Ḡ stands for the in-plane component, and Gz is a multiple of 2π/L mat z , with 2π/L mat z the thickness of the matter.

The first term is diagonal in G. It is the standard expression for the 3D-FT of the Coulomb potential.

109

The second term is a correction due to the finite thickness of the slab, proportional to 1/L mat z and depending on two G reciprocal lattice vectors. It is diagonal for the in-plane component Ḡ, which gives for the total potential the following dependence: Ṽ G1 G2 (q) = Ṽ Ḡ1 Ḡ2 ; Gz1 Gz2 (q) δ Ḡ1 , Ḡ2 . The second term also contains a phase factor ξ. It arises from the definition of the cell in the integration. When the integration is done between [-L mat z , 0] or [0, L mat z , 0], ξ = 1, while when the integration is done between [-L mat z /2, L mat z /2], ξ = (-1) n 1 +n 2 , where n i is the reduced coordinate of Gzi = 2πn L mat z . We will consider absorption spectrum and Electron Energy Loss spectra. Before presenting the results, we recall the main quantities of interest and the way they are evaluated within this formalism.

A.1.1 Absorption

The absorption spectrum is the imaginary part of the macroscopic dielectric function ε M (ω). It is calculated within TDDFT according to :

where ε M (q; ω) is given, for non local fields (NLF) and local fields (LF) cases respectively, by:

and V defined as V0 G2 = 0, ∀ G2 and V G1 G2 = Ṽ G1 G2 for G1 = 0. This leads to the following matrix representation for the slab potential:

A.1.2 EELS

The EEL spectrum is the imaginary part of the inverse dielectric matrix ε -1 G G(q; ω). For momentum transferred in the first Brillouin zone, this quantity is calculated, for NLF and NLF cases, according to: ε -1,N LF 00 (q; ω) = 1 + Ṽ00 (q) χNLF 00 (q; ω) (A.7) with χNLF 00 (q; ω) = χ 0 00 (q; ω) + χ 0 00 (q; ω) Ṽ00 (q) χNLF 00 (q; ω)

A.1.3 Remarks 1) Limits for vanishing q: lim q-→0 Ṽ00 (q) = 2πL mat z q valid for q sl < 10 -3 (A.9)

2) One should notice that with the relations Eq. (A.5) and ( A.8),

3) The slab potential, due to the presence of 1/|q|, which diverges at vanishing q, prevents the direct calculation of the component of the absorption perpendicular to the surface of the slab. Since ε M at vanishing q is a tensor, it is possible to extract this component from the linear combination of spectra calculated for in-plane and (in-plane + out-of-plane) components. 4) In the case of the EEL spectrum, ε -1 00 does not behave like a tensor at vanishing q, so it is not possible to extract a q z from any linear combination. The reason comes from the strong dependence of the spectrum with the length of q. Indeed, the dominant term of the potential is given by Eq. (A.9): it is proportional to 1/q. Since the response functions are still proportional to q 2 , it results from the definition of ε -1 00 (Eq. A.8) that the spectrum depends on the length of q. For very small q (< 10 -3 ), the amplitude of the spectrum is proportional to q. Then for larger q, the spectrum exhibit the expected dispersion, with a change of shape and energy position.

A.2 Effect of the slab potential

In order to study the behavior of the slab potential, I have studied two systems: The calculations are done using ξ = (-1] n 1 +n 2 .

A.2.1 Silicon slabs

The calculation are done on a very coarse grid (4x2x1): this explains that the amplitude of the spectra is not correct, but it does not affect our purpose. The parameters used to obtain a converged spectrum are: npwwfn = 8217 ; nbnd = 350 ; npwmat_xy = 7 ; npwmat_z = 17.

Absorption

Results are contained in files named ou4xxx.MDF:

range of q

Figure A.1 presents the Im(ε LF M ) for different length of q vector (q sl ). The black lines are spectra for q x component, the red lines spectra for q y component, the green lines are spectra for q z component (extracted from linear combination of spectra calculated with q x and q x ± q z ), and in blue spectra for q z component (extracted from linear combination of spectra calculated with q y and q y ± q z ). q sl = 10 -2 a.u. ou4MDF (y) q sl = 10 -2 a.u. ou4MDF (z) q sl = 10 -2 a.u. ou4MDF (z) q sl = 10 -3 a.u. ou4MDF (x) q sl = 10 -3 a.u. ou4MDF (y) q sl = 10 -3 a.u. ou4MDF (y) q sl = 10 -3 a.u. ou4MDF (z) q sl = 10 -4 a.u. ou4MDF (x) q sl = 10 -4 a.u. ou4MDF (y) q sl = 10 -4 a.u. ou4MDF (z) q sl = 10 -4 a.u. ou4MDF (z) q sl = 10 -5 a.u. ou4MDF (x) q sl = 10 -5 a.u. ou4MDF (y) q sl = 10 -5 a.u. ou4MDF (z) q sl = 10 -5 a.u. ou4MDF (z)

Figure A.1: Im(ε LF M ) for different q sl (ou4lf.MDF): (black) spectra for q x component; (red) spectra for q y component; (green) spectra for q z component (extracted from linear combination of spectra calculated with q x and q x ±q z ), and (blue) spectra for q z component (extracted from linear combination of spectra calculated with q y and q y ± q z ).

Results are independent of |q|, at least for the in-plane components, which allows us to conclude that χ00 is proportional to q 2 . The q z component extracted from the different linear combinations is, as expected, the same. This is an indication of the numerical stability of the procedure. The z-spectra calculated for q sl = 10 -2 a.u. is different from the one calculated for smaller values. Since the in-plane spectra exhibit no difference, it comes from the q x or y ± q z spectra. Nevertheless, q sl = 10 -2 a.u. enters in the limit of the k.p theory. ) [lines without dots (ou4nlf.MDF)] compared with Im(ε LF M ) [lines with dots (ou4lf.MDF)] for q sl = 10 -5 a.u.: (black) spectra for q x component; (red) spectra for q y component; (green) spectra for q z component (extracted from linear combination of spectra calculated with q x and q x ± q z ), and (blue) spectra for q z component (extracted from linear combination of spectra calculated with q y and q y ± q z ).

Effects of Local Fields

comparaison with surface

The surface spectra are calculated with Selected-G procedure, but using the limit L mat z -→ ∞ in Eq. (A.2), which reduced to the standard reciprocal space expression of 3D Coulomb potential:

Surface spectra are contained in ou5xxx.MDF files. ) for the slab [lines without dots (ou4lf.MDF)] and the surface [lines with dots (ou5lf.MDF)] q sl = 10 -5 a.u.: (black) spectra for q x component; (red) spectra for q y component; (green) spectra for q z component (extracted from linear combination of spectra calculated with q x and q x ± q z ), and (blue) spectra for q z component (extracted from linear combination of spectra calculated with q y and q y ± q z ).

EELS

Results are contained in files named: ou4nlf.EEL ε -1,N LF 00 (q; ω) = 1 + Ṽ00 (q) χNLF 00 (q; ω) q sl = 10 -2 a.u. ou4EEL (y) q sl = 10 -3 a.u. ou4EEL (x) q sl = 10 -3 a.u. ou4EEL (y) q sl = 10 -4 a.u. ou4EEL (x) q sl = 10 -4 a.u. ou4EEL (y) q sl = 10 -5 a.u. ou4EEL (x) q sl = 10 -5 a.u. ou4EEL (y) ) (ou4lf.EEL) for different length of q vector (q sl ): (black) spectra for q x component; (red) spectra for q y component. Bottom: same spectra normalised with 2

The plasmon of the slab is completely different from the plasmon of the bulk, where the spectrum is peaked around is around 17 eV. Here, the spectra are located in an energy range close to the one of absorption (Fig. A.4). There is a dispersion. On the contrary to the absorption case, the slab potential plays a important role. The bottom panel, where spectra are normalized shows that the EEL spectrum reached a shape similar to the absorption for vanishing q sl (see Fig. A.2). The factor has been chosen in such a way that the spectra are normalized to the absorption spectrum. Since lim q-→0 Ṽ00 (q) = 4πL mat z q , and the term which multiply χ00 is 2π |q| 2 , the factor is equal to q sl = 10 -3 a.u. ou4lf.EEL (x) q sl = 10 -4 a.u. ou4lf.EEL (x) q sl = 10 -5 a.u. ou4lf.EEL (x) q sl = 10 -2 a.u. ou4nlf.EEL (x) q sl = 10 -3 a.u. ou4nlf.EEL (x) q sl = 10 -4 a.u. ou4nlf.EEL (x) q sl = 10 -3 a.u. ou4lf.EEL (x) q sl = 10 -4 a.u. ou4lf.EEL (x) q sl = 10 -5 a.u. ou4lf.EEL (x) q sl = 10 -2 a.u. ou4nlf.EEL (x) q sl = 10 -3 a.u. ou4nlf.EEL (x) q sl = 10 -4 a.u. ou4nlf.EEL (x) q sl = 10 -5 a.u. ou4nlf.EEL (x) * 4.871821 10 -2 / q sl ) (ou4nlf.EEL) for different length of q vector (q sl ): (black) spectra with LF for q x component; (red) spectra without LF for q x component. Bottom: same spectra normalised with 2

Effects of Local Fields

Local field effects are small but visible (Fig. A.5). The extra-terms in the potential play a role. In the NLF calculations, only the head of the slab potential in used. In the LF calculations, the terms of the slab potential among the head of the matrix potential enter in the χ inversion, and in the summation for ε -1 00 . To disantangle these two contributions, I plot for each q sl , the different quantities ( Nevertheless, in the optical limit, all these contributions become small, and it seems that spectra have the shape given by χ 0 00 (furthermore normalized by 4π q 2 or Ṽ00 ). To check this hypothesis, figure A.7 shows the different χ: χ 0 00 (red); the green line is χLF 00 , which enters into the EELS (first term); the black line is χ00 , which enters into the absorption.

q sl = 10 -3 a.u. In fact, neither EELS ( χLF 00 ) nor absorption ( χ00 ) are equal to χ 0 00 . The EEL spectrum is similar to the absorption for very small q sl . The effect seen on the green curve comes from the Ṽ00 term, since it is already visible on Fig. A.7 for q sl = 10 -2 a.u..

When q increases, the black curve does not change, in agreement with the fact that the absorption spectrum is independent of the size of q. It is still quite close to χ 0 00 , since the LF are small for in-plane component. ) with local fields for the slab (ou4lf.EEL) and the surface (ou5lf.EEL) for different q sl : (black) spectra of slab for q x component; (red) spectra of slab for q y component, (green) spectra of surface for q x component; (blue) spectra of surface for q y component.

Comparison with surface

The EELS of slab is clearly different from the EELS of the surface, even for in-plane component (Fig. A.8). The collective oscillation of electrons at plasmon frequency is affect by the thickness of the slab. The surface plasmon is independent of q in the range considered. The plasmon of surface is the same that the bulk one (here we only plot the pre-edge).

A.2.2 Graphene layers' slabs

The calculation are done on a grid 20x20x1 shifted grid, which is good enough to define the structure of the spectra, but not enough to avoid spurious oscillations. The parameters used to obtain a converged spectrum are: npwwfn = 6531 ; nbnd = 280 ; npwmat_xy = 19 ; npwmat_z = 49.

Absorption

comparaison with surface

The surface spectra are calculated with Selected-G procedure, but using the standard reciprocal space expression of 3D Coulomb potential, which appears as the limit L mat z -→ ∞ of the slab potential in Eq. (A.2) :

Surface spectra are contained in ou5xxx.MDF files. ) for the slab [lines without dots (ou4lf.MDF)] and the surface [lines with dots (ou5lf.MDF)] q sl = 10 -5 a.u.: (black) spectra for q x component; (red) spectra for q y component; (green) spectra for q z component (extracted from linear combination of spectra calculated with q x and q x ± q z ), and (blue) spectra for q z component (extracted from linear combination of spectra calculated with q y and q y ± q z ). G have a negligible effect when we solve the corresponding Dyson equation. It is slightly different from the Si slab, where the spectra calculated with the surface and the slab potentials were not the same for the q z component (Fig. A.3). The spectra are independent of the length of q.

EELS

Results are contained in files named: ou4nlf.EEL ε -1,N LF 00 (q; ω) = 1 + Ṽ00 (q) χNLF 00 (q; ω) q sl = 10 -2 a.u. ou4EEL (y) q sl = 10 -3 a.u. ou4EEL (x) q sl = 10 -3 a.u. ou4EEL (y) q sl = 10 -4 a.u. ou4EEL (x) q sl = 10 -4 a.u. ou4EEL (y) q sl = 10 -5 a.u. ou4EEL (x) q sl = 10 q sl = 10 -2 a.u. ou4EEL (y) q sl = 10 -3 a.u. ou4EEL (x) q sl = 10 -3 a.u. ou4EEL (y) q sl = 10 -4 a.u. ou4EEL (x) q sl = 10 -4 a.u. ou4EEL (y) q sl = 10 -5 a.u. ou4EEL (x) q sl = 10 -5 a.u. ou4EEL (y) * 3.9720368 10 -2 / q sl Figure A.12: Top: -Im(ε -1,LF 00 ) (ou4lf.EEL) for different length of q vector (q sl ): (black) spectra for q x component; (red) spectra for q y component. Bottom: same spectra normalized with 2 |q|L mat z .

The plasmon of the slab (Fig. A.12) is different from the plasmon of the bulk, but the effect is much less important than for the case of silicon. It comes from the fact the due to the two main structures of graphitic systems, arising from the so-called π and π + σ bands, the corresponding features in the EEL and absorption spectra are more or less in the same energy range. The left panel, where spectra are normalized by 2 |q|L mat z (the origin of the factor as been explained for the silicon case) shows that the dispersion is only visible for large q sl ; the EEL spectrum reaches a shape similar to the absorption for q sl < 10 -3 . It was not reach so rapidly for silicon slabs. ou4nlf.EEL : ε -1,N LF 00 (q; ω) = 1 + Ṽ00 (q) χNLF 00 (q; ω) with χNLF 00 (q; ω) = χ 0 00 (q; ω) + χ 0 00 (q; ω) Ṽ00 (q) χNLF 00 (q; ω) ou4lf.EEL : ε -1,LF 00 (q; ω) = 1 + G Ṽ0 G(q) χ G0 (q; ω) with χ G G (q; ω) = χ 0 G G (q; ω) + G1 G2 χ 0 G G1 (q; ω) Ṽ G1 G2 (q) χ G2 G (q; ω) The EELS is calculated with the same expression as ours:

with χ GG (q; ω) = χ 0 GG (q; ω)

Absorption is deduced from ε -1,LF Bibliography